Projected Impact of Dengue Vaccination in Yucatán, Mexico

Dengue vaccines will soon provide a new tool for reducing dengue disease, but the effectiveness of widespread vaccination campaigns has not yet been determined. We developed an agent-based dengue model representing movement of and transmission dynamics among people and mosquitoes in Yucatán, Mexico, and simulated various vaccine scenarios to evaluate effectiveness under those conditions. This model includes detailed spatial representation of the Yucatán population, including the location and movement of 1.8 million people between 375,000 households and 100,000 workplaces and schools. Where possible, we designed the model to use data sources with international coverage, to simplify re-parameterization for other regions. The simulation and analysis integrate 35 years of mild and severe case data (including dengue serotype when available), results of a seroprevalence survey, satellite imagery, and climatological, census, and economic data. To fit model parameters that are not directly informed by available data, such as disease reporting rates and dengue transmission parameters, we developed a parameter estimation toolkit called AbcSmc, which we have made publicly available. After fitting the simulation model to dengue case data, we forecasted transmission and assessed the relative effectiveness of several vaccination strategies over a 20 year period. Vaccine efficacy is based on phase III trial results for the Sanofi-Pasteur vaccine, Dengvaxia. We consider routine vaccination of 2, 9, or 16 year-olds, with and without a one-time catch-up campaign to age 30. Because the durability of Dengvaxia is not yet established, we consider hypothetical vaccines that confer either durable or waning immunity, and we evaluate the use of booster doses to counter waning. We find that plausible vaccination scenarios with a durable vaccine reduce annual dengue incidence by as much as 80% within five years. However, if vaccine efficacy wanes after administration, we find that there can be years with larger epidemics than would occur without any vaccination, and that vaccine booster doses are necessary to prevent this outcome.

symptomatic outcomes also affect dynamics in that symptomatic people may withdraw 7581 (stay home) instead of going to work or school. During each day of the symptomatic 7582 period individuals withdraw with probability 0.5. Once individuals have withdrawn, 7583 they remain home until no longer symptomatic. 7584 S1.3: Pathogenicity and Severity Pathogenicity is the probability an infection 7585 produces symptoms. In our model, the reference pathogenicity is for secondary 7586 infections of with DENV1, and this is a fitted parameter. For other infections, we use 7587 relative-risk estimates for serotype-specific factors and number of past exposures to 7588 dengue.

7589
For the serotype-specific factors, we used estimates of reported cases per 1000 7590 infections (medians for fixed-duration cross-protection in Table 2 of supplement in [2]), 7591 and the observed proportions of mild and severe disease by serotype (Table 4 of [3], 7592 using all-years; DF versus sum of all DHF and deaths). For serotype i, the reported 7593 cases, C i , relate to the pathogenicity, ρ i , the mild fraction, f i , and the mild and severe 7594 expansion factors, EF m and EF s by Using these values, we calculated the relative-risk compared to infection by DENV1, 7596  For previous exposure factors, we divide infections into primary infections (no 7599 previous dengue exposure), secondary infections, and post-secondary infections. For 7600 secondary infections, pathogenicity is a fitted parameter. For post-secondary infections, 7601 the relative-risk of pathogenicity is fixed at 0.1. Primary infections, however, depend on 7602 age. 7603 For primary infections, newborns (age 0 individuals) may have maternal antibodies. 7604 Whenever newborns are exposed to an infectious bite, we randomly draw an individual 7605  [2], and observed proportion of mild disease, f i [3], by serotype.
from the females in their household of child-bearing age  years old) to identify a 7606 "mother". If that mother has a past dengue infection, then the infant may have maternal 7607 antibodies and we draw an age in months for the infant uniformly between 0 and 11. 7608 We treat maternal antibodies as cross-immunizing in months 0-2, disease-enhancing 7609 (i.e., like a second infection) in months 3-8, and as having no effect at 9+ months. For 7610 all other ages, we use a step-wise approximation of the curve from [4] for pathogenicity 7611 relative-risk. In our model, age is only a factor for pathogenicity of primary infections, 7612 in contrast to the model described in [5]. 7613 To determine case severity, we fit a parameter for the probability of severe (i.e.

7614
DHF, DSS, or death) disease in secondary cases, fit relative-risk of primary cases 7615 (compared to secondary cases) that is constrained to be less than 1, and assume a fixed 7616 relative-risk (to secondary cases) of 0.2 for post-secondary cases. 7617 S1.4: Dengue Introductions Empirically and in the model, local dengue 7618 transmission generally ceases during the winter and spring, which is relatively dry and 7619 cool (Fig. A, panel D). It is possible for transmission chains to persist, but improbable 7620 because R 0 < 1 from late fall to early spring. To continue to have outbreaks in our 7621 model, we re-introduce dengue from an unmodeled external source (representing, e.g., 7622 travelers that import it from a distinct human population, zoonotic introductions), 7623 similar to what may happen in the real system. This external source may also present 7624 novel serotypes to the population. We model these introductions as random exposures 7625 in the human population, and they may be resisted if an individual is vaccinated or has 7626 been previously infected.

