stats report

attachment_3.pdf

Home >Mathematics homework help >Statistics homework help >stats report

Epidemic Models

SIRS Model

In this set of notes, we will talk about SIR models, which are the standard for of model for infectious diseases. We will follow the steps outlined in the Introduction to Modeling notes:

1. Understand the question.

2. Translate the biological question into a mathematical question via assumptions.

3. Formulate the model that will answer the question.

4. Run/analyze the model.

5. Use the model results to answer the question.

6. Evaluate the reliability of the results.

What do we want to know about epidemics?

We will focus on three questions:

1. What characteristics of a disease make it spread widely and become an epidemic? What is different between epidemic diseases and non-epidemic ones?

2. Why do epidemics often exhibit cycles in number of infected individuals?

3. When do epidemics end?

Assumptions of the SIR model

In an SIR model, we assume that everyone in the population is in one of three categories: suspectible to the disease, infected by the disease, or recovered from the disease. Susceptibles are people who have not had the disease and who can be infected by it. Infecteds are people who currently have the disease. Depending on the disease, there may be a risk of mortality while infected. Recovered means that the person had the disease and have now recovered from it. The term SIR refers to Susceptible-Infected-Recovered. For many diseases, a person has immunity once they have recovered and thus cannot get the disease again. In this case, a person who has recovered will not ever return to the susceptible status. However, with some diseases the immunity can disappear with time and thus people can return to the susceptible status. In this case, the model is called SIRS (Susceptible-Infected-Recovered-Susceptible)

Image from institutefordiseasemodeling.github.io

In order to formulate a mathematical model, we need be more precise about how we think the disease will work. Take S as the number of susceptibles, I as the number of infecteds, R as the number of recovereds, and N as the total population size (i.e. N = S + I + R). We will start with the following assumptions:

1. All population members can be placed into one of three categories: susceptible to the disease, infected by the disease, or recovered from the disease.

2. Susceptible individuals become infected by contact with infected ones. The rate of infection is propor- tional to the rate of contact between susceptibles and infecteds. That is βSI, where β is a constant that determines how infectious the disease is.

3. Infected individuals recover at rate γ. Individuals in the recovered status cannot become infected.

4. Recovered individuals become susceptible again at rate ξ.

5. For now, we will assume no vital dynamics. That is, no birth or death. This means that the disease does not cause mortality and the time scale is short enough that changes in population size are not significant.

Formulate the Model

The basic SIR model is known as a compartment model. A compartment model describes how materials of some kind are transmitted between different compartments or states. In an SIR model, the materials are people (or other organisms) in the population. The compartments are the disease states of susceptible, infected, and recovered. The model describes how people move through these states over time.

We will use system of differential equations to formulate the model:

dt = −βSI + ξR

dt = βSI −γI

dt = γI − ξR

dS dt , dI dt dR dt

are the rate of change in the number of susceptible, infected, and recovered individuals, respectively. The right hand side of each equation quantifies how the value changes depending on the other values. That is,

(Change in Number of Susceptibles)=(loss of susceptibles to infection)+(gain from recovereds that turn susceptible)

(Change in Number of infecteds)=(gain from susceptiblels due to infection)+(loss from infecteds that recover)

(Change in Number of recovereds)=(gain of recovereds from infecteds)+(loss from recovereds that turn susceptible)

Question: Consider the total population size N = S + I + R. How does it change with time? How can we tell?

Analyze the Model

We will have to take a long detour to learn how to analyze the model. This is a set of ordinary differential equations (ODE).These are nonlinear ODEs because the right hand side is not linear in the variables.

Ideally, we would like to solve these equations and get expressions for S(t), I(t), and R(t).

Unfortunately, systems of non-linear ODEs do not generally have analytic solutions (that is, solutions that we can write down in explicit mathematical form). Some special cases can be solved, but we usually have to solve them numerically - that is, use a computer to approximate solutions. We will first talk about equilibrium solutions to the model and then talk about using an R package to solve them numerically.

Equilibrium behavior of the SIRS Model

A common strategy for gaining insight into ODE models is finding the equilibrium behavior. Typically, an ODE model will have short term transient behavior that depends on the initial conditions and will eventually settle into long-term steady state behavior.

In the steady state, the numbers of infecteds, recovereds, and susceptibles are not changing. Thus, we can find this steady state behavior by setting the derivatives equal to zero and solving for S, I, and R:

dt = 0 ⇒−βSI + ξR = 0 ⇒

R = ξ

dt = 0 ⇒ (βS −γ)I = 0

dt = γI − ξR ⇒ I = R

Note also that

S + I + R = N

This gives equations in the three unknowns (S,I,R) that we can solve.

The second equation is zero either if I = 0 or S = γ β In the I=0 case we can see from the third equation that

R=0 also. Then, we have

S̄ = N − I −R = N

In this case, the disease has died out. There are no infected individuals and no recovered ones. The entire population is susceptible.

In the other solution, we have S̄ =

β .

The third equation gives

R̄ = Ī ξ

Then, we have S̄ + Ī + R̄ = N

or γ

β + Ī + Ī

γ = N

Solving for Ī, we get

Ī = γ N − γ

γ + ξ and

R̄ = Ī ξ

γ = ξ

N − γ β

γ + ξ

To summarize, we have

S̄ = γ

Ī = γ N − γ

γ + ξ

R̄ = ξ N − γ

γ + ξ

In this solution, the epidemic goes to an equilibrium with a constant number of susceptibles, infecteds, and recovereds in the population. This does not mean that the individuals are constant in one of the states. Rather, there are always people moving between the three states (S ⇒ I ⇒ R), but the numbers in each state are in equilibrium with the rate of flow of individuals between the states.

To summarize, we have seen that this model has two steady states:

(S̄1, Ī1) = (N, 0)

(S̄2, Ī2) = ( γ

β , ξ ( N − γ

)( γ + ξ)

) ) One steady state in which the disease dies out and one in which it becomes endemic in the population. In order for the second steady state to exist (i.e. be non-negative), we must have

N > γ

If this condition is not met and N < γ β , then only the steady state (N, 0) exists. In this case, the disease

