ME 406S-I-R Model of EpidemicsPart 2Vital Dynamics IncludedsysidMathematica 6.0.3, DynPac 11.01, 1ê13ê20091. Introduction‡Description of the ModelIn this notebook, we include births and deaths in the standard S-I-R model for epidemics. The resultingmodel will allow us to look at events of longer duration. For the integration of the nonlinear differential equa-tions, we use the package DynPac. Although some familiarity with DynPac is assumed, brief descriptions ofsome DynPac commands are given as they are used. The integrations in this notebook also can be done easilywith the Mathematica function NDSolve.The S-I-R model was introduced by W.O. Kermack and A.G. McKendrick ("A Contribution to theMathematical Theory of Epidemics," Proc. Roy. Soc. London A ,115, 700-721, 1927), and has played a majorrole in mathematical epidemiology. In later work, Kermack and McKendrick added the effects of vital dynamics(Proc. Roy. Soc. London A, 138 55-83, 1932, and 141 94-122, 1933). A summary of the model and its uses isgiven by Herbert W. Hethcote (Applied Mathematical Ecology, ed. S. A. Levin, T.G. Hallam, and L.J. Gross,Springer-Verlag, 1989, pp. 119-144). In the model, a population is divided into three groups: the susceptiblesS, the infectives I, and the recovered R, with numbers s, i, and r respectively. The total population is (1) n = s + i + r . The susceptibles are those who are not infected and not immune, the infectives are those who are infected andcan transmit the disease, and the recovered are those who have been infected, have recovered and are perma-nently immune. We will include in the model the natural birth and death rates, although with simplifyingassumptions. We assume that all births are into the susceptibles. We assume that the death rate is equal formembers of all three classes, and we assume that the birth and death rates are equal so that the total population isstationary. Finally, we assume that this is a non-lethal disease so that the recovered are truly recovered and notdead. ‡Differential Equations and Parameters of the ModelThe basic differential equations are as follows: (2)dsdt= -asi + m Hs + i + rL - ms ,(3)didt= asi - bi - mi ,(4)drdt= bi - mr.These equations describe the transitions of individuals from S to I to R. The new parameter is the birth anddeath rate m. By adding the three equations, we show easily that the total population n is constant. As before,the parameter a is the transmissivity and the parameter b the recovery rate. See Part 1 of this notebook andespecially Part 2 of the SEIR notebook for a more detailed discussion of the transmissivity a. In equation (2), we may replace the sum s + i + r by the constant total population n. Then equations (2)and (3) constitute a set of two equations for s and i. After they are solved, we may calculate r from equation (1). 2. Defining the Equations for DynPacIn much of what follows, we will be integrating equations (2) and (3) numerically. In this section, wedefine the equations for DynPac. We start by resetting integration and plot options. deintreset;plotreset; imsize = 250;Now we tell DynPac that the state variables are i and s.setstate@8s, i<D;We tell DynPac that the"official" system parameters are the transmissivity a, the recovery rate b, the birth anddeath rate m, and the total population n.setparm@8a, b, m, n<D;Finally, we define the slope vector of the equations for DynPac.slopevec = 8-a * s * i + m * n - m * s, a * s * i - Hb + mL* i<;We will define integration parameters later. Our first step will be to look at the equilibrium states andtheir stability. 3. Equilibrium and StabilityWe begin by looking for equilibrium states. The DynPac command findpolyeq finds equilibrium statesfor systems with polynomial slope vectors. eq = findpolyeq:8n, 0<, :b + ma, -m H-n a + b + mLa Hb + mL>>The first state is one with no infectives, which we call eq1, and the second is an endemic state which we call eq2.2 sir2.nbeq1 = eq@@1DD8n, 0<eq2 = eq@@2DD:b + ma, -m H-n a + b + mLa Hb + mL>It is worth noting that although we can have a disease-free state with s and r both positive and i zero, it isnot strictly speaking an equilibrium, because the recovered gradually die out and are not replaced. Only when allof the recovered are gone can the population distribution be stationary.Let's look at the stability of the no-infectives state. We first evaluate the derivative matrix at that state.dermat1 = dermat ê. Thread@statevec -> eq1D88-m, -n a<, 80, n a - b - m<<Eigensystem@dermat1D:8n a - b - m, -m<, ::-n an a - b, 1>, 81, 0<>>The second eigenvalue -m is always negative, hence the possibility of an epidemic hinges on the first eigenvalue.We see that the equilibrium is unstable and there will be an epidemic if(5)na - b - m > 0 . This is essentially our old condition that the initial reproductive ratio must exceed 1, only the recovery rate in theearlier formula is now replaced by the b + m, the sum of the recovery rate and the natural death rate. If (5) issatisfied and we introduce an infective, there will be an epidemic. In this case, the equilibrium is a saddle point,and the linear stable manifold, with eigenvalue m, is the s axis -- that is, any perturbation in population with noinfectives will go back to the equilibrium population with an e-folding time of 1/m. The linear unstable manifoldpoints in the direction of increasing i and decreasing s. (See the phase plane plot for an example in section 5below.) The endemic state is relevant only if the equilibrium values of i and s are positive. We see from theexpression for eq2 above that this happens if and only if the condition (5) is satisfied. Therefore we get theendemic state only when the zero infective state is unstable. We would expect the endemic state to be stable inthat case. Let's see if we can show that. We calculate the derivative matrix at eq2, and assign it to dermat2. dermat2 = dermat ê. Thread@statevec -> eq2D::-m +m H-n a + b + mLb + m, -b - m>, :-m H-n a + b + mLb + m, 0>>Now we introduce some simplifying notation:(6)q = b + m , p = (na - (b + m ))/(b + m) .