[PYTHON] Introduction of mathematical prediction model for infectious diseases (SIR model)

1. Introduction

The number of people infected with the new coronavirus (Convid-19) that occurred in Wuhan, China in December 2019 is increasing in Japan. Understanding how infectious diseases such as influenza, AIDS, SARS, etc. spread in human populations is to confirm the effects of health policies such as vaccination settings and isolation of infected individuals. Is also important. Here, I will introduce the SIR model, which is the basis of the infectious disease mathematical prediction model, and introduce the process of calculating this model in Python.

Addendum: Although it has nothing to do with programming, a paper predicting the epidemic of the new coronavirus was published in the 2015 Nature poem. If you are interested, please read the article linked below as it will give you a hint as to how the situation of this new coronavirus will converge.

** Introduction of papers predicting the epidemic of new coronavirus in 2015 ** https://qiita.com/kotai2003/private/a44ca921314d17cc62e3

2. Mathematical prediction model of infectious diseases (SIR model)

The SIR model is a basic model that describes the epidemic process of infectious diseases. The name of the model comes from the acronym Susceptible, Infectious, Recovered for the variable names in the model. First introduced in a paper by W.O. Kermack and A.G. McKendrick in 1927.

In the SIR model, the total population is classified into the following groups, and the population increase / decrease of each group with respect to time is expressed by a differential equation. --S (Susceptible): No immunity to infectious diseases. Infectable person --I Infectious): A person who becomes ill due to infection and can be transmitted to a person who can be infected (S). Infected person --R (Removed): Those who have recovered from the onset of illness and acquired immunity. Or those who died without being able to recover from the onset of illness. (It is called Removed because it is excluded from the system of this model.)

The population increase / decrease of each group is expressed by the following differential equation.

\begin{align}

\frac{dS}{dt} &= -\beta SI \\
\frac{dI}{dt} &=  \beta SI -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

\end{align} 

\begin{align}
S &:Infectable person\quad \text{(Susceptible)}　\\
I &:Infected person\quad \text{(Infectious)}　\\
R &:Those who died after infection or who acquired immunity\quad \text{(Removed)}　\\
\beta &:Infection rate\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:Exclusion rate\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

Prerequisites for this SIR model --A person who has acquired immunity will never be infected and will not lose immunity. --There is no inflow or outflow from the outside in the total population. No one died due to causes other than birth and infection.

3. SIR model calculation

The SIR model is solved using numerical integration. Here, we use the odeint function that solves the Runge-Kutta equation of Python's Scipy module.

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
%matplotlib inline

Next, write the differential equation of the SIR model in a form that can be calculated by odeint. The code v [0], v [1], v [2] correspond to S, I, R respectively.

#Define differential equation of SIR model

'''
dS/dt = -beta * S * I
dI/dt = beta * S * I - gamma * I
dR/dt = gamma * I

[v[0], v[1], v[2]]=[S, I, R]

dv[0]/dt = -beta * v[0] * v[1]
dv[1]/dt = beta * v[0] * v[1] - gamma * v[1]
dv[2]/dt = gamma * v[1]

'''
def SIR_EQ(v, t, beta, gamma):
    return [-beta*v[0]*v[1], beta * v[0] * v[1] - gamma * v[1], gamma * v[1]]

Then, define each parameter and initial conditions required for numerical integration. Here we assume that for a population of 1,000, there was an initial infection. Set the infection rate to 0.2 / 1000 and the removal rate to 0.1.

\begin{align}
\beta &=0.2/1000 \\
\gamma &=0.1 \\
\end{align} 

Finally, calculate R0 of Basic Reproduction Number. If R0> 1, this infectious disease does not converge, but spreads.

#parameters
t_max = 160
dt = 0.01

beta_const = 0.2/1000
gamma_const = 0.1


#initial_state
S_0=999
I_0=1
R_0=0
ini_state = [S_0,I_0,R_0] #[S[0], I[0], R[0]]

#numerical integration
times =np.arange(0,t_max, dt)
args  =(beta_const, gamma_const)

#R0
N_total = S_0+I_0+R_0
R0 = N_total*beta_const *(1/gamma_const)
print(R0)

The numerical integration result is stored in result, and this result is plotted on the time axis.

#Numerical Solution using scipy.integrate
#Solver SIR model
result = odeint(SIR_EQ, ini_state, times, args)
#plot
plt.plot(times,result)
plt.legend(['Susceptible','Infectious', 'Removed'])

What can be said from this graph is that if one infected person (Infectious (0)) occurs in a population of 1,000 people, about 800 people (Removed (160)) will eventually experience the infection. And it can be seen that the peak of infected people (Infectious) occurs around 70 days, and then the infection converges. In this way, it is possible to evaluate the effects of infectious diseases on the population by using the parameters related to infectious diseases.

4. Reference materials

-Coronavirus COVID-19 Global Cases https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html?fbclid=IwAR2FwOf4Cjm4okosqFt0Ddr4k0dUgswM28oAqYkkVY6QT6BBCZQ1NlfPDXk#/bda7594740fd40299423467b48e9ecf6

https://qiita.com/Student-M/items/4e3e286bf08b7320b665

https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/

https://www.math.unm.edu/~sulsky/mathcamp/ApplyData.pdf