最近コロナウイルス関連でSIRモデルという言葉をよく聞くようになりました。SIRモデルはインフルエンザなどの感染者のグラフによく当てはまることが知られていて、今回のコロナウイルスの対策にも役立てられています。
自分はたまたま関連する授業を受けていることもあり、100番煎じくらいかもしれませんが中学生でもわかるように簡単にまとめてみようかと思います。最後にPythonのコードでシミュレーションしてみます。
#SIRモデルとは
SIRモデルとは、伝染病が広がる様子を表現する数理モデルです。Sは未感染者、Iは感染者(のべ感染者数ではなく現在感染している人)、Rは回復者を表してます。SIRモデルでは回復者は免疫を持っているのでもう感染はしないと考えます。微積を使わないでやさしくかくと次の三式で表されます。
\Delta S = - \beta S I
\Delta I = \beta S I - \nu I
\Delta R = \nu I
化学を知っている人なら反応速度の式(質量作用の法則)と同じようなものだと思ってください。$\Delta S$、$\Delta I$、$\Delta R$はそれぞれ$S$、$I$、$R$の変化分を表しています。具体的に言うと、明日どれだけ$S$、$I$、$R$がどれだけ増えたり減ったりするかと言うことです。例えば、明日の$S$は$S+\Delta S$になります。$\beta$は感染率、$\nu$は回復率というパラメータでどんな病気かによって決まります。$1/\nu$はおよそ回復までにかかる日数になります。
ここで、この三式の両辺を足してみましょう。
\Delta S + \Delta I + \Delta R = 0
となっていることが確認できます。$\Delta S + \Delta I + \Delta R$は明日増えたり減ったりする$S$、$I$、$R$の合計を表します。それが$0$だということは、このモデルが総人口が増えたり減ったりしてしまうような、へんなモデルではないということです。今回、$S+I+R=1$とします。例えば$I=0.1$のときは、その国の住民の10%が$I$(感染者)だということです。
#SIRモデルの意味
次に式の意味を考えてみましょう。一つ目の式は次のようでした。
\Delta S = - \beta S I
この式は、明日、S(未感染者)が$\beta S I$だけ減るということを表しています。減るということは感染者に変わるということです。減る量が$S$と$I$の積で決まる理由は、手っ取り早くは具体的な値を代入してみればイメージが湧きます。
全員が感染していて未感染者がいない国であれば、$S=0$となるので$SI=0$となり、同様に感染者がいない国($I=0$)でも$SI=0$です。未感染者と感染者が半々の国では$SI=0.5\times 0.5=0.25$となりもっとも大きくなります。つまり、$SI$は感染者と未感染者の遭遇率のようなものを表しています。
以上をまとめると、この式は遭遇率と感染率$\beta$の積の分だけ新たに感染する、ということ表しているということがわかります。次に三つ目の式をみてみましょう。
\Delta R = \nu I
この式は回復者が感染者の数に回復率をかけた分だけの人が回復するという意味です。回復数が感染者の割合に依存するというのは少しへんな気がしますが、SIRモデルではこう仮定します。最後に二つ目の式をみましょう。
\Delta I = \beta S I - \nu I
一つ目の式と三つ目の式の意味を踏まえて考えると、(新たに感染する人)−(新たに回復する人)の分だけ感染者が増えるという意味だとわかりますね。(新たに感染する人)<(新たに回復する人)のときこれは負になるので、感染者の数は減少することがわかります。
#基本再生産数
これから暇なとき書きます
#Pythonでのシミュレーション
説明はこれから暇なとき書きます
import matplotlib.pyplot as plt
import numpy as np
beta = 0.5
nu = 0.3
period = 100
def next_state(s, i, r):
delta_s = - beta*s*i
delta_i = beta*s*i - nu*i
s = s + delta_s
i = i + delta_i
if s < 0:
s = 0
if i > 1:
i = 1
r = 1 - s - i
return s, i, r
def main():
results = []
i = 0.01
r = 0
s = 1 - i - r
l_s = [s]
l_i = [i]
l_r = [r]
results.append([s,i,r])
r0 = beta*s/nu
print(s, i, r)
for t in range(period):
s, i, r = next_state(s, i, r)
results.append([s,i,r])
l_s.append(s)
l_i.append(i)
l_r.append(r)
print("R0:{}".format(r0))
plt.figure()
plt.title("R0:{}".format(format(r0,".3f")))
plt.xlabel('Time')
plt.ylabel('Rate')
x = np.linspace(0, 100, period+1)
plt.plot(x, l_s, label="S")
plt.plot(x, l_i, label="I")
plt.plot(x, l_r, label="R")
plt.legend()
plt.show()
main()
実行結果
