Age-Structured Epidemic Models




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

 




12.1 Introduction


Chronological age is perhaps one of the most important factors distinguishing individuals in a population that needs to be incorporated in population and epidemic models. Undoubtedly, vital characteristics such as birth and death rates differ markedly among the individuals of various ages. Age is also a key to capturing important mixing patterns in epidemic models. For instance, in childhood diseases, children predominantly mix with other children in similar age groups as well as with the individuals of the age groups of their parents and grandparents. Children are at greatest risk for contracting malaria and exhibiting strongest symptoms and highest death rate, yet malaria affects all age groups. Incidence of HIV is highest in the age groups from age twenty to age forty-five. Endemic models that incorporate births and deaths should also preferably incorporate age structure, since with time, the age profile of the population may change, and that may affect the dynamics of the disease.

Age-structured epidemic models are built on the basis of age-structured population models. There are several excellent introductory texts for age-structured population models [78, 46]. For completeness, we will introduce here first the linear age-structured population model. Most epidemic models use the linear model as a baseline population model.


12.2 Linear Age-Structured Population Model


The Malthusian model, which is also linear, considers a homogeneous population in which individuals are not distinguished by age [104]. The linear age-structured model is a strict analogue of the Malthusian model but allows for variability in age. As in the Malthusian model, the age-structured model considers a single population, not stratified by sex, which is isolated (no emigration and immigration) and lives in an invariant habitat, that is, the birth and death rates can be assumed time-independent. The only characteristics by which individuals in the population differ is age.

Models that structure the population by age can be continuous or discrete. In this chapter, we will consider the continuous age-structured model. Because age and time are considered two independent variables, the resulting continuous models are cast as partial differential equation models. It was the physician A.G. McKendrick who first considered an age-structured PDE model for the dynamics of a one-sex population [115]. Earlier age-structured models were developed by Sharpe and Lotka [144, 3] but they were formulated as integral equations of Volterra type. We will first introduce McKendrick’s age-structured PDE model, and then derive Lotka’s model.


12.2.1 Derivation of the Age-Structured Model


McKendrick model introduces two independent variables: age a and time t. The distribution of the female population can be described by a function u(a, t) that denotes the age-density of individuals at age a at time t. We assume a ∈ [0, a †], t ≥ 0, where $$a_{\dag }$$ denotes the maximal age, which may be assumed infinite. If $$a_{\dag } = \infty$$ , then we assume that u(a, t) = 0 for all sufficiently large values of a. The function gives “density” rather than numbers, because the number of individuals of ages between a and a +Δ A, where Δ a is a small increment, at time t is approximately u(a, t)Δ a. Thus, the total number of individuals of ages in the interval [a 1, a 2] at time t is given by the integral


$$\displaystyle{\int _{a_{1}}^{a_{2} }u(a,t)da.}$$
Respectively, the total number of individuals in the population or the total population size at time t is given by


$$\displaystyle{P(t) =\int _{ 0}^{a_{\dag } }u(a,t)da.}$$
Consider a cohort of individuals, that is, a group of individuals of age in an interval of length Δ a. The number of individuals in that cohort is u(a, t)Δ a. If a small interval of time Δ t elapses, then those individuals who where of age a at time t will be of age a +Δ t, and time is t +Δ t. The number of these individuals is u(a +Δ t, t +Δ t)Δ a. This is the same cohort. Their new number, however, is smaller than their original number, since some of these individuals will have died. We assume that the members of the population leave the population only by death. We denote the age-specific per capita mortality rate by μ(a). Then the number of individuals that die at age [a, a +Δ a] at time t is μ(a)u(a, t)Δ a. For the whole interval of time Δ t, that number would be $$\mu (a)u(a,t)\varDelta a\varDelta t$$ . The balance law can be written as


$$\displaystyle{ u(a +\varDelta t,t +\varDelta t)\varDelta a - u(a,t)\varDelta a = -\mu (a)\varDelta tu(a,t)\varDelta a. }$$

