Behavior–Disease Coupled Models

Classic compartmental models fix the transmission rate $\beta$ as a constant, but people are not passive substrate for a pathogen: they wash hands, mask, distance, and vaccinate when they perceive risk, and behavior–disease coupled models close the loop by letting the epidemic shape behavior and behavior reshape the epidemic. The result is a feedback system whose dynamics — flattened peaks, delays, and recurrent waves — differ qualitatively from any model with a frozen $\beta$ .

Epidemic prevalence over time with and without prevalence-dependent protective behavior, showing behavior flattens and delays the peak.

Why a fixed transmission rate fails

Transmission is a product of contact and per-contact risk, and both fall when a population reacts to an outbreak. As prevalence climbs, awareness rises, protective behavior spreads, and the effective transmission rate drops well below the value estimated early in an epidemic. A model with constant $\beta$ therefore overshoots and mistimes the true peak, treating a behaving population as if it never learned the disease was there; the systematic review of Verelst et al., 2016, J. R. Soc. Interface catalogues how widely this assumption has been relaxed.

The feedback loop

The coupling is a closed causal cycle. More infection raises perceived risk and awareness, which increases adoption of protective behavior, which lowers transmission, which reduces new infection — after which risk falls, behavior relaxes, and transmission climbs again.

$\text{infection} \;\to\; \text{awareness/behavior} \;\to\; \downarrow \text{transmission} \;\to\; \downarrow \text{infection} \;\to\; \text{relaxation} \;\to\; \cdots$

This is negative feedback with a delay — the classic recipe for oscillation — and the review by Weston et al., 2018, BMC Public Health surveys how infection-prevention behaviors have been written into this loop.

Ways behavior is encoded

Modelers close the loop in several distinct ways, ordered roughly from phenomenological to mechanistic.

Prevalence-dependent transmission. Make $\beta$ a decreasing function of current prevalence $I$ , so risk directly throttles contact without tracking individual decisions.

Information or awareness compartments. Split the population into unaware and aware classes, with awareness spreading like a second contagion and aware individuals transmitting less.

Imitation and game-theoretic decisions. Individuals weigh the cost of protecting against the perceived cost of infection and copy successful neighbors, as in the economic and evolutionary-game framing shared with evolutionary game theory and the evolution of cooperation.

Behavior-change theories. Embed a psychological model such as the Health Belief Model, so adoption depends on perceived susceptibility, severity, benefits, and barriers; Ryan et al., 2024, J. R. Soc. Interface fit exactly this into an SIRS transmission model.

The broad landscape of these couplings is mapped by Reitenbach et al., 2024, Rep. Prog. Phys..

A prevalence-dependent SIRS model

We take the simplest coupling: an SIRS system in which a fraction $P$ of the population adopts a protective behavior that scales transmission by $(1 - cP)$ . Here $c \in [0,1]$ is the maximum protection the behavior can confer, and adoption $P$ responds to current prevalence $I$ through a saturating response

$P(I) = \frac{I}{I + k},$

so that adoption is negligible when $I \ll k$ and approaches full uptake when $I \gg k$ ; the half-saturation constant $k$ sets how much disease it takes to move people. Writing the force of infection as $\lambda = \beta\,(1 - cP(I))\,I$ , the model is

$\begin{aligned} \frac{dS}{dt} &= -\beta\,(1 - cP(I))\,S I + \omega R, \\ \frac{dI}{dt} &= \beta\,(1 - cP(I))\,S I - \gamma I, \\ \frac{dR}{dt} &= \gamma I - \omega R, \end{aligned}$

with $S + I + R = 1$ . Setting $c = 0$ recovers the ordinary SIRS model; setting $c > 0$ engages the feedback loop.

A worked example

Take $\beta = 0.6$ , $\gamma = 0.2$ (so the unmitigated $R_0 = \beta/\gamma = 3$ ), waning immunity $\omega = 0.02$ , half-saturation $k = 0.02$ , and initial state $S_0 = 0.999,\; I_0 = 0.001,\; R_0 = 0$ . Integrating with behavior off ( $c = 0$ ) and then on ( $c = 0.7$ ) isolates the coupling, since every other parameter is held fixed: as prevalence rises past $k$ , the term $(1 - cP)$ falls toward $1 - c = 0.3$ , cutting effective transmission and pulling the peak down and later.

In code

R

library(deSolve)
rhs <- function(t, y, p) {
  P <- y["I"] / (y["I"] + p$k)          # saturating adoption
$  lam <- p $beta * (1 - p$ c * P) * y["I"]
  list(c(S = -lam*y["S"] + p$omega*y["R"],
$         I =  lam*y["S"] - p$gamma*y["I"],
$         R =  p $gamma*y["I"] - p$ omega*y["R"]))
}
peak <- function(c) {
  p <- list(beta=0.6, gamma=0.2, omega=0.02, k=0.02, c=c)
  out <- ode(c(S=0.999, I=0.001, R=0), seq(0,200,0.05), rhs, p)
  j <- which.max(out[,"I"]); c(out[j,"I"], out[j,"time"])
}
peak(0.0)   # peak prevalence and time, behavior off
peak(0.7)   # behavior on: lower and later

Python

import numpy as np
from scipy.integrate import solve_ivp

beta, gamma, omega, k = 0.6, 0.2, 0.02, 0.02
S0, I0, R0 = 0.999, 0.001, 0.0
t_end, t_eval = 200.0, np.linspace(0, 200, 4001)

def rhs(t, y, c):
    S, I, R = y
    P = I / (I + k)                  # prevalence-dependent adoption
    lam = beta * (1 - c * P) * I     # effective force of infection
    return [-lam * S + omega * R,
            lam * S - gamma * I,
            gamma * I - omega * R]

def peak(c):
    sol = solve_ivp(rhs, (0, t_end), [S0, I0, R0], args=(c,),
                    t_eval=t_eval, rtol=1e-8, atol=1e-10)
    I = sol.y[1]
    j = int(np.argmax(I))
    return I[j], sol.t[j]

for label, c in [("behavior off (c=0.0)", 0.0), ("behavior on  (c=0.7)", 0.7)]:
    pk, tp = peak(c)
    print(f"{label}: peak I = {pk:.4f} at t = {tp:.1f}")

behavior off (c=0.0): peak I = 0.3106 at t = 19.4
behavior on  (c=0.7): peak I = 0.0705 at t = 31.2

Julia

using DifferentialEquations
beta, gamma, omega, k = 0.6, 0.2, 0.02, 0.02
function rhs!(du, u, c, t)
    S, I, R = u
    P = I / (I + k)                  # saturating adoption
    lam = beta * (1 - c * P) * I
    du[1] = -lam*S + omega*R
    du[2] =  lam*S - gamma*I
    du[3] =  gamma*I - omega*R
end
function peak(c)
    sol = solve(ODEProblem(rhs!, [0.999, 0.001, 0.0], (0.0, 200.0), c), saveat=0.05)
    I = getindex.(sol.u, 2); j = argmax(I)
    (I[j], sol.t[j])
end
peak(0.0); peak(0.7)   # behavior on lowers and delays the peak

The printed peaks show the behavior off run cresting near $I \approx 0.31$ around $t \approx 19$ , while the behavior on run flattens to $I \approx 0.07$ and delays to $t \approx 31$ — from the coupling alone, with all epidemiological parameters unchanged.

Why it matters

Because the feedback is negative and delayed, coupled models naturally generate the recurrent waves and damped oscillations that constant- $\beta$ models can only produce by hand: behavior that relaxes as prevalence falls sets up the next wave, so multiple peaks emerge endogenously rather than being assumed. This matters for forecasting and for policy design, since an intervention’s effect depends on how people will respond to it, and integrating social and behavioral factors is a recognized priority for outbreak modeling (Bedson et al., 2021, Nat. Hum. Behav.). The central open problem is empirical: adoption functions like $P(I)$ are rarely fit to behavioral data, and this validation gap — models rich in mechanism but thin on measurement — is the recurring caution across the reviews cited above.