will die out of the population. More formal methods (beyond the scope of this class) demonstrate that this condition is key to the dynamics of the model.

If the condition is not met, typical dynamics will look like:

0 20 40 60 80 100

1 2

0 1

3 0

1 4

0 1

5 0

time

0 20 40 60 80 100

0 1

0 2

0 3

0 4

0 5

time

If the condition N > γ β is met, then disease will persist. This condition is commonly expressed as

R0 = Nβ

γ > 1

R0 is called the intrinsic reproductive rate of the disease. This is the average number of individuals that will be infected by an infected individual. β is the infection rate and γ is the average time to recovery from the disease. Thus, β

γ is the fraction of the population that is infected by an infected individual during the period

of infection. If this greater than one, then the disease will persist.

In this case, both steady states exist, but only the upper one is stable. That the (N, 0) steady state exists means that if there is no infected individuals in the population then the system will remain in that state of having no infecteds. That this steady state is unstable means that if any infecteds are introduced into the population (that is, we move a short distance off of the steady state), then the disease will spread until the equilibrium is reached at the other steady state. Typical dynamics look like

0 50 100 150 200

1 6

0 2

0 0

2 4

0 2

8 0

time

0 50 100 150 200 2

5 3

0 3

5 4

0 4

5 5

time

Answering our first question of interest

This gives us the key to answering our first question of interest:

What characteristics of a disease make it spread widely and become an epidemic? What is different between epidemic diseases and non-epidemic ones?

R0 is the key quantity in determining how widely a disease will spread in the population in its early stages.

R0 = Nβ

Here are estimated values of R0 for several pandemics:

image from https://www.the-scientist.com/features/why-r0-is-problematic-for-predicting- covid-19-spread-67690

The determinants of R0 are

β: New individuals become infected at rate βSI. S and I are the number of susceptible and infected individuals. β is impacted by the behavior of the hosts and the properties of the disease: Behavior: The more contact between susceptible and infected individuals, the higher the rate of transmission. Disease infectiousness: The mechanics of disease transmission determine how readily the disease transmits in a given situation.

γ: The longer the infectious period, the more disease transmission will occur.

Evaluate the reliability of the results.

Class discussion: How reliable do you think that our model is? Can you identity any problems with it?

Submit your answer to eLC.

https://www.the-scientist.com/features/why-r0-is-problematic-for-predicting-covid-19-spread-67690

Numerically Solving Systems of Ordinary Differential Equations:

Consider our SIR Model:

dt = −βSI + ξR

dt = βSI −γI

dt = γI − ξR

Ideally, we would like to solve this find functions S(t),I(t),R(t) that would tell us the number of Susceptibles, Infecteds, and Recovereds at any time t. Unfortunately, it is not usually possible to solve non-linear systems of ODEs. We can solve linear system, but usually not non-linear ones. In a linear system, the right hand side would involve only linear functions of the variables (that is, all terms would either be a constant or a constant multiplied by a single variable)

There are two common approaches to nonlinear ODES: numerical solutions andqualitative analysis. A numerical solution means using a computer program to give approximate solutions over specified times ranges. Qualitative analysis means using mathematical techniques to find out the big picture of how the solutions behave, without actually explicitly solving the ODEs.

We are not going to talk about the mathematics of how numerical ODE solvers work. Rather, we are simply going to apply an ODE solve in R and look at the results. We will use an R package deSolve.

We will need to install the deSolve package in order to use it. You can do this by going to the Tools menu at the top of the R studio window. Click Tools>Install Packages..., then enter deSolve under Packages and click Install. This will install the package on your computer. You then must use the library command to load it.

Within the deSolve package, there is a function called odethat will use to solve the ODEs. Instructions for using the package are at

https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf

Read this to see how to use the package. library(deSolve)

parameters<-c(beta=0.002,xi=0.2,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population

SIRmodel<-function(t,state,parameters) {

with(as.list(c(state,parameters)),{

dS=-beta*S*I+xi*R dI=beta*S*I-gamma*I dR=gamma*I-xi*R

return(list(c(dS,dI,dR)))

}) }

https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf

times<-seq(from=0,to=50,by=0.01)

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout) mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 10 20 30 40 50

2 0

0 1

0 0

time

0 10 20 30 40 50

0 4

0 0

time

0 10 20 30 40 50

0 4

0 0

time

SIR model: beta= 0.002 , xi= 0.2 , gamma= 0.4 , R0= 10

The population starts with 990 susceptibles and 10 infecteds. It goes through some short term dynamics and then settles into a steady state with 200 susceptibles, 267 infecteds, and 533 recovereds.

Recall the condition for persistance of the disease:

R0 = Nβ

γ > 1

In this example, we have (R0=1000*0.002/0.2)

## [1] 10

This is greater than 1, so the disease persists in the population. Let’s do another example where this isn’t true. We will lower the rate of infection to β = 0.0001. Now, we have (R0=0.0001/0.2)

## [1] 5e-04

Now, the intrinsic reproductive rate of the disease is less than 1. That is, each infected individual will on average infect less than one other person. We expect the disease to die out: parameters<-c(beta=0.0001,xi=0.2,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout) mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 10 20 30 40 50

9 9

0 1

0 0

time

0 10 20 30 40 50

0 6

time

0 10 20 30 40 50

0 3

time

SIR model: beta= 1e−04 , xi= 0.2 , gamma= 0.4 , R0= 0.5

As expected, the disease dies out and the population is entirely composed of susceptibles. Let’s look at a few more examples: parameters<-c(beta=0.0004,xi=0.2,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population times<-seq(from=0,to=5000,by=0.01)

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout) mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 1000 3000 5000

9 7

5 1

0 0

time

0 1000 3000 5000

0 6

time

0 1000 3000 5000

0 1

time

SIR model: beta= 4e−04 , xi= 0.2 , gamma= 0.4 , R0= 2

In this example, R0 is only a small distance above 1. The disease does persist, but at a very low level. In fact, we can’t even tell that is is non-zero in the figure. To be sure, let’s look at the actual numerical output tail(SIRout) #this shows the last few lines of SIRout

