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.

Examples

Example 1: Demo case 1

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.111	986.8	2.824	9.266	1.064
2.222	983.8	4.957	9.123	2.081
3.333	980.8	6.681	9.396	3.107
4.444	977.7	8.181	9.977	4.18
5.556	974.3	9.58	10.8	5.333
6.667	970.6	10.96	11.83	6.588
7.778	966.6	12.38	13.04	7.967
8.889	962.2	13.87	14.43	9.492
10	957.3	15.48	16	11.18

Example 2: Demo case 2

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.111	986.8	2.824	9.266	1.064
2.222	983.8	4.957	9.123	2.081
3.333	980.8	6.681	9.396	3.107
4.444	977.7	8.181	9.977	4.18
5.556	974.3	9.58	10.8	5.333
6.667	970.6	10.96	11.83	6.588
7.778	966.6	12.38	13.04	7.967
8.889	962.2	13.87	14.43	9.492
10	957.3	15.48	16	11.18

Example 3: Demo case 3

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.556	974.3	9.58	10.8	5.333
11.11	952	17.21	17.75	13.05
16.67	915.8	28.49	29.77	25.99
22.22	858.7	45.14	48.71	47.41
27.78	774.7	67.17	76.37	81.73
33.33	662.6	91.2	112.4	133.8
38.89	531.7	109.6	151.5	207.2
44.44	401.4	114.1	183.6	300.9
50	291.2	102.6	198.3	407.9

Example 4: Demo case 4

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.111	934.1	27.89	31.28	6.765
2.222	917.2	34.32	34.39	14.04
3.333	898.8	40.29	38.73	22.14
4.444	878.5	46.25	43.98	31.32
5.556	855.9	52.41	50	41.74
6.667	830.9	58.84	56.68	53.58
7.778	803.5	65.49	64	66.98
8.889	773.8	72.28	71.87	82.07
10	741.8	79.05	80.21	98.96

Python 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 "Invalid input: all initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if any(x < 0 for x in [s0, e0, i0, r0, beta_, gamma_, sigma_]):
        return "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"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"ODE integration error: {e}"
    if not sol.success:
        return f"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 "Invalid output: nan or inf detected."
        result.append([t_val, s_val, e_val, i_val, r_val])
    return result

Overview

Excel Usage

Examples

Python Code

Online Calculator