+1443 776-2705 panelessays@gmail.com
Select Page

Population Dynamics

Next, we will consider mathematical modeling of population dynamics. These models are heavily used in the
science of ecology and help us understand how populations interact with their environment and change over
time.

Take N(t) as the size of a bacterial population at time t. Suppose that the bacteria divide and produce new
bacterial cells at rate r per unit time. Then the population size at time t + ∆t is

N(t + ∆t) = N(t) + rN(t)∆t

Or

N(t + ∆t) −N(t)
∆t

= rN(r)

We take the limit as ∆t → 0 and get

dN

dt
= rN(t)

Next, we want to solve this model. Multiply both sides by dt and divide by N:

dN

N
= rdt

Integrate both sides

∫ t
0

dN

N
=
∫ t

0
rdt

or
ln N(t)|t0 = rt

ln N(t) − ln N(0) = rt

ln N(t) = ln N(0) + rt

N(t) = N(0)ert

This describes the growth of a population over time. N(0) is the population size at the initial time 0, r is the
population growth rate, and N(t) is the population at time t. Consider a population that starts at size 2 and
has a growth rate of 0.69 per year (the amount needed to double once per year):

1

N0=2
r<-c(1)
t=seq(from=0,to=100,by=0.1)
par(mfrow=c(1,1))
matplot(t,N0*exp(t(outer(r,t))),type=”l”,xlab=”time (years)”,ylab=”population size N(t)”,col=1:3,lty=1:3, main=paste(“exponential growth, r=”,r))

0 20 40 60 80 100

0
e

+
0

0
2

e
+

4
3

4
e

+
4

3

exponential growth, r= 1

time (years)

p
o

p
u

la
tio

n
s

iz
e

N
(t

)

#legend(“topleft”,c(paste(“r=”,r[1]),paste(“r=”,r[2]),paste(“r=”,r[3])),lty=1:3,col=1:3)

In 100 years, the population grows to size 1040. If this were a population of rabbits (who can easily reproduce
much faster than this), the mass of 1040 rabbits would be more than a billion times the mass of the sun.

Of course, populations cannot really grow infinitely. Instead, populations are regulated by resource limitations.

Discussion Questions

1. What resources put limits on human population growth? How about trees?

2. Make a plot of how you would expect the population growth rate to vary with population size.

2

Logistic Growth Model

A more realistic model of population growth must include resource limitations. One commonly used model is
the logistic growth model:

dN

dt
= rN(1 −

N

K
)

where

N is population size
r is the maximum population growth rate
K is the carrying capacity of the population. This is the maximum population that can be sustained given
the available resources.

The population growth rate in this model is

growth rate = r(1 −
N

K
)

When the population size N is small compared to the carrying capacity K, then we have N
K
≈ 0 and the

growth rate is approximately r:

dN

dt
≈ rN

In this case, the population grows approximately exponentially. As the population size increases, then the
population growth rate r(1 − N

K
) decreases.

r=1
K=1000
N=1:2000
growthrate=r*(1-N/K)

plot(N,growthrate,type=”l”,xlab=”population size N”, ylab=”population growth rate”, main=”population growth rate under logistic growth model”)

3

0 500 1000 1500 2000

1

.0

0
.5

0
.0

0
.5

1
.0

population growth rate under logistic growth model

population size N

p
o

p
u

la
tio

n
g

ro
w

th
r

a
te

We see that the population growth rate decreases linearly as the size of the population increases.

Discussion Question Make a plot of how you think that population will vary with time under the logistic
growth model starting at a low population size.

4

Solving the Logistic Growth Model

Unlike most population models, the logistic growth model actually has a solution:

N =
K

1 + K−N0
N0

e−rt

The plot below
r=0.69
K=1000
t=seq(from=0,to=20,by=0.1)

Nt=K/(1+(K-N0)/N0*exp(-r*t))

plot(t,Nt,type=”l”,xlab=”time (years)”, ylab=”population size N”, main=”logistic growth model with K=1000″)

0 5 10 15 20

0
2

0
0

4
0

0
6

0
0

8
0

0
1

0
0

0

logistic growth model with K=1000

time (years)

p
o

p
u

la
tio

n
s

iz
e

N

We see that the population grows smoothly to the carrying capacity. At this point, the population growth
rate is zero and so the population remains at the carrying capacity. At this point, birth and death are
balanced and the population is stably living within the available resources.

Let’s also solve this using the ODE solver:
library(deSolve)

## Warning: package ‘deSolve’ was built under R version 4.0.5

5

parameters<-c(r=0.69,K=1000) #specify model parameters
state<-c(N=1) #initial state of the population

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

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

dN=r*N*(1-N/K)
return(list(c(dN)))

})
}

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

logGrowthout<-ode(y=state,times=times,func=LogisticGrowthmodel,parms=parameters)

par(oma=c(0,0,3,0))
plot(logGrowthout,xlab=”time”,ylab=”population size N”, main=”logistc growth model from ODE solver”)

0 5 10 15 20 25 30

0
2

0
0

6
0

0
1

0
0

0

logistc growth model from ODE solver

time

p
o

p
u

la
tio

n
s

iz
e

N

We see that the dynamics look identical to the explicit solution, as we expect.

6

Human Carrying Capacity

The carrying capacity is the maximum population level that can be sustained by available resources.

If you google “human carrying capacity”, you can find many interesting resources. A good book on this topic
is How Many People Can the Earth Support by Joel Cohen. Here is an excerpt from Wikipedia:

Several estimates of the carrying capacity have been made with a wide range of population
numbers. A 2001 UN report said that two-thirds of the estimates fall in the range of 4
billion to 16 billion with unspecified standard errors, with a median of about 10 billion.[5]
More recent estimates are much lower, particularly if non-renewable resource depletion and
increased consumption are considered.[6][7] Changes in habitat quality or human behavior at
any time might increase or reduce carrying capacity. In the view of Paul and Anne Ehrlich,
“for earth as a whole (including those parts of it we call Australia and the United States),
human beings are far above carrying capacity today.”

The application of the concept of carrying capacity for the human population has been crit-
icized for not successfully capturing the multi-layered processes between humans and the
environment, which have a nature of fluidity and non-equilibrium, and for sometimes being
employed in a blame-the-victim framework. Supporters of the concept argue that the idea of
a limited carrying capacity is just as valid applied to humans as when applied to any other
species. Animal population size, living standards, and resource depletion vary, but the concept
of carrying capacity still applies. The number of people is not the only factor in the carrying
capacity of Earth. Waste and over-consumption, especially by wealthy and near-wealthy peo-
ple and nations, are also putting significant strain on the environment together with human
overpopulation. Population and consumption together appear to be at the core of many hu-
man problems.Some of these issues have been studied by computer simulation models such as
World. When scientists talk of global change today, they are usually referring to human-caused
changes in the environment of sufficient magnitude eventually to reduce the carrying capacity
of much of Earth (as opposed to local or regional areas) to support organisms, especially Homo
sapiens.

The carrying capacity is the population size at which the population is in balance with its resources and
thus can be sustained indefinitely. For humans, this is not a fixed value. It depends on technology, life styles,
etc. Consider the carrying capacity for the contemporary US versus England in the Middle Ages. We in
the US consume vastly more resources than someone living in the middle ages did. At the same time, our
technology and economic system is far more efficient and advanced. We are able to grow vastly more food on
the same amount of land and we are extremely good at extracting energy from our enivronment. Given this,
our currenty carrying capacity appears to be much higher than what it would have been 1000 years ago.

