Impact of age-specific immunity on the timing and burden of the next Zika virus outbreak

The 2015–2017 epidemics of Zika virus (ZIKV) in the Americas caused widespread infection, followed by protective immunity. The timing and burden of the next Zika virus outbreak remains unclear. We used an agent-based model to simulate the dynamics of age-specific immunity to ZIKV, and predict the future age-specific risk using data from Managua, Nicaragua. We also investigated the potential impact of a ZIKV vaccine. Assuming lifelong immunity, the risk of a ZIKV outbreak will remain low until 2035 and rise above 50% in 2047. The imbalance in age-specific immunity implies that people in the 15–29 age range will be at highest risk of infection during the next ZIKV outbreak, increasing the expected number of congenital abnormalities. ZIKV vaccine development and licensure are urgent to attain the maximum benefit in reducing the population-level risk of infection and the risk of adverse congenital outcomes. This urgency increases if immunity is not lifelong.


Introduction
outbreak seroprevalence in adults is in line with findings from seroprevalence studies from French Polynesia, Brazil, and Bolivia [25][26][27].
In this study, we used published data from the 2016 ZIKV epidemic in Managua and developed an agent-based model (ABM) to predict the evolution of age-specific protective immunity to ZIKV infection in the population of Managua, Nicaragua during the period 2017-2097. We assessed: 1) the risk of a future ZIKV outbreak; 2) the consequences of a future ZIKV outbreak on women of reproductive age; 3) the influence of loss of immunity on future attack rates; and 4) how vaccination could prevent future ZIKV outbreaks.

Modelling strategy
We assessed the consequences of future outbreaks of ZIKV infection in Managua, Nicaragua using a stochastic ABM. The model follows a basic susceptible-infected-recovered (SIR) framework and integrates processes related to ZIKV transmission, immunity, demography, adverse congenital outcomes and vaccination (Table 1). We parameterized the model based on published estimates or inferences from data about the 2016 ZIKV epidemic (Table 1, Supporting information S1). We considered different scenarios about the duration of immunity, the timing and scale of ZIKV reintroductions in the population, and the timing and scale of a hypothetical vaccination program targeted towards 15 year old girls. Table 1. Parameterization of the agent-based model.

Parameter
Comment Source

ZIKV epidemic parameters
Transmission rate a Inferred from the 2016 epidemic b [23] Recovery rate Inferred from the 2016 epidemic b [23] ZIKV immunity Initial immunity a Inferred from the 2016 epidemic b [23] Duration of immunity Lifelong or decaying with time 5 scenarios c

Exposure
Proportion of women in the first semester of pregnancy [28] Risk of microcephaly Upon infection during exposure time (3 levels of risk) [10,11]

