Type 2 poliovirus detection after global withdrawal of trivalent oral vaccine

Background Mass campaigns with oral poliovirus vaccine (OPV) have brought the world close to wild-type poliovirus eradication. However, to complete eradication, OPV must itself be withdrawn to prevent vaccine-derived poliovirus outbreaks. Synchronized global withdrawal of OPV began with serotype-2 (OPV2) in April 2016, presenting the first test of the feasibility of eradicating all polioviruses. Methods We analysed global surveillance data reporting serotype-2 vaccine poliovirus (Sabin-2) and vaccine- derived poliovirus (VDPV2) detection in stool collected during 1 January 2013 through 08 August 2017 from 431,429 children with acute flaccid paralysis in 112 countries and 5485 environmental samples from 4 high-risk countries. We used Bayesian spatiotemporal smoothing and logistic regression to identify and map risk-factors for persistent Sabin-2 and VDPV2 detection. Results Detection of Sabin-2 declined rapidly from 3.9% in stool [95% Confidence Interval: 3.5%-4.3%] and 71% [61%-80%] in sewage at the time of OPV2 withdrawal to 0.16% [0.09%-2.7%] and 13% [8%- 20%] at 2 months. However, at 12 months Sabin-2 continued to be detected (0.05% [0.01%-0.13%] in stool, 8% [5%-13%] in sewage) due to OPV2 use in response to VDPV2 outbreaks. Five outbreaks were reported after OPV2 withdrawal, associated with low routine immunisation coverage and population immunity (Odds Ratios 2.59 [1.26-6.33] and 4.65 [1.71-15.28] per 10% absolute decrease). Discussion High population immunity has facilitated rapid decline of Sabin-2 after OPV2 withdrawal and restricted circulation of VDPV2 to known areas at high risk of transmission. It will be critical to control these remaining VDPV2 outbreaks before the growth of significant cohorts of susceptible children.


Supplementary methods:
Surveillance data Children aged 0-14 years with acute flaccid paralysis (AFP) are reported in each country through a network of healthcare providers as part of routine polio surveillance. 1 Information collected from each child includes the dates of onset and notification of paralysis, the date of stool collection, and the province and district of residence. Two stool samples collected 24 hours apart within 14 days of onset of AFP are recommended to be collected from each AFP case and are analysed to determine the presence of wild, vaccine-derived or Sabin polioviruses using Global Polio Laboratory Network protocols 2, 3 . The vaccination history of these children and detection of vaccine polioviruses in their stool samples can therefore provide information on vaccine coverage and vaccine poliovirus circulation in the general population 4 . Here we analysed AFP data from countries in the African,

Eastern Mediterranean, Southeast Asian and European regions obtained through the Polio Information
System maintained by the GPEI. The Polio Information System also records vaccination campaign timing, location and vaccine used (trivalent, bivalent or monovalent OPVs, IPV). Data from AFP cases whose stool samples were collected between 1 January 2013 and 08 August were analysed.
Environmental surveillance (ES), the systematic collection of sewage samples from a given site to test for the presence of polioviruses, is currently performed in over 30 countries to supplement AFP surveillance. 5,6 Here we analyse ES data from four high-risk countries in the Eastern-Mediterranean Region (Afghanistan, Pakistan) and Africa (Kenya, Nigeria) collected between 1 January 2013 and 08 August 2017. The number and spatial distribution of collection sites is shown in Fig S1. Most sites have samples collected monthly. ES samples were collected through the 'grab' method after which samples are concentrated using the two-phase separation method 5,7,8 , while additional samples in Pakistan were also collected through the 'Bag Mediated Filtration System' 9 . Samples were then tested for the presence of polioviruses (wild, VDPV and Sabin viruses) using standard WHO protocols 8 .
Geographical information system data for boundaries of first and second administrative areas were obtained from the World Health Organization. We refer to first administrative areas as provinces and second administrative areas as districts.

