Short- and long-term impact of vaccination against cytomegalovirus: a modeling study

Background Infection with cytomegalovirus (CMV) is highly prevalent worldwide and can cause severe disease in immunocompromised persons and congenitally infected infants. The disease burden caused by congenital CMV infection is high, especially in resource-limited countries. Vaccines are currently under development for various target groups. Methods We evaluated the impact of vaccination strategies and hygiene intervention using transmission models. Model parameters were estimated from a cross-sectional serological population study (n=5179) and a retrospective birth cohort (n=31,484), providing information on the age- and sex-specific CMV prevalence and on the birth prevalence of congenital CMV (cCMV). Results The analyses show that vertical transmission and infectious reactivation are the main drivers of transmission. Vaccination strategies aimed at reducing transmission from mother to child (vaccinating pregnant women or women of reproductive age) can yield substantial reductions of cCMV in 20 years (31.7–71.4% if 70% of women are effectively vaccinated). Alternatively, hygiene intervention aimed at preventing CMV infection and re-infection of women of reproductive age from young children is expected to reduce cCMV by less than 2%. The effects of large-scale vaccination on CMV prevalence can be substantial, owing to the moderate transmissibility of CMV at the population level. However, as CMV causes lifelong infection, the timescale on which reductions in CMV prevalence are expected is in the order of several decades. Elimination of CMV infection in the long run is only feasible for a vaccine with a long duration of protection and high vaccination coverage. Conclusions Vaccination is an effective intervention to reduce the birth prevalence of cCMV. Population-level reductions in CMV prevalence can only be achieved on a long timescale. Our results stress the value of vaccinating pregnant women and women of childbearing age and provide support for the development of CMV vaccines and early planning of vaccination scenarios and rollouts.