The human carrying capacity of around 10 billion cited in the Wikipedia entry above is primarly based
on food production. However, remember that carrying capacity means the population number that can be
sustained indefinitely. Our current prosper is based heavily on energy production. Enery production is based
heavily on fossil fuels. Fossil fuels are a non-renewable resource that is being exhausted quickly. The fossil fuel
supply may last 50 years or it may last 200 years, but will definitely run out. Our current farming practices
are energy dependent. Thus, much depends on whether we can transition to renewable energy recourses.

7

Stabiltiy Analysis of the Logistic Growth Model

We saw when used the ODE solver with the logistic growth that the population went smoothly to the
equilibrium with is available resources. Will this always happen?

Consider the model again:

dN

dt
= rN(1 −

N

K
)

We see that this has a steady states atN = 0 and N = K:

dN

dt
= 0 ⇒ rN̂(1 −

K
) = 0 ⇒ N̂ = 0 or N̂ = K

If the population moves slightly above N̂ = 0, then

dN

dt
> 0

Thus, the population increases away from 0. This is an unstable steady state.

If the population is slightly below the carrying capacity, then (1 − N
K

) > 0 and dN
dt

> 0.

If the population is slightly below the carrying capacity, then (1 − N
K

) < 0 and dN
dt

< 0.
That is, if we move slightly away from the steady state in either direction, we move back to the steady state.
Thus, it is a stable equilibrium.

The condition for stability is that the system return to the steady state point when we move away. That is,
dN
dt

> 0 for N < N̂
dN
dt

< 0 for N > N̂

In other words, the condition for stability is that the slope of the right hand side of dN
dt

is negative at N̂:

if dN
dt

= f(N), then the steady state N̂ is stable if df
dN

< 0.

In the case of the logistic growth model, we have

f(N) = rN
(
1 −

N

K

)
df

dN
= r
(
1 −

N

K

)
−rN

( 1
K

)
= r
(
1 −

2N
K

)
At N̂ = 0, we have

df

dN
= r
(
1 −

2 ∗ 0
K

)
= r > 0

This is always positive and so the steady state at N̂ is unstable. This is simple to understand. When N = 0,
then the population is at a steady state because there is no one in the population and thus it can’t grow. If
move away from the steady state, this means that there are now some individuals in the population. Because
the population is low, then the reproductive rate is high and the population grows exponentially away from
N = 0.

At N̂ = K, we have
df

dN
= r
(
1 −

2 ∗K
K

)
= −r < 0

8

This is always negative and thus the steady state is always stable.

Discussion Question

1. Under the logistic growth model, the population always goes smoothly to its carrying capacity. Do you
think that it this is realistic? That is, will populations always grow smoothly to equilibrium with their
resources?

2. Do you think that the model is a realistic depiction of a growing population? List the assumptions of
the model. What shortcomings does it have?

9

Discrete Time Models

In order to better understand what is going on, we will introduce discrete time model. Humans have
overlapping generations. That is, there are people many different ages alive and reproducing at the same
time. Thus, population growth is a continuous process. In contrast, many species have approximately
non-overlapping generations. For example, in many insect species, most of the population emerges over a
short time frame in the spring time, reproduces,and then dies. This produces a new generation, which are all
born around the same, reproduce all at about the same time, and then die. This cycle continues until winter.
For such species, it is reasonable to view time as progressing in discrete increments t = 1, 2, 3, ….

The discrete time logistic growth model is

N(t + 1) = rN(t)
(
1 −

N(t)
K

)
That is, each individual in generation t produces r

(
1 − N(t)

K

)
offspring.

Discrete time models are easier to handle computationally than continous time ones:
endtime=100 #length of time to simulate
N<-vector(length=endtime)
N[1]=20
r=1.2
K=1000

for(gen in 2:endtime) #loop over generations
{N[gen]=r*N[gen-1]*(1-N[gen-1]/K)
}

plot(1:endtime,N,type=”l”,xlab=”generations”, ylab=”population size”,main=”discrete logistic growth”)

10

0 20 40 60 80 100

5
0

1
0

0
1

5
0

discrete logistic growth

generations

p
o

p
u

la
tio

n
s

iz
e

Long-term Behavior of the Discrete Logistic Model

In the case above, with the population growth rate set to 1.2 (meaning each individual produces an average
of 1.2 offspring in the next generation), the population goes smoothly to an equilibrium.

The equilibrium occurs when N(t+ 1) = N(t) – that is, the population size doesn’t change between generations
t and t + 1. This gives

N(t) = rN(t)
(
1 −

N(t)
K

)
This is satisfied by

N̂ = 0

and

1 = r
(
1 −

K

)
or

N̂ = K(1 −
1
r

)

For the above example, we have r = 1.2 and K = 1000. This yields

11

N̂ = 1000(1 −
1

1.1
) = 166.67

If the population growth rate r < 1, then we have

N̂ = K(1 −
1
r

) < 0

This steady state does not exist if r < 1. In this case the population dies out because reproduction is below
replacement rate.
endtime=100 #length of time to simulate
N<-vector(length=endtime)
N[1]=20
r=0.5
K=1000

for(gen in 2:endtime) #loop over generations
{N[gen]=r*N[gen-1]*(1-N[gen-1]/K)
}

plot(1:endtime,N,type=”l”,xlab=”generations”, ylab=”population size”,main=”discrete logistic growth”)

0 20 40 60 80 100

0
5

1
0

1
5

2
0

discrete logistic growth

generations

p
o

p
u

la
tio

n
s

iz
e

If the population growth rate is between 2 and 3, there are dampled oscillations to the equilibrium:

12

endtime=100 #length of time to simulate
N<-vector(length=endtime)
N[1]=20
r=2.9
K=1000

for(gen in 2:endtime) #loop over generations
{N[gen]=r*N[gen-1]*(1-N[gen-1]/K)
}

plot(1:endtime,N,type=”l”,xlab=”generations”, ylab=”population size”,main=paste(“discrete logistic growth, r=”,r))

0 20 40 60 80 100

0
1

0
0

3
0

0
5

0
0

7
0

0

discrete logistic growth, r= 2.9

generations

p
o

p
u

la
tio

n
s

iz
e

If the population growth rate is above 3, we get cycles:
endtime=100 #length of time to simulate
N<-vector(length=endtime)
N[1]=20
r=3.2
K=1000

for(gen in 2:endtime) #loop over generations
{N[gen]=r*N[gen-1]*(1-N[gen-1]/K)
}

plot(1:endtime,N,type=”l”,xlab=”generations”, ylab=”population size”,main=paste(“discrete logistic growth, r=”,r))

13

0 20 40 60 80 100

0
2

0
0

4
0

0
6

0
0

8
0

0
discrete logistic growth, r= 3.2

generations

p
o

p
u

la
tio

n
s

iz
e

endtime=100 #length of time to simulate
N<-vector(length=endtime)
N[1]=20
r=4
K=1000

for(gen in 2:endtime) #loop over generations
{N[gen]=r*N[gen-1]*(1-N[gen-1]/K)
}

plot(1:endtime,N,type=”l”,xlab=”generations”, ylab=”population size”,main=paste(“discrete logistic growth, r=”,r))

14

0 20 40 60 80 100

0
2

0
0

4
0

0
6

0
0

8
0

0
1

0
0

0
discrete logistic growth, r= 4

generations

p
o

p
u

la
tio

n
s

