[PYTHON] Mathematical prediction model for infectious diseases (SIR model): Case study (1)

1. Introduction

In the previous article, I introduced the SIR model, which is a mathematical prediction model for infectious diseases.

Title: Introduction of Mathematical Prediction Model for Infectious Diseases (SIR Model) https://qiita.com/kotai2003/items/3078f4095c3e94e5325c

This time, we will list the cases of infection that actually occurred and confirm the predictive performance of this SIR model.

2. Influenza infection at a boarding school in the United Kingdom (1978)

In January 1978, an influenza outbreak occurred at a boarding school in the north of England. Of the 763 boys aged 10 to 18 years old, all but 30 were confirmed to be infected with influenza. Of these, 512 were ill during the epidemic that lasted from 22 January to 4 February. Also, the infection started with one boy.

The table shows the value of the number of people with medical conditions (Infectious infected) estimated from the figure.

Plot these data in the figure. The number of Infectious people peaked on day 6 with a value of 255.

3. Parameter estimation

Now reconfirm the differential equations of the SIR model. In order to solve this model, we need to know the initial conditions of each variable. This infection started with one of the 763 school students, so S0 = 762, I0, R0 = 0.

\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}

The problem is infection rate and removal rate. In general, the SIR model is needed early in the outbreak of infection. Quickly estimate infection and elimination parameters from early infected data and calculate SIR models to see if this infection spreads or converges. If the SIR model shows that the infection spreads (R0> 1), it goes into the process of prioritizing health policies such as quarantining infected people and restricting the movement of people to prevent this spread.

Early infections showed that one sick boy infected two more boys a day later. As a rough estimate, the number of people who can be infected is expected to decrease at a rate of 2 per day.

\frac{dS}{dt} \simeq -2 \\

Since S_0 = 762 and I_0 = 1 are known as the initial conditions, the infection rate is calculated by the following formula.

\beta = 
\frac{-dS/dt}{S_0I_0}\simeq\frac{2}{762\times  1}=0.0026

Next is the estimation of the removal rate. The removal rate is the rate at which infected individuals decrease from the number of infected individuals due to death and immunity. According to reports at the time, the boys who were hospitalized were found to recover within a day or two. Therefore, assuming recovery in two days, about half of the infected population will be removed daily.

\gamma = \frac{1}{2\text{days}}=0.5[\frac{1}{day} ]

For reference, calculate the ratio of removal rate to infection rate. Its value is 192.

\rho = \frac{\gamma}{\beta}=\frac{0.5}{0.0026}= 192

4. Comparison of model prediction and actual measurement data

Put the above conditions and solve the SIR model by numerical integration.

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
%matplotlib inline
#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]]

Enter the initial conditions S_0, I_0, R_0 and the infection rate and removal rate in this part.

#parameters
t_max = 14
dt = 0.01

beta_const = 0.0026
gamma_const = 0.5


#initial_state
S_0=762
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(1,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)

#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'])
plt.title("Influenza 1978: a=0.0026, b=0.5, $S_0$=762, $I_0$=1")
plt.xlabel('time(days)')
plt.ylabel('population')

plt.show()

Plot the actual infected person data with the SIR model.

# Real Data
data_day = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
data_infectious = [3,8,28,75,221,291,255,235,190,125,70,28,12,5]

#SIR model Plot
plt.plot(times,result)

# Real Data Plot

plt.scatter(data_day, data_infectious,s=10, c="pink", alpha=0.5, linewidths="2",
            edgecolors="red")
plt.title("Influenza 1978: a=0.0026, b=0.5, $S_0$=762, $I_0$=1")
plt.legend(['Susceptible','Infectious', 'Removed','data'])
plt.xlabel('time(days)')
plt.ylabel('population')
plt.show()

Comparing the SIR model Infectious (infected person) in the figure with the actual infected person's data (data), it can be confirmed that the trends generally match. The maximum number of infected people (Infections_Max) predicted by this SIR model was 306, and the actual maximum number of infected people was 291 people, with an error of 5%. And the last infectable person (= number of people who were not infected until the end) predicted by the SIR model was 16. In reality, there were 30 people.

5. Summary

This time, we confirmed the effectiveness of the SIR model using actual cases of infectious diseases. From the information of infected persons at the initial stage, the initial conditions and the parameters of infection rate and removal rate were estimated. Then, by comparing with the actual infected person data, it was confirmed that it was almost in agreement with the prediction result of the SIR model.

6. Reference materials

1. https://qiita.com/kotai2003/items/3078f4095c3e94e5325c
2. https://www.researchgate.net/publication/336701551_On_parameter_estimation_approaches_for_predicting_disease_transmission_through_optimization_deep_learning_and_statistical_inference_methods
3. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/pdf/brmedj00115-0064.pdf