SIR Epidemic Model
Simulate infectious disease transmission dynamics using the classical Kermack-McKendrick Susceptible-Infectious-Recovered (SIR) model solved with SciPy’s solve_ivp. Durable worksheet parameters drive the continuous numerical integration in the Python for Excel notebook, publishing outbreak timelines, peak load estimates, and herd immunity metrics directly to your spreadsheet.
Result preview
With the baseline parameters (Total Population N = 100,000, Initial Infected I0 = 100, Transmission Rate beta = 0.30, Recovery Rate gamma = 0.10), the basic reproduction number is R₀ = 3.00, indicating rapid spread with a herd immunity threshold of 66.7% (fraction 0.6667). The outbreak reaches a continuous peak of 30,080 active infections on day 38.4, eventually infecting and recovering 92.9% of the population (final attack rate fraction 0.9292) over the 90-day horizon.
What this template does
- Reads population size, initial infected cases, contact rate (\(\beta\)), and recovery rate (\(\gamma\)) from labeled worksheet cells.
- Integrates coupled non-linear differential equations using
scipy.integrate.solve_ivpwith the explicit Runge-Kutta (RK45) method. - Calculates key epidemiological indicators including the basic reproduction number (\(R_0\)), herd immunity threshold, continuous peak infection timing, and final attack rate.
- Visualizes time-series epidemic curves and infection-vs-susceptible phase portraits side-by-side in the task pane.
- Publishes summary KPIs and a daily trajectory table (days 0 to T, up to 90 days) back into Excel.
Why Python
Epidemic models rely on non-linear continuous-time feedback loops where infection rates depend on the product of susceptible and infected populations (\(S \cdot I\)). Static spreadsheet formulas or simple discrete difference equations can introduce numerical divergence, artificial delays, or inaccurate peak estimates. SciPy’s solve_ivp applies adaptive step-size error control (rtol=1e-8, atol=1e-8), delivering smooth interpolation and reliable turning-point identification without requiring manual calculus approximations in Excel formulas.
Operating workflow
Author: A quantitative modeler or public health analyst maintains the differential equations, integration tolerances, scenario test vectors, and published workbook outputs.
Workbook user: A policy analyst, healthcare planner, or educator modifies transmission assumptions and reviews the recalculated epidemic peak and daily trajectory table without writing Python.
When used as a standard forecasting or scenario evaluation tool, the workbook can open in App mode for general operators. App mode changes the presentation only; it does not change permissions, source access, or workbook trust. App mode streamlines the interface while preserving the underlying calculations and workbook model.
Download the Excel template
Inputs and assumptions
Edit population counts, transmission rate, or recovery rate on the Simulation sheet. The integrated trajectory, peak timing, and phase curves recalculate immediately.
The parameter table located at Simulation!A4:B10 contains:
| Parameter | Notation | Baseline Value | Description |
|---|---|---|---|
| Total Population | \(N\) | 100,000 | Total population size (\(S + I + R\)) |
| Initial Infected | \(I_0\) | 100 | Initial count of infectious individuals |
| Initial Recovered | \(R(0)\) | 0 | Starting recovered or immune count |
| Transmission Rate | \(\beta\) | 0.30 | Effective contact rate per person per day |
| Recovery Rate | \(\gamma\) | 0.10 | Daily recovery rate (mean infectious duration \(1/\gamma = 10\) days) |
| Simulation Days | \(T\) | 90 | Duration of the modeling horizon in days, a whole number from 1 to 90 |
The model assumes homogeneous mixing across the population, constant rates during the horizon, and durable immunity following recovery.
Notebook implementation
The Python for Excel notebook reads parameters from the workbook, integrates the SIR ODE system, and publishes summary KPIs and the trajectory table:
inputs = bf.inputs(
parameters=bf.ref("Simulation!A4:B10", headers=True),
)
solution = solve_ivp(...)
publication = bf.publish(outputs={"summary": summary, "trajectory": trajectory})How the calculation/model works
The Kermack-McKendrick SIR model is formulated as three coupled non-linear ordinary differential equations:
\[ \frac{dS}{dt} = -\frac{\beta S I}{N} \] \[ \frac{dI}{dt} = \frac{\beta S I}{N} - \gamma I \] \[ \frac{dR}{dt} = \gamma I \]
The basic reproduction number (\(R_0\)) defines the expected number of secondary infections generated by a single infected individual in a fully susceptible population:
\[ R_0 = \frac{\beta}{\gamma} \]
When \(R_0 > 1\), an epidemic occurs. When \(R_0 \le 1\), no epidemic occurs and active infections decline from the start. The peak of active infections (\(dI/dt = 0\)) occurs precisely when the susceptible population declines to:
\[ S_{\text{peak}} = \frac{N}{R_0} = \frac{\gamma N}{\beta} \]
The corresponding theoretical herd immunity threshold is:
\[ \text{HIT} = 1 - \frac{1}{R_0} \]
Using solve_ivp(..., method='RK45'), the notebook computes values for each day from 0 to the simulation horizon \(T\) (at most 90), producing both a tabular daily projection and smooth phase portrait visualizations.
Architecture
- Fixed trajectory spill contract: The
trajectoryoutput spill occupies a fixed 91-row range (Simulation!A14:F105, header plus 91 daily projection rows from day 0 to 90). - Horizon capacity guard & padding:
notebook.pycapsSimulation Daysat 90 with a runtime guard (if duration > 90: raise ValueError(...)) to prevent spill errors. Shorter simulation horizons are padded withNonerows (range(90 - duration)) so obsolete rows are cleanly wiped from the worksheet. - Continuous turning-point root finding: Beyond daily evaluations, the notebook locates peak active infection timing (\(dI/dt = 0\)) using continuous ODE interpolation (
dense_output=True) via Brent’s method (scipy.optimize.root_scalar), decoupling peak metrics from integer-day discretization.
Validation / expected results
For the canonical baseline parameters:
| Check | Expected outcome |
|---|---|
| Basic Reproduction Number (\(R_0\)) | 3.0000 |
| Herd Immunity Threshold (fraction) | 0.6667 |
| Peak Infected Population | 30,079.6 |
| Peak Infection Day | 38.4 |
| Final Susceptible (\(t=90\)) | 6,149.8 |
| Final Infected (\(t=90\)) | 925.2 |
| Final Recovered (\(t=90\)) | 92,924.9 |
| Final Attack Rate (fraction) | 0.9292 |
A scenario lowering the transmission rate \(\beta\) to 0.20 (simulating distancing interventions or masking) reduces \(R_0\) to 2.00, lowers the herd immunity threshold to 0.5000 (50.0%), delays the peak to day 67.6, and flattens peak infections to 15,392.7.
Limitations
The standard SIR model assumes a closed and homogeneous population with well-mixed contact patterns and no age stratification, seasonal variation, or vital dynamics (births and deaths). Diseases with latency periods (SEIR models), waning immunity (SIRS), or asymmetric contact networks require specialized compartmental structures.
When to use this approach
Use scipy.integrate.solve_ivp whenever analyzing continuous epidemiological dynamics, viral spread curves, customer churn and adoption diffusion, or biological rate systems that cannot be captured by static linear formulas.