iz
e

Time Delay and Overshooting the Resources

We see that as the growth rate gets higher, there are very large cycles. The population size goes far above its
carrying capacity and then crashes to very low levels. This happens because the population grows past its
available resources.

In the above run, the equilibrium is at

N̂ = 1000(1 −
1
4

) = 750

When the population reaches levels approaching 1000, then it is far past the sustainable level. Starvation
(our some other consequence, depending on the limiting resource) ensues and the population crashes.

When the population growth rate is higher, then the population is able to grow further beyond the resource
before consequences set in.

The reason why the discrete model displays this characteristic and the continuous one doesn’t is that in the
continious model the population growth rate responds instantaneously to the current resource level. In the
discrete model, there is a time delay between the population growth and the consequence of the population
growth.

Consider generations 22 to 30 in the above run
N[22:30]

## [1] 967.33704 126.38436 441.64542 986.37897 53.74198 203.41513 648.14966
## [8] 912.20671 320.34250

15

At generation 23, the population size is 126. It grows to 441 in the next generation and is still below carrying
capacity.Its growth rate is

r
(
1 −

N

K
) = 4(1 −

441
1000

)
= 2.236

It grows to 986 in the next generation and is now above carrying capacity. Now the growth rate is

r
(
1 −

N(t)
K

)
= 4
(
1 −

986
1000

)
= 0.056

That is, the population drops to 54 (5.6% of its current level) because it overshot its resources.

Discussion Question

We have seen that delay in response to resources leads to unstability cycles in population models.

1. What would cycles like the last plot mean in the context of human populations?

2. Do you think that there is time delay in human response to resources?

16

Time Delay in continuous Models

Consider the situation in generation 24 in the example above. At this point, the population is well below its
carrying capacity and resources are plentiful. For this reason, everyone in the population reproduces at a
high rate. Unfortunately, this leaves the next generation well above carrying capacity and so there is mass
starvation and population crashes.

This does not happen in the continuous logistic model because reproduction occurs continuosly and responds
instantaenously to the current resources. It is possible to build a time delay in resource reponse into a
continuous model. For example

dN

dt
= rN(t)

(
1 −

N(t− τ)
K

)
Thus, the current population growth rate depends on what the population was at τ time units in the past.
Typcially, current reproduction does not depend only on the current state. Rather, it depends on the whole
lifetime experience of the organism. For an animal, current reproduction may depend on its overall state
of health and its available energy reserves for making offspring (reproduction takes a lot of energy for the
female).

In the case of humans, we are consuming resources (e.g. energy, fresh water, soil quality) that we are barely
aware. Their depeletion will probably not affect our behavior until a crisis stage is reached. Thus, there is
clearly a time delay in our response to resourcesm, although much more complex than the simple model
above.

Time delay models are beyond the scope of this course. We can get the essential ideas from discrete models
much more simply.

17

Chaos

Note that as the population growth rate increases, the dynamics get increasingly complicated:

0 10 20 30 40 50

0
3

0
0

7
0

0

discrete logistic growth, r= 3

generations

p
o

p
u

la
tio

n
s

iz
e

0 10 20 30 40 50

0
4

0
0

discrete logistic growth, r= 3.3

generations

p
o

p
u

la
tio

n
s

iz
e

0 10 20 30 40 50

0
4

0
0

discrete logistic growth, r= 3.6

generations

p
o

p
u

la
tio

n
s

iz
e

0 10 20 30 40 50

0
4

0
0

1
0

0
0

discrete logistic growth, r= 3.9

generations

p
o

p
u

la
tio

n
s

iz
e

At r = 3.3, we have regular cycles. At r = 3.6, we can discern periodic behavior, but it is complicated. There
are four peaks in each cycle. At r = 3.9 there is not any discernable pattern. The behavior looks random.
Let’s look more closely:

18

0 10 20 30 40 50

0
2

0
0

4
0

0
6

0
0

8
0

0
1

0
0

0
discrete logistic growth, r= 3.9

generations

p
o

p
u

la
tio

n
s

iz
e

There is no periodic (repeating) dynamics. Instead, the dynamics look random. This behavior is known as
chaotic dynamics. Even though the system is completely deterministic (i.e. no randomness in the data), it
looks like random fluctuations. Consider the next plot, in which I run the system three times with slightly
different starting population sizes:

19

0 10 20 30 40 50

0
2

0
0

4
0

0
6

0
0

8
0

0
1

0
0

0
discrete logistic growth, r= 3.9

generations

p
o

p
u

la
tio

n
s

iz
e

N0=10
N0=11
N0=12

The only difference in the three lines is the starting population sizes of 10,11, and 12. We see that the
population dynamics are in sync at first, but quickly diverge. If we look at any point beyond about 20
generations, the populations are at radically different sizes even though they started out very similar.

Contrast this to the case with r = 3.3:

20

0 10 20 30 40 50

0
2

0
0

4
0

0
6

0
0

8
0

0
discrete logistic growth, r= 3.3

generations

p
o

p
u

la
tio

n
s

iz
e

N0=10
N0=11
N0=12

In this case the populations syncronize and after 20 generations that have identical dynamics.

The Butterfly effect

With r = 3.3, differences in the initial conditions dissapear with time.

With r = 3.9, differences in the initial conditions grow with time

This sensitivity to initial conditions is a hallmark of chaotic dynamics. With chaotic dynamics, tiny
differences in initial conditions lead to large differences later in. This is often called the butterfly effect: “a
butterfly flapping its wings in Chicago leads to a cyclone in the Indian Ocean”. This is an exagerration, but
the basic idea is true: in chaotic systems, tiny effects can lead to big differences in the outcome. These effects
are random for practical purposes and it is very hard to tell random dynamics from chaotic ones.

Chaotic dynamics are characterized by large, random-seeming fluctuations. The phenomenon was first
observed in models of weather and was later found to exist in many types of dynamical systems. In population
models, it tends too occur when population growth rates are large. Large population growth rates lead to
large population fluctuations and large fluctuations lead to chaos.

Does Chaos Happen in Real Populations?

Chaos readily occurs in models, but does it occur in real populations? There has been a great deal of interest
in this question. Consider Figures 1 and 2 below. These show scatter plots of estimates of the Lyapunov
number for a selection of laboratory (Figure 1) and field (Figure 2) populations.

21

Figure 1: Figure from Ellner and Turchin, American Naturalist, 1995, vol 145(3)

The only thing that you need to know about the Lyapunov number is

γ < 1 ⇒ non-chaotic dynamics

γ > 1 ⇒ chaotic dyamics

We see from the figures that very few populations appear to exhibit chaotic dyamics. However, many
populations are right at the boundary of chaos.

Evolution of Population Dynamics

Why should it be that many populations are at the boundary of chaos? The answer to this is unknown.
However, we can speculate as to the reason.

Evolution works on reproduction. When a trait leads to higher survival and reproduction, then it is favored by
evolution and tends to spread in a population. All else being equal, it is beneficial to an organism to produce
more offspring. Of course “all else” is usually not equal and there are complex tradeoffs in determining the
optimal strategy for an organism in terms of reproduction. Still, producing more offspring will often be a
beneficial strategy and thus evolution will often favor it.

22

Figure 2: Figure from Ellner and Turchin, American Naturalist, 1995, vol 145(3)

23

This leads to a problem: Evolution works on individuals. An individual’s evolutionary interests are often
best served by maximizing the number of offspring. However, high reproduction leads to unstable population
dynamics, which is bad for the population.

