Modeling the spread of polio in an IPV-vaccinated population: lessons learned from the 2013 silent outbreak in southern Israel

Background Polio eradication is an extraordinary globally coordinated health program in terms of its magnitude and reach, leading to the elimination of wild poliovirus (WPV) in most parts of the world. In 2013, a silent outbreak of WPV was detected in Israel, a country using an inactivated polio vaccine (IPV) exclusively since 2005. The outbreak was detected using environmental surveillance (ES) of sewage reservoirs. Stool surveys indicated the outbreak to be restricted mainly to children under the age of 10 in the Bedouin population of southern Israel. In order to curtail the outbreak, a nationwide vaccination campaign using oral polio vaccine (OPV) was conducted, targeting all children under 10. Methods A transmission model, fitted to the results of the stool surveys, with additional conditions set by the ES measurements, was used to evaluate the prevalence of WPV in Bedouin children and the effectiveness of the vaccination campaign. Employing the parameter estimates of the model fitting, the model was used to investigate the effect of alternative timings, coverages and dosages of the OPV campaign on the outcome of the outbreak. Results The mean estimate for the mean reproductive number was 1.77 (95 % credible interval, 1.46–2.30). With seasonal variation, the reproductive number maximum range was between zero and six. The mean estimate for the mean infectious periods was 16.8 (8.6–24.9) days. The modeling indicates the OPV campaign was effective in curtailing the outbreak. The mean estimate for the attack rate in Bedouin children under 10 at the end of 2014 was 42 % (22–65 %), whereas without the campaign the mean projected attack rate was 57 % (35–74 %). The campaign also likely shortened the duration of the outbreak by a mean estimate of 309 (2–846) days. A faster initiation of the OPV campaign could have reduced the incidence of WPV even if a lower coverage was reached, at the risk of prolonging the outbreak. Conclusions OPV campaigns are essential for interrupting WPV transmission, even in a developed country setting with a high coverage of IPV. In this setting, establishing ES of WPV circulation is particularly crucial for early detection and containment of an outbreak. Electronic supplementary material The online version of this article (doi:10.1186/s12916-016-0637-z) contains supplementary material, which is available to authorized users.


Transmission model
Following is the formulation of the deterministic, discrete-time SEIR transmission model that is described schematically in figure 1 of the main text: where d I is the mean infectious period, R ̅ is the mean reproductive number and ( ) is a periodic function with the property ( ( )) = 0 describing the seasonal variation in the reproductive number, which is modeled using the estimated seasonality in the southern US states during the pre-vaccination era (see 'Modeling Seasonality' section below).
Exposed individuals become infectious at rate 1/ ( being the mean period of latency) and infectious individuals recover at rate 1/ . 1 ( ) and 2 ( ) denote the number of individuals vaccinated with a first and second dose of OPV on day t. We used the data provided by the Israeli Ministry of Health, with a delay of one week, as this is the time it takes for the vaccines to build up protective immunity [1]. denotes the per-dose efficacy of OPV (i.e., the probability that the vaccine will build up protective intestinal immunity, as in [2][3][4]]). Finally, ℎ 1 ( ) and ℎ 2 ( ) denote the probability that a first and second dose of OPV given on day t, will be given to susceptible individual. These probabilities are calculated according to the fraction of susceptible individuals in the target population for each dose (since the first dose was only given to those who did not receive one yet and the second dose was given only to individuals who received a first dose): The initial conditions of the model for the start of the epidemic are

Modeling seasonality
In order to model the seasonal variation in the reproductive number we employ information available from a study modeling the spread of WPV in the US in the 1930s-1950s [5]. The study estimated the seasonality of WPV in the different US states. We used the estimated seasonality for ten southern US states (Alabama, Arizona, California, Florida, Georgia, Louisiana, Mississippi, New Mexico, South Carolina and Texas) that share the same latitudes with Israel (30N-33N), as the study indicated that the timing of seasonality of WPV has a latitudinal gradient. Using these data we modeled the seasonal variation ω( ) (appearing in eq. S2) as: where ω ̅(t), 1 ≤ t ≤ 365, is a function of the mean seasonality of the ten US states obtained by interpolating the monthly estimates into daily data, centering the seasonal functions so that their peaks align at the day of the mean peak time, taking the mean of the functions and normalizing it (dividing by the mean and subtracting by one) so that (ω ̅(t)) = 0 ( fig. S3a-d).
The parameters δ and Φ are model parameters that allow to control the amplitude and peak time of the seasonal variation. δ sets the amplitude so that δ = 0 means there is no seasonal variation while δ = 1 means a seasonal variation similar to the mean seasonal variation of ω ̅. ϕ  (a) Estimates of the monthly variation in the reproductive number of WPV in ten southern US states obtained by fitting epidemic data from the 1930s-1950s (data taken from [5]).
(b) The curves from (a) after interpolation into daily data.
(c) The curves from (b) after centering them so that their peaks align at the mean peak time.
The thick black curve is the mean of the ten seasonal curves.
(d) The curves from (c) after normalizing each of the curves so that the mean of each curve is zero. The thick black curve is the normalized mean seasonality (ω ̅ in eq. S4).

