小田垣さんのSIQR

朝日新聞に，九州大学名誉教授の小田垣さんの計算として「PCR検査を倍にすれば、接触「5割減」でも収束可能？」という記事がでた。twitterでも注目を集めていた。

小田垣さんのホームページにその論考があった[1]。SIQRモデルということで，感染者を2段階に分けていた。我々のSIIDR2モデルと本質的に同じではないか。感受性保持者S(t)から直接我々の重症患者I2(t)，すなわち小田垣さんの隔離感染者Q(t)に遷移する項があって，ここは違うのだが，議論が始まる前の段階でこの項を落としているので結局同じです。違うのはこちらには死亡数D(t)への遷移を含んでいることくらいである。

そこで両方のモデルで使用しているパラメタを比較してみた。左が小田垣さんのSIQR，右が我々のSIIDR2である。
\begin{equation}
\begin{aligned}
N\beta \quad (0.07) &= \beta \quad (0.4-0.6)  \\
p \quad (0.096) &= \dfrac{1}{\alpha_2} \quad (0.04) \\
\gamma \quad (0.04) &= \dfrac{1}{\alpha_1} \quad (0.16) \\
\gamma \quad (0.04) &= \dfrac{1}{\gamma_1}\quad  (0.064) \\
noparameter (0) &= \dfrac{1}{\gamma_2}\quad   (0.00267) \\
\end{aligned}
\end{equation}
感染率の$\beta$が一桁違うことがわかる。その他はファクターが異なるくらいである。I(t)からR(t)への遷移時間とQ(t)からR(t)への遷移時間が同じとしていることもやや疑問に感じる。また，小田垣氏は論考の中で，基本再生産数（実効再生産数）は$N\beta/p$でなく$N\beta/(p+\gamma)$とすべきだと主張しているが，これはどうなのだろうか。

とりあえず，感染初期のS(t)=Nのときにあてはめて，接触5割減の議論に持ち込んでいる部分は興味深いが，実際にこのパラメータを使って日本の感染データの全体像を説明できるのだろうか。我々のコードを少し修正して彼らのモデルとパラメタを使った計算を再現してみた（初期値の細かな調整は我々のモデルを前提としているがとりあえずはそこはそのままにしておく）。

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

sky = @ode_def SIQR_model begin
  du0 = 1 # u0:time
  du1 = -β*u1*u2/n # u1:Noimmunity(Susceptible)
  du2 = β*u1*u2/n -u2/α2 -u2/α1 # u2:(Infected)
  du3 = u2/α2 -u3/γ1 # u3:(Quarantined)
  du4 = u3/γ2 # u4:Dead(not in use)
  du5 = u2/α1 +u3/γ1 # u5:Recovered
  du6 = u2/α2 # u6:Accumulated Quarantined
  du7 = u3/γ1 # u7:Accululated Recovered
end n α1 α2 β γ1 γ2 λ τ

function epidm(β,ν,λ,τ,T)
n=10000.0 #total number of population
α1=1/0.04   #5.0/0.8 #latent to recovery (days/%)
α2=1/0.096  #5.0/0.2 #latent to onset (days/%)
β=0.07  #0.45  #infection rate (/day・person)
γ1=1/0.04   #15.0/0.96 #onset to recovery (days/%)
γ2=15.0/0.04   #onset to death (days/%) (not in use)
u0 = [0.0,n-11ν,4ν,2ν,0.0,5ν,ν,0.0] #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

#japan-data(start=3/1)
xj=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,
14,15,16,17,18,19,20,21,22,23,24,25,26,27,
28,29,30,31,32,33,34,35,36,37,38,39,40,41,
42,43,44,45,46,47,48,49,50,51,52,53,54,55,
56,57,58,59,60,61,62]
yj=[0.0190,0.0202,0.0213,0.0225,0.0252,0.0277,0.0324,
0.0361,0.0387,0.0408,0.0451,0.0492,0.0536,0.0568,
0.0619,0.0646,0.0658,0.0658,0.0693,0.0754,0.0790,
0.0830,0.0864,0.0895,0.0947,0.103,0.110,0.119,
0.134,0.148,0.155,0.173,0.189,0.208,0.232,
0.260,0.290,0.310,0.338,0.378,0.424,0.477,
0.536,0.576,0.607,0.643,0.681,0.728,0.777,
0.822,0.853,0.882,0.912,0.946,0.983,1.018,
1.046,1.062,1.078,1.099,1.118,1.133,1.154]
zj=[0.040,0.048,0.048,0.048,0.048,0.048,0.048,
0.048,0.056,0.071,0.095,0.119,0.151,0.167,
0.175,0.190,0.222,0.222,0.230,0.262,0.278,
0.286,0.325,0.333,0.341,0.357,0.365,0.389,
0.413,0.429,0.444,0.452,0.452,0.516,0.548,
0.556,0.579,0.635,0.643,0.675,0.698,0.746,
0.778,0.810,0.865,0.944,1.079,1.175,1.222,
1.278,1.357,1.476,2.198,2.278,2.516,2.651,
2.762,2.786,2.984,3.087,3.294,3.429,3.603]/100
plot(xj,yj,st=:scatter,label="Confirmed-japan")
#plot!(xj,zj,st=:scatter,label="Deaths-japan")

β=0.07
ν=0.01
T=60

@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))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

結果は以下の通りである。確かに初期の段階ではデータを説明しているが，その後の振る舞いは説明できない。このモデル（のパラメタ）はあまりよろしくないのかもしれない。

図　日本の新規感染数累計の推移（3/1-5/1）u6をデータと比較する

［１］新型コロナウイルスの蔓延に関する一考察（小田垣孝，2020.5.5）
［２］隔離と市中の感染者を分ける SIR モデル（佐野雅己，2020.4.29）
［３］3.11以後の科学リテラシー No. 89（牧野淳一郎，2020 科学５月号）
［４］感染症の数理シミュレーション（８）（2020.3.15）