Large fluctuations in population numbers can easily lead to population extinction. When the population
number gets very low, some chance event (disease, bad weather, etc) can kill enough of the remaining
members that the population goes extinct. Thus, natural selection on individuals may lead to extinction of
the population.

This leads to multiple hypotheses about why we see many populations at the boundary of chaos, but few over:

1. Populations that cross the boundary into chaos tend to become extinct. Thus, the distribution of
Lyapunov numbers that we observe is the result of a distribution of optimal individual reproductive
numbers. However, the high end of the distribution gets truncated by population extinction.

2. Evolution favors being at the edge of chaos. That is, it might be most beneficial to have the highest
reproductive rate that doesn’t lead to chaos. Thus, natural selection would favor reproductive rates
right at the boundary. However, although this seems intuitively appealing, it has a major problem:
Natural selection works on individuals, not groups. Suppose that we are in a situation where the average
reproductive rate in the population is just below the threshold for chaos. An individual would still
benefit by producing more offspring. If some individal has a mutation that leads to more offspring, then
that mutation will tend to be represented more in future generations and it will thus spread and lead
to higher reproductive rates. Even though the population as a whole will better off if each individual
limits reproduction, individuals gain short-term benefit by “cheating” and reproducing more. Various
researchers have suggested ways in which individual natural selection could lead to stable population
dyanamics (my very first scientific publication was on this topic), but it remains an open question.

24

Dynamics of Multiple Populations: Predator-Prey Interactions.

Figure 3 shows a famous data set in population ecology. This shows the number of hare and lynx pelts
traded by the Hudson Bay Company over a period of 100 years. We see a strong cyclic pattern in both
populations, with a period of about 10 years. The lynx cycle lags behind the hare cycle by perhaps 1-2 years.
The lynx is a major predator of the hare and the hare is a major food source for the lynx. The explanation
for the observed pattern is that 1) in absence of many lynx, the hare population increases rapidly; 2) as the
hare population increases, the lynx population also increases because there is plentiful food for the lynx; 3)
As the lynx population begins to get large, they eat too many hare and the hare population collapses; 4)
When the hare population collapses, the lynx population collapses soon after because there is not enough
food; 5) Once the lynx population collapses, the hare population then starts to increase again. 6) Repeat.
Ecologists are interested in understanding what drives these dynamics. There has been a large amount of
work in developing mathematical models for population dynamics and then fitting these models to data.

Assumptions

Take X to represent the prey (e.g. hare) population number and Y to represent the predator (e.g. Lynx)
population number.

We will assume
1. In absence of predators, the prey population follows the logistic growth model (and thus follows the
assumptions of that model).
2. The prey is the only food source for the predator. Thus, in absence of the prey, the predator population
decreases at rate r2.
3. The prey has only one predator.
4. The loss in prey biomass due to predation is given by α1XY1+βX . This will be explained below.
5. The gain in predator biomass due to predation is given by α2XY1+βX .
6. There is no randomness in the amount of reproduction, predation, etc.
7. The enivronment is constant. There are no changes in prey resource availability, mortality, etc.
8. There is no structure to either population. There is no variation due to age, sex, genetics, location, etc.

25

With these assumptions, we have

dX

dt
= r1X

(
1 −

X

K

)

α1XY

1 + βX

dY

dt
= −r2Y +

α2XY

1 + βX

r1
(
1 − X

K

)
is the growth rate of the prey population in absence of the predator (logistic growth)

−α1XY1+βX is the decrease in the prey biomass due to predation

−r2Y is the is the (negative) growth rate of the predator population in absence of the prey (they die out
with no food).
α2XY
1+βX is the increase in predator biomass due their predation on prey.

Functional Response

The term α2XY1+βX is called the functional response. This determines the amount of prey eaten as a function
of prey density.

The function that we use here is called a Holling Type II functional response. Holling was a scientist who did
early work on the concept of functional response. Plot the function:
maxX=100 #max number of pret eaten
halfmax=50 #prey number at half max eaten

beta=1/halfmax
alpha=beta*maxX

X=seq(from=0,to=1000,by=1)
plot(X,alpha*X/(1+beta*X),type=”l”, ylab=”prey eaten”, xlab=”prey number”,main=”Functional Reponse”)

26

0 200 400 600 800 1000

0
2

0
4

0
6

0
8

0
Functional Reponse

prey number

p
re

y
e

a
te

n

When the prey number is low, then the number eaten per predator increases linearly – that is if there is twice
the prey, then twice the prey is eaten. As the …

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)

1

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:

dS

dt
= −βSI + ξR

dI

dt
= βSI −γI

dR

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)

2

(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:

dS

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

SI

R
=
ξ

β

dI

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

dR

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.

3

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

4

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

S

time

0 20 40 60 80 100

0
1

0
2

0
3

0
4

0
5

0

I

time

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

R0 =

γ
> 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

5

0 50 100 150 200

1
6

0
2

0
0

2
4

0
2

8
0

S

time

0 50 100 150 200
2

5
3

0
3

5
4

0
4

5
5

0

I

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 =

γ

Here are estimated values of R0 for several pandemics:

6

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

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?

7

Numerically Solving Systems of Ordinary Differential Equations:

Consider our SIR Model:

dS

dt
= −βSI + ξR

dI

dt
= βSI −γI

dR

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

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)))

})
}

8

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

0

S

time

0 10 20 30 40 50

0
4

0
0

I

time

0 10 20 30 40 50

0
4

0
0

R

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 =

γ
> 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

9

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

0

S

time

0 10 20 30 40 50

0
6

I

time

0 10 20 30 40 50

0
3

6

R

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]))

10

0 1000 3000 5000

9
7

5
1

0
0

0

S

time

0 1000 3000 5000

0
6

I

time

0 1000 3000 5000

0
1

0

R

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)

11

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

0

S

time

0 10 20 30 40 50

5
0

2
0

0

I

time

0 10 20 30 40 50

0
2

5
0

R

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]))

12

0 10 20 30 40 50

0
8

0
0

S

time

0 10 20 30 40 50

0
6

0
0

I

time

0 10 20 30 40 50

0
4

0
0

R

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]))

13

0 10 20 30 40 50

0
8

0
0

S

time

0 10 20 30 40 50

0
6

0
0

I

time

0 10 20 30 40 50

0
6

0
0

R

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

dx

dt
= xy −y

dy

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”)

14

#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))
}

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

2

flow field for example system

x

y

x nullclines
y nullclines

15

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

dx

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

dy

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

dx

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

dy

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

dx

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

dy

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)
{

16

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

17

##
## 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

x

time

0.0 0.5 1.0 1.5 2.0

0
2

0
6

0
1

0
0

y

time

18

#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)

19

0 2 4 6 8 10

0
.0

0
.1

0
.2

0
.3

0
.4

0
.5

x

time

0 2 4 6 8 10
0

.0
0

.1
0

.2
0

.3
0

.4
0

.5

y

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)

20

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

0 2 4 6 8 10

1

0
0

0
0

6

0
0

0

2
0

0
0

x

time

0 2 4 6 8 10

0
.5

0
.6

0
.7

0
.8

0
.9

1
.0

y

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

21

dS

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

dI

dt
= βSI −γI −νI

dR

dt
= γI − ξR−νR

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

Infectious Disease Dynamics with Immunization

Age Structure