## time S I R ## [499996,] 4999.95 999.5089 0.1636002 0.3275223 ## [499997,] 4999.96 999.5089 0.1635998 0.3275217 ## [499998,] 4999.97 999.5089 0.1635995 0.3275211 ## [499999,] 4999.98 999.5089 0.1635992 0.3275204 ## [500000,] 4999.99 999.5089 0.1635989 0.3275198 ## [500001,] 5000.00 999.5089 0.1635986 0.3275191

We see that the equlibirum value of I is about 0.16. Of course, we can’t really have fractional individuals. Thus, this is shortcoming of the model. In reality, the disease would die out at this low level of infectivity.

Note that this plot has a much longer time scale than the previous ones.

Next, we will increase β so that R0 = 4: parameters<-c(beta=0.0008,xi=0.2,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population times<-seq(from=0,to=50,by=0.01)

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout)

mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 10 20 30 40 50

5 0

0 1

0 0

time

0 10 20 30 40 50

5 0

2 0

time

0 10 20 30 40 50

0 2

5 0

time

SIR model: beta= 8e−04 , xi= 0.2 , gamma= 0.4 , R0= 4

Now, the steady state has a much higher proportion of infecteds and recovereds.

Increase β to 50: parameters<-c(beta=0.01,xi=0.2,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout) mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 10 20 30 40 50

0 8

0 0

time

0 10 20 30 40 50

0 6

0 0

time

0 10 20 30 40 50

0 4

0 0

time

SIR model: beta= 0.01 , xi= 0.2 , gamma= 0.4 , R0= 50

Now, the steady state situation is that most of the population is either infected or recovered. Very few are susceptible.

Finally, let’s decrease ξ to 0.01. This is decreasing the rate at which recovereds become susceptible so that it rarely happens. Now, most of the population is recovered at any given time. This situation is typical of many viruses that we commonly encounter. That is, most of the population has previously been exposed and is no longer susceptible. parameters<-c(beta=0.01,xi=0.01,gamma=0.4) #specify model parameters state<-c(S=990,I=10,R=0) #initial state of the population

SIRout<-ode(y=state,times=times,func=SIRmodel,parms=parameters)

par(oma=c(0,0,3,0)) plot(SIRout) mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

0 10 20 30 40 50

0 8

0 0

time

0 10 20 30 40 50

0 6

0 0

time

0 10 20 30 40 50

0 6

0 0

time

SIR model: beta= 0.01 , xi= 0.01 , gamma= 0.4 , R0= 1000

Phase Plane Diagrams of ODEs

Phase plan diagrams are another useful tool for analyzing systems of ODEs. A phase plan diagram shows the trajectories of the system at different combinations of the model variables.

We will use the R package PhaseR to plot flow fields for systems of ODEs. This package is designed to be compatible with deSolve. We will consider the example system

dt = xy −y

dt = xy −x

Here is the code: library(phaseR)

## Warning: package 'phaseR' was built under R version 4.0.5

## ------------------------------------------------------------------------------- ## phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems ## ------------------------------------------------------------------------------- ## ## v.2.1: For an overview of the package's functionality enter: ?phaseR ## ## For news on the latest updates enter: news(package = "phaseR")

#function to define the system of equations:

ODE1<-function(t,y,parameters) {

x<-y[1] y<-y[2]

dy<-numeric(2) dy[1]=x*y-y dy[2]=x*y-x

return(list(dy)) }

ODEflow<-flowField(ODE1,xlim=c(-2,2),ylim=c(-2,2),NULL,add=FALSE,main="flow field for example system")

ODENullclines<-nullclines(ODE1,xlim=c(-2,2),ylim=c(-2,2),NULL)

ODETrajectory<-trajectory(ODE1,y0=matrix(c(-2,-1,-1,0,2,0,0,1.5,-1,-2,1.25,1.25,-1,-1.5,1,1.2,0.5,1,1.25,1,1,0.75,0.75,0.5,0.7,0.7),nrow=13,byrow=T),tlim=c(0,2))

## Note: col has been reset as required

−2 −1 0 1 2

− 2

− 1

0 1

flow field for example system

x nullclines y nullclines

The function flowfield creates the plot with the arrows. The arrows show the direction of flow of the system at each (X,y) point. The parameters xlim and ylim define the area of the plot.

The function nullcline adds in the nullclines.These are explained below.

The function trajectory adds curves that show how the system changes with time starting from a specified point. I have the argument y0 gives the collection of starting points for the trajectories (e.g. the first point is (-2,-1) and the second point is (-1,0)). The argument tlim specifies how much time the trajectory should be plotted for. For example, the black curve starting at (-2,-1) shows the trajectory starting from (x,y)=(-2,-1) and going for 2 time units.

The direction of flow at each point is determined by the values of the derivatives of the model at that point. Take the point (x = −1,y = −1). The derivatives are

dt = xy −y = (−1) ∗ (−1) − (−1) = 2

dt = xy −x = (−1)(−1) − (−1) = 2

Thus, the vector of derivatives is (dx dt , dy dt

) = (2, 2). The vector (2, 2) points in the direction 45 degrees. Note that this is the direction of the arrow at (−1,−1) in the plot. Similarly, at (x = 1.5,y = −1),we have

dt = xy −y = (1.5) ∗ (−1) − (−1) = −1.5 + 1 = −0.5

dt = xy −x = (1.5)(−1) − (1.5) = −3

The vector of derivatives is (dx dt , dy dt

) = (−0.5,−3). We see that the arrow at $(1.5,-1) is pointing more downwards but a little to the left.

Next, we will find the nullclines. These are lines along which one of the derivatives is equal to zero. Consider our ODEs. If we set the derivatives to zero, we get

dt = xy −y = 0 =⇒ xy = y =⇒ x = 1,y = 0

dt = xy −x = 0 =⇒ xy = x =⇒ x = 0,y = 1

That is dx dt

= 0 when x = 1 or y = 0 and dy dt

= 0 when x = 0 or y = 1. The lines corresponding to these values are shown on the plot: dark blue for x and light blue for y. Note in the plot that the arrows are horizontal along the light blue line and vertical along the dark blue line.

