The SIR Model for Spread of Disease
The SIR Model for Spread of Disease
Arnoud Buzing
Wolfram Research, Inc.
The Differential Equation Model
Introduction
Introduction
The SIR (Susceptible-Infected-Recovered) model for the spread of infectious diseases is a very simple model of three linear differential equations. It examines how an infected population spreads a disease to a susceptible population, which transforms into a recovered population.
This model assumes no one dies from the disease and that the population is also isolated and stable (no new people arrive by birth or migration during the model simulation).
The independent & dependent variables
The independent & dependent variables
In this model, the independent variable is t (time, in days):
In[]:=
Clear[t]
The dependent variables are S (number of susceptible individuals at time t), I (number of infected individuals at time t) and R (number of recovered individuals at time t).
To simplify the differential equations, we will use the scaled versions of the S, I and R variables. The scaled variables are s, i and r. The scale runs from 0 (representing a zero population) to 1 (representing the full population):
In[]:=
Clear[s,i,r]
The susceptible population starts near 1, with almost everyone being susceptible. It drops over time. The infected population starts with a positive small value near 0. The recovered population starts at 0, with no one having recovered. It increases over time.
At any given time, , because the sum represents the full scaled population, which is 1.
s+i+r=1
Each infected individual is assumed to spread their disease to b new individuals each day. If b is high, the disease will spread faster:
In[]:=
Clear[b]
The fraction k represents the fraction of the infected population that recovers each day. For example, if the disease lasts 14 days, then k is (eventually) 1/14:
In[]:=
Clear[k]
The differential equations
The differential equations
The first equation describes how the susceptible population changes over time. Because we start with a susceptible population of 1 that decreases over time, we can assume the rate of change will be negative. Also, the rate of change is proportional to the number of contacts b (more contacts decreases the susceptible population faster). Additionally, the rate of change is proportionally related to the current susceptible and infected population (but not the recovered population):
In[]:=
equationS=s'[t]-bs[t]i[t]
Out[]=
′
s
The second equation describes how the infected population changes over time. It has two terms. The first term is the same as the right-hand side of the previous equation, except the sign is now positive. It represents the increase in infected population. The second term is negative because it represents the decrease of infected cases (people who recover):
In[]:=
equationI=i'[t]bs[t]i[t]-ki[t]
Out[]=
′
i
The final equation describes how the recovered population changes over time. The right-hand side is the same as the last term of the previous equation. It represents the transition of infected cases to recovered cases:
In[]:=
equationR=r'[t]ki[t]
Out[]=
′
r
Before solving the differential equations, we need to set the parameters b and k.
Let's assume that on average, each infected person spreads their disease to one person every two days:
In[]:=
b=0.5
Out[]=
0.5
Here we assume that the average recovery time is 14 days, which means that on average, 1/14th of the infected population recovers:
In[]:=
k=1/14.0
Out[]=
0.0714286
We also need initial values at for the susceptible, infected and recovered population. We need a small nonzero value for the initial infected population, so we will set it to 0.0001. We will set the initial value for the susceptible population to 0.9999 (everyone else, who is not infected). And we will set the initial recovered population to 0.
t=0
We pass the differential equations and the initial values to NDSolve and request a solution for days:
t=100
In[]:=
solution=NDSolve[{equationS,equationI,equationR,s[0]0.9999,i[0]0.0001,r[0]0.0000},{s,r,i},{t,100}];
Next we get the individual solutions for S, I and R:
In[]:=
solutionS=First[s/.solution];solutionI=First[i/.solution];solutionR=First[r/.solution];
Visualization
Visualization
We can plot the results. Here, the blue line represents the susceptible population over time, red is infected, and green is recovered:
In[]:=
Plot[{solutionS[t],solutionI[t],solutionR[t]},{t,0,100},PlotRange{0,1.01},PlotStyle{Blue,Red,Green},PlotLegends{"S(susceptible)","I(infected)","R(recovered"}]
Out[]=
Alternatively, we can stack these plots to highlight that their sum is always 1:
In[]:=
Plot[Accumulate[{solutionS[t],solutionI[t],solutionR[t]}]//Evaluate,{t,0,100},PlotRange{0,1.01},Filling{1Axis,2{1},3{2}},PlotStyle{Blue,Red,Green},PlotLegends{"S","S+I","S+I+R"}]
Out[]=
