**4. Multi-strain models motivated by dengue fever epidemiology: a review**

Multi-strain dynamics are generally modelled with SIR-type models and have demonstrated to show critical fluctuations with power law distributions of disease cases, exemplified in meningitis and dengue epidemiology (Massad et al., 2008; Stollenwerk & Jansen, 2003; Stollenwerk et al., 2004). Dengue models including multi-strain interactions via ADE but without temporary cross immunity period e.g. (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005) have shown deterministic chaos when strong infectivity on secondary infection was assumed. The addition of the temporary cross immunity period in such models shows a new chaotic attractor in an unexpected parameter region of reduced infectivity on secondary infection (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008, 2009, 2011 a), i.e. deterministic chaos was found in a wider parameter regions. This indicates that deterministic chaos is much more important in multi-strain models than previously thought, and opens new ways to data analysis of existing dengue time series, as will be shown below. It offers a promising perspective on parameter values inference from dengue cases notifications.

The basic multi-strain model divides the population into ten classes: susceptible to both strains, 1 and 2 (*S*), primarily infected with strain one (*I*1) or strain two (*I*2), recovered from the first infection with strain one (*R*1) or strain two (*R*2), susceptible with a previous infection with strain one (*S*1) or strain two (*S*2), secondarily infected with strain one when the first 12 Will-be-set-by-IN-TECH

deadly complication that is characterized by high fever and hemorrhagic phenomenae. DHF develops rapidly, usually over a period of hours, and resolves within 12 days in patients who receive appropriate fluid resuscitation. Otherwise, it can quickly progress to shock (CDC,

Treatment of uncomplicated dengue cases is only supportive, and severe dengue cases requires careful attention to fluid management and proactive treatment of hemorrhagic symptoms (CDC, 2011; WHO, 2009). A vaccine against dengue is not yet available, since it would have to simulate a protective immune response to all four serotypes (Stephenson, 2005), although several candidates of tetravalent vaccines are at various stages of development

Mathematical models describing the transmission of dengue viruses appeared in the literature early as 1970 (Fischer & Halstead, 1970). More recently, mathematical models describing the transmission of dengue viruses have focused on the ADE effect and temporary cross immunity trying to explain the irregular behavior of dengue epidemics. Such models ultimately aim to be used as a predictive tool with the objective to guide the policies of prevention and control of the dengue virus transmission, including the implementation of vaccination programs when the candidate dengue fever vaccines will be accessible. In the literature the multi-strain interaction leading to deterministic chaos via ADE has been described previously, e.g. (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005) but neglecting temporary cross immunity. Consideration of temporary cross immunity is rather complicated and up to now not in detail analyzed. Models formulated in (Loureço & Recker, 2010; Nagao & Koelle, 2008; Recker et al., 2009; Wearing & Rohani, 2006), did not investigate closer the possible dynamical structures. In (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008, 2009, 2011 a) by including temporary cross immunity into dengue models with ADE, a rich dynamic structure including deterministic chaos was found in wider and more biologically

**4. Multi-strain models motivated by dengue fever epidemiology: a review**

Multi-strain dynamics are generally modelled with SIR-type models and have demonstrated to show critical fluctuations with power law distributions of disease cases, exemplified in meningitis and dengue epidemiology (Massad et al., 2008; Stollenwerk & Jansen, 2003; Stollenwerk et al., 2004). Dengue models including multi-strain interactions via ADE but without temporary cross immunity period e.g. (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005) have shown deterministic chaos when strong infectivity on secondary infection was assumed. The addition of the temporary cross immunity period in such models shows a new chaotic attractor in an unexpected parameter region of reduced infectivity on secondary infection (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008, 2009, 2011 a), i.e. deterministic chaos was found in a wider parameter regions. This indicates that deterministic chaos is much more important in multi-strain models than previously thought, and opens new ways to data analysis of existing dengue time series, as will be shown below. It offers a promising perspective on parameter values inference from dengue cases notifications.

The basic multi-strain model divides the population into ten classes: susceptible to both strains, 1 and 2 (*S*), primarily infected with strain one (*I*1) or strain two (*I*2), recovered from the first infection with strain one (*R*1) or strain two (*R*2), susceptible with a previous infection with strain one (*S*1) or strain two (*S*2), secondarily infected with strain one when the first

2011; WHO, 2009).

(WHO, 2011).

realistic parameter regions.

infection was caused by strain two (*I*21) or for second time infected with strain two when the first infection was caused by strain one (*I*12). Notice that infection by one serotype confers life-long immunity to that serotype. Then the individuals recover from the secondary infection (*R*).

To capture differences in primary infection by one strain and secondary infection by another strain we consider a basic two-strain SIR-type model for the host population, which is only slightly refined as opposed to previously suggested models for dengue fever (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005).

The stochastic version of the multi-strain dengue model is now in complete analogy to the previously described SIR model, and the mean field ODE system for the multi-strain dengue model can be read from the following reaction scheme (24), describing the transitions for first infection with strain 1 and secondary infection with strain 2, and for the reverse process, where the first infection is caused by strain 2 and the secondary infection is caused by strain 1, the same reaction scheme can be used to describe the transitions by just changing labels.

$$\begin{aligned} S &+ I\_1 \stackrel{\beta}{\longrightarrow} I\_1 + I\_1\\ S &+ I\_{21} \stackrel{\phi\beta}{\longrightarrow} I\_1 + I\_{21} \\ I\_1 &\stackrel{\gamma}{\longrightarrow} R\_1\\ R\_1 &\stackrel{a}{\longrightarrow} S\_1\\ S\_1 + I\_2 &\stackrel{\beta}{\longrightarrow} I\_{12} + I\_2\\ S\_1 + I\_{12} &\stackrel{\phi\beta}{\longrightarrow} I\_{12} + I\_{12} \\ I\_{12} &\stackrel{\gamma}{\longrightarrow} R \end{aligned} \tag{24}$$

The demographic transitions are *<sup>S</sup>*, *<sup>I</sup>*1, *<sup>I</sup>*2, *<sup>R</sup>*1, *<sup>R</sup>*2, *<sup>S</sup>*1, *<sup>S</sup>*2, *<sup>I</sup>*12, *<sup>I</sup>*21, *<sup>R</sup> <sup>μ</sup>* −→ *S* defining the system of two strains completely (for more information on the deterministic ODE system and its parametrization, see (Aguiar et al., 2011 a)).

The complete system of ordinary differential equations for the two strain epidemiological system is given by Eq. system (25) and the dynamics are described as follows. Susceptibles to both strains can get the first infection with strain one or strain two with force of infection *<sup>β</sup><sup>I</sup> N* when the infection is acquired via an individual in his first infection or *φβ<sup>I</sup> <sup>N</sup>* when the infection is acquired via an individual in his second infection. They recover from the first infection with a recovery rate *γ*, conferring full and life-long immunity against the strain that they were exposed to, and also a short period of temporary cross-immunity *α* against the other strain, becoming susceptible to a second infection with a different strain. The susceptible with a previous infection gets the secondary infection with force of infection *<sup>β</sup><sup>I</sup> <sup>N</sup>* or *φβ<sup>I</sup> <sup>N</sup>* depending on whom (individual on his primary or secondary infection) is transmitting the infection. Then, with recovery rate *γ*, the individuals recover and become immune against all strains. We assume no epidemiological asymmetry between strains, i.e. infections with strain one or strain two contribute in the same way to the force of infection. Here, the only relevant difference concerning disease transmissibility is that the force of infection varies accordingly to the number of previous infections the hosts have experienced. The parameter *φ* in our model, is the ratio of secondary infection contribution to the force of infection. For more