Model structure
We simulated a population of 10,000 individuals for 80 years (2017-2097). We assigned agents' age and ZIKV infection status (susceptible S, infected I or immune R). Initial conditions reflected the situation in Managua, Nicaragua in 2017, when there was no documentation of active transmission. In the outbreak-free period, we only considered demographic and immunity processes: births, deaths, ageing and, if applicable, loss of immunity and vaccination. Given the scarcity of these events at the individual level, we selected a long time-step of seven days and stochastically applied the transition probabilities at each time step for each agent. After a given time, ZIKV-infected cases were reintroduced in the population. Upon reintroduction, the time step was reduced to 0.1 days, and we evaluated the epidemic-related transition probabilities: Susceptible agents may become infected at a rate β a I/N, where β a is the agedependent transmission rate and N the total population size. Infected individuals may recover with a rate γ. We ignored the influence of the vector population and assumed that the force of infection is directly proportional to the overall proportion of infected individuals. We allowed six months for the outbreak to finish after introduction. Simulations were conducted independently for each combination of scenarios and repeated 1,000 times. In the baseline scenario, we assumed no vaccination, no loss of immunity and a reintroduction of 10 infected individuals. We implemented the model in 'Stan' version 2.18 [30] and we conducted analyses with R version 3.5.1 [31]. The Bayesian inference framework Stan permits the use of probability distributions over parameters instead of single values, allowing for the direct propagation of uncertainty. Stan models are compiled in C++, which improves the efficiency of simulations. The algorithm in Supporting information S1 describes the ABM in pseudo code. The model code and data are available from http://github.com/ZikaProject/SeroProject. Parameterization ZIKV epidemic parameters. We inferred the probability distributions for the age-specific transmission rate β a and the recovery rate γ from data on the 2016 ZIKV epidemic in Managua, Nicaragua. We used surveillance data [23], which give weekly numbers of incident ZIKV infections, confirmed by RT-PCR (dataset A, n = 1,165), and survey data on age-stratified ZIKV seroprevalence, measured among participants of pediatric and household cohort studies in Managua during weeks 5-32 of 2017 (dataset B, n = 3,740 children and 1,074 adults) [23].
We conducted statistical inference using a deterministic, ordinary differential equation (ODE)-based version of the ABM with three compartments (S, I and R) and two age classes (a 2 {1, 2} corresponding to ages 0-14 and �15): We ignored demography in this model because it covers a short time span. We recorded the overall cumulative incidence of ZIKV cases using a dummy compartment: in order to compute the weekly incidence on week t: We fitted the model to weekly incidence data A using a normal likelihood after a squareroot variance-stabilizing transformation [32]: where ρ is a reporting rate parameter and σ an error parameter. In addition, we also fitted the model to the number of individuals with anti-ZIKV antibodies at the end of the epidemic by age group B a using a binomial likelihood: where B a the number of individuals with antibodies, n a is the sample size in each age group, and p a = R a (t end )/N a (t end ) the proportion of immune at the end of the epidemic. The full likelihood was obtained by multiplying Eqs 6 and 7. We chose weakly-informative priors for all parameters and fitted the model in Stan (Table 2). We describe the calculation of the basic reproduction number R 0 in Supporting information S1. We used one thousand posterior samples for β a and γ obtained by Hamiltonian Monte Carlo in the ABM model, ensuring the propagation of uncertainty of these parameters. In Supporting information S1 we provide a schematic representation of the models and the information flow. Parameter values can translate from deterministic to agent-based versions of an epidemic model if the time step is small [33], which was the reason for using a time step of 0.1 days. ZIKV immunity. We used the deterministic model, described in the previous section, to infer the proportion of people with protective immunity within each age group at the end of the 2016 epidemicp a . We used one thousand posterior samples ofp a in the ABM to allow the propagation of uncertainty. Protective immunity to ZIKV after infection was lifelong in our first scenario, so the reduction of the overall proportion of immune individuals in the population decreased only because of population renewal. Given the absence of evidence about the duration of immunity to ZIKV, we considered four scenarios assuming exponentially distributed durations of immunity with means of 15, 30, 60, 90, or 150 years. These values correspond to a proportion of initially immune agents that loses immunity after 10 years of 55%, 28%, 15%, 11% or 6%, respectively (Supporting information S1).
Demography. We based the initial age distribution of the population on data from the World Bank [34]. We used age-dependent death rates for 2016 from the World Health Organization [29]. For births, we computed a rate based on an average birth rate in Nicaragua of 2.2 births per woman, which was uniformly distributed over the female reproductive lifespan [28]. We defined the period of reproductive age between 15 and 49 years. The ageing process was linear, increasing the age of each agent by 7 days at each 7-day time step. ZIKV reintroduction. We reintroduced ZIKV in the population after a delay of d = {1, . . ., 80} years in independent simulations. We chose this approach rather than continuous reintroductions to remove some of the stochasticity and assess more clearly the association between immunity decay and risk of an outbreak. As the probability of an extinction of the outbreak depends on the number of ZIKV cases reintroduced in the population, we considered three different values for the seed (1, 5 or 10 cases) and compared the results (Supporting information S1). Simulations using continuous reintroductions each year are presented in the Supporting information S1.
Risk of adverse congenital outcomes. The estimated number of microcephaly cases resulting from the reintroduction of ZIKV depended on the exposure, i.e. the number of pregnant women infected by ZIKV during their first trimester, to which we applied three different levels of risk, based on published estimates [10,11]. We obtained the number of ZIKV infections among women aged 15-49 years from ABM simulations. As gender was not explicitly considered in the model, we assumed that women represented 50% of the population. We assumed a uniform distribution of births during the reproductive period, and considered that the first trimester constituted a third of ongoing pregnancies at a given time. We explored three different levels of risk of microcephaly in births to pregnant woman with ZIKV infection during the first trimester, as reported by Zhang et al. (2017), based on data from French Polynesia (0.95%, called low risk) and Brazil (2.19% and 4.52%, called intermediate and high risk, respectively) [6,10,11].
Vaccination. We examined the effects of a potential ZIKV vaccine, given to 15-yearold-girls. This vaccination strategy was used for rubella virus, which also causes congenital abnormalities, before the vaccine was included in the measles, mumps and rubella vaccine given in early childhood [35]. The main objective of vaccination would be the prevention of adverse congenital outcomes, including microcephaly. We simulated this intervention in the ABM, assuming vaccine implementation starting in 2021, 2025 or 2031. From that date, half of the agents reaching age 15, representing females, could transition to immune status R regardless of their initial status, with an effective vaccination coverage ranging from 20% to 80%.

