　　　　On a Thread of the Web: 感染症の数理シミュレーション（１）

2020年2月24日月曜日

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

新型環状病毒肺炎からの続き

あの神戸大学の岩田健太郎さんが新型コロナウィルス感染症の数理シミュレーションの論文を書いていた［１］。調べてみると基本的なSEIRモデルを使った簡単な微分方程式系だ。感染症に対して免疫を持たない者（Susceptible），感染症が潜伏期間中の者（Exposed），発症者（Infectious），感染症から回復し免疫を獲得した者（Recovered）から構成され，その頭文字をとってSEIRモデルと呼ばれる。このうち集団全体の出生率と死亡率を０に置いたモデルで計算していた。

このモデルの詳細については，Compartmental models in epidemiologyに詳しいが，ここでは上の日本語版WikippediaのSEIRの記事で十分だ。早速，Juliaで追計算してみた。下の例では，n=1000人の集団に１人の感染者がいるところからはじまる。感染率βは１日に１人で，潜伏期間αは12日，発症期間γは8.5日である。このモデルでは潜伏期間には感染しない。コロナウイルスは潜伏期間にも感染するようなので，この部分は簡単にモデルを修正できる。以下のコードではδとして導入した項だ（ピークがかなり前倒しになる）。

αとγは制御できないので，感染率βを如何にコントロールするかが対策のキモとなる。潜伏期間に感染するδ=1.0のモデルでは，β=0.2で集団の９割以上が感染してしまうが，β=0.1で８割，β=0.07では１年かけて５割強，β=0.05では１年後にも１割以下となり，非常にセンシティブだ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

岩田さんが例にあげているアウトブレイクするパラメタを代入してみると，よく似たグラフが得られたが，発症者のピーク位置や感染者と発症者の山の重なりが微妙に違う。彼の論文ではパラメタの領域しか提示されていないため，厳密な比較ができない。そこで，林祐輔さんのシミュレーション［２］と比較した。こちらのほうとは，ピタリと一致したのでたぶん上のコードは正しいのだろう。