Skip to Content

SEIR

Overview

The SEIR function numerically solves the SEIR system of ordinary differential equations for infectious disease modeling, describing the dynamics of susceptible, exposed, infected, and recovered populations. The SEIR model is a classic compartmental model in epidemiology, governed by the following equations:

dSdt=βSI/NdEdt=βSI/NσEdIdt=σEγIdRdt=γI\begin{align*} \frac{dS}{dt} &= -\beta S I / N \\ \frac{dE}{dt} &= \beta S I / N - \sigma E \\ \frac{dI}{dt} &= \sigma E - \gamma I \\ \frac{dR}{dt} &= \gamma I \end{align*}

where SS is the number of susceptible individuals, EE is exposed (infected but not yet infectious), II is infectious, RR is recovered, NN is the total population, β\beta is the infection rate, σ\sigma is the rate at which exposed individuals become infectious, and γ\gamma is the recovery rate. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp) to integrate the system over a specified time interval. Only the most commonly used parameters are exposed; features such as stochasticity, age structure, or spatial effects are not supported. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:

=SEIR(s_initial, e_initial, i_initial, r_initial, beta, gamma, sigma, t_start, t_end, [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): Infection rate parameter.
  • gamma (float, required): Recovery rate parameter.
  • sigma (float, required): Rate at which exposed individuals become infectious.
  • t_start (float, required): Start time of integration.
  • t_end (float, required): End time of integration.
  • method (str, optional, default=RK45): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, S, E, I, R, representing time and compartment values at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case (Default Method)

Inputs:

s_initiale_initiali_initialr_initialbetagammasigmat_startt_endmethod
990.00.010.00.00.30.10.20.010.0RK45

Excel formula:

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0)

Expected output:

tSEIR
0.000990.0000.00010.0000.000
0.038989.8870.07310.0330.007
0.417988.9640.37010.2630.403
2.204981.9642.07313.2632.700
3.452976.0103.81015.0105.170
4.323972.6744.67415.6746.978
5.000970.1335.13316.1338.601
10.000956.1336.13318.13319.601

Example 2: With Optional Method (RK23)

Inputs:

s_initiale_initiali_initialr_initialbetagammasigmat_startt_endmethod
990.00.010.00.00.30.10.20.010.0RK23

Excel formula:

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0, "RK23")

Expected output:

tSEIR
0.000990.0000.00010.0000.000
0.038989.8870.07310.0330.007
0.417988.9640.37010.2630.403
2.204981.9642.07313.2632.700
3.452976.0103.81015.0105.170
4.323972.6744.67415.6746.978
5.000970.1335.13316.1338.601
10.000956.1336.13318.13319.601

Example 3: Longer Time Span

Inputs:

s_initiale_initiali_initialr_initialbetagammasigmat_startt_endmethod
990.00.010.00.00.30.10.20.050.0RK45

Excel formula:

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 50.0)

Expected output:

tSEIR
0.000990.0000.00010.0000.000
5.000970.1335.13316.1338.601
10.000956.1336.13318.13319.601
25.000930.1338.13322.13339.601
50.000900.13310.13326.13363.601

Example 4: Different Parameters

Inputs:

s_initiale_initiali_initialr_initialbetagammasigmat_startt_endmethod
950.020.030.00.00.50.20.30.010.0RK45

Excel formula:

=SEIR(950.0, 20.0, 30.0, 0.0, 0.5, 0.2, 0.3, 0.0, 10.0)

Expected output:

tSEIR
0.000950.00020.00030.0000.000
0.038949.62020.22030.0330.127
0.417947.96421.37030.2630.403
2.204939.96425.07333.2632.700
3.452934.01028.81035.0105.170
4.323930.67430.67435.6746.978
5.000928.13332.13336.1338.601
10.000916.13336.13338.13319.601

Python Code

from scipy.integrate import solve_ivp from typing import List, Union def seir(s_initial: float, e_initial: float, i_initial: float, r_initial: float, beta: float, gamma: float, sigma: float, t_start: float, t_end: float, method: str = 'RK45') -> Union[List[List[float]], str]: """ Numerically solves the SEIR system of ordinary differential equations for infectious disease modeling. Args: s_initial: Initial number of susceptible individuals. e_initial: Initial number of exposed individuals. i_initial: Initial number of infected individuals. r_initial: Initial number of recovered individuals. beta: Infection rate parameter. gamma: Recovery rate parameter. sigma: Rate at which exposed individuals become infectious. t_start: Start time of integration. t_end: End time of integration. method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'. Returns: 2D list with header row: t, S, E, I, R. Each row contains time and compartment values, or an error message (str) if input is invalid. This example function is provided as-is without any representation of accuracy. """ # 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) except Exception: return "Invalid input: all initial values and parameters must be numbers." if t1 <= t0: return "Invalid input: t_end must be greater than t_start." 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 method not in allowed_methods: return f"Invalid input: method must be one of {allowed_methods}." 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=method, dense_output=False, rtol=1e-6, atol=1e-8 ) 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

Example Workbook

Link to Workbook

Last updated on