Observation model
We used a binomial observation model to link the prevalence given by the transmission model to the stool survey data: where (  can be written as: In addition, as mentioned in the main text, we employed the information available from ES by

Bayesian inference
We adopted a non-informative uniform prior for the for the duration of infectiousness ( ) in the range [7,49] days based on estimates ranging from 4 to 7 weeks in naïve populations [6][7][8], and the suggestions that the period of infectiousness can be shorter in a population vaccinated with IPV [9][10][11]. For the mean reproductive number ( ̅ ) we set a uniform prior in the range [1,10], as 1 is the threshold value above which an epidemic can occur, and 10 is the high end estimate for ̅ , attributed to a situation in a country with poor sanitary conditions [12]. For the two seasonality parameters (δ and ϕ) we employed prior distributions obtained from estimates of the seasonal variation of WPV in southern US states, as described above. The prior distribution for the per-dose efficacy of OPV was set as (0.56,0.23) according to estimates given in a recent meta-analysis study [13]. We set a uniform prior for the initiation time of the epidemic ( 0 ) in the range [1,145] (which implies initiation dates between September 15, 2012 and February 6, 2013). This lower limit of this range was set based on the results of a phylogenetic study, analyzing poliovirus isolates from Israel and Egypt, that indicated the Israeli and Egyptian WPV1 lineages have diverged around mid-September 2012 [14]. The upper limit for this range was set according to the timing of the first identification of WPV1 in Israel using ES [15].
We explored the posterior distribution of the model parameters by employing Markov Chain Monte Carlo (MCMC) sampling performed using the slice sampling method [16] (implemented using the slicesample function in matlab) on the sum of log ( ) given in eq. S6 and the log of the prior distributions, with the additional limitation given in eq. S7. We discarded the first 10,000 iterations as a burn-in period and used a thinning of 1 in 10 samples as to avoid autocorrelation between adjacent samples. The remaining 10,000 samples after the burn-in period and the thinning were used to obtain the posterior distribution for the model parameters. 95% credible intervals were obtained by removing the 5% of the samples with the lowest likelihood and computing the range within the remaining samples. Figure S4 shows the traces of the parameter values obtained by the MCMC. Figure S5 shows  . S5). We also found a strong negative correlation between and the overall attack rate ( fig. S6). Since we fit prevalence data, a longer period of infection is able to help fit the same data with fewer cases, thus explaining the negative correlation between d I and the attack rate. For this reason we also obtain smaller reproductive number values for larger values of d I . A larger necessitates a larger ̅ to obtain the same number of cases, which explains the positive correlation between these two parameters.

Results of model with no seasonal variation in transmission
In order to compare the results of the model with seasonal variation in transmission to the results of a model without seasonality, we ran the MCMC procedure for a second time while fixing the seasonality parameters so that = 0 (i.e., no seasonal variation) and = ̅ (i.e., the mean of the peak time in (eq. S4)). While the latter has no effect on the results of the model fitting (with = 0 the value of does not affect the transmission rates), we did so in order to obtain the maximum value of the prior distribution for , so as to allow a comparison of the likelihoods obtained with this model and the model with seasonality (see table S2 below). Table   S1 and figure S8 show the results obtained using this model. The mean and 95% CI estimates obtained using this model for the period of infectiousness ( ) and for the reproductive number  (table 2 and table S1).      Table S2: A comparison of the two models, with and without seasonality, using deviance information criterion (DIC) [17]. ̅ was calculated as the average of the sampled values obtained using the MCMC. ̅ was calculated using the mean of the log-likelihoods of all the sampled values and 2 was calculated using the variance of the log-likelihoods of all the sampled values.
The model with seasonality, which incorporates two additional parameters compared to the non-seasonal model, manages to obtain slightly better likelihoods over all, as can be seen by the difference in the value of ̅ between the two models. However, the calculated effective number of parameters (either 1 or 2 ) is higher for this model so that the obtained DIC score (either 1 or 2 ) are only slightly better for the seasonal model, and not enough to rule out the non-seasonal model. Table S3 and