Signature
F = beta * S * I / N * dt; G = gamma * I * dt; S_next = S - F; I_next = I + F - G; R_next = R + G
| Inputs | Definition | Unit |
|---|---|---|
beta | Rate of contacts sufficient for transmission per person | per day |
S | Number susceptible at the start of the step | people |
I | Number infectious at the start of the step | people |
N | Total population, S plus I plus R, above zero | people |
dt | Length of the time step | days |
gamma | Rate at which infectious people recover, the inverse of the mean infectious period | per day |
R | Number recovered or otherwise immune at the start of the step | people |
F | Number moving from S to I during the step | people |
|---|---|---|
G | Number moving from I to R during the step | people |
S_next | Number susceptible at the end of the step | people |
I_next | Number infectious at the end of the step | people |
R_next | Number recovered or immune at the end of the step | people |
Function
Compartment flow balance function of a compartmental model
Maps the numbers of people in the compartments of a population model, the flows between compartments per unit time and any inflow from or outflow to outside the model to the rate of change of each compartment, or, after a time step of length dt, to the next set of compartment counts. The flows can be fixed rates times compartment sizes, as in a static model, or can depend on the number infectious through the force of infection, which makes the model dynamic and nonlinear. A Markov cohort update (HE-FN-MM-001, HE-FM-MM-001) applies one fixed transition matrix to a closed cohort, the agent-based update HE-FN-ABM-001 with the per-contact infection probability HE-FM-ABM-001 follows individuals instead of pools, and a per-cycle probability from a constant rate follows HE-FM-TP-001. Reproduction numbers of the SIR model and of stratified models are in HE-FN-BRN-001 (HE-FM-BRN-001 to HE-FM-BRN-005 and HE-CF-BRN-001), and the final size of an epidemic and R_0 estimated from it are in HE-FN-ATR-001, HE-CF-ATR-001 and HE-FM-ATR-005, so none of them is restated here. Notation follows the Compartmental Model article.
Try this function
Implementations
Excel
SIR Euler step from named compartment and rate cells
With the rates in TransRate and RecoveryRate, the population in PopTotal, the step in StepLength and the start counts in SusNow, InfNow and RecNow, the five cells return the new infections (NewInfections), the recoveries (NewRecoveries) and the end-of-step counts SusNext, InfNext and RecNext. Copying the block down one row per step with each row's end counts as the next row's start counts runs the model forward.
=TransRate*SusNow*InfNow/PopTotal*StepLength; =RecoveryRate*InfNow*StepLength; =SusNow-NewInfections; =InfNow+NewInfections-NewRecoveries; =RecNow+NewRecoveries
Assumptions
Closed population with lasting immunity and an exponential infectious period in the SIR step
There are no births, deaths or migration over the period, recovery gives lasting immunity, and the constant recovery rate gamma implies an exponentially distributed infectious period with mean 1/gamma, the waiting time that Hethcote attaches to terms such as gamma I. Open populations add the external flows b and d of HE-FM-CMPM-001 to each compartment.
Both SIR flows computed from start-of-step counts before any compartment is updated
F and G are computed from S and I at the start of the step and then applied together. With gamma dt or lambda dt close to 1 the step moves most of a compartment at once and can overshoot, so the step length is chosen short enough that a shorter step leaves the results of interest unchanged.
Worked examples
Article's day-1 SIR step without vaccination
The article's illustrative population of 10,000 with 9,900 susceptible and 100 infectious people, beta of 0.3 and gamma of 0.1 per day has 29.7 new infections and 10 recoveries on day 1, leaving 9,870.3 susceptible, 119.7 infectious and 10 recovered, still 10,000 in all.
beta = 0.3; gamma = 0.1; N = 10000; dt = 1; S = 9900; I = 100; R = 0; F = 29.7; G = 10; S_next = 9870.3; I_next = 119.7; R_next = 10
Article's day-1 SIR step after 2,000 susceptible people are vaccinated
A perfect vaccine given to 2,000 susceptible people moves them to the immune compartment before day 1, so R starts at 2,000. Day-1 infections fall to 23.7, 6.0 fewer than without the programme, and the compartments become 7,876.3, 113.7 and 2,010.
beta = 0.3; gamma = 0.1; N = 10000; dt = 1; S = 7900; I = 100; R = 2000; F = 23.7; G = 10; S_next = 7876.3; I_next = 113.7; R_next = 2010
Zero transmission rate leaves recoveries as the only flow
With beta set to zero no one is infected and the infectious compartment only empties, the check the article lists for verifying a compartmental model, a limiting case.
beta = 0; gamma = 0.1; N = 10000; dt = 1; S = 9900; I = 100; R = 0; F = 0; G = 10; S_next = 9900; I_next = 90; R_next = 10
Common errors
Infections for the infectious compartment taken from already updated susceptibles
Updating S first and then computing the inflow to I from the new S mixes two time points. On the article's day 1 it moves 29.6109 people into I but 29.7 out of S, leaving 9,999.9109 people in the population, and the loss repeats at every step (computed here for illustration).
Euler step too long for the speed of the epidemic
Euler steps lag the true solution while infections grow. With one-day steps the article's epidemic has 570.0 infectious people on day 10, about 13% fewer than the 653.9 given by steps of 0.0001 day, which is close to the exact solution, and two five-day steps give 389.6 (computed here for illustration). The error builds over the run and feeds every count and cost taken from it.
Sources
Hethcote on the classic SIR epidemic model with standard incidence
Hethcote HW. The mathematics of infectious diseases. SIAM Review. 2000;42(4):599-653. Section 2.3, The Classic Epidemic Model (pages 604 and 605), which gives the SIR model as dS/dt = minus beta I S/N, dI/dt = beta I S/N minus gamma I and dR/dt = gamma I in equation 2.1, with S(t) + I(t) + R(t) = N, states that it uses the standard incidence and has recovery at rate gamma I corresponding to an exponential waiting time, and that it has no births and deaths because the period is short; and section 2.1 (page 603), which states that terms such as gamma I correspond to exponentially distributed waiting times with mean 1/gamma.
Pitman and colleagues on the Euler method in spreadsheet transmission models
Pitman R, Fisman D, Zaric GS, Postma M, Kretzschmar M, Edmunds J, Brisson M. Dynamic transmission modeling. A report of the ISPOR-SMDM Modeling Good Research Practices Task Force-5. Value in Health. 2012;15(6):828-834. Section Software Options, which states that in a spreadsheet environment the Euler method is often used to project the system of differential equations forward in time and that this method is not as accurate as other numerical techniques; and section Methodological uncertainty, which states that most dynamic transmission modelling has used system dynamics, in which transition between compartments is represented by differential equations.
Canonical Identity
Stable URI · Machine-readable · Resolvable · CC BY 4.0