Outcome analysis
From the simulations, we collected 1) the evolution of the age-specific ZIKV immunity in the population; 2) the attack rate resulting from the reintroduction of ZIKV at year d; 3) the age of newly infected individuals. We fitted a binary Gaussian mixture model to dichotomize the observed attack rates into either outbreaks or non-outbreaks. We defined the outbreak threshold as the 97.5% upper bound of the lower distribution. This corresponded to a threshold of 1%, so that attack rates �1% were considered as outbreaks. The age structure of newly infected individuals was used to compute relative risks of infection by age group.

Sensitivity analysis
We explored the effect of seasonality, of changes in vector density, of migration, and of an endemic circulation of ZIKV on our predictions regarding the attack rate and the proportion of introductions that result in an outbreak. Different scenarios, methods and assumptions are provided in Supporting information S2.

Immunity and population
In our forward simulations, the expected population size increased by 42% between 2017 and 2097. Under the assumption that ZIKV infection results in lifelong protective immunity, population renewal will create an imbalance in the proportion immune in different age groups. We expect the overall proportion of the population with protective immunity to have halved (from 48% to 24%) by 2051 and to be concentrated among the older age classes (Fig 2A). The 0-14 year old age group will become entirely susceptible by 2031 and the 15-29 year old age group by 2046.

Future risk of ZIKV outbreak
Reintroductions of ZIKV in the population of Managua are unlikely to develop into sizeable outbreaks before 2035, 24 years after the 2016 epidemic, assuming lifelong immunity for individuals infected in 2016 (Fig 2B). After this point, attack rates resulting from ZIKV reintroduction will rise steeply. By 2047, we predict that ZIKV reintroductions will have a 50% probability of resulting in outbreaks with attack rates greater than 1% (Fig 2C). In 2047 the median attack rate of successful introductions is 3.6% (IQR: 2.0-6.2).

Risk of infection and microcephaly births in women of reproductive age
The differences between age groups in both immunity and transmission will result in a disproportionate burden of infection in the 15-29 year age class. The relative risk of infection in this age group ranges from 1.2 to 1.6, compared with the general population if an outbreak occurs during the period 2032-2075 ( Fig 3A). As most pregnancies occur in this age group, these women are also the most likely to experience a pregnancy with an adverse outcome. The increased risk of infection in this group implies that the number of adverse congenital  Impact of age-specific immunity on the timing and burden of the next Zika virus outbreak outcomes resulting from a ZIKV outbreak during this period is likely to be higher than expected with a homogeneous distribution of immunity across ages. Assuming different values for the added risk of microcephaly after a ZIKV infection during the first trimester, we expect the mean number of additional microcephaly cases due to ZIKV infection resulting from the reintroduction of the virus in Managua, Nicaragua to reach 1 to 5 cases per 100,000 population in 2060 (Fig 3B).

Loss of immunity
If protective immunity to ZIKV is not lifelong, the time window before observing a rise in the attack rates resulting from ZIKV reintroduction will shorten (Fig 4A). For instance, if 15% of the those who were infected in 2016 lose their immunity after 10 years (a mean duration of immunity of 60 years), the time until the risk of outbreak upon reintroduction reaches 50% would be 14 years earlier (2033) than with lifelong immunity (2047). If 55% lose their immunity after 10 years (a mean duration of immunity of 15 years), in 2024, 50% of the introductions result in an outbreak, and the attack rate in 2047 is 47%. Loss of immunity over time would reduce the relative risk in the 15-29 year old age group (Fig 4B).