*<sup>S</sup>*<sup>∗</sup> <sup>=</sup> *<sup>μ</sup><sup>N</sup>* <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)(*I*<sup>∗</sup>

� *N*

� *N*

<sup>1</sup> <sup>=</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>∗</sup>

<sup>2</sup> <sup>=</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>∗</sup>

<sup>1</sup> and *I*<sup>∗</sup>

� (*γ*+*μ*) *<sup>β</sup>* − 3 �

(*I*∗ <sup>2</sup> + *φ*<sup>2</sup> *I*<sup>∗</sup> 12)

(*I*∗ <sup>1</sup> + *φ*<sup>1</sup> *I*<sup>∗</sup> 21)

*I* ∗ <sup>21</sup> <sup>=</sup> <sup>1</sup> *φ*1

*I* ∗ <sup>12</sup> <sup>=</sup> <sup>1</sup> *φ*2

*S*∗

*S*∗

*R*∗ <sup>1</sup> <sup>=</sup> *<sup>γ</sup> α* + *μ I* ∗ 1

*R*∗ <sup>2</sup> <sup>=</sup> *<sup>γ</sup> α* + *μ I* ∗ <sup>2</sup> ,

*φ* +

<sup>1</sup> <sup>−</sup> *αγ* (*α*+*μ*)(*γ*+*μ*)

*φ* +

*I* ∗

*I* ∗ <sup>2</sup> = 0 .

∗ <sup>1</sup> + *I* ∗ <sup>2</sup> + *R*<sup>∗</sup>

by the balance equation for the total population size *N*, explicitly

*R*<sup>∗</sup> = *N* − (*S*<sup>∗</sup> + *I*

and the solution of the extinction of one of the strains is as follows

<sup>1</sup> <sup>−</sup> *αγ* (*α*+*μ*)(*γ*+*μ*)

where still the stationary values of *I*∗

⎡ ⎣

− ������ *N*<sup>2</sup> 4

*αγ* (*α*+*μ*)(*γ*+*μ*)

4 (*γ*+*μ*) *μ* �

> ⎡ ⎣

*αγ* (*α*+*μ*)(*γ*+*μ*)

2 (*γ*+*μ*) *μ* �

*I* ∗ <sup>1</sup> = *I* ∗ <sup>2</sup> = − *μ*

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 243

*<sup>β</sup>*1*S*<sup>∗</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*) <sup>−</sup> <sup>1</sup>

*<sup>β</sup>*2*S*<sup>∗</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*) <sup>−</sup> <sup>1</sup>

*N β*2

*N β*1

<sup>2</sup> have to be determined.

⎦ *N* (27)

� (*γ*+*μ*) *<sup>β</sup>* − 1 �

<sup>1</sup> <sup>−</sup> *αγ* (*α*+*μ*)(*γ*+*μ*)

*φ* �

<sup>21</sup>) . (29)

⎤ ⎦ ,

(28)

*N*2*μ*

�

2 (*γ*+*μ*)<sup>2</sup> *μ*

12

21

The solution of coexistence of both strains for *I*<sup>1</sup> = *I*<sup>2</sup> = *I*<sup>∗</sup> is given by the following expression

⎤

*φ* �

<sup>1</sup> <sup>=</sup> *<sup>μ</sup>N*(*<sup>β</sup>* <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)) (*γ* + *μ*)*β*

Finally, the stationary value of *R*∗, when hosts have been recovered from both strains, is given

<sup>1</sup> + *R*<sup>∗</sup>

The time series for *φ <* 1 shows that the total number of infected *I* := *I*<sup>1</sup> + *I*<sup>2</sup> + *I*<sup>12</sup> + *I*<sup>21</sup> stays quite away from zero, avoiding the chance of extinction in stochastic systems with reasonable

<sup>2</sup> + *S*<sup>∗</sup>

<sup>1</sup> + *S*<sup>∗</sup> <sup>2</sup> + *I* ∗ <sup>12</sup> + *I* ∗

⎤ ⎦ 2 + ⎡ ⎣

*φ* �

� (*γ*+*μ*) *<sup>β</sup>* − 3 � <sup>1</sup> + *I*<sup>∗</sup> 2 )

> � *I* ∗ 1

> � *I* ∗ 2

> > (26)

information on the parametrization of ADE and secondary dengue infection by *φ*, see (Aguiar et al., 2008; Ferguson et al., 1999). The parameter values are given in Table 1, if not otherwise explicitly stated.

$$\begin{aligned} \dot{S} &= -\frac{\beta}{N}S(I\_1 + \phi I\_{21}) - \frac{\beta}{N}S(I\_2 + \phi I\_{12}) + \mu(N - S) \\\\ \dot{I}\_1 &= \frac{\beta}{N}S(I\_1 + \phi I\_{21}) - (\gamma + \mu)I\_1 \\\\ \dot{I}\_2 &= \frac{\beta}{N}S(I\_2 + \phi I\_{12}) - (\gamma + \mu)I\_2 \\\\ \dot{R}\_1 &= \gamma I\_1 - (\alpha + \mu)R\_1 \\\\ \dot{R}\_2 &= \gamma I\_2 - (\alpha + \mu)R\_2 \\\\ \dot{S}\_1 &= -\frac{\beta}{N}S\_1(I\_2 + \phi I\_{12}) + aR\_1 - \mu S\_1 \\\\ \dot{S}\_2 &= -\frac{\beta}{N}S\_2(I\_1 + \phi I\_{21}) + aR\_2 - \mu S\_2 \\\\ \dot{I}\_{12} &= \frac{\beta}{N}S\_1(I\_2 + \phi I\_{12}) - (\gamma + \mu)I\_{12} \\\\ \dot{I}\_{21} &= \frac{\beta}{N}S\_2(I\_1 + \phi I\_{21}) - (\gamma + \mu)I\_{21} \\\\ \dot{R} &= \gamma(I\_{12} + I\_{21}) - \mu R \end{aligned} \tag{25}$$


Table 1. Parameter set, rates given in units per year, ratio without unit

The stationary states can be calculated analytically by setting the time derivatives in Eq. system (25) to zero,

14 Will-be-set-by-IN-TECH

information on the parametrization of ADE and secondary dengue infection by *φ*, see (Aguiar et al., 2008; Ferguson et al., 1999). The parameter values are given in Table 1, if not otherwise

*<sup>N</sup> <sup>S</sup>*(*I*<sup>2</sup> <sup>+</sup> *<sup>φ</sup>I*12) + *<sup>μ</sup>*(*<sup>N</sup>* <sup>−</sup> *<sup>S</sup>*)

<sup>2</sup> = *γI*<sup>2</sup> − (*α* + *μ*)*R*<sup>2</sup> (25)

*<sup>N</sup> <sup>S</sup>*(*I*<sup>1</sup> <sup>+</sup> *<sup>φ</sup>I*21) <sup>−</sup> *<sup>β</sup>*

*<sup>N</sup> <sup>S</sup>*(*I*<sup>1</sup> <sup>+</sup> *<sup>φ</sup>I*21) <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>1</sup>

*<sup>N</sup> <sup>S</sup>*(*I*<sup>2</sup> <sup>+</sup> *<sup>φ</sup>I*12) <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>2</sup>

*<sup>N</sup> <sup>S</sup>*1(*I*<sup>2</sup> <sup>+</sup> *<sup>φ</sup>I*12) + *<sup>α</sup>R*<sup>1</sup> <sup>−</sup> *<sup>μ</sup>S*<sup>1</sup>

*<sup>N</sup> <sup>S</sup>*2(*I*<sup>1</sup> <sup>+</sup> *<sup>φ</sup>I*21) + *<sup>α</sup>R*<sup>2</sup> <sup>−</sup> *<sup>μ</sup>S*<sup>2</sup>

*<sup>N</sup> <sup>S</sup>*1(*I*<sup>2</sup> <sup>+</sup> *<sup>φ</sup>I*12) <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>12</sup>

*<sup>N</sup> <sup>S</sup>*2(*I*<sup>1</sup> <sup>+</sup> *<sup>φ</sup>I*21) <sup>−</sup> (*<sup>γ</sup>* <sup>+</sup> *<sup>μ</sup>*)*I*<sup>21</sup>

*<sup>R</sup>*˙ <sup>=</sup> *<sup>γ</sup>*(*I*<sup>12</sup> <sup>+</sup> *<sup>I</sup>*21) <sup>−</sup> *<sup>μ</sup><sup>R</sup>* ,

*μ* new born susceptible rate 1/65*y* (UNWPP, 2008)

*β* infection rate 2*γ* (Ferguson et al., 1999)

*γ* recovery rate 52*y*−<sup>1</sup> (Gubler et al., 1981; WHO, 2009)

*α* temporary cross-immunity rate 2*y*−<sup>1</sup> (Matheus et al., 2005; SES, 2010)

The stationary states can be calculated analytically by setting the time derivatives in Eq.

Par. Description Values Ref *N* population size 100 —

*φ* ratio of contrib. to force of inf. variable —

