SIRDモデルによるウイルス感染症の伝播シミュレーション

作成: 日大生産工 MA 大前

  • SIRDモデルとは、ウイルス感染症のシミュレーションを行う数理モデルのことを指します。
  • 厳密なモデルではなく、SIRDモデルを簡単にしたものを扱います。簡易版であり、厳密なものではありません。
  • ここでは、pythonによりSIRDモデルを少しだけ体験し、卒業研究として仕上げることを目標とします。

これは、大前研究室の卒業研究テーマの学習資料として作成されたものです。他のゼミ生がこのページを参考にして卒業研究を行う場合は、必ず事前に相談してください。異なる研究室なのに、同じテーマが卒業研究概要集に並んでいると、問題になるかもしれません。

1. 数理モデル

SIRDモデルでは、

  • $S(t)$: $t$日目の感受性者数（いわゆる、感染したことがない健康な人たちのこと）
  • $I(t)$: $t$日目の感染者数（今まさしく感染している人の人数）
  • $R(t)$: $t$日目の回復者数（感染から回復した人数。免疫獲得により、今後は感染しない）
  • $D(t)$: $t$日目の死亡者数（感染後、回復せずに死亡した人数）

という4つの変数が定義されます。また、シミュレーションを行う集団の総人口を$N$としたとき、すべての$t$について、

\begin{eqnarray} N = S(t) + I(t) + R(t) + D(t) \end{eqnarray}

を成立させます。これは、考えるすべての人々が、必ず「感受性、感染、回復、死亡」のどれかになることを意味します。すなわち、シミュレーションの途中で、突然人が増えたり減ったりしないことを約束するものです。

ここで、$t+1$日目の感受性者数を、

\begin{eqnarray} S(t+1) = S(t) - \beta S(t)I(t) \end{eqnarray}

とします。$t+1$日目の感受性者数$S(t+1)$は、$t$日目の感受性者数$S(t)$から、$\beta S(t)I(t)$人だけ減少することを表現する数式となります。「$\beta S(t)I(t)$人」とは、「$t$日目から$t+1日目にかけて、$感染してしまった人」を表します。$S(t)$と$I(t)$の積が入っていますから、「感受性者（感染する可能性のある人）」と「感染者（感染している人）」が多ければ多いほど、次の日の感染者数が多くなる、つまり、感受性者数が減ることを意味しています（ある集団で、感染する・させる人が多ければ多いほど、感染者が増えるという構造です）。$\beta$はウイルスの感染力を表しており、これが大きいほど感染力が強いこと、小さいほど感染力が弱いことを表します。なお、ウイルス本来の感染力のみではなく、マスクをする、密集を回避するなどの対策を考慮した上での値となります。

そして、$\beta S(t)I(t)$が新たに感染する「新規感染者数」を意味するわけですから、$t+1$日目の感染者数は、

\begin{eqnarray} I(t+1) = I(t) + \beta S(t)I(t) - \gamma I(t) \end{eqnarray}

となります。$t$日目の感染者$I(t)$に対し、新規感染者数$\beta S(t)I(t)$が加えられているという構造を意味します。一方で、第3項目に$\gamma I(t)$が引き算されています。これは、感染後、回復または死亡することにより、感染者が減ることを意味しています。$\gamma$は感染者数を減らす、つまり除去するという意味合いで、除去率といいます。ややこしくなるので詳細な説明は省きますが、平均感染期間の逆数を除去率$\gamma$として採用することが普通です。例えば、ウイルスの平均感染期間を10日間にしたい場合は、$\gamma=1/10=0.1$を設定します。

続いて、$t+1$日目の回復者数$R(t+1)$と死亡者数$D(t+1)$は、

\begin{eqnarray} R(t+1) &=& R(t) + \gamma (1-\delta) I(t) \\ D(t+1) &=& D(t) + \gamma \delta I(t) \end{eqnarray}

となります。除去される感染者数$\gamma I(t)$のうち、

  • 割合$1-\delta$が回復、
  • 割合$\delta$が死亡、

するということを意味しています。すなわち、$\delta$は致死率を表します。致死率が高ければ死亡しやすくなり（回復しにくい）、低ければ死亡しにくくなります（回復しやすい）。

ここで、ギリシャ文字のパラメータを一度まとめておきますと、

  • $\beta$: ウイルスの感染力
  • $\gamma$: 除去率（平均感染期間の逆数）
  • $\delta$: 致死率

