Copyright is owned by the Author of the thesis. Permission is given for a copy to be downloaded by an individual for the purpose of research and private study only. The thesis may not be reproduced elsewhere without the permission of the Author. Modelling Repeated Epidemics with General Infection Kernels. This thesis is presented in partial fulfilment of the requirement for the Degree of Masters of Information Science in Mathematics At Massey University, Albany, New Zealand Joanne L. Mann 2003 Abstract This thesis is on mathematical modelling in epidemiology, exploring the generic characteristics of diseases in two different population structures. Integral equations are used, to model the epidemics in each generation ( of the epidemic). Difference equations are then used to model the change in the populations between epidemics. Initially, single dimension populations are modelled, where the entire population is considered to be one class. Then the population is split into two classes and a similar analysis is performed, with critical differences noted between the two structures. An analytical approach is taken, with numerical examples. The work in this thesis is not specific to one disease, the main focus is to develop a stepped process between generations of the epidemic and analyse the behaviour. ii Acknowledgements I would like to thank my supervisor, Associate Professor Mick Roberts, for his patience and persistence with me throughout the year, and without whom this topic would have been lost on me! I would also like to thank my family and friends, for without their support full time study would not have been possible. lll Contents Abstract. ................ .......... ............... .... ................. .............. ..... ..................... i Acknowledgements ......... .... ... ... .. ................ ... .. .. .... ..... ...... ..... ... .. .. ... ... ....... ii Contents ............ .... .... ............. .................... .... ... ....... ........ .. .... ...... .. .. .. .... ... iii List of Figures ... ..... .................................................... ........... .... ........... ... .. vi List of Tables ........ .. .. .. .. .... ..... ... ..... ..... .. ... ......... ....................................... vii List of Tables ........................................................................................... vii Chapter 1 Introduction ..................................................... ....... ... .. .. .. ...... ... .. 1 l . l The Model.. ....................................................................................... .................. I 1.2 Probability of Infection ................... .................................................................... 4 1.2. 1 Distribution 1 ..................................................................................... ... ....... 4 1.2.2 Distribution 2 ............................................................................................... 5 1.2.3 Distribution 3 ........................................... ..... ................ ............................... 5 1.2.4 Distribution 4 ................... ........ ............... .. .. ................................. ........ ....... . 5 1.2.5 Distribution 5 .......................................... .... ................................................. 6 1.2.6 Distribution 6 ......... ......... ............... .......... .. ....... ........... .......... ......... ........ .. ... 6 1.3 The Basic Reproduction Ratio ...... ....................................................................... 7 1.4 The Incidence of Infection .................. ... .... ... ................................ ....................... 8 1.5 The Initial Growth Rate ....................................................................................... 9 1.6 Overview ........ ..................... ... .... .. .. .... .. ... ........ .... ... .................... . ..................... 10 Chapter 2 The One Dimensional Model: Epidemic Properties ...... ............ 13 2.1 Distribution One .......................... ...................................................................... 19 2. 1. 1 The Basic Reproduction Ratio ........ ........ ........ ....... ....... .......... .... .... ...... ...... 19 2.1.2 The Initial Growth Rate ...................... ... ...... ............... ............. ......... .......... 20 2.1.3 The Final Size Equation - Small Epidemic .. ... .................. .... ...................... 22 2.2 Distribution Two ........ ............. ....................... .... ............................................... 23 2.2.1 The Basic Reproduction Ratio .. ..................... ....... ...................................... 23 2.2.2 The Initial Growth Rate .................. ... .... ........ .......... ..... ........ ...................... 24 2.2.3 The Final Size Equation - Small Epidemics ................................................ 25 2.3 Distribution Three ....... ............. .................... ............ .... ...... ........... .. .... .............. 26 iv 2.3.1 The Basic Reproduction Ratio ........ ..... ... ......... ....... ......... ..... ...... ............. .. . 26 2.3.2 The Initial Growth Rate ...... .. ............... .... ................. ............. .. .. ..... .... .. ... ... 27 2.3.3 The Final Size Equation - Small Epidemics .. ........ .... ..... .... ................. .. ...... 28 2.4 Distribution Four ..... .... ......... ...... .. ... ....... ... .................... ... ..... ................ .. .. : .... ... 28 2.4.1 The Basic Reproduction Ratio ... ....... .............. ......... ... ............................. .. . 29 2.4.2 The Initial Growth Rate ......................... ............... ......... ...................... ... .... 29 2.4.3 The Final Size Equation - Small Epidemics ..... ........... .................... .. ....... .. . 31 2.5 Distribution Five ............. ...... ........ ...... ..... ............... ... ........ ... .... ....... ..... .. ....... ... 31 2.5.1 The Basic Reproduction Ratio ............. ... ............ ........... ........ ... ... ........ ... .... 32 2.5.2 The Initial Growth Rate ...... ... .. .. ................. .. ............. ... ... ....... .. ... ... .. .......... 33 2.5.3 The Final Size Equation - Small Epidemics .... ... .... .. .... .... .... .. .. .... .. ... ......... . 35 2.6 Distribution Six ............... ................... .................. ..... ...... ....... .. .. ... ........ .. ......... . 36 2.6.1 The Basic Reproduction Ratio ..... .. ...... ..... ................ .......... .. .. .. .................. 37 2.6.2 The Initial Growth Rate ...... ..... .... .. .......... .. ..... ........ .. ... ...... .. ............ ... .. ...... 37 2.6.3 The Final Size Equation - Small Epidemics ........... ........... .................... ... ... 40 2.7 Examples ................... ........ .. ............ ..... ... .. ....... ..... .. ... .. .. ... .. .............. .... ....... ..... 41 Chapter 3 Cobwebs - The One Dimensional Model.. .... ........................... 45 3 .1 The Method .... ..... ... ....... .............. ....... ....... ..... ... ....... .... ... ............ ......... ........ ..... 46 3.2 The Stability of the Fixed Point ..... ................ .............. .. .... .. ....... .... ...... .. ... .. ...... 51 Chapter 4 The Two Dimensional Model - Epidemic Properties ... ........... . 61 Preferential Mixing ..... ........ .... ... ...... .... ...... .. ............. ..... ........ ....... ....... .. ........ .. ... 63 Sexual Mixing .... .. .. .. .... ..... .. ...... ... .. .. .... ... ....... ........ ..... ..... ... ..... ...... ... .. ....... ......... 63 Proportionate Mixing ..... .. ................. .. ............ .. .. ..... ....... ........... ........ ...... ... .. .. ... . 63 Reduced Proportionate Mixing ............. .... .. ..... ... ....... ... .. .. ... ... .. .......... ..... .. ......... . 64 4.1 Distribution One ... .. ............................. .. .................. .. ..... ... ................. ............... 65 4.1.1 Special Cases .. ............. .. ....... ................. .............. ....... ................. ............ ... 66 4.1.2 Initial Growth Rate ....................................... .. ..... ................... ............. ....... 68 4.1.3 Examples ............. ........................................ ........................ .................. ..... 69 4.2 Distribution Two ............................................ ................................................... 70 4.2.1 Special Cases .......... .............................................. ...................................... 71 4.2.2 Initial Growth Rate ..... ........ .. ....... .............. ................................................. 72 4.2.3 Examples ...................................................................... .. ............................ 73 4.3 Distribution Three ................. .. ................................................... .... ................... 74 4.3.1 Special Cases .................................... .. ........................................................ 75 4.3.2 Initial Growth Rate ........................................... .................... ...................... 77 4.3.3 Examples .. ............... ......... .... ...................................................................... 78 4.4 Distribution Four ................................. .. .......................................... ... ............... 79 4.4.1 Special Cases ............................................... ............. ....... ........................... 80 4.4.2 Initial Growth Rate ................... .... ........... ..................... .. ......... ....... .. ...... .... 82 4.4.3 Examples ...... .. ... ... ....... .... .. ... .. .... .. .... .. ... ......... ..... ......... .. .......... ...... ....... .. ... 83 4.5 Distribution Five ..... ..... ..... .. .. ...... ................. ......... ................ ......... .. .... ......... .... 84 4.5.1 Special Cases .. ... .. ..... ... ......... ... ..... .. ....... .... ........................... ... .. .. ... ....... .... . 85 4.5.2 Initial Growth Rate ........................ ...... .. .. .. ..... ....... .. ........... ... ...... ... .. .. ..... ... 87 V 4.5.3 Examples .................. .................................................................... ....... ....... 88 4.6 Distribution Six .......................................................................................... ....... 89 4.6. l Special Cases .. .................................................... ........................................ 89 4.6.2 Initial Growth Rate .. .................... .......................... ............... ............ .......... 91 4.6.3 Examples ............. .. ...................... ... ................................. .. ......................... 92 4. 7 The General Case ...... ... ...................... ... ... ...................................................... ... 93 4. 7. l Special Cases ...... ... ................................................................................. .... 94 4.8 Different A( -r) For Each Class ......................................................................... 96 Chapter 5 The Final Size Equation - Two Dimensional Model ...... ......... 99 5.1 Special Cases ...................... ..................... .................................. .............. ....... 102 Preferential Mixing ............. .................................. ... ....... .... .............................. 102 Sexual Mixing ...... .......................................... ................................................... 103 Proportionate and Reduced Proportionate Mixing .. .. ............................... .. ....... . 104 5.2 Two Dimensional Cobwebs for Sexual Mixing ... ... .... .... ....................... .......... 104 Chapter 6 Conclusion ..................................... ........................................ 107 References ...... ... ................ .... .... .. ................................. .......................... 109 VI List of Figures Figure 1.1 SIR Model - the population is divided into three compartments in relation to the infection: susceptibles (S), infectives (/)and removed (R) ................................. 2 Figure 1.2 Contact Rate/Probability Distributions, time on the x axis ............................ 7 Figure 2.1 Sample final size curves, intersection of lines gives S(oo) where S(O)=IOOO, and Ro=3 ................. ..................................... ........................ ................... .. .......... 18 Figure 2.2 Distribution I ..................................................................................... ........ 19 Figure 2.3 Example of calculating r for distribution I given Ro=3 and R0=6. Other parameter values are listed in Examples section ................................................... 21 Figure 2.4 Distribution 2 ................................................................................... .......... 23 Figure 2.5 Example of calculating r for Distribution 2, given Ro=3 and Ro=6. Parameter values listed in Examples section ......................................................................... 24 Figure 2.6 Distribution 3, a=l ..................................................................................... 26 Figure 2.7 Distribution 4, a = 0.1000 and c = 0.7071 .............. .................................... 29 Figure 2.8 Distribution 5, T1=5, S(0)=1000, a=10 and c=0.3 ....................................... 32 Figure 2.9 Example of calculating r for distribution 5 given Ro=3 and Ro=6, and S(0)=1000. Please see table in Examples for other parameter values .................... 34 Figure 2.10 Distribution 6 .................................. ......................................................... 37 Figure 2.11 Example of calculating r for Distribution 6, given Ro=3 and 6. Parameter values are listed in Examples section to follow . ................................................... 38 Figure 3.2 Cobweb plot, Ro=3, 0=0.8, S0=N=l000, .xo=l ............................................ .49 Figure 3.3 Cobweb plot, Ro=3, 0=0.6, So=N=IOOO, .xo=l ............................................ .49 Figure 3.4 Cobweb plot, Ro=6, 0=0.8, So=N=IOOO, .xo=l ............................................. 50 Figure 3.5 Cobweb plot, R0=6, 0=0.6, So=N=IOOO, .xo=l ............................................. 50 Figure 3.6 The top plot shows the value of Ro at fixed points X, given 0. The second plot show the derivative of the nonlinear equation at fixed points X, given 0 using the value of Ro found from the top plot.. .................................................................... 53 Figure 3.7 Plot to depict where the zero point off(X) lies ........................................... 57 vii List of Tables Table 2.1 Initial growth rates for the six distributions when Ro=3 ................................ 42 Table 2.2 Initial growth rates for the six distributions when Ro=6 ................................ 43 Table 4.1. NumericaJJy calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases and a general case for distribution 1. Ro""-3, T1=5.5 and T2=12.5 in each case ................................................................ 70 Table 4.2. Numerically calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases for distribution 2. R0==3 and T1=5.5 in each case ..................... ........................................................................................ 74 Table 4.3. Numerically calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases and one general case for distribution 3. Ro=,-3 Si(t)=7, S2 (t)=lO and a= l in each case .... .. ...... ............................... ........ 79 Table 4.4. Numerically calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases and one general case for distribution 4. Ro""-3.amd a= 1 in each case ............................................................................. 84 Table 4.5. NumericaJJy calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases and one general case for distribution 5. R0::::3 and a= 1 in each case ............................................................................... 88 Table 4.6. Numerically calculated values of the basic reproduction ratio and the initial growth rate for the four special mixing cases and one general case for distribution 6. Ro::::3, T,=4, T2=7, T3=11 and T4=l4 in each case ............................................. 92 1 Chapter 1 Introduction. Since 1760, when Daniel Bernoulli developed a mathematical model of the impact of vaccination against smallpox1 , there has been an increasing demand for models of various infections. Mathematical models help us achieve a better understanding of the dynamics of infections, which may then allow us to efficiently implement control techniques. Infections are continually changing, whether it is change due to drug resistance or just mutation within the infectious agent. Our population dynamics are also changing, which has a large effect on the transmission and probability of infection - and so the mathematical models must also change. This thesis will look at a simple model for an infection, with six generic examples given throughout the analysis. An integral equation technique is used to model the epidemic and then a discrete mapping system is used to model the dynamics of the population between successive epidemics. We first need to describe the assumptions and terminology used when modelling a disease. 1. 1 The Model Consider a population that can be split into three classes (in relation to an epidemic): those who are susceptible to the infection, the infectious people, and those who are removed from the epidemic (through immunity or 1 Dietz & Heesterbeek (2002) death), i.e. no one member of the population may be infected twice. In the following analysis, we assume that the total population is constant. I S(t)I ~ ·I J(t) __ 1 __ ,1 R(t) Figure 1.1 Sm Model - the population is divided into three compartments in relation to the infection: susceptibles (S), infectives (I)and removed (R). 2 The simplest model is depicted in Figure 1.1. The three compartments can represent population density or population size - as the total population size is assumed to be constant it makes no difference. Susceptibles becomes infected at a rate A resulting from contact with infectives. Contact here is very loosely defined, as the amount of contact needed to become infected will be depend on the infection being modelled. lnfectives then become part of the removed compartment at a constant rate y. As the population size is constant, we know that the change in the population will be zero, .i.e. dS di dR -+-+-=0 dt dt dt The differential equations to describe this model are: dS SI -=-/3x- dt N di SI -=/3x--rI dt N dR -=yI dt (1.1) (1.2) Where x is the rate at which susceptibles contact other members of the population and /3 is the probability of a susceptible being infected given contact with an infected member of the population. We have assumed that the population size is constant, that is: s ( t) + I ( t) + R ( t) = N , so it is easy to see that one of the above equations is redundant. The rate ,l is the force of the infection, that is, the rate that susceptibles become infected2 • We have I l=/Jx­ N (1.3) 3 We can solve equations (1.2) to find a relation between the susceptibles and the infectives in the population3 . A differential equations approach has been used for numerous mathematical models, and there is a large amount of information available for the analysis of such systems. However, with constant contact parameters, the time spent in each compartment is exponentially distributed among the members of the population - this does not fit actual results. So we turn to a slightly different way of constructing a model with the use of integral equations. Using integral equations to model an infection is more intuitive than a differential equations approach, and will match the actual data more closely. However, the down side is, there is not a lot of information 2 Anderson & May (1991) 3 See Roberts & Heesterbeek (2000). published relating to the analysis of such systems. This thesis is based on the integral equation approach that will now be defined4 • 1.2 Probability of Infection 4 The probability of a susceptible being infected depends on their contact with an infective and the probability of infection given this contact, which depends on the time since the infective was itself infected. If we let p(r) be the probability of contact and infection, and x( i-) be the contact rate with an infective, where r is the time since the infective was initially infected, we let A(i-)= p(i-)z(i-) (1.4) So the function A ( i-) represents the probability of contact and infection with an infective at infection timer (the time since infection took place). Throughout the following work, six different functions A(i-) will be used to illustrate the model, where -r ~ O. 1.2.1 Distribution 1 { o. A( i-) = ~· i- < r. :z;~i-~T2 i- > r; 4Please refer to Diekmann & Heesterbeek (2000) for further elaboration on the integral equation approach. When r is less than some specified time T1 or greater than a second specified time T2 there is no chance of a susceptible being infected when contacting an infective. When r lies between the two specified times, there is a constant probability of infection when a susceptible comes in contact with an infective. The period between time zero and T1 can be seen as a latency period in the infection. 1.2.2 Distribution 2 A(i-)={a, 0, 5 This is similar to distribution one, but now there is a constant probability of contact and infection with an infective from time zero to time T1• At any other time there is no chance of infection. 1.2.3 Distribution 3 Where a and c are positive constants. For this distribution, the probability of contact and infection decreases in a negative exponential in the time since infection. Note that this is the same as for the differential equations model, as a member of the population will spend an exponential amount of time within a compartment. 1.2.4 Distribution 4 Again, a and c are positive constants. Here, the probability of contact and infection has a similar shape to the gamma distribution (see Figure 1.2 for further clarification). 1.2.5 Distribution 5 A( i-) = ae-c(r-T,)' , 'Z";?: 0 6 As expected, a, c and T1 are positive constants. The probability of contact and infection takes the shape of a shifted normal distribution curve, but we further truncate this, as we are dealing only on a positive time scale. 1.2.6 Distribution 6 _a_(i--T.) 7;:5-r:5T2 T, -T, I ' 2 I A( i-) = a, T2<-r 1 then the epidemic will continue through the population, and the number of 8 infectives will increase while the number of susceptibles will decrease. We can see that Ra will depend on the population size, the contact rates and the probability of infection, hence: (1.5) 1.4 The Incidence of Infection. The incidence of infection i(t) is the number of new cases per unit time. So we see that it will be equal to the change in the susceptible population (as we have ignored changes in the susceptible population due to other causes). At time t, the number of new cases of the infection depends on the contacts between susceptibles and infectives - those who were infected themselves before time t. So we have: (1.6) where the i0t5 ( t) accounts for the initial introduction of the infection into the population6 • We may also rewrite this in terms of the change in the susceptible population: dS(t) I' dS(t-i-) ---=i0t5(t)-S (t) A( i-)----"-----"-di- dt O dt 6 t5 ( t) is the Dirac 's delta function (1.7) 1.5 The Initial Growth Rate The number of infectives can be modelled by an exponential during the initial phases of infection. So we let 9 i ( t) ""ke" (1.8) for some positive constant r, which we call the initial growth rate of the infection. As we said above, the change in infectives is proportional to the incidence of infection. We can then state: (1.9) As we are examining the initial growth of the infection, we do not need to include the initial introduction of the infection into our population; hence we can omit the i08(t) term. We also set the size of the susceptible population equal to its initial value 7, S ( t) = S ( O). So we solve: (1.10) It is shown in Diekmann and Heesterbeek (2000), that there is a unique real r that solves equation (1.10). Note that equation (1.10) is similar to our equation for the basic reproduction ratio ( equation ( 1.5)) The correlation between the two lead to two important facts: r > O if and only if the basic reproduction ratio is greater than one, and r < O if and only if the basic 7 Note: S ( t) » i0 , and so we let S ( o+) = S ( o-) . i0 will usually be assumed to be equal to one, i.e. there will be one initial case to introduce the infection into the susceptible population. reproduction ratio is less than one. That is, we only have initial growth of the infection if we have an epidemic. 1. 6 Overview 10 The purpose of the following exercises is to detennine Ro (the basic reproduction ratio), r (the initial growth rate) and the final size equation of an epidemic, given a function A(-r) that characterises the epidemic. The susceptible population will first be viewed as one class, and will then be split into two classes with intra-class mixing introduced. All the calculations will be based on the following relation for the incidence of the epidemic (1.11) for our six functions A( -r) and for constant and non-constant s (t) . Two methods will be used to calculate the final size of the infection. For a small epidemic, Ro< 1, we assume that S(t) is constant and equal to the initial susceptible population (as there will be no major epidemic, so the change in the population due to the infection is slower that any other change in the population). For a larger epidemic, Ra > 1 , we can not assume that the susceptible population is constant, so we use a direct method applied to equation ( 1. 7) to calculate the final size of the epidemic. Using the methods outlines above, we then construct a repeated epidemic process, where we consider epidemics on a discrete generation basis. Initially we assume that the entire population is susceptible, and we let an epidemic occur. We then calculate the final number of susceptibles and let 11 a portion of them continue on to the next epidemic generation. New susceptibles are introduced into the population to maintain a constant population. We then let another epidemic occur, calculate the final number of susceptibles from this second generation and let a proportion continue and introduce new members into the population. This is repeated, with either an epidemic occurring each generation or the infection not persisting within the population. Numerical calculations are given for this, and then a full analytic proof into the nature of the solution is given. We then repeat out analysis of the six functions A ( -c) when the population is split into two subclasses. The basic reproduction ratio will be calculated for four different mixing schemes between classes. The same methods can be used as for the one dimensional case, with slight alterations to the equations. The final size equations are calculated, and a brief introduction into applying a repeated epidemic process is given. MATLAB has been used to generate the numerical examples with this thesis, and Maple was used for some of the analytical work.