7627
The daily number of people exposed is sampled based on a fitted parameter, λ E : # 7628 exposed per circulating serotype ∼ Pois(λ E ) (see Section S5). This rate is flat over time: 7629 it does not vary by year, nor seasonally over the course of a year. However, since it is 7630 the rate per circulating serotype, the total introduction rate is higher in years with more 7631 circulating serotypes. When an exposure occurs, it is randomly assigned one of the 7632 circulating serotypes. The serotypes being introduced depends on the presence / 7633 absence time series for each serotype (see next section). If there are no serotypes 7634 circulating in a given year, then there are no introductions. 7635 S1.5: Serotype of Introductions Historical case data for Yucatán (Table B) 7636 indicates that serotypes are present (runs) and then absent (gaps) in streaks [6]. Before 7637 each simulation begins, time series are independently constructed for each of the four 7638 serotypes. For the fitting period, we introduce the serotypes in the years in which they 7639 were observed. For historical periods with no recorded data and the forecast period, 7640 time series are generated for each of the four serotypes according to separately fit 7641 geometric distributions (see Section S5). In a given year there may be any combination 7642 of serotypes 1-4 introduced, including no serotypes.  Fig. 4) [7,8]. In that study, serum samples were collected from a random 7646 sample of the population aged from one month to 65 years old who lived in Mérida, 7647 with a particular focus on children (aged 2-10). The adults were sampled from 7648 healthcare facilities in the study area and children were sampled from the general 7649 population and schools.

7650
The age-group-specific prevalence of pre-existing antibody-mediated immunity to 7651 any of the dengue viruses was defined as the proportion of the sampled individuals in 7652 the age group whose serologic specimen is positive to any serotype using Panbio Dengue 7653 IgG Indirect ELISA (titer > 0.10). We computed the exact binomial confidence 7654 intervals (95%) for these prevalence estimates.

7656
Individuals belong to a household, may travel to a school or workplace during the day, 7657 and have an age and gender. These attributes do not change, thus population 7658 demographics are static. The only aspect of the population that is not static is immune 7659 profile.

7660
For simulations of many years of dengue epidemics, which are necessary to establish 7661 long-term dynamics and population-level patterns of immunity, assuming a static 7662 population is obviously imperfect. A demographic model of changing population 7663 structure, formation of new households and shifting spatial distribution (e.g., due to 7664 urbanization) could be more realistic, but it would also be difficult to parameterize and 7665 validate, and to reproduce for other populations of interest in future research. The 7666 current approach avoids the challenges of modeling past and future demography across 7667 varying regions, for which little (or no) data are available. International (IPUMS) data [10]. IPUMS characterizes households by number of 7673 residents, their ages, genders, and employment/student status. Individual households 7674 were sampled from the IPUMS household distribution for each municipality.

7675
Using boundaries from DIVA-GIS [11], households are located by pixel (resolution of 7676 ∼500m) within each municipality, proportional to light output recorded by nighttime 7677 satellite imagery [12] (main text, Fig. 1). Within pixels, households are located 7678 uniformly randomly. The satellite images are composite images taken on nights with no 7679 moonlight. The composite images do not contain clouds, but some background noise is 7680 present, for example from reflected star light. To prevent the placement of households 7681 in regions that are unpopulated, we subtracted from the entire region the average 7682 luminosity of an unpopulated region in western Yucatán, the Petenes-ría Celestún  People are assigned to daytime locations according to the rules below. Rules are 7697 listed in order of precedence.

