SIDRモデル（１）

 エンデミックからの続き　（参考：感染症の数理シミュレーション（２））

西浦博さんが紹介している既知のSIRSモデルは，年齢構造や平均余命まで考慮した精緻なモデルだが，一般ピープルがその振る舞いを試してみるには大層なので，単純なSIDRモデルを考えてみた。本質的に，牧野淳一郎さんが最近計算したモデルと等価であり，念のため死亡者数をちょっと取り出しただけだ。

図１：SIDRの概念図

$\dfrac{du_1}{dt} = -\beta \dfrac{ u_1 \cdot u_2}{n} + \dfrac{u_3}{\alpha}$

$\dfrac{du_2}{dt} = \beta \dfrac{ u_1 \cdot u_2}{n} -\dfrac{u_2}{\gamma_1} -\dfrac{u_2}{\gamma_2}$

$\dfrac{du_3}{dt} = \dfrac{u_2}{\gamma_1} -\dfrac{u_3}{\alpha}$

$\dfrac{du_4}{dt} = \dfrac{u_2}{\gamma_2}$

ここで，$u_1+u_2+u_3+u_4=n$は定数で，全国民数である。また，$\alpha$が免疫保持期間（day），$\gamma$は感染期間（day）である。CFRを$c$とすると，$\frac{1}{\gamma}=\frac{1}{\gamma_1} + \frac{1}{\gamma_2}$が成り立ち，$\frac{1}{\gamma_1}=\frac{1-c}{\gamma}，\frac{1}{\gamma_2}=\frac{c}{\gamma}$である。

微分方程式の両辺を$n$で割って，$x=\frac{u_2}{n}，y=\frac{u_3}{n}，z=\frac{u_4}{n}$とする。また，$s=\frac{u_1}{n}$を，$s=1-(x+y+z)$で消去すると次の3変数の微分方程式系が得られる。

$\dfrac{d x}{dt} = \beta (1-x-y-z)\ x -\dfrac{x}{\gamma}$

$\dfrac{d y}{dt} = \dfrac{(1-c)\ x}{\gamma} -\dfrac{y}{\alpha}$

$\dfrac{d z}{dt} = \dfrac{ c\ x}{\gamma}$

この式で$c=0$とすれば，牧野さんの式となる。次のMathematicaコードで検算したところ，彼の結果は再現できた。ここでは，次の値で計算してみる。$R=\beta \gamma =2，\gamma=13，\alpha=130，c=0.0013，x[0]=0.01$

R = 2; \[Gamma] = 13; \[Beta] = R/\[Gamma]; \[Alpha] = 130;
sol = NDSolve[{
x'[t] == \[Beta] x[t] (1 - x[t] - y[t] - z[t]) - x[t]/\[Gamma],
y'[t] == 0.9987 * x[t]/\[Gamma] - y[t]/\[Alpha],
z'[t] == 0.0013 * x[t]/\[Gamma],
x[0] == 0.01, y[0] == 0, z[0] == 0}, {x, y, z}, {t, 0, 1000}]
fx[t_] := x[t] /. sol[[1, 1]]
fy[t_] := y[t] /. sol[[1, 2]]
fz[t_] := z[t] /. sol[[1, 3]]
Plot[fx[t], {t, 0, 1000}, PlotRange -> {0, 0.2}]

その結果は図２のようなものであり，振動の後に平衡状態となり，感染者割合は4.5%に収束する。

図２：SIDRモデルにおける感染者割合の推移

また，400日目を過ぎると，死亡数は図３のように線形で増加するが，このときの1日あたりの死亡者割合は，4.5 x 10^-6 となり，540人/日に相当するので，年間死亡数は19万人に達する。これは季節性インフルエンザの超過死亡数（1万人）の20倍のオーダーである。

図３：SIDRモデルにおける累計死亡数の推移