Journal of Theoretical Biology 364 (2015) 266–274

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Transmission dynamics of oral polio vaccine viruses and vaccine-derived polioviruses on networks Jong-Hoon Kim a,b,n, Seong-Hwan Rho b a b

International Vaccine Institute, 1 Gwanak-ro, Gwanak-gu, Seoul, Korea 151-742 Simulacre Modeling Group, 4 Baekbeom-ro 45-gil, Yongsan-gu, Seoul, Korea 140-897

H I G H L I G H T S


OPV viruses may outcompete more transmissible, emerging cVDPVs in the short term.
Efforts to reduce disease transmission may increase the risk of cVDPV outbreaks.
Increasing vaccine coverage rate can be effective at preventing cVDPV outbreaks.

art ic l e i nf o

a b s t r a c t

Article history: Received 23 June 2014 Received in revised form 16 September 2014 Accepted 17 September 2014 Available online 28 September 2014

One drawback of oral polio vaccine (OPV) is the potential reversion to more transmissible, virulent circulating vaccine-derived polioviruses (cVDPVs), which may cause outbreaks of paralytic poliomyelitis. Previous modeling studies of the transmission of cVDPVs assume an unrealistic homogeneous mixing of the population and/or ignore that OPV viruses and cVDPVs compete for susceptibles, which we show is a key to understanding the dynamics of the transmission of cVDPVs. We examined the transmission of OPV viruses and cVDPVs on heterogeneous, dynamic contact networks using differential equation-based and individual-based models. Despite the lower transmissibility, OPV viruses may outcompete more transmissible cVDPVs in the short run by spreading extensively before cVDPVs emerge. If viruses become endemic, however, cVDPVs eventually dominate and force OPV viruses to extinction. This study improves our understanding of the emergence of cVDPVs and helps develop more detailed models to plan a policy to control paralytic polio associated with the continued use of OPV in many countries. & 2014 Elsevier Ltd. All rights reserved.

Keywords: OPV VDPV Network Reversion Transmission model

1. Introduction Polio epidemics first appeared in North America and Europe in the late 18th century (Nathanson and Kew, 2010) and, when the Global Polio Eradication Initiative was launched in 1988, the disease was endemic in more than 125 countries with the estimated 350,000 cases annually (World Health Organization, 1997). Since then, combating this disease has been successful and since 2012, the number of countries with endemic polio decreased down to three: Pakistan, Afghanistan, and Nigeria (World Health Organization, 2014). At the center of the efforts to eradicate the disease is oral polio vaccine (OPV) introduced in 1960s, which provides a less expensive, easy-to-administer option compared with the previously available inactivated poliovirus vaccine (IPV). OPV contains live, attenuated viruses, which may grow and induce n Corresponding author at: International Vaccine Institute, 1 Gwanak-ro, Gwanak-gu, Seoul, Korea 151-742. Tel.: þ 82 2 881 1295; fax: þ 82 2 881 2803. E-mail addresses: kimfi[email protected] (J.-H. Kim), [email protected] (S.-H. Rho).

http://dx.doi.org/10.1016/j.jtbi.2014.09.026 0022-5193/& 2014 Elsevier Ltd. All rights reserved.

immunity in the gut of vaccine recipients and transmit to the contacts leading to secondary immunity. One drawback of OPV is that viruses may evolve to regain the properties of their parental viruses (i.e., wild polio virus (WPV)) and cause outbreaks of paralytic polio. Recent and historic experiences of circulating vaccine-derived poliovirus (cVDPV) outbreaks in Nigeria (Wassilak et al., 2011), Indonesia (Estivariz et al., 2008), Hispaniola (Kew et al., 2002), Philippines (World Health Organization 2001), Madagascar (Rousset et al., 2003), China (World Health Organization, 2004; Liang et al., 2006), and retrospective finding of potential cVDPV outbreaks in Egypt (Yang et al., 2003), Poland (Martin et al., 2001) and Belarus (Korotkova et al., 2003) illustrate that the virus species that are as transmissible and neurovirulent as WPV may emerge from OPV. OPV continues to be used in many countries in routine immunizations or supplementary immunization activities (SIAs), which makes it important to understand how effectively we can reduce the risk of cVDPV outbreaks. Evidences suggest that cVDPV outbreaks occurred in populations where the vaccine coverage has dropped after the apparent absence of circulating indigenous

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

WPVs of the same serotype, which likely caused the population immunity to decrease low enough to allow for OPV viruses to circulate and evolve to cVDPVs (Kew et al., 2004; Dowdle et al., 2003; Kew et al., 2005; Dowdle and Kew, 2006). One recent study provides an extensive review of the evidence and potential mechanisms on the emergence of more transmissible and neurovirulent viruses in the OPV recipients (Duintjer Tebbens et al., 2013b). However, as the authors acknowledge, transmission of OPV (-derived) viruses between individuals can be more critical for the cVDPV emergence and subsequent outbreaks than reversion within the vaccine recipient and knowledge on what influences the transmission dynamics of OPV viruses and the emerging cVDPVs is still lacking. Mathematical models are useful for understanding the dynamics of infectious disease transmission (Anderson and May, 1992) and several studies used mathematical models to examine the dynamics of the OPV viruses and the cVDPV transmission (Duintjer Tebbens et al., 2013b; Wagner and Earn, 2008; Duintjer Tebbens et al., 2014). One study analyzed the risk of the cVDPV outbreak as a function of OPV coverage rate, but unrealistically assumed homogeneous mixing of individuals and disregarded the secondary transmission of OPV viruses (Wagner and Earn, 2008). Other studies explicitly included the secondary transmission and the reversion of OPV viruses to model historic cVDPV outbreaks or explore a potential impact of including older children and adults in SIAs (Duintjer Tebbens et al., 2014; Duintjer Tebbens et al., 2013b). These studies also assume homogeneous mixing in a subpopulation or a “mixing site” and generally do not provide the insights on what influences transmission dynamics of OPVs and emerging cVDPVs. In this paper, we study the transmission dynamics of OPV viruses and emerging cVDPVs on dynamic, heterogeneous networks of contacts between individuals to overcome the limitation of the previous models that assume random, one-off contacts in a (sub) population or a mixing site (Duintjer Tebbens et al., 2013b; Wagner and Earn, 2008; Duintjer Tebbens et al., 2014). To model contact networks, we mainly use a differential equation-based (DEB) model with probability generating functions (PGFs) to bypass extensive simulations required for the individual-based (IB) model, which is, however, also used to explore the effects of finite population size. We evaluated the models to identify various social and biological aspects that substantially influence the spread of OPV viruses and cVDPVs. In particular, we examined the effects of the difference in transmissibility between the OPV virus and the cVDPV, the dynamics and the distribution of contacts between individuals, finite population size, and population demographics. In the Methods section, we describe the model formulations and provide rationale for the choice of input values. In the Results section, we present figures and diagrams to illustrate the relative dominance between OPV viruses and cVDPVs in infecting people under varying assumptions. Finally, we discuss the robustness of our inferences and the implications of our findings to policy decisions.

