【Python】SEIRモデルを実装してみた

はじめに

新型コロナウイルス感染症（COVID-19）の感染者数がどのように予測されているのか気になって調べてみたところ、SEIRモデルと呼ばれる感染症の流行を予測する数理モデルがあるという情報を得ました。 「数理モデルなら簡単に実装できるのでは…？」という軽いノリで実装してみたので、その覚え書きを残しておきます。

ja.wikipedia.org

※私自身は感染症に関して全くの素人です。 記事中の内容に関して間違い等がございましたら、ご指摘いただけるとありがたいです。 なお、本記事によって生じた損害については、一切責任を負いません。 あくまでも個人的な趣味の範囲内でご利用ください。

参考文献

実装にあたっては、東京大学大学院理学系研究科生物科学専攻のサイトで公開されていた以下のページを参考にしました。

「新型コロナウイルス感染症の流行予測　正しく理解し、正しく怖がり、適切な行動をとるために」（最終閲覧日：2020-04-15）

Pythonで実装

Susceptible、Exposed、Infectious、Recoveredに関する微分方程式をそれぞれ離散化して実装します。 時間方向の微分だけなので、簡単にできました。 取り敢えず、人口100万人の都市に1人の感染者がやってきたという想定で計算してみます。

# --- ライブラリの読み込み
import matplotlib.pyplot as plt
import numpy as np

# --- パラメーターを指定
n       = 1000000 # --- 人口
days    = 365     # --- 予測日数
delta_t = 1.0     # --- タイムステップ [日]
r0      = 2.5     # --- 基本再生産数
l       = 5.0     # --- 平均潜伏期間
i       = 10.0    # --- 平均発症期間

# --- 予測値を格納する配列を作成
Sus = np.full(days, np.nan) # --- 未感染者
Exp = np.full(days, np.nan) # --- 潜伏期感染者
Inf = np.full(days, np.nan) # --- 発症者
Rec = np.full(days, np.nan) # --- 免疫獲得者

# --- 初期値を設定
Sus[0] = n
Exp[0] = 0.0
Inf[0] = 1.0
Rec[0] = 0.0

# --- 翌日の値を計算
for t in range(days-1):
    Sus[t+1] = Sus[t] - ((r0 / i) * (Sus[t] / n) * Inf[t]) * delta_t
    Exp[t+1] = Exp[t] + ((r0 / i) * (Sus[t] / n) * Inf[t] - (1.0 / l) * Exp[t]) * delta_t
    Inf[t+1] = Inf[t] + ((1.0 / l) * Exp[t] - (1.0 / i) * Inf[t]) * delta_t
    Rec[t+1] = Rec[t] + ((1.0 / i) * Inf[t]) * delta_t

# --- グラフを作成
fig = plt.figure(figsize = (8, 4))
ax = fig.add_subplot(111)

ax.plot(np.arange(days), Sus, linewidth = 2, label = 'Susceptible')
ax.plot(np.arange(days), Exp, linewidth = 2, label = 'Exposed')
ax.plot(np.arange(days), Inf, linewidth = 2, label = 'Infectious')
ax.plot(np.arange(days), Rec, linewidth = 2, label = 'Recovered')

ax.set_xlabel('Lapsed Days', fontsize = 16)
ax.set_ylabel('Number of People', fontsize = 16)
ax.tick_params(labelsize = 14)
plt.tight_layout()
plt.legend(fontsize = 14)

plt.savefig('./seir.png')

出力結果

こんな感じの図が描けました。 モデルのパラメーターを変えるとそれぞれの予測時系列がどのように変化するのか、色々実験してみると面白そうです。