Continuous model formulation
Let us denote by S i (a, t) the age density of seronegative (susceptible, S) individuals of sex i ∈ { , } at time t. Similarly we denote by L i (a, t) and B i (a, t) the age densities of latently infected (latent, L) with low antibody concentrations and latently infected with high antibody concentrations (boosted, B) of sex i at time t. The age densities of infectious individuals with primary acute infection (I 1 ), and re-infection or reactivation from the L and B classes (I 2 and I 3 ) are I i 1 (a, t), I i 2 (a, t) and I i 3 (a, t), respectively. The age density of individuals of sex i at time t is given by P i (a, t) = S i (a, t) + L i (a, t) + B i (a, t) + I i 1 (a, t) + I i 2 (a, t) + I i 3 (a, t). (1) The total size of the population of sex i at time t is obtained by integrating P i (a, t) defined by Eq. (1) over all ages where M is the maximum attainable age of the population of sex i. The age-structured CMV model can then be formulated as follows ∂I i 1 (a, t) ∂t + ∂I i 1 (a, t) ∂a = λ i (a, t)S i (a, t) − γ i 1 (a) + µ i (a) I i 1 (a, t), ∂L i (a, t) ∂t + ∂L i (a, t) ∂a = γ i 1 (a)I i 1 (a, t) − ρ i (a) + zλ i (a, t) L i (a, t) − µ i (a)L i (a, t) + + (1 − p LB )γ i 2 (a)I i 2 (a, t), (5) ∂I i 2 (a, t) ∂t ∂B i (a, t) ∂t The initial and the boundary conditions are I i 1 (0, t) = q M 0 φ i (a)[I 1 (a, t) + L (a, t) + I 2 (a, t) + B (a, t) + I 3 (a, t)]da, S i (a, 0) = S i 0 (a), In the model defined by Eqs.
(3)-(20) µ i (a) denotes the age-specific death rate of individuals of sex i; φ i (a) is the age-specific fertility rate of women; q ∈ [0, 1] is the probability of vertical transmission; λ i (a, t) is the force of infection for individuals of sex i and age a at time t; ρ i (a) is age-specific reactivation rate for individuals of sex i; z ∈ [0, 1] is the reduction in susceptibility to re-infection in latently infected individuals compared to seronegative individuals; p LB ∈ [0, 1] is the probability of progression from low to high antibody concentrations; γ i 1 (a), γ i 2 (a), and γ i 3 (a) are the rates of progression from acute infectious state to latent uninfectious state.
The dynamics of the population can be determined independently from the CMV dynamics. The equation for the age density of individuals of sex i at time t is obtained by adding Eqs. (3)-(8) Adding up the boundary, Eqs. (9)-(14), and initial, Eqs. (15)-(20), conditions we obtain where P i 0 (a) = S i 0 (a) + L i 0 (a) + B i 0 (a) + I i 10 (a) + I i 20 (a) + I i 30 (a). In the demographic steady state the reproduction rate of the population is 1, which can be expressed as follows 2 Model with discrete age groups In this section we derive from the initial boundary-value problem defined by Eqs.
(3)-(20) a system of 6 × 2 × n (# of variables × # of sexes × # of age groups) ordinary differential equations for the number of individuals of different types at time t in n age groups. The age groups are defined by age intervals [a k−1 , a k ] with a 0 = 0 < a 1 < a 2 . . . < a n−1 < a n = M , k = 1, . . . , n. Let us denote by S i k (t), L i k (t), B i k (t), I i 1,k , I i 2,k and I i 3,k the number of S, L, B, I 1 , I 2 , and I 3 individuals of sex i ∈ { , } in the k-th age group [a k−1 , a k ] as follows We assume that the rates of fertility, death, reactivation, progression from acute to latent state and the force of infection are constant in each age group, i. e.
We integrate the equation for S i (a, t), Eq. (3), from a 0 = 0 to a 1 a 1 Using Eqs. (26), (33) and (38) the above equation can be rewritten in the following form Using Eqs. (26)-(32) the boundary condition given by Eq. (9) becomes We further define the progression rate from age group k to age group (k + 1) as follows Finally, the equation for the number of susceptible individuals of sex i in the first age group [0, a 1 ] becomes Integrating Eq. (3) from a k−1 to a k , where k = 2, . . . , n, we obtain The above equation simplifies to Using Eq. (42) we finally arrive to the following equation where k = 2, . . . , n.
The equations for other variables can be obtained similarly. The final system of the ordinary differential equations describing the model with discrete age groups reads as follows dS i 1 (t) dt = n k=1 φ i k {S k (t) + (1 − q)[I 1,k (t) + L k (t) + I 2,k (t) + B k (t) + I 3,k (t)]} − − m 1 + λ i where k = 2, . . . , n.

Force of infection
The force of infection for individuals of sex i and age a at time t is expressed as follows where β 1 and β 2 are proportionality parameters determining the infectivities of primary infection and re-infection/reactivation, c ij (a, a ) is the contact rate between individuals of 7 sex i and age a and individuals of sex j and age a . More precisely, c ij (a, a ) represents the number of contacts per unit of time that one individual of sex i and age a receives from individuals of sex j and age a . The contact rate c ij (a, a ) can thus be written down as the product of the number of contacts per unit of time between one individual of sex i and age a and one individual of sex j and age a and the age density of individuals of sex j at time t c ij (a, a ) =c ij (a, a )P j (a , t). (60) Note that by definitionc ij (a, a ) is symmetric while c ij (a, a ) is not. Using the above equation the expression for the force of infection simplifies to Eq. (61) is used in the continuous formulation of the model defined by Eqs.
(3)-(8). In the model with discrete age groups defined by Eqs. (47)-(58), the force of infection has to be discretized. This is done assuming that the contact rate is constant for the interactions between age groupsc ij (a, a ) =c ij kl (62) for a ∈ [a k−1 , a k ] and a ∈ [a l−1 , a l ], k, l = 1, . . . , n.
The force of infection for individuals of sex i in the k-th age group [a k−1 , a k ] becomes Eq. (63) is used in the model with discrete age groups, Eqs. (47)-(58).

Simplified model with discrete age groups
We now aim to obtain a simplified version of the model with discrete age groups, Eqs.
(47)-(58) that can be fit to the cross-sectional serological data providing information on the infection status of the Dutch population in 2006/2007, and a retrospective cohort study carried out in 2008 providing information on the birth prevalence of cCMV (i.e. the fraction of infants infected during pregnancy) [33,34]. Firstly, we assume that the male to female ratio is 1 to 1. Secondly, the population is assumed to be in the demographic equilibrium, i.e. population sizes of different age groups are constant. Finally, we model Type I mortality, namely the probability of surviving is 1 till the maximum age M and zero above that. We start by obtaining equations for the number of individuals of sex i in the k-th age group [a k−1 , a k ] where k = 2, . . . , n.
We now find the expression for the progression rate m k from age group k to age group (k + 1), k = 1, . . . , n for the demographically stable population: P i (a, t) ≡ P i (a). For this special case Eq. (21) simplifies to Solving the above equation on the interval [a k−1 , a k ] results in Using this expression we further integrate P i (a) over the interval [a k−1 , a k ] to get The progression rate m k from age group k to age group (k + 1) then becomes 9 The final expression for m k is Note that Type I mortality implies that all µ i k = 0, therefore we must find the limit of m k when µ i k → 0: Finally, we substitute µ i k = 0 into Eqs. (65)-(66) and set dN i k (t)/dt = 0, k = 1, . . . , n. We obtain the following system of equations . . .
where k = 3, . . . , n. Solving Eqs. (76)-(78) we obtain the equation for N i k in terms of N i 1 : Since the ratio of men to women is 1 to 1, we can write Eq. (75) as follows Combining the last two expressions we obtain which can be written as The above equation is the discrete version of Eq. (24). Denoting the total population size for sex i as N i and using Eq. (79) we get If the population size of sex i, N i , the maximum attainable age, M , and the age group intervals [a k−1 , a k ] are given, the population sizes in different age groups are then determined as This means that the age distribution is uniform on the interval [a 0 = 0, a n = M ].
Since by construction the population sizes of different age groups do not change we can express S i The equations for the model with discrete age groups are where k = 2, . . . , n.

In the above equations
and The steady state solutions of the model are obtained by setting the left hand sides of Eqs. (85)-(94) to zero. To calculate the disease free steady state we substitute in the resulting set of equations q = 0, λ i k (t) = 0 and L i where the bar denotes the steady-state solution.

Calculation of the effective reproduction number
In this section we describe how to calculate the basic reproduction number, R 0 , in the absence of vaccination and the effective reproduction number, R e , in the population with vaccination or hygienic intervention. We used the method known in the literature as the next-generation matrix approach. The detailed description of the method can be found in [48,49]. The calculation of R e follows the same steps, with the only difference that the starting point is the system of differential equations for the CMV model with vaccination. In the following sections we give the model equations for different vaccination scenarios and the respective disease free equilibrium.  Figure S1: Schematic of the model with universal vaccination with a vaccine preventing infection. The vaccine is assumed to protect against primary infection in seronegative individuals. Proportion p of susceptible (S) individuals in age group k is effectively vaccinated (V), where p is given by the product of vaccination coverage and vaccine efficacy. Vaccinated individuals lose protection at rate δ (average duration of protection, 1/δ), returning to the susceptible class.

Model with universal vaccination
We considered a suite of universal vaccination strategies with varying proportions of effectively vaccinated persons, ages at vaccination, sexes to be vaccinated, and durations of protection. We distinguished between scenarios in which the vaccine is assumed to protect only against primary infection in seronegative individuals ("prevention of infection") or against primary infection in seronegative individuals and re-infection/reactivation in seropositive individuals ("prevention of (re-)infection and reactivation"). The target population for vaccination was either infants in the first year of life, adolescent boys and girls at the age of 10 years, adolescent girls at the age of 10 years, or women of reproductive age (15-50 years).

Prevention of infection
In the model with universal vaccination with a vaccine preventing infection, the proportion p of susceptible (seronegative, S) individuals in age group k is effectively vaccinated (V). Henceforth, we will refer to p as the effectively vaccinated proportion, which is the product of vaccination coverage and vaccine efficacy. Vaccinated individuals lose protection at rate δ (average duration of protection, 1/δ), returning to the susceptible class. The special case of δ = 0 describes the model with life-long protection. Schematic of this model is given in Figure S1. For different scenarios, the vaccinated age groups are k = 2 (6-months-old), k = 4 (10-years-old), k = 7 (25-years-old).

6-months-old boys and girls and life-long protection
In the scenario with vaccination of infants in the first year of life, the fraction p of susceptible boys and girls (S) is vaccinated (V) at 6 months of age (k = 2). In the simplest case when the vaccine-induced protection does not wane, these individuals stay protected for the rest of their lives. Let us denote V i k (t) denote the number of vaccinated individuals of sex i in the age group k at time t. The model equations for this scenario read as follows where k = 2, . . . , n in the equations for I i 1,k (t), L i k (t), I i 2,k (t), B i k (t) and I i 3,k (t). The number of individuals of sex i in the k-th age group [a k−1 , a k ] can be written down as follows: To calculate the effective reproduction number, R e , we expressed S i and dropped the equations for S i k (t). The equations for the calculation of R e are where k = 2, . . . , n in the equations for and The disease free equilibrium is 6.1.2 6-months-old boys and girls and waning protection If the protection wanes with rate δ, the equations for susceptible and vaccinated persons are changed: The equations for calculation of the effective reproduction number are The disease free equilibrium is , k = 3, . . . , n,

Individuals in age class k and waning protection
The equations for the model with vaccination of seronegative individuals in age class k and waning protection ( Figure S1) are (109) where k = 2, . . . , n in the equations for . The model equations for the calculation of R e (without the equations for susceptibles): where k = 2, . . . , n in the equations for L i k (t), I i 2,k (t), B i k (t) and I i 3,k (t).
The disease free equilibrium for the calculation of R e is

Prevention of (re-)infection and reactivation
In the model with universal vaccination with a vaccine preventing (re-)infection and reactivation, the proportion p of seronegative (S) and seropositive individuals (L, B, I 1 , I 2 , and I 3 ) individuals in age group k is effectively vaccinated. The seronegative and seropositive vaccinated individuals transit to class V 1 and V 2 , respectively. As before, p is the effectively vaccinated proportion, which is the product of vaccination coverage and vaccine efficacy. Vaccinated individuals (V 1 and V 2 ) lose protection at rate δ (average duration of protection, 1/δ), returning from class V 1 to the susceptible class (S) and from class V 2 to the latent class with low antibody concentrations (L). Schematic of this model is given in Figure S2. For different scenarios, the vaccinated age groups are k = 2 (6-months-old), k = 4 (10-years-old), k = 7 (25-years-old).

Individuals in age class k and waning protection
The model equations are  Figure S2: Schematic of the model with universal vaccination with a vaccine preventing (re-)infection and reactivation. The vaccine is assumed to protect against primary infection in seronegative individuals and re-infection/reactivation in seropositive individuals. In this vaccination scenario, the proportion p of susceptible individuals (S) in age group k is effectively vaccinated (V 1 ). The proportion p of all seropositive persons (L, B, I 1 , I 2 , and I 3 ) in age group k is vaccinated (V 2 ) as well. After the protection is lost (average duration: 1/δ), V 1 and V 2 individuals return to the S and L class, respectively. (112) where k = 2, . . . , n in the equations for The equations for the calculation of R e (without the equations for susceptibles) are The disease free steady state for the calculation of R e is

Model outcomes
We evaluated the impact of all interventions on the percent reduction in the incidence of cCMV, primary infections and re-infections/reactivations after a time frame of 20 years. The incidence of primary infections at time t was computed as 24 The incidence of re-infections/reactivations at time t was computed as n k=1 i∈{ , } The incidence of congenital infections at time t was computed as where 1/γ is the duration of acute infection.

Parameter summary
The description of the model parameters is given in Table S1. probability of progression from low to high antibody concentrations estimated 1/γ i 1,k = 1/γ days 1/γ i 2,k = 1/γ days duration of acute infection 14 days, sensitivity analyses 1/γ i 3,k = 1/γ days β 1 1/year infectivity of primary infection estimated β 2 1/year infectivity of re-infection/reactivation estimated p effectively vaccinated proportion (vaccine efficacy × vaccination coverage) 0-100%, sensitivity analyses 1/δ years duration of protection 2.5-80 years, sensitivity analyses  pregnancy) [33,34]. Parameters were broadly estimated using methods developed earlier [41]. Briefly, assuming an endemic equilibrium and the short disease approximation, estimates were obtained using a Bayesian framework. In comparison with our earlier analyses, we extended the model by (i) allowing for multiple reactivation and re-infection events occurring over a person's life ( Figure 1 in the main text), (ii) using cubic B-splines for more flexible estimation of the age-dependent reactivation rates, and (iii) including the birth cohort data to enable estimation of the probability of cCMV. For horizontal transmission we used an age-and sex-specific contact matrix with 17 age classes. The model is fitted to the data using Hamiltonian Monte Carlo method as implemented in Stan (https://mc-stan.org). Data, model code, and Rmarkdown file are available at https://github.com/mvboven/cmv-vaccination. The Rmarkdown file with prior specifications and parameter settings is attached at the end of the Supplementary Material. Figure S3 shows an overview of all parameter estimates which are not given in the main text.
10 Results for the scenario with high reactivation rates Estimates of the reactivation rates, infectivity of re-infection and reactivation, and probability to move from the L to B class depend sensitively on the assumed prior distributions of the reactivation rates. In the main analyses, we took prior distributions such that a priori reactivation rates range from 0.005 to 0.11 per year (95% range), which seemed reasonable given what we had found earlier [41].However, with our more flexible model which allows for multiple reactivation and re-infection events in the infected classes (L and B), the reactivation rate estimates are found to depend sensitively on the prior distributions ( Figure S4). For instance, with broader and higher prior distributions for these rates, the posterior distributions of the reactivation rates are also (markedly) higher, while the infectivity and probability to move from the L to the B class are (markedly) lower ( Figure S5). Fortunately, due to the neutralizing effects of the parameter shifts, the impact of vaccination, in particular the estimated reductions in cCMV, primary infections, re-infections, and reactivations are quantitatively very similar to the main scenario and sensitivity analyses (cf. Table 1 in the main text and Table S2). Hence, while considerable uncertainty surrounds the actual magnitude of the reactivation rates and associated parameters, the impact on the effectiveness of vaccination is small.  Figure S4: Overview of the estimation results in the scenario with high reactivation rates. (A) and (B) show the age-specific prevalence of seronegative, seropositive with low antibody concentrations and seropositive with high antibody concentrations in females (A) and males (B), respectively. (C) shows the age-specific reactivation rates per year for females (black) and males (blue), respectively. The solid lines are the median estimates, and the shaded regions correspond to 95% credible intervals obtained from 1000 parameter samples from the posterior distribution. Note that the seroprevalence is estimated with high precision, and that credible intervals for the reactivation rates are quite broad.   Table S2: Impact of interventions on cCMV, primary infection and reinfection/reactivation in the scenario with high reactivation rates. The reductions are evaluated 20 years after the start of the intervention. The proportion of effectively vaccinated persons (vaccination coverage × vaccine efficacy) is 70%, and the average duration of protection is 10 years. Hygienic measures are modeled as a 70% reduction in infectious contacts between women of reproductive age (15-50 years) and young children (0-5 years). The effective reproduction number is defined as the average number of secondary infections at the start of an epidemic with one infected individual introduced in a population where 70% of persons are effectively vaccinated. This number smaller than 1 indicates that a given intervention is going to lead to the disease elimination in the long run.

