Excel VBA で SIR モデルを実装する。

##新型コロナウイルス(COVID-19）と SIR モデル
新型コロナウイルス(COVID-19）の患者数推移を表す数理モデル「SIRモデル」をよく目にするようになりました。
SIR モデルでは未感染者をS (susceptible)、感染者を I (infectious)、治癒者や死亡者を R (removed) とします。

S, I, R 各群の増加速度は以下のように表せます。

\begin{align}
\frac{dS}{dt} &= -βSI\\
\frac{dI}{dt} &= βSI-γI\\
\frac{dR}{dt} &= γI\\
\end{align}

$\ β$ は感染率、$\ γ$ は除去率（平均罹病期間の逆数）となります。

\ γ=\frac{1}{平均罹病期間}

$\ 全人口をNとすると基本再生産数\ R_0は$

\ R_0=\frac{βN}{γ}\\
\therefore β=\frac{γR_0}{N}

以上より、2つの係数 $\ γ, β$ は平均罹病期間と基本再生産数 (COVID-19 では2~3)から算出できることが分かります。

上記の微分方程式から、

\begin{align}
\Delta S&=-βSI\Delta t\\
\Delta I&=(βSI-γI)\Delta t\\
\Delta R&=γI\Delta t\\
\ S(t+\Delta t)&=S+\Delta S\\
\ I(t+\Delta t)&=I+\Delta I\\
\ R(t+\Delta t)&=R+\Delta R\\
\end{align}

##Excel VBA マクロのコード
上記から Excel VBA マクロのプログラムを作りました。
数学的に微分方程式の解を求めなくても、for ループで各パラメータを漸増させていけば曲線が得られます。
セルB2 を (S) 総人口1億2000万人
セルB3 を (I) 13852人（2020.4.29の患者数)
セルB4 を (R) 3763人（2020.4.29の治癒者＋死亡者数)
セルB5 を基本再生産数 2.5
セルB6 を平均罹病期間 14日
セルB7 を観察期間 180日とします。
$\Delta t=0.5日$ とします。
基本再生産数を色々変えてみて下さい。
以下にスクリーンショットとコードを示します。

 Sub SIR()
 If ActiveSheet.ChartObjects.Count = 1 Then ActiveSheet.ChartObjects(1).Delete
 Range(Columns(4), Columns(7)).Clear
        S = Range("B2").Value
        I = Range("B3").Value
        R = Range("B4").Value
          saiseisansu = Range("B5").Value
         ribyoukikan = Range("B6").Value
          kansatsukikan = Range("B7").Value
         Dim beta As Double
         Dim gam As Double
         Dim deltaS As Double
         Dim deltaI As Double
         Dim deltaR As Double
         Dim deltaT As Double
         Dim myCount As Integer
         gam = 1 / ribyoukikan
         beta = saiseisansu * gam / S
         deltaT = 0.5
         myCount = 0
          Cells(1, 4) = "days"
           Cells(1, 5) = "susceptible"
           Cells(1, 6) = "infetcted"
           Cells(1, 7) = "removed"
          For myTime = 0 To kansatsukikan Step deltaT
           deltaS = -beta * S * I * deltaT
           deltaI = (beta * S * I - gam * I) * deltaT
           deltaR = gam * I * deltaT
            S = S + deltaS
            I = I + deltaI
            R = R + deltaR
            Cells(2 + myCount, 4) = myTime
             Cells(2 + myCount, 5) = S
             Cells(2 + myCount, 6) = I
             Cells(2 + myCount, 7) = R
             myCount = myCount + 1
          Next myTime
          
          Set trange = Worksheets("Sheet1").Range(Cells(1, 4), Cells(1 + myCount, 7))
      Charts.Add
      With ActiveChart
            .HasLegend = True
            .Legend.Position = xlLegendPositionTop
            .ChartType = xlXYScatterLines
            .SetSourceData Source:=trange, PlotBy:=xlColumns
            .Location Where:=xlLocationAsObject, Name:="Sheet1"
      End With
      With ActiveChart
            .Axes(xlValue, xlPrimary).HasTitle = True
            .Axes(xlValue, xlPrimary).AxisTitle.Characters.Text = "value"
            .Axes(xlValue, xlPrimary).AxisTitle.AutoScaleFont = True
            .Axes(xlCategory, xlPrimary).HasTitle = True
            .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time"
            .Axes(xlCategory, xlPrimary).AxisTitle.AutoScaleFont = True
            .PlotArea.Interior.ColorIndex = 2
            .Axes(xlValue).MajorGridlines.Delete
      End With
      With ActiveSheet
            .ChartObjects(1).Top = 50
            .ChartObjects(1).Width = 350
            .ChartObjects(1).Height = 350
            .ChartObjects(1).Left = 600
      End With
          
 End Sub

COVID-19 に SIR モデルを当てはめる際の注意点 をご参考ください。
###参考サイト
稲葉 寿：感染症数理モデルとCOVID-19 | 日本医師会 COVID-19有識者会議
SIR型の感染症数理モデル