Where these lines cross, both derivatives are equal to zero and thus there are steady states. This occurs at (x = 0,y = 0), (1, 1). There are no arrows at these points.

The arrows point away from the point (1, 1). This means that this point is an unstable state. Let’s test this with the ODE solver: parameters<-c() #specify model parameters state<-c(x=1.1,y=1.1) #initial state of the population times<-seq(from=0,to=10,by=0.01)

ODEmod1<-function(t,state,parameters) {

with(as.list(c(state,parameters)),{

dx=x*y-y dy=x*y-x

return(list(c(dx,dy)))

}) }

ODEout<-ode(y=state,times=times,func=ODEmod1,parms=parameters)

## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.94531e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.94531e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.61137e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.61137e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.61137e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.28813e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.28813e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.06701e-16

## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.06701e-16 ## ## DLSODA- Warning..Internal T (=R1) and H (=R2) are ## such that in the machine, T + H = T on the next step ## (H = step size). Solver will continue anyway. ## In above message, R1 = 2.39789, R2 = 1.06701e-16 ## ## DLSODA- Above warning has been issued I1 times. ## It will not be issued again for this problem. ## In above message, I1 = 10 ## ## DLSODA- At current T (=R1), MXSTEP (=I1) steps ## taken on this call before reaching TOUT ## In above message, I1 = 5000 ## ## In above message, R1 = 2.39789 ## par(oma=c(0,0,3,0)) plot(ODEout)

0.0 0.5 1.0 1.5 2.0

0 2

0 6

0 1

0 0

time

0.0 0.5 1.0 1.5 2.0

0 2

0 6

0 1

0 0

time

#mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

We started at (x = 1.1,y = 1, 1) - slightly off of the (x = 1,y = 1) equilibrium. We see that from this point the system goes away from (1, 1) towards positive x and positive y. This is what we expect if the equilibrium is unstable. If we try any other points around (x = 1,y = 1) we will find the same thing - that the system will move away from the equilibrium.

Next, Let’s look at the point (x = 0,y = 0). First, start at (x = 0,y = 0): parameters<-c() #specify model parameters state<-c(x=0.5,y=0.5) #initial state of the population times<-seq(from=0,to=10,by=0.01)

ODEmod1<-function(t,state,parameters) {

with(as.list(c(state,parameters)),{

dx=x*y-y dy=x*y-x

return(list(c(dx,dy)))

}) }

ODEout<-ode(y=state,times=times,func=ODEmod1,parms=parameters)

par(oma=c(0,0,3,0)) plot(ODEout)

0 2 4 6 8 10

0 .0

0 .1

0 .2

0 .3

0 .4

0 .5

time

0 2 4 6 8 10 0

.0 0

.1 0

.2 0

.3 0

.4 0

time

#mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

We see that the system goes to the equilibrium at (x = 0,y = 0). Note that the arrow at (x = 0.5,y = 0.5) points towards (x = 0,y = 0). However, if we start at (x = 0,y = 0.5), the system ends up on the y = 1 nullcline and then heads to x = −∞. parameters<-c() #specify model parameters state<-c(x=0,y=0.5) #initial state of the population times<-seq(from=0,to=10,by=0.01)

ODEmod1<-function(t,state,parameters) {

with(as.list(c(state,parameters)),{

dx=x*y-y dy=x*y-x

return(list(c(dx,dy)))

}) }

ODEout<-ode(y=state,times=times,func=ODEmod1,parms=parameters)

par(oma=c(0,0,3,0)) plot(ODEout)

0 2 4 6 8 10

− 1

0 0

− 6

0 0

0 −

2 0

0 0

time

0 2 4 6 8 10

0 .5

0 .6

0 .7

0 .8

0 .9

1 .0

time

#mtext(outer=TRUE,side=3,paste("SIR model:","beta=",parameters[1],", xi=",parameters[2],", gamma=",parameters[3], ", R0=",sum(state)*parameters[1]/parameters[2]))

Adding birth and death to the model

One key factor that our model does not include is the possibility of birth and death. Ignoring birth and death is reasonable if 1) the disease does not cause mortality (or possibly lowe mortality) and 2) the time scale of the dynamics is short enough that the population size won’t change significantly during the course of the epidemic.

Now, we will introduce birth and death into the model. Let µ and ν represent the birth and death rates, respectively. ν is the “background” mortality that occurs in absence of the disease.

We assume

1. birth and death rates are balanced (µ = ν), so that the population size will remain constant.

2. All newborns are in the susceptible class. This is true for many diseases, but not all. For example, HIV can be passed from the mother to the baby in utero.

3. For now, we will assume that the disease does not cause mortality. Thus, the three classes in the model (S,I,R) have the same death rate.

Then, the model is

dt = µN −βSI + ξR−νS

dt = βSI −γI −νI

dt = γI − ξR−νR

New susceptibles enter the system by birth. Individuals in all three classes leave the system by death.

System Dynamics

Here is the model for the function ode: SIRmodelBD<-function(t,z,parameters) {#SIRS model with vital dynamics

beta=parameters[1] xi=parameters[2] gamma=parameters[3] N=parameters[4] mu=parameters[5] nu=parameters[6]

S=z[1] I=z[2]

with(as.list(c(z,parameters)),{

dS=mu*N-beta*S*I+xi*(N-S-I)-nu*S dI=beta*S*I-gamma*I-nu*I

return(list(c(dS,dI)))

}) }

See the file SteadyStateBDSIR.pdf to see the derivation of the steady states. The condition for the disease to persist in the population is

N > γ + µ β

Let’s take µ = ν = 0.01

β = 0.001

ξ = 0.01

γ = 0.2

Then, the condition for the disease to persist is

N > 0.2 + 0.01

0.001 = 210

library(deSolve)

par1<-c(beta=0.001,xi=0.01,gamma=0.2,N=300,mu=0.01,nu=0.01)

state<-c(S=290,I=10) #initial state of the population times<-seq(from=0,to=500,by=0.01)

SIRout<-ode(y=state,times=times,func=SIRmodelBD,parms=par1)

