# SIR Epidemic Model

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