Discrete Epidemic Models




(1)
Department of Mathematics, University of Florida, Gainesville, FL, USA

 




16.1 Single-Species Discrete Population Models


The continuous population models that we have considered in previous chapters model population and epidemic processes that occur continuously in time. In particular, they assume that births and deaths in the population occur continuously. This assumption is true for the human population, but many insect and plant populations have discrete, nonoverlapping generations. Such populations reproduce during specific time intervals of the year. Consequently, population censuses are taken at those specific times. As a result, modeling such populations and the distribution of disease in them should happen at discrete times. In this chapter we introduce discrete single-species population and epidemic models.


16.1.1 Simple Discrete Population Models


We assume that we measure the population at discrete, equally spaced, moments of time: 
$$ t_{0},t_{1},\ldots,t_{n},\ldots $$
, and we find that the population numbers at these moments of time are N t , where t takes the values of 
$$ t_{0},t_{1},\ldots,t_{n},\ldots $$
. For simplicity, we will set 
$$ N_{t_{n}} = N_{n} $$
. Thus, the population size is described by a sequence: 
$$ N_{1},N_{2},\ldots,N_{n},\ldots $$
. A discrete population model can be written in the following general form:



$$ \displaystyle{ N_{n+1} = \mathcal{F}(N_{n}), } $$

(16.1)
where 
$$ \mathcal{F} $$
is a specified function of N n . That is, if we know the population size at time t n , the model tells us what the populations size at time t n+1 should be. Such a model is equipped with a given initial condition: the population size N 0 at time t 0 is given. Another way to rewrite Eq. (16.1) is



$$ \displaystyle{ N_{n+1} = N_{n}f(N_{n}). } $$

(16.2)
The function f(N n ) is called a fitness function or per capita rate of population growth or net reproduction rate.


Definition 16.1.

Equations of the form (16.1) are called difference equations.

Such difference equations are of first order, because they contain only one time step. They are also autonomous, because 
$$ \mathcal{F} $$
does not depend explicitly on the time t n . The simplest discrete population model is derived under the assumption that individuals die with constant probability d. Furthermore, we assume that b individuals are born per individual in the population. The model then becomes



$$ \displaystyle{N_{n+1} = N_{n} + bN_{n} - dN_{n},} $$
that is, the number of individuals at the time step t n+1 is the number from the time step t n plus those who have been born, minus those who have died. Defining 
$$ \mathcal{R} = 1 + b - d $$
, we obtain the following linear discrete equation of population growth:



$$ \displaystyle{ N_{n+1} = \mathcal{R}N_{n}. } $$

(16.3)
The parameter 
$$ \mathcal{R} $$
is called the net reproduction number. We note that 
$$ \mathcal{R} > 0 $$
” src=”/wp-content/uploads/2016/11/A304573_1_En_16_Chapter_IEq9.gif”></SPAN>, since <SPAN class=EmphasisTypeItalic>b</SPAN> and <SPAN class=EmphasisTypeItalic>d</SPAN> are probabilities and are less than one. Model (<SPAN class=InternalRef><A href=16.3) is a discrete analogue of the Malthusian equation. Equation (16.3) is a special case of Eq. (16.2) with 
$$ f(N_{n}) = \mathcal{R} $$
. Model (16.3) can be solved. Given initial population size N 0, we have



$$ \displaystyle{ \begin{array}{l} N_{1} = \mathcal{R}N_{0}, \\ N_{2} = \mathcal{R}N_{1} = \mathcal{R}^{2}N_{0},\\ \vdots \\ N_{n} = \mathcal{R}N_{n-1} = \mathcal{R}^{n}N_{0}.\end{array} } $$

(16.4)
If 
$$ \mathcal{R} > 1 $$
” src=”/wp-content/uploads/2016/11/A304573_1_En_16_Chapter_IEq11.gif”></SPAN>, then each individual on average leaves more than one descendant, and the population grows geometrically. If <SPAN id=IEq12 class=InlineEquation><IMG alt=, then each individual leaves fewer than one descendant, and the population declines geometrically. If 
$$ \mathcal{R} = 1 $$
, the population remains constant. These model predictions are valid under the assumption that the resources are unlimited.

In practice, populations do not experience unlimited growth, so models that predict asymptotically bounded growth are more realistic. One such model is the discrete analogue of the logistic equation. To derive such an analogue, we approximate the continuous time derivative with N n+1N n , assuming that the time step is equal to one. Thus the discrete logistic equation takes the form



$$ \displaystyle{ N_{n+1} = N_{n} + rN_{n}\left (1 -\frac{N_{n}} {K} \right ). } $$

