感染症の数理シミュレーション（８）

感染症の数理シミュレーション（７）からの続き

中間まとめ
いろいろと書き散らしてきたので，ここまでに得られた結果を整理してみる。

感染症の数理シミュレーション（３）において，感染対策時間因子を修正した図１のSIIDR2モデルで，湖北省（武漢以外）の新規感染者数累計と死亡数累計を説明できるパラメタセットを得ることができた。

図１　SIIDR2モデルの概念図，$n=u_1+u_2+u_3+u_4+u_5$, $p+q=1$

\begin{equation}
\begin{aligned}
S: \dfrac{du_1}{dt} &= -\dfrac{\beta(t)}{n} u_1 u_2\\
I_1: \dfrac{du_2}{dt} &= \dfrac{\beta(t)}{n} u_1 u_2 -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2}\\
I_2: \dfrac{du_3}{dt} &= \dfrac{u_2}{\alpha_2} - \dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2}\\
D: \dfrac{du_4}{dt} &= \dfrac{u_3}{\gamma_2}\\
R: \dfrac{du_5}{dt} &= \dfrac{u_2}{\alpha_1} + \dfrac{u_3}{\gamma_1}\\
I_a: \dfrac{du_6}{dt} &= \dfrac{u_2}{\alpha_2}\\
\beta(t)&= \dfrac{\beta}{15}\{8+7 \tanh(-\frac{t-\tau}{\lambda})
\end{aligned}
\end{equation}
湖北省（武漢以外）のデータを再現するパラメタセットは以下のとおりであり，これを用いたシミュレーション結果は図２のように与えられる。
$\alpha_1=5/0.80, \alpha_2=5/0.20=25$, $\beta=0.915$, $\gamma_1=15/0.96, \gamma_2=15/0.04=375$,  $\lambda=\tau=7, \nu=0.025$，$u_2$が軽症感染者数（発症無），$u_3$が重症感染者数（発症有），$\dfrac{u_2}{\alpha_2}$が新規重症感染者数$f(t)$，$\dfrac{u_3}{\gamma_2}$が新規死亡数$h(t)$，$u_4$が死亡数累計$i(t)=\int h(t) dt$，$u_6$が新規重症感染者累計$g(t)=\int f(t) dt$，である。WHOで報告されているデータ[1]では，感染者累計（Confirmed）が$g(t)$，死亡数累計（Deaths）が$i(t)$に対応している。

図２　湖北省（武漢以外）のSIIDR2モデルによるシミュレーション

（○はConfirmedとDeathsの観測値，WHOとtencentニュースサイトのデータ）

半分パラメタフィッティングなので，$\nu,\tau,\lambda$を調整して勝手に象の鼻を描いている感は否めないのが微妙なところである。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()

sky = @ode_def SIIDR2_model begin
  du0 = 1 # u0:time
  du1 = -β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n # u1:Noimmunity(Susceptible)
  du2 = β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n -u2/α1 -u2/α2 # u2:Mild(Infected-a)
  du3 = u2/α2 -u3/γ1 -u3/γ2 # u3:Serious(Infected-b)
  du4 = u3/γ2 # u4:Dead
  du5 = u2/α1 +u3/γ1 # u5:Recovered
  du6 = u2/α2 # u6:Accumulated Infected-b
end n α1 α2 β γ1 γ2 λ τ

function epidm(β,ν,λ,τ,T)
n=10000.0 #total number of population
α1=5.0/0.80 #latent to recovery (days/%)
α2=5.0/0.20 #latent to onset (days/%)
#β=0.45 #infection rate (/day・person)
γ1=15.0/0.96 #onset to recovery (days/%)
γ2=15.0/0.04 #onset to death (days/%)
u0 = [0.0,n-11ν,4ν,2ν,0.0,5ν,ν] #initial values
p = (n,α1,α2,β,γ1,γ2,λ,τ) #parameters
tspan = (0.0,T) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end

β=0.915 #infection rate
ν=0.025 #inital value of accumulated infected-b
λ=7 #pandemic supression range (days)
τ=7 #pandemic supression start (day)
T=40 #period of simulation

xc=[0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48]
yc=[0.003,0.03,0.23,0.67,1.27,2.00,2.52,2.88,
3.44,3.58,3.71,3.69,3.71,3.71,3.71,3.71,3.71]
zc=[0.0,0.0,0.0,0.01,0.02,0.03,0.04,0.05,
0.07,0.09,0.10,0.11,0.11,0.12,0.12,0.13,0.13]

# kohoku-bukan model
# β=0.915,ν=0.025,λ=7,τ=7,α2=5.0/0.20,γ2=15.0/0.04
@time so=epidm(β,ν,λ,τ,T)
#plot(so,vars=(0,2))
plot(so,vars=(0,3))
plot!(so,vars=(0,4))
plot!(so,vars=(0,5))
plot!(so,vars=(0,6))
plot!(so,vars=(0,7))
plot!(xc,yc,st=:scatter)
plot!(xc,zc,st=:scatter)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

［１］Coronavirus disease (COVID-2019) situation reports
［２］新型環状病毒肺炎疫情動態（tencentニュースサイト）
［３］新型コロナウイルス感染症対策専門家会議(第5回)（3月2日）
［４］特設サイト：新型コロナウイルス（NHK）
［５］新型コロナウイルス感染速報（Su Wei）
［６］Coronavirus Disease (COVID-19) – Statistics and Research（Our World in DATA）
［７］Databrew's COVID-19 data explorer（Databrew）
［８］COVID-19情報共有（佐藤彰洋）
［９］時間遅れを考慮した確率 SIR モデルの安定性解析（石川昌明）
［１０］分布的時間遅れをもつ確率感染症モデルの安定性解析（石川昌明）

P. S. [6]のように，プロがシミュレーションをすると遅延SIRモデルみたいなことになるのか。感染率にステップ関数を導入しているので鋭い構造が出現しているようだ。遅延確率SIRモデルになるとちょっとついていけません(3/17/2020)。


感染症の数理シミュレーション（９）に続く