The final thing that we will do with epidemic models is to study the effect of immunization. In order to do
this, we will need a model that has age structure. Our previous models do not track age. They treat everyone
as identical in the population. Howeer, the dynamics of vaccination are highly dependent on age structure.

We can write differential equation models that include age structure. However, this would require (partial
differential equations). This a more advanced area that we won’t study in this class. Instead, we will modify
the simulation model.

Previously, the population was represented with a vector pop that stored the disease status (susceptible,
infected, recovered) for each member of the population. Now, we will represent the population with a matrix,
that stores disease status, age, age at infection, and sex:

The function below initializes the population:
#function that generates poisson(lambda) values until it gets a non-zero value
trunc_pois<-function(lambda)
{ d=0

while(d==0) d=rpois(1,lambda)
return(d)

}

FindWorkPlace<-function()
{

}

#Make random households
#F0 fanilies with kids. All have two parents, make and fenale
#S0 households adults only of random size

make_households<-function(pop,H0)
{ for(f in 1:F0)

{ ageM=runif(20,50) #qge of mother
ageD=ageM+rnorm(0,5) #age of father
nKids=trunc_pois(1)

mom=c(0,ageM,NA,1,f,)

}

1

}

initialize_pop<-function(S0,I0,maxAge)
{#initialize the population. pop0 is matrix with pop0[i,1]=susceptble/infected/recovered status

#pop0[i,2] is current age. pop0[i,3] will be assigned age at infection, pop0[i,4]=1 if female

pop0=matrix(nrow=S0+I0,ncol=4)

#assigned susceptible/infected status:
pop0[1:S0,1]=0
pop0[(S0+1):(S0+I0),1]=1
pop0[(S0+1):(S0+I0),3]=0 #assign age at infection

#assign ages. all ages have equal probability.
pop0[,2]=sample(maxAge*365,S0+I0,replace=T)

#assign sex (only females reproduce)
pop0[,4]=ifelse(runif(nrow(pop0))<=0.5,1,0) #1 if female, 0 if male

return(pop0)

}

Now, we can have important parameters (infection rates, death rates, etc) vary by age. The functions nelow
define age-specific rates. Each function takes age a as input and outputs the appropriate rate. Note that
these currently are not very realistic. I didn’t use any data in setting the rates and instead just took very
simple values. However, we can easily incorporate detailed data.
InfectionRatebyAge<-function(a)
{ #function returns the infection rate for an individual of age a

# age is in days

if(a<=5*365) { #less than 5 years old
return(0.015)
} else if(a<=10*365) { #5-10 years old

return(0.012)
} else if(a<=69*365) {#10-69 years old

return(0.01)
} else{ #> 69

return(0.02)
}

}

RecoveryRatebyAge<-function(a)
{ #function returns the infection recovery rate for an individual of age a

# age is in days

2

if(a<=5*365) {
return(0.1)

} else if(a<=10*365) {
return(0.1)

} else if(a<=69*365) {
return(0.1)

} else{
return(0.1)

}

}

BirthRate<-function(a)
{#function returns the birth rate for an individual of age a

# age is in days
#currently set so that women from 18 to 45 have same fertility rate
#if less than 18 or more than 45, fertility is zero.
#this is not very realistic. We can improve with data
if(a<=18*365){

return(0)
} else if(a<=45*365) {

return(1/27)
} else {

return(0)
}

}

BackgroundMortalitybyAge<-function(a)
{ #function returns the infection rate for an individual of age a

if(a<=5*365) {
return(0.015/365)

} else if(a<=10*365) {
return(0.015/365)

} else if(a<=69*365) {
return(0.015)

} else{
return(0.05/365)

}

}

The functions below do all do the updating of the population. The functions infect, recover,
backgroundmortality, and birth work similarly to how they did previously. However, they are now
modified to account for age structure. The functions aging, clear_dead and cap_pop are new: aging ages
the population by one day every day andclear_dead removes dead people from the population each day.
The function cap_pop puts a cap on the population. If the population exceeds maxpop at the end of the day,
then the function randomly selects people and removes them from the population. This function is required
to prevent the population from growing exponetially and making the program run very slowly. Ideally, we
would have death rates and birh rates balanced so that the population stays constant. However, this is
difficult to acheive in practice and so we force a cap on the population.

3

aging<-function(pop)
{#function to age everyone by one day

pop[,2]=pop[,2]+1

}

#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,InfectRate)
{ nI<-sum(pop[,1]==1) #number of infecteds

nS<-sum(pop[,1]==0) #number of susceptibles
pop1=pop #create copy of pop

susceptibles=which(pop[,1]==0) #make vector giving positions of suscepibles

beta=sapply(pop[susceptibles,2],InfectRate) #sapply applies the function to the vector
#produces a vector of infection probabilities.
#where beta[i] is infection prob for ith susceptible

pi=1-(1-beta)^nI #probility of susceptible being infected. Calcs seperate prob for each individual, based on
# that individuals age-specific infection prob.

#pi=ifelse(beta*nI<=1,beta*nI,1)

pop1[susceptibles,1]=ifelse(runif(nS)<pi,1,0) #check if each susceptible becomes infected

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,RecRate)
{

nI<-sum(pop[,1]==1) #number of infected
pop1=pop

infecteds=which(pop[,1]==1)

gamma=sapply(pop[infecteds,2],RecRate) #sapply applies the function to the vector
#produces a vector of infection probabilities.

pop1[infecteds,1]=ifelse(runif(nI)<gamma,2,1) #check if each infected recovers. This is the “coin flip”

return(pop1)

4

}

#function to check whether each individual dies from non-disease causes
#MortRate is a function that gives the mortality rate for each age
backgroundmortality<-function(pop,MortRate)
{

pop1=pop
n=nrow(pop1) #current population size

#note that this doesn’t distinguish between living and dead. That is, dead people can die again. This does not
#affect anything.

nu=sapply(pop[,2],MortRate) #produces a vector of mortality rates for everyone in the population
#This is speficic to their age

pop1[,1]=ifelse(runif(n)<=nu,3,pop1) #check if each person dies . If so, write 3, if not write current value

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,birthrate)
{ #function to check for births

#birthrate is a function that gives fertility rate by age.

females<-pop[pop[,4]==1,] #find the females
nf<-nrow(females) #number of females

nbaby<-sum(runif(nf)<sapply(females[,2],birthrate)) #number of babies born today

if(nbaby>0) #if any babies are born today, then add them to the population
{ babies<-matrix(nrow=nbaby,ncol=4)

babies[,1]=0 #babies born susceptible (ignoring possibility of temporary immunity from mother)
babies[,2]=0 #age 0
babies[,4]=ifelse(runif(nbaby)<=0.5,1,0) #assign sex
return(rbind(pop,babies))

}

return(pop) #if no babies, then return pop as is

5

}

{#function to remove dead from population

return(pop[-which(pop[,1]==3),])
}

cap_pop<-function(pop,maxpop)
{#this function caps the population. If the population size exceeds maxpop, it randomly selects (n-maxpop) individuals
# and removes them from the population

n<-nrow(pop)

if(n<=maxpop) return(pop)

remove<-sample(1:n,(n-maxpop),replace=F)

return(pop[-remove,])

}

Finally, here is the main program. This follows the same format as the model without age structure, but
with functions that account for age structure.
#initialize the population
pop=initialize_pop(S0,I0,max_age)

#initialize the output variables:
St=S0
It=I0
Rt=0
Mt=0
At=NULL
Vt=0
Qt=NULL
RI=NULL #will store proportion of