7698
Everyone with an employment status code that is a student code goes to a school. 7699 Everyone with an employment status code that is a work code goes to a workplace. 7700 Children under 5 years stay home.

7701
Children age 5-11 years go to school.

7704
Students attend the closest school. We use a censored gravity model for employment: 7705 workers select from the 1000 nearest workplaces with probability proportional to 7706 workplace size over Euclidean distance to the workplace squared. These schools or 7707 workplaces are static for individuals, just like their households, and they visit them 7708 every day. The only change to this behavior is if they have a symptomatic infection, in 7709 which case they may stay home (see Section S1), and they resume their daily movement 7710 once their infectious period ends. 7711 S2.3: Immunity & Aging Immunity in the population is established by simulating 7712 100 years of dynamics before assessing baseline dynamics or introducing interventions. 7713 An initial population is created with a modest level of immunity by assuming a 7714 probability of past infection corresponding to a homogeneous 1% annual attack rate per 7715 serotype. We use this approach rather than a fully susceptible initial population in 7716 order to avoid huge early epidemics that are time consuming to model and are not 7717 meaningful, as no data are available from this time period for comparison. The 7718 population is exposed to random dengue introductions for 77 years, followed by a 23 7719 year period with reduced mosquito populations and introductions (see Time), followed 7720 by 35 years of random introductions.  year younger "donates" their immune state. We considered an alternative approach 7731 using a censored gravity model: the donor was selected from the closest 100 individuals 7732 that were one year younger, with probability proportional to the inverse squared 7733 distance between recipient and donor. However, this approach produced no apparent 7734 difference in dynamics and was substantially more computationally demanding.

7735
Newborns (age 0 individuals) have no prior exposure, but may have maternal antibodies 7736 present (see next section for details) contributed from a female in their household. 7737 We "age" the entire population on day 99 of the Julian calendar, which corresponds 7738 to April 9 of non-leap years. This is approximately the middle of the inter-epidemic school or work, while some remain home. This change occurs at fixed times every 7748 simulated day. 7749 We use the female biting behavior in the rainy season in Thailand to estimate 7750 relative preference for biting while people are at home (night) versus their daytime 7751 locations: we assume 8% of biting takes place between sunrise (7 a.m.) and 9 a.m., 76% 7752 takes place between 9 a.m. and 5 p.m, and 16% takes place between 5 p.m. and 7753 sunrise [16]. We define the work/school day as 9 to 5, thus there is a 24% preference for 7754 dusk-night-dawn biting, and 76% for daytime biting.

7755
For each location, we calculate the biting-time-preference weighted total human 7756 population, n, and the preference weighted infectious human population, v, from the 7757 biting preference and time-dependent human populations. From n and v, we calculated 7758 the weighted fraction of bites that expose mosquitoes to infection, f v : We use f v with the effective biting rate, β HM , to determine to determine the number 7760 of newly infected mosquitoes: where m is the number of 7761 susceptible mosquitoes. The susceptible mosquito population at a location is the 7762 location's seasonally adjusted capacity minus the number of tracked mosquitoes (those 7763 that can potentially infect people) currently at that location. This is in contrast to the 7764 previously published model description in which the susceptible mosquito population 7765 was the location's capacity minus the number of infected mosquitoes generated at that 7766 location that are still alive [5]. The infecting serotype is drawn from the serotype 7767 frequencies at that location, also time-weighted by biting preference.

