Torus bifurcations, isolas and chaotic attractors in a simple dengue model with ADE and temporary cross immunity
We analyse an epidemiological model of competing strains of pathogens and hence differences in transmission for first versus secondary infection due to interaction of the strains with previously aquired immunities, as has been described for dengue fever (in dengue known as antibody dependent enhancement, ADE). Such models show a rich variety of dynamics through bifurcations up to deterministic chaos. Including temporary cross-immunity even enlarges the parameter range of such chaotic attractors, and also gives rise to various coexisting attractors, which are difficult to identify by standard numerical bifurcation programs using continuation methods. A combination of techniques, including classical bifurcation plots and Lyapunov exponent spectra has to be applied in comparison to get further insight into such dynamical structures. Here we present for the first time multi-parameter studies in a range of biologically plausible values for dengue. The multi-strain interaction with the immune system is expected to also have implications for the epidemiology of other diseases.
umerical bifurcation analysis, Lyapunov exponents, symmetry, coexisting attractors, antibody dependent enhancement (ADE)
Centro de Matemática e Aplicações Fundamentais Universidade de Lisboa, Portugal
Faculty of Earth and Life Sciences, Department of Theoretical Biology Vrije Universiteit Amsterdam, Nederland
Epidemic models are classically phrased in ordinary differential equation (ODE) systems for the host population divided in classes of susceptible individuals and infected ones (SIS system), or in addition, a class of recovered individuals due to immunity after an infection to the respective pathogen (SIR epidemics). The infection term includes a product of two variables, hence a non-linearity which in extended systems can cause complicated dynamics. Though these simple SIS and SIR models only show fixed points as equilibrium solutions, they already show non-trivial equilibria arising from bifurcations, and in stochastic versions of the system critical fluctuations at the threshold. Further refinements of the SIR model in terms of external forcing or distinction of infections with different strains of a pathogen, hence classes of infected with one or another strain recovered from one or another strain, infected with more than one strain etc., can induce more complicated dynamical attractors including equilibria, limit cycles, tori and chaotic attractors.
Classical examples of chaos in epidemiological models are childhood diseases with extremely high infection rates, so that a moderate seasonal forcing can generate Feigenbaum sequences of period doubling bifurcations into chaos. The success in analysing childhood diseases in terms of modelling and data comparison lies in the fact that they are just childhood diseases with such high infectivity. Otherwise host populations cannot sustain the respective pathogens. In other infectious diseases much lower forces of infection have to be considered leading to further conceptual problems with noise affecting the system more than the deterministic part, leading even to critical fluctuations with power law behaviour, when considering evolutionary processes of harmless strains of pathogens versus occasional accidents of pathogenic mutants . Only explicitly stochastic models, of which the classical ODE models are mean field versions, can capture the fluctuations observed in time series data .
More recently it has been demonstrated that the interaction of various strains on the infection of the host with eventual cross-immunities or other interactions between host immune system and multiple strains can generate complicated dynamic attractors. A prime example is dengue fever. A first infection is often mild or even asymptomatic and leads to life long immunity against this strain. However, a subsequent infection with another strain of the virus often causes clinical complications up to life threatening conditions and hospitalization, due to ADE. More on the biology of dengue and its consequences for the detailed epidemiological model structure can be found in Aguiar and Stollenwerk  including literature on previous modelling attempts, see also . On the biological evidence for ADE see e.g. . Besides the difference in the force of infection between primary and secondary infection, parametrized by a so called ADE parameter , which has been demonstrated to show chaotic attractors in a certain parameter region, another effect, the temporary cross-immunity after a first infection against all dengue virus strains, parametrized by the temporary cross-immunity rate , shows bifurcations up to chaotic attractors in a much wider and biologically more realistic parameter region. The model presented in the Appendix has been described in detail in  and has recently been analysed for a parameter value of corresponding to on average half a year of temporary cross immunity which is biologically plausible . For increasing ADE parameter first an equilibrium which bifurcates via a Hopf bifurcation into a stable limit cycle and then after further continuation the limit cycle becomes unstable in a torus bifurcation. This torus bifurcation can be located using numerical bifurcation software based on continuation methods tracking known equilibria or limit cycles up to bifurcation points . The continuation techniques and the theory behind it are described e.g. in Kuznetsov . Complementary methods like Lyapunov exponent spectra can also characterize chaotic attractor [9, 10], and led ultimately to the detection of coexisting attractors to the main limit cycles and tori originated from the analytically accessible fixed point for small . Such coexisting structures are often missed in bifurcation analysis of higher dimensional dynamical systems but are demonstrated to be crucial at times in understanding qualitatively the real world data, as for example demonstrated previously in a childhood disease study . In such a study first the understanding of the deterministic system’s attractor structure is needed, and then eventually the interplay between attractors mediated by population noise in the stochastic version of the system gives the full understanding of the data. Here we present for the first time extended results of the bifurcation structure for various parameter values of the temporary cross immunity in the region of biological relevance and multi-parameter bifurcation analysis. This reveals besides the torus bifurcation route to chaos also the classical Feigenbaum period doubling sequence and the origin of so called isola solutions. The symmetry of the different strains leads to symmerty breaking bifurcations of limit cycles, which are rarely described in the epidemiological literature but well known in the biochemical literature, e.g for coupled identical cells. The interplay between different numerical procedures and basic analytic insight in terms of symmetries help to understand the attractor structure of multi-strain interactions in the present case of dengue fever, and will contribute to the final understanding of dengue epidemiology including the observed fluctuations in real world data. In the literature the multi-strain interaction leading to deterministic chaos via ADE has been described previously, e.g. [12, 13] but neglecting temporary cross immunity and hence getting stuck in rather unbiological parameter regions, whereas more recently the first considerations of temporary cross immunity in rather complicated and up to now not in detail analysed models including all kinds of interations have appeared [14, 15], in this case failing to investigate closer the possible dynamical structures.
2 Dynamical system
The multistrain model under investigation can be given as an ODE system
for the state vector of the epidemiological host classes and besides other fixed parameters which are biologically undisputed the parameter vector of varied parameters . For a detailed description of the biological content of state variables and parameters see . The ODE equations and fixed parameter values are given in the appendix. The equilibrium values are given by the equilibrium condition , respectively for limit cycles with period . For chaotic attractors the trajectory of the dynamical system reaches in the time limit of infinity the attractor trajectory , equally for tori with irrational winding ratios. In all cases the stability can be analysed considering small perturbations around the attractor trajectories
Here, any attractor is notified by , be it an equilibrium, periodic orbit or chaotic attractor. In this ODE system the linearized dynamics is given with the Jacobian matrix of the ODE system Eq. (1) evaluated at the trajectory points given in notation of . The Jacobian matrix is analyzed for equilibria in terms of eigenvalues to determine stability and the loss of it at bifurcation points, negative real part indicating stability. For the stability and loss of it for limit cylces Floquet multipliers are more common (essentially the exponentials of eigenvalues), multipliers inside the unit circle indicating stability, and where they leave eventually the unit circle determining the type of limit cycle bifurcations. And for chaotic systems Lyapunov exponents are determined from the Jacobian around the trajectory, positive largest exponents showing deterministic chaos, zero largest showing limit cycles including tori, largest smaller zero indicating fixed points.
To investigate the bifurcation structure of the system under investigation
we first observe the symmetries due to the multi-strain structure of the
model. This becomes important for the time being for
we have the following symmetry:
with equilibrium values or limit cycle for all times . For the right hand side of the ODE system Eq. (1) the kind of symmetry found above is called -symmetry when the following equivariance condition holds
with a matrix that obeys and , where is the unit matrix. Observe that besides also satisfies (5). The symmetry transformation matrix in Eq. (3) fulfills these requirements. It is easy to verify that the -equivariance conditions Eq. (5) and the properties of are satisfied for our ODE system. In Seydel  a simplified version of the famous Brusselator that shows this type of symmetry is discussed. There, an equilibrium and also a limit cycle show a pitchfork bifurcation with symmetry breaking.
An equilibrium is called fixed when (see ). Two equilibria where , are called -conjugate if their corresponding solutions satisfy (and because also ). For limit cycles a similar terminology is introduced. A periodic solution is called fixed when and the associated limit cycles are also called fixed . There is another type of periodic solution that is not fixed but called symmetric when
where is the period. Again the associated limit cycles are also called symmetric. Both types of limit cycles are -invariant as curves : . That is, in the phase-plane where time parameterizes the orbit, the cycle and the transformed cycle are equal. A -invariant cycle is either fixed or symmetric. Two noninvariant limit cycles () are called -conjugate if their corresponding periodic solutions satisfy . The properties of the symmetric systems and the introduced terminology are used below with the interpretation of the numerical bifurcation analysis results. We refer to  for an overview of the possible bifurcations of equilibria and limit cycles of -equivariant systems.
3 Bifurcation diagrams for various values
We show the results of the bifurcation analysis in bifurcation diagrams for several values, varying continuously. Besides the previously investigated case of , we show also a case of smaller and a case of larger value, obtaining more information on the bifurcations possible in the model as a whole. The above mentioned symmetries help in understanding the present bifurcation structure.
3.1 Bifurcation diagram for
For the one-parameter bifurcation diagram is shown in Fig. 1 a). Starting with there is a stable fixed equilibrium, fixed in the above mentioned notion for symmetric systems. This equilibrium becomes unstable at a Hopf bifurcation at . A stable fixed limit cycle originates at this Hopf bifurcation. This limit cycle shows a supercritical pitch-fork bifurcation , i.e. a bifurcation of a limit cycle with Floquet multiplier 1, splitting the original limit cycle into two new ones. Besides the now unstable branch two new branches originate for the pair of conjugated limit cycles. The branches merge again at another supercritical pitch-fork bifurcation , after which the limit cycle is stable again for higher -values. The pair of -conjugate limit cycles become unstable at a torus bifurcation at .
Besides this main bifurcation pattern we found two isolas, that is an isolated solution branch of limit cycles . These isola cycles are not -invariant, that is . Isolas consisting of isolated limit cycles exist between two tangent bifurcations. One isola consists of a stable and an unstable branch. The other shows more complex bifurcation patterns. There is no full stable branch. For at the tangent bifurcation a stable and an unstable limit cycle collide. The stable branch becomes unstable via a flip bifurcation or periodic doubling bifurcation , with Floquet multiplier , at which is also pitchfork bifurcation for the period-two limit cycles. At the other end of that branch at the tangent bifurcation at both colliding limit cycles are unstable. Close to this point at one branch there is a torus bifurcation , also called Neimark-Sacker bifurcation, at and a flip bifurcation at which is again a pitchfork bifurcation for the period-two limit cycles. Contiuation of the stable branch originating for the flip bifurcation at gives another flip bifurcation at and one closed to the other end at , namely at . These results suggest that for this isola two classical routes to chaos can exist, namely via the torus or Neimark-Sacker bifurcation where the dynamics on the originating torus is chaotic, and the cascade of period doubling route to chaos.
3.2 Bifurcation diagram for
For the one-parameter bifurcation diagram is shown in Fig. 1 b). The stable fixed equilibrium becomes unstable at a supercritical Hopf bifurcation at where a stable fixed limit cycle originates. This stable limit cycle becomes unstable at a superciritcal pitchfork bifurcation point at for a limit cycle. This point marks the origin of a pair of -conjugate stable limit cycles besides the now unstable fixed limit cycle. Here one has to consider the two infected subpopulations and to distinguish the conjugate limit cycles. Because the two variables and are interchangeable this can also be interpreted as the stable limit cycles for the single variable say . The fixed stable equilibrium below the Hopf bifurcation where we have , , and is a fixed equilibrium. For the fixed limit cycle in the parameter interval between the Hopf bifurcation and the pitchfork bifurcation we have , , and . This means that at the Hopf bifurcation the stable fixed equilibrium becomes an unstable fixed equilibrium. In the parameter interval between the two pitchfork bifurcations at and subcritical at , two stable limit cycles coexist and these limit cycles are -conjugate. At the pitchfork bifurcation points the fixed limit cycle becomes unstable and remains fixed, and two stable -conjugate limit cycles originate (see [8, Theorem 7.7]). The invariant plane forms the separatrix between the pair of stable -conjugate limit cycles and . The initial values of the two state variables and together with the point on the invariant plane, determine to which limit cycle the system converges. Continuation of the stable symmetric limit cycle gives a torus or Neimark-Sacker bifurcation at point denoted by at . At his point the limit cycles become unstable because a pair of complex-conjugate multipliers crosses the unit circle. Observe that at this point in the time series plot [3, there Fig. 12] the chaotic region starts. In  the following route to chaos, namely the sequence of Neimark-Sacker bifurcations into chaos, is mentioned. Increasing the bifurcation parameter along the now unstable pair of -conjugate limit cycles leads to a tangent bifurcation at where a pair of two unstable limit cycles collide. This branch terminates at the second pitchfork bifurcation point denoted by at . Because the first fold point gave rise to a stable limit cycle and this fold point to an unstable limit cycle we call the first pitchfork bifurcation supercritical and the latter pitchfork bifurcation subcritical. These results agree very well with the simulation results shown in the bifurcation diagram for the maxima and minima of the overall infected [3, there Fig. 15]. Notice that AUTO  calculates only the global extrema during a cycle, not the local extrema. Fig. 1 b) shows also two isolas similar to those for in Fig. 1 a).
3.3 Bifurcation diagram for
For the bifurcation diagram is shown in Fig 1 c). In the lower parameter range there is bistability of two limit cycles in an interval bounded by two tangent bifurcations . The stable manifold of the intermediate saddle limit cycle acts as a separatrix. Inceasing the stable limit cycles become unstable at the pitchfork bifurcation at . Following the unstable primary branch, for larger values of we observe an open loop bounded by two tangent bifurcations . The extreme value for is at . Then lowering there is a pitchfork bifurcation at . Later we will return to the description of this point. Lowering further the limit cycle becomes stable again at the tangent bifurcations at . Increasing this limit cycle becomes unstable again at the pitchfork bifurcation at .
Continuation of the secondary branch of the two -conjugated limit cycles from this point reveals that the stable limit cycle becomes unstable at a torus bifurcation at . The simulation results depicted in [3, Fig. 13] show that there is chaos beyond this point. The secondary pair of -conjugate limit cycles that originate from pitchfork bifurcation at becomes unstable at a flip bifurcation . Increasing further it becomes stable again at a flip bifurcation . Below we return to the interval between these two flip bifurcations. The stable part becomes unstable at a tangent bifurcation , then continuing, after a tangent bifurcation and a Neimark-Sacker bifurcation . This bifurcation can lead to a sequence of Neimark-Sacker bifurcations into chaos. The unstable limit cycles terminates via a tangent bifurcation where the primary limit cycle possesses a pitchfork bifurcation at . At the flip bifurcation the cycle becomes unstable and a new stable limit cycle with double period emanates. The stable branch becomes unstable at a flip bifurcation again. We conclude that there is a cascade of period doubling route to chaos. Similarly this happens in reversed order ending at the flip bifurcation where the secondary branch becomes stable again.
Fig. 2 a) gives the results for the interval where only the minima are show. In this plot also a “period three” limit cycle is shown. In a small region it is stable and coexists together with the “period one” limit cycle. The cycles are shown in Fig. 2 b) and c) for . The one in c) looks like a period-3 limit cycle. In Fig. 2 continuation of the limit cycle gives a closed graph bounded at the two ends by trangent bifurcations where a stable and an unstable limit cycle collide. The intervals where the limit cycle is stable, are on the other end bounded by flip bifurcations . One unstable part intersects the higher period cycles that originate via the cascade of period doubling between the period-1 limit cycle flip bifurcations at and . This suggest that the period-3 limit cycle is associated with a “period-3 window” of the chaotic attractor. We conjecture that this interval is bounded by two homoclinic bifurcations for a period-3 limit cycle (see [19, 20, 21, 22]). The bifurcation diagram shown in [3, there Fig. 13] shows the point where the chaotic attractor disappears abruptly, possible at one of the two homoclinic bifurcations. In that region the two conjugated limit cycles that originate at the pitchfork bifurcation at are the attractors. These results suggest that there are chaotic attractors associated with the period-1 limit cycle, one occurs via a cascade of flip bifurcations originating from the two ends at and and one via a Neimark-Sacker bifurcation at .
4 Two-parameter diagram
We will now link the three studies of the different values by investigating a two-parameter diagram for and , concentrating especially on the creation of isolated limit cycles, which sometimes lead to further bifurcations inside the isola region. Fig. 3 gives a two-parameter bifurcation diagram where and are the free parameters. For low -values there is the Hopf bifurcation and all other curves are tangent bifurcation curves.
Isolas appear or disappears upon crossing an isola variety. At an elliptic isola point an isolated solution branch is born, while at a hyperbolic isola point an isolated solution branch vanishes by coalescence with another branch . From Fig. 3 we see that at two values of isolas are born. Furthermore, period doubling bifurcations appear for lower values, indicating the Feigenbaum route to chaos. However, only the calculation of Lyapunov exponents, which are discussed in the next section, can clearly indicate chaos.
5 Lyapunov spectra for various values
The Lyapunov exponents are the logarithms of the eigenvalues of the Jacobian matrix along the integrated trajectories, Eq. (2), in the limit of large integration times. Besides for very simple iterated maps no analytic expressions for chaotic systems can be given for the Lyapunov exponents. For the calculation of the iterated Jacobian matrix and its eigenvalues, we use the QR decomposition algorithm .
In Fig. 4 we show for various values the four largest Lyapunov exponents in the range between zero and one. For in Fig. 4 a) we see for small values fixed point behaviour indicated by a negative largest Lyapunov exponent up to around . There, at the Hopf bifurcation point, the largest Lyapunov exponent becomes zero, indicating limit cycle behaviour for the whole range of , apart from the final bit before , where a small spike with positive Lyapunov exponent might be present, but difficult to distinguish from the noisy numerical background.
For in Fig. 4 b) however, we see a large window with positive largest Lyapunov exponent, well separated from the second largest being zero. This is s clear sign of deterministically chaotic attractors present for this range. Just a few windows with periodic attractors, indicated by the zero largest Lyapunov exponent are visible in the region of . For smaller values we observe qualitatively the same behaviour as already seen for . For the smaller value of in Fig. 4 c) the chaotic window is even larger than for . Hence deterministic chaos is present for temporary cross immunity in the range around in the range of between zero and one.
We have presented a detailed bifurcation analysis for a multi-strain dengue fever model in terms of the ADE parameter , in the previously not well investigated region between zero and one, and a parameter for the temporary cross immunity . The symmetries implied by the strain structure, are taken into account in the analysis. Many of the possible bifurcations of equilibria and limit cycles of -equivariant systems can be distinguished. Using AUTO  the different dynamical structures were calculated. Future time series analysis of epidemiological data has good chances to give insight into the relevant parameter values purely on topological information of the dynamics, rather than classical parameter estimation of which application is in general restricted to farely simple dynamical scenarios.
This work has been supported by the European Union under the Marie Curie grant MEXT-CT-2004-14338. We thank Gabriela Gomes and Luis Sanchez, Lisbon, for scientific support.
Appendix A Appendix: Epidemic model equations
The complete system of ordinary differential equations for a two strain epidemiological system allowing for differences in primary versus secondary infection and temporary cross immunity is given by
For two different strains, and , we label the SIR classes for the hosts that have seen the individual strains. Susceptibles to both strains (S) get infected with strain () or strain (), with infection rate . They recover from infection with strain (becoming temporary cross-immune ) or from strain (becoming ), with recovery rate etc.. With rate , the and enter again in the susceptible classes ( being immune against strain 1 but susceptible to 2, respectively ), where the index represents the first infection strain. Now, can be reinfected with strain (becoming ), meeting with infection rate or meeting with infection rate , secondary infected contributing differently to the force of infection than primary infected, etc..
We include demography of the host population denoting the birth and death rate by . For constant population size we have for the immune to all strains and therefore we only need to consider the first 9 equations of Eq. (7), giving 9 Lyapunov exponents. In our numerical studies we take the population size equal to so that numbers of susceptibles, infected etc. are given in percentage. As fixed parameter values we take , , . The parameters and are varied.
- Equilibria are often called fixed points in dynamical systems theory, here we try to avoid this term, since in symmetry the term fixed is used in a more specific way, see below.
- N. Stollenwerk and V.A.A. Jansen, Evolution towards criticality in an epidemiological model for meningococcal disease, Physics Letters A 317 (2003) 87–96.
- N. Stollenwerk, M.C.J. Maiden and V.A.A. Jansen, Diversity in pathogenicity can cause outbreaks of menigococcal disease, Proc. Natl. Acad. Sci. USA 101 (2004) 10229–10234.
- M. Aguiar and N. Stollenwerk, A new chaotic attractor in a basic multi-strain epidemiological model with temporary cross-immunity, arXiv:0704.3174v1 [nlin.CD] (2007) (accessible electronically at http://arxive.org).
- E. Massad, M. Chen, S. Ma, C. J. Struchiner, N. Stollenwerk, and M. Aguiar, Scale-free Network for a Dengue Epidemic, Applied Mathematics and Computation 159 (2008) 376–381.
- S.B. Halstead, Neutralization and antibody-dependent enhancement of dengue viruses, Advances in Virus Research 60 (2003) 421–67.
- M. Aguiar, B.W. Kooi, and N. Stollenwerk, Epidemiology of dengue fever: A model with temporary cross-immunity and possible secondary infection shows bifurcations and chaotic behaviour in wide parameter regions, submitted (2008).
- E.J. Doedel, R.C. Paffenroth, A.R. Champneys, T.F. Fairgrieve, Y.A. Kusnetsov, B. Sandstede, B. Oldeman, X.J. Wang, and C. Zhang, AUTO 07P – Continuation and bifurcation software for ordinary differential equations, Technical Report: Concordia University, Montreal, Canada (2007) (accessible electronically at http://indy.cs.concordia.ca/auto/).
- Y.A. Kuznetsov, Elements of Applied Bifurcation Theory Applied Mathematical Sciences 112, Springer-Verlag, 3 edition, New York, 2004.
- D. Ruelle,Chaotic Evolution and Strange Attractors, Cambridge University Press, Cambridge, 1989.
- E. Ott, Chaos in Dynamical Systems, Cambridge University Press, Cambridge, 2002.
- F.R. Drepper, R. Engbert and N. Stollenwerk, Nonlinear time series analysis of empirical population dynamics, Ecological Modelling 75/76 (1994) 171–181.
- N. Ferguson, R. Anderson, and S. Grupta, The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens, Proc. Natl. Acad. Sci. USA 96 (1999) 790–94.
- L. Billings, B.I. Schwartz, B.L. Shaw, M. McCrary, D.S. Burke and T.A.D. Cummings, Instabilities in multiserotype disease models with antibody-dependent enhancement, Journal of Theoretical Biology 246 (2007) 18–27.
- H.J. Wearing and P. Rohani, Ecological and immunological determinants of dengue epidemics, Proc. Natl. Acad. Sci. USA 103 (2006) 11802–11807.
- Y. Nagao and K. Koelle, Decreases in dengue transmission may act to increase the incidence of dengue hemorrhagic fever, Proc. Natl. Acad. Sci 105 (2008) 2238–2243.
- R. Seydel, Practical bifurcation and stability analysis-from equilibrium to chaos, Springer-Verlag, New York, 1994.
- M. Golubitsky and D.G. Schaeffer, Singularities and groups in bifurcation theory, Springer, New York, 1985.
- D. Albers and J. Sprott, Routes to chaos in high-dimensional dynamical systems: A qualitative numerical study, Physica D 223 (2006) 194–207.
- M.P. Boer,B.W. Kooi and S.A.L.M. Kooijman, Homoclinic and heteroclinic orbits in a tri-trophic food chain, Journal of Mathematical Biology 39 (1999) 19–38.
- M.P. Boer, B.W. Kooi and S.A.L.M. Kooijman, Multiple attractors and boundary crises in a tri-trophic food chain, Mathematical Biosciences 169 (2001) 109–128.
- B.W. Kooi and M.P. Boer, Chaotic behaviour of a predator-prey system, Dynamics of Continuous, Discrete and Impulsive Systems, Series B: Applications and Algorithms 10 (2003) 259–272.
- B.W. Kooi, L.D.J. Kuijper and S.A.L.M. Kooijman, Consequence of symbiosis for food web dynamics, Journal of Mathematical Biology 49 (2004) 227–271.
- J.P. Eckmann, S. Oliffson-Kamphorst, D. Ruelle and S. Ciliberto, Liapunov exponents from time series, Phys. Rev. A 34 (1986) 4971–9.