#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)
{ if(time/(10*365)==floor(time/(10*365))) print(time/365) #print time every ten years to track progress

pop=infect(pop,InfectionRatebyAge)
pop=recover(pop,RecoveryRatebyAge)
pop=backgroundmortality(pop,BackgroundMortalitybyAge)
pop=diseasemortality(pop,m)
pop=birth(pop,BirthRate)
pop=aging(pop)

6

St=c(St,sum(pop[,1]==0))
It=c(It,sum(pop[,1]==1))
Rt=c(Rt,sum(pop[,1]==2))
Mt=c(Mt,sum(pop[,1]==3))
At=c(At,mean(pop[,3],na.rm=T))
Vt=c(Vt,sum(pop[,5]==1))

#Qt=c(Qt,St*mean(sapply(pop[,2],InfectionRatebyAge))/mean(sapply(pop[,2],RecoveryRatebyAge)))

pop=cap_pop(pop,maxpop)

if(time/(10*365)==floor(time/(10*365))) if(sum(pop[,1]==1)<1) pop[1,1]=1 #if infected drops below 1 add 1 every ten years.

}

Here is a run of the model with age structure. We see that it goes through large amplitude cycles, with
disease outbreaks about once every 20 years. In each outbreak, most of the population is infected (we can
tell this because the number of susceptibles gets close to zero). The cycles are irregular both in amplitude
and period, suggesting the possibility of chaotic dyanmics (we will discuss chaos in the context of population
dynamics, our next topic). They don’t show any sign up damping over the 1000 years plotted (in fact, I ran
it for 2000 years total and there was not sign of damping over that range either). I am not sure why the
cycles are so much bigger with the age stuctured population.

7

0 200 400 600 800 1000

0
5

0
0

1
5

0
0

2
5

0
0

time (years)

n
u

m
b

e
r

Susceptible
Infected
Recovered

Immunization

Next, we will consider the effect of immunization.

From Wikipedia:

A vaccine is a biological preparation that provides active acquired immunity to a particular disease. A vaccine
typically contains an agent that resembles a disease-causing microorganism and is often made from weakened
or killed forms of the microbe, its toxins, or one of its surface proteins. The agent stimulates the body’s
immune system to recognize the agent as a threat, destroy it, and to further recognize and destroy any of the
microorganisms associated with that agent that it may encounter in the future. Vaccines can be prophylactic
(example: to prevent or ameliorate the effects of a future infection by a natural or “wild” pathogen), or
therapeutic (e.g., vaccines against cancer are being investigated).

The administration of vaccines is called vaccination. Vaccination is the most effective method of preventing
infectious diseases; widespread immunity due to vaccination is largely responsible for the worldwide eradication
of smallpox and the restriction of diseases such as polio, measles, and tetanus from much of the world. The
effectiveness of vaccination has been widely studied and verified; for example, the influenza vaccine, the HPV
vaccine, and the chicken pox vaccine. The World Health Organization (WHO) reports that licensed vaccines
are currently available for twenty-five different preventable infections.

The widespread adoption of vaccines has greatly reduced the incidence of many diseases. Smallpox, formerly
one of the largest killers of humans, has been completely eradicated from the world (there are samples at two
locations in the world: the CDC in Atlanta, and a government lab in Russia).

In the case of measles, there were an estimated 2-3 million deaths world wide each year before the vaccine was
introduced. That number has dropped to under 100,000 per year (e.g. an estimated 90,000 measles deaths in

8

2016). These deaths are predominantly in developing countries where there are not widespread immunization
programs.

In the US, there were 700,000 measles cases in the US in 1958 (before the adoption of the measles vaccine),
resulting in 558 deaths. Since the vaccine came into wide use, the number of cases dropped to under 100 per
year (Figure 1).

The table in Figure 2 summarizes measles cases in 2008:

9

reference:https://www.cdc.gov/mmwr/preview/mmwrhtml/mm5718a5.htm

We see that most cases were from unvaccinated cases. Although measles is very rare in the US, it is still
fairly common in some developing countries. All or nearly all outbreaks in the US start with people traveling
from other countries.

Figure 4 below shows measles in recent years (source=CDC). Unfortunately, measles has been on the increase
because of the unfortunate trend of parents not vaccinating their children.

Vaccine Schedule

In order to maximize protection, it is recommended that children be vaccinated as soon their immune
system is sufficiently developed to respond to a particular vaccine. This has led to a complicated schedule of
recommended ages for vaccinations. The link below shows the schedule for the US.

10

Children in the US typically receive the vaccines during their regularly checkups at the pediatrician. The
pediatrician cannot force children to get vaccines. However, all 50 states have laws requiring children to be
vaccinated before entering kindergarten in public schools.

Unfortunately, there has been a trend towards parents not vaccinating their children in the US and other
western countries. This is based on 1) mistaken beliefs about the safety of vaccines and 2) people not under-
standing how devastating these diseases once were for children before vaccines were introduced. Fortunately,
most parents still get their children vaccinated and these diseases are largely kept in check. However, there
have been outbreaks of measles, whooping cough, and other diseases in parts of the country where large
numbers of parents don’t vaccinate their kids.

We will apply immunization using the function below. For simplicity, we assume that children are vaccinated
exactly on their 5th birthday. We define a vaccination rate ir. On each child’s 5th birthday, we “flip the
coin” to determine if they are vaccinated. If so, there disease status is changed to 2. The status of 2 formerly
represented recovered. Now it represents immune whether by recovering from the disease or from vaccination.
We assume that the vaccine is 100% effective and that immunity is never lost.
immunize<-function(pop,ir)
{ #function to immunize. Assumes that kids are immunized exactly on their five-year-old birthday

#ir is immunization rate

immunize_today<-which(pop[,2]==365*5) #find kids with 5-year-old birthday today
nb<-length(immunize_today) #number of kids with 5-year old birthday today

if(nb>0) #if there are any with 5-year birthday today, then check whether immunuized.
{ pop[immunize_today,]=ifelse(runif(nb)<ir,2,pop[immunize_today,]) #if vaccinated, then change status to immunized
} #else leave as is

return(pop)

}

Here is a run with 90% of children immunized at age 5. The model runs for 100 years before the immunization
program begins. Because only 5-year-olds are vaccinated, it takes decades before a large portion of the
population is immunized. Once this happens, the disease outbreaks are greatly reduced. After the introduction
of the vaccine, there are three outbreaks in 200 years. All three are of much smaller magnitude than what
was observed before the vaccine was introduced.

11

0 50 100 150 200 250 300

0
5

0
0

1
0

0
0

2
0

0
0

90% immunization at age 5

time (years)

n
u

m
b

e
r Susceptible

Infected
Recovered/Immune
vaccinated

Next, let’s see what happens as the proportion vaccinated is decreased. The next plot shows the same run,
but with vaccination rates of 80% rather than 90%. Now we see that the disease outbreaks are larger in
magnitude and possibly more frequent. It is easier to see the magnitude of the outbreak in the changes in the
number of susceptibles than in the number infected. There are four large outbreaks after the introduction of
the vaccine. These are much smaller in magnitude than the pre-vaccine (Where nearly everyone was infected),
but still there are large numbers of people getting the disease.

12

0 50 100 150 200 250 300

0
5

0
0

1
0

0
0

2
0

0
0

80% immunization at age 5

time (years)

n
u

m
b

e
r Susceptible

Infected
Recovered/Immune
vaccinated

The next plot shows the case when the vaccination rate is only 50%. Now, the vaccine fails to control the
disease. The size of the outbreaks is slightly reduced, but the effect is minimal.

13

0 50 100 150 200 250 300

0
5

0
0

1
0

0
0

2
0

0
0

50% immunization at age 5

time (years)

n
u

m
b

e
r Susceptible

Infected
Recovered/Immune
vaccinated

Herd Immunity

Look back at the first plot. Note that after vaccination has become common that the proportion of Susceptibles
and Recovered/Immune was approximately the same as before the epidemic (e.g.the number of susceptible
after the vaccine introduction is similar to the average across cycles prior to the introduction). The manner
in which people gain immunity changes: before vaccination they gained immunity through getting the disease.
After the vaccine, they gain immunity by being vaccinated. However, the numbers are similar.

The key point: Even though the number of susceptibles in the population is about the same
as pre-vaccine, the susceptibles mostly don’t catch the disease anymore. This is because the
incidence of the disease is much lower in the population. Even though people are susceptible,
they don’t get the disease because they don’t encounter infected people. This phenomenon is
known as herd immunity

Here is the definition if herd immunity from Wikipedia:

Herd immunity (also called herd effect, community immunity, population immunity, or social immunity)
is a form of indirect protection from infectious disease that occurs when a large percentage of a population
has become immune to an infection, thereby providing a measure of protection for individuals who are not
immune.In a population in which a large number of individuals are immune, chains of infection are likely
to be disrupted, which stops or slows the spread of disease. The greater the proportion of individuals in a
community who are immune, the smaller the probability that those who are not immune will come into contact
with an infectious individual.

Herd immunity is a key part of the benefit of vaccination programs. Because of this, parents choosing to not
vaccinate their children has impact beyond their own children. There are many people in the population
who cannot be vaccinated. This includes kids that are too young to receive the vaccine and people who with

14

compromised immune systems who are unable to be vaccinated. When vaccination rates are high in the
rest of the population, then these people are protected by herd immunity. However, when vaccination rates
decrease, then herd immunity is lessened and disease outbreaks can occur. We see in the simulations above
that the size of disease outbreaks was much higher for 80% vaccination rate than for 90%. That is, even if
only 10% of the population decide to not vaccinate their children then many more people will get the disease.
Typically, babies and young children are hardest hit because they are too young for vaccination.

Mathematical models have played a large role in epidemiology and particularly in the area of vaccination
policy. As we have seen the dynamics are complicated and often non-intuitive. The use of models allows
us to explore the ramifications of different vaccination strategies without having to try them out in a real
population.

In order to eradicate a disease. the reproductive rate must be reduced to less than one. Suppose that before
the vaccination campaign, that the equilibrium number of susceptibles is S∗. Then, the disease effective
intrinsic reproductive rate is

Re =
S∗β

γ + µ

If the vaccination rate is v, then the number of susceptibles is reduced by the proportion v and Re becomes

Re =
(1 − p)S∗β
γ + µ

If Re is reduced below one after the vaccination campaign, then the disease will be eradicated.

Consider the disease simulated in the last three figures. Before the vaccination begins, there are outbreaks
when the number of susceptibles reaches about 1800.

The mean infection rate is about 0.05/365 (it varies by age, but is near 0.05 per year) and the recovery rate is
0.1 for all age groups. The background mortality rate is 0.015 except for individuals over 70. Thus, we have

Re =
S∗β

γ + µ

1800 ∗ 0.05/365
0.1 + 0.015

= 2.14

Children are vaccinated at 5 years old and are susceptible up until that point. The proportion of children
sum(pop[,2]<5*365)/3000

## [1] 0.1473333

Thus, the total proportion of vaccinated people in the population is

proportion over 5 ∗ vaccination rate

With a 90% vaccination rate for 5-year olds, we have

Re = 2.14 ∗ (1 − 0.9 ∗ 0.85) = 0.5

With an 80% vaccination rate, we have

Re = 2.15 ∗ (1 − 0.8 ∗ 0.85) = 0.42

15

With an 50% vaccination rate, we have

Re = 2.15 ∗ (1 − 0.5 ∗ 0.85) = 1.23

Re is below one for 90% and 80% vaccination rates. In both of these cases, the disease was largely eliminated
except for small outbreaks coming from new input of infecteds. With a vaccination rate of 50%, Re > 1 and
we see that the disease persists in the population with major outbreaks continuing

When do epidemics end?

An epidemic ends when Re drops below one. That is, when infected people on average infect less than one
other person.

Recall this expression for Re:

Re =
S

N

γ + µ
=

γ + µ

Key points:

1. The effective disease reproductive rate Re is decreased when susceptibles become a lower proportion of
the population.

2. The proportion of people who are susceptible drops by people becoming immune.

3. People become immune either by getting the disease or by being vaccinated.

The epidemic ends when the number of people who become immune by getting the disease
plus the number of people who become immune by being vaccinated becomes large enough
that infected people will not encounter large numbers of susceptible people.

The exact percentage of immune people that is required depends on the details of the disease.

Another key point:

If there is not vaccination (and nothing else changes) then there are likely to be further
outbreaks. This is because there will be an inflow of new susceptibles by birth. If they are
not vaccinated, then the number of susceptibles will eventually become high enough to fuel a
new outbreak. Alternatively, the disease will become endemic: always in the population at a

If we don’t want this, then

GET VACCINATED

Age at Infection

The impacts of vaccination are heavily impacted by the age at which infections typically occur. In order to
be effective, the vaccine has to be administered at an earlier age than most people will get the disease. With
many of the historically deadly diseases (e.g. measles), people would most often catch them in childhood.
This is not necessarily because children are more susceptible. If infection rates are high, then most people
will catch it within a few years of exposure. For example, suppose that the probability of infection is 20% per
year for people of all ages. Then most kids will catch the disease before they are five – not because they are
more vulnerable than anyone else to infection, but because it is not possible last longer than a few years
without catching it. For a disease with lower infection rates, the average age of infection can be much older.

16

Consider the plot below. I set the infection rate to 5% per year for all ages and the recory rate is 10% per
day. The bottom plot shows the average age at infection over the population. We see that it varies between
10 and 18 years of age as the intensity of the epidemic varies. When the number of susceptibles is high, then
children are catching the disease around 10-12 years old. After an outbreak when the number of suscepibles
is low, the average age of infection increases to 16-18. In all cases, it is predominant children who catch the
disease, even though the probability of disease transmission is the same for all ages. This is because older
people predominantly already caught the disease when they were young and thus are immune.

200 250 300 350 400 450 500

0
1

5
0

0

disease dynamics, beta=0.05/365 for all ages

time (years)

n
u

m
b

e
r

200 250 300 350 400 450 500

1
0

1
4

1
8

Average Age at Infection

a
ve

ra
g

e
a

g
e

(
ye

a
rs

)

Vaccination and Age of Infection

Now, let’s consider what happens when we introduce vaccination:

The first plot below has vaccination at age 5. In this case, the average age at infection drops after the
vaccination campaign begins. This is because most infections happen in children who are not vaccinated yet.

17

200 300 400 500 600

0
1

5
0

0
vaccinated at age 5

time (years)

n
u

m
b

e
r

200 300 400 500 600

0
1

5
3

0

Average Age at Infection

a
ve

ra
g

e
a

g
e

(
ye

a
rs

)

In the next simulation, children are vaccianted at birth. Before the vaccination campaign, the average age
of infection is around six. After the vaccination campaign, the average age of infection increases to the
range 14-18. This happens because of herd immunity. Because there are many fewer infected people in the
population than before the vaccination, those who are not vaccinated encounter many fewer infected people.
Thus, they are typically teenagers when they get the disease.

18

200 300 400 500 600

0
1

5
0

0
vaccinated at birth

time (years)

n
u

m
b

e
r

200 300 400 500 600

6
1

2
1

8

Average Age at Infection

a
ve

ra
g

e
a

g
e

(
ye

a
rs

)

Vaccines can strongly impact the age at which people tend to get the disease and frequently
will increase the average age at which someone gets the disease

Is this necessarily a good thing? Consider the case of rubella:

Case Study: Rubella

From Wikipedia:

Rubella, also known as German measles or three-day measles, is an infection caused by the rubella virus
This disease is often mild with half of people not realizing that they are infected. A rash may start around
two weeks after exposure and last for three days. It usually starts on the face and spreads to the rest of the
body.The rash is sometimes itchy and is not as bright as that of measles.Swollen lymph nodes are common
and may last a few weeks. A fever, sore throat, and fatigue may also occur. In adults joint pain is common.
Complications may include bleeding problems, testicular swelling, and inflammation of nerves. Infection
during early pregnancy may result in a child born with congenital rubella syndrome (CRS) or miscarriage.
Symptoms of CRS include problems with the eyes such as cataracts, ears such as deafness, heart, and brain.
Problems are rare after the 20th week of pregnancy. Rubella is usually spread through the air via coughs of
people who are infected.People are infectious during the week before and after the appearance of the rash.
Babies with CRS may spread the virus for more than a year.[1] Only humans are infected. Insects do not
spread the disease. Once recovered, people are immune to future infections. Testing is available that can
verify immunity. Diagnosis is confirmed by finding the virus in the blood, throat, or urine. Testing the blood
for antibodies may also be useful.[1] Rubella is preventable with the rubella vaccine with a single dose being
more than 95Rubella is a common infection in many areas of the world. Each year about 100,000 cases of
congenital rubella syndrome occur. Rates of disease have decreased in many areas as a result of vaccination.
There are ongoing efforts to eliminate the disease globally. In April 2015 the World Health Organization

19

declared the Americas free of rubella transmission. The name “rubella” is from Latin and means little red.
It was first described as a separate disease by German physicians in 1814 resulting in the name “German
measles.”

Children in the US receive the MMR (measles-mumps-rubella) vaccine. The first dose is at around one year
old and the second dose is in the range 4-6 years old.

A key characteristic of rubella is that it is usually mild, except when pregnant mothers transmit it to their
babies in early pregnancy. The resulting congenital rubella syndrome can be devastating in the newborn
baby. Thus, the most important goal in controlling rubella is preventing women of child-bearing age from
becoming infected.

Rubella vaccine camapigns commenced in the 1970’s. In a seminal 1983 paper, Anderson and May (Vaccination
against rubella and measles: quantitative investigations of different policies) explored the implications of
different vaccination policies for rubella.

In the US, the initial rubella vaccination campaign took the strategy of vaccinating all children 15 months or
older. The …

Project #1 – STAT 4660
Fall 2022

Due Friday March 4 by midnight

In this this project you will conduct a simulation study of different COVID mitigation strategies and
use it to make a case for what you think the best strategy would have been. You will simulate a town of
population approximately 10,000 people for several years (or however long it takes for the epidemic to
run its course). In your town, people should live in households with other people and go to work or
school. This will allow you to model disease transmission dynamics in a more realistic way than the
previous models that we have studied. You should have a moderately realistic structure to your town in
terms of household composition and number and size of businesses. You introduce the disease into the
town in a small number of individuals and then see how the disease spreads. You should test at least the
following scenarios:

1) No mitigation strategies are attempted.
2) Total lockdown with no one leaving home.
3) Lockdowns with e.g. 90%, 80%, …50%,… of people staying home.
4) Measures such as masking and social distancing. You can model these simply as reductions in

transmission probabilities.
5) School closing only. Businesses stay open.
6) Businesses only closing. Schools stay open.
7) You should consider these scenarios with and without a vaccine.

You should consider metrics such as total deaths, total cases, number of people with long
COVID, etc. It might make sense to consider not just deaths, but years of life lost (i.e. a child dying is
many more years of life lost than an 80-year-old and this is relevant).
The various parameters in your model should be based on real data. This includes:

– Household composition (i.e. number of households with one person, two people, two people
with kids, how many kids, etc).

– Infection rates in different settings: home, school, work.
– Mortality rates by age.
– Length of time that someone is infectious.
– Probability of long COVID.
You should write a report that explains what you did, your results with figures/tables, and your

conclusions about the best COVID strategy.
In order to get a good grade, you must do the following:
– Clearly explain your assumptions (i.e. the rules of how your simulation works).
– Make a reasonable effort to have realistic assumptions.
– Have a thorough discussion of caveats and shortcomings of your model and results.
– Have your model parameters based on data and reference the source of that data.
– Have code that is well organized and documented.
– Have appropriate figures and tables.
– Make a good case for your preferred mitigation strategy based on model output. We are

modeling the disease side of things in detail, but not modeling things like economic impacts, lost

education opportunities, mental health impacts, etc. Thus, you can make qualitative arguments

Hints:
1. Your program should follow the same basic form as the IBM model in the lecture notes. That is,

there should be a matrix or data frame that tracks the state of the population, with rows
corresponding to members of the population and columns corresponding to variables that you
want to track. There should be a loop over time and in each time step various functions are
applied to the population that implement processes such as infection, recovery, etc.

2. My population data frame has columns of disease status, age, age at infection, sex, household
index, workplace/school, employee index/school grade, day of recovery if infected, day of
death if they die from disease, indicator long COVID. My program does not actually use age at
infection or sex, but I kept these from the previous program in case I want them in the future.
You do not need to have exactly the same columns. I am telling you this to give an example of
what I did.

3. You should make a function to construct your population. You will need to do this many times
and you will want to experiment with different population compositions. Having a function will
make this more convenient. My function is of form make_households(F0,nw,sb,lb,nsb,nlb),
where F0 is the number of households, nw is the number of people who don’t work, sb is the
number of small businesses with nsb employees, and lb is the number of large businesses with
nlb employees. The function loops over desired households. For each one, it randomly
determines the household composition (single individual, couple, couple with kids, etc) and the
number of children (if there are any) based on real data for the US. I have another function that
takes the values nw,sb,lb,nsb, and nlb and constructs the available jobs. It then randomly
assigns adults to jobs. It does similarly with schools for children. It calculates grade in school
based on age.

4. Households, jobs, and schools only matter in the model because they determine which
members of the population interact with each other and therefore can infect each other. Rates
of infection may vary between different settings.

5. I have separate functions for disease transmission in households, businesses, schools, and in the
community. I also have a functions that implement mitigation measures such as lockdowns.
Putting things in functions allows you to easily change them in the future. For example, if I
decide that I want to have household disease transmission work in a different way, it just means
writing a new function. Provided that the input and output is the same, I don’t need to change
anything else in the program.

6. There are many different pieces to this program, but work on one thing at a time. Keep testing
as you go along. Check at every step to make sure that your code is doing what it is supposed to.

7. The first thing that you probably should do is the make the function that constructs your
population because you can’t test anything else until you have a correct population data frame.

8. You will get most of the data for model parameters from scientific papers. One way to find these