Asymmetry in the presence of migration stabilizes multistrain disease outbreaks

In this work we study the stability of a multistrain disease model in two coupled populations. Multistrain diseases are diseases with multiple coexisting strains, such as influenza (Andreasen et al.,, 1997), HIV (Hu et al.,, 1996), and dengue (Ferguson et al., 1999b, ). We consider two kinds of strain interactions: cross immunity and antibody-dependent enhancement. When the disease infects an individual, his or her immune system creates serotype-specific antibodies, which will protect the individual against that serotype. However, antibodies may also give some cross-protection to the other serotypes (Halstead,, 2007). This reduced susceptibility to the other serotypes may be temporary. When the temporary cross immunity wanes, heterologous secondary infections are possible. Low level antibodies developed from primary infections may form complexes with the virus so that more cells are affected, and viral load is increased (Vaughn et al.,, 2000). This effect is called antibody-dependent enhancement (ADE), and it has been observed in vitro in diseases such as ebola (Takada et al.,, 2003) and dengue (Halstead and O’Rourke,, 1977). Throughout this work we make the hypothesis that ADE increases the infectiousness of secondary infective cases due to the higher viral load (Cummings et al.,, 2005; Schwartz et al.,, 2005). Alternative views of ADE as an increase in mortality associated with secondary infectives can be considered (Kawaguchi et al.,, 2003). We focus in this paper on the multistrain disease dengue, which is believed to exhibit both temporary cross immunity and antibody-dependent enhancement (Halstead,, 2007). Dengue is a subtropical mosquito-borne disease that exhibits up to four serotypes. It is widespread in tropical regions of southeast Asia, Africa, and the Americas, infecting an estimated 50 to 100 million people every year (World Health Organization Website,, 2006). Primary infections may be asymptomatic, while secondary infections are more severe, with about $5\%$ of secondary infections leading to dengue hemorrhagic fever (DHF) or dengue shock syndrome (DSS), the potentially fatal forms of the disease (World Health Organization Website,, 2006). An effective vaccine against dengue is very difficult to achieve. Because of ADE, infection with an unvaccinated strain following a single strain vaccination may lead to the more severe symptoms associated with secondary infections (Halstead and Deen,, 2002). Therefore, an effective vaccine must protect against all four serotypes simultaneously. A recent theoretical study (Billings et al.,, 2008) has shown that eradication of dengue using only single-strain vaccines is unlikely. Because a tetravalent vaccine is not immediately forthcoming, deepening our understanding of the dynamics of dengue in more realistic models is of great importance.

The dynamics of dengue in a single, well mixed population has been studied in several recent publications, including but not limited to Ferguson et al., 1999a ; Schwartz et al., (2005); Cummings et al., (2005); Wearing and Rohani, (2006); Bianco et al., (2009). ADE and cross immunity have been shown to play a fundamental role in the spreading of dengue, causing instability, desynchronization of serotypes, and chaotic outbreaks (Schwartz et al.,, 2005; Cummings et al.,, 2005; Bianco et al.,, 2009). However, real populations may be spatially heterogeneous. To gain further insight into the wide spectrum of possible dengue dynamics, we relax the assumption of a well mixed population. The division of a population into spatially distinct patches simulates the potentially heterogeneous environment in which the disease spreads. Spatial heterogeneity has been invoked in the past to account for the persistent character of some infectious diseases (Lloyd and May,, 1996; Hagenaars et al.,, 2004) and to explain in-phase and out-of-phase dynamics of diseases (He and Stone,, 2003). A multipatch model with age structure has been used to explain bi- and tri-annual oscillations in the spread of measles (Bolker and Grenfell,, 1995).

Including human mobility patterns has a significant influence on the spreading of infectious disease. We will assume here that the coupling between patches is via migration, movement of individuals from one patch into another (as in Liebovitch and Schwartz, (2004); Ruan et al., (2006); Sattenspiel and Dietz, (1995); Grais et al., (2003)). An alternative coupling strategy is mass action coupling, in which susceptibles in one patch are assumed to interact with infectives in another patch, introducing nonlinear coupling terms (Bolker and Grenfell,, 1995; Lloyd and May,, 1996; Rohani et al.,, 1999). Stochastic coupling has also been modeled (see Keeling and Rohani, (2002) and references therein, and Colizza and Vespignani, (2008)).

The purpose of this work is to analyze the stability of a model for multistrain diseases with interacting strains, using dengue as an example, in a system of two coupled patches. Since chaotic outbreaks are likely to produce a higher number of infected individuals, understanding the stability properties may play an important role in public health. The case of non-identical parameters in the two patches will be of particular interest, as this serves as a model for spatial heterogeneity. The paper is divided as follows. In Section II we introduce the epidemic model. Section III summarizes the bifurcation structure for a single patch. In Section IV we present numerical results for bifurcations in two coupled patches. Section V motivates these results via analysis of a simple, lower dimensional model, and Section VI concludes.

