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 987.151 2.5789 9.30939 0.960608
2 984.437 4.5703 9.1143 1.87832
3 981.732 6.19384 9.27836 2.79539
4 978.946 7.59875 9.71272 3.74297
5 976.007 8.88733 10.3605 4.74502
6 972.862 10.1306 11.1866 5.821
7 969.463 11.379 12.1706 6.9876
8 965.769 12.6688 13.302 8.26003
9 961.742 14.0274 14.5775 9.6528
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 987.151 2.5789 9.30939 0.960608
2 984.437 4.5703 9.1143 1.87832
3 981.732 6.19384 9.27836 2.79539
4 978.946 7.59875 9.71271 3.74297
5 976.007 8.88734 10.3605 4.74503
6 972.862 10.1307 11.1866 5.821
7 969.463 11.379 12.1706 6.9876
8 965.769 12.6688 13.302 8.26003
9 961.742 14.0274 14.5775 9.6528
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 976.007 8.88733 10.3605 4.74502
10 957.345 15.476 15.9986 11.1804
15 928.483 24.5936 25.5388 21.3841
20 884.488 37.7811 40.1657 37.5654
25 820.36 55.5943 61.3916 62.6537
30 732.997 76.8898 89.9255 100.188
35 624.565 97.7305 124.155 153.55
40 504.911 111.785 158.863 224.441
45 389.268 113.599 185.988 311.145
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 935.686 27.1893 31.052 6.07294
2 920.717 33.0956 33.6638 12.5238
3 904.554 38.5115 37.3265 19.608
4 886.867 43.8462 41.7801 27.5069
5 867.445 49.3 46.8908 36.3637
6 846.159 54.9541 52.5852 46.302
7 822.936 60.8162 58.8145 57.4335
8 797.762 66.844 65.5332 69.8605
9 770.679 72.9576 72.6872 83.6758
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.
    """
    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)

      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}."

      t_eval = np.linspace(t0, t1, ntp + 1)

      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]

      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
      )
      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)):
        row = [float(sol.t[i]), float(sol.y[0][i]), float(sol.y[1][i]), float(sol.y[2][i]), float(sol.y[3][i])]
        if not all(np.isfinite(v) for v in row):
          return "Error: Invalid output: nan or inf detected."
        result.append(row)
      return result
    except Exception as e:
      return f"Error: {str(e)}"

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.