03
Writing the deterministic SIR equations
This lesson develops confidence in working with model equations. We begin with biological statements, identify flows into and out of each compartment, form the differential equations, substitute numerical values and finally translate the same calculations into Python.
Biological scenario
An epidemic has reached the state
\[S=900,\qquad I=80,\qquad R=20.\]
The total population is \(N=1000\). We use \(\beta=0.3\) and \(\gamma=0.1\). We want to determine the instantaneous direction and rate of change of every compartment.
The general balance principle
A compartment changes because quantities enter or leave it. This gives a principle used throughout mathematical biology:
We apply this principle separately to \(S\), \(I\) and \(R\).
1. Construct the susceptible equation
dS_dt = -infection_flowThe minus sign is essential. Infection is an outgoing flow from \(S\), so it makes the susceptible population decrease.
2. Construct the recovered equation
dR_dt = recovery_flowThe recovery flow is incoming to \(R\), so it has a positive sign.
3. Construct the infectious equation
The infectious compartment has both an incoming flow and an outgoing flow.
Words first
rate of change of infectious people = infection flow in − recovery flow out
\[\frac{dI}{dt}=\beta\frac{SI}{N}-\gamma I.\]
Read this equation as:
“The rate of change of the infectious population with respect to time equals the infection flow into the infectious compartment minus the recovery flow out of it.”
Meaning of the infectious equation
| Expression | Meaning |
|---|---|
| \(\displaystyle \frac{dI}{dt}\) | Instantaneous rate of change of the infectious population |
| \(\beta\) | Transmission parameter |
| \(SI/N\) | Susceptible–infectious interaction adjusted for population size |
| \(\displaystyle \beta\frac{SI}{N}\) | Infection flow into \(I\) |
| \(\gamma\) | Recovery-rate parameter |
| \(\gamma I\) | Recovery flow out of \(I\) |
| minus sign | Recovery removes people from the infectious compartment |
The complete SIR system
\[\frac{dS}{dt}=-\beta\frac{SI}{N},\]
\[\frac{dI}{dt}=\beta\frac{SI}{N}-\gamma I,\]
\[\frac{dR}{dt}=\gamma I.\]
These equations are connected. The same infection flow leaves \(S\) and enters \(I\). The same recovery flow leaves \(I\) and enters \(R\).
Substitute the numerical values by hand
Substitution means replacing symbols with their given numerical values.
Infection flow
\[\beta\frac{SI}{N}=0.3\times\frac{900\times80}{1000}=21.6.\]
Recovery flow
\[\gamma I=0.1\times80=8.\]
\[\frac{dS}{dt}=-21.6,\qquad \frac{dI}{dt}=21.6-8=13.6,\qquad \frac{dR}{dt}=8.\]
Translate the calculation into Python
Step 1: Write the mathematical quantities
S = 900
I = 80
R = 20
N = S + I + R
beta = 0.3
gamma = 0.1Step 2: Write the flow formulas
infection_flow = beta * S * I / N
recovery_flow = gamma * IPython uses * for multiplication and / for division.
Step 3: Write the equations
dS_dt = -infection_flow
dI_dt = infection_flow - recovery_flow
dR_dt = recovery_flowThe variable names dS_dt, dI_dt and dR_dt are readable Python versions of the derivative notation.
Run the equations
Output
Run the code to see the result.
Interpret the signs
| Derivative | Value | Interpretation at this state |
|---|---|---|
| \(dS/dt\) | −21.60 | − Susceptible population is decreasing |
| \(dI/dt\) | +13.60 | + Infectious population is increasing |
| \(dR/dt\) | +8.00 | + Recovered population is increasing |
The signs describe the direction of change at the current state. The infection flow exceeds the recovery flow, so \(I\) is increasing. These rates will change as the compartment sizes change.
Why the total population remains constant
Adding the equations makes the repeated flows cancel:
\[\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=-\beta\frac{SI}{N}+\beta\frac{SI}{N}-\gamma I+\gamma I=0.\]
Python verifies the same result using total_rate.
Understanding the conservation check
tolerance = 1e-12
assert abs(total_rate) < toleranceMathematically, total_rate should equal exactly zero. However, computers store many decimal numbers approximately, so a correct calculation can sometimes produce an extremely small value near zero.
- 1e-12 means \(1\times10^{-12}=0.000000000001\). This is the accepted numerical tolerance.
- abs(total_rate) gives the distance of total_rate from zero, whether the small value is positive or negative.
- < tolerance asks whether that distance is sufficiently small to be treated as zero.
- assert checks that the condition is true. If it is false, Python stops with an AssertionError.
We do not normally use total_rate == 0 for decimal calculations because a harmless computer-rounding difference could make that exact comparison false.
Place the equations in a reusable function
After the individual equations are understood, we can group them:
def sir_rates(S, I, R, beta, gamma):
N = S + I + R
infection_flow = beta * S * I / N
recovery_flow = gamma * I
dS_dt = -infection_flow
dI_dt = infection_flow - recovery_flow
dR_dt = recovery_flow
return dS_dt, dI_dt, dR_dtThe function does not change the mathematics. It packages the same equations so that a numerical method can evaluate them repeatedly at different epidemic states.
Try it yourself
Change beta = 0.3 to beta = 0.05.
- Substitute the new value into the infection-flow formula by hand.
- Predict the sign of \(dI/dt\).
- Run the Python code and compare it with your calculation.
- Explain why the total rate is still zero.