We use a compartmental model for multistrain disease spread with cross immunity and antibody-dependent enhancement, which was previously studied in a single well-mixed population by Bianco et al., (2009). In this model, individuals may develop a primary infection with any of the serotypes. Immediately after recovering, the individual experiences a period of temporary partial cross immunity to all other serotypes. When the cross immunity wears off, immunity to the primary infecting strain is retained, but the individual may develop a secondary infection with a different serotype. Infectiousness of the secondary infectives is increased due to antibody-dependent enhancement. After the secondary infection, complete immunity to all serotypes is assumed. A flow diagram for the single patch model with two serotypes is shown in Figure 1 for simplicity, but we present results here for all four serotypes. We extend the model to two spatially distinct patches, which are coupled by linear migration terms. For two patches (indexed by $q$) and $n$ strains ($n$ arbitrary), the model is as follows:

where the variables are $s_{q}$, the fraction of susceptibles in patch $q$; $x_{q,i}$, the fraction of primary infectives with strain $i$ in patch $q$; $c_{q,i}$, the fraction of individuals in patch $q$ with partial cross immunity to strain $i$; $r_{q,i}$, the fraction of individuals in patch $q$ that are recovered from a primary infection with strain $i$ and no longer have cross immunity to the other strains; and $x_{q,ij}$, the fraction of individuals in patch $q$ that were previously infected with strain $i$ and currently infected with strain $j$. The parameters are the number of strains $n$, the contact rate in patch $q$ $\beta_{q}$, the recovery rate $\sigma$, the ADE factor $\phi$, the strength of cross immunity $\epsilon$, the rate $\theta$ for cross immunity to wear off, the birth and mortality rate in patch $q$ $\mu_{q}$, and the migration rate between patches $\nu$. A list of parameters appears in Table 1.

For simplicity, the birth and death rates in a patch are set equal to each other so that the total population of each patch is constant. The model of Eqs. 1-5 allows for one reinfection. The parameter $\epsilon$ determines how susceptible the cross immune compartments $c_{i}$ are to a secondary infection, where $\epsilon=0$ means no cross immunity (the infection rate is identical for compartments $c_{i}$ and $r_{i}$) and $\epsilon=1$ confers complete cross immunity (cross immunes are immune to a secondary infection for an average time $\theta^{-1}$). We allow $\epsilon$ to take any value between $0$ and $1$. The ADE factor $\phi$ is the enhancement in the infectiousness of secondary infectives. $\phi=1$ means that secondary infectives are as infectious as primary infectives, while $\phi=2$ means that secondary infectives are twice as infectious as primary, and so forth. In contrast to Aguiar and Stollenwerk, (2007), we will not consider values of $\phi$ smaller than $1$. The migration rate from patch $q$ to $q^{\prime}$, and from patch $q^{\prime}$ to $q$, is $\nu$. All individuals migrate with equal probability, independent of their infection status. The migration rate is assumed to be slow compared to the infection spread. For convenience, we put it on the same order as the birth/death rate. We assume that the social parameters, which are the contact rate $\beta_{q}$ and birth/death rate $\mu_{q}$, may vary between patches, while the epidemic parameters are the same in all regions. Because the social parameters depend on human factors (and in the case of the contact rate also include mosquito levels, which are weather-dependent), these parameters are the most likely to be different in adjacent regions. We use parameter values compatible with dengue fever, which are summarized in Table 1. Our contact rate $\beta$ corresponds to a reproductive rate of infection $R_{0}$ of $3.2-4.8$, which is consistent with previous estimates (Ferguson et al., 1999b, ; Nagao and Koelle,, 2008).

A similar model to the one of Eqs. 1-5 has recently been used to analyze the dynamics of dengue fever in a single, well mixed, population (Bianco et al.,, 2009). The ADE $\phi$ and cross immunity strength $\epsilon$ were varied as bifurcation parameters. In the absence of cross immunity ($\epsilon=0$), ADE alone generates instability, desynchronization, and ultimately chaotic outbreaks (Cummings et al.,, 2005; Schwartz et al.,, 2005). A Hopf bifurcation is observed for a critical value of the ADE factor $\phi$, above which oscillatory solutions are obtained. Weak cross immunity stabilizes the system, while strong cross immunity triggers instability and chaos even in the absence of ADE (Bianco et al.,, 2009). In the latter case, destabilization occurs via a Hopf bifurcation for a critical value of the cross immunity strength $\epsilon$. At the bifurcation, three identical complex pairs of eigenvalues of the Jacobian simultaneously become unstable. Although Bianco et al., (2009) discusses the full two parameter bifurcation structure (in $\epsilon$ and $\phi$), we will consider the cross immunity and ADE effects separately in the present work.

