あの神戸大学の岩田健太郎さんが新型コロナウィルス感染症の数理シミュレーションの論文を書いていた[1]。調べてみると基本的なSEIRモデルを使った簡単な微分方程式系だ。感染症に対して免疫を持たない者(Susceptible),感染症が潜伏期間中の者(Exposed),発症者(Infectious),感染症から回復し免疫を獲得した者(Recovered)から構成され,その頭文字をとってSEIRモデルと呼ばれる。このうち集団全体の出生率と死亡率を0に置いたモデルで計算していた。
このモデルの詳細については,Compartmental models in epidemiologyに詳しいが,ここでは上の日本語版WikippediaのSEIRの記事で十分だ。早速,Juliaで追計算してみた。下の例では,n=1000人の集団に1人の感染者がいるところからはじまる。感染率βは1日に1人で,潜伏期間αは12日,発症期間γは8.5日である。このモデルでは潜伏期間には感染しない。コロナウイルスは潜伏期間にも感染するようなので,この部分は簡単にモデルを修正できる。以下のコードではδとして導入した項だ(ピークがかなり前倒しになる)。
αとγは制御できないので,感染率βを如何にコントロールするかが対策のキモとなる。潜伏期間に感染するδ=1.0のモデルでは,β=0.2で集団の9割以上が感染してしまうが,β=0.1で8割,β=0.07では1年かけて5割強,β=0.05では1年後にも1割以下となり,非常にセンシティブだ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()
sky = @ode_def SEIR_model begin
du1 = -β*u1*(δ*u2+u3)/n # u1:Susceptible
du2 = β*u1*(δ*u2+u3)/n -u2/α # u2:Exposed
du3 = u2/α -u3/γ # u3:Infectious
du4 = u3/γ # u4:Recovered
end n α β γ δ
function dif()
# hayashi-model
#n=1000.0 #total number of population
#α=5.5 #latent period (days)
#β=1.0 #disease transmission rate
#γ=7.0 #infectious period (days)
#δ=0.0 #latent transmission rate
# iwata-model
n=1000.0 #total number of population
α=12.0 #latent period (days)
β=1.0 #disease transmission rate
γ=8.5 #infectious period (days)
δ=0.0 #latent transmission rate
u0 = [999.0,1.0,0,0] #initial values
p = (n,α,β,γ,δ) #parameters
tspan = (0.0,100.0) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
plot(sol,vars=(0,1))
plot!(sol,vars=(0,2))
plot!(sol,vars=(0,3))
plot!(sol,vars=(0,4))
end
@time dif()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
岩田さんが例にあげているアウトブレイクするパラメタを代入してみると,よく似たグラフが得られたが,発症者のピーク位置や感染者と発症者の山の重なりが微妙に違う。彼の論文ではパラメタの領域しか提示されていないため,厳密な比較ができない。そこで,林祐輔さんのシミュレーション[2]と比較した。こちらのほうとは,ピタリと一致したのでたぶん上のコードは正しいのだろう。
図 JuliaによるSEIRモデルの確認
[1]A Simulation on Potential Secondary Spread of Novel Coronavirus in an Exported Country Using a Stochastic Epidemic SEIR Model(K. Iwata, C. Miyakoshi)
[2]COVID-19/stochasticSEIRmodel.ipynb(Y. Hayashi)
[3]面向新冠疫情的数据可视化分析与模拟预测(陈宝权 北京大学前沿计算研究中心)
[4]corona_virus_model fitting_by optina(H. Maruyama)
感染症の数理シミュレーション(2)に続く