4H Mathematical Biology Computer Lab 1
Epidemic models
An SIR model with partial immunity
The basic SIR model can be extended to model a disease for which immunity is only tem-
porary, in that individuals in the recovered class can be reinfected at a slower rate than the
susceptibles who have never been infected.
Consider the modified SIR model :
ds
dT
= µn� �si� µs
di
dT
= +�si� �i+ ��ir � µi
dr
dT
= �i� ��ir � µr
where all the parameters are positive constants and 0 � 1.
1. Explain the terms in the model. Let s(T ) + i(T ) + r(T ) = N(T ) and N(0)=n, show
that N(T ) = n for all time, T � 0. Why it is su�cient to consider only equations for
s and i?
2. Nondimensionalise the model by setting
s = Sn, i = In, T = t/(µ+ �),
where S, I, t are scaled variables, susceptibles, infectives and time respectively. Show
the non-dimensionalised model is given by
dS
dt
= e�R0SI � eS
dI
dt
= +R0SI � I + �R0I(1� S � I)
(1)
Give the expressions for R0 and e in terms of the original model parameters. Explain
why the parameter R0 is the basic reproduction number for this model.
3. Letting
e = 0.0012, R0 = 3.5, � = 0.25
and initial conditions S(0) = 0.99, I(0) = 0.01. Write a matlab program to solve
equations (1) numerically using Matlab’s ODE solver ODE45. You should be creating
two files: a function file which contains the equations and a script file that calls the
ODE solver which will access your equation file. Plot your solution, I(t) against t, S(t)
against t and R(t) (= 1 � S(t) � I(t)) against t (on the same graph) for 0 t 50.
Now vary R0. What is the e↵ect of R0 on the graphs, pay particular attention to
features such as peak infection, whether the infection dies out, and the endemic steady
state. Why does R0 a↵ect the graph in this way?
4H Mathematical Biology Computer Lab 1
4. For R0 > 0 the number of infectives in the unique endemic steady state (I 6= 0) is
given by
I =
e
2R0�
2
4
s✓
1�R0�
e
+ �
◆2
+ 4(R0 � 1)
�
e
�
✓
1�R0�
e
+ �
◆3
5
Write a matlab program to plot the steady state value of I given in the equation above
versus R0, for 1 R0 6. Set � = 0.25 and e = 0.0012.
5. Try changing the value of �, how does this change the graph, try some other values
and study the e↵ect of partial immunity on the endemic steady state (remember 0
� 1.). Consider your observations in terms of the e↵ect of partial immunity on the
terms in the model.
6. Now we modify the model so that a fraction v of the population are vaccinated at birth.
Introduce a new compartment for individuals who are vaccinated and are immune to
the disease. Modify your equations so that a proportion v go directly into this class and
the rest become susceptible (vaccinated individuals will still die with rate µ). Setting
R0 = 5 and using the same parameters as in part (c), explore the e↵ect of v on the
dynamics. In particular think about when infection dies out and when it persists. You
might find it helpful to pass a parameter into the ODE solver using the syntax:
% function
function dy=myfunction(t,y,R0)
…
end
%script
[T,Y] = ode45(@(t,y)myfunction(t,y,R0),[0:1:500],[0.99 0.01]);
then run the solver within a loop.
Some useful matlab commands
hold on If you use this command directly following a plot it will allow any subsequent plots
to appear on the same figure (it will not overwrite them).
size size(Y) gives the size of the matrix Y . For example if Y is an 7 ⇥ 4 matrix then
size(Y)=[7 4].