となります。この値と、各状態の初期人数、

  • $S(0)$: シミュレーション開始時点の感受性者数
  • $I(0)$: シミュレーション開始時点の感染者数
  • $R(0)$: シミュレーション開始時点の回復者数
  • $D(0)$: シミュレーション開始時点の死亡者数

を自分で決定すれば、シミュレーションが開始します。

2. シミュレーション

  • 今回は、合計1万人、初期感染者数が10人、感受性者数は9990人とします。最初は回復者や死亡者が居ないものとします。
  • また、感染力は0.00003、除去率は0.1（平均感染期間10日間の逆数）、致死率0.01（1%）とします。
  • この条件で、SIRDモデルの数式を実装し、計算します。
In [1]:
Tmax = 100 # シミュレーション日数

# 各状態の初期値をセット（N = 10000人でシミュレーション）
S, I, R, D = [], [], [], [] # 各状態を保存するリストを用意
S.append(9990) # 感受性者は9990人
I.append(10) # 感染者は10人
R.append(0) # 回復者はいない
D.append(0) # 死亡者はいない

# パラメータをセット
beta = 0.00003 # 感染力
gam = 0.1 # 除去率（平均感染期間を10日間とし、1/10=0.1を採用）
delta = 0.01 # 致死率（1%=0.01を採用）

# シミュレーション開始
for t in range(Tmax):
    S.append( S[t] - beta*S[t]*I[t] )
    I.append( I[t] + beta*S[t]*I[t] - gam*I[t] )
    R.append( R[t] + gam*(1-delta)*I[t] )
    D.append( D[t] + gam*delta*I[t])

3. 結果の可視化

得られた結果を可視化してみましょう。可視化には、plt.plot関数を利用します。グラフ中の英語はそれぞれ、

  • Susceptible: 状態S、感受性者
  • Infected: 状態I、感染者
  • Recovered: 状態R、回復者
  • Dead: 状態D、死亡者
  • Elapsed day: 経過日数
  • Number of persons: 人数

を表します。また、pltの引数は、

  • ls: 線の形状
  • lw: 線の太さ
  • c: 線の色

を意味します。

In [2]:
import matplotlib.pyplot as plt

#plt.figure(figsize=(20, 5)) 
plt.plot(S, label="Susceptible", lw=2, c="blue") # S: 感受性者数
plt.plot(I, label="Infected", ls="--", lw=2, c="red") # I: 感染者数
plt.plot(R, label="Recovered", ls="-.", lw=2, c="green") # 回復者数
plt.plot(D, label="Dead", ls=":", lw=3, c="k") # 死亡者数
plt.legend()
plt.grid()
plt.xlabel("Elapsed day")
plt.ylabel("Number of persons")
plt.show()

まとめてグラフにすると見にくい場合は、subplot機能を利用して別々に書くのもアリです。特に感染力が低い場合は、バラバラに描かないと、とても見にくいです。

In [3]:
plt.figure(figsize=(20, 4))  # 横サイズ、縦サイズ
plt.subplots_adjust(wspace=0.25, hspace=0.6)

plt.subplot(1, 4, 1)  #縦1横4枠の、1枚目
plt.plot(S, label="Susceptible", lw=2, c="blue") # S: 感受性者数
plt.legend()
plt.grid()
plt.xlabel("Elapsed day")
plt.ylabel("Number of persons")

plt.subplot(1, 4, 2) #縦1横4枠の、2枚目
plt.plot(I, label="Infected", ls="--", lw=2, c="red") # I: 感染者数
plt.legend()
plt.grid()
plt.xlabel("Elapsed day")
plt.ylabel("Number of persons")

plt.subplot(1, 4, 3) #縦1横4枠の、3枚目
plt.plot(R, label="Recovered", ls="-.", lw=2, c="green") # 回復者数
plt.legend()
plt.grid()
plt.xlabel("Elapsed day")
plt.ylabel("Number of persons")

plt.subplot(1, 4, 4) #縦1横4枠の、4枚目
plt.plot(D, label="Dead", ls=":", lw=3, c="k") # 死亡者数
plt.legend()
plt.grid()
plt.xlabel("Elapsed day")
plt.ylabel("Number of persons")
plt.show()

