We consider a Leslie-type matrix approach to an SIR epidemic in discrete time. We give examples of the population of susceptibles, infectives, and removals for different birth rates and two different infection rates. Finally, when the infection rate depe

Gani, Joseph

Stals, Linda

English

The Australian National University

A Leslie Matrix Approach to an Age-StructuredEpidemicJoe Gani∗ Linda Stals∗June 28, 2012AbstractWe consider a Leslie-type matrix approach to an SIR epidemic in dis-crete time. We give examples of the population of susceptibles, infectivesand removals for different birthrates and two different infection rates. Fi-nally, when the infection rate depends on the number of infectives, wederive conditions for a steady state.1 IntroductionAn interesting problem in epidemiology is the modelling of age-structured epi-demics. Several such studies in continuous time have been carried out in recentyears, among them [2, 3, 5, 7].For example, the authors of [2] derive a set of partial differential equations forthe time-dependent age-specific densities of susceptible, infective and removedindividuals in an SIR epidemic. Among other results, the authors obtain con-ditions for the existence of steady states.Their calculations are lengthy and result in a complicated system of equa-tions, from which, it is not easy to get an intuitive insight. One is led to considerif a simpler model based on a discrete time Leslie-type matrix [4], might serveequally well as a model. Models of this type have already been used by Zhou, Maand Brauer [8] for the transmission of SARS in China, and by Powell, Slapnicarand van der Werf [6] for the spread of plant pathogens.The purpose of this paper is not to develop new mathematical methods, butrather to indicate that discrete time age-structured epidemic models can char-acterise SIR epidemics very adequately. We propose to examine this alternative,and find conditions under which steady-state solutions exist.∗Mathematical Sciences Institute, Australian National University, Canberra ACT 0200,Australia12 The simplest Leslie-type matrixSuppose that at time t = 0, 1, 2, · · · we have a population consisting of N + 1groups of size s(0, t), s(1, t), · · · , s(N, t) of susceptibles s(a, t) of ages 0 ≤ a ≤ N(from a to a + 1 − 0), with death occurring when a = N + 1. There are i(t)infectives who spread the disease, while birth rates f(a) ≥ 0 apply to each ofthese age groups, with some of the f(a) = 0. If we assume that there is aprobability 0 < q < 1 of non-infection, independent of the number of infectives,an infection rate p = 1− q and a removal rate r > 0 of infectives in an SIR typeepidemic process, then we can writey(t+ 1) = Ay(t), (1)where y(t+ 1) ∈ RN+2, y(t)a = s(a, t) for 0 ≤ a ≤ N and y(t)N+1 = i(t) andA =f(0) f(1) f(2) · · · f(N − 2) f(N − 1) f(N) 0q 0 0 · · · 0 0 0 00 q 0 · · · 0 0 0 00 0 q · · · 0 0 0 0....... . ..........0 0 0 · · · q 0 0 00 0 0 · · · 0 q 0 0p p p · · · p p p 1− r. (2)We see that the total number of births in (t, t+ 1) is the sum∑Na=0 f(a)s(a, t)while the number of infectives remaining after removal isi(t+ 1) = (1− r)i(t) + pN∑a=0s(a, t). (3)We now present some numerical examples of such populations.Example 2.1 Let us suppose that we have a population of insects of ages a =0, 1 and 2, where only the group aged 1 can give birth at the rate of f(1) = 0.85.We may start at t = 0 with s(0, 0) = 412, s(1, 0) = 492, s(2, 0) = 465 withdeath occurring when a = 3. At t = 0, i(0) = 1 so we begin with one infectivein the population, with the probability 0.9 of non-infection and a removal rater = 0.5.Thus, for the matrixA =0 0.85 0 00.9 0 0 00 0.9 0 00.1 0.1 0.1 0.5 , (4)we find as t = 0, 1, 2, 3, 4, 5 the population evolves as shown in Table 1 with thetotal population size decreasing from 1370 to 896. The calculations were done in2Table 1: Change in the population size for Example 2.1t 0 1 2 3 4 5s(0, t) 412 418 315 320 241 245s(1, t) 492 371 376 284 288 217s(2, t) 465 443 334 339 255 259i(t) 1 137 192 198 193 175double precision, but the s(a, t) and i(t) presented in the table are rounded off tothe nearest integer. The continual decrease in the population in the long-termis shown in figures 1 and 2.0 5 10 15 20 25 30 35 40 45 50050100150200250300350400450500timesusceptiblesChange in susceptibles  s(0, t)s(1, t)s(2, t)Figure 1: Change in the number of susceptibles in each age-group over time forExample 2.1.One may well ask if such models can lead to the gradual extinction of theinfectives i(t). From our simple insect model, we can see that i(t+ 1) < i(t) ifp (s(0, t) + s(1, t) + s(2, t)) + (1− r)i(t) < i(t),orprS(t) < i(t),where S(t) =∑Na=0 s(a, t) and p/r = 1/5 in the example given.For Example 2.1, with a birth rate of 0.85, the size of the population rapidlydecreases, as does also the number i(t) of infectives. For t = 15 the populationis y(t) =[65 58 69 48.19]Tor S(t)/5 = 192/5 = 38.4 so that a furtherdecline of the infectives is expected. In fact, the decline in both susceptiblesand infectives is steady, by t = 30, S(30) = 28, i(30) = 7.41.30 5 10 15 20 25 30 35 40 45 500200400600800100012001400timeChange in population  ∑a s(a, t)i (t)Figure 2: Change in the population size over time for Example 2.1.Table 2: Change in the population size for Example 2.2t 0 1 2 3 4 5s(0, t) 412 566 426 586 441 606s(1, t) 492 371 509 384 527 397s(2, t) 465 443 334 458 345 474i(t) 1 137 207 230 258 260Example 2.2 In the next example we increase f(1) to 1.15. All other param-eters are the same as Example 2.1.In which caseA =0 1.15 0 00.9 0 0 00 0.9 0 00.1 0.1 0.1 0.5. . (5)The change in each age-group for t = 0, 1, 2, 3, 4, 5 is given in Table 2. Theapparent growth evident in the table is confirmed over the long-term in Figure3.Note that if at any time t > 0, i(t) > S(t)p/r then i(t + 1) < i(t). FromEquation (3),i(t+ 1)− i(t) = pS(t)− ri(t) = pS(t)− r(S(t)pr+ ε), ε > 0,so that for this t, i(t+ 1)− i(t) = −rε < 0. In Section 3, we state in detail theconditions under which both S(t) and i(t) decrease as t increases.40 5 10 15 20 25 30 35 40 45 500500100015002000250030003500timeChange in population  ∑a s(a, t)i (t)Figure 3: Change in the population size over time for Example 2.2.3 Some consequences of the matrix structureConsider the Leslie-type matrix given in (2). We can write the matrix in blockform asA =[L 0xT 1− r],where L is the standard Leslie matrixL =f(0) f(1) f(2) · · · f(N − 2) f(N − 1) f(N)q 0 0 · · · 0 0 00 q 0 · · · 0 0 00 0 q · · · 0 0 0....... . .......0 0 0 · · · q 0 00 0 0 · · · 0 q 0, (6)and x is a vector with xa = p (0 ≤ a ≤ N). Consequently, by cofactor expansion,|A− λI2| = (1− r − λ) |L− λI1| ,and hence the eigenvalues of A are 1 − r together with the eigenvalues of L.The matrices I2 ∈ R(N+2)×(N+2) and I1 ∈ R(N+1)×(N+1) are identity matrices.The standard Leslie matrix is well known and techniques for calculating itseigenvalues appear in elementary linear algebra text books such as [1].Consider the eigenvalue λN+1 = 1− r. Let uN+1 be a vector such that the(N + 2)th entry is one, and all remaining entries are zero. Then(A− λN+1I2)uN+1 =[L− λN+1I1 0xT 0]uN+1 = 0.5So, uN+1 is an eigenvector.We now focus on λa, (0 ≤ a ≤ N), where λa is an eigenvalue of L withcorresponding eigenvector va. Let ua =[vTa α]Tfor some scalar term α.Then(A− λaI2)[vaα]=[L− λaI1 0xT 0] [vaα]=[(L− λaI)vaxTva + α(1− r − λa)]=[0xTva + α(1− r − λa)].Therefore, if α = −xTva/(1− r − λa), ua is an eigenvector.To be more specific, let us look at the examples in which the Leslie matrixis of the formL =f(0) f(1) 0q 0 00 q 0 .The eigenvalues of L areλ1 =f(0)−√f2(0) + 4f(1)q2, λ2 =f(0) +√f2(0) + 4f(1)q2and λ3 = 0with corresponding eigenvectorsv1 =f(1)/(λ1 − f(0))1q/λ1 , v2 =f(1)/(λ2 − f(0))1q/λ2 and v3 =001 .So, the eigenvalues ofA =f(0) f(1) 0 0q 0 0 00 q 0 0p p p 1− rare λ1, λ2 and λ3, as listed above, with λ4 = 1−r and corresponding eigenvectorsu1 =f(1)λ1−f(0)1qλ1α1 , u2 =f(1)λ2−f(0)1qλ2α2 ,u3 =001α3 and u4 =0001 ,whereα1 =−p [q(λ1 − f(0)) + λ1(λ1 − f(0)) + f(1)λ1]λ1(λ1 − f(0))(1− r − λ1),α2 =−p [q(λ2 − f(0)) + λ2(λ2 − f(0)) + f(1)λ2]λ2(λ2 − f(0))(1− r − λ2),α3 =−p1− r.6Once we obtain the eigenvalues and eigenvectors we can determine the long-term behaviour of the system.Consider examples 2.1 and 2.2 from Section 2. In this case f(0) = 0, sothe eigenvalues are λ1 = −√fq, λ2 =√fq, λ3 = 0 and λ4 = 1 − r, wheref(1) = f . Given that the eigenvectors u1, u2, u3 and u4 form a basis, ourinitial population distribution can be written asy(0) = c1u1 + c2u2 + c3u3 + c4u4,for some scalar terms c1, c2, c3 and c4. Consequently,y(t) = c1λt1u1 + c2λt2u2 + c3λt3u4 + c4λt4u4= c1(−√fq)tu1 + c2(√fq)tu2 + c4(1− r)tu4.As (1− r) < 1, the last term will always approach zero. If√fq < 1, or fq < 1,the overall population will decrease with time, while if√fq > 1 it will increasewith time. If√fq = 1, the population will achieve a steady-state solution thatoscillates between c1u1 + c2u2 and −c1u1 + c2u2. With Example 2.1, where Ais defined in (4), fq = 0.9 × 0.85 = 0.765 < 1 and we see that the populationdecreases with time. However, with Example 2.2, , where A is defined in (5),fq = 0.9× 1.15 = 1.035 > 1 and the population increases with time.As t grows, we expect y(t) to approach(√fq)t(c1(−1)tu1 + c2u2).For the example given in (4), this means we expect y(t) to approach(0.8746)tc1(−1)t−0.97181.0000−1.02900.0728+ c20.97181.00001.02900.8010 .In this case there is no single dominant eigenvalue and the long-term y(t)depends on the initial distribution y(0). However, we can still determine thequalitative behaviour of the different age-groups. For example, the (−1)t termindicates that s(a, t) will oscillate. We would also expect s(0, t), s(1, t) ands(2, t) to reach roughly similar percentages of the total population, as theircorresponding entries in the eigenvector all have about the same absolute value(their absolute values are between 0.97 and 1.03). Finally, we can see that theinfectives should represent a smaller percentage of the overall population. UsingMATLAB to calculate the distribution of each group gives the results in Figure 4,where such qualitative behaviour is evident.Example 3.1 Instead of setting f(0) = 0, let us try f(0) = 0.56. All of theother parameters are the same as Example 2.1.70 5 10 15 20 25 30 35 40 45 500510152025303540time%Change in population percentage  s(0, t)s(1, t)s(2, t)i (t)Figure 4: Change in the population ratio over time for Example 2.1 when A isgiven by (4).So that,A =0.56 0.85 0 00.9 0 0 00 0.9 0 00.1 0.1 0.1 0.5 . (7)Then, we have a single dominant eigenvalueλ =f(0) +√f2(0) + 4f(1)q2= 1.1984,with corresponding eigenvectoru =1.33151.00000.75100.4414 .Consequently, we expect the long-term solution to approach some multiple ofu. Figure 5 shows the distribution of each age-group when A is given by (7).80 5 10 15 20 25 30 35 40 45 50051015202530354045time%Change in population percentage  s(0, t)s(1, t)s(2, t)i (t)Figure 5: Change in the population ratio over time for Example 3.1 when A isgiven by (7).4 A more realistic Leslie-type matrixIf instead of a probability q of non-infection, we replace this by qi(t), for i(t)independent infectives, then the Leslie-type matrix becomesA(t) =f(0) f(1) f(2) · · · f(N − 2) f(N − 1) f(N) 0qi(t) 0 0 · · · 0 0 0 00 qi(t) 0 · · · 0 0 0 00 0 qi(t) · · · 0 0 0 0....... . ..........0 0 0 · · · qi(t) 0 0 00 0 0 · · · 0 qi(t) 0 01− qi(t) 1− qi(t) 1− qi(t) · · · 1− qi(t) 1− qi(t) 1− qi(t) 1− r(8)and, once again, we can write the SIR process as y(t+ 1) = A(t)y(t).Example 4.1 This starts with the same initial populations, s(0, 0) = 412,s(1, 0) = 492 and s(2, 0) = 465 as for Example 2.1, but with the Leslie-typematrixA =0 0.85 0 00.9i(t) 0 0 00 0.9i(t) 0 01− 0.9i(t) 1− 0.9i(t) 1− 0.9i(t) 0.5 .Example 4.2 Similarly we apply the Leslie-type matrix (8) to Example 2.2where f(1) = 1.15.9Example 4.3 The final example applies the Leslie-type matrix (8) to Example3.1 where f(1) = 0.56 and f(1) = 0.85.1 2 3 4 5 6 7 8 9 1002004006008001000120014001600timesusceptiblesChange in susceptibles and infectives  s(0, t)s(1, t)s(2, t)i (t)Figure 6: Change in the number of susceptibles and infectives in each age-groupover time for Example 4.1.The results for the SIR process y(t+ 1) = A(t)y(t) using the examples 4.1,4.2 and 4.3 are shown in figures 6, 7 and 8. We can see that the results are nowmarkedly different from those shown in Sections 2 and 3. In these cases, as inthe more realistic cases, most of the insects end up as infectives.Is it possible to find examples where the steady-state solution does not haveall susceptibles removed from the system? In other words, can we find y(t+ 1)such that y(t+ 1)− y(t) = 0, but y(t+ 1)a = s(a, t+ 1) 6= 0 (for 0 ≤ a ≤ N)?Now, y(t+ 1)− y(t) = 0 if y(t) is in the null space of(A(t)− I2) =[L(t)− I1 0x(t)T −r].where Ii ∈ R(N+i)×(N+i), andL(t) =f(0) f(1) f(2) · · · f(N − 2) f(N − 1) f(N)qi(t) 0 0 · · · 0 0 00 qi(t) 0 · · · 0 0 00 0 qi(t) · · · 0 0 0....... . .......0 0 0 · · · qi(t) 0 00 0 0 · · · 0 qi(t) 0, (9)101 2 3 4 5 6 7 8 9 1002004006008001000120014001600timesusceptiblesChange in susceptibles and infectives  s(0, t)s(1, t)s(2, t)i (t)Figure 7: Change in the number of susceptibles and infectives in each age-groupover time for Example 4.2.1 2 3 4 5 6 7 8 9 1002004006008001000120014001600timesusceptiblesChange in susceptibles and infectives  s(0, t)s(1, t)s(2, t)i (t)Figure 8: Change in the number of susceptibles and infectives in each age-groupover time for Example 4.3.11with x(t)a = 1 − qi(t). Let us focus on L(t) − I1. A vector s ∈ RN+1 is in thenull space of L(t)− I1 if0 = (L− I1)s=f(0)− 1 f(1) f(2) · · · f(N − 2) f(N − 1) f(N)qi(t) −1 0 · · · 0 0 00 qi(t) −1 · · · 0 0 00 0 qi(t) · · · 0 0 0....... . .......0 0 0 · · · qi(t) −1 00 0 0 · · · 0 qi(t) −1s0s1s2...sN−1sN.Using back substitution, we getqi(t)sN−1 − sN = 0⇒ sN−1 = q−i(t)sN ,qi(t)sN−2 − sN−1 = 0⇒ sN−2 = q−i(t)sN−1 = q−2i(t)sN ,and so forth, until finallyqi(t)s0 − s1 = 0⇒ s0 = q−i(t)s1 = q−Ni(t)sN .But we still need to take the first row into account. From the first row,0 = (f(0)− 1)s0 + f(1)s1 + f(2)s2 + · · ·+ f(N)sN={(f(0)− 1)q−Ni(t) + f(1)q−(N−1)i(t) + · · ·+ f(N)}sN .Recall that we are looking for the steady-state solutions, so we are assumingi(t + 1) = i(t) = i for all t > T , at some point T in time. Consequently, if thefollowing constraint0 = (f(0)− 1)q−Ni + f(1)q−(N−1)i + · · ·+ f(N) (10)is satisfied, the non-trivial null space of L(t)− I1 is spanned by,w =[q−Ni q−(N−1)i . . . q−i 1]T. (11)If Constraint (10) is not satisfied, the null space of L(t)− I1 only contains thezero vector.In our example problems we have f(0) = f(2) = f(3) = 0. So, for Constraint(10) to be satisfied, we must have f(1) = q−i.Suppose we have s that belongs to the null space of L(t)− I1, then s = cwfor some scalar term c. Furthermore, y =[sT i]T ∈ RN+2 is in the null space12of A(t)− I2 if xTs− ri = 0. In other words,i =1rxTs=1r(1− qi) N∑a=0sa=cr(1− qi) N∑a=0q−ai=cr(q−Ni − qi).Rearranging the terms in the above equations gives the following additionalconstraint on i0 =rci+ qi(1− q−(N+1)i). (12)To find an example of a non-zero steady-state solution we set q = 0.9 andfound an example of i that satisfied Constraint (12) with r/c = 0.5. That is,we set i = 5.969091. We then set f(1) = q−i and used Equation (11) to finds, with c = 1.5. As shown in Figure 9, the solution does not change with time.This appears to be a stable equilibrium point. We tried changing i to i = 2.4691(set f(1) = q−i and used Equation (11) to find s), where the resulting systemshown in Figure 10 is periodic.0 5 10 15 20 25 30 35 40 45 500123456timeAn exmaple of a steady−state solution (i = 5.969091)  s(0, t)s(1, t)s(2, t)i (t)Figure 9: An example of a non-zero steady-state solution130 5 10 15 20 25 30 35 40 45 500123456timeAn exmaple of a periodic solution (i =  2.4691)  s(0, t)s(1, t)s(2, t)i (t)Figure 10: An example of a periodic equilibrium solution5 ConclusionOur discrete time approach proves simpler than the continuous time analysispresented in [2]. It offers readily calculable results, such as those of Table 2and the graphs in figures 2, 4, 5, and 6. From the view point of the practisingepidemiologist, this simple applicability of the Leslie-type matrix approach is adistinct asset.References[1] H. Anton and C. Rorres. Elementary Linear Algebra. Wiley, 10th edition,2010.[2] Y. Cha, M. Iannelli, and F. A. Milner. Existence and uniqueness of endemicstates for the age structured S-I-R epidemic model. Math. Biosci., 150:177–190, 1998.[3] J. M. Hyman and J. Li. Infection age-structured epidemic models withbehavior change or treatment. J. Biol. Dynamics, 1:109–131, 2007.[4] P. H. Leslie. On the use of matrices in certain population mathematics.Biometrika, 33:183–212, 1945.[5] X.-Z. Li and J.-X. Liu. An age-structured two-strain epidemic model withsuper-infection. Math. Biosci. Eng., 7:123–147, 2010.[6] J. A. Powell, I. Slapnicar, and W. van der Werf. Epidemic spread of alesion-forming plant pathogen–analysis of a mechanistic model with infiniteage structure. Lin. Alg. Applns, 398:117–140, 2005.14[7] E. Shim, Z. Feng, M. Martcheva, and C. Castillo-Chavez. An age-structuredepidemic model of rotavirus with vaccination. J. Math. Biol, 53:719–746,2006.[8] Y. Zhou, Z. Ma, and F. Brauer. A discrete epidemic model for SARS trans-mission and control in China. Math. Comput. Modelling, 40(13):1491–1506,2004.15

A Leslie matrix approach to an age-structured epidemic

https://openresearch-repository.anu.edu.au/bitstream/1885/73112/9/01_Gani_A_Leslie_Matrix_2013.pdf

A Leslie matrix approach to an age-structured epidemic

Abstract

Similar works

Full text

Available Versions

The Australian National University