COVID-19, Ordinary Differential Equations & Power of Numerical Analysis

A quarantined graduate student is no different from a normal one. While there is nothing else to do in the lockdown, you'd expect the procrastination levels to drop. But no! we always find new ways to procrastinate. Instead of working on things I should be working on, I decided to develop and implement a classic "S-I-R Epidemic Model[1]" in MATLAB® to simulate the COVID-19 outbreak and compare the model predictions with the data [2] from first the 20 days of the outbreak in Islamic Republic of Iran.

The classic S-I-R model of epidemics is,

\[ \frac{dS}{dt} = -\epsilon \frac{S}{N} I, \ \ \ \ \ \ \ \ \ \ -(1) \]                      

\[\frac{dI}{dt} = \epsilon \frac{S}{N} I -\rho I, \ \ \ \ \ \ \ \ \ \ -(2) \]                     

\[\frac{dR}{dt} = \rho I \ \ \ \ \ \ \ \ \ \ -(3) \]                   

where, S, I and R are the susceptible, infected and the recovered number of people in a population of size N. 

\(\epsilon \) is the effective contact rate,

\(\rho \) is the recovery rate, 

\(1/\rho \) is the mean period of  time during which an infected individual can pass it on.

Also, basic reproduction number, 

\(R_0 =  \epsilon/\rho \) is the average number of people infected from one infected person.

The model is based on a few simple assumptions, which are:

1. We are dealing with a closed population of size 'N'.
2. The rates of transmission and recovery are assumed to be constant.
3. As the population is closed, so no new births/deaths are considered.
4. Well mixed population where any infected individual has a probability of contacting any susceptible individual.

The governing equations (Equations. 1 through 3), which are ordinary differential equations governing the growth of each set of S, I or R are integrated in time using 4th order Runge-Kutta method [3]. The implementation in MATLAB® is, #include iostream

%-----------------------------------------------------------------------\
% S-I-R Model Implementation for Epidemic Spread in Matlab |
% Integrator for rate ODEs: 4th Order Runge-Kutta |
%-----------------------------------------------------------------------/
% Shujaut H. Bader
% Ph.D. Candidate, Aerospace Engineering Department, Iowa State University
% -----------------------------------------------------------------------
clc

clear

ti = 0.0; % starting time

tf = 120; % final time, let's say, days!

Ndays = 120; % total no. of days

dt = (tf-ti)/(Ndays-1); % step size for integrator (a day)

t = ti:dt:tf; % array for time (days)

N = 15000; % started with this sample size.. you can use a different one

%epsilon and rho are given for Iran:

epsilon =0.5961; % 0.00003974*N % effective contact rate, infected person comes in contact with \epsilon times S

rho = 0.08023309; % recovery rate, takes 1/\rho days for someone infected to pass on the infection

% fprime is a variable name given to store RHSs of the rate equations

fprime = @(t,vec)[-epsilon*(vec(1)/N)*vec(2); epsilon*(vec(1)/N)*vec(2) - rho*vec(2) ; rho*vec(2)];

% where 't' is time in days, and 'vec' is vector of RHSs, vec(1) = S, vec(2)= I , vec(3) = R

% intial values of new variables S, I, R

% initially S = N, all the population is susceptible

% I = 2

% Zero recovered cases initially

vec(:,1) = [N;2;0];

% 4th-Order Runge-Kutta Iterations

for i = 1:Ndays

  s1 = dt*fprime(t(i),vec(:,i));

  s2 = dt*fprime(t(i)+0.5*dt, vec(:,i)+0.5*s1);

  s3 = dt*fprime(t(i)+0.5*dt, vec(:,i)+0.5*s2);

  s4 = dt*fprime(t(i)+dt, vec(:,i)+s3);

  vec(:,i+1) = vec(:,i) + s1/6 + s2/3 + s3/3 + s4/6;

end


plot(t(1:i),vec(1,1:i),'k','LineWidth',2) % plotting S

hold on

plot(t(1:i),vec(2,1:i),'r','LineWidth',2) % plotting I

hold on

plot(t(1:i),vec(3,1:i),'g','LineWidth',2) % plotting R

% hold on

% iran = load("irandata.xy"); % should've a file by that name with data

% plot(iran(1:20,1),iran(1:20,2),'o-');

xlabel("days")

ylabel("No. of people")

legend({'Susceptible, S','Infected, I','Recovered, R', 'Data(Iran)'},'NumColumns',1, 'Location', 'northeast')

title({'S-I-R Epidemic model predictions compared', 'with data the for first 20 days of the COVID-19 outbreak in Iran'})

Below, the results are shown for Reproduction number, contact rate and recovery suitable to Iran's situation in the first two weeks of the outbreak [4]. 

S-I-R model predictions

The real curve in Iran has indeed shot up and not in agreement with what this model will predict after 20 days. That is because of the assumptions we have made. This is a very rudimentary way of formulating how an epidemic can spread. In some simpler cases, it will work well- where all the assumptions are met. You can copy the source code into your MATLAB® editor and run it with different values of contact and recovery rates, simulating different scenarios.

Cite as

Bader, Shujaut H., “COVID-19, Ordinary Differential Equations & Power of Numerical Analysis.” Backscatter, April 1, 2020, www.backscatterblog.blogspot.com/2020/04/covid-19-ordinary-differential.html