はじめに
昨今はテレビ等でも「新型コロナ肺炎を収束させるためには8割の行動規制が必要である」と言われています。この8割の数字の根拠はSIRモデルと呼ばれる古典的なモデル方程式だそうです。興味本位ですが、SIRモデルを実装して自粛率の効果を確認してみました。
※筆者はSIRモデルについてはWikipedia程度の知識しかないため、もしかしたら間違っている部分があるかもしれません。本記事のソースコードは、あくまでも、行動規制でどの程度感染を抑えられるかを直感的に理解するためのものです。ここでのモデルは非常に簡易的なので、実際の感染者数の予測・対策は難しいと思います。
SIRモデル
SIRモデル（エスアイアールモデル）は、感染症の短期的な流行過程を決定論的に記述する古典的なモデル方程式である。名称はモデルが対象とする
- 感受性保持者（Susceptible）
- 感染者（Infected）
- 免疫保持者（Recovered、あるいは隔離者 Removed）
の頭文字にちなむ。原型となるモデルは、W・O・ケルマック（英語版）とA・G・マッケンドリック（英語版）の1927年の論文で提案された[1]。
ただし、β > 0 は感染率、γ > 0 は回復率（隔離率）を表す（逆数 1/γ は平均感染期間を表す）。
出展 : wikipedia
SIRモデルはS,R,Iの3元連立微分方程式ですが、
$$
S(t)+I(t)+R(t)=const.
$$
であることから、2変数の方程式に置き換えられます。さらに、S>>Iと仮定し、γ=1、R=βSと置くことで、Iのみの方程式に置き換えられます。ここでは、以下の式を支配方程式としました。
$$
\frac{dI}{dt}=I(t)(R-1)
$$
ここでRは基本再生産数と呼ばれる値で、COVID-19については2～3くらいとされているそうです。
文献 → https://www.biorxiv.org/content/10.1101/2020.01.25.919787v1
上式を前進差分で書き下すと、
$$
I(t+Δt)=I(t)+I(t)(R-1)Δt
$$
となります。以降のソースコードではこの式を元にシミュレーションしていきます。
ソースコード
インポート
グラフの可視化のためにmatplotlibを使用します。
from matplotlib import pyplot
パラメータ
R = 2.5
activity = 0.2
x0 = 200
infection_days = 14
dt = 1/infection_days
stop_day = 30
start_day = 150
total_day = 180
| パラメータ | 説明 |
|---|---|
| R | 基本再生産数（通常活動時の感染者一人当たりの新規感染者数） |
| activity | 自粛期間中の行動率（1が通常時） 例えば、8割自粛なら1-0.8=0.2になります。 |
| x0 | 初期感染者数 |
| infection_days | 隔離or回復までの日数 |
| dt | 計算上のタイムステップ |
| stop_day | 自粛開始日（この日以降行動率が1→activityになる） |
| start_day | 自粛を終了して通常の活動を再開する日 |
| total_day | 総シミュレーション日数 |
変数の初期化
x = [x0]
new = [0]
| 変数名 | 説明 |
|---|---|
| x[day] | t=dayにおける感染者数 |
| new[day] | t=dayにおける新規感染者数 |
変数の更新
for day in range(1,total_day):
p = 1
if day >= stop_day:
p = activity
if day >= start_day:
p = 1
x.append(x[day-1]+x[day-1]*(R*p-1.0)*dt)
new.append(x[day-1]*R*p*dt)
pは行動率を表し、dayがstop_dayになるまでは1、stop_dayからstart_dayまではactivity、それ以降は1になります。この値をRに掛けることで、自粛時の基本再生産数を表しています。例えば、p=0.5であればR×p=2.5×0.5=1.25となります。
グラフの表示
pyplot.plot(x,label='Infected persons')
pyplot.plot(new,label='New infected persons')
pyplot.xlabel('days')
pyplot.ylabel('Persons')
pyplot.legend()
pyplot.show()
行動8割減時（activity=0.2）の感染者推移
自粛開始日（stop_day）から感染者数が急速に低下しましたが、完全に撲滅する前に活動開始日（start_day）を迎えてしまい、その日以降また徐々に感染者数が増加しています。感染者数が減ってきても普段してはダメだとわかりますね。
行動6割減時（activity=0.4）の感染者推移
グラフが見辛くなるので、start_day=180（再活動なし）としました。テレビとかでよく見るグラフですね。activity=0.4とするとR×p=1となり、感染者一人当たり平均一人に感染させることになります。R=2.5の前提下では、6割減で現状維持、不確実性を考慮して確実に感染者数を減らすには8割減が必要ということですね。
さいごに
繰り返しになりますが、本記事のソースコードはあくまでも、行動規制でどの程度感染を抑えられるかを直感的に理解するためのものです。定量的な評価をするにはモデルが簡易的過ぎてきついかなと思います。
APPENDIX コード全容
from matplotlib import pyplot
#基本再生産数（通常活動時の感染者一人当たりの新規感染者数）
R = 2.5
#自粛期間中の行動率（1が通常時）
activity = 0.2
#初期感染者数
x0 = 200
#隔離or回復までの日数
infection_days = 14
#計算上のタイムステップ
dt = 1/infection_days
#自粛開始日（この日以降行動率が1→activityになる）
stop_day = 30
#自粛を終了して通常の活動を再開する日
start_day = 150
#総シミュレーション日数
total_day = 180
#現在の感染者数
x = [x0]
#新規感染者数
new = [0]
#値の更新（差分法）
for day in range(1,total_day):
p = 1
if day >= stop_day:
p = activity
if day >= start_day:
p = 1
x.append(x[day-1]+x[day-1]*(R*p-1.0)*dt)
new.append(x[day-1]*R*p*dt)
#グラフ表示
pyplot.plot(x,label='Infected persons')
pyplot.plot(new,label='New infected persons')
pyplot.xlabel('days')
pyplot.ylabel('Persons')
pyplot.legend()
pyplot.show()