Targeted vaccination
The implementation of a vaccination program targeted towards 15 year old girls between 2021 and 2031 would reduce the risk of infection in women aged 15-29 years and would also indirectly reduce the overall risk of a ZIKV outbreak in the population (Fig 5). If effective vaccine coverage is 60-80% amongst 15 year old girls, the prolongation of herd immunity could effectively mitigate the overall risk of a ZIKV outbreak in the population. The reduction in the number of microcephaly cases would then exceed what would be expected by considering only the direct protection granted by a vaccine to future mothers. A later implementation of the intervention would be less effective, as it becomes more difficult to maintain the herd immunity ( Fig 5B).

Sensitivity analysis
We considered additional model features that may impact our predictions in a sensitivity analysis (Table 3 and Supporting information S2). Accounting for seasonality or for a future  increase in vector abundance would result in higher transmission rates. This would lead to a shorter time window until a rise in the risk of ZIKV outbreak, and higher overall attack rates. A future diminution of vector abundance would have the opposite effects. Human migration from rural areas to Managua, Nicaragua would lead to a sharper decline of protective immunity in the population, also lowering the time window before the next ZIKV outbreak. Finally, a continuous endemic circulation of ZIKV in the region would increase the probability of an outbreak early on and lead to more stochasticity.

Discussion
In this mathematical modelling study, we show that a new ZIKV outbreak in Nicaragua would affect proportionally more women in the young reproductive age range (15-29 years) than the general population, owing to the age-dependent infection pattern and population renewal. The risk of a new ZIKV outbreak in Nicaragua, after reintroduction, will remain low before 2035 because of herd immunity, then rise to 50% in 2047. If protective immunity to ZIKV decays with time, ZIKV recurrence could occur sooner. Timely introduction of targeted vaccination, focusing on females aged 15 years would both reduce the risk of adverse congenital outcomes and extend herd immunity, mitigating the overall risk of an outbreak and resulting in lower attack rates if an outbreak occurs.

Strengths and limitations of the study
A strength of our approach is that it allows for the propagation of uncertainty from the initial data into the risk assessment, by transferring the posterior distributions of the parameters from the deterministic model fitted to surveillance and seroprevalence data on the 2016 epidemic into the ABM used for simulations. Roche et al. showed that, when a sufficiently small time step was chosen, stochastic and deterministic models using the same parameter values led to similar results [36]. Additionally, we benefited from the availability of high quality data from population-based surveys that included participants from age 2 to 80 years in Managua, Nicaragua. The age-stratified seroprevalence data allowed us to investigate the risk in different age groups and better assess the evolution of the age-specific immunity, which is crucial when studying adverse congenital events caused by ZIKV infection during pregnancy. We chose a simple approach based on an SIR structure, similar to the model used by Netto et al. (2017), to focus on the dynamics of infection and immunity in the human population [26]. We did not model vector populations and behavior explicitly, as in some other studies [14,16,37]. This simplification limits the mechanistic interpretation of the epidemic parameters, but provides a phenomenological description of the transmission dynamics. We believe that this approach is appropriate because our main objective was to determine the risk of an outbreak after reintroduction of ZIKV, which is mostly influenced by the level of protective immunity in the human population. We acknowledge that the future occurrence of ZIKV in the area also depends on the presence of a competent vector. Our choice is supported by sensitivity analyses that show that more complex model structures (delayed SIR and Ross-MacDonald-type models) were not superior to a simple SIR structure in describing the 2016 ZIKV epidemic of Managua (Supporting information S1). Similarly, Pandey et al. (2013) showed that additional model complexity does not result in a better description of the dynamics of transmission of dengue virus (another Aedes-borne virus) in a human population compared with a SIR model [38]. In our model, the transmission rate (β a ) captures both human-mosquito and mosquito-human transmission; we assumed a constant transmission rate, as observed in the 2016 outbreak.
Despite having modeled the effect of migration on our predictions, uncertainty remains; factors such as the political instability in Nicaragua could drive migration and influence disease transmission, as we currently observe in Venezuela and bordering countries [39].

