[PYTHON] I tried to predict the behavior of the new coronavirus with the SEIR model.

1.First of all

The number of people infected with the new coronavirus pneumonia (Covid-19) that occurred in Wuhan, Hubei Province, China in December 2019 is increasing in Japan as well. Understanding the process by which infectious diseases such as influenza, AIDS, SARS, etc. spread within the human population is to confirm the effects of health policies such as vaccination settings and isolation of infected individuals. Is also important. In the previous article, I explained the SIR model, but I will introduce the SEIR model that considers the incubation period of infectious diseases. I will introduce the process of calculating this model with Python and the calculation result of how the new coronavirus spreads using this SEIR model.

Addendum: I tried to make it GUI based on the contents of this article.

** GUI simulation of new coronavirus (SEIR model) ** https://qiita.com/kotai2003/items/f6cf36e9c22c3e776dee

2. Infectious disease prediction model (SEIR model)

In the SEIR 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 --E (Exposed): Those who have an infectious disease during the incubation period --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 \frac{SI}{N} \\
\frac{dE}{dt} &=  \beta \frac{SI}{N}  -\epsilon E \\
\frac{dI}{dt} &=  \epsilon E -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

\end{align} 

\begin{align}
S &:What can be infected, can be infected without immunity\quad \text{(Susceptible)}　\\
E &:Infectious disease during the incubation period\quad \text{(Infectious)}　\\
I &:Those who have an infectious disease, those who can be infected by contact(S)Infects the disease\quad \text{(Infectious)}　\\
R &:Those who died after infection or who acquired immunity\quad 
\text{(Removed)}　\\
N &:Total population, S+E+I+R
\end{align}

\begin{align}
\beta &:Infection rate\quad \text{(The infectious rate)}  \\
\epsilon &:Rate of getting infection after exposure\quad \text{(The rate at which an exposed person becomes infective)} \quad [1/day] \\
\gamma &:Exclusion rate\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

\begin{align}
\ l_p &:Infection waiting time\text{(latency period [day])}\quad \epsilon= \frac{1}{l_p} \\
\ i_p &:Infection period\text{(Infectious period [day])}\quad \gamma= \frac{1}{i_p} \\
\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 from any cause other than birth and infection.

3. SEIR 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
import matplotlib.pyplot as plt

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

# Define differential equation of SEIR model

'''
dS/dt = -beta * S * I / N
dE/dt = beta* S * I / N - epsilon * E
dI/dt = epsilon * E - gamma * I
dR/dt = gamma * I

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

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

'''


def SEIR_EQ(v, t, beta, epsilon, gamma, N ):
    return [-beta * v[0] * v[2] / N ,beta * v[0] * v[2] / N - epsilon * v[1],
            epsilon * v[1] - gamma * v[2],gamma * v[2]]

Then, define each parameter and initial conditions required for numerical integration. Let's assume that for a population of 1,000, there is one person with an early infection during the incubation period. Set the infection rate to 1, the infection waiting time to 2 days, and the infection period to 7.4 days.

# parameters
t_max = 100 #days
dt = 0.01

# initial_state
S_0 = 99
E_0 = 1
I_0 = 0
R_0 = 0
N_pop = S_0 + E_0 + I_0 + R_0
ini_state = [S_0, E_0, I_0, R_0]  # [S[0],E,[0], I[0], R[0]]


#Infection rate
beta_const = 1 #Infection rate

#Rate of getting infection after exposure
latency_period = 2 #days
epsilon_const = 1/latency_period

#Recovery rate and quarantine rate
infectious_period = 7.4 #days
gamma_const = 1/infectious_period

Store the numerical integration result in result and plot this result on the time axis.

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

# Numerical Solution using scipy.integrate
# Solver SEIR model
result = odeint(SEIR_EQ, ini_state, times, args)
# plot
plt.plot(times, result)
plt.legend(['Susceptible', 'Exposed', 'Infectious', 'Removed'])
plt.title("SEIR model  COVID-19")
plt.xlabel('time(days)')
plt.ylabel('population')
plt.grid()

plt.show()

From this graph, it can be said that in a population of 100 people, one person (Exposed (0)) has an infection during the incubation period, and finally 100 people (Removed (100)) experience the infection. That is. It can be seen that the peak of infected people (Infectious) occurs around 18 days, after which the infection converges. In this way, it is possible to evaluate the effects of infectious diseases on the population with parameters related to infectious diseases.

4. Coronavirus prediction

Currently, many research papers have been published to estimate SEIR parameters from the cases of new coronavirus cases. This time, I will calculate the SEIR model with the parameter estimates published in the paper published on February 16. (Reference 2)

The calculation result is shown. Note that each has a different population figure on the vertical axis. The first graph is a forecast of the spread of infection in mainland China, excluding Hubei.

The second graph is a forecast of the spread of infection in Hubei Province, excluding Wuhan.

The third graph is a prediction of the spread of Wuhan infection.

With this parameter, the number of infected people peaks in 30 to 40 days, and eventually almost 100% of the population will experience infection (Removed (100)). I was really surprized. Probably, the prediction result is based on the worst case. I think it is necessary to pay attention to each step in this coronavirus.

5. References

1. SEIR and SEIR models https://institutefordiseasemodeling.github.io/Documentation/general/model-seir.html

2. Epidemic analysis of COVID-19 in China by dynamical modeling https://arxiv.org/abs/2002.06563

3. Introduction of papers that predicted the epidemic of coronavirus in 2015 https://qiita.com/kotai2003/private/a44ca921314d17cc62e3