(16.5)
First we factor N n and r + 1. Furthermore, we make the following changes in dependent variables and parameters:



$$ \displaystyle{y_{n} = \frac{r} {r + 1} \frac{N_{n}} {K} \qquad a = r + 1.} $$
We obtain a classical form for the discrete logistic equation:



$$ \displaystyle{y_{n+1} = ay_{n}(1 - y_{n}).} $$
This method for producing discrete equations is not foolproof, however. The discrete logistic equation above is not well posed, in the sense that its solutions can become negative. This is not hard to see. Suppose we start from y 0 = 0. 5 and a = 6. Then y 1 = 1. 5. Consequently, y 2 < 0. Thus, the logistic equation is not a very good discrete population model.

We can derive a discrete version of the simplified logistic model. Suppose the population increases in each time interval by a constant amount 
$$ \varLambda $$
, and that γ ≤ 1 is the probability for survival of individuals to the next time period. Then the simplified logistic model takes the form



$$ \displaystyle{ N_{n+1} =\varLambda +\gamma N_{n}. } $$

(16.6)
This model can also be solved explicitly:



$$ \displaystyle{ \begin{array}{l} N_{1} =\varLambda +\gamma N_{0}, \\ N_{2} =\varLambda +\gamma (\varLambda +\gamma N_{0}) =\gamma ^{2}N_{0} + (1+\gamma )\varLambda,\\ \vdots \\ N_{n} =\gamma ^{n}N_{0} + (1 +\gamma +\ldots +\gamma ^{n-1})\varLambda.\end{array} } $$

(16.7)
Hence,



