[SOLVED] matlab 4H Mathematical Biology Computer Lab 1

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 sucient 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
= eR0SI 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 Matlabs 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 eect 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 aect 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

1R0
e

+

2
+ 4(R0 1)

e

1R0

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 eect of partial immunity on the endemic steady state (remember 0
1.). Consider your observations in terms of the eect 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 eect 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].