The simplest ordinary differential equations that describe the spread of disease in a population are:
\begin{align*}
\dot{S} = -\beta I \frac{S}{N} \\
\dot{I} = \beta I \frac{S}{N} - \frac{I}{\tau} \\
\dot{R} = \frac{I}{\tau}
\end{align*}
where:
\begin{equation*}
N = S + I + R
\end{equation*}
This is the so called SIR model. The time varying variables are:
- \(S\): number of people susceptible to infection
- \(I\): number of people infected
- \(R\): number of people recovered (i.e. no longer susceptible)
- \(N\): total number of people
The constant parameters are:
- \(\beta\): the time rate of infectious transmissions
- \(\tau\): the mean infectious period
- \(R_0 = \beta\tau\): the basic reproduction number
The right hand side of the ordinary differential equations are evaluated with the following function eval_SIR_rhs.m:
function xdot = eval_SIR_rhs(t, x, c)
% Inputs:
% t : 1x1 scalar time
% x : 3x1 state vector [S; I; R]
% c : structure containing 2 constant parameters
% Outputs:
% xdot : 3x1 state derivative vector [S'; I', R']
S = x(1); % susceptible [person]
I = x(2); % infected [person]
R = x(3); % recovered [person]
beta = c.beta; % infectious contacts per time [person/day]
tau = c.tau; % mean infectious period [day]
N = S + I + R;
Sdot = -beta*I*S/N;
Idot = beta*I*S/N - I/tau;
Rdot = I/tau;
xdot = [Sdot; Idot; Rdot];
end
The ordinary differential equations are integrated (simulated) with the following script simulate_SIR_model.m:
% constant parameters
c.beta = 1/7; % infectious contacts per day [person/days]
c.tau = 14; % days infectious [day]
% intial population values
S0 = 39.56E6; % population of California [person]
I0 = 1000; % infected [person]
R0 = 200; % recovered [person]
x0 = [S0; I0; R0];
num_months = 8;
num_days_per_month = 30;
ts = 0:num_months*num_days_per_month; % units of days
% integrate the equations
f = @(t, x) eval_SIR_rhs(t, x, c);
[ts, xs] = ode45(f, ts, x0);
% plot the results
hold on
title(sprintf('R0 = %1.2f', c.beta*c.tau))
plot(ts, xs(:, 1)) % S(t)
plot(ts, xs(:, 2)) % I(t)
plot(ts, xs(:, 3)) % R(t)
xlabel('Time [days]')
ylabel('Number of people')
legend('S', 'I', 'R')
hold off