$$ \displaystyle{ N_{n} = \left \{\begin{array}{ll} N_{0} +\varLambda n, &\qquad \gamma = 1, \\ \gamma ^{n}\left (N_{ 0} - \frac{\varLambda } {1-\gamma }\right ) + \frac{\varLambda } {1-\gamma },&\qquad \gamma < 1. \end{array} \right. } $$

(16.8)

Other discrete population models have been proposed that guarantee that the population remains positive for all times. One such model, proposed by Bill Ricker [138], is the Ricker model:



$$ \displaystyle{ N_{n+1} = N_{n}e^{r\left (1-\frac{N_{n}} {K} \right )}. } $$

(16.9)

Another model also widely used is the Beverton–Holt model [23], also called the Verhulst equation:



$$ \displaystyle{ N_{n+1} = \frac{rN_{n}} {A + N_{n}}. } $$

(16.10)

A generalization of the Beverton–Holt model can be made that is known as the Hassell equation [72]:



$$ \displaystyle{ N_{n+1} = \frac{rN_{n}} {(A + N_{n})^{b}}, } $$

(16.11)
where b > 0 is a positive parameter.


16.1.2 Analysis of Single-Species Discrete Models


Difference equations also have solutions that do not depend on time, called equilibria. Since the solution does not depend on time, all members of the sequence have the same value, that is, we have



$$ \displaystyle{N_{n} = N^{{\ast}}\qquad \qquad \mathrm{for}\quad \mathrm{all}\quad n \geq 0.} $$
Consequently, equilibria of the difference equation (16.1) must satisfy 
$$ N^{{\ast}} = \mathcal{F}(N^{{\ast}}) $$
.


Definition 16.2.

A value N that satisfies



$$ \displaystyle{N^{{\ast}} = \mathcal{F}(N^{{\ast}})} $$
is called a fixed point of the function 
$$ \mathcal{F} $$
.


Example 16.1.

Consider the equilibria of the logistic equation



$$ \displaystyle{N^{{\ast}} = N^{{\ast}} + rN^{{\ast}}\left (1 -\frac{N^{{\ast}}} {K} \right ).} $$
The solutions of this equation are N 1  = 0 and N 2  = K, that is, the equilibria in the discrete case are exactly the same as in the continuous case. The equilibrium N 1  = 0 is called a trivial equilibrium, while the equilibrium N 2  = K is called a nontrivial equilibrium.

To describe the behavior of the solutions near an equilibrium, we use again a process called linearization. Let N be the equilibrium, and u n the perturbation of the solution from the equilibrium, that is,



$$ \displaystyle{N_{n} = N^{{\ast}} + u_{ n}.} $$
Substituting this equation into Eq. (16.1), we have 
$$ u_{n+1} + N^{{\ast}} = \mathcal{F}(u_{n} + N^{{\ast}}) $$
. Expanding 
$$ \mathcal{F} $$
in a Taylor series and neglecting all terms containing powers of u n greater than one, we obtain



$$ \displaystyle{u_{n+1} + N^{{\ast}} = \mathcal{F}(N^{{\ast}}) + \mathcal{F}'(N^{{\ast}})u_{ n}.} $$
Recall that since N is an equilibrium, we have 
$$ N^{{\ast}} = \mathcal{F}(N^{{\ast}}) $$
. Hence, we obtain the following linearized equation:



$$ \displaystyle{ u_{n+1} = \mathcal{F}'(N^{{\ast}})u_{ n}. } $$

(16.12)
We note that 
$$ \mathcal{F}'(N^{{\ast}}) $$
is a fixed number, which may be positive or negative. If we consider



$$ \displaystyle{ u_{n+1} = \vert \mathcal{F}'(N^{{\ast}})\vert u_{ n}, } $$

(16.13)

then Eq. (16.13) is exactly the discrete Malthus equation. Consequently, we have the following:

1.

If 
$$ \vert \mathcal{F}'(N^{{\ast}})\vert < 1 $$
, then u n  → 0. Hence, N n N  → 0 and N n  → N . This is the case if N 0 is close enough to N , that is, this result is local. In this case, we call N  locally asymptotically stable.

 

2.

If 
$$ \vert \mathcal{F}'(N^{{\ast}})\vert > 1 $$
” src=”/wp-content/uploads/2016/11/A304573_1_En_16_Chapter_IEq22.gif”></SPAN>, then <SPAN id=IEq23 class=InlineEquation><IMG alt=. Hence 
$$ N_{n} - N^{{\ast}}\rightarrow \infty  $$
, and N n diverges from N . This is the case if N 0 is close enough to N . In this case, we call N  unstable.

 

We note that if 
$$ \vert \mathcal{F}'(N^{{\ast}})\vert = 1 $$
, we cannot draw conclusions from the local analysis.

We summarize the above discussion in the following theorem:


Theorem 16.1.

The equilibrium N  of the discrete equation(16.1)is locally asymptotically stable if and only if 
$$ \vert \mathcal{F}'(N^{{\ast}})\vert < 1 $$
. The equilibrium N  of the discrete equation (16.1) is unstable if and only if 
$$ \vert \mathcal{F}'(N^{{\ast}})\vert > 1 $$
” src=”/wp-content/uploads/2016/11/A304573_1_En_16_Chapter_IEq27.gif”></SPAN> <SPAN class=EmphasisTypeItalic>.</SPAN> </DIV></DIV><br />
<DIV class=Para>To illustrate the use of the theorem above, we consider the local stability of the equilibria of the logistic equation.</DIV><br />
<DIV id=FPar5 class=
Example 16.2.

In the case of the logistic equation (16.5), the function 
$$ \mathcal{F} $$
is given by



$$ \displaystyle{\mathcal{F}(N) = N + rN\left (1 -\frac{N} {K}\right ).} $$
The derivative is given by



$$ \displaystyle{\mathcal{F}'(N) = 1 + r\left (1 -\frac{N} {K}\right ) - \frac{r} {K}N.} $$
In the case of the trivial equilibrium N  = 0, we have



$$ \displaystyle{\mathcal{F}'(0) = 1 + r > 1.} $$
” src=”/wp-content/uploads/2016/11/A304573_1_En_16_Chapter_Equk.gif”></DIV></DIV></DIV>Consequently, the trivial equilibrium is always unstable. Now we consider the nontrivial equilibrium <SPAN class=EmphasisTypeItalic>N</SPAN> <SUP>∗</SUP> = <SPAN class=EmphasisTypeItalic>K</SPAN>. We have<br />
<DIV id=Equl class=Equation><br />
<DIV class=EquationContent><br />
<DIV class=MediaObject><IMG alt=
So if | 1 − r |  < 1, or equivalently, if 0 < r < 2, then the nontrivial equilibrium is locally asymptotically stable.

When r > 2, simulations suggests that the logistic equation can experience very complex behavior. To investigate this behavior through simulations, we will study the nondimensionalized version of the logistic equation:



$$ \displaystyle{ y_{n+1} =\rho y_{n}(1 - y_{n}). } $$

(16.14)
Recall that ρ = 1 + r, so we can expect complex behavior for ρ > 3. We notice that the corresponding equilibria of the nondimensional logistic model are y  = 0 and y  = 1. The first complexity that appears is a 2-cycle.


Definition 16.3.

A 2-cycle of model (16.1) is a system of two solutions y 1 and y 2 such that



$$ \displaystyle{ \begin{array}{l} y_{1} = \mathcal{F}(y_{2}), \\ y_{2} = \mathcal{F}(y_{1}).\\ \end{array} } $$

(16.15)

In model (16.14), 
$$ \mathcal{F}(y) =\rho y(1 - y) $$
. As ρ increases, the system experiences a process, called period-doubling, to a 4-cycle. Similarly, a 4-cycle of model (16.14) is a system of four solutions y 1, y 2, y 2, y 4 such that



$$ \displaystyle{ \begin{array}{l} y_{1} = \mathcal{F}(y_{4}), \\ y_{2} = \mathcal{F}(y_{1}), \\ y_{3} = \mathcal{F}(y_{2}), \\ y_{4} = \mathcal{F}(y_{3}). \end{array} } $$

(16.16)
Further period-doubling occurs to an 8-cycle. The period-doubling continues until the system begins to exhibit chaos. We illustrate period-doubling and chaos in Fig. 16.1.

A304573_1_En_16_Fig1_HTML.gif


Fig. 16.1
The figure shows time series of the model y n+1 = ρ y n (1 − y n ) for different values of ρ. The first figure shows a two-cycle with ρ = 3. 43. The second figure on the top line shows a 4-cycle with ρ = 3. 47. The left figure on the bottom row shows an 8-cycle with ρ = 3. 58. The right figure on the bottom row shows chaos with ρ = 3. 7

We need single-species discrete population models to capture the demographic processes in epidemic models. Many books focus on single-species discrete models and provide an excellent introduction to these models (for instance, see [27, 90]).


16.2 Discrete Epidemic Models


Just like single-species population models, discrete epidemic models can also be obtained from a discretization of the continuous epidemic models. However, this approach results in models that have issues like those of the discrete logistic equation. To avoid these problems, a modeling approach specific to discrete models should be taken. We follow here the approach of Castillo-Chavez and Yakubu [39].


16.2.1 A Discrete SIS Epidemic Model


We begin with a general population model



$$ \displaystyle{ N_{n+1} = f(N_{n}) +\gamma N_{n}, } $$

(16.17)
where γ < 1 is the probability of survival to the next time period, and f(N n ) is a recruitment function. We assume that the disease does not affect the population dynamics, that is, we assume that the disease is nonfatal and does not affect the birth process. We will build an SIS epidemic process on top of the demographic process. We denote by S n and I n the susceptible and infected individuals at time t n . Individuals survive with probability γ < 1 (die with probability 1 −γ) in each generation. Infected individuals recover with probability 
$$ 1-\sigma $$
(do not recover with probability 
$$ \sigma < 1 $$
) in each generation. In each generation, susceptible individuals become infected with probability 1 − G (remain susceptible with probability G). The function G is a function of the prevalence I n N n , which is weighted with coefficient α. The model assumes a sequential process: at each generation, γ S n susceptibles survive, and the surviving susceptibles become infected with probability 1 − G. Similarly, γ I n infected individuals survive, and the surviving ones recover with probability 
$$ (1-\sigma ) $$
:



$$ \displaystyle\begin{array}{rcl} S_{n+1}& =& f(N_{n}) +\gamma S_{n}G\left ( \frac{\alpha I_{n}} {N_{n}}\right ) +\gamma (1-\sigma )I_{n}, \\ I_{n+1}& =& \gamma S_{n}\left (1 - G\left ( \frac{\alpha I_{n}} {N_{n}}\right )\right ) +\gamma \sigma I_{n}. {}\end{array} $$

(16.18)
The function G must satisfy the following conditions:

1.


$$ G: [0,\infty ) \rightarrow [0,1] $$
.

 

2.

G(0) = 1.

 

3.

G is a monotone decreasing function with G′(x) < 0 and G″(x) ≥ 0.

 

An example of such a function that we will use is G(x) = e x . Another example is G(x) = A∕(x + A). Adding the two equations in system (16.18) gives Eq. (16.17).


16.2.2 Analysis of Multidimensional Discrete Models


In this subsection, we introduce the techniques that help us analyze systems of discrete equations. Suppose that we are given the following system:



$$ \displaystyle{ \mathbf{x}_{n+1} = \mathcal{F}(\mathbf{x}_{n}), } $$

(16.19)
where x is an M-dimensional vector of variables. As before, an equilibrium of system (16.19) is the solution of the problem



$$ \displaystyle{\mathbf{x}^{{\ast}} = \mathcal{F}(\mathbf{x}^{{\ast}}).} $$
To find the behavior of the solutions near an equilibrium, we use linearization. We set x n  = u n +x . We obtain the following linear system:



$$ \displaystyle{ \mathbf{u}_{n+1} = \mathbf{J}(\mathbf{x}^{{\ast}})\mathbf{u}_{ n}, } $$

(16.20)
where J is the Jacobian of the system, that is,



$$ \displaystyle{ \mathbf{J}(\mathbf{x}^{{\ast}}) = \left (\begin{array}{ccc} \frac{\partial \mathcal{F}_{1}} {\partial x_{1}} & \ldots \frac{\partial \mathcal{F}_{1}} {\partial x_{M}}\\ & \vdots & \\ \frac{\partial \mathcal{F}_{M}} {\partial x_{1}} & \ldots \frac{\partial \mathcal{F}_{M}} {\partial x_{M}} \end{array} \right )\vert _{\mathbf{x}=\mathbf{x}^{{\ast}}}. } $$

(16.21)


Definition 16.4.

An equilibrium point x is said to be locally asymptotically stable if there exists a neighborhood U of x such that for each starting value x 0 ∈ U, we get



$$ \displaystyle{ \lim _{n\rightarrow \infty }\mathbf{x}_{n} = \mathbf{x}^{{\ast}}. } $$

(16.22)
The equilibrium point x is called unstable if x is not (locally asymptotically) stable.

The limit (16.22) holds if for system (16.20), we have 
$$ \lim _{n\rightarrow \infty }\mathbf{u}_{n} = 0 $$
. The following theorem gives the conditions for convergence of solutions of the linear system (16.20) to zero:


Theorem 16.2.

Let J be an M × M matrix with ρ( J ) < 1, where 



$$ \displaystyle{\rho (\mathbf{J}) =\max \{ \vert \lambda \vert:\lambda \mathrm{ \,\,is\,\,an\,\,eigenvalue\,\,of}\,\,\mathbf{J}\}.} $$
 Then every solution of ( 16.20 ) satisfies



$$ \displaystyle{\lim _{n\rightarrow \infty }\mathbf{u}_{n} = 0.} $$
If ρ( J ) > 1, then there are solutions that go to infinity.

This implies the following criterion for stability of an equilibrium x of system (16.19).


Theorem 16.3.

Consider the nonlinear autonomous system ( 16.19 ). Suppose 
$$ \mathcal{F}: \mathcal{D}\rightarrow \mathcal{D} $$
, where 
$$ \mathcal{D}\subset \mathbf{R}^{M} $$
 and 
$$ \mathcal{D} $$
 is an open set. Suppose 
$$ \mathcal{F} $$
 is twice continuously differentiable in some neighborhood of a fixed point 
$$ \mathbf{x}^{{\ast}} \in \mathcal{D} $$
. Let J ( x ) be the Jacobian matrix of 
$$ \mathcal{F} $$
 evaluated at x . Then the following hold:

1.

x is locally asymptotically stable if all eigenvalues of J ( x ) have magnitude less than one.

 

2.

x is unstable if at least one eigenvalue of J ( x ) has magnitude greater than one.

 

The Routh–Hurwitz criterion will not be helpful here in determining which matrices are stable, since Routh–Hurwitz identifies matrices whose eigenvalues lie in the left half of the complex plane. However, there is an analogous criterion that can help determine whether the spectral radius of a matrix is smaller than one. This criterion is called the Jury conditions. Let



$$ \displaystyle{p(\lambda ) = \vert \mathbf{J} -\lambda I\vert = a_{M}\lambda ^{M} +\ldots +a_{ 1}\lambda + a_{0},} $$
where a M  = 1. To introduce the Jury conditions, we first have to introduce the Jury array. The Jury array is composed as follows: First we write out a row of the coefficients, and then we write out another row with the same coefficients in reverse order. The first two rows of the Jury array are composed of the coefficients of the polynomial 
$$ p(\lambda ) $$
above. Once we have the first two rows of the a coefficients, the next two rows are of the b coefficients, and so on. We obtain the array of Table 16.1, where the b coefficients, c coefficients, etc., are composed as follows:


Table 16.1
Jury array









































Number

Coeff.

Coeff.

Coeff.

Coeff.

Coeff.

(1)

a 0

a 1


a M−1

a M

(2)

a M

a M−1


a 1

a 0

(3)

b 0

b 1


b M−1
 

(4)

Only gold members can continue reading. Log In or Register to continue

Stay updated, free articles. Join our Telegram channel

Nov 20, 2016 | Posted by in PUBLIC HEALTH AND EPIDEMIOLOGY | Comments Off on Discrete Epidemic Models

Full access? Get Clinical Tree

Get Clinical Tree app for offline access