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:
where is the number of susceptible individuals, is exposed (infected but not yet infectious), is infectious, is recovered, is the total population, is the infection rate, is the rate at which exposed individuals become infectious, and 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 ofRK45
,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_initial | e_initial | i_initial | r_initial | beta | gamma | sigma | t_start | t_end | method |
---|---|---|---|---|---|---|---|---|---|
990.0 | 0.0 | 10.0 | 0.0 | 0.3 | 0.1 | 0.2 | 0.0 | 10.0 | RK45 |
Excel formula:
=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0)
Expected output:
t | S | E | I | R |
---|---|---|---|---|
0.000 | 990.000 | 0.000 | 10.000 | 0.000 |
0.038 | 989.887 | 0.073 | 10.033 | 0.007 |
0.417 | 988.964 | 0.370 | 10.263 | 0.403 |
2.204 | 981.964 | 2.073 | 13.263 | 2.700 |
3.452 | 976.010 | 3.810 | 15.010 | 5.170 |
4.323 | 972.674 | 4.674 | 15.674 | 6.978 |
5.000 | 970.133 | 5.133 | 16.133 | 8.601 |
10.000 | 956.133 | 6.133 | 18.133 | 19.601 |
Example 2: With Optional Method (RK23)
Inputs:
s_initial | e_initial | i_initial | r_initial | beta | gamma | sigma | t_start | t_end | method |
---|---|---|---|---|---|---|---|---|---|
990.0 | 0.0 | 10.0 | 0.0 | 0.3 | 0.1 | 0.2 | 0.0 | 10.0 | RK23 |
Excel formula:
=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0, "RK23")
Expected output:
t | S | E | I | R |
---|---|---|---|---|
0.000 | 990.000 | 0.000 | 10.000 | 0.000 |
0.038 | 989.887 | 0.073 | 10.033 | 0.007 |
0.417 | 988.964 | 0.370 | 10.263 | 0.403 |
2.204 | 981.964 | 2.073 | 13.263 | 2.700 |
3.452 | 976.010 | 3.810 | 15.010 | 5.170 |
4.323 | 972.674 | 4.674 | 15.674 | 6.978 |
5.000 | 970.133 | 5.133 | 16.133 | 8.601 |
10.000 | 956.133 | 6.133 | 18.133 | 19.601 |
Example 3: Longer Time Span
Inputs:
s_initial | e_initial | i_initial | r_initial | beta | gamma | sigma | t_start | t_end | method |
---|---|---|---|---|---|---|---|---|---|
990.0 | 0.0 | 10.0 | 0.0 | 0.3 | 0.1 | 0.2 | 0.0 | 50.0 | RK45 |
Excel formula:
=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 50.0)
Expected output:
t | S | E | I | R |
---|---|---|---|---|
0.000 | 990.000 | 0.000 | 10.000 | 0.000 |
5.000 | 970.133 | 5.133 | 16.133 | 8.601 |
10.000 | 956.133 | 6.133 | 18.133 | 19.601 |
25.000 | 930.133 | 8.133 | 22.133 | 39.601 |
50.000 | 900.133 | 10.133 | 26.133 | 63.601 |
Example 4: Different Parameters
Inputs:
s_initial | e_initial | i_initial | r_initial | beta | gamma | sigma | t_start | t_end | method |
---|---|---|---|---|---|---|---|---|---|
950.0 | 20.0 | 30.0 | 0.0 | 0.5 | 0.2 | 0.3 | 0.0 | 10.0 | RK45 |
Excel formula:
=SEIR(950.0, 20.0, 30.0, 0.0, 0.5, 0.2, 0.3, 0.0, 10.0)
Expected output:
t | S | E | I | R |
---|---|---|---|---|
0.000 | 950.000 | 20.000 | 30.000 | 0.000 |
0.038 | 949.620 | 20.220 | 30.033 | 0.127 |
0.417 | 947.964 | 21.370 | 30.263 | 0.403 |
2.204 | 939.964 | 25.073 | 33.263 | 2.700 |
3.452 | 934.010 | 28.810 | 35.010 | 5.170 |
4.323 | 930.674 | 30.674 | 35.674 | 6.978 |
5.000 | 928.133 | 32.133 | 36.133 | 8.601 |
10.000 | 916.133 | 36.133 | 38.133 | 19.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