Interpretation in comparison with other studies
This study shows that the lower attack rate of ZIKV in children than in adults will hasten the emergence of a population that will be fully susceptible to infection, especially if immunity is not lifelong. The advantage of our approach is that we used the age-specific attack rates to model the processes of ageing in relation to protective immunity to ZIKV explicitly. Even with lifelong immunity, our model predicts that children aged 0-14 years will become entirely susceptible by 2031 and 15-29 year olds by 2046. In future outbreaks, the attack rate will then be highest amongst 15-29 year olds, including women who will be at risk of ZIKV infection in pregnancy. If immunity wanes, the time until the next ZIKV outbreak will be reduced and, in that case, the distribution of infection risk would be more equal across age groups (Fig 4). Several authors have studied the time to a next ZIKV outbreak, but none studied the effect of the loss of immunity over time in relation to age. Assuming lifelong immunity, our estimates of the time until the risk increases are similar to the 12-20 years before re-emergence estimated for French Polynesia [16]. Netto et al. (2017) used an SEIR model to show that in Salvador, Brazil, the effective reproduction number was insufficient to cause a new outbreak during the "subsequent years" [26]. Lourenço et al. (2017) showed the same for the whole of Brazil: herd immunity should protect the population from a new outbreak in the coming years [40]. Ferguson et al. (2016) concluded that the age distribution of future ZIKV outbreaks will likely differ and that a new large epidemic will be delayed for "at least a decade" [14].
Other ZIKV vaccination studies confirm our findings. However, they do not show the effect in risk groups nor assume herd immunity from previous outbreaks as we did; Durham et al. (2018) showed that immunizing females aged 9 to 49 years with a 75% effective vaccine and a coverage of 90%, would reduce the incidence of prenatal infections by at least 94% [41]. Similarly, Bartsch et al. (2018) showed that women of childbearing age or young adults would be an ideal target group for vaccination [42]. Valega-Mackenzie et al. (2018) formulated a vaccination model for ZIKV transmission that included mosquito and sexual transmission [43]. They found that vaccination works if high coverage is achieved, both when sexual transmission or vector-borne transmission is most important.

Implications and future research
Our finding that people in the 15-29 year age range are more at risk of infection implies that we expect a higher number of congenital abnormalities due to ZIKV infection. Thus, vaccine development efforts should be increased. Our conclusions are drawn based on data from Managua, Nicaragua, but should be relevant to many regions in the Americas and the Pacific that have documented high post-epidemic levels of seropositivity [25][26][27]. In regions where ZIKV has not yet caused an epidemic but competent vectors are present, vaccination would be in place as well. Further age-stratified seroprevalence studies, using sensitive and specific tests and with longitudinal follow-up, are needed to improve our understanding of ZIKV antibody distribution in populations and to quantify the duration of immunity. This information will provide important information to improve mathematical modeling of ZIKV risk.
ZIKV vaccine development faces considerable hurdles. First, the evaluation of vaccine efficacy has stalled because the reduced circulation of ZIKV has reduced the visibility of ZIKVassociated disease [22]. Second, it remains unclear if neutralizing antibodies induced by vaccination are sufficient to protect women against vertical transmission and congenital abnormalities [44]. Third, it is not clear whether or how vaccine-induced antibodies against ZIKV will cross-react with other flaviviruses. To move vaccine development forward, we need to find regions where disease will occur to be able to conduct trials. This requires identifying populations that are at risk, and implementing surveillance there. These can either be regions where ZIKV is endemic, or where ZIKV outbreaks are likely to occur; throughout the Americas, there might be regions that did not experience an outbreak, but do have suitable conditions such as competent vectors. Conducting vaccine trials in disease outbreaks is complex, but there are tools to facilitate planning [45]. ZIKV in an endemic setting, such as in Africa and Asia, could prove a suitable setting as well. However, ZIKV circulation in endemic setting is not well described and the occurrence of adverse outcomes in this context is less documented [9]. Further research in understanding the transmission of the virus in an endemic context is therefore needed. Similarly, we need to increase the understanding of changes over time in vector abundance and population composition, since these influence the risk of new outbreaks.

Conclusion
Preparedness is vital; the time until the next outbreak gives us the opportunity to be prepared. The next sizeable ZIKV outbreak in Nicaragua will likely not occur before 2035 but the probability of outbreaks will increase. Young women of reproductive age will be at highest risk of infection during the next ZIKV outbreak. Vaccination targeted to young women could curb the risk of a large outbreak and extend herd immunity. ZIKV vaccine development and licensure are urgent to attain the maximum benefit in reducing the population-level risk of infection and the risk of adverse congenital outcomes. The urgency of ZIKV vaccine development increases if immunity is not lifelong.