71日目 SIRモデルでこの自粛がいつまで続くか予測してみました

#SIRモデルを計算してみました。

新型コロナの流行はいつまで続くんだろう。
昨日のNHKスペシャル[「新型コロナウイルス 瀬戸際の攻防 〜感染拡大阻止最前線からの報告〜」](https://www.nhk.jp/p/special/ts/2NY2QQLPM3/episode/te/BLNZP53V36/)を参考に計算してみました。

##[SIRモデル](https://ja.wikipedia.org/wiki/SIR%E3%83%A2%E3%83%87%E3%83%AB)とは
感染症の流行〜収束を計算するのに使われる方程式です。

##PythonでSIRモデルを実装する
Pythonで実装したコードが[scipython.com](https://scipython.com/)にありましたので、使わせていただきました。


[The SIR epidemic model](https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/)

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

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 1, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.2, 1./10 
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

# Initial conditions vector
y0 = S0, I0, R0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, axis_bgcolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()
```

#インフルエンザの場合
まず、例として設定されていたインフルエンザをみてみます。

##インフルエンザ(免疫なし）
```
N  =1000 #人口、1000人というと学校や会社、町内会あたり？
I0 =1    #感染者。最初は１人。
R0 =0    #回復して免疫を獲得した人。最初は0人
S0 =N - I0 - R0 #免疫をもたず感染しやすい人。最初はN-I0-R0人
beta =0.2       #1日に感染させる人数
gamma =1./10    #回復率
t =160          #期間(日数) 
```

ひとりの感染者が1日0.2人に感染させます。
感染期間は10日間。
回復して免疫を獲得すると、この後は接触しても感染しません。
免疫ありが半分以上になると収束に向かいます。

![i0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/b8458f6e-8180-0db6-968e-b3d11a7bf95c.png)

ひとりのキャリアからピークまで67日、15.3%、153名が感染し収束となりました。
インフルエンザで学級閉鎖するとき、こんな感じですよね。

##インフルエンザ(免疫20%)
インフルエンザのシーズン前にワクチンを打って免疫を獲得したケースを考えてみます。
まず、２割の人が免疫を獲得したところにインフルエンザが流行した場合です。

```
# Initial number of infected and recovered individuals, I0 and R0.
#I0, R0 = 1, 0
I0, R0 = 1, (N*0.2)
```

![i2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/d7c677f9-f8b1-eee8-971c-cc0d606bf47f.png)

２割の人が免疫を持つだけで、流行がずいぶん抑えられました。ワクチンすごい。

##インフルエンザ(免疫50%)

2019年度の[インフルエンザワクチンの定期接種率](https://www.jmedj.co.jp/journal/paper/detail.php?id=13945)は50%でした。

```
# Initial number of infected and recovered individuals, I0 and R0.
#I0, R0 = 1, 0
I0, R0 = 1, (N*0.5)
```

![i1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/18ddb184-6d30-afcb-5e1d-d24fea783eb3.png)

完全に抑えられています。
こうなれば安心ですね。

#COVID-19の場合

新型コロナの場合を考えてみます。
東京都、武漢市など1000万人の都市を想定、感染期間は３週間(21日)としました。
その他の性質はインフルエンザと同じとしました。

##COVID-19(免疫なし）
```
N  =10000000 #人口1000万人(東京、武漢を想定)|
I0 =1        #感染者。最初は１人。
R0 =0        #回復して免疫を獲得した人。最初は0人
S0 =N - I0 - R0 #免疫をもたず感染しやすい人。最初はN-I0-R0人
beta =0.2       #1日に感染させる人数 
gamma =1./21    #回復に要する日数の逆数 
t =160          #期間(日数) 
```

![c0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/d1f79adb-6f40-7558-7552-28a83efe67cc.png)


感染期間が21日になるだけで、ピーク時の感染者数が42%、420万人と凶悪なウイルスに進化しました。
このうち2割が発症し、1%が死亡となると・・・
84万人が発症、4万2千人が死亡、228日で収束します。


地獄です。

このグラフを見た後、数日ほど滅入ってしまいました。

しかしこのツイートを見かけて、ひょっとしたら希望はあるかもと思いました。

> 衝撃の調査結果‥
ドイツで最初に新型コロナの集団感染が起こったハインスベルク郡で、免疫調査をまず約500人に行ったところ、何と約15％の人が感染済みだったと‥。
ハインスベルクのPCR検査での人口10万人あたりの感染者数は604.4人、判明していた感性率はたった0.6%です‥
https://twitter.com/R3000C/status/1248594913589002240

##COVID-19(免疫15%）
今の日本も15%の人が感染済みとします。
15%の人が感染済みになるまで89日。

89日目でロックダウンして流行を止めます。

１ヶ月後に全ての人が回復し15%の人が免疫を持つとします。
ロックダウンを解除。流行が再燃するのですが、、、

![c1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/232b6c0e-e392-c243-311a-a8c199c4b598.png)

第二波のカーブは15%の免疫によりゆるくなりました。

再びピークの手前、感染率が15%でロックダウンすることにします。
今度は119日目でした。

１ヶ月後に全員が回復し免疫を獲得、合計30%の人が免疫を持ちます。

##COVID-19(免疫30%）

![c2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/30132458-c566-adc6-9794-1c82ff112361.png)

第三波はピークまで178日。最大20%の感染率です。
15％の感染率まで149日でした。

##COVID-19(免疫45%）
第四波はピークまで252日、最大11%の感染率です。

![c3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/356263/314f2318-3f18-3831-2d45-dc07b0002329.png)

ここまでくればゆるい自粛で乗り越えられそうです。

#この自粛はいつまで続く？
15%感染→ロックダウン→回復→解除を３セットで
(89日+30日)＋(119日＋30日)+(149日＋30日)=447日

447日間ひきこもって生き延びる工夫が必要になりますが、犠牲者は少なくすみます。

経済を回しながらガンガン感染させる作戦だと、99%命に別状はなく228日でケリがつくわけで。
悪魔の誘惑に耐えられない人がいるかも・・・しれませんね。

武漢市は4月7日に封鎖を解除しました。
第二波が来ても今回ほどの規模にはならないはず。
免疫獲得率が気になります。
徹底した管理によって、ぐっと流行を抑えられる可能性もあります。
そうなると世界で最初に復興するのが武漢市、ここから先はSFでしょうか。

そしてワクチン。これがあれば一気に状況がひっくり返ります。
手洗い、うがい、ソーシャルディスタンス。ぐっとカーブが抑えられます。
それから検査。人口の半分が免疫を獲得できれば完全に抑え込めるのです。

長い戦いになりそうです。