2022年8月22日月曜日

SIDRモデル(1)

 エンデミックからの続き (参考:感染症の数理シミュレーション(2)

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

図1: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}]

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



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

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



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

0 件のコメント:

コメントを投稿