Constructing a stochastic SIS model

A stochastic SIS SDE is not created by adding arbitrary noise. Its drift and diffusion are derived from the infection and recovery events of the underlying epidemic process.

Scenario and state

A closed SIS population contains \(N\) people. The state is the infectious count \(I(t)\), and \(S(t)=N-I(t)\). Recovery returns a person immediately to susceptibility.

Begin with biological events

Event	Change in \(I\)	Rate
Infection	+1	\(a_1(I)=\beta(N-I)I/N\)
Recovery	−1	\(a_2(I)=\gamma I\)

The SDE approximation must retain both the average net effect and the variability generated by these two event streams.

Construct the drift

\[f(I)=(+1)a_1(I)+(-1)a_2(I)=a_1(I)-a_2(I).\]

Therefore

\[f(I)=\beta\frac{(N-I)I}{N}-\gamma I.\]

Construct the diffusion variance

For independent event counts, variances add after multiplying by the squared state changes:

\[g^2(I)=(+1)^2a_1(I)+(-1)^2a_2(I)=a_1(I)+a_2(I).\]

\[g(I)=\sqrt{\beta\frac{(N-I)I}{N}+\gamma I}.\]

The recovery sign disappears after squaring. Recoveries reduce the drift but still increase uncertainty.

The stochastic SIS model

\[ dI=\left[\beta\frac{(N-I)I}{N}-\gamma I\right]dt +\sqrt{\beta\frac{(N-I)I}{N}+\gamma I}\,dW. \]

At \(I=0\), both coefficients vanish, matching the disease-free absorbing state. The SDE is a continuous approximation: \(I\) need not remain integer.

Interactive Python laboratory

Run one Euler–Maruyama SIS path beside its deterministic solution. A second graph shows how drift direction and diffusion strength change across the state space.

Interactive PythonStochastic SIS model

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

def infection_rate(I, N, beta):
    return beta * (N - I) * I / N

def recovery_rate(I, gamma):
    return gamma * I

def sis_drift(I, N, beta, gamma):
    return (
        infection_rate(I, N, beta)
        - recovery_rate(I, gamma)
    )

def sis_diffusion(I, N, beta, gamma):
    total_event_rate = (
        infection_rate(I, N, beta)
        + recovery_rate(I, gamma)
    )
    return np.sqrt(max(total_event_rate, 0))

N = 500
I0 = 10.0
beta = 0.25
gamma = 0.15
T = 100
dt = 0.01
seed = 44

if not 0 <= I0 <= N:
    raise ValueError("I0 must lie between 0 and N.")

if dt <= 0:
    raise ValueError("dt must be positive.")

steps = int(round(T / dt))
times = np.linspace(0, T, steps + 1)
rng = np.random.default_rng(seed)

I_sde = np.zeros(steps + 1)
I_ode = np.zeros(steps + 1)
I_sde[0] = I0
I_ode[0] = I0

boundary_corrections = 0

for n in range(steps):
    current_I = I_sde[n]

if current_I == 0:
        I_sde[n + 1] = 0
    else:
        dW = np.sqrt(dt) * rng.normal(0, 1)
        proposed_I = (
            current_I
            + sis_drift(
                current_I, N, beta, gamma
            ) * dt
            + sis_diffusion(
                current_I, N, beta, gamma
            ) * dW
        )

# Explicitly stated projection rule for this approximation
        I_sde[n + 1] = np.clip(proposed_I, 0, N)

if I_sde[n + 1] != proposed_I:
            boundary_corrections += 1

I_ode[n + 1] = (
        I_ode[n]
        + sis_drift(I_ode[n], N, beta, gamma) * dt
    )

summary = pd.DataFrame({
    "quantity": [
        "R0 = beta/gamma",
        "Final deterministic I",
        "Final stochastic I",
        "Minimum stochastic I",
        "Maximum stochastic I",
        "Boundary corrections"
    ],
    "value": [
        beta / gamma,
        I_ode[-1],
        I_sde[-1],
        I_sde.min(),
        I_sde.max(),
        boundary_corrections
    ]
})

print(summary.round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, I_ode, color="#174c34",
         linewidth=2.2, label="Deterministic SIS")
plt.plot(times, I_sde, color="#b42318",
         linewidth=1, label="Stochastic SIS SDE")
plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title("One stochastic SIS path and deterministic reference")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

state_grid = np.linspace(0, N, 501)
drift_values = np.array([
    sis_drift(i, N, beta, gamma)
    for i in state_grid
])
diffusion_values = np.array([
    sis_diffusion(i, N, beta, gamma)
    for i in state_grid
])

fig, left_axis = plt.subplots(figsize=(9.5, 5.2))
right_axis = left_axis.twinx()

left_axis.plot(
    state_grid, drift_values,
    color="#174c34", label="Drift"
)
right_axis.plot(
    state_grid, diffusion_values,
    color="#b42318", label="Diffusion"
)
left_axis.axhline(0, color="black", linewidth=0.8)
left_axis.set_xlabel("Current infectious population, I")
left_axis.set_ylabel("Drift, people per day", color="#174c34")
right_axis.set_ylabel(
    "Diffusion, people per √day", color="#b42318"
)
left_axis.set_title("SIS coefficients depend on the current state")
left_axis.grid(alpha=0.25)
fig.tight_layout()

Output

Run the code to see the result.

Why `max(total_event_rate, 0)` appears

Mathematically, the total rate is non-negative on \(0\le I\le N\). Tiny floating-point errors or an invalid numerical state could otherwise make a square root fail. This guard does not make an invalid biological model valid; the state and parameter checks remain necessary.

The boundary projection is a declared numerical choice

np.clip(proposed_I, 0, N) projects an Euler–Maruyama proposal back into the biological interval. It is simple, but it changes the raw numerical scheme. The program reports how often it was used. Later lessons compare boundary treatments and time-step effects rather than hiding them.

SDE versus CTMC SIS

CTMC	SDE approximation
Integer state and individual jumps	Continuous state and continuous path
Exact extinction event at \(I=0\)	Boundary behaviour depends on model and numerical treatment
Best for small counts and event detail	Often efficient for larger populations and approximate fluctuations

What this lesson adds

You can now derive a one-dimensional stochastic SIS model from its infection and recovery events, explain why rates subtract in drift but add in diffusion variance, program the state-dependent coefficients, and state the limitations of boundary projection.