# 全系横断演習 SIRモデル
## 今後の方針
* [08/05 15:32] Miyazaki Hayato今後の方針
1) 稲葉先生のPDFのワクチン接種を反映したモデルの数値計算をできるようになる。
2) 次を行う:
a) どの感染症を扱うか
b) ワクチンの効果について調べる
c) 対象地域を決定し、ワクチン接種前の感染者、感染予備軍、免疫獲得者の割合を決定する。
## 仕様
* 人口の増加(出生,移住),減少(死亡など)を考慮
| 変数 | 意味 | 備考 |
| - | -| -|
| t | 時間 | day |
| S | 感染予備軍 | |
| I | 感染者 | |
| R | Recover | |
## link
### SIRモデル参考
* [SIRモデルの基本的なシミュレーション](https://wagtail.cds.tohoku.ac.jp/coda/python/p-6-application-sup-ode-sir-model.html)
* [緊急事態宣言は早すぎた?SIRモデルで検証(おもしろいby瀧本)](https://qiita.com/mashino-satoshi/items/ad59057dab52308eebdd)
### 人口
* [一日に起きる出来事一覧(出生,死亡など)](https://www.mhlw.go.jp/wp/hakusyo/kousei/11-3/dl/1jp.pdf)
## 目標
* ワクチンによる感染が収まるために、「どの程度の時間がかかるか」をおおよそわかるようにする
a) ワクチン接種の割合をかえる
b) ($S$, $I$, $R$) の比を変えてシミュレーションしてみる
c) 赤ん坊にのみワクチンを打つ
etc...
## 課題
* ワクチンの接種割合(?)v をどういう風にしていくか(指数関数, 二次関数...)
* どの地域について絞るか
## 意見箱
新型コロナ版SIRモデルwith vaccine
$$
\begin{align}
\frac{dS(t)}{dt} ={}& (1-v(t))b - \mu S(t) - \underbrace{\beta(t)}_{\text{ワクチンの効果を入れる?}} S(t)I(t) \\
\frac{dI(t)}{dt} ={}& \underbrace{\beta(t)}_{\text{ワクチンの効果を入れる?}} S(t)I(t) - (\mu + \gamma )I(t) \\
\frac{dR(t)}{dt} ={}& vb - \mu R(t) + \gamma I(t)
\end{align}
$$
ワクチンの効果を入れるなら、$v$と$\beta$をいじる?
時間経過の要素を入れるなら、$t$の関数として。
例:$v(t)=t/1000$. $\beta$ を $\beta-v(t)$に置き換え。
## 次回までの課題
適当な$v$を見つけてくる。
* $\beta$ は定数と考えてv(t)は$\beta$に近づいていく関数がイイかも
* $v(t)=t/1000$のような関数として、$\beta$になったらそれ以上増えない関数を設定する?もしくは、現実的に$v(t)$が$\beta$になるという意味は、全員がワクチン接種したという意味なので、そこでシミュレーションを終える。
## 出生と死亡を考慮したSIRモデルについて
出生と死亡を考慮したSIRモデルは以下のように表せる:
$$
\begin{align}
\begin{aligned}
\frac{dS(t)}{dt} ={}& bN(t) - \frac{\beta}{N(t)} S(t)I(t) - \mu S(t), \\
\frac{dI(t)}{dt} ={}& \frac{\beta}{N(t)} S(t)I(t) - (\mu + \gamma )I(t), \\
\frac{dR(t)}{dt} ={}& \gamma I(t) -\mu R(t).
\end{aligned}
\tag{SIR1}
\end{align}
$$
ここで $N(t)$ は総人口であり, $N=S+I+R$ である. また, 各定数の意味は以下である:
| 定数 | 意味 | 備考 |
| - | -| -|
| $b$ | 出生率 | /day |
| $\beta$ | 感染率 | /day |
| $\mu$ | 自然死亡率 | /day |
| $\gamma$ | 回復率 | /day |
多くの論文では, 閉鎖人口を仮定 ($N(t)$が一定) したり, 出生率と死亡率が等しいことを仮定 ($b=\mu$) したりしている. これは以下のように考えれば, 自然な仮定である:
未知関数を
$$
\tilde{S} = \frac{S}{N}, \quad \tilde{I} = \frac{I}{N}, \quad \tilde{R} = \frac{R}{N}
$$
と取り直すと, (SIR1)は
$$
\begin{align}
\begin{aligned}
\frac{d\tilde{S}(t)}{dt} ={}& b - \beta \tilde{S}(t)\tilde{I}(t) - b \tilde{S}(t), \\
\frac{d\tilde{I}(t)}{dt} ={}& \beta \tilde{S}(t)\tilde{I}(t) - (b + \gamma )\tilde{I}(t), \\
\frac{d\tilde{R}(t)}{dt} ={}& \gamma \tilde{I}(t) -b \tilde{R}(t).
\end{aligned}
\tag{SIR2}
\end{align}
$$
と書き直すことができる. 詳しい計算は [Trawicki(2017)](https://www.mdpi.com/2227-7390/5/1/7/htm) のAppendix Aを参照. ざっくり言うと, 積の微分と, (SIR1)の3つの方程式を足すと得られる総人口に関する方程式
$$
\frac{dN(t)}{dt} = (b-\mu)N(t) \tag{TP}
$$
を用いると(SIR2)を導出できる. つまり総人口$N(t)$が時間によって変化する場合でも
$$
\begin{align}
\tilde{S}(t) + \tilde{I}(t) + \tilde{R}(t) = 1 \tag{TP1}
\end{align}
$$
と総人口を1, つまり一定値とすることができ, 死亡率$\mu$も出生率$b$と一致する.
また, (TP)はすぐに解くことができて
$$
\begin{align}
N(t) = N(0)e^{(b-\mu)t}. \tag{TP2}
\end{align}
$$
したがって, ワクチンを考慮したり数値計算したりするときは(SIR2)を扱い, 必要に応じて(TP2)で総人口を求めればよい. ワクチンを考慮したSEIRSモデルが [Trawicki(2017)](https://www.mdpi.com/2227-7390/5/1/7/htm) で考察されているので, これを簡略化すればSIRモデルの場合も出る. また, $\mu$は自然死亡率なので, [高須夫悟先生のスライド](http://gi.ics.nara-wu.ac.jp/~takasu/lecture/global/20-MathBiol-9.pdf)の21ページのように感染症による死亡率を入れればよいと思う. 21ページのスライドで出生率が$\mu K$で自然死亡率とずれているのは, 時間が経つと総人口が$K$に近づくように設定しているから. このように設定したい場合, (TP1)を忘れて
$$
\hat{S} = K\tilde{S}, \quad \hat{I} = K \tilde{I}, \quad \hat{R} = K \tilde{R}
$$
とおき, (SIR2)に代入すると
$$
\begin{align}
\frac{d\hat{S}(t)}{dt} ={}& bK - \tilde{\beta} \hat{S}(t)\hat{I}(t) - b \hat{S}(t), \\
\frac{d\hat{I}(t)}{dt} ={}& \tilde{\beta} \hat{S}(t)\hat{I}(t) - (b + \gamma )\hat{I}(t), \\
\frac{d\hat{R}(t)}{dt} ={}& \gamma \hat{I}(t) -b \hat{R}(t)
\end{align}
$$
となる. ここで $\tilde{\beta} = \frac{\beta}{K}$である. 実際, $\hat{N}(t) = \hat{S}(t)+\hat{I}(t)+\hat{R}(t)$ とおき, 上の3つの方程式を足し合わせると
$$
\frac{d\hat{N}(t)}{dt} = b K - b \hat{N}(t)
$$
が得られるので, 変数分離で解くと
$$
\hat{N}(t) = K - (K - \hat{N}(0))e^{-b t}
$$
となる. また, $N=K$ は $\frac{d\hat{N}(t)}{dt}=0$ を満たす. つまり, たとえ(TP1)を忘れて総人口が一定でないとしても, 結局時間が経てば不動点$K$に近づいていく($K$は安定な不動点であるという)ことがわかる.