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 형 감염증 수리 모델