**2. The SIR epidemic model**

The SIR epidemic model divides the population into three classes: susceptible (*S*), Infected (*I*) and Recovered (*R*). It can be applied to infectious diseases where waning immunity can happen, and assuming that the transmission of the disease is contagious from person to person, the susceptibles become infected and infectious, are cured and become recovered. After a waning immunity period, the recovered individual can become susceptible again. This model was for the first time proposed by William Ogilvy Kermack and Anderson Gray McKendrick in 1927 (Weisstein, 2010). The model was brought back to prominence after decades of neglect by Anderson and May (Anderson & May, 1979).

In the simple SIR epidemics without strain structure of the pathogens we have the following reaction scheme for the possible transitions from one to another disease related state, susceptibles *S*, infected *I* and recovered *R*,

2 Will-be-set-by-IN-TECH

of disease, the duration of infectivity, the infection rate, the waning immunity, and so forth. In such a way, differential equation models are a simplified representation of reality, which are designed to facilitate prediction and calculation of rates of change as functions of the

There are two common approaches in modeling, the deterministic and the stochastic one. In the first case, the model is one in which the variable states are uniquely determined by parameters in the model and by sets of previous states of these variables. In mathematics, a deterministic system is a system in which no randomness is involved in the development of future states of the system. In a stochastic model, randomness is present, and variable states are not described by unique values, but rather by probability distributions. Stochastic epidemic models are appropriate stochastic processes that can be used to model disease propagation. Disease propagation is an inherently stochastic phenomenon and there are a number of reasons why one should use stochastic models to capture the transmission process. Real life epidemics, in the absence of intervention from outside, can either go extinct with a limited number of individuals getting ultimately infected, or end up with a significant proportion of the population having contracted the disease in question. It is only stochastic, as opposed to deterministic, models that can capture this behavior and the probability of each

Only few stochastic processes can be solved explicitly. The simplest and most thoroughly studied stochastic model of epidemics are based on the assumption of homogeneous mixing, i.e. individuals interact randomly at a certain rate. The mean field approximation is a good approximation to be used in order to understand better the behavior of the stochastic systems in certain parameter regions, where the dynamics of the mean quantities are approximated by neglecting correlations, giving closed ordinary differential equations (ODE) systems, hence

In the following section of this chapter we present 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. The goal is to introduce notation, terminology, and results that will be generalized in later sections on more advanced models motivated by dengue fever

The SIR epidemic model divides the population into three classes: susceptible (*S*), Infected (*I*) and Recovered (*R*). It can be applied to infectious diseases where waning immunity can happen, and assuming that the transmission of the disease is contagious from person to person, the susceptibles become infected and infectious, are cured and become recovered. After a waning immunity period, the recovered individual can become susceptible again. This model was for the first time proposed by William Ogilvy Kermack and Anderson Gray McKendrick in 1927 (Weisstein, 2010). The model was brought back to prominence after

In the simple SIR epidemics without strain structure of the pathogens we have the following reaction scheme for the possible transitions from one to another disease related state,

mathematically deterministic systems which are easier to analyze.

decades of neglect by Anderson and May (Anderson & May, 1979).

epidemiology as an example of multi-strain systems.

**2. The SIR epidemic model**

conditions or the components of the system.

event taking place.

$$\begin{array}{c} \mathcal{S} + I \stackrel{\beta}{\longrightarrow} I + I \\ I \stackrel{\gamma}{\longrightarrow} \mathcal{R} \\ \mathcal{R} \stackrel{\alpha}{\longrightarrow} \mathcal{S} \end{array}$$

for a host population of *N* individuals, with contact and infection rate *β*, recovery rate *γ* and waning immunity rate *α*. The dynamic model in terms of ordinary differential equations (ODE) reads,

$$\dot{S} = -\frac{\beta}{N}IS + \mathfrak{a}(N - S - I) \tag{1}$$

$$\mathbf{M} = \frac{\beta}{N} \mathbf{I} \mathbf{S} - \gamma \mathbf{I} \quad , \tag{2}$$