(12.1)

Dividing both sides by Δ a Δ t, we obtain a difference quotient on the left-hand side:


$$\displaystyle{ \frac{u(a +\varDelta t,t +\varDelta t) - u(a,t)} {\varDelta t} = -\mu (a)u(a,t). }$$

(12.2)
We take the limit as t → 0. If u(a, t) is a differentiable function with respect to each variable, we have


$$\displaystyle\begin{array}{rcl} & & \lim _{\varDelta t\rightarrow 0}\frac{u(a +\varDelta t,t +\varDelta t) - u(a,t)} {\varDelta t} \\ & & \quad =\lim _{\varDelta t\rightarrow 0}\frac{u(a +\varDelta t,t +\varDelta t) - u(a,t +\varDelta t)} {\varDelta t} \\ & & \qquad +\lim _{\varDelta t\rightarrow 0}\frac{u(a,t +\varDelta t) - u(a,t)} {\varDelta t} \\ & & \quad = u_{a}(a,t) + u_{t}(a,t). {}\end{array}$$

(12.3)
We obtain McKendrick’s equation, also referred to as the McKendrick–von Foerster equation, since it was later derived by von Foerster (1959) in the context of cell biology [63]:


$$\displaystyle{ u_{a}(a,t) + u_{t}(a,t) = -\mu (a)u(a,t). }$$

(12.4)

The mortality rate is connected to the probability of survival (Fig. 12.1). Let π(a) be the probability of survival from birth to age a. Suppose we start from a cohort of newborns of size P. Then π(a)P gives the number of individuals from the cohort who are of age a. Similarly, π(a +Δ a)P gives the number of individuals from the cohort who are of age a +Δ a. Then π(a +Δ a)P −π(a)P is the number of individuals who have died in the age interval Δ a. We have


$$\displaystyle{\pi (a +\varDelta a)P -\pi (a)P = -\mu (a)\pi (a)P\varDelta a.}$$
Canceling P and dividing by Δ a, we have


$$\displaystyle{\frac{\pi (a +\varDelta a) -\pi (a)} {\varDelta a} = -\mu (a)\pi (a).}$$
Taking the limit as Δ a → 0 and assuming that the probability of survival π(a) is a differentiable function, we obtain the following ordinary differential equation for π(a):

A304573_1_En_12_Fig1_HTML.gif


Fig. 12.1
Data on age-specific probability of survival in the United States in 2007. The data are taken from the life tables in the U.S. Vital Statistics Reports [15]. Data are interpolated to show the continuous curve




