Graphing The SIR Model With Python
Last Updated on July 21, 2023 by Editorial Team
Author(s): Joaquin de Castro
Originally published on Towards AI.
Data Visualization, Programming
Graphing and solving simultaneous differential equations to model COVID-19 spread
If one good thing has come out of the COVID-19 pandemic, itβs the vast amount of data we have acquired. In light of technological advancements, we have access to more information and computing power, which we can use to predict and curb the spread of the virus. One of the simplest ways to do this is through the SIR model.
The SIR is a compartmental model that categorizes a constant population into three groups, namely the susceptible, infected, and recovered. These can all be expressed as functions that take time as an argument.
S(t) β The number of people who are susceptible to the disease
I(t) β The number of people infected with the disease
R(t) β Those incapable of infecting others; either recovered or diseased (hence a better name for this group may be βremovedβ)
Of note, S(t) + I(t) + R(t) = N at any time t, where N is the total population.
Importing the needed Python libraries
from scipy.integrate import odeintimport numpy as npimport matplotlib.pyplot as plt
There are several libraries that we can use to smoothen out the calculations and graphing process. The SciPy library contains calculation methods that are βuser-friendly and efficient.β We will be using their numerical integration function for this program. We will also be using Numpy for float step values and Matplotlib for graphing.
Finding the rate of change of each function
We canβt directly find an equation for each function. But we can derive the rate of change of each function at time t. In other words, it's a derivative. The amount of susceptible people generally starts off close to the total population. This number then decreases with time as susceptible people become infected. The number of newly infected people is a percentage of the possible interactions between susceptible and infected individuals. We can call this infection rate as βaβ, while the possible interactions being the product of S(t) and I(t). The change in susceptible people is therefore
Sβ(t) = -a*S(t)*I(t)
The decrease in the number of susceptible people is the same as the increase in the number of infected people. To find the entire derivative of I(t), we must also consider those who have recovered or died after being infected. That is simply the recovery rate multiplied by the current number of infected individuals. With the recovery rate as βbβ, we then have
Iβ(t) = a*S(t)*I(t) β b*I(t)
Calculating the derivative of R(t) is then a simple matter, as it is just the second term of Iβ(t). In the SIR model, recovery (or more aptly removal) only increases with time. The increase in R(t) is the product of the recovery rate and the infected population:
Rβ(t) = b*I(t)
We can now use these derivatives to solve the system of ordinary differential equations through the SciPy library.
Defining the necessary constants
From our equations above, we see that there are two constants that need to be defined: the transmission rate and the recovery rate. For now, weβll set the transmission rate to be 100% and the recovery rate to be 10%.
a = 1 # infection rate
b = 0.1 # recovery rate
Creating function f(y, t) to calculate derivatives
# FUNCTION TO RETURN DERIVATIVES AT Tdef f(y,t):
S, I, R = y d0 = -a*S*I # derivative of S(t) d1 = a*S*I β b*I # derivative of I(t) d2 = b*I # derivative of R(t) return [d0, d1, d2]
Next, we have to define a function to return the derivatives of S(t), I(t), and R(t) for a given time t. Remember that we actually solved for these already and are just encoding the equations in a function. In the following lines of code, d0, d1, and d2 are the derivatives of S, I, and R, respectively.
d0 = -a*S*I # derivative of S(t)d1 = a*S*I β b*I # derivative of I(t)d2 = b*I # derivative of R(t)
Note that the values of S, I, and R are not yet defined here (although we donβt need R(t) to find its derivative). Since these functions are dependent on each other, we will first get the previous values of S(t), I(t), and R(t) to calculate their derivatives.
S, I, R = y# orS = y[0]
I = y[1]
R = y[2]
Replacing the variable R with an underscore _ is also acceptable, but itβs best to be explicit and descriptive.
Defining the necessary initial values
Before we can calculate the values of the functions at time t, we must first find the initial values. The Philippines will be used as a sample population, t=0 on March 1, 2020. The initial number of susceptible people is the total population, which is around 109,581,078. Based on the Department of Healthβs COVID-19 tracker, the initial cases on March 1, 2020, is 50 people. And of course, we can set the total recovered from being 0. It would make things cleaner to keep the values between 0 and 1. Such can be accomplished by dividing all the values by the population.
S_0 = 1
I_0 = 50/109_581_078
R_0 = 0
It would be helpful to store these values in a list, and weβll see why this is important later.
y_0 = [S_0,I_0,R_0]
Letβs also not forget to define the domain or the time range. Weβll do this with Numpyβs linspace so as to incorporate decimal values for steps if desired.
t = np.linspace(start=1,stop=100,num=100)
Solving the differential equations with ODEInt
Now that we have defined all the needed variables, constants, and functions, we can now solve the system of ordinary differential equations. We will be solving each differential equation for a range of 100 days as specified in the time range. Doing so is actually simple with all the parameters set; the code is just one line:
y = odeint(f,y_0,t)
f is the function f(y, t), which we defined earlier, y_0 is the list of our initial values, and t is the list of equally spaced time values.
The variable y is actually a Numpy ndarray, or an N-dimensional array.
>>> print(type(y))
<class βnumpy.ndarrayβ>
In this case the odeint() function returns a 2 dimensional array. A row is a list of values for a given time t, while a column represents the values for S, I, or R. The values for S(t) would be found in the first column, I(t) in the second, and R(t) in the third. We can, therefore, access them as follows:
An element of such an array located in the nth row, mth column can be indexed as
element = y[n,m]
And we can use list splicing to view all the values in a column.
S = y[:,0]
I = y[:,1]
R = y[:,2]
Graphing the values for each function
plt.plot(t,S,βrβ,label='S(t)')plt.plot(t,I,βbβ,label='I(t)')plt.plot(t,R,βgβ,label='R(t)')plt.legend()plt.show()
The first 2 arguments in the plot() function indicate the x and y values. We can then specify the color and label of each line. The legend() function places a legend on the graph, referencing the βlabelβ keyword argument for each line.
Finally, we use the show() function to actually generate the graph.
Further exploration
It is worth pointing out that we can increase the precision of the graph by decreasing the specified step values. (In this case, by increasing the βnumβ keyword argument).
t = np.linspace(start=1,stop=100,num=200)
t = np.linspace(start=1,stop=100,num=500)
t = np.linspace(start=1,stop=100,num=1000)
Moreover, making the parameters such as the transmission and recovery rate closer to the actual data would make the model more accurate.
All that said and done, the SIR model demonstrates the invaluable role technology and mathematics plays in dealing with real-world issues.
You can find the code for this article here:
https://github.com/JoaquindeCastro/SIR-Model/blob/main/SIR-ODE-Integrate.py
References and resources
The SIR Model for Spread of Disease – The Differential Equation Model
As the first step in the modeling process, we identify the independent and dependent variables. The independentβ¦
www.maa.org
Pyplot tutorial – Matplotlib 3.3.2 documentation
An introduction to the pyplot interface. is a collection of functions that make matplotlib work like MATLAB. Eachβ¦
matplotlib.org
The N-dimensional array (ndarray) – NumPy v1.19 Manual
An is a (usually fixed-size) multidimensional container of items of the same type and size. The number of dimensionsβ¦
numpy.org
COVID-19 Tracker U+007C Department of Health website
Edit description
www.doh.gov.ph
A SIR model assumption for the spread of COVID-19 in different communities
In this paper, we study the effectiveness of the modeling approach on the pandemic due to the spreading of the novelβ¦
www.sciencedirect.com
Join thousands of data leaders on the AI newsletter. Join over 80,000 subscribers and keep up to date with the latest developments in AI. From research to projects and ideas. If you are building an AI startup, an AI-related product, or a service, we invite you to consider becoming aΒ sponsor.
Published via Towards AI