一昨日のニュース。2020年2月26日，基本方針が全く具体的でないと厳しく指摘されてあせる政府の要請により，PerfumeやEXILEのコンサートが中止になった。また，3月11日まで国立劇場の公演は延期や中止，国立博物館・美術館が閉館となり，社会の空気が相転移した。一方，大阪でただ1名の感染報告者のガイドさんが，回復・退院後再び症状が現れてしまった。さらに，昨日は安倍官邸が全国の小中学校の休校を要請するに至った。SNSでは根拠が不確かな情報がひろがっているということなので，このblogもあまり信用しないほうがいい（ソースコードは提供しており反証は可能だ）。
さて，パラメタが４つあれば象の絵が描けるといわれているので，５つもパラメタを含んだモデルで遊ぶのは慎んだほうがよいのだ。しかし，精緻な理論化に取り組む前に，現象の本質につながるトイモデルであれこれ試すのは，物理学あるいはコンピュータ・シンキング（プログラミング的思考じゃないダス）のキモである。前回のモデルにおけるパラメタセットの問題点を指摘し，それを除く修正モデルを提案して今後の感染状況の推移の全体像の可能性を探りたい。
前回のSIIDRモデルには６つのパラメタが現れる。このうち集団の人口$n$は固定されており，これは問題ない。感染率（１人１日あたりの感染人数）$\beta$や軽症感染者からの平均離脱日$\alpha_1, \alpha_2$，重症感染者からの平均離脱日$\gamma_1, \gamma_2$の５つの範囲をどのように合理的に推定するかがポイントだ。
東京大学の黒田玲子産業医が個人的な見解としてまとめたメモによれば，新型コロナウイルス感染症の潜伏期は1〜12.5日（多くは5〜6日）で，感染力は1人→2〜3人，予後はほぼ治癒だが，ICU入院が5%弱で致死率が0.7%とある。したがって，軽症感染者の平均離脱日を7日とするのはまあまあだ。また，重症感染者への移行率は5%とした。黒田メモの分母は感染者だが，我々のモデルの感染者には発症していないものも含む（5%より小さくなる）。その一方でICUまで至らない重症感染者も考えられる（5%より大きくなる）ので，とりあえず5%を採用する。そこで，$\alpha_1=7/0.95, \alpha_2=7/0.05$とする。また，重症感染者についても同じ7日を仮定し，重症感染者からの死亡率を5%とした。つまり，$\gamma_1=7/0.95, \gamma_2=7/0.05$である。その結果，このモデルの感染者致死率は0.25%となるが，ここでの軽症感染者には発症していないものを含めているので，上記の致死率0.7%とは矛盾しない。こうして４つのパラメタは合理的な範囲で収められる。
最も重要な不確定因子は感染率$\beta$である。これは，ウイルスとホストの生物学的な特性でほぼ決まるが，同時に，ホストである人間集団の社会構造（人口密度，経済活動，生活習慣，文化的要因）にも強く影響されるので，想定する集団によって大きく変わりうる。一方で，感染力（基本再生産数 Ro）は2.2と報告されているため，軽症感染者の平均離脱日数7日で割れば，$\beta=0.3$あたりがこの値の想定範囲となる。
とりあえず，前回のモデルをダイヤモンドプリンセスに当てはめて計算してみた。初期値として5%が感染したところから初めた（2/28/2020日経朝刊の対談より）。$\beta=0.3$のときに3700人あたり4人の死亡数（30日目）が再現できる。このときの累積重症感染者は3700人に換算して590人であり，報告されている700人とほぼ対応している。
次に，このパラメタセットで中国の湖北省の状況を説明できるかどうかを確認したい。モデルの軽症感染者には感染力を持つ非発症者を含んでいるため，報告されている感染者数とは直接比較できない。一方，軽症感染者数の時間構造と重症感染者数の時間構造は１週間ほどのタイムラグがあるがほぼ同形である。また，死亡者数についてもあいまいさが少ない。そこで，我々は，(1)モデルの重症感染者数（1万人当り）の感染者数飽和時期と(2) 現時点の死亡数（1万人当り）の２要素を湖北省のデータと比較する。2019年12月8日に最初の報告があってから，ほぼ３ヶ月たった現時点で湖北省の感染者数は飽和しそうである。また，湖北省の死亡数は約2700人だ（湖北省総人口は5900万人で人口密度は日本と同程度。また，死亡数の8割が集中する武漢は人口が1100万人で人口密度は愛知県なみ）。湖北省の感染者数飽和時の1万人当りの死亡数は0.5人のオーダーと推定される（特異的な武漢を除く湖北省では0.1人）。
図１ β=0.3(左の山/丘)とβ=0.2(右の山/丘)の重症感染者数(u4)と死亡数(u4)
図１のように，上記のパラメタで3ヶ月後(90日目)に感染者数飽和を与える$\beta$は0.30前後だが，このときの死亡数は20人となり1桁以上違う。一方，$\beta=0.20$として，3ヶ月後の死亡数を0.5人あたりに近づけると感染者数飽和時期（右の山の累積値の飽和）が7ヶ月後になり，このモデルでは湖北省データをうまく説明できない。その原因は，中国の強力な対策によって感染率$\beta$を制御できたことにあるのではないかと考え，$\beta$が時間$t$の関数として減少する図２のSIIDR2モデルを導入する。パラメタ数をなるべく減らしたいので，ここでは，$\beta(t)=\beta \exp(-t / \lambda)$という１パラメタの指数関数を導入する。$\lambda$は感染対策時間因子であり１〜数ヶ月のオーダーだ。短期間の変化は１次関数で近似される。
図２ 感染率がβ(t)として時間に依存するSIIDR2モデル
図３のように，$\beta=0.30, \lambda=65$で，湖北省データと矛盾しない結果が得られる。これは感染初期に3ヶ月で感染率を1/4まで減らしたことに対応し，実際中国は武漢封鎖などの強力な措置によってこれを達成したのかもしれない。このときの飽和感染者数は70人（1万人当り）であり，武漢に換算すると8万人となる（非発症者を含む）。なお，この感染対策時間因子が$\lambda=100$にとどまれば，重症感染者数の飽和数は数倍に跳ね上がる。
図３ 小さな山/丘(λ=65)が湖北省データに対応，大きな山/丘はλ=100の場合
この湖北省で得られたパラメタをこのまま日本に当てはめると，30日目の現在の重症感染者数は0.65人/1万人，0.05人/1万人である。つまり現時点における全国の重症感染者数はu3=6500人，死亡数はu4=500人に達することになる。しかし，たいへん残念なことに日本には正確なデータがない。通常の肺炎の死亡者は，年間10万人（日平均300人）なので，埋もれて観測されていない可能性も排除できないがそれはないと思いたい。例えば，$\beta=0.20, \lambda=100$とすると（武漢を除く湖北省に近いケース），30日目の現時点においてそれぞれu3=1000人とu4=100人となり（観測されていないものがあるという考え），よりマイルドな結果を与える。この場合，3〜4ヶ月後に飽和する重症感染者数は全国で7〜8万人，最終死亡数は500人のオーダーとなる（パラメタ次第で，結果を象の鼻を動かすように変えられるのがいたい・・・orz）。
P. S. 下記の参考文献[6][7][8]のように，中国からはよりまともな遅延ダイナミカル・シミュレーションの結果が報告されている。日本では誰か計算していないのだろうか・・・
［１］緊急寄稿新型コロナウィルス(COVID-19)感染症（川名明彦・三笠桂一・泉川公一）
［２］2019-nCoVについてのメモとリンク（中澤港）
［３］人口と感染症の数理（稲葉寿）
［４］数理モデルによる感染症流行制御の検証（大森亮介）
［５］新興感染症の国際的伝播を予測する数理モデル（西浦博・木下諒）
［６］CoVID-19 in Japan: What could happen in the future?（medRxivから）
［７］The reproductive number R0 of COVID-19 based on estimate of a statistical time delay dynamical system（medRxivから）
［８］A Time Delay Dynamical Model for Outbreak of 2019-nCoV and the Parameter Identification（arxivから）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()
sky = @ode_def SIIDR2_model begin
du0 = 1 # u0:time
du1 = -β*exp(-u0/λ)*u1*u2/n # u1:Noimmunity(Susceptible)
du2 = β*exp(-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 = u3 # u6:Accumulated Infection
end n α1 α2 β γ1 γ2 λ
n=10000.0 #total number of population
α1=7.0/0.95 #latent to recovery (days/%)
α2=7.0/0.05 #latent to onset (days/%)
β=0.3 #infection rate
γ1=7.0/0.95 #onset to recovery (days/%)
γ2=7.0/0.05 #onset to death (days/%)
λ=65 #pandemic supression time (days)
function epidm(n,α1,α2,β,γ1,γ2,λ)
u0 = [0.0,9999.0,1.0,0.0,0.0,0.0,0.0] #initial values
p = (n,α1,α2,β,γ1,γ2,λ) #parameters
tspan = (0.0,30.0) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end
β=0.20
λ=100
@time so=epidm(n,α1,α2,β,γ1,γ2,λ)
plot(so,vars=(0,4))
plot!(so,vars=(0,5))
#plot!(so,vars=(0,7))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
0 件のコメント:
コメントを投稿