Description of statistical analysis
Statistical model of Sabin-2 detection i) Detection from non-polio AFP The probability of Sabin-2 detection from non-polio AFP as a function of time since the last tOPV campaign was estimated through fitting a logistic regression model with a (logit link) offset to the proportion of non-polio AFP cases positive for Sabin-2 in stool, using data from before OPV2 withdrawal and conditioning on cases with stool collected within six months of a campaign. The model was fitted to non-polio AFP data with stool samples collected between 1 January 2013 and 15 April 2016. If Ni,t is the total number non-polio AFP in population i, t days since the last tOPV campaign in the respective province and Xi,t is the number of observed non-polio AFP with Sabin-2 detected in stool then ,~( , , , ) . The parameter , is the probability of detecting Sabin-2 and is defined as , = + (1 − ) , where ( , (1 − , ) ⁄ ) = + . The parameter qi allows the probability of Sabin-2 detection to asymptotically decline to a non-zero level following a tOPV campaign to account for the fact that routine immunisation with tOPV leads to a low background prevalence of Sabin-2 detection. The background level also accounts for detection of  Table S1) to allow estimating a separate value of qi per population. A summary of model parameters is given in Table S3. The fit of the model was assessed through plotting the observed data against the posterior predictive distribution, aggregating at the monthly level ( Figure 1). In addition, the distribution of the binned residuals (observedposterior predictive distribution) by month was plotted ( Figure S9), as described by 11 . No trends in the binned residuals occurred.
We tested the correlation across populations between the estimated time (days) needed for the prevalence of Sabin-2 to reach background levels (+ 0.1%) following a tOPV campaign and serotype-2 population immunity. Serotype-2 population immunity was estimated at subnational levels over 6-month time periods from vaccination histories of non-polio AFP cases (<36 months old) and serotypespecific estimates of OPV efficacy against serotype-2 poliomyelitis using a spatiotemporal randomeffects model as previously described 10 . Each non-polio AFP case in each population was assigned the immunity estimate for six-month period the stool sample was collected in and the average population immunity was calculated as the median estimate across AFP cases in each population (where populations are defined as in Table S1). To assess general trends in immunity over time across all considered countries the population weighted mean and standard deviation were calculated using population size data from Worldpop 12 .
The statistical model of Sabin-2 detection was then used to predict the expected probability (στ,t,i) of Sabin-2 detection following OPV2 withdrawal: , , = , + (1 − , ) , where , to allow for a decline in Sabin-2 detection from the final children immunised with tOPV through routine immunisation (compared to a constant value above), with τ being the number of days since OPV2 withdrawal (1 May 2016), and t the number of days since the last tOPV or mOPV2 campaign. The model assumes that the rate in the decline of Sabin-2 from the final children immunised with tOPV through routine immunisation occurs at the same rate as the decline following tOPV campaigns (hence is multiplied by in the linear predictor of , ). The model also assumes that the change in prevalence of Sabin-2 following an outbreak response mOPV2 campaign after OPV2 withdrawal is similar to that following a tOPV campaign prior to OPV2 withdrawal. A summary of additional model parameters is given in Table S4. Latin hypercube sampling 13 was performed to efficiently sample , , , and from their posterior marginal distributions as it was not possible to jointly sample these within the R-INLA framework.
The predictions were performed from 1000 samples and the 2.5% and 97.5% quantiles of the simulations were compared to the observed proportion of non-polio AFP cases with Sabin-2 detection following the OPV2 withdrawal. Unexpected Sabin-2 detections were defined when the observations were outside the 95% prediction interval.
ii) Detection from ES samples A similar statistical model was fitted to the ES data whereby , , represents the number of ES samples with Sabin-2 detected from area i, t days since the last tOPV campaign from site j, and , , the number of samples tested. The term area refers to either Afghanistan, Kenya, North Nigeria, South Nigeria, Punjab and Islamabad provinces in Pakistan or the remainder of Pakistan. In addition, a persite random effect was included in this model to account for site-specific variation. Therefore where was assumed to be Normally distributed with mean 0 and estimated variance σ 2 . Samples that were scheduled but not collected were not included in the analysis. As described above, the model assumes independence between provinces. The fit of the model was assessed through plotting the observed data against the posterior predictive distribution, aggregating at the monthly level ( Figure 1). In addition, the distribution of the residuals (observedposterior predictive distribution) by month was plotted ( Figure S10). When comparing the observed decline of Sabin-2 to the predicted decline, ES sites that began sampling after OPV2 withdrawal were not included (Fig 2).

Risk factors of cVDPV2 cases post OPV2 withdrawal
The emergence and spread of cVDPV2 has been previously found to occur in populations with poor routine immunisation coverage, low population immunity and high birth rates 14 . In this work we explore whether provinces that report at least one cVDPV2 case after OPV2 withdrawal occurred in places that would be considered high risk (the number of cases by province is given in Table S5).
Provinces, , , in each country, j, were classified as either reporting or not reporting cVDPV2 cases post OPV2 withdrawal and the log odds of reporting cases was assumed to be a function of previously The method to estimate serotype-2 population immunity against poliomyelitis is referred to above.
The province estimates for the first half of 2016 were used in the regression to reflect immunity that would have allowed VDPV2 that emerged pre OPV2 withdrawal to circulate into the post-withdrawal time period. Province estimates for Nigeria and Pakistan were obtained by taking the mean across relevant districts.
Sub-national estimates of the number of annual births were not available for all the considered countries (i.e. Syrian Arab Republic) and therefore we use population size as a proxy measure for this risk factor. For each country the 2015 estimated population size from spatially smoothed surfaces 12,18 were aggregated the province level using the R 'velox' package. 17 Note that the population estimates for Syria are based upon census data collected prior to the recent conflict in this country.
Univariable analyses were performed followed by multivariable model selection using the 'widely applicable information criteria' (WAIC) 19 . All analyses were performed using the R-INLA package. 20 Supplementary Results

Risk Factors for cVDPV2 cases reported after OPV2 withdrawal
A map of the routine immunisation coverage per province is given in Figure S11. Routine immunisation coverage was estimated to be variable both between and within countries. The areas with the lowest routine immunisation coverage were northern Syria, Western Pakistan, central DRC and northern Nigeria. The total population size also was variable across provinces, with the highest population located in eastern Pakistan ( Figure S12).
Multivariable model selection was performed to determine the most parsimonious yet adequate model (Table S6). The model which included routine immunisation coverage, population immunity and population size resulted in the lowest WAIC. The estimated variance for the country level random effect in the final model was 0.000 (95% credible interval: 0.000 -0.002).