par(oma=c(0,0,3,0)) plot(SIRout)

0 100 300 500

1 8

0 2

2 0

2 6

time

0 100 300 500

5 1

0 1

5 2

0 I

time

In this case, the disease settles in an endemic state where it persists at a low level in the population. If we take N < 210, the disease dies out: library(deSolve)

par1<-c(beta=0.001,xi=0.01,gamma=0.2,N=209,mu=0.01,nu=0.01)

state<-c(S=200,I=9) #initial state of the population times<-seq(from=0,to=500,by=0.01)

SIRout<-ode(y=state,times=times,func=SIRmodelBD,parms=par1)

par(oma=c(0,0,3,0)) plot(SIRout)

0 100 300 500

1 7

5 1

8 5

1 9

5 2

0 5

time

0 100 300 500

0 2

4 6

time

Disease dynamics: Susceptibles must be regenerated for the disease to persist

Now consider what happens if we set ξ = 0. This means that recovered individuals don’t ever become susceptible again:

0 100 300 500

1 0

0 2

0 0

3 0

0 4

0 0

5 0

time

0 100 300 500

0 2

0 4

0 6

0 8

0 1

2 0

time

no return to suseptible state

The disease is still able to persist in the population. Now, I will also set

µ = ν = 0

(no birth and death):

0 100 300 500

1 0

0 2

0 0

3 0

0 4

0 0

5 0

time

0 100 300 500 0

2 0

6 0

1 0

time

no return to susceptible state and no birth/death

We see that the disease dies out. This reveals an important principle:

in order for the disease to persist, there must be a way to renew susceptibles.

There are two ways that this can happen:

1. Loss of immunity causes a gain in susceptibles from recovereds.

2. Gain of new susceptibles via new births.

Disease Dynamics: Cycles in Infections disease

We are now in a position to answer our second question:

Why do epidemics often exhibit cycles in number of infected individuals?

Consider the case below. In this case, we will consider a longer time scale disease. We will take the time units as years. Then, set

γ = 10

. This means that the average time to recovery is 1/10 of a year.

µ = 1/70

This means that the average life expectancy is 70 years.

β = 0.01

This means that there is one infection per year per 100 susceptible-infectious contacts

N = 5000

The population size is 5000 (e.g. a small town)

ξ = 0

This means that recovered people are always immune and do not become susceptible again.

The basic disease reproductive rate is

R0 = Nβ

γ + µ =

5000 ∗ 0.01 10 + 170

= 4.95

Here are the dynamics:

0 20 40 60 80 100

0 .0

0 .8

time (years)

su sc

e p

tib le

0 20 40 60 80 100

0 .0

0 .2

0 .4

time (years)

in fe

ct e

d s

I have scaled the y-axis to show the proportion of the population. We see that about half of the population is infected in the initial outbreak of the disease. This happens over a period of a few months. At this point there are very few susceptibles and so the number of infecteds crashes to a very low level. After that, the scale of the dynamics is much slower. It takes about 30 years for the susceptible to build up enough to get another disease outbreak. This epidemic is less severe because most of the population is still immune. Each subsequent outbreak is quicker and of smaller amplitude. Eventually, epidemics stop and the disease settles into an endemic state at level 1

R0 .

Stages of Disease Cycles: 1. At the initial disease outbreak, the entire population is susceptible. 2. The disease spreads through the population. Susceptibles become infected and then recover. 3. The susceptible population drops quickly. Soon, there are few people who are susceptible to the disease

and so the number of infecteds drops to a low level as well. Now, the disease is present at only a very low level in the population. 4. Gradually, the number of susceptibles builds up again through either birth or loss of immunity. 5. When the number of susceptibles gets high enough again, then the a new disease outbreak will occur and the cycle begins again.

Nature of the Disease Equilibrium

The equilibrium occurs at the point where the number of new infections each year is balanced by the number of new susceptibles. Again take R0 as the intrinsic reproductive rate of the disease. In the case of the model with birth and death, this is given by

R0 = Nβ

γ + µ This is the average number of individuals infected by an infected individual. However, this is only valid at the beginning of the outbreak when the whole population is susceptible. If the number of susceptibles is S, then the disease reproductive rate is

Re = Sβ

γ + µ

Where Re is called the effective intrinsic reproductive rate. In other words, the number of new infections per infected is reduced by the fraction of the population that is susceptible.

Re = S

Nβ

γ + µ =

Sβ

γ + µ

The disease will be at equilibrium when Re = 1. That is, when each infected individual infects an average of one new individual. That is, equilibrium will occur at

Re = 1 ⇒ Sβ

γ + µ = 1 ⇒ S =

γ + µ β

Which we already derived through other means (See the Maple file).

Overshooting the Equilibrium Leads to Cycles

Cycles will tend to happen when the disease “overshoots” the equilibrium. That is, when then infection spreads so rapidly that the number of susceptibles drops below the point at which Re = 1. Then, Re < 1 and the number of susceptibles will increase. However, the number of infecteds will lag and thus the number susceptibles overshoots the equilibrium again, but in the other direction (too high rather than too low). Now, Re > 1) and thus the number of susceptibles drops. The amplitude of the cycles decreases because it overshoots by less and less each cycle.

0 50 100 150 200

0 .2

0 .8

time (years)

su sc

e p

tib le

0 50 100 150 200

0 .0

0 0

.1 0

0 .2

time (years)

in fe

ct e

d s

Individually Based Models

Now, rather than continuing with an ODE model, we are going to introduce a new type of model. This is called an individually based model (IBM) or an Agent Based Model.This is a computer simulation model in which we track individuals in a population. We can make much more sophisticated and realistic models in an IBM framework.

The program below implements an IBM epidemic model. I will discuss how it works below. I will first use it to study the impact of disease mortality on disease dynamics. I will start with a baseline model that does not have disease mortality: #This program implements an individually based SIR model.

set.seed(1222) S0=1000 #initital number of susceptibles I0=10 #initial number of infecteds R0=0 #initital number of recovereds M0=0 beta=0.001 #infection rate per contact per time step gamma=0.1 #recovery rate per time step m=0.0 #disease mortality rate per time step endtime=1000 b=0.02 #birthrate per time step nu=0.01 #background mortality rate per time step.