system (25) to zero,

Table 1. Parameter set, rates given in units per year, ratio without unit

<sup>1</sup> = *γI*<sup>1</sup> − (*α* + *μ*)*R*<sup>1</sup>

explicitly stated.

*<sup>S</sup>*˙ <sup>=</sup> <sup>−</sup> *<sup>β</sup>*

˙ *<sup>I</sup>*<sup>1</sup> <sup>=</sup> *<sup>β</sup>*

˙ *<sup>I</sup>*<sup>2</sup> <sup>=</sup> *<sup>β</sup>*

*R*˙

*R*˙

*S*˙ <sup>1</sup> <sup>=</sup> <sup>−</sup> *<sup>β</sup>*

*S*˙ <sup>2</sup> <sup>=</sup> <sup>−</sup> *<sup>β</sup>*

˙*I*<sup>12</sup> <sup>=</sup> *<sup>β</sup>*

˙*I*<sup>21</sup> <sup>=</sup> *<sup>β</sup>*

$$\begin{aligned} S^\* &= \frac{\mu N - (\gamma + \mu)(I\_1^\* + I\_2^\*)}{\mu} \\ I\_{21}^\* &= \frac{1}{\Phi\_1} \left( \frac{N}{\beta\_1 S^\*} (\gamma + \mu) - 1 \right) I\_1^\* \\ I\_{12}^\* &= \frac{1}{\Phi\_2} \left( \frac{N}{\beta\_2 S^\*} (\gamma + \mu) - 1 \right) I\_2^\* \\ S\_1^\* &= \frac{(\gamma + \mu) I\_{12}^\*}{(I\_2^\* + \Phi\_2 I\_{12}^\*)} \frac{N}{\beta\_2} \\ S\_2^\* &= \frac{(\gamma + \mu) I\_{21}^\*}{(I\_1^\* + \Phi\_1 I\_{21}^\*)} \frac{N}{\beta\_1} \\ R\_1^\* &= \frac{\gamma}{\alpha + \mu} I\_1^\* \\ R\_2^\* &= \frac{\gamma}{\alpha + \mu} I\_2^\* \end{aligned} \tag{26}$$

where still the stationary values of *I*∗ <sup>1</sup> and *I*<sup>∗</sup> <sup>2</sup> have to be determined.

The solution of coexistence of both strains for *I*<sup>1</sup> = *I*<sup>2</sup> = *I*<sup>∗</sup> is given by the following expression

$$I\_1^\* = I\_2^\* = -\left[\frac{\frac{a\gamma}{(a+\mu)(\gamma+\mu)}\phi + \left(\frac{(\gamma+\mu)}{\beta} - 3\right)}{4\frac{(\gamma+\mu)}{\mu}\left(1 - \frac{a\gamma}{(a+\mu)(\gamma+\mu)}\phi\right)}\right]N\tag{27}$$

$$1 - \sqrt{\frac{N^2}{4} \left[ \frac{\frac{a\gamma}{(a+\mu)\left(\gamma+\mu\right)} \phi + \left(\frac{(\gamma+\mu)}{\beta} - 3\right)}{2\frac{(\gamma+\mu)}{\mu}\left(1 - \frac{a\gamma}{(a+\mu)\left(\gamma+\mu\right)}\phi\right)} \right]^2 + \left[\frac{N^2\mu\left(\frac{(\gamma+\mu)}{\beta} - 1\right)}{2\frac{(\gamma+\mu)^2}{\mu}\left(1 - \frac{a\gamma}{(a+\mu)\left(\gamma+\mu\right)}\phi\right)} \right]} \right]$$

and the solution of the extinction of one of the strains is as follows

$$\begin{aligned} I\_1^\* &= \frac{\mu N(\beta - (\gamma + \mu))}{(\gamma + \mu)\beta} \\\\ I\_2^\* &= 0 \quad . \end{aligned} \tag{28}$$

Finally, the stationary value of *R*∗, when hosts have been recovered from both strains, is given by the balance equation for the total population size *N*, explicitly

$$R^\* = N - \left(\mathcal{S}^\* + I\_1^\* + I\_2^\* + R\_1^\* + R\_2^\* + S\_1^\* + S\_2^\* + I\_{12}^\* + I\_{21}^\*\right) \tag{29}$$

The time series for *φ <* 1 shows that the total number of infected *I* := *I*<sup>1</sup> + *I*<sup>2</sup> + *I*<sup>12</sup> + *I*<sup>21</sup> stays quite away from zero, avoiding the chance of extinction in stochastic systems with reasonable

a)

as people with first infection.

λ<sup>i</sup>


unstable equilibria or limit cycles.

a)


ln(I)

0 0.5 1 1.5 2 2.5 3

parameter *φ*. In a) *α* = 2 (six month) and in b) *α* = 52 (one week).

0 0.2 0.4 0.6 0.8 1 1.2

<sup>φ</sup> b)

Fig. 13. In a) spectrum of the four largest Lyapunov exponents with changing parameter *φ* and fixed *α* = 2. In b) we show the one-parameter bifurcation diagram with temporary cross-immunity rate *α* = 2 and varying the ratio of secondary infection contribution to the force of infection *φ*. Solid lines denote stable equilibria or limit cycles, and dashed lines

We quantify the attractor structure, fixed point, limit cycle or chaotic attractor etc., by calculating Lyapunov exponents (Ott, 1993; Ruelle, 1989), which were calculated using an iterated technique along a trajectory using the QR decomposition algorithm via Householder matrices (see (Aguiar et al., 2008; Holzfuss & Lauterborn, 1989; Holzfuss & Parlitz, 1991)). Lyapunov exponents are essentially a generalization of eigenvalues determining stability versus instability along trajectories. A negative largest Lyapunov exponent indicates a stable fixed point as attractor, a zero largest Lyapunov exponent indicates a stable limit cycle and a positive largest Lyapunov exponent indicates a chaotic attractor. Fig. 13a) shows the


ln(I)

0 0.2 0.4 0.6 0.8 1

φ

H P TR P T

<sup>φ</sup> b)

Fig. 12. Bifurcation diagram for the local extrema of the overall infected with changing

attractor found first by Aguiar et al. (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008).

period, i.e. by putting *α* → ∞. The recovered individuals can be immediately infected with another strain, whereas consideration of temporary cross-immunity brings a new chaotic

This finding encourages to look closer to the parameter region of *φ <* 1, when dengue patients in a secondary infection evolving to severe disease because of the ADE phenomenon contribute less to the force of infection, and not more, as previous models suggested. This assumption is likely to be more realistic for dengue fever since the possible severity of a secondary infection may hospitalize people, not contributing to the force of infections as much


ln(I)

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 245

0 0.5 1 1.5 2 2.5 3

φ

system size (see Fig. 10 a)). The parameter region previously considered to model ADE effects on dengue epidemiology, i.e. *φ >* 1, leads to rather low troughs for the total number of infected giving unrealistically low numbers of infected. In Fig. 10 b) the logarithm of total number of infected goes as low as −70 for *φ* = 2.7 in the chaotic region of *φ >* 1. Population fluctuations would in this case drive almost surely the system to extinction.

Fig. 10. Time series of the logarithm of the overall infected (ln(*I*)) comparison: a) simulation for *φ* = 2.7 and b) simulation for *φ* = 0.6 for the same time interval.

The state space plots in terms of the variables *S* and the logarithm of the total number of infected *I* show a rich dynamical behavior with bifurcations from fixed point to limit cycles, till completely irregular behavior (see Fig. 11).

Fig. 11. Attractors for various values of *φ <* 1: a) fixed point for *φ* = 0.1, and b) limit cycle for *φ* = 0.4, and c) chaotic attractor for *φ* = 0.6.

Looking for higher values of *φ*, the chaotic attractor becomes unstable, just leaving simple limit cycles as attractors for large parameter regions beyond *φ* = 1 (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008). Only for much higher values of *φ >>* 1, another chaotic attractor appears, the classical "ADE chaotic attractor" (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008; Ferguson et al., 1999).

