SIRコンパートメントモデルで感染症の流行予測

数学セミナーの９月号を買いました。今月号の特集は「新型コロナウイルスと闘うために数学ができること」。感染症予測に関する記事が７つ載っています。

f:id:z_logger:20200915013732j:plain — 数学セミナー9月号

今回は、これらの記事を参考に実際にプログラムを作って感染者の時間変化を計算してみようと思います。

注意：ただの素人が雑誌の記事を参考に計算しているだけなので、間違いがある可能性があります。詳しく知りたい方は、数学セミナーの９月号を買って読んでください。

f:id:z_logger:20200915013830j:plain — 特に参考にした記事は「国内の流行初期のデータによる予測とその評価」

SIRコンパートメントモデル

数学セミナーでは、SIR（SEIR）コンパートメントモデルという感染症のモデルが解説されています。これは、全人口を

感受性人口（感染する可能性のある人たち）S
潜伏期 E
感染 I
回復 R

に分けてそれぞれの人口の経時変化を計算するものです。

具体的には、単位時間内に一人のIが一人のSと接触する確率をβ、発症率をε、単位時間内に一人のIが治癒してRに移行する確率をγとして、各人口の時間変化は、

$$ \begin{eqnarray} \frac{ dS }{ dt } &=& -\beta S(t) I(t) \\ \frac{ dE }{ dt } &=& \beta S(t) I(t) - \epsilon E(t) \\ \frac{ dI }{ dt } &=& \epsilon E(t) - \gamma I(t) \\ \frac{ dR }{ dt } &=& \gamma I(t) \end{eqnarray} $$

で計算できるそうです。実際は感染した全員を検出できるわけではないので、検出率やデータの誤差、再感染を考慮したモデルもあるようですが、ここではこの式だけで簡易的に計算してみます。また、基本再生産数R₀（一人の感染者が起こす二次感染の総数）はβ/γで計算できます。

今回はCOVID-19を例に計算してみるので、潜伏期間と感染期間は記事の値を参考に、ε=0.2, γ=0.1（潜伏期間５日、感染期間１０日）に固定して計算します。

基本再生産数R₀を変えて計算してみる

R₀=2.5（感染者１人が平均2.5人の人に感染させる）の場合、

f:id:z_logger:20200914234719p:plain — R₀=2.5の場合。S(0)=0.9999, I(0)=0.0001で計算。

赤が感染者数の推移を表しています。I(0)=0.0001としたので、日本だとするとだいたい1000人ぐらいの人が感染したところを0日としています。感染者数は約７０日目でピークに達して、最終的には約90%の人が感染します。実際の基本再生産数の値は一定ではなく、常に変化するので非現実的な計算ですが、実際の感染者数のグラフを見ると非常によく似ていて興味深いです。

f:id:z_logger:20200915014120p:plain — 日本の新型コロナウイルス感染者数の時間変化。SmartNewsより引用。

R₀=5.0とした場合の結果が以下。R₀=2.5のときと違って30日で感染のピークに達して、二か月ぐらいで完全に収束しますね。ただ流行のピーク時は全人口の半分が感染することになるので大変なことになりそう。

f:id:z_logger:20200915012848p:plain — R₀=5の場合

R₀=1.5とした場合。感染は１００日以上かけてとても緩やかに進みます。１年ぐらいかけて６０％ぐらいの人が感染することになりますね。

f:id:z_logger:20200915012922p:plain — R₀=1.5の場合

Stay Homeの効果

R₀=2.5を基準にしてStay Homeの効果を検証してみます。感染爆発の初期（50日目）でStayHomeした場合と、感染者数がほぼピークに達した時(70日目)からStay Homeした場合で計算してみます。なお、本記事で行った計算では、簡単にStayHomeするとβが半分になるとして計算を行いました。