Programming the SEIR model with Euler’s method

The SEIR model introduces an exposed compartment. A newly infected person enters \(E\) before becoming infectious, allowing the model to represent a delay between infection and infectiousness.

Biological scenario

A closed population contains 1,000 people. Initially, 990 are susceptible, 5 are exposed, 5 are infectious and none have recovered.

We use \(\beta=0.3\) per day, \(\sigma=0.2\) per day and \(\gamma=0.1\) per day. We want to simulate the epidemic for 200 days and compare the exposed and infectious peaks.

The new exposed compartment

SSusceptible

→infection

EExposed

→progression

IInfectious

→recovery

RRecovered

In this basic SEIR model, exposed people are infected but are not yet infectious. They progress from \(E\) to \(I\) at rate \(\sigma E\).

Three biological flows

Infection flow\[\beta\frac{SI}{N}\]

Progression flow\[\sigma E\]

Recovery flow\[\gamma I\]

The progression parameter is related to the mean exposed period:

\[\text{mean exposed period}=\frac{1}{\sigma}.\]

For \(\sigma=0.2\) per day, the mean exposed period is \(1/0.2=5\) days.

Construct the four equations

For each compartment, use “flow in minus flow out.”

Susceptible\[\frac{dS}{dt}=-\beta\frac{SI}{N}\]

Exposed\[\frac{dE}{dt}=\beta\frac{SI}{N}-\sigma E\]

Infectious\[\frac{dI}{dt}=\sigma E-\gamma I\]

Recovered\[\frac{dR}{dt}=\gamma I\]

Compartment	Flow in	Flow out	Rate of change
\(S\)	None	Infection	− infection flow
\(E\)	Infection	Progression	infection − progression
\(I\)	Progression	Recovery	progression − recovery
\(R\)	Recovery	None	+ recovery flow

SIR and SEIR are different

SIR	SEIR
Newly infected people enter \(I\) immediately.	Newly infected people first enter \(E\).
The infection flow appears directly in \(dI/dt\).	The infection flow appears in \(dE/dt\); progression enters \(dI/dt\).
No explicit delay before infectiousness.	The exposed period creates a delay before infectiousness.
Three compartment arrays.	Four compartment arrays.

Step 1: Write the SEIR rate function

def seir_rates(S, E, I, R, beta, sigma, gamma):
    N = S + E + I + R

    infection_flow = beta * S * I / N
    progression_flow = sigma * E
    recovery_flow = gamma * I

    dS_dt = -infection_flow
    dE_dt = infection_flow - progression_flow
    dI_dt = progression_flow - recovery_flow
    dR_dt = recovery_flow

    return dS_dt, dE_dt, dI_dt, dR_dt

The function returns four derivatives in the same order as the compartments.

Step 2: Add exposed storage

S = np.zeros(number_of_steps + 1)
E = np.zeros(number_of_steps + 1)
I = np.zeros(number_of_steps + 1)
R = np.zeros(number_of_steps + 1)

S[0], E[0], I[0], R[0] = S0, E0, I0, R0

The arrays have matching lengths. Index \(n\) represents one common epidemic state \((S_n,E_n,I_n,R_n)\).

Step 3: Update all four compartments

\[S_{n+1}=S_n+\left(\frac{dS}{dt}\right)_n\Delta t\]

\[E_{n+1}=E_n+\left(\frac{dE}{dt}\right)_n\Delta t\]

\[I_{n+1}=I_n+\left(\frac{dI}{dt}\right)_n\Delta t\]

\[R_{n+1}=R_n+\left(\frac{dR}{dt}\right)_n\Delta t\]

As in the SIR program, calculate all derivatives from the same current state before storing any next value.

Run the complete SEIR simulation

Interactive PythonComplete SEIR Euler simulation

import numpy as np
import matplotlib.pyplot as plt

def seir_rates(S, E, I, R, beta, sigma, gamma):
    """Return the deterministic SEIR rates at one state."""
    
    N = S + E + I + R
    
    infection_flow = beta * S * I / N
    progression_flow = sigma * E
    recovery_flow = gamma * I
    
    dS_dt = -infection_flow
    dE_dt = infection_flow - progression_flow
    dI_dt = progression_flow - recovery_flow
    dR_dt = recovery_flow
    
    return dS_dt, dE_dt, dI_dt, dR_dt

# Initial conditions
S0 = 990
E0 = 5
I0 = 5
R0 = 0
N = S0 + E0 + I0 + R0

# Model parameters
beta = 0.3
sigma = 0.2
gamma = 0.1

# Numerical settings
end_time = 200
dt = 0.1

# Validate settings
if min(S0, E0, I0, R0, beta, sigma, gamma) < 0:
    raise ValueError("Compartments and parameters cannot be negative.")