$$\displaystyle{\pi '(a) = -\mu (a)\pi (a).}$$
Since π(0) is the probability of survival until age 0, we may assume that π(0) = 1. Solving the initial value problem, we have the following constitutive form for the probability of survival until age a:


$$\displaystyle{\pi (a) = e^{-\int _{0}^{a}\mu (s)ds }.}$$
We further assume that the probability of survival until age a 1 is independent of the probability of survival until age a 2. If a 2 > a 1, we have


$$\displaystyle{ \pi (a_{2}) =\pi (a_{1})e^{-\int _{a_{1}}^{a_{2}}\mu (s)ds }. }$$

(12.5)
Thus, $$e^{-\int _{a_{1}}^{a_{2}}\mu (s)ds }$$ is the probability that an individual will survive from age a 1 to age a 2. Age-dependent mortality μ(a) can be found in the national vital statistics files [121], but it is advisable that it be estimated from the life tables and the probability of survival. From the formula (12.5), we have


$$\displaystyle{\int _{a_{1}}^{a_{2} }\mu (s)ds =\ln \left (\frac{\pi (a_{1})} {\pi (a_{2})}\right ).}$$
Hence approximately [12],


$$\displaystyle{\mu (a_{1}) = \frac{1} {a_{2} - a_{1}}\ln \left (\frac{\pi (a_{1})} {\pi (a_{2})}\right ).}$$
Computing the discrete values of μ(a), we may obtain a typical form of the age-specific mortality function. The age-specific mortality data for the United States in 2007 and a function μ(a) that fits them are given in Fig. 12.2.

A304573_1_En_12_Fig2_HTML.gif


Fig. 12.2
Data on age-specific mortality in the United States in 2007. The data are taken from life tables of the the U.S. Vital Statistics Reports [15] as explained in the text. The age-specific function fitted to the data is $$\mu (a) = \frac{0.748123a} {(110 - a)^{2}}$$


Remark 12.1.

If $$a_{\dag } <\infty$$ , then we must have


$$\displaystyle{\lim _{a\rightarrow a_{\dag }}\pi (a) = 0,}$$
which means that no one survives until the maximal age. That, in turn, implies that


$$\displaystyle{\lim _{a\rightarrow a_{\dag }}\mu (a) = \infty.}$$

Equation (12.4) is a first-order linear partial differential equation. It is defined on the domain


$$\displaystyle{\mathcal{D} =\{ (a,t): a \geq 0,t \geq 0\},}$$
that is, in the first quadrant. We will need conditions on the boundary of the domain to complete the model. In particular, we need to specify the age density for the initial population distribution at time t = 0:


$$\displaystyle{u(a,0) = u_{0}(a),}$$
where u 0(a) is a given function. The function u 0(a) is called the initial population density. It is assumed that u 0(a) ≥ 0. Furthermore,


$$\displaystyle{P_{0} =\int _{ 0}^{a_{\dag } }u_{0}(a)da <\infty,}$$
where P 0 is the initial total population size.

The value u(0, t) is the number of newborns at time t. To model the birth process, we introduce the age-specific per capita birth rate β(a). Since there are u(a, t)Δ a females in the population with ages in the interval [a, a +Δ a], it follows that β(a)u(a, t)Δ a gives the number of births to females of age [a, a +Δ a] at time t. Therefore, the total number of births at time t is the sum of all births at all ages: $$\sum _{i}\beta (a_{i})u(a_{i},t)\varDelta a$$ . Taking the limit as Δ a → 0, we obtain an integral in place of the sum. Thus, the total birth rate is given by


$$\displaystyle{B(t) =\int _{ 0}^{a_{\dag } }\beta (a)u(a,t)da,}$$
which also gives the total number of newborns at time t. That is, we have


$$\displaystyle{u(0,t) =\int _{ 0}^{a_{\dag } }\beta (a)u(a,t)da.}$$

Data on fertility for the United States in 2010 and a typical form of birth rate are given in Fig. 12.3.

A304573_1_En_12_Fig3_HTML.gif


Fig. 12.3
Data on age-specific fertility rate in the United States in 2010. The data are taken from the U.S. Vital Statistics Reports [70]. The age-specific function fitted to the data is $$\beta (a) = 8.151 {\ast} 10^{-7}(a - 4.48268)^{10}e^{-0.459a}$$

The full McKendrick–von Foerster age-structured population model takes the form


$$\displaystyle{ \left \{\begin{array}{l} u_{a}(a,t) + u_{t}(a,t) = -\mu (a)u(a,t), \\ u(0,t) =\int _{ 0}^{a_{\dag } }\beta (a)u(a,t)da, \\ u(a,0) = u_{0}(a).\end{array} \right. }$$

(12.6)
The condition specified on the boundary a = 0, that is, u(0, t), is called a boundary condition. Since the boundary condition is not specified through a given

function but through an equation that depends on the entire unknown function u(a, t), this boundary condition is referred to as a nonlocal boundary condition. The condition giving a value to u(a, 0) is called an initial condition.


12.2.2 Reformulation of the Model Through the Method of the Characteristics. The Renewal Equation


Model (12.6) can be recast in an integral equation form, a form first derived by Sharpe and Lotka [3, 144]. To obtain the integral equation form, we must integrate the partial differential equation in model (12.6). That can be done through the method of characteristics. The method of characteristics typically applies to first-order hyperbolic partial differential equations. The method identifies curves, called characteristics, along which the partial differential equation reduces to an ordinary differential equation. Then the ordinary differential equation can be integrated from some initial data given on a suitable curve. Finally, the solution of the ordinary differential equation can be transformed into a solution for the original PDE.

An important step in the method of characteristics is identifying the characteristics of the PDE. This is relatively simple for almost linear PDEs. For an almost linear PDE of the form


$$\displaystyle{c_{1}(x,y)\frac{\partial u} {\partial x} + c_{2}(x,y)\frac{\partial u} {\partial y} = c_{0}(x,y,u)}$$
in parameterized form, the characteristics are given by


$$\displaystyle\begin{array}{rcl} \frac{dx} {ds}& =& c_{1}(x,y), \\ \frac{dy} {ds}& =& c_{2}(x,y).{}\end{array}$$

(12.7)
Therefore, in our case we have


$$\displaystyle{\frac{da} {ds} = 1,\qquad \frac{dt} {ds} = 1.}$$
Dividing the one of the equations by the other, we have dt∕da = 1. Hence t = a + c, where c is an arbitrary constant. This implies that the characteristics are lines of slope 1, called characteristic lines of the PDE. What is the practical significance of the characteristic lines? The value of u(a, t) along the characteristic lines models one cohort of individuals. Thus the value of u(a, t) is determined by previous values of u(a, t) along the characteristics, and in particular by the value of u where the characteristic line crosses the boundary of $$\mathcal{D}$$ .

The integral formulation can be derived from the partial differential equation through a rigorous procedure, called integration along the characteristic lines. To apply the procedure, we fix a point (a 0, t 0) in the first quadrant. We parameterize the characteristic line that goes through that point. If we denote by s the parameter for the fixed point (a 0, t 0) and a variable s, then u(a 0 + s, t 0 + s) gives the value of u along the characteristic line. Define v(s) = u(a 0 + s, t 0 + s). Then the derivative along the characteristic line can be written as


$$\displaystyle{\frac{dv} {ds} = u_{a} + u_{t}.}$$
We let also $$\bar{\mu }(s) =\mu (a_{0} + s)$$ . Then the partial differential equation in model (12.6) can be written as the following ordinary differential equation along the characteristic line:


$$\displaystyle{\frac{dv} {ds} = -\bar{\mu }(s)v(s).}$$
This ODE can easily be solved to give


$$\displaystyle{ v(s) = v(0)e^{-\int _{0}^{s}\bar{\mu }(\tau )d\tau }. }$$

(12.8)
Now we have to interpret this solution in terms of u and two variables a and t. We consider two cases:

1.

a 0 ≥ t 0. The solution (12.8) gives, in the original variables,


$$\displaystyle{ u(a_{0} + s,t_{0} + s) = u(a_{0},t_{0})e^{-\int _{0}^{s}\mu (a_{ 0}+\tau )d\tau }. }$$

(12.9)
Now we specify our original point to have been fixed on one of the boundaries. Since we are in the case a 0 ≥ t 0, we must take t 0 = 0. Then s = t and a 0 = a − t. The choice of s and a 0 is made in such a way that a 0 + s = a and t 0 + s = t. With these specifications, from (12.9) we obtain the solution in the case a 0 ≥ t 0, that is, the solution when a ≥ t:


$$\displaystyle\begin{array}{rcl} u(a,t)& =& u(a - t,0)e^{-\int _{0}^{t}\mu (a-t+\tau )d\tau }, \\ & =& u(a - t,0)e^{-\int _{a-t}^{a}\mu (\sigma )d\sigma }, \\ & =& u_{0}(a - t) \frac{\pi (a)} {\pi (a - t)}.{}\end{array}$$

(12.10)

 

2.

a 0 < t 0. Again


$$\displaystyle{ u(a_{0} + s,t_{0} + s) = u(a_{0},t_{0})e^{-\int _{0}^{s}\mu (a_{ 0}+\tau )d\tau }. }$$

(12.11)
Now we set a 0 = 0. Then s = a, t 0 = t − a. Substituting above, we get for a < t,


$$\displaystyle{u(a,t) = u(0,t - a)e^{-\int _{0}^{a}\mu (\tau )d\tau } = B(t - a)\pi (a).}$$

 

So finally, the solution of the PDE can be written in the form




$$\displaystyle{ u(a,t) = \left \{\begin{array}{ll} u_{0}(a - t) \frac{\pi (a)} {\pi (a - t)} &~~~~~~a \geq t, \\ B(t - a)\pi (a) &~~~~~~a <t. \end{array} \right. }$$

(12.12)

This would have been an explicit solution of the partial differential equation if B(t) were a given function. However, B(t) depends on u(a, t). Thus, this representation is another equation for u(a, t). This representation is not easy to work with, since it is an equation for a function of two variables. It is customary to rewrite this representation in the form of an equation for the function B(t). To do that, we substitute the expression for u(a, t) from (12.12) into the formula for B(t). We have two cases:

Case $$a_{\dag } <\infty$$ :

In this case, we have


$$\displaystyle{ B(t) = \left \{\begin{array}{ll} \int _{0}^{t}\beta (a)\pi (a)B(t - a)da +\int _{ t}^{a_{\dag } }\beta (a) \frac{\pi (a)} {\pi (a - t)}u_{0}(a - t)da \\ \qquad t \leq a_{\dag }, \\ \int _{0}^{a_{\dag } }\beta (a)\pi (a)B(t - a)da t> a_{\dag }.\end{array} \right. }$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_12_Chapter_Equ13.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(12.13)</DIV></DIV></DIV></DIV></DIV><br /> <DIV class=DefinitionListEntry><SPAN class=Term>Case <SPAN id=IEq13 class=InlineEquation><IMG alt=:

In this case, we have only $$t <\infty$$ . Therefore, we have the following equation for B(t):




$$\displaystyle{ B(t) =\int _{ 0}^{t}\beta (a)\pi (a)B(t - a)da + F(t), }$$

(12.14)

where F(t) is a given function:


$$\displaystyle{F(t) =\int _{ t}^{\infty }\beta (a) \frac{\pi (a)} {\pi (a - t)}u_{0}(a - t)da =\int _{ 0}^{\infty }\beta (a + t)\frac{\pi (a + t)} {\pi (a)} u_{0}(a)da.}$$

Equation (12.14) is a linear Volterra integral equation of convolution type with kernel K(a) = β(a)π(a). The function K(a) is sometimes referred to as a maternity function [78]. Equation (12.14) is called the renewal equation or Lotka equation. If we can solve (12.14) and determine B(t) from it, then we can obtain u(a, t) from (12.12).


12.2.3 Separable Solutions. Asymptotic Behavior


The McKendrick–von Foerster model (12.6) is a linear model. Thus its solutions are not necessarily bounded but grow or decay to zero nearly exponentially in time. It can be shown that all solutions approach in time a separable solution. Separable solutions for partial differential equation are solutions that can be written as a product of two functions: a function of age and a function of time, that is, the solution can be written in the form u(a, t) = T(t)ϕ(a). It can be shown (see Problem 12.4) that the time-dependent function T(t) is actually an exponential function. Hence, the separable solutions of the model (12.6) are of the form


$$\displaystyle{u(a,t) = e^{\lambda t}\phi (a).}$$
To find the separable solutions, we have to find the number $$\lambda$$ and the function ϕ(a). To find $$\lambda$$ and ϕ(a), we substitute in the differential equation of model (12.6). Therefore, we have


$$\displaystyle{e^{\lambda t}\phi '(a) +\lambda e^{\lambda t}\phi (a) = -\mu (a)e^{\lambda t}\phi (a),}$$
where the prime denotes a derivative with respect to a. After canceling $$e^{\lambda t}$$ , we obtain the differential equation for the function ϕ(a):


$$\displaystyle{\phi '(a) +\lambda \phi (a) = -\mu (a)\phi (a).}$$
This is a first-order linear ordinary differential equation, which can be explicitly solved. The solution is given by


$$\displaystyle{\phi (a) =\phi _{0}e^{-\lambda a}\pi (a),}$$
where the constant ϕ 0 = ϕ(0) is unknown and will be specified later. Hence, we have the following form for the solution u(a, t):


$$\displaystyle{u(a,t) =\phi _{0}e^{\lambda t}e^{-\lambda a}\pi (a).}$$
We seek to satisfy the boundary condition in model (12.6):


$$\displaystyle{\phi _{0}e^{\lambda t} =\phi _{ 0}e^{\lambda t}\int _{ 0}^{a_{\dag } }\beta (a)e^{-\lambda a}\pi (a)da.}$$
Therefore, $$\lambda$$ should be identified in such a way that the equation below is satisfied.

The characteristic equation of the McKendrick–von Foerster model is


$$\displaystyle{ \int _{0}^{a_{\dag } }\beta (a)e^{-\lambda a}\pi (a)da = 1. }$$

(12.15)

This is a transcendental equation, and it can have multiple solutions for $$\lambda$$ that are real and complex. It turns out that the Eq. (12.15) has a unique real solution if β(a) is positive on some positive interval:


Lemma 12.1.

Assume $$\beta (a) \geq \hat{\beta }> 0$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_12_Chapter_IEq20.gif”></SPAN> <SPAN class=EmphasisTypeItalic>for a ∈ [a</SPAN> <SUB>1</SUB> <SPAN class=EmphasisTypeItalic>,a</SPAN> <SUB>2</SUB> <SPAN class=EmphasisTypeItalic>]. Then there is a unique real solution</SPAN> <SPAN id=IEq21 class=InlineEquation><IMG alt= such that Eq.  (12.15) holds.

The unique real solution $$\lambda ^{{\ast}}$$ of Eq. (12.15) gives the growth rate of the population for the age-structured model (12.6). It is similar to the population growth rate derived from the Malthusian model.


Definition 12.1.

The parameter $$\lambda ^{{\ast}}$$ is called a Malthusian parameter or growth rate of the population.


Proof.

Define


$$\displaystyle{f(\lambda ) =\int _{ 0}^{a_{\dag } }\beta (a)e^{-\lambda a}\pi (a)da.}$$
Because $$\beta (a)\geqslant \hat{\beta }> 0$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_12_Chapter_IEq24.gif”></SPAN> on a nonzero interval, it follows that <SPAN id=IEq25 class=InlineEquation><IMG alt= is a strictly decreasing function of $$\lambda$$ . The function $$f(\lambda )$$ may not be defined for all values of $$\lambda$$ , particularly if $$a_{\dag } = \infty$$ . Assume that $$f(\lambda )$$ is defined and continuous on the interval $$(-L,\infty )$$ . We have


$$\displaystyle{\lim _{\lambda \rightarrow -L}f(\lambda ) = \infty \qquad \qquad \lim _{\lambda \rightarrow \infty }f(\lambda ) = 0.}$$
Hence, there is a unique real value $$\lambda ^{{\ast}}$$ that satisfies $$f(\lambda ^{{\ast}}) = 1$$ . □ 

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

Related posts:

  1. Ecological Context of Epidemiology
  2. Class-Age Structured Epidemic Models
  3. Techniques for Computing ℛ 0 $$\mathcal{R}_{0}$$
  4. The SIR Model with Demography: General Properties of Planar Systems

Stay updated, free articles. Join our Telegram channel
Tags:
Nov 20, 2016 | Posted by in PUBLIC HEALTH AND EPIDEMIOLOGY | Comments Off on Age-Structured Epidemic Models

Full access? Get Clinical Tree
Get Clinical Tree app for offline access