pop=c(rep(0,S0),rep(1,I0),rep(2,R0)) #pop[i] gives the infected status of subject i. #0 indicates susceptible, 1 indicates infected, 2 indicates recovered, 3 indicates dead

#Function to make new infections: #"flips coin" for each population individual to see whether they get the disease. #if so, their status is changed from susceptible (0) to infected (1). infect<-function(pop,beta) { nI<-sum(pop==1) #number 1of infecteds

nS<-sum(pop==0) #number of susceptibles pop1=pop

pi=1-(1-beta)^nI #probility of susceptible being infected pop1[pop==0]=ifelse(runif(nS)<pi,1,0) #check if each susceptible becomes infected #change to loop to make easier to understand? return(pop1)

}

#function to make individuals recover #"flips coin" for each infected individual to see whether they recover from the disease. #if so, their status is changed from infected (1) to recovered. recover<-function(pop,gamma) {

nI<-sum(pop==1) #number of infected pop1=pop pop1[pop==1]=ifelse(runif(nI)<gamma,2,1) #check if each infected recovers. This is the "coin flip"

return(pop1)

}

#function to check whether each individual dies from non-disease causes backgroundmortality<-function(pop,nu) {

pop1=pop n=length(pop1) pop1=ifelse(runif(n)<nu,3,pop1) #check if . If so, write 3, if not write 0

return(pop1)

}

#function to check whether each individual dies from disease diseasemortality<-function(pop,m) {

nI<-sum(pop==1) #number of infected pop1=pop pop1[pop==1]=ifelse(runif(nI)<m,3,1) #check if each infected dies

return(pop1) }

birth<-function(pop,b) {

n<-sum(pop!=3) #population size nf=floor(n/2) #number of females = half total pop births=rbinom(1,nf,b) #number of births pop1=c(pop,rep(0,births))

return(pop1) }

#initialize the population: St=S0 It=I0 Rt=R0 Mt=M0

#this is the main part of the program. This loops over time. Each time step is one day. On each day, it checks whether each #person has caught the disease, recovered from the disease, died from non-disease causes, died from disease, or had a baby. for (time in 1:endtime) { pop=infect(pop,beta)

pop=recover(pop,gamma) pop=backgroundmortality(pop,nu) pop=diseasemortality(pop,m) pop=birth(pop,b) St=c(St,sum(pop==0)) It=c(It,sum(pop==1)) Rt=c(Rt,sum(pop==2))

Mt=c(Mt,sum(pop==3)) }

#time=1000 #matplot(1:time,cbind(St,It,Rt,Mt)[1:time,],type="l",xlab="time",ylab="number",lty=1:4,col=1:4) #legend("right",c("Susceptible","Infected","Recovered","Died"),lty=1:4,col=1:4)

time=1000 par(mfrow=c(1,1)) matplot(1:time,cbind(St,It,Rt)[1:time,],type="l",xlab="time",ylab="number",lty=c(1,2,4),col=2:4) legend("right",c("Susceptible","Infected","Recovered"),lty=c(1,2,4),col=2:4)

0 200 400 600 800 1000

0 2

0 0

4 0

0 6

0 0

8 0

0 1

0 0

time

n u

m b

e r Susceptible

Infected Recovered

We see that the disease persists in the population at a low level. The fluctuations up and down are from random fluctuations from day to day in what happens.

Note that I balanced the birth and death rate (birth=2xdeath because only females reproduce) so that the population will remain approximately constant.

Now, I will add in disease mortality. I take a mortality rate of 0.001 per unit time. Since the disease lasts on average ten time units (because the recovery rate is 10%), this means that the disease mortality rate is 10*0.001=0.01. That is, the disease has a 1% mortality rate.

Note that I don’t need to rerun the function definitions (infect, recover, etc) set.seed(1984) S0=1000 #initital number of susceptibles I0=10 #initial number of infecteds R0=0 #initital number of recovereds

M0=0 #initial number dead beta=0.001 #infection rate per contact per time step gamma=0.1 #recovery rate per time step m=0.001 #disease mortality rate per time step endtime=1000 b=0.02 #birthrate per time step nu=0.01 #background mortality rate per time step.

pop=c(rep(0,S0),rep(1,I0),rep(2,R0)) #pop[i] gives the infected status of subject i. #0 indicates susceptible, 1 indicates infected, 2 indicates recovered, 3 indicates dead

#initialize the population: St=S0 It=I0 Rt=R0 Mt=M0

pop=recover(pop,gamma) pop=backgroundmortality(pop,nu) pop=diseasemortality(pop,m) pop=birth(pop,b) St=c(St,sum(pop==0)) It=c(It,sum(pop==1)) Rt=c(Rt,sum(pop==2)) Mt=c(Mt,sum(pop==3))

}

#time=1000 #matplot(1:time,cbind(St,It,Rt,Mt)[1:time,],type="l",xlab="time",ylab="number",lty=1:4,col=1:4) #legend("right",c("Susceptible","Infected","Recovered","Died"),lty=1:4,col=1:4)

time=1000 matplot(1:time,cbind(St,It,Rt)[1:time,],type="l",xlab="time",ylab="number",lty=c(1,2,4),col=2:4, main="disease mortality=0.01") legend("right",c("Susceptible","Infected","Recovered"),lty=c(1,2,4),col=2:4)

0 200 400 600 800 1000

0 2

0 0

4 0

0 6

0 0

8 0

0 1

0 0

0 disease mortality=0.01

time

n u

m b

e r Susceptible

Infected Recovered

This doesn’t change things very much because the disease morality is low. We do see a slight downward trend in the number of susceptibles and therefore the total population size. This is because total deaths now exceed total births. Birth and death was balanced before the disease was introduced. Even though the disease morality is low, it is enough to cause an overall decrease in population size.

Now I will increase mortality:

0 200 400 600 800 1000

0 2

0 0

4 0

0 6

0 0

8 0

0 1

0 0

0 disease mortality=0.1

time

n u

m b

e r

Susceptible Infected Recovered

