Why nonlinear models produce unclosed moments

Nonlinear biological interactions make lower-moment equations depend on higher moments, creating an infinite hierarchy.

Nonlinear SIS infection rate

\[b(I)=\beta\frac{(N-I)I}{N}=\beta I-\frac{\beta}{N}I^2,\qquad d(I)=\gamma I.\]

The product \((N-I)I\) represents susceptible–infectious contact and makes the rate nonlinear.

First moment requires the second moment

\[ \frac{dm_1}{dt} = (\beta-\gamma)m_1-\frac{\beta}{N}m_2. \]

The first-moment equation is not closed because knowing \(m_1\) alone does not determine \(m_2\).

Second moment requires the third moment

\[ \frac{dm_2}{dt} = \left[2(\beta-\gamma)-\frac{\beta}{N}\right]m_2 -\frac{2\beta}{N}m_3 +(\beta+\gamma)m_1. \]

Deriving \(m_3'\) introduces \(m_4\), and so on. This is the moment hierarchy.

Why replacing \(\mathbb E[I^2]\) by \(\mathbb E[I]^2\) is an assumption

\[\mathbb E[I^2]=\operatorname{Var}(I)+\mathbb E[I]^2.\]

Therefore \(\mathbb E[I^2]=\mathbb E[I]^2\) holds only when variance is zero. Using it in a stochastic epidemic discards distributional spread.

Interactive Python laboratory

Interactive PythonDemonstrating the closure gap

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

N, I0 = 100, 5
beta, gamma = 0.25, 0.15
grid = np.linspace(0, 30, 61)
M = 4000
rng = np.random.default_rng(270)

def path_on_grid():
    I = I0
    time = 0.0
    event_times = [0.0]
    values = [I]
    while I > 0 and time < grid[-1]:
        b = beta * (N - I) * I / N
        d = gamma * I
        wait = rng.exponential(1 / (b + d))
        if time + wait > grid[-1]:
            break
        time += wait
        I += 1 if rng.random() < b / (b + d) else -1
        event_times.append(time)
        values.append(I)
    positions = np.searchsorted(
        event_times, grid, side="right"
    ) - 1
    return np.array(values)[positions]

paths = np.array([path_on_grid() for _ in range(M)])
m1 = paths.mean(axis=0)
m2 = np.mean(paths ** 2, axis=0)
variance = paths.var(axis=0, ddof=0)

exact_first_derivative = (
    (beta - gamma) * m1 - beta * m2 / N
)
mean_field_derivative = (
    (beta - gamma) * m1 - beta * m1 ** 2 / N
)

selected = [0, 10, 20, 30]
idx = [np.argmin(np.abs(grid - t)) for t in selected]
table = pd.DataFrame({
    "time": selected,
    "m1": m1[idx],
    "m2": m2[idx],
    "m1 squared": (m1 ** 2)[idx],
    "variance": variance[idx],
    "exact moment drift": exact_first_derivative[idx],
    "mean-field drift": mean_field_derivative[idx]
})
print(table.round(3).to_string(index=False))

plt.figure(figsize=(9.5, 5.1))
plt.plot(grid, m2, label="m2 = E[I²]",
         color="#b42318")
plt.plot(grid, m1 ** 2, "--",
         label="m1² = E[I]²", color="#174c34")
plt.fill_between(grid, m1 ** 2, m2,
                 color="#f2c7c2", alpha=0.55,
                 label="Difference = variance")
plt.xlabel("Time")
plt.ylabel("Moment value")
plt.title("Non-zero variance creates a closure gap")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Exact hierarchy versus approximation

The unclosed equations are exact statements about the CTMC moments. The difficulty is that a finite ODE system cannot contain infinitely many moments. Closing the system requires an additional distributional or factorisation assumption, introduced next.

What this lesson adds

You can now identify how nonlinear rates raise moment order, write the first two SIS moment equations, recognise the infinite hierarchy and explain why factorising a higher moment is an approximation.