7768
Infected mosquitoes are represented as mobile, individual agents that are generated 7769 in locations with infected humans. Infectious mosquitoes make an infectious bite (which 7770 may be resisted via immunity) with probability β MH each day; this parameter 7771 represents the product of biting frequency and mosquito-to-human transmissibility, and 7772 β MH can be understood as the daily probability of transmission per infectious mosquito 7773 in a fully susceptible human population. Infectious bites occur either during the day or 7774 night, proportional to the population during those times weighted by biting preference: 7775 EIP, or time between a mosquito taking up virus from a host until the mosquito is 7783 able to transmit the virus, is highly sensitive to ambient temperature [1]. Chan and 7784 Johansson (2012) found that the distribution of EIP (in days) for dengue at a constant 7785 temperature (in Celsius) is best fit by a Log-Normal (ln N ) distribution: However, the real environment has a seasonally varying temperature, not a constant 7787 one, so we must integrate to obtain the correct expected EIP as a function of time of 7788 year. We can do so by calculating the incubation rate associated with the expected EIP 7789 (EIP = e µ(T )+ σ 2 2 ) in a given interval δt i with temperature T i : Then the EIP forward from time interval k is We used an input time series of hourly (i.e. ∆t i = 1/24) temperatures for Mérida, 7792 and linearly pro-rated the final interval to obtain EIP k = 1 exactly. This time series 7793 was not complete for the fitting interval (missing ∼8.6% of hourly temperature 7794 readings), so we estimated missing values using temperature data from Miami, USA, 7795 which is climatically similar to Mérida and is more complete (∼0.02% missing). We 7796 related the temperatures between Miami and Mérida using ordinary least squares 7797 regression, which in turn identified high standardized residual (|SR| > 5) in a modest 7798 number (73 out ∼334k non-missing points, or 0.02%) of points in the Mérida time series. 7799 We discarded these points, refit, and used the resulting model to fill in missing 7800 temperatures. We used the hourly temperatures to calculate hourly incubation rates, 7801 which we then integrated over as described above to obtain the EIPs on an hourly basis. 7802 We found the weighted average EIP for each day using the day versus night biting 7803 preferences assumptions (see previous section). In the simulation, this expected EIP is 7804 used to calculate µ for the Log-Normal distribution on a particular day, which we draw 7805 from when infecting mosquitoes.

7806
Mosquito population size also varies seasonally. The number of susceptible 7807 mosquitoes per location is dependent on a daily scaling factor, and a location-specific 7808 baseline capacity.

7809
Daily scaling factors are based on precipitation history (Fig. A, panel B) in Mérida. 7810 Using a climatological dataset from NOAA, we calculated the proportion of years that 7811 each day had precipitation. We then fit a cubic smoothing spline, lagged by 7 days to 7812 accommodate the pre-adult life stages of Ae. aegypti [17]. We normalize the spline to a 7813 maximum of 1.

7814
At the beginning of the simulation, each location is assigned a maximum capacity 7815 sampled from an exponential distribution with a fitted mean, 1/M peak (see Section S5). 7816 We multiply the sampled peak for a given location by the seasonal scaling series to 7817 obtain the daily capacities for that location.  Adjacency among the approximately 480k locations is determined using a filtered 7825 Delaunay triangulation [19]. We preclude long distance movement by removing edges with a haversine length greater than 1 km. As a consequence of traveling along this 7827 network, mosquitoes will only travel short distances in densely populated areas and 7828 longer distance in sparser areas, which is consistent with an empirical study [20].  incidence. We can maintain the correct expected dengue incidence, however, because we 7853 are fitting the number of biting mosquitoes. In this example, the mosquito population 7854 size parameter could increase to maintain a realistic force of infection.

11/38
A B q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  trial results for the Sanofi-Pasteur vaccine in Latin America [22,23], where the overall 7864 VE by serotype is from the intent-to-treat analysis for the trial. In that trial, the 7865 estimated overall VE , ignoring serotype, was twice as high for those who were 7866 antibody-primed for any dengue serotype (i.e., those previously infected), compared to 7867 those who were antibody-naïve (never previously infected). These VE estimates for 7868 antibody-primed and antibody-naïve individuals per serotype are not available from the 7869 vaccine trial results, but we produced the approximate values (see main text, Table 1), 7870 by assuming the antibody positive rate in the vaccine trial participants was 7871 ρ sero+ = 60% when they were vaccinated. Specifically, we used the following relations to 7872 calculate VE by serotype: Vaccinees retain perfect protective immunity against serotypes they were previously 7874 infected with. Among individuals with an infection history ("primed") who are 7875 challenged with a novel serotype, the modeled vaccine is highly effective against