2. Methods Fig. 1 shows a schematic of the models analyzed. Individuals may exist in four mutually exclusive states; susceptible (S), infected with the OPV virus (V), infected with the cVDPV (I), and recovered from infection (R). Susceptibles may be infected with either OPV viruses or cVDPVs, and OPV viruses revert to cVDPVs at a constant rate. Individuals become infectious as soon as they become infected and remain so until recovery, which occurs at a constant rate. In the DEB model, susceptible population becomes infected with OPV viruses or cVDPVs at the rate defined by a set

267

Fig. 1. Schematic for the models analyzed. Susceptibles (S) may be infected with either cVDPVs (I) or OPV viruses (V). OPV viruses revert to cVDPVs (i.e., V-I) at rate δ and infected people recover at rate γ. People die at rate μ independent of infection status. In the DEB model, S becomes V or I at the rate of rpIθg0 (θ) and ϕrpVθg0 (θ), respectively (see Appendix). In the IB model, each susceptible individual becomes infected with OPV viruses and cVDPVs depending on the number and infection status of the neighbors. Explanations about the variables and parameters appear in Tables 1 and 2.

of nonlinear, ordinary differential equations with PGFs whereas in the IB model, each susceptible individual has varying rates of infection depending on the number and infection status of the contacts. 2.1. Model inputs The model inputs that are applied to both the DEB and the IB models (Table 1) include the transmission rate of OPV viruses and cVDPVs per day (rϕ and r, respectively), the recovery rate of infected people (γ), reversion rate of OPV viruses to cVDPVs (δ), swapping rate for a contact (ρ), and natural birth (death) rate (μ). The model assumes pre-determined distributions for contact networks (i.e., Poisson, exponential, or power law), which remain constant over time although the identities of contacts may change. We explore these inputs over a wide range since our main purpose is to examine the model that is not limited to any one specific situation, except although baseline values are based on the available data. Details about the history of development of OPV and the difference among OPV viruses, WPVs, and cVDPVs appear elsewhere (Kew et al., 2004; Duintjer Tebbens et al., 2013b; Fine and Carneiro, 1999; Sutter et al., 2008; Agol, 2006). In brief, Albert Sabin produced OPV through the successive passages of WPVs in animal cells (e.g., monkey kidney cells) whereby WPVs were gradually attenuated both in transmissibility and neurovirulence. Being alive, however, OPV viruses can regain the properties of the WPV in OPV-infected individuals, although much uncertainty exists regarding whether OPV viruses fully regain properties of the WPV and if so, how quickly the full reversion occurs. In reality, OPV isolates exist in varying degrees of reversion, measured as the differentiation of the VP1 region of the virus genome (Kew et al., 2005). For simplicity, however, we disregard OPV-derived viruses with varying degrees of reversion and assume only two types of viruses: the OPV virus and the cVDPV. We simplify the process of reversion using a single parameter δ that represents the rate at which OPV viruses revert to cVDPVs within the host. We assume that co-infection with OPV viruses and cVDPVs do not occur although we acknowledge that biological processes underlying infection will likely be more complicated. Individuals recover at a rate γ regardless of virus strains with which they are infected. A review of literature seems to suggest a slightly shorter period of infection for OPV viruses than cVDPVs (Alexander et al., 1997), but the difference may not be significant, mainly because of the lack of data (Duintjer Tebbens et al., 2013a; Duintjer Tebbens et al., 2011). Individuals infected with OPV viruses and cVDPVs may transmit infection to their contacts at ϕr and r rates, respectively. Here, 0 rϕr1 represents the relative transmissibility of the OPV virus compared with the cVDPV, which will likely vary by serotype of polioviruses (Fine and Carneiro, 1999). We assume that once they recover from infection,

268

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

Table 1 Model inputs. Symbol Meaning

Range

References

δ

The rate at which OPV viruses revert to cVDPVs per day

1/60–1/730

ra ϕ

Transmission rate of the cVDPV per day Relative transmission rate of the OPV virus compared with the cVDPV Rate of recovery from infection with the OPV virus or the cVDPV per day Natural birth or death rate per day

6–14 for R* 0.25–0.6

(Estivariz et al., 2008; Kew et al., 2002; Dowdle et al., 2003; Dowdle and Kew, 2006; Duintjer Tebbens et al., 2013b; Agol, 2006; Minor and Dunn, 1988; Dunn et al., 1990; Contreras et al., 1992; Dane et al., 1961) (Duintjer Tebbens et al., 2005; Patriarca et al., 1997) (Duintjer Tebbens et al., 2013b; Eichner and Dietz, 1996; Duintjer Tebbens et al., 2013)

1/35

(Duintjer Tebbens et al., 2005)

1/(57  365)– 1/(80  365) 0–1

(World Health Organization)

μ ν ρ

Rate at which the contact exchanges per day a

NA

