Graphing The SIR Model With Python

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