The bifurcation diagram was obtained plotting the local extrema of *ln*(*I*) over the varying parameter *φ* (see Fig. 12). Fixed points appear as one dot per parameter value, limit cycles appear as two dots, double-limit cycles as four dots, more complicated limit cycles as more dots, and chaotic attractors as continuously distributed dots for a single *φ* value (Ruelle, 1989). We observe two chaotic windows, one for *φ <* 1, where this dynamical behavior has never been described before, and also another one for *φ >* 1 (see Fig. 12a)) where the minimal values go to very low numbers of infected, which already has been described (see Fig. 12b)) in previous publications (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005). In Fig. 12b) the dynamical behavior was obtained when neglecting the temporary cross-immunity 16 Will-be-set-by-IN-TECH

system size (see Fig. 10 a)). The parameter region previously considered to model ADE effects on dengue epidemiology, i.e. *φ >* 1, leads to rather low troughs for the total number of infected giving unrealistically low numbers of infected. In Fig. 10 b) the logarithm of total number of infected goes as low as −70 for *φ* = 2.7 in the chaotic region of *φ >* 1. Population

> -70 -60 -50 -40 -30 -20 -10 0 10

ln(I1+I2+I12+I21)

0 20 40 60 80 100 120 140 160 180 200

t


ln(I1+I2+I12+I21) (t)

36 38 40 42 44 46 48 50

S(t)

fluctuations would in this case drive almost surely the system to extinction.

<sup>t</sup> b)

Fig. 10. Time series of the logarithm of the overall infected (ln(*I*)) comparison: a) simulation

The state space plots in terms of the variables *S* and the logarithm of the total number of infected *I* show a rich dynamical behavior with bifurcations from fixed point to limit cycles,

36 38 40 42 44 46 48 50

Fig. 11. Attractors for various values of *φ <* 1: a) fixed point for *φ* = 0.1, and b) limit cycle for

Looking for higher values of *φ*, the chaotic attractor becomes unstable, just leaving simple limit cycles as attractors for large parameter regions beyond *φ* = 1 (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008). Only for much higher values of *φ >>* 1, another chaotic attractor appears, the classical "ADE chaotic attractor" (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008;

The bifurcation diagram was obtained plotting the local extrema of *ln*(*I*) over the varying parameter *φ* (see Fig. 12). Fixed points appear as one dot per parameter value, limit cycles appear as two dots, double-limit cycles as four dots, more complicated limit cycles as more dots, and chaotic attractors as continuously distributed dots for a single *φ* value (Ruelle, 1989). We observe two chaotic windows, one for *φ <* 1, where this dynamical behavior has never been described before, and also another one for *φ >* 1 (see Fig. 12a)) where the minimal values go to very low numbers of infected, which already has been described (see Fig. 12b)) in previous publications (Billings et al., 2007; Ferguson et al., 1999; Schwartz et al, 2005). In Fig. 12b) the dynamical behavior was obtained when neglecting the temporary cross-immunity

S(t) c)

0 20 40 60 80 100 120 140 160 180 200

for *φ* = 2.7 and b) simulation for *φ* = 0.6 for the same time interval.


ln(I1+I2+I12+I21) (t)

a)

a)


Ferguson et al., 1999).

ln(I1+I2+I12+I21) (t)


till completely irregular behavior (see Fig. 11).

36 38 40 42 44 46 48 50

*φ* = 0.4, and c) chaotic attractor for *φ* = 0.6.

S(t) b)

ln(I1+I2+I12+I21)

Fig. 12. Bifurcation diagram for the local extrema of the overall infected with changing parameter *φ*. In a) *α* = 2 (six month) and in b) *α* = 52 (one week).

period, i.e. by putting *α* → ∞. The recovered individuals can be immediately infected with another strain, whereas consideration of temporary cross-immunity brings a new chaotic attractor found first by Aguiar et al. (Aguiar & Stollenwerk, 2007; Aguiar et al., 2008).

This finding encourages to look closer to the parameter region of *φ <* 1, when dengue patients in a secondary infection evolving to severe disease because of the ADE phenomenon contribute less to the force of infection, and not more, as previous models suggested. This assumption is likely to be more realistic for dengue fever since the possible severity of a secondary infection may hospitalize people, not contributing to the force of infections as much as people with first infection.

Fig. 13. In a) spectrum of the four largest Lyapunov exponents with changing parameter *φ* and fixed *α* = 2. In b) we show the one-parameter bifurcation diagram with temporary cross-immunity rate *α* = 2 and varying the ratio of secondary infection contribution to the force of infection *φ*. Solid lines denote stable equilibria or limit cycles, and dashed lines unstable equilibria or limit cycles.

We quantify the attractor structure, fixed point, limit cycle or chaotic attractor etc., by calculating Lyapunov exponents (Ott, 1993; Ruelle, 1989), which were calculated using an iterated technique along a trajectory using the QR decomposition algorithm via Householder matrices (see (Aguiar et al., 2008; Holzfuss & Lauterborn, 1989; Holzfuss & Parlitz, 1991)). Lyapunov exponents are essentially a generalization of eigenvalues determining stability versus instability along trajectories. A negative largest Lyapunov exponent indicates a stable fixed point as attractor, a zero largest Lyapunov exponent indicates a stable limit cycle and a positive largest Lyapunov exponent indicates a chaotic attractor. Fig. 13a) shows the

Fig. 15. The state flow diagram for the seasonal multi-strain model. The transition rate *μ*

The seasonal multi-strain model is represented in Fig. 15 by using a state flow diagram where the boxes represent the disease related stages and the arrows indicate the transition rates. In the same manner as described for the non-seasonal, the population is divided into ten classes, with constant size *N* = *S* + *I*<sup>1</sup> + *I*<sup>2</sup> + *R*<sup>1</sup> + *R*<sup>2</sup> + *S*<sup>1</sup> + *S*<sup>2</sup> + *I*<sup>12</sup> + *I*<sup>21</sup> + *R*. The complete system of ordinary differential equations for the seasonal multi-strain epidemiological can be written as shown in system (25), with the difference that now the parameter *β* takes the seasonal forcing

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 247

where *β*<sup>0</sup> is the infection rate, and *η* is the degree of seasonality. In this model, a susceptible individual can become infected also by meeting an infected individual from an external population (hence (*β*/*N* · *S* · *I*) goes to (*β*/*N* · *S* · (*I* + *ρ* · *N*))) contributing to the force of

The parameters are fixed, temporary cross immunity rate *α* = 2*y*−1, recovery rate *γ* = 52*y*−1, infection rate *<sup>β</sup>*<sup>0</sup> <sup>=</sup> <sup>2</sup> · *<sup>γ</sup>*, seasonality *<sup>η</sup>* <sup>=</sup> 0.35, import factor *<sup>ρ</sup>* <sup>=</sup> <sup>10</sup>−10, birth and death rate *μ* = 1/65*y* and the ratio of secondary infection contribution to the force of infection *φ* = 0.9. In Fig. 16a) the time series simulation results for the total number of infected (*I*<sup>1</sup> + *I*<sup>2</sup> + *I*<sup>12</sup> + *I*21) in the non-seasonal system, previously studied in (Aguiar et al., 2008), is shown. Besides showing an irregular pattern of outbreaks that happens every 5 years, the non-seasonal system and its time series are not able to represent dengue fever epidemiology that is characterized as a yearly cycle of incidences. By adding low seasonality into the system, the epidemic outbreaks appear every year (see (Aguiar et al., 2011 a)). However, between two large outbreaks there is a very low number of cases in subsequent years, which is also not data alike. In Fig. 7c), the time series simulation in the high seasonal system with a low import of infected contributing to the force of infection is shown. The addition of import into the

*β*(*t*) = *β*<sup>0</sup> · (1 + *η* · cos(*ω* · *t*)) , (30)

coming out of the class *R* represents the death rates of all classes, *S*, *I*1, *I*2, *R*1, *R*2, *S*1, *S*2, *I*12, *I*21, *R*, getting into the class *S* as a birth rate.

into account as a cosine function given explicitly by

infection with an import parameter *ρ*.

largest four Lyapunov exponents as a function of *φ*. We observe that for small *φ* up to 0.1 all four Lyapunov exponents are negative, indicating the stable fixed point solution. Then follows a region up to *φ* = 0.5 where the largest Lyapunov exponent is zero, characteristic for stable limit cycles. Above *φ* = 0.5 a positive Lyapunov exponent, clearly separated from the second largest Lyapunov exponent being zero, indicates deterministically chaotic attractors. In the chaotic window between *φ* = 0.5 and *φ* = 1 also periodic windows appear, giving a zero largest Lyapunov exponent. These findings are in good agreement with the numerical bifurcation diagram, Fig. 13b).