配置が気にくわない場合は、

  • plt.subplot(2, 2, n): 縦2横2枚の枠
  • plt.subplot(4, 1, n): 縦4横1枚の枠

などの工夫で、描画方法を変更できます。

4. ピーク感染者数、累計感染者数、死亡者数の取得

グラフを見ることで大雑把な傾向がわかります。ここで、社会運営を行う上で重要となる指標、

  • ピーク感染者数（あるいは、最大感染者数）
  • 累計感染者数（一度でも感染した人の総数）
  • 累計死亡者数

を算出してみます。ピーク感染者数とは、「最も感染が拡大しているときに、何人の感染者がいるか」という値になります。この値が一定値を超えると、病院が捌ける数を超え、医療崩壊の状態となります。

In [4]:
MaxInfectors = int(max(I)) # リストIの最大値（を整数で表現）
TotalInfectors = int(I[-1]+R[-1]+D[-t]) # シミュレーション終了時点の、S以外の数
TotalDeath = int(D[-1]) # リストDの最後の値（を整数で表現）

print("ピーク感染者数: ", MaxInfectors, "人")
print("累計感染者数: ", TotalInfectors, "人")
print("累計死者数: ", TotalDeath, "人")

ピーク感染者数:  3127 人
累計感染者数:  9366 人
累計死者数:  94 人

グラフと合わせてみると、40日目付近で、最大3127人もの感染者が出るそうです。総人口が1万人ですから、とても感染力が強いことがわかります。また、累計感染者数が9366人ですから、人口のほとんどが感染してしまったことも見て取れます。一方で、累計死者数は94名となりました。感染者数が多い割に、死者数はそこまで多くないのかもしれません（致死率を1%としたことを思い出すと、だいたい死亡者は1%となっていることがわかります）。

5. 課題

(1) 上のシミュレーションの感染力、致死率をベースとします。ここで、変異株が登場し、感染力が$n$倍、致死率が$m$倍のウイルスが登場した場合、ピーク感染者数、累計感染者数、死亡者数はどう変化するでしょうか？

(2) 人々がマスクや3密の回避を徹底することで、感染力が$1/n$倍になりました。この際、ピーク感染者数、累計感染者数、死亡者数はどう変化するでしょうか？

(3) その他、他に色々パラメータを変化させて考察しましょう。例えば、

  • 平均感染期間が長いウイルスが出た場合は？（20日間にしたいなら、除去率$\gamma=1/20$）

などです。感染力、除去率、致死率を変更すると、色々なウイルスを想定することができます。

6. 卒業研究概要集への手引き

  • このテーマを卒業研究のテーマにしたい場合は、ここの対応をお願いします（必ず事前に相談してください）。
  • 以下、すべて含めて、A4用紙2枚にまとめてください。
  • 厳密に守る必要はありませんが、文章量の比率を目安にしてください（15%とは、A4用紙2枚分の15%という意味）。

1章 はじめに:（文章量の比率: 15%）

  • 感染伝播シミュレーションの必要性を論じてください。
  • SIRモデルを用いたどのような先行研究があるか、まとめてください（SIRDではなく、SIRでok）。
  • 先行研究は、簡単にで良いので、2, 3個紹介してください。

2章 SIRDモデル:（文章量の比率: 30%）

  • SIRDモデルの概要を説明してください。
  • 数式を書く自信があるなら書く。自信がないなら、日本語のみでうまく説明する。

3章 シミュレーション:（文章量の比率: 40%）

3.1節 概要:

  • なんのために、どんなパラメータを採用してシミュレーションを行うのか、説明
  • 1パターンだけではなくて、複数パターン必要です（結果同士を比較して、こっちよりこっちの方がいいなど、考察が必要なため）。

3.2節 結果と考察:

  • S, I, R, Dのグラフ、ピーク感染者数・累計感染者数・死亡者数を表で記載し、説明
  • そこから何が読み取れるのか、考察

4章 おわりに:（文章量の比率: 15%）

  • 今回の分析に対する感想などを記載してください。

参考文献:

参考にした資料を、2〜3件記載してください。以下、書き方の例です。

  • [1] 大前ほか, XXXに関する分析, XXX学会論文誌, 2020.
  • [2] XXXに関する情報, http：//xxx.ddd.ttt.com, 2020年4月20日閲覧

本文中で引用する場合は「大前らはXXXを実施している [1]。また、〜」のように、どこで引用したのか明白にしてください。