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

[Geoff Quickfall]

Interesting info Angel.

I will get it installed and start playing with it.

Under COVID a new definition of R (recovered) has been proposed. Actually it has been an alternative for years.

R = removed (not recovered) as a death affects R in the same manner as recovered; both removing and no longer infected or susceptible.

[Ángel Martin]

(04-08-2020 03:56 PM)Geoff Quickfall Wrote:  Under COVID a new definition of R (recovered) has been proposed. Actually it has been an alternative for years.

R = removed (not recovered) as a death affects R in the same manner as recovered; both removing and no longer infected or susceptible.

Yes, that's exactly how it's labeled in this model, see text:
"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."

and:
"Besides, the Removed section includes both the recovered (cured) and the dead patients."

The attached RAW file should facilitate things.

[SlideRule]

"                                             An Epidemic Model
The SIR model is a simple model, due to Kermack and McKendrick, of an epidemic of an infectious
disease in a large population. We assume the population consists of three types of individuals, whose
numbers are denoted by the letters S, I and R (which is why this is called an SIR model). All these
are functions of the time t, and they change according to a system of differential equations.

• S is the number of susceptibles, who are not infected but could become infected.
• I is the number of infectives. These individuals have the disease and can transmit it to the
susceptibles.
• R is the number of removed individuals. These may or may not have the disease, but they can't
become infected and they can't transmit the disease to others. They may have a natural immunity, or
they may have recovered from the disease and are immune from getting it again, or they may have the
disease but are incapable of transmitting it (e.g. because they may have been placed in isolation), or
they may have died. The mathematical model doesn't distinguish among those possibilities.

The model we will consider assumes a time scale short enough that births and deaths (other than deaths
from this disease) can be neglected."

[attachment=8324]
a visualization

Might I recommend The Mathematical Theory of Infectious Diseases and its application, 2e (ISBN 0 85264 231 8)

BEST!
SlideRule

... je souhaite seulement que dans une question qui regarde de si près
le bien de l'humanité, on ne décide rien qu'avec toute la connoissance
de cause qu'un peu d'analyse & de calcul peut fournir.
                                                                  Daniel Bernoulli, 1760

[Ángel Martin]

(04-09-2020 12:24 AM)Ángel Martin Wrote:  The attached RAW file should facilitate things.

RAW file removed - it's been replaced by a new ROM with the SIR model using ODEs plus a direct method using a combined SOLVE/INTEG for numerical evaluation...

Edited: Here it is, make sure both MOD files are plugged to use the direct method (LBL "R<T>") as both FINTG and FROOT are used. Most recommendable to use TURBO mode in V41 (or a CL w/ Turbo50) due to the slow numerical process involved.

The ODE-based method is accessed with LBL "SIR" (general population size N) and LBL "SIR2+" (unitary population size N=1, i.e. ratios for S, I, and R).

Thanks to Jean-Marc for his suggestions and contributions to this module.

[Ángel Martin]

Updated MOD file is attached, it now offers two versions for the direct method:

1. LBL "R=RT" - using FOCAL routines for Solve & Integ, self-contained
2. LBL "R<T>", using FINTG and FROOT in the Solve & Integrate ROM

The SIR model also has two versions as follows:

1. LBL "SIR", population N
2. LBL "SIR2+", population N=1

A few additional routines are thrown in as well, such as EASTER - written by Kari Passanen, the EKG Calculations , the LDL Cholesterol formula and the totally unrelated Audio Tape Counter programs from the forum.