please follow instructions attached lectures
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
d
N
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.
6
9
per year (the amount needed to double once per year):
1
N0=2
r<-c(1)
t=seq(from=0,to=
10
0,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
4
0 60
8
0 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
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
5
00 1000
15
00 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
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 paramete
rs
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
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
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 −
N̂
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
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 −
N̂
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 1
26
.38436 441.64542 986.37897 53.74
19
8 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
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
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
d
X
dt
= r1X
(
1 −
X
K
)
−
α1X
Y
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 number of prey keeps on increasing, then the number eaten per
predator begins to flatten out. This is because the predators are getting enough to and so don’t eat more as
the number of prey increases.
We assume that the functional reponse is the same between predators and prey except for the constants α1
and α2. The ratio α2α1 determines how efficiently prey biomass is converted into predator biomass).
Dynamics of the Model
Let’s run the model using the ODE solver:
library(deSolve)
parameters<-c(r1=1,K=500,alpha1=0.01,beta=0.005,r2=1,alpha2=0.01) #specify model parameters state<-c(X=200,Y=10) #initial state of the population
PredatorPreymodel<-function(t,state,parameters) {
X=state[1]
Y=state[2]
with(as.list(c(state,parameters)),{
dX=r1*X*(1-X/K)-alpha1*X*Y/(1+beta*X)
dY=-r2*Y+alpha2*X*Y/(1+beta*X)
27
return(list(c(dX,dY)))
})
}
times<-seq(from=0,to=200,by=0.01)
PredPreyout<-ode(y=state,times=times,func=PredatorPreymodel,parms=parameters)
par(mfrow=c(1,1))
matplot(PredPreyout[,1],PredPreyout[,2:3],type=”l”,xlab=”time”, ylab=”population numbers”, main=”Predator-Prey Model”)
0 50 100 150 200
0
1
0
0
2
0
0
3
0
0
4
0
0
Predator−Prey Model
time
p
o
p
u
la
tio
n
n
u
m
b
e
rs
In this case, the predator and prey populations undergo damped oscillations to an equilibrium. If we increase
the carrying capacity to K = 625, the we get cycles:
library(deSolve)
parameters<-c(r1=1,K=625,alpha1=0.01,beta=0.005,r2=1,alpha2=0.01) #specify model parameters state<-c(X=200,Y=10) #initial state of the population
PredatorPreymodel<-function(t,state,parameters) {
28
X=state[1]
Y=state[2]
with(as.list(c(state,parameters)),{
dX=r1*X*(1-X/K)-alpha1*X*Y/(1+beta*X)
dY=-r2*Y+alpha2*X*Y/(1+beta*X)
return(list(c(dX,dY)))
})
}
times<-seq(from=0,to=200,by=0.01)
PredPreyout<-ode(y=state,times=times,func=PredatorPreymodel,parms=parameters)
par(mfrow=c(1,1))
matplot(PredPreyout[,1],PredPreyout[,2:3],type="l",xlab="time", ylab="population numbers", main="Predator-Prey Model")
0 50 100 150 200
0
1
0
0
2
0
0
3
0
0
4
0
0
5
0
0
Predator−Prey Model
time
p
o
p
u
la
tio
n
n
u
m
b
e
rs
This shows the phase plane analysis:
29
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”)
lim=500
PredPreyflow<-flowField(PredatorPreymodel,xlim=c(0,lim), ylim=c(0,lim), parameters=parameters,add=FALSE,xlab="X",ylab="Y")
PredPreyNullclines<-nullclines(PredatorPreymodel,xlim=c(0,lim), ylim=c(0,lim), parameters=parameters,state.names = c("X","Y"))
PredPreyTrajectory<-trajectory(PredatorPreymodel,y0=matrix(c(25,200,150,150,450,200,150,50,400,400,200,25),nrow=6,byrow=T),parameters=parameters,tlim=c(0,20))
## Note: col has been reset as required
0 100 200 300 400 500
0
1
0
0
2
0
0
3
0
0
4
0
0
5
0
0
X
Y
X nullclines
Y nullclines
We see that with these parameter values there is a single stable equilibrium. The system will osciallate around
this equilibrium wherever it starts from (except if either population starts at zero).
Using models such as this one, we can explore how the dynamics change as we vary different aspects of the
model. In particular, ecologists are interested in what characteristics will cause cycles versus more stable
dynamics and how well the correspondance is between the models and real populations.
30
Now we will turn out attention to a more specialized example of predation: large scale fishing. Let’s start
with the concept of a fishery. Here is a definition from https://www.sciencedaily.com/terms/fishery.htm
Fishery A fishery (plural: fisheries) is an organized effort by humans to catch fish or other aquatic species, an
activity known as fishing. Generally, a fishery exists for the purpose of providing human food, although other
aims are possible (such as sport or recreational fishing), or obtaining ornamental fish or fish products such
as fish oil. Regardless of purpose, however, the term fishery generally refers to a fishing effort centered on
either a particular ecoregion or a particular species or type of fish or aquatic animal, and usually fisheries are
differentiated by both criteria. Examples would be the salmon fishery of Alaska, the cod fishery off the Lofoten
islands or the tuna fishery of the Eastern Pacific. Most fisheries are marine, rather than freshwater; most
marine fisheries are based near the coast.This is not only because harvesting from relatively shallow waters is
easier than in the open ocean, but also because fish are much more abundant near the coastal shelf, due to
coastal upwelling and the abundance of nutrients available there.
Many fisheries have collapsed due to overfishing.
Stock Peak Catch (year) Catch in 1981
Antartic Blue Whale 29,000 whales (19
31
) Nil
Antartic fin whale 27,000 whales (1938) Nil
Hokkaido herring 850,000 tons (1913) Nil
Peruvian anchoveta 12.3 million tons(1970) 0.3 million tons
SW African pilchard 1.4 million tons (1962) Nil
Nort Sea herring 1.5 million tons (1962) negligible
California sardine 640,000 tons (1936) Nil
Georges Bank herring 374,000 tons (1968) Nil
Japanse sardine 2.3 million tons (1939) 17,000 tons
When a fishery collapes, it is a disaster both environmentally (for the fish) and economically (for the fishermen
and consumers of the fish)
Even when fisheries do not collapse, fishing is ofen on economically marginal occupation. In other words,
most fishermen are poor.
Goal: Manage fisheries in way that is environmentally and economically sustainable, provides
good work for fishermen, and provides a good source to consumers.
Mathematical models have played a large role in fisheries management.
###Modeling of Fisheries
31
https://www.sciencedaily.com/terms/fishery.htm
- Logistic Growth Model
Solving the Logistic Growth Model
Human Carrying Capacity
Stabiltiy Analysis of the Logistic Growth Model
Discrete Time Models
Long-term Behavior of the Discrete Logistic Model
Time Delay and Overshooting the Resources
Time Delay in continuous Models
Chaos
The Butterfly effect
Does Chaos Happen in Real Populations?
Evolution of Population Dynamics
Dynamics of Multiple Populations: Predator-Prey Interactions.
Assumptions
Functional Response
Dynamics of the Model
Dynamics of Large Scale Fishing
Epidemic Model
s
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
.
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?
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.
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:
d
S
dt
= −βSI + ξ
R
d
I
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?
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 eventuall
y
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̄ = Ī
ξ
γ
Then, we have
S̄ + Ī + R̄ = N
or
γ
β
+ Ī + Ī
ξ
γ
= N
Solving for Ī, we get
Ī = γ
N − γ
β
γ + ξ
and
R̄ = Ī
ξ
γ
= ξ
N − γ
β
γ + ξ
To summarize, we have
S̄ =
γ
β
Ī = γ
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, Ī1) = (N, 0)
(S̄2, Ī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
8
0
10
0
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 =
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
5
0 50 100
15
0 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
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:
6
image from https://www.the-scientist.com/features/why-r0-is-problematic-for-predicting-
covid-1
9
-spread-6
7
690
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.
7
https://www.the-scientist.com/features/why-r0-is-problematic-for-predicting-covid-
19
-spread-67690
https://www.the-scientist.com/features/why-r0-is-problematic-for-predicting-covid-19-spread-67690
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
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
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
https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve
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,
26
7 infecteds, and 5
33
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
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.
16
36
002 0.3
27
5
22
3
## [499997,] 4999.96 999.5089 0.16
35
998 0.
32
752
17
## [499998,] 4999.97 999.5089 0.1635995 0.32752
11
## [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 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
d
x
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))
}
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
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.
39
789, R2 = 1.945
31
e-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.611
37
e-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.
28
813e-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]))
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 leave the system by death.
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 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
.
22
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=
29
0,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
0
S
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)
23
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
S
time
0 100 300 500
0
2
4
6
8
I
time
24
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
0
S
time
0 100 300 500
0
2
0
4
0
6
0
8
0
1
2
0
I
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):
25
0 100 300 500
1
0
0
2
0
0
3
0
0
4
0
0
5
0
0
S
time
0 100 300 500
0
2
0
6
0
1
0
0
I
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.
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
26
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
s
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
27
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
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.
28
0 50 100 150 200
0
.2
0
.8
time (years)
su
sc
e
p
tib
le
s
0 50 100 150 200
0
.0
0
0
.1
0
0
.2
0
time (years)
in
fe
ct
e
d
s
29
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)
}
#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) return(pop1) } #function to check whether each individual dies from non-disease causes pop1=pop return(pop1) #function to check whether each individual dies from disease nI<-sum(pop==1) #number of infected
pop1=pop
pop1[pop==1]=ifelse(runif(nI) 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: #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 pop=recover(pop,gamma) pop=backgroundmortality(pop,nu) pop=birth(pop,b) 31 Mt=c(Mt,sum(pop==3)) #time=1000 #legend(“right”,c(“Susceptible”,”Infected”,”Recovered”,”Died”),lty=1:4,col=1:4)
time=1000 0 200 400 600 800 1000
0 8 0 n m e Infected We see that the disease persists in the population at a low level. The fluctuations up and down are from Note that I balanced the birth and death rate (birth=2xdeath because only females reproduce) so that the Now, I will add in disease mortality. I take a mortality rate of 0.001 per unit time. Since the disease lasts Note that I don’t need to rerun the function definitions (infect, recover, etc) 32 M0=0 #initial number dead pop=c(rep(0,S0),rep(1,I0),rep(2,R0)) pop=recover(pop,gamma) St=c(St,sum(pop==0)) } time=1000 33 0 200 400 600 800 1000 0 time This doesn’t change things very much because the disease morality is low. We do see a slight downward Now I will increase mortality:
34 0 200 400 600 800 1000 0 time e Susceptible The disease mortality rate is now 0.1 per time unit. Thus, the overall disease morality rate is close to 100%. If we make the birth rate higher, we can have the disease persist with high mortality because new susceptibles 35 0 200 400 600 800 1000 0 time Note that the recovered population is growing. This is because births are greater than deaths. The disease 36 0 200 400 600 800 1000 Coding Details of the Simulation.
In order to understand better how the simulation program works, lets look at the infect function: nS<-sum(pop==0) #number of susceptibles
pop1=pop
pi=1-(1-beta)^nI #probility of susceptible being infected } In order to demonstrate the code, we will initialize the population and also define beta. I am making the The first two lines in the function are These count the number of infecteds and susceptibles. The code pop==1 creates a vector of TRUE and FALSE
37 values: TRUE where the value in pop is equal to one and FALSE where it is not equal to one: ## [1] 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 ## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE We call a vector of TRUE and FALSE values a logical vector. Applying the sum command counts the number of ## [1] 10 Next, we calculate the probability that a susceptible is infected: ## [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 ## [1] 0.33223562 0.01786395 0.65783307 0.970 43 755 0. 38 175202 0.59197094 The command set.seed makes it so that the same random number sequene is generated each time. This is Taking runif(nS) ## [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 We use the ifelse command to modify the pop vector. The command ifelse(condition,A,B) takes a 38 set.seed(3332) ## [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: We first take pop1=pop. This creates a copy of pop.
Then pop1[pop==0] extracts the 0 elements from pop1: ## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE ## [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 Now, pop1 has been modified so that newly infected individuals have been changed from 0 to 1. The other 1. Find out the number of individuals of a relevant type (e.g. suscepibles and infecteds in function 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. The main part of the program is pop=recover(pop,gamma) This is a loop over time. That is, for (time in 1:endtime) steps through times and executes the commands In each time step, a function for each possible event is applied to the population: 39 pop=backgroundmortality(pop,nu) Then, the current state of the population is recorded: 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 0 20 40 60 80 100 fr ct n f o u tio susceptible 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 40 0.0 0.5 1.0 1.5 2.0 2.5
0 We see that the infection dies out. Why is this different from the ODE model? Consider detailed output ## time S I 42 4.6930 3.884276e-05 41 ## [22,] 7.20 438.3983 2.585547e-05 We see that the number of infecteds drops very low (less than $10ˆ{-4}). However, there cannot really be In order for a disease to persist, the population must be sufficiently large to hold a resevoir If the population is large, then that tiny fraction of the population will be small number of infected individuals How do we deal with this problem in the simulation? One solution would be to make the population size 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. #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 pop=infect(pop,beta) 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)
42
#legend(“right”,c(“Susceptible”,”Infected”,”Recovered”,”Died”),lty=1:4,col=1:4) #SIR2=cbind(St,It,Rt)
#save(SIR2,file=”workspaces/SIR2.RData”)
load(“workspaces/SIR2.RData”)
endtime=200*365 0 50 100 150 200 We now see cycles like we expect. This requirement for a minimum population is not something that we can We will use this model more in our next section of notes.
43 Infectious Disease Dynamics with The final thing that we will do with epidemic models is to study the effect of immunization. In order to do We can write differential equation models that include age structure. However, this would require (partial Previously, the population was represented with a vector pop that stored the disease status (susceptible, The function below initializes the population: )
{ d= 0
while(d==0) d=rpois( 1 ,lambda)
}
FindWorkPlace<-function()
{
} #Make random households make_households<-function(pop,H0)
{ for(f in 1:F0)
{ ageM=runif( 2 0, 5 0) #qge of mother mom=c(0,ageM,NA,1,f,)
} } 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: #assign ages. all ages have equal probability. 6 5,S0+I0,replace=T)
#assign sex (only females reproduce) return(pop0)
} Now, we can have important parameters (infection rates, death rates, etc) vary by age. The functions nelow # age is in days
if(a<=5*365) { #less than 5 years old
return(0.0
15 ) 10 *365) { #5-10 years old
return(0.0 12 ) 9 *365) {#10-69 years old
return(0.01) return(0.02) } RecoveryRatebyAge<-function(a)
{ #function returns the infection recovery rate for an individual of age a
# age is in days if(a<=5*365) {
return(0.1)
} else if(a<=10*365) {
return(0.1)
} else if(a<=69*365) {
return(0.1)
} else{ } BirthRate<-function(a)
{#function returns the birth rate for an individual of age a
# age is in days 8 to 45 have same fertility rate 18 or more than 45, fertility is zero. return(0) return(1/2 7 ) 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{ } The functions below do all do the updating of the population. The functions infect, recover, 3 aging<-function(pop)
{#function to age everyone by one day
pop[,2]=pop[,2]+1
} #Function to make new infections: 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 pi=1-(1-beta)^nI #probility of susceptible being infected. Calcs seperate prob for each individual, based on #pi=ifelse(beta*nI<=1,beta*nI,1)
pop1[susceptibles,1]=ifelse(runif(nS) return(pop1)
}
#function to make individuals recover 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 pop1[infecteds,1]=ifelse(runif(nI) return(pop1) } #function to check whether each individual dies from non-disease causes pop1=pop #note that this doesn’t distinguish between living and dead. That is, dead people can die again. This does not nu=sapply(pop[,2],MortRate) #produces a vector of mortality rates for everyone in the population 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 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) if(nbaby>0) #if any babies are born today, then add them to the population babies[,1]=0 #babies born susceptible (ignoring possibility of temporary immunity from mother) } return(pop) #if no babies, then return pop as is
5 } clear_dead<-function(pop)
{#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 #initialize the output variables: #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 pop=infect(pop,InfectionRatebyAge) 6 St=c(St,sum(pop[,1]==0)) #Qt=c(Qt,St*mean(sapply(pop[,2],InfectionRatebyAge))/mean(sapply(pop[,2],RecoveryRatebyAge)))
pop=clear_dead(pop) 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 20 years. In each outbreak, most of the population is infected (we can 7 0 200 400 600 800 1000
0
0 1 0 2 0 time (years)
n m e Susceptible 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 The administration of vaccines is called vaccination. Vaccination is the most effective method of preventing The widespread adoption of vaccines has greatly reduced the incidence of many diseases. Smallpox, formerly In the case of measles, there were an estimated 2-3 million deaths world wide each year before the vaccine was 8
20 16 ). These deaths are predominantly in developing countries where there are not widespread immunization In the US, there were 700,000 measles cases in the US in 19 58 (before the adoption of the measles vaccine), 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 Figure 4 below shows measles in recent years (source=CDC). Unfortunately, measles has been on the increase Vaccine Schedule
In order to maximize protection, it is recommended that children be vaccinated as soon their immune 10 https://www.cdc.gov/mmwr/preview/mmwrhtml/mm5718a5.htm https://www.cdc.gov/vaccines/schedules/hcp/imz/child-adolescent.html
Children in the US typically receive the vaccines during their regularly checkups at the pediatrician. The Unfortunately, there has been a trend towards parents not vaccinating their children in the US and other We will apply immunization using the function below. For simplicity, we assume that children are vaccinated #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. return(pop) Here is a run with 90% of children immunized at age 5. The model runs for 100 years before the immunization 11 https://www.cdc.gov/vaccines/schedules/hcp/imz/child-adolescent.html 0 50 100 150 200 250 300
0 1 0 2 0 90% immunization at age 5
time (years) e Infected Next, let’s see what happens as the proportion vaccinated is decreased. The next plot shows the same run, 12 0 50 100 150 200 250 300 80% immunization at age 5
time (years) The next plot shows the case when the vaccination rate is only 50%. Now, the vaccine fails to control the 13 0 50 100 150 200 250 300 50% immunization at age 5
time (years) Look back at the first plot. Note that after vaccination has become common that the proportion of Susceptibles The key point: Even though the number of susceptibles in the population is about the same Here is the definition if herd immunity from Wikipedia:
Herd immunity (also called herd effect, community immunity, population immunity, or social immunity) Herd immunity is a key part of the benefit of vaccination programs. Because of this, parents choosing to not 14 compromised immune systems who are unable to be vaccinated. When vaccination rates are high in the Mathematical models have played a large role in epidemiology and particularly in the area of vaccination In order to eradicate a disease. the reproductive rate must be reduced to less than one. Suppose that before Re = γ + µ
If the vaccination rate is v, then the number of susceptibles is reduced by the proportion v and Re becomes
Re = 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 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 Re = γ + µ 1800 ∗ 0.05/365 = 2.14
Children are vaccinated at 5 years old and are susceptible up until that point. The proportion of children ## [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 An epidemic ends when Re drops below one. That is, when infected people on average infect less than one Recall this expression for Re:
Re = N
Nβ
γ + µ Sβ
γ + µ Key points:
1. The effective disease reproductive rate Re is decreased when susceptibles become a lower proportion of 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 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 If we don’t want this, then
GET VACCINATED
The impacts of vaccination are heavily impacted by the age at which infections typically occur. In order to 16 Consider the plot below. I set the infection rate to 5% per year for all ages and the recory rate is 10% per 200 250 300 350 400 450 500
0 5 0 disease dynamics, beta=0.05/365 for all ages
time (years) 1 1 Average Age at Infection
a ra e g ( a ) 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 17 200 300 400 500 600
0 0 time (years) 5 0 In the next simulation, children are vaccianted at birth. Before the vaccination campaign, the average age 18 200 300 400 500 600 0 time (years) 6 2 8 Vaccines can strongly impact the age at which people tend to get the disease and frequently Is this necessarily a good thing? Consider the case of rubella:
From Wikipedia: Rubella, also known as German measles or three-day measles, is an infection caused by the rubella virus 19 declared the Americas free of rubella transmission. The name “rubella” is from Latin and means little red. Children in the US receive the MMR (measles-mumps-rubella) vaccine. The first dose is at around one year A key characteristic of rubella is that it is usually mild, except when pregnant mothers transmit it to their Rubella vaccine camapigns commenced in the 1970’s. In a seminal 1983 paper, Anderson and May (Vaccination In the US, the initial rubella vaccination campaign took the strategy of vaccinating all children 15 months or Contrast this to the United Kingdom. There, the vaccination strategy was to vaccinate only girls between 20 Why did this happen? It is because of the impact on age at infection of the two policies. The table shown in In many countries the average age of infection before vaccination was less than ten years old. This is well 21 lower end of the child-bearing years. Of course, this is average age at infection – some people will get it later. The UK policy targeted only girls shortly before they would reach the age of bearing children. This has the Note the very young age of infection it the Gambia (2-3 years). Ironically, Gambia had an advantage over more Anderson and May used mathematical models to study the impact of vaccination strategies on the incidence 22 If the average age of infection is low (which will happen when the disease incidence is high) and if the Decreasing the incidence of disease increases the average age of infection. Because rubella In contrast, the UK strategy of vaccinating at age 15 has no downside. Since rubella usually has minimal 23 Many other disease (e.g. measles and chicken pox) are also worse in adults than in children. However, unlike The Anderson and May paper had a large impact on understanding of vaccine policy. This is a an example 24 Project #1 – STAT 4660 Due Friday March 4 by midnight In this this project you will conduct a simulation study of different COVID mitigation strategies and 1) No mitigation strategies are attempted. transmission probabilities. You should consider metrics such as total deaths, total cases, number of people with long – Household composition (i.e. number of households with one person, two people, two people – Infection rates in different settings: home, school, work. conclusions about the best COVID strategy. 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: there should be a matrix or data frame that tracks the state of the population, with rows 2. My population data frame has columns of disease status, age, age at infection, sex, household 3. You should make a function to construct your population. You will need to do this many times 4. Households, jobs, and schools only matter in the model because they determine which 5. I have separate functions for disease transmission in households, businesses, schools, and in the 6. There are many different pieces to this program, but work on one thing at a time. Keep testing 7. The first thing that you probably should do is the make the function that constructs your 8. You will get most of the data for model parameters from scientific papers. One way to find these 9. One advantage of using functions is that you can start with parameter values that you make up 10. There is partial credit! Your program does not have to be perfect to get points.
30
backgroundmortality<-function(pop,nu)
{
n=length(pop1)
pop1=ifelse(runif(n)
}
diseasemortality<-function(pop,m)
{
}
}
St=S0
It=I0
Rt=R0
Mt=M0
#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=diseasemortality(pop,m)
St=c(St,sum(pop==0))
It=c(It,sum(pop==1))
Rt=c(Rt,sum(pop==2))
}
#matplot(1:time,cbind(St,It,Rt,Mt)[1:time,],type=”l”,xlab=”time”,ylab=”number”,lty=1:4,col=1:4)
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)
2
0
0
4
0
0
6
0
0
0
1
0
0
0
time
u
b
r Susceptible
Recovered
random fluctuations from day to day in what happens.
population will remain approximately constant.
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.
set.seed(1984)
S0=1000 #initital number of susceptibles
I0=10 #initial number of infecteds
R0=0 #initital number of recovereds
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[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
#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=backgroundmortality(pop,nu)
pop=diseasemortality(pop,m)
pop=birth(pop,b)
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)
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
2
0
0
4
0
0
6
0
0
8
0
0
1
0
0
disease mortality=0.01
n
u
m
b
e
r Susceptible
Infected
Recovered
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.
0
2
0
0
4
0
0
6
0
0
8
0
0
1
0
0
disease mortality=0.1
n
u
m
b
r
Infected
Recovered
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.
keep coming in to the population:
0
2
0
0
4
0
0
6
0
0
8
0
0
1
0
0
disease mortality=0.1, birth rate=0.035
n
u
m
b
e
r Susceptible
Infected
Recovered
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
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
infect<-function(pop,beta)
{ nI<-sum(pop==1) #number of infecteds
pop1[pop==0]=ifelse(runif(nS)
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
nI<-sum(pop==1) #number 1of infecteds
nS<-sum(pop==0) #number of susceptibles
pop
pop==1
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
TRUE values:
sum(pop==1)
pi=1-(1-beta)^nI
pi
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)
## [7] 0.39927833 0.68927095 0.44402178 0.50352531
for illustrative purposes for these notes. We do not have this in the program.
determine whether an individual is infected. In this example, the second individual is infected and the rest
are not.
logical vector condition and returns A for every element where Condition is TRUE and B for every element
whereCondition is FALSE:
ifelse(runif(nS)
pop1=pop
pop1[pop==0]=ifelse(runif(nS)
pop1=pop
pop==0
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
pop1[pop==0]
occupied by 0’s.
functions all work in a similar way. They each follow the steps
infect).
4. Modify pop for the individuals for whom the event occured.
for (time in 1:endtime)
{ pop=infect(pop,beta)
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))
}
inside the brackets for each time step.
pop=infect(pop,beta)
pop=recover(pop,gamma)
pop=diseasemortality(pop,m)
St=c(St,sum(pop==0))
It=c(It,sum(pop==1))
Rt=c(Rt,sum(pop==2))
Mt=c(Mt,sum(pop==3))
run below:
0
.0
0
.2
0
.4
0
.6
0
.8
1
.0
time (years)
a
io
o
p
p
la
n
infected
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).
1
0
0
0
3
0
0
0
5
0
0
0
time (years)
n
u
m
b
e
r Susceptible
Infected
Recovered
from the ODE model:
SIRout[700:725,]
## [1,] 6.99
## [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
## [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
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:
of infecteds when the disease reaches the low point of its cycle.
and then the disease can persist in small pockets. This is one reason why cities are often necessary for
epidemics to occur.
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.
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
#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)
pop=recover(pop,gamma)
pop=backgroundmortality(pop,nu)
pop=diseasemortality(pop,m)
pop=birth(pop,b)
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)
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
1
0
0
0
3
0
0
0
5
0
0
0
time (years)
n
u
m
b
e
r Susceptible
Infected
Recovered
learn from the ODE model.
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
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.
differential equations). This a more advanced area that we won’t study in this class. Instead, we will modify
the simulation model.
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:
#function that generates poisson(lambda) values until it gets a non-zero value
trunc_pois<-function(lambda
return(d)
#F0 fanilies with kids. All have two parents, make and fenale
#S0 households adults only of random size
ageD=ageM+rnorm(0,5) #age of father
nKids=trunc_pois(1)
1
pop0[1:S0,1]=0
pop0[(S0+1):(S0+I0),1]=1
pop0[(S0+1):(S0+I0),3]=0 #assign age at infection
pop0[,2]=sample(maxAge*3
pop0[,4]=ifelse(runif(nrow(pop0))<=0.5,1,0) #1 if female, 0 if male
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
} else if(a<=
} else if(a<=6
} else{ #> 69
}
2
return(0.1)
}
#currently set so that women from 1
#if less than
#this is not very realistic. We can improve with data
if(a<=18*365){
} else if(a<=45*365) {
} else {
}
return(0.05/365)
}
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.
#”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
#produces a vector of infection probabilities.
#where beta[i] is infection prob for ith susceptible
# that individuals age-specific infection prob.
#”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)
{
#produces a vector of infection probabilities.
4
#MortRate is a function that gives the mortality rate for each age
backgroundmortality<-function(pop,MortRate)
{
n=nrow(pop1) #current population size
#affect anything.
#This is speficic to their age
}
diseasemortality<-function(pop,m)
{
}
{ babies<-matrix(nrow=nbaby,ncol=4)
babies[,2]=0 #age 0
babies[,4]=ifelse(runif(nbaby)<=0.5,1,0) #assign sex
return(rbind(pop,babies))
}
with functions that account for age structure.
#initialize the population
pop=initialize_pop(S0,I0,max_age)
St=S0
It=I0
Rt=0
Mt=0
At=NULL
Vt=0
Qt=NULL
RI=NULL #will store proportion of
#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=recover(pop,RecoveryRatebyAge)
pop=backgroundmortality(pop,BackgroundMortalitybyAge)
pop=diseasemortality(pop,m)
pop=birth(pop,BirthRate)
pop=aging(pop)
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))
pop=cap_pop(pop,maxpop)
disease outbreaks about once every
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.
5
0
5
0
5
0
u
b
r
Infected
Recovered
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).
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.
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).
introduced. That number has dropped to under 100,000 per year (e.g. an estimated 90,000 measles deaths in
programs.
resulting in 558 deaths. Since the vaccine came into wide use, the number of cases dropped to under 100 per
year (Figure 1).
fairly common in some developing countries. All or nearly all outbreaks in the US start with people traveling
from other countries.
because of the unfortunate trend of parents not vaccinating their children.
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.
pediatrician cannot force children to get vaccines. However, all 50 states have laws requiring children to be
vaccinated before entering kindergarten in public schools.
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.
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
{ pop[immunize_today,]=ifelse(runif(nb)
}
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.
5
0
0
0
0
0
0
n
u
m
b
r Susceptible
Recovered/Immune
vaccinated
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.
0
5
0
0
1
0
0
0
2
0
0
0
n
u
m
b
e
r Susceptible
Infected
Recovered/Immune
vaccinated
disease. The size of the outbreaks is slightly reduced, but the effect is minimal.
0
5
0
0
1
0
0
0
2
0
0
0
n
u
m
b
e
r Susceptible
Infected
Recovered/Immune
vaccinated
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.
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
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.
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
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.
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.
the vaccination campaign, that the equilibrium number of susceptibles is S∗. Then, the disease effective
intrinsic reproductive rate is
S∗β
(1 − p)S∗β
γ + µ
when the number of susceptibles reaches about 1800.
0.1 for all age groups. The background mortality rate is 0.015 except for individuals over 70. Thus, we have
S∗β
≈
0.1 + 0.015
under 5 is about 15%:
sum(pop[,2]<5*365)/3000
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
other person.
S
=
the population.
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.
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
steady low level.
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.
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.
1
0
n
u
m
b
e
r
200 250 300 350 400 450 500
1
0
4
8
ve
g
a
e
ye
rs
vaccination campaign begins. This is because most infections happen in children who are not vaccinated yet.
1
5
0
vaccinated at age 5
n
u
m
b
e
r
200 300 400 500 600
0
1
3
Average Age at Infection
a
ve
ra
g
e
a
g
e
(
ye
a
rs
)
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.
0
1
5
0
vaccinated at birth
n
u
m
b
e
r
200 300 400 500 600
1
1
Average Age at Infection
a
ve
ra
g
e
a
g
e
(
ye
a
rs
)
will increase the average age at which someone gets the disease
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
It was first described as a separate disease by German physicians in 1814 resulting in the name “German
measles.”
old and the second dose is in the range 4-6 years old.
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.
against rubella and measles: quantitative investigations of different policies) explored the implications of
different vaccination policies for rubella.
older. The trend over time in cases of rubella and CRS are shown in Figure 4. We see that the number of
rubella cases declined after the introduction of the vaccine, as we would expect. However, the number of
CRS cases actually increased as the vaccination campaign took hold even as the number of rubella cases
dropped very low:
ages 10-15. Figure 5 shows the trend in cases of CRS in the UK. We see that the number of cases decreased.
Figure 6 gives the average age of rubella infection for various countries.
before child-bearing age, and so it would not lead to CRS. In the US, the average age was 10-11 before
vaccination. After vaccination, the average age of infection increased to 13-16 years. This is approaching the
Thus, the US policy greatly reduced the incidence of rubella, but the people who did get it tended to get it
later. This was problematic because it meant more pregnant women were getting rubella and transmitting it
to their children.
advantage that it doesn’t greatly impact the force of disease infection. That is, there are still many people
with rubella in the population and thus those who get rubella will still tend to get it young. Note that that
the average age of infection was 9-10 in the UK even after the vaccination campaign.
devloped countries: because the disease was so prevalent, nearly everyone caught it very young. Therefore,
CRS was very rare.
of CRS. The figure below is from the Anderson and May paper and is based on their mathematical model.
The solid lines show the ratio of CRS cases after and before vaccination under the US policy of vaccination
versus the proportion p who are vaccinated. The different curves correspond to different values of A, the
average age of infection before vaccination. The dashed line shows the same ratio under the UK vaccination
policy.
proportion vaccinated is not high (e.g. 80%) then the US vaccination policy can actually increase the incidence
of CRS:
is most harmful in pregnant women (by causing CRS), the US vaccination strategy will be
counterproductive if 1) The disease prevalence if high before vaccination and 2) the vaccination
rate doesn’t approach 100%
effect on children, there is not a large cost to vaccinating later. Vaccinating later has the big benefit that it
doesn’t change the “force of infection” and therefore the age at which people get rubella. Currently, nearly
all children in the US get the MMR (measles-mumps-rubella) vaccine and so the vaccination rate approaches
100% and the policy is effective.
rubella, many such diseases are also harmful in children and so there is a cost to vaccinating later.
of the usefulness of mathematical models. The dynamics are too complex to work out without the use of
models and we don’t want to use real populations as guinea pigs for vaccination strategies. Epidemiological
models have played a large role in our understanding of infectious diseases and play a large role in public
health policy.
Age Structure
Immunization
Vaccine Schedule
Herd Immunity
Eradicating Diseases
When do epidemics end?
Age at Infection
Vaccination and Age of Infection
Case Study: Rubella
Fall 2022
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:
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
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.
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:
with kids, how many kids, etc).
– 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
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 code match your stated assumptions.
– 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
about these things.
1. Your program should follow the same basic form as the IBM model in the lecture notes. That is,
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.
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.
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.
members of the population interact with each other and therefore can infect each other. Rates
of infection may vary between different settings.
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.
as you go along. Check at every step to make sure that your code is doing what it is supposed to.
population because you can’t test anything else until you have a correct population data frame.
is using Google Scholar. Search “scholar” on Google and it will give you Google Scholar. E.g.
searching “us household composition” led me to the data on which I based my
make_households function.
for testing and then easily plug in better values later after doing research.