The disease mortality rate is now 0.1 per time unit. Thus, the overall disease morality rate is close to 100%. Now, the disease dies quickly out of the population because it kills too many people. Because they die so quickly, they don’t have time to infect many other people.

If we make the birth rate higher, we can have the disease persist with high mortality because new susceptibles keep coming in to the population:

0 200 400 600 800 1000

0 2

0 0

4 0

0 6

0 0

8 0

0 1

0 0

0 disease mortality=0.1, birth rate=0.035

time

n u

m b

e r Susceptible

Infected Recovered

Note that the recovered population is growing. This is because births are greater than deaths. The disease is able to persist at 10% mortality because there is a constant flow of new births. In these conditions, the disease settles into a state of long term high mortality. That is, 5-10% of the population is infected at any give time and most infected individuals die. The disease persists because there is a constant inflow of susceptible newborns. This is essentially how many of the terrible historical killers of humans like measles and tuberculosis worked. For centuries, before modern vaccination campaigns, they were ALWAYS there killing large numbers of people (mostly children in the case of measles) year after year.

0 200 400 600 800 1000

0 2

0 0

4 0

0 6

0 0

8 0

0 1

0 0

0 disease mortality=0.1, birth rate=0.035

time

n u

m b

e r Susceptible

Infected Recovered

Coding Details of the Simulation.

In order to understand better how the simulation program works, lets look at the infect function: infect<-function(pop,beta) { nI<-sum(pop==1) #number of infecteds

nS<-sum(pop==0) #number of susceptibles pop1=pop

pi=1-(1-beta)^nI #probility of susceptible being infected pop1[pop==0]=ifelse(runif(nS)<pi,1,0) #check if each susceptible becomes infected #change to loop to make easier to understand? return(pop1)

}

In order to demonstrate the code, we will initialize the population and also define beta. I am making the population small (just 20) so that we can more easily see what is going on. pop=c(rep(0,10),rep(1,10),rep(2,0)) beta=0.01

The first two lines in the function are nI<-sum(pop==1) #number 1of infecteds nS<-sum(pop==0) #number of susceptibles

These count the number of infecteds and susceptibles. The code pop==1 creates a vector of TRUE and FALSE

values: TRUE where the value in pop is equal to one and FALSE where it is not equal to one: pop

## [1] 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 pop==1

## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE ## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

We call a vector of TRUE and FALSE values a logical vector. Applying the sum command counts the number of TRUE values: sum(pop==1)

## [1] 10

Next, we calculate the probability that a susceptible is infected: pi=1-(1-beta)^nI pi

## [1] 0.09561792

1. The probability that a suspectible is infected by contact with one infected individual is β.

2. The probability that they are not infected is (1 −β).

3. The probability that they are not infected by any of the nI infecteds is (1 −β)nI .

4. The probability that they are infected by at least one is 1 − (1 −β)nI

Next, we find out who gets infected. We do this by generating a uniform random number between 0 and 1 for each susceptible and comparing to the infection probability pi. runif(nS) generates nS values from a ‘uniform(0,1) distribution: set.seed(3332) runif(nS)

## [1] 0.33223562 0.01786395 0.65783307 0.97043755 0.38175202 0.59197094 ## [7] 0.39927833 0.68927095 0.44402178 0.50352531

The command set.seed makes it so that the same random number sequene is generated each time. This is for illustrative purposes for these notes. We do not have this in the program.

Taking runif(nS)<pi checks whether each random value is below pi. If so, a TRUE is returned. If not a FALSE is returned: set.seed(3332) runif(nS)<pi

## [1] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

We see that the second value is less than pi=0.096. Thus, TRUE is returned. This is the “coin flip” to determine whether an individual is infected. In this example, the second individual is infected and the rest are not.

We use the ifelse command to modify the pop vector. The command ifelse(condition,A,B) takes a logical vector condition and returns A for every element where Condition is TRUE and B for every element whereCondition is FALSE:

set.seed(3332) ifelse(runif(nS)<pi,1,0)

## [1] 0 1 0 0 0 0 0 0 0 0

Thus, it produces 1 for the individual who was infected and 0 for everyone else.

Next, we need to take pop and change the values for the infected individuals from 0 to 1: pop1=pop pop1[pop==0]=ifelse(runif(nS)<pi,1,0)

We first take pop1=pop. This creates a copy of pop.

Then pop1[pop==0] extracts the 0 elements from pop1: pop1=pop pop==0

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE ## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE pop1[pop==0]

## [1] 0 0 0 0 0 0 0 0 0 0

Then, we set this equal to the ifelse. This writes the results of the ifelse into the locations in pop currently occupied by 0’s.

Now, pop1 has been modified so that newly infected individuals have been changed from 0 to 1. The other functions all work in a similar way. They each follow the steps

1. Find out the number of individuals of a relevant type (e.g. suscepibles and infecteds in function infect).

2. Calculate the probability that an event (infection, recover, birth, death) occcurs.

3. Generate a uniform(0,1) random number for each relevant individual to determine if the event happens. 4. Modify pop for the individuals for whom the event occured.

The main part of the program is for (time in 1:endtime) { pop=infect(pop,beta)

pop=recover(pop,gamma) pop=backgroundmortality(pop,nu) pop=diseasemortality(pop,m) pop=birth(pop,b) St=c(St,sum(pop==0)) It=c(It,sum(pop==1)) Rt=c(Rt,sum(pop==2)) Mt=c(Mt,sum(pop==3))

}

This is a loop over time. That is, for (time in 1:endtime) steps through times and executes the commands inside the brackets for each time step.

In each time step, a function for each possible event is applied to the population: pop=infect(pop,beta) pop=recover(pop,gamma)

pop=backgroundmortality(pop,nu) pop=diseasemortality(pop,m)

Then, the current state of the population is recorded: St=c(St,sum(pop==0)) It=c(It,sum(pop==1)) Rt=c(Rt,sum(pop==2)) Mt=c(Mt,sum(pop==3))

Minimum Population Size for Disease Spread

What can we learn from the simulation model that we can’t learn from the ODE model? Consider the ODE run below:

0 20 40 60 80 100