A further analysis of the bifurcation structure, in the region of interest of *φ <* 1, was performed using the numerical software AUTO (AUTO, 2009). Various bifurcations were found: Hopf bifurcation *H*(*φ* = 0.11326), pitchfork bifurcations *P*(*φ* = 0.41145, 0.99214), torus bifurcation *TR*(*φ* = 0.55069) and tangent bifurcations *T*(*φ* = 0.4.9406, 0.53874, 0.93103, 0.97825, 1.05242). In addition to this main bifurcation pattern we found two isolas, consisting of isolated limit cycles existing between two tangent bifurcations (see Fig. 13b), for more information on the isolas see (Aguiar et al., 2008, 2009). These results agree very well with the simulation results shown in the bifurcation diagram for the maxima and minima of the overall infected in Fig. 12a).

t

Fig. 14. Time series of DHF incidence in Thailand.

Dengue fever epidemiology is characterized as a yearly cycle of incidences (see Fig. **??**), therefore in order to be able to reproduce the yearly cycle in dengue incidence seasonal forcing and a low import of infected have to be included in the models. The first recorded epidemic of DHF in Thailand (population of approximately 66 million people (Wikipedia, 2011)) was in 1958 (WHO, 2009). The co-circulation of all four dengue serotypes and their capacity to produce severe dengue disease was demonstrated as early as 1960 in Bangkok, Thailand (Halstead et al., 1969). DHF occurred first only in Bangkok, but was disseminated to the whole region during the 1970*s* (Chareonsook et al., 1999; Gubler, 2002; Halstead et al., 1969). Physicians in Thailand are trained to recognize and treat dengue fever and practically all cases of DHF and DSS are hospitalized. A system for reporting communicable diseases including DHF/DSS was considered fully installed in 1974 and the data bank of DHF and DSS is available at the Ministry of Public Health, Bangkok (Chareonsook et al., 1999).

We extend the previously studied non-seasonal model by adding seasonal forcing, mimicking the vectorial dynamics, and a low import of infected individuals, which is realistic in the dynamics of infectious diseases, in order to get a more realistic pattern of dengue fever epidemics, with irregular, yearly and smooth outbreaks.

18 Will-be-set-by-IN-TECH

largest four Lyapunov exponents as a function of *φ*. We observe that for small *φ* up to 0.1 all four Lyapunov exponents are negative, indicating the stable fixed point solution. Then follows a region up to *φ* = 0.5 where the largest Lyapunov exponent is zero, characteristic for stable limit cycles. Above *φ* = 0.5 a positive Lyapunov exponent, clearly separated from the second largest Lyapunov exponent being zero, indicates deterministically chaotic attractors. In the chaotic window between *φ* = 0.5 and *φ* = 1 also periodic windows appear, giving a zero largest Lyapunov exponent. These findings are in good agreement with the numerical

A further analysis of the bifurcation structure, in the region of interest of *φ <* 1, was performed using the numerical software AUTO (AUTO, 2009). Various bifurcations were found: Hopf bifurcation *H*(*φ* = 0.11326), pitchfork bifurcations *P*(*φ* = 0.41145, 0.99214), torus bifurcation *TR*(*φ* = 0.55069) and tangent bifurcations *T*(*φ* = 0.4.9406, 0.53874, 0.93103, 0.97825, 1.05242). In addition to this main bifurcation pattern we found two isolas, consisting of isolated limit cycles existing between two tangent bifurcations (see Fig. 13b), for more information on the isolas see (Aguiar et al., 2008, 2009). These results agree very well with the simulation results shown in the bifurcation diagram for the maxima and minima of the overall infected in Fig.

1983 1986 1989 1992 1995 1998 2001 2004

t

Dengue fever epidemiology is characterized as a yearly cycle of incidences (see Fig. **??**), therefore in order to be able to reproduce the yearly cycle in dengue incidence seasonal forcing and a low import of infected have to be included in the models. The first recorded epidemic of DHF in Thailand (population of approximately 66 million people (Wikipedia, 2011)) was in 1958 (WHO, 2009). The co-circulation of all four dengue serotypes and their capacity to produce severe dengue disease was demonstrated as early as 1960 in Bangkok, Thailand (Halstead et al., 1969). DHF occurred first only in Bangkok, but was disseminated to the whole region during the 1970*s* (Chareonsook et al., 1999; Gubler, 2002; Halstead et al., 1969). Physicians in Thailand are trained to recognize and treat dengue fever and practically all cases of DHF and DSS are hospitalized. A system for reporting communicable diseases including DHF/DSS was considered fully installed in 1974 and the data bank of DHF and

DSS is available at the Ministry of Public Health, Bangkok (Chareonsook et al., 1999).

We extend the previously studied non-seasonal model by adding seasonal forcing, mimicking the vectorial dynamics, and a low import of infected individuals, which is realistic in the dynamics of infectious diseases, in order to get a more realistic pattern of dengue fever

Thailand

bifurcation diagram, Fig. 13b).

epidemics, with irregular, yearly and smooth outbreaks.

DHF incidence

Fig. 14. Time series of DHF incidence in Thailand.

12a).

Fig. 15. The state flow diagram for the seasonal multi-strain model. The transition rate *μ* coming out of the class *R* represents the death rates of all classes, *S*, *I*1, *I*2, *R*1, *R*2, *S*1, *S*2, *I*12, *I*21, *R*, getting into the class *S* as a birth rate.

The seasonal multi-strain model is represented in Fig. 15 by using a state flow diagram where the boxes represent the disease related stages and the arrows indicate the transition rates. In the same manner as described for the non-seasonal, the population is divided into ten classes, with constant size *N* = *S* + *I*<sup>1</sup> + *I*<sup>2</sup> + *R*<sup>1</sup> + *R*<sup>2</sup> + *S*<sup>1</sup> + *S*<sup>2</sup> + *I*<sup>12</sup> + *I*<sup>21</sup> + *R*. The complete system of ordinary differential equations for the seasonal multi-strain epidemiological can be written as shown in system (25), with the difference that now the parameter *β* takes the seasonal forcing into account as a cosine function given explicitly by

$$\beta(t) = \beta\_0 \cdot (1 + \eta \cdot \cos(\omega \cdot t)) \quad , \tag{30}$$

where *β*<sup>0</sup> is the infection rate, and *η* is the degree of seasonality. In this model, a susceptible individual can become infected also by meeting an infected individual from an external population (hence (*β*/*N* · *S* · *I*) goes to (*β*/*N* · *S* · (*I* + *ρ* · *N*))) contributing to the force of infection with an import parameter *ρ*.

The parameters are fixed, temporary cross immunity rate *α* = 2*y*−1, recovery rate *γ* = 52*y*−1, infection rate *<sup>β</sup>*<sup>0</sup> <sup>=</sup> <sup>2</sup> · *<sup>γ</sup>*, seasonality *<sup>η</sup>* <sup>=</sup> 0.35, import factor *<sup>ρ</sup>* <sup>=</sup> <sup>10</sup>−10, birth and death rate *μ* = 1/65*y* and the ratio of secondary infection contribution to the force of infection *φ* = 0.9.

In Fig. 16a) the time series simulation results for the total number of infected (*I*<sup>1</sup> + *I*<sup>2</sup> + *I*<sup>12</sup> + *I*21) in the non-seasonal system, previously studied in (Aguiar et al., 2008), is shown. Besides showing an irregular pattern of outbreaks that happens every 5 years, the non-seasonal system and its time series are not able to represent dengue fever epidemiology that is characterized as a yearly cycle of incidences. By adding low seasonality into the system, the epidemic outbreaks appear every year (see (Aguiar et al., 2011 a)). However, between two large outbreaks there is a very low number of cases in subsequent years, which is also not data alike. In Fig. 7c), the time series simulation in the high seasonal system with a low import of infected contributing to the force of infection is shown. The addition of import into the

a)

c)



λ<sup>i</sup>

λ<sup>i</sup>

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1 1.2

<sup>φ</sup> b)

<sup>φ</sup> d)

spectrum and in d) the time series for the seasonal model with import.

with *�* = 0.001. (for details on the perturbed system see (Aguiar et al., 2011 a)).

a well defined maximum each year of irregular hight for the Northern Provinces.

