Population sensitivity of acute flaccid paralysis and environmental surveillance for serotype 1 poliovirus in Pakistan: an observational study

Background To support poliomyelitis eradication in Pakistan, environmental surveillance (ES) of wastewater has been expanded alongside surveillance for acute flaccid paralysis (AFP). ES is a relatively new method of surveillance, and the population sensitivity of detecting poliovirus within endemic settings requires estimation. Methods Data for wild serotype 1 poliovirus from AFP and ES from January 2011 to September 2015 from 14 districts in Pakistan were analysed using a multi-state model framework. This framework was used to estimate the sensitivity of poliovirus detection from each surveillance source and parameters such as the duration of infection within a community. Results The location and timing of poliomyelitis cases showed spatial and temporal variability. The sensitivity of AFP surveillance to detect serotype 1 poliovirus infection in a district and its neighbours per month was on average 30.0% (95% CI 24.8–35.8) and increased with the incidence of poliomyelitis cases. The average population sensitivity of a single environmental sample was 59.4% (95% CI 55.4–63.0), with significant variation in site-specific estimates (median varied from 33.3–79.2%). The combined population sensitivity of environmental and AFP surveillance in a given month was on average 98.1% (95% CI 97.2–98.7), assuming four samples per month for each site. Conclusions ES can be a highly sensitive supplement to AFP surveillance in areas with converging sewage systems. As ES for poliovirus is expanded, it will be important to identify factors associated with variation in site sensitivity, leading to improved site selection and surveillance system performance. Electronic supplementary material The online version of this article (10.1186/s12879-018-3070-4) contains supplementary material, which is available to authorized users.


S1. Description of the multistate modelling framework and parameter estimation
Within each district we assumed the population was either free from poliovirus circulation (uninfected) or infected with poliovirus. These assumptions correspond to a two-state Markov process with constant hazard rates λ and γ of infection and recovery, respectively. A Markov process is one in which the future state of the system depends only on the present state of the system and not on the total history. Using the steps outlined in Gentleman et al. [1], the transition probability matrix can be derived for this and more complex systems, and is essentially a function of the hazards λ and γ. A transition probability matrix P (t) is used to represent the probability of moving from one state to another, where p denotes the transition intensity of moving from uninfected to infected and 1 − q refers to moving from infected to uninfected. For a two-stage Markov model the transition intensities can be easily calculated by solving the differential equations that describe the system; The state of the system within district i(i = 1, ..., M ) at time t(t = 1, ..., T ) is denoted by z it , where a value z it = 0 denotes an uninfected district and a value z it = 1 denotes an infected district. The initial conditions for z (ie. z it=0 ) are estimated. It is important to note that z it is an augmented (latent) parameter, as the true infection status is not directly observed but only estimated through the surveillance data. The AFP data are denoted by x it where if x it = 0 then no poliomyelitis cases were reported in the district and if x it = 1 then at least 1 poliomyelitis case was reported within the time period and the sensitivity of one AFP month is denoted by α, in other words the probability of a positive result in an infected district, p(x it = 1|z it = 1), is equal to α. The environmental data consist of y sit positive samples from s sites within each district, from a total of n sit samples tested that month, with sensitivity ω per sample. In the first instance we assume that ω does not vary and consequently y it = S s=1 y sit and n it = S s=1 n sit . In this way the Markov model utilizes data on the infection status whilst allowing for error in reporting due to suboptimal sensitivity of each surveillance system. Within a Bayesian framework the posterior density of the parameters are a function of the likelihood of the model and the priors of the parameters and augmented data; The model parameters were estimated using Markov chain Monte Carlo (MCMC) methods [2]. Where the surveillance sensitivity parameters (α and ω) take one value they were estimated using a Gibbs sampler where priors were assumed to be beta distributed with parameters α 1 = 2 and β 1 = 2 in the case of α and α 2 = 2 and β 2 = 2 in the case of ω. The augmented data z were updated using a Gibbs sampler based on the current values of the transition intensities, the data and values of α and ω. The posterior estimate of z i was used in estimation of the false omission rate. The parameters λ and γ were estimated using the Metropolis-Hastings step and the priors for λ and γ were log-normally distributed with mean (λ m = −2 and γ m = −2) and variance (λ v = 0.1 and γ v = 0.1).
The model described above was extended the test hypotheses of variation in the sensitivity of each surveillance system. District variation in sensitivity of AFP surveillance was included by allowing α to vary by district, first assuming independence of each district and sampling estimates via a Gibbs sampler. A linear relationship between the sensitivity of α i and the incidence of poliomyelitis k i (on a log 10 scale) was modelled assuming log α i 1−α i = g + hk i where g and h were estimated using the Metropolis-Hastings step. For environmental surveillance a similar approach was taken, first assuming independence and then exploring additional functional forms. A mixed effects model was assumed where the sensitivity at each site took the form log ω si 1−ω si = β + b i + si where β was the average sensitivity and was modelled assuming β ∼ Normal(µ β , σ 2 β ), b i was the district random effects which were modelled assuming b i ∼ Normal(0, σ 2 b ), and si were the site random effects which were modelled assuming si ∼ Normal(0, σ 2 ). Each parameter (β, σ 2 b and σ 2 ) was estimated using a Gibbs sampler where the priors on the variances were inverse gamma distributed. In another model sensitivity was assumed to be a function of catchment size m si (on a log 10 scale). A cubic relationship was assumed between catchment size and sensitivity to allow for a non-linear relationship (ie. log ω si 1−ω si = dm 3 si + cm 2 si + bm si + a), where a, b, c and d were estimated, and simpler models were tested by assuming each parameter was equal to zero (ie. for a quadratic model d = 0). The fit of these alternative models to the data were compared using Bayes factors.

