{'fraction_non_alphanumeric': 0.03200930052315443, 'fraction_numerical': 0.01747723309436156, 'mean_word_length': 4.603908794788274, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymmetry in the presence of migration stabilizes multistrain disease outbreaks

We study the effect of migration between coupled populations, or patches, on the stability properties of multistrain disease dynamics. The epidemic model used in this work displays a Hopf bifurcation to oscillations in a single well mixed population. It is shown numerically that migration between two non-identical patches stabilizes the endemic steady state, delaying the onset of large amplitude outbreaks and reducing the total number of infections. This result is motivated by analyzing generic Hopf bifurcations with different frequencies and with diffusive coupling between them. Stabilization of the steady state is again seen, indicating that our observation in the full multistrain model is based on qualitative characteristics of the dynamics rather than on details of the disease model.


Introduction
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. 1 However, there is evidence that antibodies also give some cross-protection to the other serotypes (Halstead 2007). This reduced susceptibility to the other serotypes is temporary. When the temporary cross immunity wanes, heterologous secondary infections are possible. Low level antibodies developed from primary infections are believed to form complexes with the virus so that more cells are infected, 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 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 are sometimes 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 could 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 on a 2 strains system (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, Wearing and Rohani (2006), Bianco et al. (2009). ADE and cross immunity have been shown to play a fundamental role in mathematical models for the spreading of dengue, causing instability, desynchronization of serotypes, and chaotic outbreaks 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).
Human mobility patterns have 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 can be modeled (see Keeling andRohani 2002 andreferences therein, andColizza andVespignani 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 Sect. 2 we introduce the epidemic model. Section 3 summarizes the bifurcation structure for a single patch. In Sect. 4 we present numerical results for bifurcations in two coupled patches. Section 5 motivates these results via analysis of a simple, lower dimensional model, and Sect. 6 concludes.

Model
We use a compartmental model for multistrain disease spread with cross immunity and antibody-dependent enhancement, previously studied in a single, well-mixed population (Bianco et al. 2009). In this model, individuals can 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 can 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 Fig. 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 recovered from strain i and currently infected with strain j . The parameters are the number of strains n, the contact rate in patch q β q , the recovery rate σ , the ADE factor φ, the strength of cross immunity , the rate θ for cross immunity to wear off, the birth and mortality rate in patch q μ q , and the migration rate between patches ν. A list of parameters appears in Table 1.
For clarity, we describe in detail the terms appearing in (5) describing the time evolution of the secondary infectives x q,ij . These are individuals who were previously infected with strain i and are currently infected with strain j . Individuals recovered from an infection with i and no longer retaining cross immunity (r q,i ) can become infected by any infective with strain j . We use mass action contact between infectives and non-invectives, contributing a term β q r q,i (x q,j + φ k =j x q,kj ) (where the secondary infectives have enhanced contagion given by the ADE factor φ). Cross immune individuals c q,i become infected with a reduced probability (1 − ). Linear terms −σ x q,ij and −μ q x q,ij appear as secondary infectives recover with rate σ or are lost to natural death with rate μ q . Finally, individuals are exchanged between patches q and q with rate ν, yielding migration terms −νx q,ij + νx q ,ij .
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 (1)-(5) allows for one reinfection. Tertiary infections are not considered (Nisalak et al. 2003). The parameter determines how susceptible the cross immune compartments c i are to  (Bianco et al. 2009). Note the reduction of susceptibility to a secondary infection through the cross immunity factor (1 − ) and the enhancement of secondary infectiousness due to the ADE factor φ. Mortality terms for each compartment are not included in the diagram for ease of reading  Ferguson et al. (1999b) β, transmission coefficient, years −1 ∼200 Ferguson et al. (1999b) σ , recovery rate, years −1 50 Rigau-Perez et al. (1998) θ, rate to leave the cross 2 Wearing and Rohani (2006) immune compartment, years −1 φ, ADE factor ≥1 S c h w a r t z e t a l . ( 2005) , strength of cross immunity 0-1 Bianco et al. (2009) ν, migration rate, years −1 0-0.05 n, number of strains 4 a secondary infection, where = 0 means no cross immunity (the infection rate is identical for compartments c i and r i ) and = 1 confers complete cross immunity (cross immunes are immune to a secondary infection for an average time θ −1 ). We allow to take any value between 0 and 1. The ADE factor φ is the enhancement in the infectiousness of secondary infectives. φ = 1 means that secondary infectives are as infectious as primary infectives, while φ = 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 φ smaller than 1. The migration rate from patch q to q , and from patch q to q, is ν. For simplicity, we assume that all individuals migrate with equal probability, independent of their infection status. This assumption may be relaxed in a future study. 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 β q and birth/death rate μ 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 β 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).

Single Patch Bifurcation Structure
A similar model to the one of (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 φ and cross immunity strength were varied as bifurcation parameters. In the absence of cross immunity ( = 0), ADE alone generates instability, desynchronization, and ultimately chaotic outbreaks Schwartz et al. 2005). A Hopf bifurcation is observed for a critical value of the ADE factor φ, above which oscillatory solutions are obtained. Weak cross immunity stabilizes the system, while strong cross immunity triggers instability and chaos even in  Table 1 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 . 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 and φ), 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 β 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 β increases, and the dependence is approximately linear. The period of the periodic orbit, shown in Fig. 2b for a fixed ADE value, also varies with β. (Note that multistrain models with ADE can display subcritical Hopf bifurcations, Billings et al. 2007, so the periodic orbit shown in Fig. 2b exists throughout the range of β values shown.) Similarly, varying the contact rate β in the absence of ADE affects the location of the Hopf bifurcation in cross immunity and its characteristic frequency of oscillation. Likewise, varying the other social parameter, the birth rate μ, affects the location and frequency of the Hopf bifurcations in φ and (data not shown).