7889
Allowing coverage to vary while keeping the number of vaccinations fixed makes it 7890 easier to attribute the cause of effectiveness differences between strategies. In particular, 7891 it becomes possible to distinguish between the effects of better vaccine efficacy due to 7892 vaccinating older, antibody-positive individuals, versus administering more doses.  study the effects of a vaccine with waning efficacy, we use four different linear waning 7902 models: no waning, and waning with half-lives of 10, 5, and 2 years (Fig. D). When a 7903 vaccinated (but otherwise susceptible) person is bitten by an infectious mosquito,  Quantifying what is meant by saying simulated results are similar to observed data 7925 can be challenging. For some models and observations-e.g. a disease with stable, 7926 predictable annual outbreaks-we might be satisfied reproducing only the mean epidemic 7927 size and similarity would reduce to the absolute error in this result. For reproducing 7928 more complicated observed data, such as matching different types of data (e.g. epidemic 7929 size distribution moments and autocorrelation patterns), how to score similarity is less 7930 obvious. That is the case here, as dengue in Yucatán has complex interannual dynamics 7931 that we would like to be able to reproduce, and we have different types of observations, 7932 e.g. quantiles on total reported cases, annual autocorrelation, severe case trends, and 7933 seroprevalence in a particular population at a particular time. How to get a single 7934 distance between the observed data and simulated results is less obvious for this case. 7935 A simple but flawed approach to matching multiple outcomes, or metrics, is to useful information in others. Partial least squares regression (PLS) addresses these 7943 problems [27,28]. 7944 We do not know a priori the correlation between model parameters and simulated 7945 observed metrics, but we can use PLS to simultaneously orthonormalize both, into a   A more precise description of our ABC-SMC algorithm is the following, based on an 7971 approach described in [26], where N is the number of parameter combinations sampled, 7972 ρ is the fraction retained in the predictive prior/posterior, and n = ρ * N is the size of 7973 the predictive prior/posterior: ii. Simulate data x to sum to 1.

7997
The weight function described in step 6(b)i favors parameter vectors that have 7998 relatively greater density in the prior distribution and lesser density in the most recent 7999 predictive prior, thus offsetting a sampling bias that would otherwise be introduced.

8000
The algorithm above differs from that in Beaumont (2010) in that the best 8001 n = ρ * N parameter combinations is retained, rather than an absolute number that 8002 score better than some threshold that decreases with each SMC step. Beaumont's 8003 approach involves sampling parameter space an indefinite number of times until a fixed 8004 number of samples have resulted in a distance less than . Because we necessarily use an 8005 updated PLS model for each SMC step, the meaning of the numerical value of changes 8006 from step to step. Indeed, the dimensionality of the space that would be applied to is 8007 changing. It is also impossible to confidently estimate the running time required for 8008 Beaumont's algorithm, which is a common requirement in supercomputing 8009 environments that use a job scheduler. If the selected threshold is too strict, it is 8010 possible that no parameter vectors would ever be selected, and the algorithm would run 8011 indefinitely. Specifying but setting an arbitrary time limit could result in an undefined 8012 predictive prior.