30
To compute the disease burden caused by cCMV in the Netherlands, we used the yearly incidence of cCMV predicted by the model. From the yearly incidence we computed numbers of cases of cCMV per year using a population size of 17 million for the Netherlands in 2017. Furthermore, we used an estimate of the disability-adjusted life years (DALYs) per case of 3.034 (95%CrI: 1.202-6.105) from the meta-analysis for Belgium [13].The DALY per case result is not reported in the published paper, but was provided to us by Brecht Devleesschauwer (personal communication). For every vaccination scenario, we computed the yearly incidence for a period of 20 years, calculated the number of cCMV cases based on the population size of 17 million, and computed the number of DALYs per year by multiplying with the DALY per case estimate. We then calculated the number of DALYs prevented in subsequent years compared to reference year zero and added those prevented DALYs for the period of 20 years to arrive at the number of DALYs prevented over that time period. We did the computation for the median incidence and the 95% credible interval of the incidences (which is based on the uncertainty of model parameters, but not the uncertainty in DALY per case estimate). Table S3 shows the computed numbers for one scenario, namely the vaccination of pregnant women. Figure S6 shows the DALYs prevented for all scenarios together with the 95% credible intervals based on the uncertainty in incidences.
Year  T  h  e  a  n  a  l  y  s  e  s  a  n  d  r  e  s  u  l  t  s  p  r  e  s  e  n  t  e  d  h  e  r  e  f  o  r  m  t  h  e  s  t  a  r  t  i  n  g  p  o  i  n  t  o  f  a  m  a  n  u  s  c  r  i  p  t  t  h  a  t  e  v  a  l  u  a  t  e  s  t  h  e  i  m  p  a  c  t  o  f  v  a  r  i  o  u  s  v  a  c  c  i  n  a  t  i  o  n  s  t  r  a  t  e  g  i  e  s  a  g  a  i  n  s  t  C  M  V  .  A  u  t  h  o  r  s  o  f  t  h  e  m  a  n  u  s  c  r  i  p  t  a  r  e  G  a  n  n  a  R  o  z  h  n  o  v  a  ,  M  i  r  j  a  m  K  r  e  t  z  s  c  h  m  a  r  ,  F  i  o  n  a  v  a  n  d  e  r  K  l  i  s  ,  D  e  b  b  i  e  v  a  n  B  a  a  r  l  e  ,  M  a  r  j  o  l  e  i  n  K  o  r  n  d  e  w  a  l  ,  A  n  n  V  o  s  s r  e  c  o  n  s  t  a  n  t  o  n  t  h  e  i  n  t  e  r  v  a  l  s  .  T  h  e  m  e  s  h  c  a  n  b  e  m  a  d  e  e  v  e  n  m  o  r  e  f  i  n  e  g  r  a  i  n  e  d  ,  o  r  i  n  t  e  r  p  o  l  a  t  i  o  n  s  c  a  n  b  e  u  s  e  d  t  o  m  a  k  e  t  h  e  l  i  k  e  l  i  h  o  o  d  c  o  n  t  r  i  b  u  t  i  o  n  s  o  f  t  h  e  s  e  r  o  l  o  g  i  c  a  l  d  a  t  a  m  o  r  e  p  r  e  c  i  s  e  .  I  n  o  u  r  e  x  p  e  r  i  e  n  c  e  ,  t  h  e  r  e  i  s  l  i  t  t  l  e  a  d  d  i  t  i  o  n  a  l  p  r  e  c  i  s  i  o  n  t  o  b  e  g  a  i  n  e  d  b  y  s  u  c  h  a  p  p  r  o  a  c  h  e  s  .  S  e  e  b  e  l  o  w  f  o  r  c  o  d  e  r  e  l  a  t  i  n  g  t  o  t  h  e  f  i  r  s  t  t  w  o  s  t  e  p s .