h   i ð1Þ ρ ρ Rn ¼ r þ ρr þ μ g″ g0 ð1Þ 1 þ μ þ μ (Volz and Meyers, 2009).

individuals remain immune to infection for simplicity although studies show that the immunity wanes in people vaccinated with OPV viruses and in people vaccinated with IPV as people age (Nishio et al., 1984; Bottiger, 1990). 2.2. Neighbor exchange model To incorporate heterogeneous, dynamic contact networks while still keeping the model simple, we use the neighbor exchange (NE) model, of which details appear elsewhere (Volz and Meyers, 2007; Volz, 2008; Kamp, 2010). Briefly, the NE model assumes a pre-determined distribution of contacts (“arcs” in network terms) across the individuals (“nodes” in network terms) using a PGF. Identities of contacts may change over time, but the number of arcs per node and thus the overall distribution of contacts remain constant. For example, two nodes, ego1 and ego2, with distinct contacts (ego1, alter1) and (ego2, alter2) may exchange contacts such that resulting contacts are (ego1, alter2) and (ego2, alter1). Any given contact (ego1, alter1) will be reassigned to (ego1, alter2) at a constant rate ρ. In other words, two contacts are rearranged at a rate ρ/2. The model comprises dynamic variables (Table 2) indicating the fraction of susceptibles with degree k ¼1 (θ) and the probability that the contacts of a susceptible are infected with cVDPVs (pI) or infected with OPV viruses (pV). In addition, we use such variables as MXY that indicates the fraction of contacts between ego from the set X and alter from the set Y, where X,Y ¼S, V, I, R. (see Appendix A1) As in the previous studies (Volz and Meyers, 2007; Volz, 2008), we assume that each contact of a susceptible node (ego, alter) where ego A S has a uniform probability pI that alter A I, and a probability pV that alter A V. This implies that a susceptible node has an expected number of kpI contacts with nodes A I and kpV contacts with nodes A V. Thus, the instantaneous rate of infection for a degree k node is λk ðt Þ ¼ rkpI þ rϕkpV :

Table 2 Dynamic variables in the model. Notation

Meaning

θ S V I R g(x) g0 (x)

Fraction of susceptibles with degree k ¼ 1 Fraction of susceptible nodes Fraction of nodes infected with OPV viruses Fraction of nodes infected with cVDPV viruses Fraction of nodes that recovered from infection Probability generating function, g(x) ¼ a0 þ a1x þ a2x2 þa3x3 þ... g 0 (x) ¼ a1 þ 2a2x þ3a3x2 þ...

distribution ak will be generated by the PGF g ðxÞ ¼ a0 þ a1 x þa2 x2 þ a3 x3 þ⋯

The dummy variable x serves as placeholder for dynamic variables and the mean of each distribution is given by g0 (1). Poisson degree distributions, ak ¼ zk e  z =k!; k Z 0, are generated by gðxÞ ¼ ezðx  1Þ α and power-law with cut-off, ak ¼ k =∑κi ¼ 1 i  α ; 1 rk rκ, is genα k α κ erated by gðxÞ ¼ ∑k k x =∑i ¼ 1 i , where κ denotes the maximum degree. Let θ¼u1(t) be the fraction of degree k ¼1 nodes that remain susceptible at time t. Then, uk(t) ¼θk and the fraction of nodes that remain susceptible at time t is S ¼ a0 þ a1 θ þ a2 θ2 þ a3 θ3 þ ⋯ ¼ ∑ ak θk ¼ g ðθÞ

To incorporate population turnover, we assume that people are born and die at the same rate and the number of contacts of newborns is a random number from the pre-determined distribution for the population, and newborns connect with existing individuals uniformly at random. Therefore, the expected population size and the expected distribution of contacts remain constant, which makes PGF constant against time. The following set of equations describes dynamics of individuals in S, V, I, and R states:

ð1Þ

_ 0 ðθÞ  μðS  NÞ; S_ ¼ θg V_ ¼ S_ V  ðμ þ δ þ γ ÞV; _I ¼ S_ I þ δV  ðμ þ γ ÞI;

¼ exp Z ¼ exp

τ¼0 t

τ¼0 t τ¼0

   rkpI ðτÞ þ rϕkpV ðτÞ dτ    rpI ðτÞ þrϕpV ðτÞ dτ



k ð2Þ

Given uk(t), we can compute the fraction of nodes remaining susceptible at time t for various degree distributions. A degree

ð4Þ

k

Let uk(t) denote the fraction of degree k nodes remaining susceptible at time t. Then, uk(t) is given as follows: Z t  uk ðtÞ ¼ exp  λk ðτÞdτ Z

ð3Þ

R_ ¼ γ ðV þ I Þ  μR;

ð5Þ

where S_ V and S_ I denote new infection rates of susceptible individuals by the OPV viruses and the cVDPVs, respectively. Additional equations including those for S_ V and S_ I that are needed to close the system appear in the Appendix A1. To explore the effects of finite population size, we use an individual-based (IB), discrete-time, Markov model. We generate the contact networks using the configuration model, as described

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

