SEIRモデルで新型コロナウイルスの挙動を予測してみた。

#1.はじめに
2019年12月に中国の湖北省の武漢で発生した新型コロナ肺炎(Covid-19)が、日本にも感染者数が増加しています。インフルエンザ, AIDS, SARS,などの感染病がどのプロセスで人間集団の中で拡大していくことを理解することは、予防接種の設定、感染者の隔離などの保健政策の効果を確認するためにも重要であります。
前回の記事では、SIRモデルを説明しましたが、感染病の潜伏期間を考慮したSEIRモデルを紹介します。
このモデルをPythonで計算する過程と、このSEIRモデルを利用し新型コロナウイルスがどう拡散するか計算結果を紹介します。

追記：
この記事の内容をベースにGUI化してみました。

新型コロナウイルスのGUIシミュレーション　(SEIRモデル）
https://qiita.com/kotai2003/items/f6cf36e9c22c3e776dee

#2.伝染病の予測モデル(SEIRモデル)

SEIRモデルでは、全体人口を次の集団に分類し、時間に関する各集団の人口増減を微分方程式で表す。　

• S(Susceptible) : 伝染病に対する免疫がない。感染可能者
• E(Exposed)：感染症が潜伏期間中の者
• I(Infectious)：感染により発病し、感染可能者(S)への伝染可能な者。感染者
• R (Removed) : 発病から回復し、免疫を獲得したもの。あるいは発病から回復できず死亡した者。（このモデルのシステムから排除されるため、Removedと言う。）

各集団の人口増減を次の微分方程式で表す。

\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 &: 感染しうるもの、免疫を持たず感染可能 \quad \text{(Susceptible)}　\\
E &: 感染症が潜伏期間中のもの \quad \text{(Infectious)}　\\
I &: 感染症が発症しているもの、接触した感染しうるもの(S)に病気を感染 \quad \text{(Infectious)}　\\
R &: 感染後死亡者、もしくは免疫を獲得した者 \quad 
\text{(Removed)}　\\
N &: 全人口, S+E+I+R
\end{align}

\begin{align}
\beta &: 感染率\quad \text{(The infectious rate)}  \\
\epsilon &: 暴露後に感染症を得る率\quad \text{(The rate at which an exposed person becomes infective)} \quad [1/day] \\
\gamma &:除去率\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

\begin{align}
\ l_p &: 感染待ち時間 \text{(latency period [day])}\quad \epsilon= \frac{1}{l_p} \\
\ i_p &: 感染期間 \text{(Infectious period [day])}\quad \gamma= \frac{1}{i_p} \\
\end{align}

このSIRモデルの前提条件

• 一度免疫を獲得した者は2度と感染することなく、免疫を失うこともありません。
• 全体人口で外部からの流入及び流出はありません。なお出生者及び感染以外の原因で死亡する者もいません。

#3.SEIRモデルの計算

SIRモデルは、数値積分を用いて解きます。ここでは、PythonのScipyモジュールのRunge-Kutta方程式を解くodeint関数を利用します。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

次に、SIRモデルの微分方程式をodeintで計算できる形に記述します。ここで、v[0], v[1], v[2],v[3]がそれぞれS,E,I,Rに対応します。

# 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]]

そして、数値積分に必要な各パラメータ及び初期条件を定義します。ここでは、人口1,000名の集団に対し、初期の感染症が潜伏期間中の者が1名がいると仮定します。
感染率を1に、感染待ち時間を2日、感染期間を7.4日と設定します。

# 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]]


#感染率
beta_const = 1 #感染率

#暴露後に感染症を得る率
latency_period = 2 #days
epsilon_const = 1/latency_period

#回復率や隔離率
infectious_period = 7.4 #days
gamma_const = 1/infectious_period

数値積分結果をresultに格納し、この結果を時間軸に対しプロットします。

# 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()

このグラフから言えることは、100名の集団に感染症が潜伏期間中の者が1名(Exposed(0))発生することで、最終的に100名(Removed(100))が感染を経験することである。そして、感染者(Infectious)のピークは約18日辺りで発生し、その後は感染が収束することがわかる。
このように感染病に関するパラメータをもって、感染病が人口集団に及ぼす影響を評価することが可能である。

#4.コロナウイルスの予測

現在、新型コロナウイルスの発症事例より、SEIRのパラメータを推定する研究論文が多数発表されています。今回は、2月16日に発表された論文に掲載されたパラメータ推定値で、SEIRモデルを計算してみます。（参考文献　2)

Parameter	中国本土（湖北省除く）	湖北省（武漢除く）	武漢
人口数 N (million)	1340	45	14
感染率　[beta]	1.0	1.0	1.0
Latency period (days)	2	2	2
infectious_period (days)	6.6	7.2	7.4
E_0	696	592	318
I_0	652	515	389

計算結果を示します。それぞれ縦軸の人口の数値が異なることに注意してください。
1番目のグラフは、湖北省を除く中国本土の感染拡大の予測です。

2番目のグラフは、武漢を除く湖北省の感染拡大の予測です。

3番目のグラフは、武漢の感染拡大の予測です。

このパラメータでの予測だと、30日～40日で感染者数はピークを迎えて、最終的にはほぼ100%の人口が感染を経験する（Removed(100)）ことになります。正直、驚きました。
恐らくWorst Caseを想定した予測結果になったと思います。今回のコロナウイルスには各段に注意が必要だと思います。

#5. 参考文献

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. 2015年にコロナウイルスの流行を予測した論文の紹介 https://qiita.com/kotai2003/private/a44ca921314d17cc62e3