8013
For fitting of the serotype introduction patterns (see below), we use N = 10 6 and 8014 select the best performing 1% of parameter combinations. More thorough sampling is 8015 possible than with the epidemic model parameters, because simulation of the serotype 8016 patterns is orders of magnitude faster. All four dengue serotypes have been observed in Yucatán, but they are present with 8038 varying runs and gaps: for example, DENV1 was observed from 1987 to 1997 (an 11 8039 year run), and then was not observed from 1998 to 2001 (a 4 year gap). A plausible 8040 biological explanation is that an introduced serotype is able to persist regionally until a 8041 threshold of herd immunity is reached, at which point dengue levels may decrease until 8042 another serotype is introduced or enough new susceptible people enter the population, 8043 via birth or immigration, for the first serotype to circulate again. Consistent with that 8044 interpretation-i.e. gaps are driven by herd immunity-we assume that the initial delay 8045 in appearance of a serotype is not a gap. 8046 We generate random sequences of alternating runs and gaps for each serotype by separate distributions: for runs and gaps, and for each of the four serotypes. We fit the 8050 distribution parameters p run , the probability that a run will end in a given year, and 8051 p gap , the probability that a gap will end, separately for each serotype using ABC-SMC. 8052 Because the observed data are incomplete (right censored after 2013, missing for Fitting these parameters at the same time as the others might improve the dengue 8070 introduction match to observations, since there is an interplay between introduced 8071 serotypes and which are observed at what time, as well as inter-serotype dynamics.

8072
However, it is not practical because of how much it would increase the complexity of 8073 parameter space, and how long it takes to run the epidemic simulation once.

8074
Furthermore, the data available from the fitting period may not be particularly 8075 applicable to the past and forecast introduction series, given that these may be largely 8076 driven by external phenomena that are insufficiently reflected in the fitting period, so a 8077 more detailed approach is unlikely to improve forecasts.

8078
The actual observed, simulated observed, and distribution means are reported in E. 8079  Table E. Mean duration (years) of presence (runs) and absence (gaps) for each serotype. Simulated observed means were generated using the distribution means and then censored as the observed data were censored.

S5.4: Epidemic Model Parameter Estimation
The remaining model parameters 8080 were estimated using the ABC-SMC procedure described above to identify parameter 8081 combinations that produce dynamics similar to the observed data. For "similarity", we 8082 consider the following metrics: Priors for these parameters are described in Table D. 8116 We use AbcSmc to fit these parameters. We ran AbcSmc for a total of 14 SMC sets, 8117 with a sample size of 10,000 parameter combinations for each set. To reduce the previously (see Section S2). The last 7 SMC sets used the full Yucatán population. The 8124 metrics we use are calculated on a per capita basis, so the interpretation of the metrics 8125 is insensitive to the size of the modeled population, but the size and spatial structure of 8126 the population (e.g., one large city versus that city and the surrounding towns and rural 8127 areas) may affect epidemic dynamics. Using a smaller population initially allowed us to 8128 focus in on promising regions of parameter space with faster-running simulations, which 8129 we then explored more thoroughly with the full population. 8130 We observed during the Yucatán SMC sets that the predictive priors were 8131 fluctuating noisily rather than converging: see, for example, the introduction rate λ E in 8132 Fig. E. This is may be due to inadequate sample size; however, sampling substantially 8133 more than 10,000 parameter combinations per set would be prohibitive, as each The last column, All-Y is not a normal SMC set, but instead represents all Yucatán parameter combinations simulated in Y1-7 (dark grey) and the global 100 best performing combinations (light grey).  year [30]. We start by inverting the previously calculated daily EIPs into T eff (day) (the 8153 constant temperature corresponding with that day's mean EIP), then create a 8154 forecasting time series of T eff ((day)) + 0.02 * y, where y is the forecast year. We then 8155 recalculate EIPs using the previously described approach, but with the daily T eff s 8156 instead.

8157
In general, we expect dengue incidence to increase with warmer climate. Because we 8158 assume that EIP decreases with increasing temperature, the likelihood that an infected 8159 mosquito will spread dengue before dying increases in turn. Indeed, when long-term 8160 warming is included in the model, we see that dengue incidence increases in both 8161 baseline and intervention scenarios (Fig. G). However, because the increases in incidence 8162 under climate change are approximately proportional to incidence without climate 8163 change, vaccine effectiveness remains essentially unaffected.

8164
While we believe the effectiveness results to be relatively robust to alternate climate 8165 forecasts, we emphasize that the incidence forecasts will be more sensitive. For example, 8166 since temperature has a non-linear effect on EIP, warmer nighttime lows may have a life-cycle factors influenced by temperature [31,32]. Similarly, we did not include any 8171 changes in precipitation. An international consensus model of climate change indicates 8172 the region will become drier [33], which may have a moderating effect on any increase in 8173 force of infection.  Cumulative V eff is computed similarly, but from the start of the intervention up the 8190 relevant year, rather than annually.

8191
Annual V eff shows short-term vaccination scenario performance over moving one-year 8192 windows of an intervention, which corresponds to thinking about annual snapshots of 8193 the program. This is the during-an-outbreak view: for example, what the general 8194 populace or officials in a region undertaking a vaccine program might consider their 8195 outcomes compared to neighboring regions without vaccination. For effective programs, 8196 they will see lower disease burden in most years. However, in some particular year, an 8197 outbreak might occur in the region with vaccination and not in a non-vaccinating 8198 neighboring area, simply due to the stochastic nature of introduction of the pathogen, 8199 disease spread, and (particularly for dengue) manifestation of disease after infection. In 8200 the long view, for an effective intervention, these relatively bad years should be fewer 8201 and outweighed by the benefit in most years, so it is also useful to consider cumulative 8202 V eff .

8203
When aggregating over many simulated samples of these perspective, robustness is 8204 an important consideration in choice of statistic (e.g. mean versus median effectiveness), 8205 since the distribution of V eff can be highly skewed: it has an upper bound of 1 (at most, 8206 all cases can be avoided), is typically between 0 and 1 (i.e. the intervention is effective), 8207 but can go to negative infinity (i.e. when comparing any number of cases in an 8208 intervention to 0 cases in the baseline). It can also produce large negative V eff by chance 8209 when case counts are generally low, or (in the case of annual effectiveness) when an 8210 intervention postpones-but does not prevent-an epidemic from one year into another. 8211 Even though the resulting epidemic may be smaller, it may be compared to a small 8212 number of cases from the baseline where many people have cross-immunity from the 8213 previous year's large epidemic. When these delays are misaligned between many 8214 samples, large negative mean V eff appears in many years because the individual large 8215 negative V eff observations appear at many different times across different runs.

8216
The net result of these concerns is that mean V eff is not a robust statistic. If we 8217 calculated V eff on mean cases (instead of using the cases from each sample and then 8218 taking the mean of V eff ), this is equivalent to calculating a weighted mean V eff where 8219 the weights are the relative number of cases in the baseline. This ignores the magnitude 8220 of the large negative V eff values associated with low baseline case counts, and likely 8221 corresponds with an intuitive perception of those events in many cases (e.g. the 8222 non-vaccinating, neighboring region reporting no cases when we reported five is not a 8223 huge public health failure). However, it does not properly account for large, negative 8224 annual V eff values associated with delayed epidemics in an intervention scenario when 8225 the baseline merely had its large epidemic earlier.

8226
Median and quantile V eff values, however, are robust, because they are not sensitive 8227 to the skewness of the V eff distribution. The sensitivity to negative V eff also tends to be 8228 eliminated when considering cumulative V eff s, since the synchronization of annual 8229 epidemics between scenarios becomes irrelevant: only total relative case counts for the 8230 considered years matters. Similarly, small outbreak years (and their potential for 8231 unimportant, yet large negative V eff ) are aggregated with large epidemic years, thus 8232 tending to make even net negative V eff small unless the intervention actually harms the 8233 population. This corresponds to our intuition about how good (or bad) an intervention 8234 is. This robustness is also why forecasted cumulative V eff has decreasing uncertainty: 8235 including more years means more time for the trend to emerge above inter-annual 8236 variation. intervention and with routine, durable vaccination, across target ages with and without 8239 catchup campaigns), then for annual and cumulative V eff for the routine vaccine, then for 8240 annual and cumulative V eff for a waning vaccine with and without a booster program. 8241 Figure H. Annual total cases per capita 95% (lightest) and 68% (darker) prediction intervals, and medians (solid line) for interventions with a durable vaccine. The panels are arranged by routine age of vaccination as well as non-intervention (rows) and routine versus routine-and-catchup campaigns (columns). Color is replicated from main text figures.
Figure I. Cumulative total cases per capita 95% (lightest) and 68% (darker) prediction intervals, and medians (solid line) for interventions with a durable vaccine. The panels are arranged by routine age of vaccination as well as non-intervention (rows) and routine versus routine-and-catchup campaigns (columns). Color is replicated from main text figures.
Figure J. Annual VE 95% (lightest) and 68% (darker) prediction intervals, and medians (solid line) for interventions with a durable vaccine. The panels are faceted by routine age of vaccination (rows) and routine versus routine and catchup campaigns (columns). Color is replicated from main text figures. Note that the lower forecast bound (2.5%) for the annual view is always negative, though it trends slowly upward (left out of frame), similar to the low bound (16%) for routine only strategies. For strategies with catchup campaigns, however, the resulting annual VE becomes less certain as the large initial coverage of seropositive vaccinees declines due to mortality.

29/38
Figure K. Cumulative VE 95% (lightest) and 68% (darker) prediction intervals, and medians (solid line) for interventions with a durable vaccine. The panels are faceted by routine age of vaccination (rows) and routine versus routine and catchup campaigns (columns). Color is replicated from main text figures.
Figure L. Annual VE 95% (lightest) and 68% (darker) prediction intervals, and medians (solid line) for interventions with a waning vaccine and no booster vaccination campaigns. Routine age of vaccination is 9 years old. The panels are faceted by routine versus routine and catchup campaigns (rows) and waning half-life (columns). Color is replicated from main text figures. As noted in main text, a waning vaccine leads to limited VE, and when combined with a catchup campaign, leads to an expected temporary negative VE, which is not seen in any other scenario.