SEIR

Overview

The SEIR function numerically solves the SEIR compartmental model, a system of ordinary differential equations used in epidemiology to simulate the spread of infectious diseases with a latent incubation period. The model extends the classic SIR model by adding an Exposed (E) compartment for individuals who have been infected but are not yet infectious themselves.

The SEIR model divides a population into four compartments representing disease progression:

  • S (Susceptible): Individuals who can contract the disease
  • E (Exposed): Infected individuals in the latent period who are not yet infectious
  • I (Infectious): Individuals actively spreading the disease
  • R (Recovered): Individuals who have recovered and gained immunity

The model is governed by the following system of differential equations:

\begin{aligned} \frac{dS}{dt} &= -\beta \frac{SI}{N} \\ \frac{dE}{dt} &= \beta \frac{SI}{N} - \sigma E \\ \frac{dI}{dt} &= \sigma E - \gamma I \\ \frac{dR}{dt} &= \gamma I \end{aligned}

where N = S + E + I + R is the total population, \beta is the transmission rate, \sigma is the incubation rate (inverse of the average latent period), and \gamma is the recovery rate (inverse of the average infectious period).

This implementation uses SciPy’s solve_ivp function for numerical integration, which provides multiple solver methods including Runge-Kutta (RK45, RK23, DOP853) and implicit methods (Radau, BDF, LSODA). The SEIR model is particularly useful for diseases with significant incubation periods such as measles, influenza, and COVID-19, where the delay between infection and infectiousness meaningfully affects outbreak dynamics.

The basic reproduction number R_0 for the SEIR model is given by:

R_0 = \frac{\beta}{\gamma}

When R_0 > 1, an epidemic outbreak occurs; when R_0 \leq 1, the disease will eventually die out. The original formulation of compartmental epidemic models was developed by Kermack and McKendrick in 1927 and remains foundational to modern epidemiological modeling.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=SEIR(s_initial, e_initial, i_initial, r_initial, beta, gamma, sigma, t_start, t_end, timesteps, solve_ivp_method)
  • s_initial (float, required): Initial number of susceptible individuals.
  • e_initial (float, required): Initial number of exposed individuals.
  • i_initial (float, required): Initial number of infected individuals.
  • r_initial (float, required): Initial number of recovered individuals.
  • beta (float, required): Transmission rate (contacts per time).
  • gamma (float, required): Recovery rate (1/time).
  • sigma (float, required): Incubation rate (1/time).
  • t_start (float, required): Start time of integration.
  • t_end (float, required): End time of integration.
  • timesteps (int, optional, default: 10): Number of timesteps to evaluate.
  • solve_ivp_method (str, optional, default: “RK45”): Integration method for solve_ivp.

Returns (list[list]): 2D list with header [t, S, E, I, R], or error message string.

Example 1: SEIR epidemic with incubation period over 10 days

Inputs:

s_initial e_initial i_initial r_initial beta gamma sigma t_start t_end timesteps
990 0 10 0 0.3 0.1 0.2 0 10 10

Excel formula:

=SEIR(990, 0, 10, 0, 0.3, 0.1, 0.2, 0, 10, 10)

Expected output:

t S E I R
0 990 0 10 0
1.11111 986.846 2.8244 9.26625 1.0638
2.22222 983.839 4.95671 9.12347 2.08093
3.33333 980.816 6.68059 9.39646 3.10656
4.44444 977.662 8.181 9.97678 4.18036
5.55556 974.288 9.58005 10.7991 5.33255
6.66667 970.627 10.9598 11.8259 6.58771
7.77778 966.617 12.3771 13.0381 7.96737
8.88889 962.207 13.8723 14.4286 9.49166
10 957.345 15.476 15.9986 11.1804
Example 2: SEIR integration using RK23 method

Inputs:

s_initial e_initial i_initial r_initial beta gamma sigma t_start t_end timesteps solve_ivp_method
990 0 10 0 0.3 0.1 0.2 0 10 10 RK23

Excel formula:

=SEIR(990, 0, 10, 0, 0.3, 0.1, 0.2, 0, 10, 10, "RK23")

Expected output:

t S E I R
0 990 0 10 0
1.11111 986.846 2.8244 9.26625 1.0638
2.22222 983.839 4.95672 9.12346 2.08093
3.33333 980.816 6.6806 9.39646 3.10656
4.44444 977.662 8.18101 9.97678 4.18036
5.55556 974.288 9.58006 10.7991 5.33255
6.66667 970.627 10.9599 11.8259 6.58771
7.77778 966.617 12.3771 13.0381 7.96737
8.88889 962.207 13.8724 14.4286 9.49165
10 957.345 15.476 15.9986 11.1804
Example 3: Full epidemic curve over 50 days

Inputs:

s_initial e_initial i_initial r_initial beta gamma sigma t_start t_end timesteps
990 0 10 0 0.3 0.1 0.2 0 50 10

Excel formula:

=SEIR(990, 0, 10, 0, 0.3, 0.1, 0.2, 0, 50, 10)

Expected output:

t S E I R
0 990 0 10 0
5.55556 974.288 9.58005 10.7991 5.33255
11.1111 951.98 17.2121 17.7546 13.0538
16.6667 915.756 28.4877 29.7715 25.9851
22.2222 858.735 45.1413 48.7089 47.4149
27.7778 774.719 67.1748 76.3711 81.7347
33.3333 662.602 91.1995 112.355 133.844
38.8889 531.721 109.582 151.502 207.195
44.4444 401.445 114.072 183.605 300.878
50 291.178 102.642 198.257 407.923
Example 4: Higher transmission rate with initial exposed population

Inputs:

s_initial e_initial i_initial r_initial beta gamma sigma t_start t_end timesteps
950 20 30 0 0.5 0.2 0.3 0 10 10

Excel formula:

=SEIR(950, 20, 30, 0, 0.5, 0.2, 0.3, 0, 10, 10)

Expected output:

t S E I R
0 950 20 30 0
1.11111 934.067 27.8898 31.2776 6.76546
2.22222 917.242 34.3227 34.3988 14.0361
3.33333 898.84 40.2857 38.7316 22.1428
4.44444 878.459 46.2485 43.9757 31.3173
5.55556 855.856 52.4146 49.9857 41.7441
6.66667 830.894 58.8408 56.6813 53.5839
7.77778 803.524 65.4937 64.0004 66.982
8.88889 773.78 72.2772 71.8728 82.0696
10 741.787 79.0468 80.2061 98.9597

Python Code

Show Code 
from scipy.integrate import solve_ivp
import numpy as np

def seir(s_initial, e_initial, i_initial, r_initial, beta, gamma, sigma, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the SEIR system of ordinary differential equations for infectious disease modeling using scipy.integrate.solve_ivp.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        s_initial (float): Initial number of susceptible individuals.
        e_initial (float): Initial number of exposed individuals.
        i_initial (float): Initial number of infected individuals.
        r_initial (float): Initial number of recovered individuals.
        beta (float): Transmission rate (contacts per time).
        gamma (float): Recovery rate (1/time).
        sigma (float): Incubation rate (1/time).
        t_start (float): Start time of integration.
        t_end (float): End time of integration.
        timesteps (int, optional): Number of timesteps to evaluate. Default is 10.
        solve_ivp_method (str, optional): Integration method for solve_ivp. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, S, E, I, R], or error message string.
    """
    # Validate inputs
    try:
        s0 = float(s_initial)
        e0 = float(e_initial)
        i0 = float(i_initial)
        r0 = float(r_initial)
        beta_ = float(beta)
        gamma_ = float(gamma)
        sigma_ = float(sigma)
        t0 = float(t_start)
        t1 = float(t_end)
        ntp = int(timesteps)
    except Exception:
        return "Error: Invalid input: all initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Error: Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Error: Invalid input: timesteps must be a positive integer."
    if any(x < 0 for x in [s0, e0, i0, r0, beta_, gamma_, sigma_]):
        return "Error: Invalid input: all initial values and parameters must be non-negative."
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in allowed_methods:
        return f"Error: Invalid input: solve_ivp_method must be one of {allowed_methods}."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)

    N = s0 + e0 + i0 + r0
    def seir_ode(t, y):
        S, E, I, R = y
        dSdt = -beta_ * S * I / N
        dEdt = beta_ * S * I / N - sigma_ * E
        dIdt = sigma_ * E - gamma_ * I
        dRdt = gamma_ * I
        return [dSdt, dEdt, dIdt, dRdt]
    try:
        sol = solve_ivp(
            seir_ode,
            [t0, t1],
            [s0, e0, i0, r0],
            method=solve_ivp_method,
            dense_output=False,
            rtol=1e-6,
            atol=1e-8,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Error: ODE integration error: {e}"
    if not sol.success:
        return f"Error: ODE integration failed: {sol.message}"
    result = [["t", "S", "E", "I", "R"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        s_val = float(sol.y[0][i])
        e_val = float(sol.y[1][i])
        i_val = float(sol.y[2][i])
        r_val = float(sol.y[3][i])
        # Disallow nan/inf
        if any([
            not isinstance(t_val, float) or not isinstance(s_val, float) or not isinstance(e_val, float) or not isinstance(i_val, float) or not isinstance(r_val, float),
            t_val != t_val or s_val != s_val or e_val != e_val or i_val != i_val or r_val != r_val,
            abs(t_val) == float('inf') or abs(s_val) == float('inf') or abs(e_val) == float('inf') or abs(i_val) == float('inf') or abs(r_val) == float('inf')
        ]):
            return "Error: Invalid output: nan or inf detected."
        result.append([t_val, s_val, e_val, i_val, r_val])
    return result

Online Calculator

Initial number of susceptible individuals.
Initial number of exposed individuals.
Initial number of infected individuals.
Initial number of recovered individuals.
Transmission rate (contacts per time).
Recovery rate (1/time).
Incubation rate (1/time).
Start time of integration.
End time of integration.
Number of timesteps to evaluate.
Integration method for solve_ivp.