Figure 2a shows the bifurcation behavior of the single patch model with no cross immunity as the contact rate $\beta$ is varied. The bifurcation structure was computed using a continuation routine (Doedel et al.,, 1997). The ADE value at which the bifurcation occurs increases slightly as $\beta$ increases, and the dependence is approximately linear. The period of the periodic orbit, shown in Figure 2b for a fixed ADE value, also varies with $\beta$. (Note that multistrain models with ADE can display subcritical Hopf bifurcations (Billings et al.,, 2007), so the periodic orbit shown in Figure 2b exists throughout the range of $\beta$ values shown.) Similarly, varying the contact rate $\beta$ in the absence of ADE affects the location of the Hopf bifurcation in cross immunity $\epsilon$ and its characteristic frequency of oscillation. Likewise, varying the other social parameter, the birth rate $\mu$, affects the location and frequency of the Hopf bifurcations in $\phi$ and $\epsilon$ (data not shown).

We next examine the effect of coupling between distinct regions on the dynamics previously observed in the single patch model. In particular, we consider how migration between nonidentical patches affects the stability of the steady state. We investigate the dynamics of the coupled systems by numerically integrating Eqs. 1-5 and by tracking the bifurcations using a continuation routine (Doedel et al.,, 1997). As mentioned in the previous Section, the model has a Hopf bifurcation at critical values of the bifurcation parameters $\phi$, the ADE factor, and $\epsilon$, the cross immunity strength. We study the effect of coupling and asymmetry on the stability by observing their effect on the location of the Hopf bifurcations. We consider cross immunity and ADE separately; that is, we analyze the system either with ADE and no cross immunity or with cross immunity and no ADE. Including both cross immunity and ADE together leads to qualitatively similar behavior to that reported here, at least locally when the asymmetry is not too large.

We proceed by fixing the patch-specific parameters in patch 2 and varying a social parameter ($\beta_{1}$ or $\mu_{1}$) in patch 1 to observe the effect of increasing asymmetry on the dynamics. As previously mentioned, changing the social parameters modifies the natural frequency of the system. For each value of asymmetry, we look for the critical points at which the coupled system loses stability. We also consider several migration rates $\nu$.

The effect of asymmetry in the contact rate is shown in Figs. 3(a) and 3(b) for two migration rates, namely $\nu=0.02$ (black) and $\nu=0.05$ (gray). The value of the critical parameter $\epsilon$ (Fig. 3(a)) or $\phi$ (Fig. 3(b)) at which the Hopf bifurcation occurs is plotted against the varying contact rate. The symmetric case is at the bottom of the curve. We see that two identical patches have the same bifurcation point as a single, well mixed population even when migration is present (cf. Fig. 2a). However, a striking difference in the dynamics appears when we make the two patches weakly asymmetric. The values of $\phi$ and $\epsilon$ at which the Hopf bifurcation occurs are dramatically different from the symmetric case. The steady state stability persists well into the parameter regime where a single patch would display chaotic oscillations (Billings et al.,, 2007; Bianco et al.,, 2009). The location of the single patch Hopf point does not depend strongly on $\beta$, and decreasing the contact rate is actually destabilizing (Figure 2a), so the stabilization observed here is clearly a result of the coupling between asymmetric systems. Slightly increasing the migration rate (from $\nu=0.02$ to $\nu=0.05$ in Figure 3) further increases the stability of the system.

Similar results are obtained if the asymmetry occurs in the other social parameter, the birth rate. Again, for asymmetric patches either stronger cross immunity or stronger ADE is needed to destabilize the steady state through a Hopf bifurcation (figure in Supplementary Material).

We now motivate the results of the preceding Section using a simple lower dimensional model and show that the qualitative features depend only on the bifurcation structure and characteristic frequencies rather than on details of the epidemic model. In this Section we study the general behavior of two coupled systems, each displaying a Hopf bifurcation, but with different characteristic frequencies.

The generic form for a Hopf bifurcation is the following, in polar coordinates (Strogatz,, 2001):

A Hopf bifurcation occurs at $\lambda=0$, where periodic oscillations with frequency $\omega$ are observed. We now couple two systems of the sort in Eq. 6, each with a potentially different frequency $\omega_{q}$, via linear migration terms with migration rate $\nu$. To lowest order (not displaying cubic terms), the coupled system in cartesian coordinates is