0 .0

0 .2

0 .4

0 .6

0 .8

1 .0

time (years)

fr a

ct io

n o

f p

o p

u la

tio n

susceptible infected

The disease has damped oscillations to a stable equilibrium where the disease persists.

Now, I will run the simulation with the same parameters. I have set the time step in the simulation to be days, whereas the time unit is years in the ODE model. Unlike the ODE model, the time scale matters in the simulation model. If we only did infection, recovery, etc once per year, then we would miss most of the action in the model (which occurs on a much faster time scale).

0.0 0.5 1.0 1.5 2.0 2.5

0 1

0 0

0 3

0 0

0 5

0 0

time (years)

n u

m b

e r Susceptible

Infected Recovered

We see that the infection dies out. Why is this different from the ODE model? Consider detailed output from the ODE model: SIRout[700:725,]

## time S I ## [1,] 6.99 424.6930 3.884276e-05 ## [2,] 7.00 425.3465 3.808725e-05 ## [3,] 7.01 426.0000 3.734742e-05 ## [4,] 7.02 426.6534 3.662291e-05 ## [5,] 7.03 427.3067 3.591339e-05 ## [6,] 7.04 427.9599 3.521854e-05 ## [7,] 7.05 428.6130 3.453804e-05 ## [8,] 7.06 429.2660 3.387157e-05 ## [9,] 7.07 429.9189 3.321883e-05 ## [10,] 7.08 430.5717 3.257952e-05 ## [11,] 7.09 431.2244 3.195335e-05 ## [12,] 7.10 431.8771 3.134003e-05 ## [13,] 7.11 432.5296 3.073928e-05 ## [14,] 7.12 433.1821 3.015084e-05 ## [15,] 7.13 433.8344 2.957443e-05 ## [16,] 7.14 434.4867 2.900980e-05 ## [17,] 7.15 435.1388 2.845670e-05 ## [18,] 7.16 435.7909 2.791486e-05 ## [19,] 7.17 436.4429 2.738406e-05 ## [20,] 7.18 437.0948 2.686405e-05 ## [21,] 7.19 437.7466 2.635460e-05

## [22,] 7.20 438.3983 2.585547e-05 ## [23,] 7.21 439.0499 2.536649e-05 ## [24,] 7.22 439.7014 2.488737e-05 ## [25,] 7.23 440.3528 2.441794e-05 ## [26,] 7.24 441.0042 2.395799e-05

We see that the number of infecteds drops very low (less than $10ˆ{-4}). However, there cannot really be fractional individuals. In a real population the disease would die out. This is what happens in the simulation because it tracks individuals. This reveals another important principle:

In order for a disease to persist, the population must be sufficiently large to hold a resevoir of infecteds when the disease reaches the low point of its cycle.

If the population is large, then that tiny fraction of the population will be small number of infected individuals and then the disease can persist in small pockets. This is one reason why cities are often necessary for epidemics to occur.

How do we deal with this problem in the simulation? One solution would be to make the population size much bigger. Unfortunately, this is not very feasible because the time to run an individually based simulation model increases with the population size. If I increase the population size to e.g. 100 thousand, then it will take hours to run on my computer. I could do it on the UGA computing cluster, but this more effort than we want to do in this class.

Here is a kludge solution: if(sum(pop==1)<1) pop[1]=1 #if infected drops below 1, maintain at 1.

If we add this line to the code, it ensures that there is always at least one infected individual in the population. Thus, when the susceptible get numerous enough again, we can start a new epidemic cycle. #initialize the population: St=S0 It=I0 Rt=R0 Mt=M0

pop=infect(pop,beta) pop=recover(pop,gamma) pop=backgroundmortality(pop,nu) pop=diseasemortality(pop,m) pop=birth(pop,b)

if(sum(pop==1)<1) pop[1]=1 #if infected drops below 1, maintain at 1.

St=c(St,sum(pop==0)) It=c(It,sum(pop==1)) Rt=c(Rt,sum(pop==2)) Mt=c(Mt,sum(pop==3))

}

#time=1000 #matplot(1:time,cbind(St,It,Rt,Mt)[1:time,],type="l",xlab="time",ylab="number",lty=1:4,col=1:4)

#legend("right",c("Susceptible","Infected","Recovered","Died"),lty=1:4,col=1:4) par(mfrow=c(1,1)) time1=1 #74*365 time2=endtime #75*365 matplot((time1:time2)/365,cbind(St,It,Rt)[time1:time2,],type="l",xlab="time (years)",ylab="number",lty=2:4,col=2:4) legend("right",c("Susceptible","Infected","Recovered"),lty=2:4,col=2:4)

#SIR2=cbind(St,It,Rt)

#save(SIR2,file="workspaces/SIR2.RData")

load("workspaces/SIR2.RData")

endtime=200*365 par(mfrow=c(1,1)) time1=1 #74*365 time2=endtime #75*365 matplot((time1:time2)/365,SIR2[time1:time2,],type="l",xlab="time (years)",ylab="number",lty=c(1,2,4),col=2:4) legend("right",c("Susceptible","Infected","Recovered"),lty=c(1,2,4),col=2:4)

0 50 100 150 200

0 1

0 0

0 3

0 0

0 5

0 0

time (years)

n u

m b

e r Susceptible

Infected Recovered

We now see cycles like we expect. This requirement for a minimum population is not something that we can learn from the ODE model.

We will use this model more in our next section of notes.

SIRS Model
What do we want to know about epidemics?
Assumptions of the SIR model
Formulate the Model
Analyze the Model

Equilibrium behavior of the SIRS Model

Answering our first question of interest
Evaluate the reliability of the results.

Submit your answer to eLC.

Numerically Solving Systems of Ordinary Differential Equations:
Phase Plane Diagrams of ODEs
Adding birth and death to the model
System Dynamics

Disease dynamics: Susceptibles must be regenerated for the disease to persist

Disease Dynamics: Cycles in Infections disease

Nature of the Disease Equilibrium
Overshooting the Equilibrium Leads to Cycles

Individually Based Models

Coding Details of the Simulation.
Minimum Population Size for Disease Spread