2020年3月8日日曜日

感染症の数理シミュレーション(6)

感染症の数理ミュレーション(5)からの続き

湖北省(武漢以外)のデータでは,新規/累計感染者数 f(t), g(t) と新規/累計死亡数 h(t), i(t) を簡単な近似関数で表現した。これを用いると,感染症の数理シミュレーション(3)のSIIDR2モデルにおけるパラメータを絞り込むことができるのではないだろうか。

そこで,前回示したSIIDR2の方程式系と今回の実測データ近似関数の関係を整理してみよう。
\begin{equation}
\begin{aligned}
S: \dfrac{du_1}{dt} &= -\dfrac{\beta}{n} u_1 u_2 = -\dfrac{\beta}{n}\alpha_2 f u_1 \\
I_1: \dfrac{du_2}{dt} &= \dfrac{\beta}{n} u_1 u_2 -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2} = \dfrac{\beta}{n}\alpha_2 f u_1 - (1 + \dfrac{\alpha_2}{\alpha_1}) f \\
I_2: \dfrac{du_3}{dt} &= \dfrac{u_2}{\alpha_2} - \dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2} = f - (1 + \dfrac{\gamma_2}{\gamma_1}) h \\
D: \dfrac{du_4}{dt} &= \dfrac{u_3}{\gamma_2} = h\\
R: \dfrac{du_5}{dt} &= \dfrac{u_2}{\alpha_1} + \dfrac{u_3}{\gamma_1} = \dfrac{\alpha_2}{\alpha_1} f + \dfrac{\gamma_2}{\gamma_1} h\\
I_a: \dfrac{du_6}{dt} &= \dfrac{u_2}{\alpha_2} = f
\end{aligned}
\end{equation}
これから次の関係式が得られる。
\begin{equation}
\begin{aligned}
u_1 &= n - u_2 - u_3 - u_4 - u_5\\
u_2 &= \alpha_2 f\\
u_3 &= \gamma_2 h\\
u_4 &= \int h dt = i\\
u_5 &= \dfrac{\alpha_2}{\alpha_1} g + \dfrac{\gamma_2}{\gamma_1}i \\
u_6 &= \int f dt = g
\end{aligned}
\end{equation}
さらに,次の前提条件が成り立つとする。
①感染者の2割が発症する[1]。
②致命率(発症者の死亡率)は0.7%である[2]。
③発症までの平均潜伏期間は$\alpha=5$とする$(1/\alpha=1/\alpha_1+1/\alpha_2)$ [3]。
④基本再生産数$R_o$=2.2(1.4〜3.9)より,$\beta=R_o/\lambda=0.46\pm0.12$  である [2]。
⑤湖北省(武漢以外)の新規死亡数h(t)は20日のピークに0.006/日,40日目の累計で0.12人(1万人当り)[4]。
⑥同上の新規感染者数f(t)は12日のピークに0.25/日,40日目の累計で3.70人(1万人当り)[4]。

$I_2$の$u_3$に対する微分方程式$du_3/dt + u_3/\gamma = f$を解いて求めた$u_3$の関数形が$\gamma_2 h$に近づくように定めると,$\gamma=5,\   \gamma_1=5/0.965,  \gamma_2=5/0.035$が得られる。ただし,$(1/\gamma=1/\gamma_1+1/\gamma_2)$である訂正有り 3/13/2020)→感染症の数理シミュレーション(7)

各関数のT日目の値がわかれば,非感染者の減少分$\Delta u_1$は次式で与えられる。
\begin{equation}
\Delta u_1=n-u_1(T)=u_2(t)+u_3(T)+u_4(T)+u_5(T)
\end{equation}
さらに湖北省(武漢以外)のT=40日のデータを代入して,$u_2(T)=\alpha_2 f(T)≒0$, $u_3(T)=\gamma_2 h(T)≒0$より,$\Delta u_1= \dfrac{\alpha_2}{\alpha_1} g(T) + \dfrac{\gamma_2}{\gamma_1}i(T)+i(T) =14.8+3.4=18.2 $である。ここで,$\alpha_1=5/0.8, \alpha_2=5/0.2$という①の仮定を用いた。

$u_1(t)$において,$\beta$が定数であると仮定して微分方程式を解けば,$u_1(t)=n \exp (-\frac{\beta}{n} \alpha_2 g(t))$となることから,$\beta$と$\alpha_2$を与えれば,$u_1(T)=9982$が得られる。先ほどの式と組み合わせると$\beta=0.2$となる。

まとめると,湖北省(武漢以外)のデータを再現するSIIDR2モデルのパラメタのセットは次のようになる。訂正有り 3/13/2020)→感染症の数理シミュレーション(7)
\begin{equation}
\alpha_1=\dfrac{5}{0.8},\  \alpha_2=\dfrac{5}{0.2},\  \beta=0.2,\  \gamma_1=\dfrac{5}{0.965},\   \gamma_2=\dfrac{5}{0.035}
\end{equation}

SIIDR2モデルでは,感染率を時間の関数$\beta(t)$として感染対策が行われることを表している。上記の$\beta$はそれを平均したものと考えることができる。ところで,上記の$I_2$の$u_2$に対する微分方程式は,$du_2/dt + u_2/\alpha = \beta u_2$と近似できる($\therefore u_1/n ≒ 1$)ことから,次のように与えられる。
\begin{equation}
\beta(t) = \dfrac{du_2/dt}{u_2} + \dfrac{1}{\alpha}
\end{equation}
$u_2(t)=\alpha_2 f(t)$を上式に代入したものと,感染症の数理シミュレーション(5)で考えたモデル感染対策時間因子で$\beta=0.91, \lambda=7, \tau=7$と置いたものを下記の図に示す。


図:湖北省(武漢以外)のデータからSIIDR2で導いたβ(t)(下)とモデル化したβ(t)(上)

また,このモデル化された感染対策時間因子$\beta(t)$をt=0からt=40まで積分したものを40日で割って平均すると($T=40, \beta=0.91, \lambda=7, \tau=7$),
\begin{equation}
\bar{\beta} = \dfrac{1}{T} \int_0^T \dfrac{\beta}{15} (8 + 7 \tanh (-(t-\tau)/\lambda ) dt = 0.24
\end{equation}
となって,先ほど湖北省(武漢以外)のデータからSIIDR2モデルを経由して導いた$\beta=0.2$と矛盾しないような感染対策時間因子のパラメタセットが存在することがわかる。

[1]新型コロナウイルス感染症対策の見解(厚生労働省)
[2]新型コロナウイルス(COVID-19)感染症(川名・三笠・泉川)
[3]新型コロナウイルス肺炎COVID-19(下畑享良)
[4]tenetニュースサイト
[5]Coronavirus disease (COVID-2019) situation reports(WHO)
[6]COVID-19(奥村晴彦)

感染症の数理シミュレーション(7)へ続く

0 件のコメント: