Virus Epidemics SIR Model

[Ángel Martin]

Chances are you're reading this from your "shelter in place" Covid-19 confinement, and thus will resonate with the concepts involved. This short program is a direct application of Jean-Marc Baillard's ODE routine "RK4C" from the Differential Equations module, available at TOS and the CL Library (id$ = "DIFF").

The SIR model calculates the values of Susceptible, Infected and Removed groups of population in a total population of N individuals. N = S + I + R.

The virus is modeled with two parameters that indicate the infection rate "a" and the recovery rate "b". These are the crux of the model, as they need to be expressed in the same units used by the ODEs to show the individual results - i.e. a=2.3 infected people per person and DAY, and b = 0.3 people recovered per DAY if we want to look at daily numbers.

Supposedly your local Covid-19 statistics could be used to estimate the parameters, but this is tricky since the reported infection cases are much lower than the actual ones. Besides, the Removed section includes both the recovered (cured) and the dead patients.

The ODEs are normalized by the total population size N, but this is transparent to the user and it's done by the routine itself. The ODE's are implemented in the LBL "d/dt" subroutine:

dS/dt = -a. I(t).S(t)
dI/dt = a. S(t). I(t) - b.I(t)
dR/dt = b.I(t)

So here you have it, play with your a,b parameters (XEQ A) and define your initial conditions N, So, Io (XEQ B) to start getting the results for the S,I,R sectors of population. The model assumes Ro=0. Remember that N=S+I+R at all times.

Program Listing:

Code:

 9:41AM 04/08
 01*[b]LBL "SIR"[/b]
 02*[u]LBL A[/u]
 03 RCL 12
 04 "a="
 05 ARCL X
 06 >"?"
 07 PROMPT
 08 STO 12
 09 RCL 13
 10 "b="
 11 ARCL X
 12 >"?"
 13 PROMPT
 14 STO 13
 15*[u]LBL B[/u]
 16 RCL 14
 17 "N="
 18 ARCL X
 19 >"?"
 20 PROMPT
 21 STO 14
 22 RCL 02
 23 *
 24 "S0="
 25 ARCL X
 26 >"?"
 27 PROMPT
 28 RCL 14
 29 /
 30 STO 02
 31 RCL 03
 32 RCL 14
 33 *
 34 "I0="
 35 ARCL X
 36 >"?"
 37 PROMPT
 38 RCL 14
 39 /
 40 STO 03
 41 0
 42 STO 04
 43 STO 01
 44 ,1
 45 STO 05
 46 10
 47 STO 06
 48*LBL C
 49 RCL 04
 50 RCL 14
 51 *
 52 RCL 03
 53 RCL 14
 54 *
 55 RCL 02
 56 RCL 14
 57 *
 58 RCL 01
 59 "S"
 60 ARCLI
 61 "`="
 62 ARCL Y
 63 PROMPT
 64 "I"
 65 ARCLI
 66 >"="
 67 ARCL Z
 68 PROMPT
 69 "R"
 70 ARCLI
 71 >"="
 72 ARCL T
 73 PROMPT
 74 "d/dT"
 75 ASTO 00
 76 [color=#0000CD][b]XROM "RK4C"[/b][/color]
 77 GTO C
 78 RTN
 79*[b]LBL "d/dT"[/b]
 80 RDN
 81 X<>Y
 82 *
 83 LASTX
 84 RCL 13
 85 *
 86 X<>Y
 87 RCL 12
 88 *
 89 ENTER^
 90 CHS
 91 X<>Y
 92 RCL Z
 93 ST- Y
 94 END

Note: ARCLI is in the AMC_OS/X Module. It appends the integer part o X to ALPHA.

PS. Looks that the formatting codes don't work within a CODE block?

PSS. - SIR material abounds on the Web, here's one video that was useful to me to grasp some details:
https://www.youtube.com/watch?v=k6nLfCbAzgo&t=867s