in (Molloy and Reed, 1995). In brief, we assign the number of “halfedges” to each individual by drawing a random number from a pre-determined distribution. If the total number of half-edges is not even, we make it so by choosing an individual uniformly at random among those with one or more half-edges and reducing the number of the half-edges of the individual by one. We then connect half-edges randomly until no half-edges are left. We restrict self-edges (i.e., both ends of an edge are connected to the same individual) and multiple edges between a pair of individuals although their frequency becomes negligible as population size becomes large. We then randomly select a fraction ε (e.g., 0.1%) of population and set them initially infected with the OPV viruses. Then, following events occur at each time step with Δ indicating the length of a time step in days: 1. OPV viruses revert to cVDPVs with probability Δδ 2. Individuals infected with OPV viruses and cVDPVs infect susceptible contacts with probability Δrϕ and Δr, respectively 3. Infected individuals recover with probability Δγ 4. Contacts swap with probability Δρ 5. Susceptible births occur with probability Δμ per person 6. People die with probability Δμ 3. Results Fig. 2 illustrates that model behavior may be in one of four qualitatively different regimes depending on basic reproductive ratio of the cVDPV on network (Rn) and the relative transmission rate of the OPV virus compared with the cVDPV (ϕ). First of all, epidemics are possible only if Rn is over one (i.e., solid vertical line). Rn is different from the more commonly known quantity, basic reproductive ratio (R0) and represents the expected number of secondary infections from an individual who is infected early in the epidemic, but not the very first case in the population (Volz and Meyers, 2009). Until Rn or ϕ increases further such that basic reproductive ratio of the OPV virus on network is also over threshold value of one, only cVDPVs can cause epidemics (i.e., cVDPVs-only regime). Above this threshold, both OPV viruses and cVDPVs may cause epidemics, but more transmissible cVDPVs may be dominant (i.e., cVDPVs-dominant regime). If Rn or ϕ increases further, the system may be in a regime where OPV viruses are transmissible enough to infect a larger portion of the population than cVDPVs (i.e., OPV-dominant regime). Despite the lower

Fig. 2. Model behavior may be in one of four qualitatively different regimes depending on basic reproductive ratio of cVDPV on network (R*) and the relative transmission rate of the OPV virus compared with the cVDPV (ϕ). “Dominant” means a larger cumulative number of infections based on results of the DEB model that assume Poisson-degree distributions with mean degree z¼ 20. γ¼1/35, δ¼ 1/365, μ¼ 0, S(0)¼ 0.999, V(0)¼ 0.001, I(0)¼ R(0) ¼0. The vertical arrow indicates probable values for ϕ based on the literature (i.e., ϕ ¼0.25–0.6).

269

transmissibility, OPV viruses can be dominant because they can spread extensively before a substantial fraction of OPVs revert to cVDPVs. In addition to Rn and ϕ, other model inputs can influence transition across different regimes which we cover later in Fig. 4. Fig. 3(a)–(c) illustrate the system behavior explained in Fig. 2. We assume that the distribution of contacts follows Poisson distribution with mean of 20 and contacts do not change over time (i.e., ρ¼0). We assume that the mean degree of the network is arbitrarily large enough to explore the wide range of values for Rn. Fig. 3(a) shows the cVDPV-only regime. cVDPVs emerge through the reversion of OPV viruses and infect a significant fraction of the population whereas OPV viruses soon die out because their transmissibility is under the threshold. Fig. 3(b) shows that increasing Rn and ϕ make transmissibility of both OPV viruses and cVDPVs above the threshold, but cVDPVs produce a larger epidemic. In this regime, OPV viruses spread first, but the extent to which OPV viruses spread is not large enough to limit the spread of cVDPVs. Accordingly, later emerging cVDPVs that are more transmissible produce a larger epidemic. Fig. 3(c) shows that increasing Rn and ϕ even further makes OPV viruses dominant. In this regime, the rate of spread of OPV viruses is fast enough to infect a majority of the population before cVDPVs emerge and spread significantly. Comparison of Figs. 3(b) and (c) shows that cVDPV epidemics may be greater in populations where viruses are less transmissible (e.g., populations with good hygiene) than in populations where viruses are more transmissible (e.g., populations with poor hygiene). This is because a population with good hygiene will benefit less from secondary and tertiary, etc. transmissions of OPV viruses and thus a larger proportion of population will remain susceptible which emerging, more transmissible cVDPVs can infect afterwards. Fig. 4(a)–(c) illustrate additional model inputs that influence the dynamics of the transmission of OPV viruses and cVDPVs. Fig. 4(a) shows that increasing the reversion rate (δ) may shift the dominance from OPV viruses to cVDPVs. The cVDPVs emerge from OPV viruses more rapidly as δ increases and can overtake OPV viruses. Fig. 4(b) shows that increasing the rate of contact swapping ρ speeds up infection spread and thus enhances the transmission of OPV viruses before cVDPVs emerge and spread, which implies that the models assuming one-off contacts will likely underestimate the spread of cVDPVs. Fig. 4(c) shows that changing the contact distribution from Poisson to exponential to power law with the same mean may shift the dominance from cVDPVs to OPV viruses. Prior studies showed that as the degree distribution changes from Poisson to exponential to power law, time to the peak incidence shortens while the epidemic size remains largest for Poisson distribution (Volz, 2008; Barthelemy et al., 2005). Shorter time to the peak (i.e., more rapid transmission) imply that OPV viruses can spread more extensively before cVDPVs emerge. Accordingly, OPV viruses may be dominant in a population where contacts between people follow exponential or power law distributions. Fig. 4(d) shows that increasing death rate μ (i.e., decreasing life expectancy) increases the transmission of cVDPVs relatively more than that of OPV viruses. This holds even when OPV viruses dominant in the short term (not shown). Increased death rate raises the chance of the endemic circulation of viruses, which eventually makes cVDPVs dominant. This is in agreement with the consensus that pathogens with the higher R0 will force the other pathogens to extinction in the long run from the studies of competing pathogens (Keeling and Rohani, 2008). Abrupt increases in the population infected with cVDPVs in Fig. 4(d) arise because smaller cVDPV epidemics following the first largest epidemic that we showed in Fig. 3 occurs with increasing frequency as μ increases. For example, as the life expectancy, 1/μ/ 365, becomes shorter than 50 years, one additional cVDPV epidemic occurs at the end of 10 simulated years for which we

270

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

Fig. 3. (a)-(c). The fraction of population with OPV viruses (solid line) and cVDPVs (dashed line) based on the DEB models change in response to the change in R* or ϕ. S(0)¼ 0.999, V(0)¼ 0.001, I(0) ¼R(0) ¼ 0, μ¼ 0, γ ¼1/35, δ ¼ 1/365. Contact networks are Poisson-distributed with mean z of 20. (a) R* ¼3, ϕ ¼0.25, (b) R* ¼ 6, ϕ ¼0.6, (c) R* ¼ 14, ϕ¼ 0.6.