Stability is determined by evaluating the Jacobian of Eq. 7 at the steady state $(x_{1},y_{1},x_{2},y_{2})=0$. At the Hopf bifurcation, the real part of the largest eigenvalue crosses zero, with nonzero imaginary part. We obtain the Hopf point directly by setting the real part of the eigenvalue to zero and solving for the critical value of $\lambda$ at the bifurcation. Details may be found in the Supplementary Material. The Hopf bifurcation occurs at

In Fig. 4 we show the location of the Hopf bifurcation given by Eq. 8 as a function of the asymmetry between the two systems for $\nu=0.05$. Comparison with Fig. 3 (and Fig. 1 of the Supplementary Material) shows the qualitative agreement of the theoretical results for the lower dimensional system with the full multistrain system. Increasing the migration rate $\nu$ increases the value of the bifurcation parameter to which the Hopf point saturates (not shown), as in Fig. 3. It is also worth noticing that, in the case of identical frequencies $\omega_{1}=\omega_{2}$, the Hopf bifurcation occurs at $\lambda_{c}=0$, the location in the absence of coupling. This is consistent with the numerical results for the multistrain system in the case of symmetric patches.

We have studied the endemic steady state stability properties for a multistrain epidemic model on two migration-coupled patches. Interactions between strains in the model were governed by temporary partial cross immunity and antibody-dependent enhancement. In the absence of coupling, the system displayed Hopf bifurcations in two epidemic parameters. Coupling between patches with non-identical parameters, which gave them non-identical characteristic frequencies of oscillation, was shown to shift the Hopf bifurcations, stabilizing the steady state. This behavior was observed for the Hopf bifurcation obtained by sweeping the cross immunity in the absence of ADE and the bifurcation obtained by sweeping the ADE in the absence of cross immunity. It occurred for asymmetry in either of our two social parameters, the birth rate and the contact rate.

To motivate this result, we diffusively coupled two low dimensional Hopf bifurcations with different characteristic frequencies and analyzed the stability of the steady state. We again saw that coupling between asymmetric systems led to stabilization. This indicates that the stabilization in the epidemic model is a result of the underlying dynamics, rather than the details of the model. The stabilization may occur because when the frequencies of oscillation in two coupled systems are different, the oscillations tend to cancel each other because of phase differences. This topic will be studied in more detail in a future work.

Bifurcations from steady state to oscillatory behavior can be associated with an increased number of infection cases, particularly if chaotic oscillations occur, as in previous dengue models (Billings et al.,, 2007; Bianco et al.,, 2009). Our results suggest that if control strategies in one region are able to generate enough asymmetry, this could lead to a stabilization of the outbreaks, which would have a positive effect on adjacent regions. Asymmetry could be generated in the effective contact rate by mosquito control, which could include reducing mosquito breeding sites (Slosek,, 1986), or through new genetic controls which are under development (Barbazan et al.,, 2008). Asymmetry in the birth rate could be generated by lowering the effective birth rate through vaccination of new susceptibles once a vaccine is available. However, because the bifurcation point saturates rather than increasing indefinitely as asymmetry increases, such a strategy would be successful only if the epidemic parameters in the real system are moderately close to the bifurcation. (Extreme asymmetry may even be destabilizing, so this strategy works best locally when the asymmetry is not too strong.) In addition, the role of seasonality in exciting oscillations should not be ignored. Seasonal variations in the contact rate have been included in previous dengue models (Schwartz et al.,, 2005). The interplay of seasonality and coupling is a topic for future study.

The results of the present work may also be useful in lower bound estimation of the parameter values for ADE and cross immunity. ADE especially is difficult to measure in vivo and must be estimated by other means. Since recent publications have suggested that epidemiological data from Thailand show chaotic outbreaks (Schwartz et al.,, 2005), and given the certainty of human migration and asymmetry between adjacent regions of the country, it is possible that the actual parameter values for ADE and cross immunity are higher than the ones estimated by studying a single well mixed population.

Finally, the work discussed here shows a potential effect of human movement between heterogeneous regions. As spatial effects are further studied in epidemic models, it remains to be seen how this phenomenon will extend to more complicated spatial geometries, including more patches and perhaps non-symmetric coupling terms. Because the migration-induced stabilization depends only on the existence of a Hopf bifurcation in the model, it is expected that the stabilization will be observed in other population models that also contain Hopf bifurcations (e.g., Fussmann et al., (2000); Greenhalgh et al., (2004)).

This work was supported in part by the Jeffress Memorial Trust. The authors wish to thank Ira Schwartz for helpful discussions.