The inspection of the available DHF incidence data in Thailand shows a smooth behavior with

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 249

In Fig. 18 the Lyapunov spectrum for both non-seasonal model and the seasonal model with import are shown and compared concerning the prediction horizon of the monthly peaks in the multi-strain dengue model time series. We take as an example the Dominant Lyapunov Exponent (DLE) for *φ* = 0.9 in the region where the system is chaotic (positive DLE). For the non-seasonal system, the DLE = 0.04 giving around 25 years of prediction horizon in the monthly time series (see Fig. 18b)), whereas for the seasonal system with import,the DLE= 0.118 giving around 8.5 years of prediction horizon in the monthly time series. It is clear that the addition of seasonal forcing into the system by itself decreases the practical predictability, however, the addition of a low import into the seasonally forced system helps to get a more complex dynamics and a better prediction horizon in the monthly time series. In order to get a qualitative insight into the predictability in the monthly sampled time series, i.e. to show how the original system behaves under a small perturbation we plot two different trajectories of the same system (for the non-seasonal model in Fig. 18b), and for the seasonal model with import of infected in Fig. 18d)), where the perturbed system (black line) is compared with the original model simulation (red line). To get the trajectory of the perturbed system, we keep the last point of the transient of the original system and use those values as starting values to compute the new and perturbed trajectory. The perturbation is given by *Sp* = *S* + *R* · *�* and *Rp* = *R* · (1.0 − *�*), where *Sp* is the susceptibles perturbed and *Rp* is the recovered perturbed

Secondary Infections

Secondary Infections

Fig. 18. Qualitative insight into the predictability in the monthly time series. In a) the Lyapunov spectrum and in b) the time series for the non-seasonal model. In c) the Lyapunov

1985 1990 1995 2000 2005 2010 2015 2020 2025 2030

1985 1990 1995 2000 2005 2010

t

 attractor perturbed

 attractor perturbed

t

Fig. 16. Time series simulations. In a) time series simulation for the non-seasonal model (*η* = 0). In b) time series simulation for the seasonal model with a low import of infected. Here, the degree of seasonality is *η* = 0.35 and the import of infected *ρ* = 10−10.

Fig. 17. Bifurcation diagram for the seasonal model with import. Here, the degree of seasonality *η* = 0.35 and the import factor *ρ* = 10−10.

seasonal system gives a much more realistic pattern of dengue fever epidemics, with irregular, yearly and smooth outbreaks. The system has a reasonable size (the number of infected stays quite away from zero), avoiding the chance of extinction in stochastic systems. For detailed analysis on the attractors in state space for the seasonal model, see (Aguiar et al., 2011 a). The bifurcation diagram for the seasonal model with import is shown in Fig. 17.

For the seasonal model with import AUTO predicted a torus bifurcation *TR* at *φ* = 0.13, and at *φ* = 0.522 which are also predicted very well when comparing with the results given by the Lyapunov exponent calculation. In the limiting case where the amplitude of the seasonal forcing is zero, the torus bifurcation *TR* of the seasonally forced system coincides with the Hopf bifurcation *H* of the non-seasonal system, as was shown in (Aguiar et al., 2011 a).

20 Will-be-set-by-IN-TECH

Fig. 16. Time series simulations. In a) time series simulation for the non-seasonal model (*η* = 0). In b) time series simulation for the seasonal model with a low import of infected.

0 0.2 0.4 0.6 0.8 1

seasonal system gives a much more realistic pattern of dengue fever epidemics, with irregular, yearly and smooth outbreaks. The system has a reasonable size (the number of infected stays quite away from zero), avoiding the chance of extinction in stochastic systems. For detailed analysis on the attractors in state space for the seasonal model, see (Aguiar et al., 2011 a). The

For the seasonal model with import AUTO predicted a torus bifurcation *TR* at *φ* = 0.13, and at *φ* = 0.522 which are also predicted very well when comparing with the results given by the Lyapunov exponent calculation. In the limiting case where the amplitude of the seasonal forcing is zero, the torus bifurcation *TR* of the seasonally forced system coincides with the Hopf bifurcation *H* of the non-seasonal system, as was shown in (Aguiar et al., 2011 a).

Fig. 17. Bifurcation diagram for the seasonal model with import. Here, the degree of

φ

Here, the degree of seasonality is *η* = 0.35 and the import of infected *ρ* = 10−10.


bifurcation diagram for the seasonal model with import is shown in Fig. 17.

ln(I)

a)

b)

seasonality *η* = 0.35 and the import factor *ρ* = 10−10.

Fig. 18. Qualitative insight into the predictability in the monthly time series. In a) the Lyapunov spectrum and in b) the time series for the non-seasonal model. In c) the Lyapunov spectrum and in d) the time series for the seasonal model with import.

In Fig. 18 the Lyapunov spectrum for both non-seasonal model and the seasonal model with import are shown and compared concerning the prediction horizon of the monthly peaks in the multi-strain dengue model time series. We take as an example the Dominant Lyapunov Exponent (DLE) for *φ* = 0.9 in the region where the system is chaotic (positive DLE). For the non-seasonal system, the DLE = 0.04 giving around 25 years of prediction horizon in the monthly time series (see Fig. 18b)), whereas for the seasonal system with import,the DLE= 0.118 giving around 8.5 years of prediction horizon in the monthly time series. It is clear that the addition of seasonal forcing into the system by itself decreases the practical predictability, however, the addition of a low import into the seasonally forced system helps to get a more complex dynamics and a better prediction horizon in the monthly time series. In order to get a qualitative insight into the predictability in the monthly sampled time series, i.e. to show how the original system behaves under a small perturbation we plot two different trajectories of the same system (for the non-seasonal model in Fig. 18b), and for the seasonal model with import of infected in Fig. 18d)), where the perturbed system (black line) is compared with the original model simulation (red line). To get the trajectory of the perturbed system, we keep the last point of the transient of the original system and use those values as starting values to compute the new and perturbed trajectory. The perturbation is given by *Sp* = *S* + *R* · *�* and *Rp* = *R* · (1.0 − *�*), where *Sp* is the susceptibles perturbed and *Rp* is the recovered perturbed with *�* = 0.001. (for details on the perturbed system see (Aguiar et al., 2011 a)).

The inspection of the available DHF incidence data in Thailand shows a smooth behavior with a well defined maximum each year of irregular hight for the Northern Provinces.

parameter estimation of which application is in general restricted to fairly simple dynamical

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 251

Being able to predict future outbreaks of dengue in the absence of human interventions is a major goal if one wants to understand the effects of control measures. Even after a dengue virus vaccine has become accessible, this holds true for the implementation of a vaccination program. For example, to perform a vaccine trial in a year with normally low numbers of cases would make statistical tests of vaccine efficacy much more difficult than when it was performed in a year with naturally high numbers of cases. Thus predictability of the next season's hight of the dengue peak on the basis of deterministic balance of infected and

The research presented in this book chapter has been supported by the Fundação para a Ciência e a Tecnologia, FCT grant SFRH/BD/43236/2008, and has been further supported by the Portuguese FCT project PTDC/MAT/115168/2009 and by the EU project EPIWORK under Framework Program 7. The work was carried out at Centro de Matemática e Aplicações

Aguiar, M. & Stollenwerk, N. (2007). A new chaotic attractor in a basic multi-strain

Aguiar, M., Kooi, B. W. & Stollenwerk, N. (2008). Epidemiology of Dengue Fever: A

Aguiar, M., Stollenwerk, N. & Kooi, B. W. (2009). Torus bifurcations, isolas and chaotic

*International Journal of Computer Mathematics*, 86, 1867–1877, ISSN 0020-7160. Aguiar, M., et al. (2011 a). The role of seasonality and import in a minimalistic multi-strain

Aguiar, M., Kooi B. W., & Stollenwerk N. (2011 b). Scaling of stochasticity in DHF epidemics.

Alonso, D., McKane, A. & Pascual, M. (2006). Stochastic Amplification in Epidemics. *Journal*

Anderson, R. M. & May, R. M. (1979). Population biology of infectious diseases: Part I. *Nature*,

Doedel J. E. and Oldeman, B. (2009). AUTO 07P âA ¸ ˘ S Continuation and bifurcation software

for ordinary differential equations. *Technical Report: Concordia University, Montreal,*

*of the Royal Society Interface*, 4, 575–582, ISSN 1742-5689.