Coupled Patches
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 non-identical patches affects the stability of the steady state. We investigate the dynamics of the coupled systems by numerically integrating (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 φ, the ADE factor, and , the cross immunity strength. We Fig. 3 Critical values of the parameters (a) and (b) φ at which the Hopf bifurcation occurs for coupled patches, as a function of the contact rate in patch 1, for two values of the migration rate, ν = 0.02 (solid black) and ν = 0.05 (dashed gray). The contact rate in patch 2 is fixed at β 2 = 200. In (a), φ = 1 (no ADE). In (b), = 0 (no cross immunity). μ 1 = μ 2 = 0.02 and other parameters are as in Table 1 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 (β 1 or μ 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 ν.
The effect of asymmetry in the contact rate is shown in Figs. 3(a) and 3(b) for two migration rates, namely ν = 0.02 (black) and ν = 0.05 (gray). The value of the critical parameter ( Fig. 3(a)) or φ (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 φ and 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 β, and decreasing the contact rate is actually destabilizing (Fig. 2a), so the stabilization observed here is clearly a result of the coupling between asymmetric systems. Slightly increasing the migration rate (from ν = 0.02 to ν = 0.05 in Fig. 3) further increases the stability of the system. Fig. 4 Critical values of the parameters (a) and (b) φ at which the Hopf bifurcation occurs for coupled patches, as a function of the birth/death rate in patch 1, for ν = 0.02. The birth/death rate in patch 2 is fixed at μ 2 = 0.02. In (a), φ = 1 (no ADE). In (b), = 0 (no cross immunity). β 1 = β 2 = 200 and other parameters are as in Table 1 Similar results are obtained if the asymmetry occurs in the other social parameter, the birth rate, as depicted in Fig. 4. Again, for asymmetric patches either stronger cross immunity or stronger ADE is needed to destabilize the steady state through a Hopf bifurcation.

Analytical Motivation
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 α = 0, where periodic oscillations with frequency ω are observed. We now couple two systems of the sort in (6), each with a potentially different frequency ω q , via linear migration terms with migration rate ν. To lowest order (not displaying cubic terms), the coupled system in cartesian coordinates iṡ x 1 = αx 1 − ω 1 y 1 − νx 1 + νx 2 y 1 = αy 1 + ω 1 x 1 − νy 1 + νy 2 x 2 = αx 2 − ω 2 y 2 − νx 2 + νx 1 y 2 = αy 2 + ω 2 x 2 − νy 2 + νy 1 Stability is determined by evaluating the Jacobian of (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. The roots of the characteristic polynomial f (λ) of the Jacobian are the eigenvalues {λ}. The four eigenvalues are but it is not easy to see from these expressions where the Hopf bifurcation occurs. Instead, we solve directly for the critical value of α at the Hopf point. At a Hopf bifurcation, a pair of eigenvalues has zero real part. Thus to obtain the Hopf point, we can set λ = ib, where b is real. If ib is an eigenvalue, we require and If α = ν, (11) has nonzero roots b = ± 1 2 4α 2 − 8αν + 2ω 2 1 + 2ω 2 2 . When this b is substituted into (10), we obtain four potential roots for α at the Hopf point. Two are complex, which is unphysical and can be ignored. The one corresponding to the Hopf bifurcation is α c = ν − 1 2 4ν 2 − (ω 1 − ω 2 ) 2 . (The fourth root is larger and corresponds to loss of stability of the second pair of eigenvalues.) This α c is real and thus gives the Hopf point location when |ω 1 − ω 2 | ≤ 2ν. When ω 1 − ω 2 = 0 and ν > 0, α c is positive, indicating stabilization due to the asymmetry and coupling.
Summarizing, the Hopf bifurcation occurs at In Fig. 5 we show the location of the Hopf bifurcation given by (12) as a function of the asymmetry between the two systems for ν = 0.05. Comparison with Figs. 3 and 4 shows the qualitative agreement of the theoretical results for the lower dimensional system with the full multistrain system. Increasing the migration rate ν increases the value of the bifurcation parameter to which the Hopf point saturates, as in Fig. 3. It is also worth noticing that, in the case of identical frequencies ω 1 = ω 2 , the Hopf bifurcation occurs at α 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.

Conclusions and Discussion
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. We suggest that the stabilization may occur as a result of the two different coupled frequencies generating oscillations that 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.) Real systems exhibiting oscillatory behavior may already be well into an unstable region of parameters. For such systems, our suggested control strategy may not apply since asymmetry between adjacent regions is already likely, thus reducing that impact that an increase in the asymmetry could have. On the other hand, in the case of dengue fever a recent study (Cummings et al. 2009) has highlighted a decreasing trend in both the force of infection and birth and death rates in the population of Thailand during a time span of about 30 years. This phenomenon could potentially drive the system closer to the bifurcation point and make it more amenable to control.
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 . The interplay of seasonality and coupling is a topic for future study.
The results of the present work may also be useful in 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 , 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. This work represents a first step towards understanding the role of migration and spatial heterogeneity in dynamical properties of dengue observed in epidemiological data, such as traveling waves of infection in Thailand (Cummings et al. 2004). Furthermore, 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).