Fig. 4. (a)-(d). Various model inputs may change the fraction of OPV-infected (solid line) and cVDPV-infected (dashed line) population shifting the dominance between them based on the DEB models. S(0)¼ 0.999, V(0)¼ 0.001, I(0) ¼R(0)¼ 0, γ ¼1/35, δ¼ 1/365, μ¼0 (for (a)-(c)). Contact networks are Poisson-distributed with mean z¼ 20 for (a), (b), and (d). R* ¼ 10, ϕ¼ 0.6. For (b), R* ¼ 10 for ρ ¼ 0. Note that R* increases as ρ increases. For (c), λ ¼ 20.5 for exponential distribution, and α ¼0.976, κ ¼ 100 for power law distribution with exponential cut off. Cumulative fraction of population infected is computed for a 1-year period for (a)-(c), and for a 10-year period for (d).

have computed the incidence. Similarly, if the life expectancy becomes shorter than 30 years, two additional cVDPV epidemics with decreasing size occur for the same simulation period. A

similar trend continues to occur as the life expectancy decreases and almost stable endemic equilibrium is reached if the life expectancy is around 10 years.

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

271

Fig. 5. (a)-(d). Dynamics of OPV viruses and cVDPVs under routine immunization or SIAs. S(0)¼ 0.999, V(0) ¼0.001, I(0) ¼ R(0)¼ 0, μ¼1/(365  57), γ ¼ 1/35, δ ¼ 1/365, R* ¼10, ϕ¼ 0.6. Contact networks are Poisson-distributed with mean z ¼20. Routine vaccinations (vaccinating 50% of births) increase the level of OPV-infected people (a), but decrease the level of cVDPV-infected people (b). Given that the cumulative number of vaccine doses is the same as routine vaccinations, SIAs also reduce the level of cVDPVinfected people (c), and can be more effective than routine vaccinations (d). SIAs that correspond to complete routine vaccinations (empty triangles) are missing because it is impossible to vaccinate the same number of people using SIAs due to slightly more transmissions of cVDPVs if the vaccination program is implemented from the second year.

In Fig. 5, we explore the effects of simple vaccination strategies on the transmission of OPV and cVDPVs (see Appendix A1 for models with vaccination programs). Fig. 5(a) shows that vaccinating 50% of births (i.e., routine vaccination) increases the level of OPV-infected people, but decreases the level of cVDPV-infected people. In particular, as shown in Fig. 5(b), routine vaccination can prevent the cVDPV outbreaks that otherwise would occur around 10–25 years. Fig. 5(c) shows that SIAs that last a week at the beginning of each year reduces the size of the first largest cVDPV epidemic significantly, and prevents subsequent smaller epidemics as do routine vaccinations. Here, the total number of vaccine doses is the same for both routine vaccinations and SIAs. Fig. 5(d) shows that as the vaccine coverage increases, the number of cumulative cVDPV infections decreases. Comparison of the solid line and the filled triangles shows that SIAs are more effective at reducing the transmission of cVDPVs than routine vaccinations. This supports the effectiveness shown in historic SIAs in developing countries (Olive et al., 1997; Sabin, 1991; Birmingham et al., 1997; Sutter and Maher, 2006) and also the current recommendation by WHO (World Health Organization, 2014). The difference between routine immunizations and SIAs arises because SIAs, but not routine vaccinations, reduce the first largest epidemic significantly. If we assume that the vaccination programs are implemented from the second year (i.e., missing the first epidemic), then the difference in their effects between routine vaccinations and SIAs are negligible, as seen by the comparison between the dashed line and the empty triangles in Fig. 5(d). Fig. 6 shows that finite population size favors the spread of OPV viruses over emerging cVDPVs, in particular, in a small population (e.g., r10,000). This arises mainly because reduced population size implies a smaller size of the population infected with OPV

viruses, which in turn leads to an even smaller number of emerging cVDPVs that will likely experience stochastic die-outs more frequently. As population size increases (e.g., 100,000), the IB model produces almost the same results as the DEB model.

4. Discussion We have illustrated that the transmission dynamics of OPV viruses and cVDPVs are influenced by various social and biological factors. In particular, OPV viruses may outcompete more transmissible cVDPVs by spreading further before cVDPVs emerge. This happens if conditions favor the transmission of infectious disease such as low hygiene or frequent change in contacts or increased proportion of individuals with a large number of contacts. This implies that ironically, efforts to reduce the transmission of a disease in general (e.g., improving sanitation) may lead to a higher risk of a cVDPV outbreak if those efforts limit the transmission of OPV viruses, but are not strong enough to prevent the spread of cVDPVs. Despite that, conventional wisdom holds that increasing OPV coverage rate or population immunity is effective at preventing cVDPV outbreaks (Fig. A1 in the Appendix). In addition, realistic epidemic die-out in a small population will likely reduce the spread of emerging cVDPVs. To prevent cVDPV outbreaks, it is important to reduce the time during which viruses circulate in the population because upon continued transmission, more transmissible cVDPVs will eventually dominate and force OPV viruses to extinction even if OPV viruses are dominant during the early stages of the epidemic. This is in agreement with the consensus that the disease with the higher basic reproduction ratio, R0, will

272

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