*Canada*, Retrieved from http://indy.cs.concordia.ca/auto/

epidemiological model with temporary cross-immunity. *arXiv:0704.3174v1 [nlin.CD]*.

Model with Temporary Cross-Immunity and Possible Secondary Infection Shows Bifurcations and Chaotic Behaviour in Wide Parameter Regions. *Math. Model. Nat.*

attractors in a simple dengue model with ADE and temporary cross immunity.

dengue model capturing differences between primary and secondary infections: complex dynamics and its implications for data analysis. *Journal of Theoretical Biology*,

scenarios.

**6. Conclusions**

**7. Acknowledgments**

**8. References**

susceptible would be of major practical use.

Fundamentais, Science Faculty, Lisbon University, Portugal.

*Phenom.*, 4, 48–70, ISSN 0973-5348.

289, 181–196, ISSN 0022-5193.

280, 361–67, ISSN 0028-0836.

*Under review*.

Fig. 19. Empirical DHF incidence data matched with the model simulation.

We take the Province of Chiang Mai as a case study where the empirical DHF incidence data and the time series simulation for the seasonal model with import are compared (see Fig. 19)). The seasonal model with import shows complex dynamics and qualitatively a very good result when comparing empirical DHF and simulations. However, the extended model needs to be parametrized on data referring to incidence of severe disease.

#### **5. Discussions**

In this chapter we presented the properties of the basic SIR epidemic model for infectious diseases with a summary of the analysis of the dynamics, identifying the thresholds and equilibrium points in order to introduce notation, terminology. The results that were generalized to more advanced models motivated by dengue fever epidemiology.

The epidemiology of dengue fever was described presenting the relevant biological features that are taken into the modeling process. Then, multi-strain models previously described in the literature were presented. We focused in a minimal model motivated by dengue fever epidemiology, formulated first by Aguiar et al. (Aguiar & Stollenwerk, 2007), where the notion of at least two different strains is needed to describe differences between primary infections, often asymptomatic, and secondary infection, associated with the severe form of the disease. We discussed the role of seasonal forcing and the import of infected individuals in such systems, the biological relevance and its implications for the analysis of the available dengue data. The extended model (Aguiar et al., 2011 a) shows complex dynamics and qualitatively a good agreement between empirical DHF monitoring data and the obtained model simulation. This suggests that the used parameter set can be the starting set for a more detailed parameter estimation procedure. Such a technical parameter estimation is notoriously difficult for chaotic time series but temporally local approaches are possible (He et al., 2010; Ionides et al., 2006). At the moment only such minimalistic models have a chance to be qualitatively understood well and eventually tested against existing data.

The introduction of stochasticity is needed to explain the fluctuations observed in some of the available data sets, revealing a scenario where noise and complex deterministic skeleton strongly interact (Aguiar et al., 2011 b). For large enough population size, the stochastic system can be well described by the deterministic skeleton gaining 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 fairly simple dynamical scenarios.

### **6. Conclusions**

22 Will-be-set-by-IN-TECH

1985 1990 1995 2000 2005 2010 2015 2020 2025 2030

 Simulated Chiang Mai

t

We take the Province of Chiang Mai as a case study where the empirical DHF incidence data and the time series simulation for the seasonal model with import are compared (see Fig. 19)). The seasonal model with import shows complex dynamics and qualitatively a very good result when comparing empirical DHF and simulations. However, the extended model needs

In this chapter we presented the properties of the basic SIR epidemic model for infectious diseases with a summary of the analysis of the dynamics, identifying the thresholds and equilibrium points in order to introduce notation, terminology. The results that were

The epidemiology of dengue fever was described presenting the relevant biological features that are taken into the modeling process. Then, multi-strain models previously described in the literature were presented. We focused in a minimal model motivated by dengue fever epidemiology, formulated first by Aguiar et al. (Aguiar & Stollenwerk, 2007), where the notion of at least two different strains is needed to describe differences between primary infections, often asymptomatic, and secondary infection, associated with the severe form of the disease. We discussed the role of seasonal forcing and the import of infected individuals in such systems, the biological relevance and its implications for the analysis of the available dengue data. The extended model (Aguiar et al., 2011 a) shows complex dynamics and qualitatively a good agreement between empirical DHF monitoring data and the obtained model simulation. This suggests that the used parameter set can be the starting set for a more detailed parameter estimation procedure. Such a technical parameter estimation is notoriously difficult for chaotic time series but temporally local approaches are possible (He et al., 2010; Ionides et al., 2006). At the moment only such minimalistic models have a chance

generalized to more advanced models motivated by dengue fever epidemiology.

to be qualitatively understood well and eventually tested against existing data.

The introduction of stochasticity is needed to explain the fluctuations observed in some of the available data sets, revealing a scenario where noise and complex deterministic skeleton strongly interact (Aguiar et al., 2011 b). For large enough population size, the stochastic system can be well described by the deterministic skeleton gaining insight into the relevant parameter values purely on topological information of the dynamics, rather than classical

Fig. 19. Empirical DHF incidence data matched with the model simulation.

to be parametrized on data referring to incidence of severe disease.

Secondary Infections

**5. Discussions**

Being able to predict future outbreaks of dengue in the absence of human interventions is a major goal if one wants to understand the effects of control measures. Even after a dengue virus vaccine has become accessible, this holds true for the implementation of a vaccination program. For example, to perform a vaccine trial in a year with normally low numbers of cases would make statistical tests of vaccine efficacy much more difficult than when it was performed in a year with naturally high numbers of cases. Thus predictability of the next season's hight of the dengue peak on the basis of deterministic balance of infected and susceptible would be of major practical use.

#### **7. Acknowledgments**

The research presented in this book chapter has been supported by the Fundação para a Ciência e a Tecnologia, FCT grant SFRH/BD/43236/2008, and has been further supported by the Portuguese FCT project PTDC/MAT/115168/2009 and by the EU project EPIWORK under Framework Program 7. The work was carried out at Centro de Matemática e Aplicações Fundamentais, Science Faculty, Lisbon University, Portugal.

#### **8. References**


Holzfuss, J. & Lauterborn, W. (1989). Liapunov exponents from a time series of acoustic chaos.

Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study 253

Holzfuss, J. & Parlitz, U. (1991). Lyapunov exponents from time series. *Lecture Notes in*

Ionides, E., Breto, C., & King, A. A. (2006). Inference for nonlinear dynamical systems. *Proc.*

Lourenço, J. & Recker, M. (2010). Viral and Epidemiological Determinants of the Invasion

Mackenzie, J. S., Gubler, D. J. & Petersen, L. R. (2004). Emerging flaviviruses: the spread and

Massad, E.,et al. (2008). Scale-free network of a dengue epidemic. *Applied Mathematics and*

Matheus, S. et al. (2005). Discrimination between Primary and Secondary Dengue Virus

Monath T. P., (1994). Dengue: The risk to developed and developing countries. *Proc. Natl.*

Nagao, Y. & Koelle, K.(2008). Decreases in dengue transmission may act to increase the

Nisalak, A. et al. (2003). Serotype-specific dengue virus circulation and dengue disease in

