Accounting for herd immunity in projections of the population at risk

To provide a more realistic estimate of the population at risk of ZIKV in Asia, we applied a method that accounts for the protective effect of herd immunity and also considered how pre-existing immunity to ZIKV among Asian populations may impact these projections. Having previously used this method to make projections for the Americas,16 here we made analogous projections for 15 Asian countries with a significant history of DENV transmission.3 This method combines classic results from epidemiological theory with empirical estimates of infection attack rates (IAR) from past epidemics of ZIKV and CHIKV and highly spatially resolved data on drivers of transmission to make location-specific projections of IAR in the event of an epidemic in an immunologically naïve population.

Data to inform projections of the population at risk

We first assembled a suite of spatial data layers obtained from different sources to serve as covariates in our spatially explicit model that focused on ZIKV transmission by the Ae. aegypti mosquito, given its dominant role in recent outbreaks and its high potential to contribute to future transmission in Asia.17 18 For demography, we downloaded high-resolution gridded data from WorldPop (http://worldpop.org.uk) on total population and total births in 2015 (at 100 m x 100 m resolution) covering 15 countries: Bangladesh, Brunei, Cambodia, India, Indonesia, Laos, Malaysia, Myanmar, Pakistan, Philippines, Singapore, Sri Lanka, Timor-Leste, Thailand, and Vietnam.19 For five countries that did not have gridded data on births (Brunei, Malaysia, Singapore, Sri Lanka, Thailand), we imputed births at each grid cell by scaling total population in each grid cell by national-level birth rates.20 For temperature, we downloaded WorldClim (www.worldclim.org) gridded mean monthly temperature at a resolution of 5 km x 5 km covering all 15 countries.21 We then calculated two location-specific transmission parameters as a function of temperature: daily mortality of Ae. aegypti 22 and extrinsic incubation period based on studies of DENV infections in Ae. aegypti.23 Similarly, we obtained a set of 100 gridded Ae. aegypti occurrence probability layers at a resolution of 5 km x 5 km, with each layer reflecting a unique sample from the posterior distribution of parameters.24 To account for heterogeneity in mosquito-human contact due to socioeconomic factors,25 we incorporated estimates of per capita gross cell product (GCP) adjusted for inflation at a resolution of one degree.26 Missing data for two countries (Timor-Leste and Myanmar) were uniformly filled with inflation-adjusted GCP figures obtained from the United States Central Intelligence Agency World Factbook.27

Modelling to inform projections of the population at risk

We used the aforementioned spatial data layers to calculate the basic reproduction number R ₀ as a function of temperature T and the GCP economic index based on the classic Ross-Macdonald form,28

(1)

where a is mosquito biting rate, m is the ratio of mosquitoes to humans, α(GCP) is a scaling factor that incorporates the effect of GCP on mosquito-human contact, b and c are human-to-mosquito and mosquito-to-human infection probabilities, μ(T) is the daily mortality rate of adult mosquitoes at temperature T, 1/r is the average duration of infectiousness in humans, and n(T) is the extrinsic incubation period at temperature T. Because of our assumption of a common dominant vector in Ae. aegypti, we assumed that this model for R ₀ was a reasonable approximation for any of ZIKV, DENV, or CHIKV in settings where those viruses are transmitted predominantly by Ae. aegypti, although there may be important ways in which these viruses differ.29

Based on a theoretical one-to-one relationship between the basic reproduction number, R ₀, and the first-wave infection attack rate, IAR, under an SIR transmission model,14 we generated location-specific projections of IAR. Lesser known parameters in the formulas for R ₀ and IAR (ie, the effect of GCP and the scaling relationship between R ₀ and IAR) were calibrated by comparing model-generated IAR with 13 published estimates of IAR experienced during invasions of CHIKV and ZIKV into naïve populations, as determined by virus-specific serological tests in 12/13 cases and an extrapolation of clinical attack rate in 1/13 cases (table 1).30–37 Although there are limitations associated with serological testing for ZIKV,38 these data presently reflect the best estimates available of the epidemic size of ZIKV or related viruses in naïve populations in areas with high-quality gridded spatial data on drivers of transmission. To generate projections of the overall population and annual births at risk of ZIKV infection, we multiplied location-specific infection attack rate projections, IAR, by the corresponding population and births in each grid cell.

Because they do not account for spatial processes such as ZIKV importation and the relative isolation of some locations, our IAR projections are best interpreted as location-specific projections conditional on an epidemic occurring in the first place. Our projections are also not temporally specific, meaning that they do not take into account the timing of ZIKV introduction in relation to seasonal transmission patterns, and they are not specific about the timeframe over which an epidemic will subside. They also focus exclusively on the potential for ZIKV transmission by Ae. aegypti mosquitoes and may underestimate the contributions to ZIKV transmission by other mosquito vectors or by other transmission routes.17 Despite these limitations, our model-based approach is one of the few modelling tools currently available for making projections of ZIKV epidemic size.12

Validating model projections of the population at risk

Since our initial publication16 describing the method used here, four additional studies39–42 have been published that offer an opportunity to validate our method against empirical estimates of IAR for CHIKV and ZIKV (table 2). One estimate of IAR for CHIKV following an epidemic in Puerto Rico was based on IgM and IgG ELISAs performed on blood donations before and after the CHIKV epidemic there.41 Another estimate of IAR for CHIKV following an epidemic in Chapada, Brazil was also based on IgM and IgG ELISAs and was performed on a random sample of residents.39 One estimate of IAR for ZIKV over the course of an epidemic in Guayaqil, Ecuador was based on biweekly PCR assays for ZIKV from control subjects in a case-control study of pregnant women.42 Another estimate of IAR for ZIKV over the course of an epidemic in Recife, Brazil was also based on subjects in a case-control study of pregnant women but used plaque-reduction neutralization tests (PRNT).40 For each of these four studies, we obtained a probabilistic estimate of IAR under a Bayesian framework by assuming a non-informative beta prior probability distribution and a binomial likelihood, which implied a beta posterior probability distribution with parameters and .43

We compared our location-specific IAR projections to these four empirical estimates of CHIKV and ZIKV IAR that were independent from the data used to calibrate the model. To make IAR projections comparable to these estimates, we first compiled data on temperature, Ae. aegypti occurrence probability, and the GCP economic index for each site, and we then computed a distribution of IAR projections for each site by applying 1000 random draws of model parameters to equation 1. We found that the 95% posterior predictive intervals of three recent empirical estimates of CHIKV and ZIKV IAR intersected with the range of IAR projections (figure 1, Guayaquil, Puerto Rico, Chapada). However, the 95% posterior predictive interval of one recent empirical estimate of ZIKV IAR40 was noticeably lower than the corresponding empirical estimate (figure 1, Recife). This resulted in location-specific IAR projections averaged across the four new sites being somewhat lower than empirical estimates of IAR averaged across those same four sites (mean of projected IARs=0.31, mean of empirical IARs=0.40) as compared to the original 13 sites (mean of projected  IARs=0.27, mean of empirical IARs=0.23). Based on this exercise, we conclude that our projections may be relatively accurate in some areas but could represent an underestimate in other areas, such as Recife, where the epidemic was particularly explosive.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Comparison of model projections of infection attack rate (IAR) (blue and red) and empirical estimates of IAR (black) for each of 17 sites. Coloured violin plots show the distribution of 1,000 Monte Carlo samples of model projections of IAR at the 13 sites used to calibrate the model (blue) and at the four sites used to independently validate the model (red). Black lines show the 95% posterior predictive intervals of empirical estimates of IAR of either CHIKV or ZIKV at those sites, and the points shows the respective medians.

An upper limit for the projected population at risk

Although it is important to account for the limiting effects of pre-existing immunity on ZIKV projections in Asia, we first describe projections made under the assumption of no pre-existing immunity. This is necessary for two reasons. First, although these projections are almost certainly an overestimate across the region, they nonetheless provide an upper limit for the number of possible ZIKV infections across the region and are necessarily much lower than estimates by Bogoch et al 13 due to our acknowledgement of herd immunity. Second, given that historical patterns of ZIKV transmission in Asia may have been heterogeneous in time and space,4 11 the assumption of little to no pre-existing immunity in our location-specific projections could be reasonable for some locations.44

Our upper-limit projections indicate that as many as 785 (range: 730–992) million people could become infected by ZIKV in the event of a relatively rapid first-wave epidemic in an immunologically naïve population, were an epidemic to occur in Asia under circumstances similar to the epidemic in the Americas (figure 2). Because as many as 80% of ZIKV infections may be asymptomatic,37 we note that our implied upper-limit projection for the number of clinical episodes of Zika is considerably lower (ie, 157 million assuming that 20% of infections are symptomatic). At the country level, we project that India has the potential for the largest number of ZIKV infections by more than six-fold that of any other country, mainly due its large population and relatively high suitability for ZIKV transmission. The country-level, per capita probability of infection ranges from 0.21 in Laos to 0.40 in Cambodia and Bangladesh. Among childbearing women across the region as a whole, there could be as many as 15.3 (range: 14.2–19.6) million infected under this worst-case scenario. India would likewise be projected to have the highest number of births affected by Zika under this scenario, but countries such as Pakistan and the Philippines would have disproportionately more infections among childbearing women due to those countries’ relatively high birth rates.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Model-based projections of ZIKV epidemic potential in Asia. (A) Median projections of first-wave ZIKV infection attack rates at 5 km x 5 km resolution across Bangladesh, Brunei, Cambodia, India, Indonesia, Laos, Malaysia, Myanmar, Pakistan, Philippines, Singapore, Sri Lanka, Timor-Leste, Thailand, and Vietnam, using the approach by Perkins et al.16 These projected infection attack rates were combined with spatial demographic data19 to project (D) expected numbers of infections for men and women of all ages and then (B, C) summed nationally. These location-specific projections are shown in greater detail for (E) western Java, (F) southern Vietnam, (G) Bangladesh, and (H) Singapore and Kuala Lumpur. Tables indicate median values and the total range of values from 1,000 Monte Carlo samples that reflect uncertainty in model parameters.

Compared with our projections for the Americas,16 these totals are nearly one order of magnitude greater due to a combination of larger populations and higher projected infection attack rates. India alone has more than twice the population of the region of the Americas we considered in our previous analysis, yet it is projected to have a five-fold larger number of infections in the event of a first-wave ZIKV epidemic in an immunologically naïve population. In addition, country-level averages of projected infection attack rates are also higher in several cases in Asia than they are in the Americas.

Putting the projections into context

To understand the extent to which our projections may overestimate the true population at risk of ZIKV infection due to the fact that our projections do not account for pre-existing immunity, it is important to understand how a given level of pre-existing immunity would impact our projections. Specifically, for a given level 1-s of pre-existing immunity in a population with a projected infection attack rate of IAR, what does this imply about infection attack rates IAR _i in the presence of pre-existing immunity? Pre-existing immunity would mean that the effective reproduction number R assumes a fraction S of its value R ₀ in a fully susceptible population,9 the latter being what our projections of IAR are based on. Such a linear effect of S on R implies a nonlinear effect of S on IAR _i . Specifically, in regions with the highest R ₀ we would expect IAR _i > (1-s) IAR among the susceptible population, whereas in regions with the lowest R ₀ we would expect IAR _i < (1-s) IAR among the susceptible population, possibly even 0 (figure 3). As a result, pre-existing immunity in high-transmission settings in Asia may result in less indirect protection of susceptibles than the simple notion of a linear protective effect of pre-existing immunity might suggest.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

The effect of pre-existing immunity on infection attack rates (IAR). The nonlinear relationship between the basic reproduction number R ₀ and the infection attack rate IAR implies that the same level of pre-existing immunity 1-s could have different impacts on future epidemics of ZIKV in Asia, depending on location-specific values of R ₀ or, in other words, the overall potential for transmission in a given area. Consider a population with 50% pre-existing immunity. The reproduction number R = (1-s) R ₀=R ₀/2. For locations with high R ₀, this reduces IAR to IAR_i by a relatively small fraction (eg, 19%). For locations with low R ₀, this reduces IAR to IAR_i by a relatively large fraction (eg, 92%). In both scenarios, IAR_i is the projected infection attack rate among individuals susceptible at the onset of the epidemic in question.

Although it is clear that ZIKV transmission has been documented numerous times in Asia both before and after its emergence in the Americas,44 it is not clear how different the situation in Asia may be now than it was before 2015. Travel histories of recently documented ZIKV infections in Asia indicate that there has been a mixture of ZIKV importation from the Americas and ZIKV infection that has occurred locally in Asia,44 and travel volumes between these continents appear sufficient for international spread to continue for as long as ZIKV is present in either Asia or the Americas.45 Virus evolution may also play a role, with one recent study suggesting that the American strain of ZIKV may be transmitted more efficiently by Asian Ae. aegypti than is the Asian strain of ZIKV.18 There are also questions about whether large-scale climatic forces such as El Niño may have played a role in driving recent patterns of ZIKV transmission on both continents in 2015–2016.46 Definitive answers to these questions are still beyond the reach of current analyses given currently available data,47 but resolving these uncertainties nonetheless remains important for clarifying the future of ZIKV in Asia and for the manner in which health systems should prepare for future ZIKV activity.