force the other strain to extinction from the studies of competing pathogens (Keeling and Rohani, 2008). This study has many limitations. First, the model assumes that the OPV virus reverts to the cVDPV as a one-step process although in reality reversion will likely occur through accumulated mutations (Kew et al., 2004; Kew et al., 2005; Agol, 2006; Agol, 2006). This simplifying assumption leads to that the waiting time before reversion is exponentially distributed and thus a substantial fraction of OPV viruses revert to cVDPVs almost instantaneously, which in turn will likely increase the chance of cVDPV transmission and thus overestimate the risk of the cVDPV outbreak. Second, we modeled a single arbitrary serotype of poliovirus although in reality three different serotypes (types 1, 2, and 3) exist. OPV usually contains all of three serotypes of attenuated polioviruses with serotypes presenting different transmissibility and neurovirulence relative to their progenitors. For example, it seems that relative transmissibility of OPV viruses is highest for type 2 (Fine and Carneiro, 1999) and reversion rate is higher for type 2 and type 3 than for type 1 (Kew et al., 2005). However, there seems to be little cross protection between three serotypes (Sutter et al., 2008; Burnet and Macnamara, 1931), which implies we can model each serotype separately without loss of generality. Third, we assumed only simple networks such as Poisson or power-law degree distributions although real-world contact networks will likely be more complicated. For example, human contact networks are often clustered, i.e., there is a high probability that “the friend of my friend is also my friend” (Watts and Strogatz, 1998; Newman, 2003). A recent study shows that clustering tends to slow the spread of infection and interacts with the degree distributions to influence the dynamics of infection transmission (Volz et al., 2011). If clustering slows the spread of polioviruses, this may mean that clustering will favor the dominance of cVDPVs over OPV viruses, at least in the short term, since more cVDPVs will arise from OPV viruses before OPV viruses spread extensively. Also, our network structure does not specify any physical locations (or mixing sites) although one study, using the IB model (Rahmandad et al., 2011), shows that different network structures including mixing site assumptions influence the dynamics of the transmission of polioviruses. Assuming simple networks, however, facilitates the analysis of the model and the neighbor exchange model allows us to assume various types of networks in the context of the DEB model and thereby bypass extensive simulations required for the IB models. Fourth, we assumed that people have the same immunity to the infection of OPV viruses or cVDPVs although prior studies suggest that people will likely have different susceptibility to infection, the length

of infectious period, etc. depending on the exposure history (Duintjer Tebbens et al., 2013b; Duintjer Tebbens et al., 2005). This issue of population immunity is partly addressed in Fig. A1(a), which shows that increasing the proportion of population that is immune to polioviruses reduces the transmission of both OPV viruses and cVDPVs. Fifth, we ignored that some of the people who harbor cVDPVs are originally infected with OPV viruses, who may present different clinical outcomes from the people who are infected with cVDPVs in the first place. If clinical outcomes are milder in the people who harbor cVDPVs but were originally infected with OPV viruses, the magnitude of the simulated cVDPV outbreaks will be an overestimation. Finally, we assume that only OPV vaccination as an option of the vaccination program to focus on the dynamics of the competition between OPV viruses and cVDPVs. In reality, IPV is also a viable option and seems to have differential effects on the transmission of OPV viruses and cVDPVs (Fine and Carneiro, 1999). Despite these limitations, we believe that this study provides insights on the cVDPV emergence and has implications for the control of the potential cVDPV outbreaks. This model can serve as a good starting point from which more detailed models can be developed to plan a policy to control paralytic poliomyelitis.

Acknowledgment J-HK and S-HR thank anonymous reviewers as well as Justin Im from International Vaccine Institute for their helpful and insightful comments.

Appendix. Effects of initial population immunity and initial OPV coverage rate Given that polioviruses have been circulating for a long time, it is natural to assume that some portion of the population is already immune when OPV is introduced. Also, the proportion of population to whom OPV is given (i.e., OPV coverage rate) can vary and will likely be much higher than we have assumed in the main text (i.e., 0.1%). These realities will influence the transmission dynamics of OPV viruses and cVDPVs. Fig. A1(a) shows that the fraction of population infected with of OPV viruses or cVDPVs decreases with the increased initial proportion of population that is immune, R(0). Fig. A1(b) shows that increased initial proportion of population infected with the OPV viruses, V(0), decreases the proportion of population infected with cVDPVs.

Fig. 6. Fraction of population infected with (a) OPV viruses and (b) cVDPVs from IB models deviate from those from DEB model as population size decreases. Empty triangles, circles, and squares with bars represent means and 95% confidence intervals of 100 simulation runs of IB models with the population size of 103, 104, and 105, respectively. Solid lines indicate the results from DEB models. S(0) ¼0.999, V(0)¼ 0.001, I(0)¼ R(0) ¼0, μ¼ 0, γ ¼ 1/35, δ ¼ 1/365, ϕ¼ 0.6, and R* ¼ 6. Contact networks are Poisson-distributed with mean z ¼ 20. Solid lines indicate the results from DEB models with the same input values.

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

273

Fig. A1. The fraction of population infected with OPV viruses or cVDPVs from DEB models decreases as the initial proportion of population that are immune, R(0), or the proportion of the population to whom OPVs are given, V(0), increases. For (a), V(0) ¼ 0.001, R(0) ¼ 0–0.6 and for (b), V(0)¼ 0.001–0.9, R(0)¼ 0. For both (a) and (b), I(0)¼ 0, S (0) ¼1–V(0)-I(0)-R(0), μ¼ 0, γ ¼ 1/35, δ ¼1/365, R* ¼10, ϕ ¼0.25 or 0.6. Contact networks are Poisson-distributed with mean z¼ 20.

Table A1 Parameters and dynamic variables in the model. Notation Meaning α β MX MXY pX

Fraction of OPV-vaccinated births Rate at which susceptibles are vaccinated with OPV Number of contacts where ego A X where X ¼S, V, I, R Number of contacts where ego A X and alter AY, for X, Y ¼ S, V, I, R The probability that an arc with a susceptible ego has an alter in set X ¼ S, V, I, R (i.e., pX ¼ MSX/MS), Average excess degree of S nodes with at least one infected partner that are connected to S nodes Average excess degree of S nodes with at least one infected partner that are connected to X nodes, for X¼ S, V, I

πS πS,X