where we use the time derivative *S*˙ = *dS*/*dt* with time *t* for a constant population size of *N* = *S* + *I* + *R* individuals. The solution of *R*(*t*) is given by *R*(*t*) = *N* − *I*(*t*) − *S*(*t*) which can be calculated using the solution of the ODEs. The susceptible individuals become infected with infection rate *β*, recover from the infection with recovery rate *γ* and become susceptible again after waning immunity rate *α*.

In Fig. 1 we show the dynamical behavior of the susceptible, infected and recovered individuals in a given population *N*, when solving the above ODE system.

Fig. 1. Time dependent solution simulation for the SIR epidemic model. With a population *N* = 100, and starting values *I* = 40, *S* = 60 and *R* = 0, we fixed *β* = 2.5, *α* = 0.1, and *γ* = 1. In green the dynamics for the susceptibles *S*(*t*), in pink the dynamics for the infected *I*(*t*) and in blue the dynamics of the recovered *R*(*t*). Note that *N* = 100 allows for the interpretation for the class abundances in percentages.

The basic SIR model has only fixed points as possible stationary solutions, that can be calculated setting the rates of change *S*˙ and ˙ *I* to zero. For the disease free equilibrium state, the solution is given by

$$I\_1^\* = 0 \tag{3}$$

$$S\_1^\* = N \tag{4}$$

**2.1 The disease free equilibrium state** For the disease free equilibrium state (*I*∗

by

and *γ* = 1.0.

a)

d)


 0 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 0.001

Infected (β < γ)

Infected (β > γ)

the fixed point *I*∗

99.99 99.992 99.994 99.996 99.998 100

eigendirection uλ1 eigendirection uλ2

Susceptible

eigendirection uλ1 eigendirection uλ2

99.99 99.992 99.994 99.996 99.998 100

Susceptible

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

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

0 0.5 1 1.5 2 2.5

β

99.99 99.992 99.994 99.996 99.998 100

eigendirection uλ1 eigendirection uλ2

Susceptible

eigendirection uλ1 eigendirection uλ2

99.99 99.992 99.994 99.996 99.998 100

Susceptible

<sup>1</sup> is stable. The eigenvalues *λ*<sup>1</sup> and *λ*<sup>2</sup> are equal at the point *β* = *γ* − *α* and

c)

Infected (β > γ)

f)


 0 0.0005 0.001 0.0015

Infected (β < γ)

99.99 99.992 99.994 99.996 99.998 100

eigendirection uλ1 eigendirection uλ2

Susceptible

eigendirection uλ1 eigendirection uλ2

99.9918 99.9936 99.9954 99.9972 99.999

Susceptible

<sup>1</sup> = 0 is unstable. The

Fig. 3. Eigenvalues for the disease free equilibrium state as functions of *β* when fixing *α* = 0.1


b)

e)

for *β > γ*, *λ*<sup>1</sup> is positive and *λ*<sup>2</sup> is negative, therefore the fixed point *I*<sup>∗</sup>

*β* = 0.999, in d) *β* = 1.001, in e) *β* = 1.1, and in f) *β* = 1.3.



Fig. 4. eigenvectors for the disease free equilibrium state in function of *β*. For a population *N* = 100, where *I* = 0.0001, *S* = 99.99, and *R* = *N* − *I* − *S*, in a) *β* = 0.9, in b) *β* = 0.97, in c)

When looking at Eq. (10) and Eq. (11) we see that for *β < γ* both eigenvalues are negative, i.e.

stability of the system changes when one of the eigenvalues of the system becomes zero. At

Infected (β > γ)

Infected (β < γ)

λ1 λ2

Eigenvalues

<sup>1</sup>), Eq.(3) and Eq. (4), the eigenvalues are given

*λ*<sup>1</sup> = *β* − *γ* (10) *λ*<sup>2</sup> = −*α* . (11)

and for the disease endemic equilibrium state, the solution is

$$I\_2^\* = N \left( 1 - \frac{\gamma}{\beta} \right) \frac{\alpha}{(\alpha + \gamma)} \tag{5}$$

$$S\_2^\* = N \frac{\gamma}{\beta} \quad . \tag{6}$$

The epidemic dynamic as a function of the parameter *β* shows the spread of the epidemic when *β > γ* (see Fig. 2 a)), and its extinction when *β < γ* (see Fig. 2 b)).

Fig. 2. Epidemic dynamics as a function of *β*. With the same initial values as used in Fig. 1, we plot time dependent solutions *I*(*t*) for several *β* values. In a)*β* ∈ [1.5, 2.5], with a resolution Δ*β* = 0.1 and in b) *β* ∈ [0, 0.9] where Δ*β* = 0.2.

In order to analyze the stability of the equilibrium states, we look at the Jacobian matrix and its eigenvalues. Let the dynamics for the state *x* := (*S*, *I*) be *f*(*x*), hence *<sup>d</sup> dt x* = *f*(*x*) which explicitly gives Δ*x* := *x*(*t*) − *x*<sup>∗</sup> as a small perturbation around the fixed point *x*∗. We linearize the dynamic *<sup>d</sup> dt*Δ*<sup>x</sup>* <sup>=</sup> *<sup>d</sup> dt*(*x*(*t*) − *x*∗) applying Taylor's expansion

$$
\underline{f}(\underline{\mathbf{x}}^\* + \Delta \underline{\mathbf{x}}) = \underline{f}(\underline{\mathbf{x}}^\*) + \left. \frac{d\underline{f}}{d\underline{\mathbf{x}}} \right|\_{\underline{\mathbf{x}}^\*} \cdot (\Delta \underline{\mathbf{x}}) + \mathcal{O}((\Delta \underline{\mathbf{x}})^2) \tag{7}
$$

with *f*(*x*∗) = 0 for the fixed point and neglecting higher order terms. For our system we have the following linear differential equation system

$$\frac{d}{dt}\begin{pmatrix} S(t) - S^\* \\ I(t) - I^\* \end{pmatrix} = \begin{pmatrix} \frac{\partial f}{\partial S} \ \frac{\partial f}{\partial I} \\\\ \frac{\partial g}{\partial S} \ \frac{\partial g}{\partial I} \end{pmatrix} \Bigg|\_{\begin{pmatrix} S \\ I \end{pmatrix} = \begin{pmatrix} S^\* \\ I^\* \end{pmatrix}} \cdot \begin{pmatrix} S - S^\* \\ I - I^\* \end{pmatrix} \tag{8}$$

where *f* := (*f* , *g*) and the Jacobian matrix is explicitly given by

$$\begin{pmatrix} -\frac{\mathcal{\beta}}{N}I^\*-a & -\frac{\mathcal{\beta}}{N}S^\*-a\\\\ \frac{\mathcal{\beta}}{N}I^\* & \frac{\mathcal{\beta}}{N}S^\*-\gamma \end{pmatrix} =: A \tag{9}$$

where we have to insert for *S*∗ and *I*∗, the respective steady states. In order to decoupled the linear differential equation system, we diagonalize the matrix *A*, (9), with the eigenvalue decomposition *A u* = *λ u*, *u* is an eigenvector of *A*, and *λ* is an eigenvalue of *A* corresponding to the eigenvector *u*.

The eigenvalues can be calculated setting the determinant of [*A* − *λ* **1**] equal zero.

#### **2.1 The disease free equilibrium state**

4 Will-be-set-by-IN-TECH

The epidemic dynamic as a function of the parameter *β* shows the spread of the epidemic

� *α*

I(t,β)

(*<sup>α</sup>* <sup>+</sup> *<sup>γ</sup>*) (5)

*<sup>β</sup>* . (6)

0 5 10 15 20 25

t

· (Δ*x*) + O((Δ*x*)

� ·

� *<sup>S</sup>* <sup>−</sup> *<sup>S</sup>*<sup>∗</sup> *I* − *I*<sup>∗</sup>

�

⎟⎠ <sup>=</sup>: *<sup>A</sup>* (9)

*dt x* = *f*(*x*) which

(8)

<sup>2</sup>) (7)

� <sup>1</sup> <sup>−</sup> *<sup>γ</sup> β*

and for the disease endemic equilibrium state, the solution is

a)

the dynamic *<sup>d</sup>*

*dt*Δ*<sup>x</sup>* <sup>=</sup> *<sup>d</sup>*

*d dt*

corresponding to the eigenvector *u*.

the following linear differential equation system

� *<sup>S</sup>*(*t*) <sup>−</sup> *<sup>S</sup>*<sup>∗</sup> *I*(*t*) − *I*<sup>∗</sup>

I(t,β)

*I* ∗ <sup>2</sup> = *N*

*S*∗ <sup>2</sup> <sup>=</sup> *<sup>N</sup> <sup>γ</sup>*

when *β > γ* (see Fig. 2 a)), and its extinction when *β < γ* (see Fig. 2 b)).

0 10 20 30 40 50

resolution Δ*β* = 0.1 and in b) *β* ∈ [0, 0.9] where Δ*β* = 0.2.

<sup>t</sup> b)

we plot time dependent solutions *I*(*t*) for several *β* values. In a)*β* ∈ [1.5, 2.5], with a

*dt*(*x*(*t*) − *x*∗) applying Taylor's expansion

its eigenvalues. Let the dynamics for the state *x* := (*S*, *I*) be *f*(*x*), hence *<sup>d</sup>*

⎛

*∂ f ∂S ∂ f ∂I*

*∂g ∂S ∂g ∂I*

*<sup>N</sup> <sup>I</sup>*<sup>∗</sup> <sup>−</sup> *<sup>α</sup>* <sup>−</sup> *<sup>β</sup>*

*<sup>N</sup> <sup>I</sup>*<sup>∗</sup> *<sup>β</sup>*

The eigenvalues can be calculated setting the determinant of [*A* − *λ* **1**] equal zero.

⎜⎝

*<sup>f</sup>*(*x*<sup>∗</sup> <sup>+</sup> <sup>Δ</sup>*x*) = *<sup>f</sup>*(*x*∗) + *d f*

� =

where *f* := (*f* , *g*) and the Jacobian matrix is explicitly given by ⎛

⎜⎝ − *β*

*β*

Fig. 2. Epidemic dynamics as a function of *β*. With the same initial values as used in Fig. 1,

In order to analyze the stability of the equilibrium states, we look at the Jacobian matrix and

explicitly gives Δ*x* := *x*(*t*) − *x*<sup>∗</sup> as a small perturbation around the fixed point *x*∗. We linearize

*dx*

with *f*(*x*∗) = 0 for the fixed point and neglecting higher order terms. For our system we have

⎞

� � � � � � � � *S I* � = � *S*<sup>∗</sup> *I*∗

*<sup>N</sup> S*<sup>∗</sup> − *α*

⎞

*<sup>N</sup> S*<sup>∗</sup> − *γ*

where we have to insert for *S*∗ and *I*∗, the respective steady states. In order to decoupled the linear differential equation system, we diagonalize the matrix *A*, (9), with the eigenvalue decomposition *A u* = *λ u*, *u* is an eigenvector of *A*, and *λ* is an eigenvalue of *A*

⎟⎠

� � � � � *x*∗ For the disease free equilibrium state (*I*∗ <sup>1</sup> and *S*<sup>∗</sup> <sup>1</sup>), Eq.(3) and Eq. (4), the eigenvalues are given by

$$
\lambda\_1 = \beta - \gamma \tag{10}
$$

$$
\lambda\_2 = -\mathfrak{a} \quad . \tag{11}
$$

Fig. 3. Eigenvalues for the disease free equilibrium state as functions of *β* when fixing *α* = 0.1 and *γ* = 1.0.

Fig. 4. eigenvectors for the disease free equilibrium state in function of *β*. For a population *N* = 100, where *I* = 0.0001, *S* = 99.99, and *R* = *N* − *I* − *S*, in a) *β* = 0.9, in b) *β* = 0.97, in c) *β* = 0.999, in d) *β* = 1.001, in e) *β* = 1.1, and in f) *β* = 1.3.

When looking at Eq. (10) and Eq. (11) we see that for *β < γ* both eigenvalues are negative, i.e. the fixed point *I*∗ <sup>1</sup> is stable. The eigenvalues *λ*<sup>1</sup> and *λ*<sup>2</sup> are equal at the point *β* = *γ* − *α* and for *β > γ*, *λ*<sup>1</sup> is positive and *λ*<sup>2</sup> is negative, therefore the fixed point *I*<sup>∗</sup> <sup>1</sup> = 0 is unstable. The stability of the system changes when one of the eigenvalues of the system becomes zero. At

part *i*

solution

eigenvalues. For *β < γ* the fixed *I*∗

For *β > γ*, the fixed point *I*∗

a)


correspondent eigenvector *u*<sup>1</sup> is given by

Eigenvalues

become complex, where the real part *a* gives the contraction or expansion, and the imaginary

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

changes stability when *β<sup>c</sup>* = *γ* and becomes complex when *λ*<sup>1</sup> = *λ*<sup>2</sup> (see Fig. 5).

0.9 0.95 1 1.05 1.1

λ1 λ2

β

*<sup>u</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> � 1 + � *<sup>a</sup>*<sup>−</sup> √ *b γ*+*α*

�|*b*<sup>|</sup> gives the frequency of oscillations of the trajectories spiraling into the fixed point as is shown in Fig. 5. The parabola curve shows the contraction and expansion of the

b)

Fig. 5. Plot of the eigenvalues as functions of *β* when fixing *α* = 0.1 and *γ* = 1.0. In a) the real part *a* of the eigenvalues gives the contraction and expansion of the trajectories. For *β < γ* the fixed point is unstable (*λ*<sup>1</sup> *>* 0). The system changes stability when *β<sup>c</sup>* = *γ* and becomes complex when *λ*<sup>1</sup> = *λ*2. Here, the straight green line represents only the real parts of the complex eigenvalues obtained putting √−*<sup>b</sup>* <sup>=</sup> 0. In b) the imaginary part of the eigenvalues gives the frequency of the oscillatory behavior on the trajectory toward at the fixed point.

The correspondent eigenvectors of the disease endemic equilibrium state can be calculated for the eigenvalues in the same manner as it was shown above. For the first eigenvalue *λ*1, the

⎛

⎜⎜⎝

*a*− √ *b γ*+*α* · �

� ·

In Fig. 6 we show the eigenvectors for the disease endemic equilibrium state in function of *β*.

when *λ*<sup>1</sup> �= *λ*2. By including the respective eigenvalues, Eq. (10) and Eq.(11), we get as a

⎛

1

*u*2

−*αt*

*a*+ √*b γ*+*α*

⎜⎝

*u*<sup>1</sup> + *C*2*eλ*<sup>2</sup> *<sup>t</sup>*

*u*<sup>1</sup> + *C*2*e*

�2 ·

and for the second eigenvalue *λ*2, the correspondent eigenvector *u*<sup>2</sup> is given by

For the real eigenvalue the general solution of the linearized system is given by

*x*(*t*) = *C*1*eλ*<sup>1</sup> *<sup>t</sup>*

*x*(*t*) = *C*1*e*(*β*−*γ*)*<sup>t</sup>*

*<sup>u</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> � 1 + � *<sup>a</sup>*+ √*b γ*+*α*

Frequency

 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014

1

1 1+ � *<sup>a</sup>*<sup>−</sup> √ *b γ*+*α* �2

⎞

<sup>2</sup> point is unstable, with one positive eigenvalue.

<sup>2</sup> becomes stable with both eigenvalues negative. The system

0.9 0.95 1 1.05 1.1

β

⎞

⎟⎟⎠ . (18)

⎟⎠ . (19)

*u*<sup>2</sup> (20)

this point, *β<sup>c</sup>* = *γ*, when *I*∗ <sup>1</sup> becomes unstable and *I*<sup>∗</sup> <sup>2</sup> stable. Fig. 3 shows the eigenvalues for the disease free equilibrium state as functions of *β*.

To calculate the corresponding eigenvectors we use (*A* − *λ***1**) *u* = 0, with *λ<sup>i</sup>* and *ui* =: *u*1*<sup>i</sup> u*2*<sup>i</sup>* . For the first eigenvalue, *λ*1, the correspondent eigenvector *u*<sup>1</sup> is giving by

$$\underline{\mu}\_1 = \frac{1}{\sqrt{1 + \left(\frac{\gamma - \beta - \alpha}{\beta + \alpha}\right)^2}} \cdot \left(\left(\frac{1}{\frac{\gamma - \beta - \alpha}{\beta + \alpha}}\right)^{\gamma}\right) \quad \gamma$$

and for *λ*2, the correspondent eigenvector *u*<sup>2</sup> is is giving by

$$
\underline{u}\_2 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}
$$

.

In Fig. 4 we show the eigenvectors for the disease free equilibrium state as functions of *β*, when fixing *α* = 0.1, and *γ* = 1. We plot the eigenvectors, *u*<sup>1</sup> (blue line) and *u*<sup>2</sup> (green line) on top of the trajectory of the infected individuals (red line). Note that *λ*<sup>1</sup> = *λ*<sup>2</sup> at *β* = (*γ* − *α*), i.e. the eigenvectors *u*<sup>1</sup> = *u*<sup>2</sup> (see Fig. 4a)). By increasing *β* toward the critical value *β<sup>c</sup>* = *γ* the trajectory needs longer time to hit the fixed point (see Fig. 4b) and 4c)). For *β > γ*, the trajectory goes toward the other fixed point *I*∗ <sup>2</sup> (see Fig. 4d) and 4e)).

#### **2.2 The disease endemic equilibrium state**

For the disease endemic equilibrium state (*I*∗ <sup>2</sup> and *S*<sup>∗</sup> <sup>2</sup>), Eq. (5) and Eq. (6), the eigenvalues are giving by

$$\lambda\_1 = -\frac{a}{2}\left(1 + \frac{\beta - \gamma}{a + \gamma}\right) + \sqrt{\left[\frac{a}{2}\left(1 + \frac{\beta - \gamma}{a + \gamma}\right)\right]^2 - (\beta - \gamma)a} \tag{12}$$

$$\lambda\_2 = -\frac{a}{2}\left(1 + \frac{\beta - \gamma}{a + \gamma}\right) - \sqrt{\left[\frac{a}{2}\left(1 + \frac{\beta - \gamma}{a + \gamma}\right)\right]^2 - (\beta - \gamma)a} \quad . \tag{13}$$

In order to simplify the notation, let <sup>−</sup> *<sup>α</sup>* 2 1 + *<sup>β</sup>*−*<sup>γ</sup> α*+*γ* <sup>=</sup>: *<sup>a</sup>* and *α* 2 1 + *<sup>β</sup>*−*<sup>γ</sup> α*+*γ* <sup>2</sup> − (*β* − *γ*)*α* =: *b*. If *b >* 0 the eigenvalues are real numbers, giving the contraction or expansion of the trajectories near to the considered fixed point, and can be written as

$$
\lambda\_1 = a + \sqrt{b} \tag{14}
$$

$$
\lambda\_2 = a - \sqrt{b} \quad . \tag{15}
$$

If *b <* 0, the eigenvalues

$$
\lambda\_1 = a + i\sqrt{|b|} \tag{16}
$$

$$
\lambda\_2 = a - i\sqrt{|b|} \tag{17}
$$

6 Will-be-set-by-IN-TECH

<sup>2</sup> stable. Fig. 3 shows the eigenvalues for

<sup>2</sup>), Eq. (5) and Eq. (6), the eigenvalues are

*b* (14)

*b* . (15)



− (*β* − *γ*)*α* (12)

− (*β* − *γ*)*α* . (13)

− (*β* − *γ*)*α* =: *b*.

 *u*1*<sup>i</sup> u*2*<sup>i</sup>* .

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

For the first eigenvalue, *λ*1, the correspondent eigenvector *u*<sup>1</sup> is giving by

*<sup>u</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> 1 +

and for *λ*2, the correspondent eigenvector *u*<sup>2</sup> is is giving by

<sup>1</sup> <sup>+</sup> *<sup>β</sup>* <sup>−</sup> *<sup>γ</sup> α* + *γ*

<sup>1</sup> <sup>+</sup> *<sup>β</sup>* <sup>−</sup> *<sup>γ</sup> α* + *γ*

trajectories near to the considered fixed point, and can be written as

 + *α* 2 

 − *α* 2 

2 1 + *<sup>β</sup>*−*<sup>γ</sup> α*+*γ* 

trajectory goes toward the other fixed point *I*∗

**2.2 The disease endemic equilibrium state** For the disease endemic equilibrium state (*I*∗

> *<sup>λ</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup> *<sup>α</sup>* 2

> *<sup>λ</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup> *<sup>α</sup>* 2

In order to simplify the notation, let <sup>−</sup> *<sup>α</sup>*

If *b <* 0, the eigenvalues

giving by

To calculate the corresponding eigenvectors we use (*A* − *λ***1**) *u* = 0, with *λ<sup>i</sup>* and *ui* =:

 *<sup>γ</sup>*−*β*−*<sup>α</sup> β*+*α*

*u*<sup>2</sup> =

2 ·

 1 0 .

In Fig. 4 we show the eigenvectors for the disease free equilibrium state as functions of *β*, when fixing *α* = 0.1, and *γ* = 1. We plot the eigenvectors, *u*<sup>1</sup> (blue line) and *u*<sup>2</sup> (green line) on top of the trajectory of the infected individuals (red line). Note that *λ*<sup>1</sup> = *λ*<sup>2</sup> at *β* = (*γ* − *α*), i.e. the eigenvectors *u*<sup>1</sup> = *u*<sup>2</sup> (see Fig. 4a)). By increasing *β* toward the critical value *β<sup>c</sup>* = *γ* the trajectory needs longer time to hit the fixed point (see Fig. 4b) and 4c)). For *β > γ*, the

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

If *b >* 0 the eigenvalues are real numbers, giving the contraction or expansion of the

*<sup>λ</sup>*<sup>1</sup> <sup>=</sup> *<sup>a</sup>* <sup>+</sup> <sup>√</sup>

*<sup>λ</sup>*<sup>2</sup> <sup>=</sup> *<sup>a</sup>* <sup>−</sup> <sup>√</sup>

*λ*<sup>1</sup> = *a* + *i*

*λ*<sup>2</sup> = *a* − *i*

 1 *<sup>γ</sup>*−*β*−*<sup>α</sup> β*+*α* ,

<sup>2</sup> (see Fig. 4d) and 4e)).

<sup>1</sup> <sup>+</sup> *<sup>β</sup>* <sup>−</sup> *<sup>γ</sup> α* + *γ*

<sup>1</sup> <sup>+</sup> *<sup>β</sup>* <sup>−</sup> *<sup>γ</sup> α* + *γ*

=: *a* and

<sup>2</sup>

<sup>2</sup>

 *α* 2 1 + *<sup>β</sup>*−*<sup>γ</sup> α*+*γ* <sup>2</sup>

this point, *β<sup>c</sup>* = *γ*, when *I*∗

the disease free equilibrium state as functions of *β*.

become complex, where the real part *a* gives the contraction or expansion, and the imaginary part *i* �|*b*<sup>|</sup> gives the frequency of oscillations of the trajectories spiraling into the fixed point as is shown in Fig. 5. The parabola curve shows the contraction and expansion of the eigenvalues. For *β < γ* the fixed *I*∗ <sup>2</sup> point is unstable, with one positive eigenvalue.

For *β > γ*, the fixed point *I*∗ <sup>2</sup> becomes stable with both eigenvalues negative. The system changes stability when *β<sup>c</sup>* = *γ* and becomes complex when *λ*<sup>1</sup> = *λ*<sup>2</sup> (see Fig. 5).

Fig. 5. Plot of the eigenvalues as functions of *β* when fixing *α* = 0.1 and *γ* = 1.0. In a) the real part *a* of the eigenvalues gives the contraction and expansion of the trajectories. For *β < γ* the fixed point is unstable (*λ*<sup>1</sup> *>* 0). The system changes stability when *β<sup>c</sup>* = *γ* and becomes complex when *λ*<sup>1</sup> = *λ*2. Here, the straight green line represents only the real parts of the complex eigenvalues obtained putting √−*<sup>b</sup>* <sup>=</sup> 0. In b) the imaginary part of the eigenvalues gives the frequency of the oscillatory behavior on the trajectory toward at the fixed point.

The correspondent eigenvectors of the disease endemic equilibrium state can be calculated for the eigenvalues in the same manner as it was shown above. For the first eigenvalue *λ*1, the correspondent eigenvector *u*<sup>1</sup> is given by

$$\underline{\mu}\_{1} = \frac{1}{\sqrt{1 + \left(\frac{a-\sqrt{b}}{\gamma+a}\right)^2}} \cdot \begin{pmatrix} 1\\ \frac{a-\sqrt{b}}{\gamma+a} \cdot \frac{1}{\sqrt{1 + \left(\frac{a-\sqrt{b}}{\gamma+a}\right)^2}} \end{pmatrix} \tag{18}$$

and for the second eigenvalue *λ*2, the correspondent eigenvector *u*<sup>2</sup> is given by

$$\underline{u}\_2 = \frac{1}{\sqrt{1 + \left(\frac{a + \sqrt{b}}{\gamma + a}\right)}} \cdot \begin{pmatrix} 1\\ \frac{a + \sqrt{b}}{\gamma + a} \end{pmatrix} \quad . \tag{19}$$

In Fig. 6 we show the eigenvectors for the disease endemic equilibrium state in function of *β*. For the real eigenvalue the general solution of the linearized system is given by

$$\underline{\mathbf{x}}(t) = \mathbf{C}\_1 e^{\lambda\_1 t} \underline{u}\_1 + \mathbf{C}\_2 e^{\lambda\_2 t} \underline{u}\_2$$

when *λ*<sup>1</sup> �= *λ*2. By including the respective eigenvalues, Eq. (10) and Eq.(11), we get as a solution

$$\underline{\mathbf{x}}(t) = \mathbf{C}\_1 e^{(\beta - \gamma)t} \underline{u\_1} + \mathbf{C}\_2 e^{-\mathbf{a}t} \underline{u\_2} \tag{20}$$

a)

I(t)

0 10 20 30 40 50

t

*N* = 1000 and in c) *N* = 100000.

mean values obtaining

b)

For mean values of infected �*I*� and susceptibles �*S*�, defined as e.g.

�*I*�(*t*) :=

*d*

*d dt*�*I*� <sup>=</sup> *<sup>β</sup>*

show here continued oscillations (Alonso et al., 2006).

I(t)

0 10 20 30 40 50

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

t

Fig. 7. Stochastic simulations for the basic SIR epidemic model. Here 10 realizations are plotted. We fixed *α* = 0.1, *γ* = 1 and *β* = 2.5. The deterministic trajectory is shown (pink line) top of the stochastic realizations for different population size *N* a. In a) *N* = 100, in b)

> *N* ∑ *S*=0

*N* ∑ *I*=0

one can calculate the dynamics by inserting the master equation into the definition of the

with �*R*� = *N* − �*S*�−�*I*�. For more details on the calculations see e.g. (Stollenwerk & Jansen, 2010). These equations for the mean dynamics include now due to the nonlinear transition rates in the master equation also higher moments �*S* · *I*�. The simplest approximation to obtain a closed ODE system is to neglect cross-correlations �*S* · *I*�−�*S*�·�*I*� ≈ 0, the so-called mean field approximation (originally introduced for spatially extended systems in statistical physics (Stollenwerk et al., 2010)). Hence, the equation system (23) gives with identifying the higher moment �*S* · *I*� = �*S*�·�*I*� by a product of simple moments gives again the ODE system for SIR system, as it was just presented above. For certain parameter regions the mean field approximation describes the system well in terms of its mean dynamics and only small fluctuations around it. Then the previously shown analysis of the system is appropriate. However, noise can stabilize transients, a feature which becomes important in parameter regions where in the deterministic description a fixed point is reached via decreasing oscillations, as we have observed them in the SIR system. The noisy system would

In Fig. 7 we compare the deterministic and stochastic dynamics and we see that the magnitude of stochastic fluctuations decreases when the population size increases. However, the good approximation (see Fig. 7c)) is only achieved when the population size is large enough (see

Almost all mathematical models of diseases start from the same basic premise: that the population can be subdivided into a set of distinct classes. The most commonly used

Fig. 7a) where most simulations die out very quickly for small population size).

*<sup>N</sup>* �*S I*�

*<sup>N</sup>* �*S I*� − *<sup>γ</sup>*�*I*�

*dt*�*S*� <sup>=</sup> *<sup>α</sup>*�*R*� − *<sup>β</sup>*

c)

I(t)

*I p*(*S*, *I*, *t*) . (22)

0 10 20 30 40 50

t

(23)

Fig. 6. eigenvectors for the disease endemic equilibrium state in function of *β*. For a population *N* = 100, where *I* = 0.001 and *S* = 99, we fixed *α* = 0.1, *γ* = 1 and vary *β*. We plot the eigenvectors, *u*<sup>1</sup> (blue line) and *u*<sup>2</sup> (greenline), on top of the trajectory of the infected individuals toward the second fixed point *I*∗ <sup>2</sup> . In a) *β* = 1.01, and in b) *β* = 1.025. In c) we show the oscillatory trajectory toward the fixed point when *β* = 2.5

where the eigenvector *u*<sup>1</sup> is the driving force for *t* → ∞, since *λ*<sup>1</sup> *> λ*<sup>2</sup> (see Fig. 4b) to 4f)). When *λ*<sup>1</sup> *< λ*<sup>2</sup> the eigenvector *u*<sup>2</sup> is the driving force for *t* → ∞.

For the special case where *λ*<sup>1</sup> = *λ*<sup>2</sup> =: *λ* and therefore the eigenvectors *u*<sup>1</sup> = *u*<sup>2</sup> =: *u* the general solution is then given by

$$\underline{\mathfrak{x}}(t) = \mathsf{C}\_1 e^{\lambda t} \underline{\mathfrak{u}} + \mathsf{C}\_2 (t e^{\lambda t} \underline{\mathfrak{u}} + e^{\lambda t} \underline{\mathfrak{w}}),$$

where *w* is the so called generalized eigenvector, satisfying (*A* − *λ***1**) *w* = *u*. In this case, for *t* → ∞, the eigenvector *u* is again the driving force (see Fig. 4a)).

For the complex eigenvalues, where the real part *a* gives the contraction or expansion, and the imaginary part *i* <sup>|</sup>*b*<sup>|</sup> gives the frequency of oscillations of the trajectories spiraling into the fixed point, the general solution of the linearized system is given by

$$\underline{\mathbf{x}}(t) = 2e^{\mathrm{d}t} \left( \left[ \mathbf{C}\_1 \cos \left( \sqrt{b}t \right) - \mathbf{C}\_2 \sin \left( \sqrt{b}t \right) \right] \underline{u}\_1 - \left[ \mathbf{C}\_1 \sin \left( \sqrt{b}t \right) + \mathbf{C}\_2 \cos \left( \sqrt{b}t \right) \right] \underline{u}\_2 \right) \right.$$

where *C*<sup>1</sup> and *C*<sup>2</sup> depend on the initial conditions and *ui*, respectively the real and imaginary parts of the complex eigenvector. This expression shows that the stability of the fixed point depends on the sign of *a*. For detailed information on the solution of a linear two dimensional ODE system, see (Mattheij & Molenaar, 1996).

The stochastic SIR epidemic is modeled as a time-continuous Markov process to capture population noise. The dynamics of the probability of integer infected and integer susceptibles, while the recovered follow from this due to constant population size, can be give as a master equation (van Kampen, 1992) in the following form

$$\frac{dp(\mathcal{S}, I, t)}{dt} = \frac{\mathcal{B}}{N}(\mathcal{S} + 1)(I - 1) \quad p(\mathcal{S} + 1, I - 1, t) + \gamma(I + 1) \quad p(\mathcal{S}, I + 1, t) \tag{21}$$

$$\mathcal{L}\_{\mathcal{L}}(\mathcal{N} \quad (\mathcal{L} \quad 1) \quad \mathcal{N} \quad \pi(\mathcal{L} \quad \mathcal{I} \quad 1, t)) \quad \left(\mathcal{S} \quad \bot \quad \bot \quad \neg \mathcal{N}(\mathcal{I} \quad \mathcal{I} \quad \bot)\right) \tag{22}$$

$$\begin{aligned} \ &+ \ \mathfrak{a}(N - (\mathcal{S} - 1) - I) \quad p(\mathcal{S} - 1, I, t) - \left(\frac{\mathcal{B}}{N} + \gamma I + \mathfrak{a}(N - \mathcal{S} - I)\right) \quad p(\mathcal{S}, I, t) \quad \dots \end{aligned}$$

This process can be simulated by the Gillespie algorithm giving stochastic realizations of infected and susceptibles in time (Gillespie, 1976, 1978).

8 Will-be-set-by-IN-TECH

97 97.5 98 98.5 99 99.5 100

eigendirection uλ1 eigendirection uλ2

Susceptible

where the eigenvector *u*<sup>1</sup> is the driving force for *t* → ∞, since *λ*<sup>1</sup> *> λ*<sup>2</sup> (see Fig. 4b) to 4f)).

For the special case where *λ*<sup>1</sup> = *λ*<sup>2</sup> =: *λ* and therefore the eigenvectors *u*<sup>1</sup> = *u*<sup>2</sup> =: *u* the

*u* + *C*2(*teλ<sup>t</sup>*

where *w* is the so called generalized eigenvector, satisfying (*A* − *λ***1**) *w* = *u*. In this case, for

For the complex eigenvalues, where the real part *a* gives the contraction or expansion, and the

where *C*<sup>1</sup> and *C*<sup>2</sup> depend on the initial conditions and *ui*, respectively the real and imaginary parts of the complex eigenvector. This expression shows that the stability of the fixed point depends on the sign of *a*. For detailed information on the solution of a linear two dimensional

The stochastic SIR epidemic is modeled as a time-continuous Markov process to capture population noise. The dynamics of the probability of integer infected and integer susceptibles, while the recovered follow from this due to constant population size, can be give as a master

This process can be simulated by the Gillespie algorithm giving stochastic realizations of

*u* + *eλ<sup>t</sup>*

<sup>|</sup>*b*<sup>|</sup> gives the frequency of oscillations of the trajectories spiraling into the

*<sup>N</sup>* (*<sup>S</sup>* <sup>+</sup> <sup>1</sup>)(*<sup>I</sup>* <sup>−</sup> <sup>1</sup>) *<sup>p</sup>*(*<sup>S</sup>* <sup>+</sup> 1, *<sup>I</sup>* <sup>−</sup> 1, *<sup>t</sup>*) + *<sup>γ</sup>*(*<sup>I</sup>* <sup>+</sup> <sup>1</sup>) *<sup>p</sup>*(*S*, *<sup>I</sup>* <sup>+</sup> 1, *<sup>t</sup>*) (21)

*<sup>N</sup>* <sup>+</sup> *<sup>γ</sup><sup>I</sup>* <sup>+</sup> *<sup>α</sup>*(*<sup>N</sup>* <sup>−</sup> *<sup>S</sup>* <sup>−</sup> *<sup>I</sup>*)

*β*

*w*)

√ *bt* 

+ *C*<sup>2</sup> cos

√ *bt u*2 ,

*p*(*S*, *I*,*t*) .

c)

<sup>2</sup> . In a) *β* = 1.01, and in b) *β* = 1.025. In c) we

Infected (β > γ)

55 60 65 70 75 80 85 90 95 100

Susceptible

a)

 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14

imaginary part *i*

*<sup>x</sup>*(*t*) = <sup>2</sup>*eat*

*dp*(*S*, *I*, *t*)

*dt* <sup>=</sup> *<sup>β</sup>*

Infected (β > γ)

99 99.2 99.4 99.6 99.8 100

eigendirection uλ1 eigendirection uλ2

Susceptible

general solution is then given by

*C*<sup>1</sup> cos

√ *bt* 

ODE system, see (Mattheij & Molenaar, 1996).

equation (van Kampen, 1992) in the following form

+ *α*(*N* − (*S* − 1) − *I*) *p*(*S* − 1, *I*,*t*) −

infected and susceptibles in time (Gillespie, 1976, 1978).

individuals toward the second fixed point *I*∗

b)

show the oscillatory trajectory toward the fixed point when *β* = 2.5

*x*(*t*) = *C*1*e*

When *λ*<sup>1</sup> *< λ*<sup>2</sup> the eigenvector *u*<sup>2</sup> is the driving force for *t* → ∞.

*t* → ∞, the eigenvector *u* is again the driving force (see Fig. 4a)).

fixed point, the general solution of the linearized system is given by

− *C*<sup>2</sup> sin

√ *bt u*<sup>1</sup> − *C*<sup>1</sup> sin

 0 0.05 0.1 0.15 0.2 0.25

Fig. 6. eigenvectors for the disease endemic equilibrium state in function of *β*. For a population *N* = 100, where *I* = 0.001 and *S* = 99, we fixed *α* = 0.1, *γ* = 1 and vary *β*. We plot the eigenvectors, *u*<sup>1</sup> (blue line) and *u*<sup>2</sup> (greenline), on top of the trajectory of the infected

*λt*

Infected (β > γ)

Fig. 7. Stochastic simulations for the basic SIR epidemic model. Here 10 realizations are plotted. We fixed *α* = 0.1, *γ* = 1 and *β* = 2.5. The deterministic trajectory is shown (pink line) top of the stochastic realizations for different population size *N* a. In a) *N* = 100, in b) *N* = 1000 and in c) *N* = 100000.

For mean values of infected �*I*� and susceptibles �*S*�, defined as e.g.

$$\langle I \rangle (t) := \sum\_{\mathcal{S}=0}^{N} \sum\_{I=0}^{N} I \ p(\mathcal{S}, I, t) \quad . \tag{22}$$

one can calculate the dynamics by inserting the master equation into the definition of the mean values obtaining

$$\begin{aligned} \frac{d}{dt} \langle S \rangle &= \alpha \langle R \rangle - \frac{\beta}{N} \, \langle SI \rangle \\\\ \frac{d}{dt} \langle I \rangle &= \frac{\beta}{N} \, \langle SI \rangle - \gamma \langle I \rangle \end{aligned} \tag{23}$$

with �*R*� = *N* − �*S*�−�*I*�. For more details on the calculations see e.g. (Stollenwerk & Jansen, 2010). These equations for the mean dynamics include now due to the nonlinear transition rates in the master equation also higher moments �*S* · *I*�. The simplest approximation to obtain a closed ODE system is to neglect cross-correlations �*S* · *I*�−�*S*�·�*I*� ≈ 0, the so-called mean field approximation (originally introduced for spatially extended systems in statistical physics (Stollenwerk et al., 2010)). Hence, the equation system (23) gives with identifying the higher moment �*S* · *I*� = �*S*�·�*I*� by a product of simple moments gives again the ODE system for SIR system, as it was just presented above. For certain parameter regions the mean field approximation describes the system well in terms of its mean dynamics and only small fluctuations around it. Then the previously shown analysis of the system is appropriate. However, noise can stabilize transients, a feature which becomes important in parameter regions where in the deterministic description a fixed point is reached via decreasing oscillations, as we have observed them in the SIR system. The noisy system would show here continued oscillations (Alonso et al., 2006).

In Fig. 7 we compare the deterministic and stochastic dynamics and we see that the magnitude of stochastic fluctuations decreases when the population size increases. However, the good approximation (see Fig. 7c)) is only achieved when the population size is large enough (see Fig. 7a) where most simulations die out very quickly for small population size).

Almost all mathematical models of diseases start from the same basic premise: that the population can be subdivided into a set of distinct classes. The most commonly used

of illness, and dengue hemorrhagic fever (DHF), which may evolve toward a severe form

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

Epidemiological studies support the association of DHF with secondary dengue infection (Guzmán et al., 2000; Halstead, 1982, 2003; Nisalak et al., 2003; Vaughn, 2000), and there is good evidence that sequential infection increases the risk of developing DHF, due to a process described as antibody-dependent enhancement (ADE), where the pre-existing antibodies to

Fig. 9. Scheme of the immunological response on recurrent dengue infections. In (a.) the first

antibodies (IgG class, the so called memory antibodies). In (e.) the temporary cross immunity period, that lasts between 3-9 months. After that period, the individual can get infected again with another dengue virus serotype (f.). In (g.) the IgG from the previous dengue infection binds to the new virus but do not inactivate them. In (h.) the complex antibody-virus enhances the new infection (i.). In (j.) the production of antibodies (IgM class) which is then able to inactivate the new viruses, leading to (l.), an enhanced immune response, such that

In the first dengue infection virus particles will be captured and processed by so-called antigen presenting cells. These viruses will be presented to T-cells causing them to become activated. And likewise B-cells will encounter their antigen free floating and become activated. B-cells produce antibodies that are used to tag the viruses to encourage their uptake by macrophages and inactivate them. In a secondary infection the antibodies from the first infection will attach to the virus particles but will not inactivate them. The antibody-virus complex suppresses innate immune responses, increasing intracellular infection and generating inflammatory citokines and chemokines that, collectively, result in enhanced disease (Dejnirattisai et al., 2010; Guzmán et al., 2010; Halstead, 1982, 1994, 2003; Mackenzie et al., 2004; WHO, 2009). Fig.9 is an scheme to illustrate the immunological response on recurrent dengue infections. DF is characterized by headache, retro-orbital pain, myalgia, arthralgia, rash, leukopenia, and mild thrombocytopenia. The symptoms resolve after 27 days. DHF is a potentially

infection with a given dengue virus serotype, in (b.) production of antibodies (Immunoglobulin M (IgM)), in (c.) inactivation of the virus and in (d.) production of

hemorrhagic symptoms are observed. In (m.) production of IgG antibodies.

previous dengue infection cannot neutralize but rather enhance the new infection.

known as dengue shock syndrome (DSS).

framework for epidemiological systems, is still the SIR type model, a good and simple model for many infectious diseases. However, different extensions of the classical single-strain SIR model show a rich dynamic behavior, e.g. (Stone et al., 2007) in measles, or in generalized multi-strain SIR type models to describe the epidemiology of dengue fever (Aguiar et al., 2008).