Ott, E. (1993). *Chaos in Dynamical Systems*, ISBN 0-52143-799-7. (Cambridge University Press,

Pediatric Dengue Vaccine Initiative. International Vaccine Institute (IVI). *Global Burden of Dengue*. Retrieved from http://www.pdvi.org/about\_dengue/GBD.asp Recker, M. et al. (2009). Immunological serotype interactions and their effect on the

Rigau-Pérez, J. G., et al. (1998). Dengue and dengue haemorrhagic fever. *The Lancet*, 352,

Rosen, L. et al. (1983). Transovarial transmission of dengue viruses by mosquitoes: A. albopictus and A. aegypti. *Am. J. Trop. Med. Hyg.*, 32, 1108–19, ISSN: 0002-9637. D. Ruelle. (1989). *Chaotic Evolution and Strange Attractors*, ISBN 978-0-52136-272-6 (Cambridge

Pers comm.: Francisco Lemos, Secretaria de Estado de Saúde de Minas Gerais, Brazil; Sônia

Schwartz, I. B., et al. (2005). Chaotic desynchronization of multi-strain diseases. *Physical*

Stephenson, J. R., (2005). Understanding dengue pathogenesis: implications for vaccine

design. *Bull. World Health Organ.*, 83, 308–14, ISSN 0042-9686.

Sample. *Journal of Clinical Microbiology*, 43, 2793–2797, ISSN: 0095-1137. Mattheij, R. M. M. and Molenaar, J. (1996). *Ordinary differential equations in theory and practice.*,

Dynamics of Novel Dengue Genotypes. *PLoS Negl. Trop. Dis.*, 4, e894, ISSN

resurgence of Japanese encephalitis, West Nile and dengue viruses. *Nature Medicine*

Infection by an Immunoglobulin G Aviditnoy Test Using a Single Acute-Phase Serum

incidence of dengue hemorrhagic fever. *Proc. Natl. Acad. Sci. USA*, 105, 2238–2243,

Bangkok, Thailand from 1973 to 1999. *Am. J. Trop. Med. Hyg.*, 68, 191–202, ISSN:

epidemiological pattern of dengue. *Proc. R. Soc. B.*, 276, 2541–2548, ISSN: 1471-2954.

Diniz, Fundação Ezequiel Dias, Minas Gerais, Brazil and Scott Halstead, Pedriatic

*Physical Review A*, 39, 2146–2152, ISSN 1943-2879.

*Natl. Acad. Sci. USA*, 103, 18438–18443, ISSN 0027-8424.

ISBN 0-89871-531-8, (Wiley, Chichester and New York).

*Acad. Sci. U.S.A.*, 91, 2395–2400, ISSN 0027-8424.

*Mathematics*, 1486, 263–270, ISSN: 0075-8434.

*Review*, 12, S98–S109, ISSN : 1078-8956.

*Computation*, 195, 376–381, ISSN: 0096-3003.

1935-2735.

ISSN 0027-8424.

Cambridge, 2nd edition).

971–77, ISSN: 0140-6736.

University Press, Cambridge).

Dengue Vaccine Initiative, Maryland, USA.

*Review*, E 72, 066201–6, ISSN 1943-2879.

0002-9637.


24 Will-be-set-by-IN-TECH

Billings, L., et al. (2007). Instabilities in multiserotype disease models with

Centers for Disease Control and Prevention. (2011). *Dengue*. Retrieved from

Chareonsook, O.et al. (1999). Changing epidemiology of dengue hemorrhagic fever in

Dejnirattisai, W. et al. (2010). Cross-Reacting Antibodies Enhance Dengue Virus Infection in

Favier, C., et al. (2005). Influence of spatial heterogeneity on an emerging infectious disease: the case of dengue epidemics. *Proc. Biol. Sci.*, 272, 1171–7, ISSN 0962-8452. Ferguson, N., Anderson, R. and Gupta, S. (1999). The effect of antibody-dependent

Fischer, D. B. & Halstead, S. B. (1970). Observations related to pathogenesis of dengue

Gillespie, D.T. (1976). A general method for numerically simulating the stochastic time

Gillespie, D.T. (1978). Monte Carlo simulation of random walks with residence time

Gubler, J. D., Suharyono, W., Tan, R., Abidin, M. and Sie, A. (1981). Viraemia in patients with

Gubler D. J., (2002). Epidemic dengue/dengue hemorrhagic fever as a public health, social

Guzmán, M.G. et al. (2000). Epidemiologic Studies on Dengue in Santiago de Cuba, 1997. *Am.*

Guzmán, M.G. et al. (2010). Dengue: a continuing global threat. *Nature Reviews Microbiology*,

Halstead S. B., et al. (1969). Dengue and chikungunya virus infection in man in Thailand,

Halstead, S.B. (1982). Immune enhancement of viral infection. *Progress in Allergy*, 31,

Halstead, S. B. (1994). Antibody-dependent Enhancement of Infection: A Mechanism for

He, D., Ionides, E. L., King, A. A. (2010). Plug-and-play inference for disease dynamics:

493–516, ISBN 0-87969-429-7. (Cold Spring Harbor Laboratory Press). Halstead, S.B. (2003). Neutralization and antibody-dependent enhancement of dengue

viruses. *Advances in Virus Research*, 60, 421–467, ISSN 0065-3527.

pathogens. *Proc. Natl. Acad. Sci. USA*, 96, 790–94, ISSN 0027-8424.

enhancement on the transmission dynamics and persistence of multiple-strain

hemorrhagic fever. V. Examination of age specific sequential infection rates using

evolution of coupled chemical reactions. *Journal of Computational Physics*, 22, 403–434,

dependent transition probability rates. *Journal of Computational Physics*, 28, 395–407,

naturally acquired dengue infection. *Bull. World Health Organ.*, 59, 623–630, ISSN

and economic problem in the 21st century. *Trends in Microbiology*, 10, 100–103, ISSN:

1962–1964. V. Epidemiologic observations outside Bangkok. *Am. J. Trop. Med. Hyg.*

Indirect Virus Entry into Cells. *Cellular Receptors for Animal Viruses*, 28, Chapter 25,

measles in large and small populations as a case study. *J. R. Soc. Interface*, 7, 271–283,

Thailand. *Epidemiol. Infect.*, 122, 161–166, ISSN 0950-2688.

Humans. *Science*, 328, 745–748, ISSN 0036-8075.

a mathematical model. *J. Biol. Med.*, 42, 329–49.

*J. Epidemiol.*, 152, 793–799, ISSN 0002-9262.

8, S7–S16, ISSN : 1740-1526.

18, 1022–33, ISSN: 0002-9637.

301–364,ISSN 0079-6034.

ISSN 1742-5689.

ISSN 0022-5193.

ISSN 0021-9991.

ISSN 0021-9991.

0042-9686.

0966-842X.

http://www.cdc.gov/dengue/

antibody-dependent enhancement. *Journal of Theoretical Biology*, 246, 18–27,


**13** 

 *Russia* 

B. Lapin and M. Chikobava

*Research Institute of Medical Primatology of the Russian Academy of Medical Sciences, Sochi* 

**Epidemiology of Simian Polyomavirus SV40 in** 

DNA-containing polyomavirus SV40 was isolated in 1961 from poliovaccine prepared on kidney cell culture of M.mulatta naturally infected with this virus. Main part of poliovaccines (IPV and OPV) was contaminated with SV40. From 1955 till 1963 hundred millions of people in USA, Russia (USSR), and in some other countries were immunized with contaminated SV40 Salk and Sabin vaccine and a lot of people became carriers of SV40. This virus was detected also in some human tumors. Epidemiological researches were not carried out in Russia before 2009. By present time possibility and ways of SV40 excretion into environment and possibility of its horizontal spreading is established. From the beginning of 60s (in Russia a little bit later) vaccines were free from SV40 contamination. In this article the data about SV40 infection of population in Moscow, St. Petersburg, in some cities of the Black Sea coast of Caucasus, in Novosibirsk and in one area of Krasnoyarsk

region are presented., Analysis of the source of SV40 people infection is performed.

it was established later (Eddy et al., 1961; Sweet & Hilleman, 1960).

Small size DNA virus was isolated from mice in 1953. The virus was called "polyomavirus" SV40 (Gross, 1951; Gross, 1953). because it was able to produce multiple malignant tumors

Simian SV40 was the second discovered polyomavirus. It was isolated in 1961 from poliovaccine propagated on M.mulatta kidney cell cultures naturally infected with SV40 as

SV 40 is a small DNA virus which does not have envelope; its capsid consisted of 72 capsomers. The size of virion is 45 nm; the viral genome is represented by two-stranded DNA aprox. 5000 bp long. On character of the changes caused in cell cultures, the virus also has received the name "vacuolating virus"(Sweet & Hilleman, 1960; Butel et al., 1999).

During next years polyomaviruses were isolated from different animals including different monkey species and from humans with different pathology as well (Simon,2008; Gjoerup,

In 1957 baboon polyomavirus was discovered (Gardner et al., 1989); polyovavirus of African green monkeys was discovered in 1979 (zur Hausen & Giosmann, 1979), after that polyomaviruses of M.fascicularis and chimpanzee were discovered (but not isolated). Two

**1. Introduction** 

in newborn rodents.

2010; Voevodin, 2009).).

**Different Areas of Russian Federation (RF)** 


http://mathworld.wolfram.com/Kermack-McKendrickModel.html


http://www.who.int/vaccine\_research/diseases/vector/en/index1.html#virology