(Table A1). We define β as follows: ( βn if t ðmod ωÞ r k β¼ 0 otherwise where t, ω, and k represent time, period of SIAs, and the length of SIAs, respectively. That is, β takes the value of βn if the remainder on dividing t by ω is less than or equal to k and is zero otherwise. For example, ω¼365 and k ¼7 imply that susceptibles are vaccinated at rate βn for a week every year. For given α, ω, and k, we numerically identified βn that produces the same cumulative number of vaccinated people.

References A1. equations for the model with the vaccination program

_ 0 ðθÞ μðS  ð1  αÞN Þ βS; S_ ¼ θg V_ ¼ S_ V  ðμþ δ þ γ ÞV þ αμN þβS; _I ¼ S_ I þ δV  ðμþ γ ÞI; R_ ¼ γ ðV þ I Þ μR;   θ_ ¼  rθ pI þ ϕpV ; S_ V ¼ rϕpV θg 0 ðθÞ; S_ I ¼ rp θg 0 ðθÞ; I

_ SV M _ S M SV   M p_ V ¼  þρ M V  pV ; MS MS 2 _ SI M _ S M SI   M  þ ρ M I  pI ; 2 MS MS     _ S ¼  r ϕp þ p M S þ μðN  SÞ g 0 ð1Þ M V I   _ SV ¼  r ϕp þ p M S π S;V þ rϕp M S π S;S  ðrϕ þ γ þ δÞM SV M V I V p_ I ¼

þ μg 0 ð1Þðð1  αÞV þ αS  2M SV Þ þ ðM SS  M SV Þβ;   _ SI ¼ r ϕp þp M S π S;I þ rp M S π S;S  ðr þ γ ÞM SI þ πM SV M V I I þ μðg 0 ð1Þð1  αÞI 2M SI Þ  βM SI ;   _ M SS ¼  2r ϕpV þ pI M S π S;S þ νð2g 0 ð1Þð1  αÞS 2M SS Þ  βM SS ; π S ¼ θg″ðθÞ=g 0 ðθÞ π S;X ¼ δS pX for X ¼ S; V; I: In the above equations, α and β represent the fraction of vaccinated births and the rate at which the susceptible people are vaccinated

Agol, V.I., 2006. Vaccine-derived polioviruses. Biologicals 34 (2), 103–108. Agol, V.I., 2006. Molecular mechanisms of poliovirus variation and evolution. Curr. Top. Microbiol 299, 211–259. Alexander Jr., J.P., Gary Jr., H.E., Pallansch, M.A., 1997. Duration of poliovirus excretion and its implications for acute flaccid paralysis surveillance: a review of the literature. J. Infect. Dis 175 (1), S176–S182. Anderson, R.M., May, R.M., 1992. Infectious Diseases of Humans: Dynamics and Control,. Oxford University Press, Oxford. World Health Organization, 2001. Acute flaccid paralysis associated with circulating vaccine-derived poliovirus, Philippines, 2001. Wkly. Epidemiol. Rec 76, 319–320. Barthelemy, M., et al., 2005. Dynamical patterns of epidemic outbreaks in complex heterogeneous networks. J. Theor. Biol. 235 (2), 275–288. Birmingham, M.E., et al., 1997. National immunization days: state of the art. J. Infect. Dis. 175 (1), S183–S188. Bottiger, M., 1990. Polio immunity to killed vaccine: an 18-year follow-up. Vaccine 8 (5), 443–445. Burnet, F., Macnamara, J., 1931. Immunological differences between strains of poliomyelitic virus. Brit. J. Exp. Pathol 12 (2), 57. Contreras, G., et al., 1992. Genetic characterization of Sabin types 1 and 3 poliovaccine virus following serial passage in the human intestinal tract. Biologicals 20 (1), 15–26. Dane, D., et al., 1961. Vaccination against poliomyelitis with live virus vaccines: 6. Changes in sabin type II oral vaccine virus after human passage. Brit. Med. J 2 (5247), 259–265. Dowdle, W., Kew, O., 2006. Vaccine-derived polioviruses: is it time to stop using the word rare? J. Infect. Dis 194 (5), 539–541. Dowdle, W.R., et al., 2003. Polio eradication: the OPV paradox. Rev. Med. Virol 13 (5), 277–291. Duintjer Tebbens, R.J., et al., 2005. A dynamic model of poliomyelitis outbreaks: learning from the past to help inform the future. Am. J. Epidemiol 162 (4), 358–372. Duintjer Tebbens, R.J., et al., Expert review on poliovirus immunity and transmission. Submitted, 2011. Duintjer Tebbens, R.J., et al., 2013a. Review and assessment of poliovirus immunity and transmission: synthesis of knowledge gaps and identification of research needs. Risk Anal 33 (4), 606–646. Duintjer Tebbens, R.J., et al., 2013b. Oral poliovirus vaccine evolution and insights relevant to modeling the risks of circulating vaccine-derived polioviruses (cVDPVs). Risk Anal. 33 (4), 680–702.

274

J.-H. Kim, S.-H. Rho / Journal of Theoretical Biology 364 (2015) 266–274

Duintjer Tebbens, R.J., et al., 2013. Characterizing poliovirus transmission and evolution: insights from modeling experiences with wild and vaccine-related polioviruses. Risk Anal. 33 (4), 703–749. Duintjer Tebbens, R.J., et al., 2013. Review and assessment of poliovirus immunity and transmission: synthesis of knowledge gaps and identification of research needs. Risk Anal. 33 (4), 606–646. Duintjer Tebbens, R.J., et al., 2014. The potential impact of expanding target age groups for polio immunization campaigns. BMC Inf. Dis 14, 45. Dunn, G., et al., 1990. Virus excretion and mutation by infants following primary vaccination with live oral poliovaccine from two sources. J. Med. Virol. 32 (2), 92–95. Eichner, M., Dietz, K., 1996. Eradication of poliomyelitis: when can one be sure that polio virus transmission has been terminated? Am. J. Epidemiol. 143 (8), 816–822. Estivariz, C.F., et al., 2008. A large vaccine-derived poliovirus outbreak on Madura Island–Indonesia, 2005. J. Inf. Dis 197 (3), 347–354. Fine, P.E., Carneiro, I.A., 1999. Transmissibility and persistence of oral polio vaccine viruses: implications for the global poliomyelitis eradication initiative. Am. J. Epidemiol. 150 (10), 1001–1021. Fine, P.E., Carneiro, I.A., 1999. Transmissibility and persistence of oral polio vaccine viruses: implications for the global poliomyelitis eradication initiative. Am. J. Epidemiol 150 (10), 1001–1021. Kamp, C., 2010. Untangling the interplay between epidemic spread and transmission network dynamics. PLoS Comput. Biol 6 (11), e1000984. Keeling, M.J., Rohani, P., 2008. Modeling Infectious Diseases in Humans and Animals, 408. Princeton University Press,, Princeton,. Keeling, M.J., Rohani, P., 2008. Modeling Infectious Diseases in Humans and Animals,. Princeton University Press, Princeton, pp. 107–108. Kew, O., et al., 2002. Outbreak of poliomyelitis in Hispaniola associated with circulating type 1 vaccine-derived poliovirus. Science 296 (5566), 356. Kew, O., et al., 2004. Circulating vaccine-derived polioviruses: current state of knowledge. B World Health Organ 82, 16–23. Kew, O., et al., 2005. Vaccine-derived polioviruses and the endgame strategy for global polio eradication. Ann. Rev. Microbiol 59, 587–635. Kew, O.M., et al., 2005. Vaccine-derived polioviruses and the endgame strategy for global polio eradication. Annu. Rev. Microbiol. 59, 587–635. Korotkova, E.A., et al., 2003. Retrospective analysis of a local cessation of vaccination against poliomyelitis: a possible scenario for the future. J. Virol 77 (23), 12460–12465. Liang, X., et al., 2006. An outbreak of poliomyelitis caused by type 1 vaccine-derived poliovirus in China. J. Inf. Dis 194 (5), 545–551. Martin, J., et al., 2001. Risks of reintroduction of polio after eradication: the vaccine origin of an outbreak of type 3 poliomyelitis. Dev. Biol 105, 83–92. Minor, P.D.,, Dunn, G., 1988. The effect of sequences in the 50 non-coding region on the replication of polioviruses in the human gut. J. Gen. Virol. 69 (Pt 5), 1091–1096. Molloy, M.,, Reed, B., 1995. A critical point for random graphs with a given degree sequence. Random Struct. Algor 6 (2-3), 161–180. Nathanson, N., Kew, O., 2010. From emergence to eradication: the epidemiology of poliomyelitis deconstructed. Am. J. Epidemiol 172 (11), 1213–1229.

Newman, M.E., 2003. Properties of highly clustered networks. Phys. Rev. E Stat. Nonlin. Soft Matter. Phys. 68 (2 Pt 2), 026121. Nishio, O., et al., 1984. The trend of acquired immunity with live poliovirus vaccine and the effect of revaccination: follow-up of vaccinees for ten years. J. Biol. Stand 12 (1), 1–10. Olive, J.M., Risi Jr., J.B., de Quadros, C.A., 1997. National immunization days: experience in Latin America. J Infect. Dis 175 (1), S189–S193. Patriarca, P.A., Sutter, R.W., Oostvogel, P.M., 1997. Outbreaks of paralytic poliomyelitis, 1976-1995. J. Infect Dis. 175 (1), S165–S172. Rahmandad, H., et al., 2011. Development of an individual-based model for polioviruses: implications of the selection of network type and outcome metrics. Epidemiol. Infect. 139 (6), 836–848. Rousset, D., et al., 2003. Recombinant vaccine-derived poliovirus in Madagascar. Emerg. Infect. Dis 9 (7), 885–886. Sabin, A.B., 1991. Perspectives on rapid elimination and ultimate global eradication of paralytic poliomyelitis caused by polioviruses. Eur. J. Epidemiol. 7 (2), 95–120. Sutter, R.W., Maher, C., 2006. Mass vaccination campaigns for polio eradication: an essential strategy for success. Curr. Top Microbiol. Immunol. 304, 195–220. Sutter, R.W., Kew, O.M., Cochi, S.L., 2008. Poliovirus vaccine—live. In: P. S.A., O. W.A., O.P.A. (Eds.), Vaccines. Saunders Elsevier, pp. 631–686. Volz, E., 2008. SIR dynamics in random networks with heterogeneous connectivity. J. Math. Biol 56 (3), 293–310. Volz, E., Meyers, L.A., 2007. Susceptible-infected-recovered epidemics in dynamic contact networks. P. Roy. Soc. B Bio 274 (1628), 2925–2933. Volz, E., Meyers, L.A., 2009. Epidemic thresholds in dynamic contact networks. J. Roy. Soc. Interf 6 (32), 233–241. Volz, E.M., et al., 2011. Effects of heterogeneous and clustered contact patterns on infectious disease dynamics. PLoS Comput. Biol. 7 (6), e1002042. Wagner, B.G., Earn, D.J., 2008. Circulating vaccine derived polio viruses and their impact on global polio eradication. B. Math. Biol 70 (1), 253–280. Wassilak, S., et al., 2011. Outbreak of type 2 vaccine-derived poliovirus in Nigeria: emergence and widespread circulation in an underimmunized population. J. Infect. Dis 203 (7), 898–909. Watts, D.J., Strogatz, S.H., 1998. Collective dynamics of ‘small-world’ networks. Nature 393 (6684), 440–442. World Health Organization, Life expectancy at birth. 〈http://www.who.int/gho/ mortality_burden_disease/life_tables/situation_trends_text/en/index.html〉. World Health Organization AU: Please provide the full details for this reference., Polio: The beginning of the end. 1997. World Health Organization, 2004. Laboratory surveillance for wild and vaccinederived polioviruses, January 2003–June 2004. Wkly Epidemiol. Rec 79, 393–398. World Health Organization, 2014. Polio vaccines: WHO position paper, January 2014. Wkly. Epidemiol. Rec. 89 (9), 73–92. Yang, C., et al., 2003. Circulation of endemic type 2 vaccine-derived poliovirus in Egypt from 1983 to 1993. J. Virol 77 (15), 8366.