S2. Estimation of catchment area size from digital elevation modelling
During wastewater sample collection, GPS coordinates are taken of the sampling location, where in subsequent visits the same location is sampled. The catchment area of sampling locations has been estimated using a digital elevation modelling approach, where the resolution of the topography is 30x30m [3]. The topography of the landscape and river flows are used to identify land areas where wastewater is likely to flow towards the environmental sampling site. The estimated polygon is overlaid onto resolute population size estimates of Pakistan, available at a 100x100m resolution [4]. The estimates of catchment areas do not include any information regarding sewer networks in Pakistan, largely because the information is unavailable. A summary of the estimated catchment sizes (as of March 2016) are shown in S1 Table. S3. Estimation of Bayes factors using thermodynamic integration Bayes factors were estimated using thermodynamic integration (TI) methods, which enabled comparison of the different multistate models [5]. These methods lend themselves well to estimation of model fit in a framework with augmented data such as multistate models, because they are not reliant upon knowing the number of parameters within a model (unlike the better-known deviance information criteria). Given a model M , with parameter vector θ and applied on dataset D, the posterior probability distribution is given by Bayes Theorem; the gradient of the log of the power-posterior, log(p(D|M )) is equivalent to the expected loglikelihood of θ, where the expectation is taken over the power posterior. This can be approximated by selecting r values of β ranging from 0 to 1, estimating the power posterior and using Simpson's Rule to estimate the model evidence (equivalent to the area under the curve). Note that this estimate can become computationally intensive as each value of β requires an MCMC output of sufficient length to limit estimation error, and r needs to be sufficiently large to reduce discretization error from the values of β. In this application r = 50 and 100,000 iterations of the MCMC chain were used. Methods to improve the efficiency of TI approximations are an important area of research.
When two models M 0 and M 1 are being compared, the Bayes Factor refers to which on a log-scale is log(p(D|M 1 )) − log(p(D|M 0 )). Values greater than 1 favor M 1 , while values less than 1 favor M 2 , although the cut-off chosen may vary between applications. Here simulations (described below) indicated that a cut-off of 1.00 may limit type II errors.

S4. Simulations to test the model framework
Simulated data were used to test whether the model framework was able to accurately estimate parameter values and identify the correct model. Values of ω and α were varied and 100 simulations were generated for each parameter set. The model parameters were re-estimated and the median and 95% CI were compared to the simulated values; simulations correctly re-estimating parameters included the values within the 95% CI. Variation in ω was tested using Bayes Factors between the simple model and a model that included variation in ω, if the Bayes Factors were >1.00 then there was evidence for the more complex model being a better fit to the data. These simulations were repeated 10 times for each parameter set. Figure S1 illustrates that estimates of surveillance sensitivity were both accurate and precise for outputs from the model without covariates. A difference in environmental sensitivity (ω) between districts could be detected if the difference was >10% and if AFP sensitivity was near 50% (Table  S2). For example when AFP sensitivity (α) was 50% and the difference in environmental sensitivity was greater than 20% all simulations selected the model with variation in sensitivity of environmental sampling, when the simulated difference in environmental sensitivity was less pronounced model selection was less accurate. When α = 0.1 model selection was still able to detect a difference in sensitivity when the true difference in ES sensitivity was more than 20% but at values less than 20% there was insufficient discriminatory power to correctly identify variation in ω. S1 Figure. Comparison of simulated (x-axis) and estimated (y-axis) values to compare the accuracy and precision of (A) AFP and (B) environmental surveillance sensitivity. For clarity the simulated values are jittered across the x-axis. The grey lines indicate the 95% credible intervals of the posterior distribution for the sensitivity parameters.
1 Tables S1 Table. Environmental sampling sites in Pakistan. Catchment size estimates were obtained from the Novel-t website. Estimates of sensitivity are from the model assuming a mixed effects structure.  Figure. Prior and posterior distributions for each of the parameters in the best-fitting model A) The infection rate (λ) B) the recovery rate (γ) C) Sensitivity of environmental surveillance (ω) for each sampled site D) Sensitivity of AFP surveillance (α) for each district.