if N <= 0:
    raise ValueError("The population must be positive.")
if dt <= 0:
    raise ValueError("dt must be positive.")

# Prepare time and storage
number_of_steps = int(round(end_time / dt))
time = np.linspace(0, end_time, number_of_steps + 1)

S = np.zeros(number_of_steps + 1)
E = np.zeros(number_of_steps + 1)
I = np.zeros(number_of_steps + 1)
R = np.zeros(number_of_steps + 1)

S[0], E[0], I[0], R[0] = S0, E0, I0, R0

# Forward Euler simulation
for n in range(number_of_steps):
    dS_dt, dE_dt, dI_dt, dR_dt = seir_rates(
        S[n], E[n], I[n], R[n],
        beta, sigma, gamma
    )
    
    S[n + 1] = S[n] + dS_dt * dt
    E[n + 1] = E[n] + dE_dt * dt
    I[n + 1] = I[n] + dI_dt * dt
    R[n + 1] = R[n] + dR_dt * dt
    
    if min(S[n + 1], E[n + 1], I[n + 1], R[n + 1]) < -1e-10:
        raise ValueError("A negative compartment occurred. Reduce dt.")

# Calculate exposed and infectious peaks
exposed_peak_index = np.argmax(E)
infectious_peak_index = np.argmax(I)

exposed_peak = E[exposed_peak_index]
infectious_peak = I[infectious_peak_index]

exposed_peak_time = time[exposed_peak_index]
infectious_peak_time = time[infectious_peak_index]
peak_delay = infectious_peak_time - exposed_peak_time

# Check population conservation
population = S + E + I + R
conservation_error = np.max(np.abs(population - N))

# Display summary
print("SEIR epidemic summary")
print("-----------------------------------------")
print(f"Peak exposed:          {exposed_peak:.2f}")
print(f"Exposed peak time:     {exposed_peak_time:.1f} days")
print(f"Peak infectious:       {infectious_peak:.2f}")
print(f"Infectious peak time:  {infectious_peak_time:.1f} days")
print(f"Delay between peaks:   {peak_delay:.1f} days")
print(f"Final recovered:       {R[-1]:.2f}")
print(f"Maximum total error:   {conservation_error:.2e}")

# Display selected days
print("\nDay       S         E         I         R")
print("------------------------------------------------")
for day in [0, 20, 40, 60, 80, 120, 160, 200]:
    index = int(round(day / dt))
    print(
        f"{day:3d}  {S[index]:8.2f}  {E[index]:8.2f}  "
        f"{I[index]:8.2f}  {R[index]:8.2f}"
    )

# Plot the epidemic
plt.figure(figsize=(8, 4.8))
plt.plot(time, S, label="Susceptible")
plt.plot(time, E, label="Exposed")
plt.plot(time, I, label="Infectious")
plt.plot(time, R, label="Recovered")
plt.axvline(
    infectious_peak_time,
    color="grey",
    linestyle=":",
    label=f"Infectious peak: day {infectious_peak_time:.1f}"
)
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("Deterministic SEIR model: forward Euler method")
plt.grid(alpha=0.25)
plt.legend()
plt.tight_layout()

Output

Run the code to see the result.

New outcomes in the SEIR model

Exposed peaklargest value reached by \(E(t)\)

Infectious peaklargest value reached by \(I(t)\)

Peak delayinfectious peak time minus exposed peak time

Code	Meaning
`np.argmax(E)`	Index at which the exposed population is largest
`np.argmax(I)`	Index at which the infectious population is largest
`time[infectious_peak_index]`	Time corresponding to the infectious peak
`infectious_peak_time - exposed_peak_time`	Numerical delay between the two peaks

Biological interpretation

New infections increase the exposed population first. Only after progression do they increase the infectious population. Consequently, the exposed curve usually responds earlier, while the infectious curve is delayed and smoothed by the exposed period.

The parameter \(\sigma\) controls this delay. Increasing \(\sigma\) shortens the mean exposed period and moves people into \(I\) more rapidly. Decreasing \(\sigma\) lengthens the delay.

Important model and numerical conditions

The basic model assumes exposed people do not transmit infection.
All initial compartments and parameters must be non-negative.
The closed population satisfies \(S+E+I+R=N\).
\(\Delta t\) must be positive and sufficiently small.
All four Euler derivatives must be evaluated at the same current state.
The simulation must be long enough to include both peaks and the epidemic decline.
The interpretation \(1/\sigma\) assumes a constant per-person progression rate.

Change the exposed period

Run the model with sigma = 0.5 and then with sigma = 0.1.

Calculate the corresponding mean exposed periods.
Compare the infectious peak times.
Compare the delay between the exposed and infectious peaks.
Explain the graphical differences biologically.