Variance of the $SIS$ Epidemic on Networks: A Diffusion Approximation

Abstract: Functional laws of large numbers (FLLNs) describe the mean-field trajectory of epidemics on networks, but say nothing about the fluctuations around it. These fluctuations are governed by moments of the degree distribution not relevant at the level of the mean. A rigorous functional central limit theorem (FCLT) exists for the susceptible--infected ($SI$) process on configuration-model graphs, but no analogue exists for $SIS$, where recovery reintroduces vertices into the susceptible pool with partially known neighborhoods, breaking the clean neighborhood distribution the $SI$ derivation relies on. We develop a tractable variance approximation for Markovian $SIS$ on configuration-model graphs, combining Gleeson's approximate master equation (AME) framework with a van Kampen system-size expansion in the spirit of the $SI$ FCLT. We derive a closed drift and diffusion matrix for a reduced susceptible/$SI$-edge/$SS$-edge count vector and obtain the time-dependent covariance via the associated Langevin/Lyapunov equation. Validation against Gillespie simulation across Poisson, regular, and power-law networks shows close agreement, with deviations near the epidemic threshold and in strongly heterogeneous networks.

• The paper introduces a diffusion approximation framework that extends fluctuation analysis to SIS epidemic dynamics on complex networks.
• It employs system-size expansion, AME reduction, and hypergeometric closure to derive closed-form drift and diffusion terms.
• Numerical validations via Gillespie simulations highlight the model's practical implications for risk quantification and intervention planning.

Diffusion Approximation for SIS Fluctuations on Configuration-Model Networks

Introduction and Motivation

The deterministic trajectory of an epidemic on a network, as described by functional laws of large numbers (FLLN), is insufficient for quantifying stochastic variations—essential in predicting extinction probabilities, peak sizes, and risk-based interventions. While central limit theorems (CLTs) for the SI process yield elegant expressions for epidemic variance governed by higher-order moments of the degree distribution, there has been no comparable analytic framework for SIS dynamics on networks, due largely to the reinfection and memory effects introduced by recovery. This paper presents a tractable, time-dependent variance approximation for Markovian SIS dynamics on configuration-model graphs, extending the van Kampen system-size expansion combined with Gleeson’s approximate master equation (AME) formalism to yield a closed drift and diffusion description.

Model and Diffusion Approximation Framework

The SIS process is constructed as a continuous-time Markov chain on a network $G = (V, E)$, where each node toggles between susceptible ($S$) and infected ($I$), with infection and recovery rates $\beta$ and $\gamma$, respectively. The mean-field behavior is described at the degree-resolved level using AME [Gleeson 2013], tracking the fractions of susceptible/infected vertices of degree $k$ with $m$ infected neighbors. This extensive set of ODEs captures the evolution of compartmental statistics, e.g., the count vector $(X_{S}, X_{SI}, X_{SS})$ (susceptible, SI edges, SS edges).

The system-size expansion decomposes stochastic fluctuations $$\boldsymbol{X}(t) = N \boldsymbol{\phi}(t) + \sqrt{N}\,\boldsymbol{\zeta}(t)$$, linearizing around the AME trajectory. The covariance matrix $$C(t) = \text{Cov}(\boldsymbol{\zeta}, \boldsymbol{\zeta}^\intercal)$$ evolves according to a Lyapunov differential equation,

$S$0

where $S$1 is the Jacobian of the AME drift evaluated along the deterministic solution, and $S$2 is the event-wise noise matrix, encoding the effect of both infection and recovery.

Approximate Closure for High-Order Motifs

Exact closure requires explicit knowledge of the degree and neighbor-state distribution around each susceptible, which is disrupted in SIS by recovery events. To address this, the neighborhood of a degree-$S$3 susceptible is approximated as a hypergeometric draw from the current pool of SI and SS half-edges, leading to

$S$4

aligning with the SI configuration model limits for rapid enough renewal ($S$5) or high average degree. This facilitates the closure of wedge motifs such as $S$6, $S$7—crucial for expressing the drift, Jacobian, and noise in terms of only $S$8.

The reduction introduces structural factors

$S$9

encoding normalized $I$0-star motif counts (conditional on susceptibility), making explicit the dependence of variance on higher moments, as predicted by the SI FCLT [KhudaBukhsh et al. 2022]. For regular graphs, $I$1 are constants; for heterogeneous graphs, they capture the shift due to degree skewness and variation.

Numerical Results and Regimes of Validity

The validity of the diffusion approximation is systematically tested against extensive Gillespie simulations, across varying degree distributions and epidemic parameters.

Figure 2: Comparison of the AME-Langevin predictions for Var(s) with the empirical variance measured from Gillespie simulations as a function of time $I$2 on Poisson configuration-model networks; agreement is strong except near threshold and in highly heterogeneous regimes.

Agreement is robust for Poisson random graphs away from threshold, accurately capturing both transient variance peaks and stationary values. The largest discrepancies occur near the epidemic threshold, where finite-size corrections and branching-time to equilibrium are not fully captured by the diffusion expansion.

Figure 1: Numerical validation of variance estimates across different network sizes demonstrates the $I$3 scaling of amplitude and increasing accuracy with system size.

The approximation is also validated for varying $I$4. The $I$5 scaling of variance is preserved, and the mean-field approximation improves monotonically with increasing network size, in line with central-limit scaling.

Figure 3: Validation across $I$6-regular networks with varying degree, showing best correspondence for higher $I$7 (i.e., more neighbors per vertex).

For $I$8-regular networks, the method is exact in the limit $I$9 and performs well for moderate $\beta$0. For low degree (sparse networks), discrepancies arise, especially at the variance peak. This is due to edge-wise correlations that cannot be averaged out in small neighborhoods, causing the hypergeometric closure to underestimate or overestimate higher moments.

Figure 4: Application to configuration-model networks with a truncated power-law degree distribution. The timing of the variance peak is captured, but the amplitude is systematically underestimated for strong heterogeneity ($\beta$1).

On scale-free (power-law) networks, agreement is restricted: while the AME-Langevin system localizes the peak in prevalence variance, the amplitude is misestimated due to the repeated reinfection into and out of influential hub nodes—correlations not well-approximated by moment closure. The analysis reveals a fundamental breakdown of the mean-field assumption in the presence of very large degree heterogeneity, tracing the error directly to inapplicability of the susceptible neighbor sampling closure.

Theoretical and Practical Implications

This approach extends the classic SI functional CLT reasoning to the richer, technically challenging SIS setting. The time-dependent variance is governed by degree distribution moments beyond the mean and variance, exposing the critical role of degree skewness and motif counts in transient and stationary epidemic fluctuations. The explicit encapsulation of network structure through factors like $\beta$2 and $\beta$3 enables rapid computational exploration of epidemic uncertainty, bypassing the need for large-scale stochastic simulation. For practical risk quantification, this framework gives immediate access to prevalence uncertainty curves—vital for scenario planning and control strategies sensitive to large deviations.

Computational tractability is a central advantage: the reduced $\beta$4 ODE system for the covariance matrix requires minimal resources compared to the thousands of stochastic realizations needed for empirical variances, making it suitable for parameter sweeps and sensitivity analyses.

Limitations and Future Directions

Leading error sources are present in the closure of motif statistics for strongly non-regular graphs and in the neglect of derivatives of the wedge structural factors with respect to the reduced variables in the Jacobian computation. The current method treats $\beta$5 as time-dependent functions input from the AME solution and omits their derivatives; exactness is restored for regular or weakly heterogeneous graphs. To address the strong correlations in hub-dominated networks, future work may introduce degree-resolved or neighborhood-resolved Langevin expansions, or dynamic corrections based on local motif accumulation.

A more systematic inclusion of the time-dependence and derivatives of the wedge factors in the Jacobian could refine predictions for heterogeneous graphs. The methodology is adaptable to other binary-state network processes (SIR, SIRS, contact processes, etc.), and may inform optimal intervention design, early warning, and forecast uncertainty quantification, especially under incomplete network information.

Conclusion

A diffusion approximation for the SIS epidemic on configuration-model networks, constructed via AME reduction and system-size expansion, provides an analytically tractable estimate for the time-dependent prevalence variance. This work demonstrates both the feasibility and the practical utility of extending fluctuation analyses from irreversible to recurrent dynamics, clarifying the network structural dependencies of macroscopic uncertainty in finite systems, and motivating refined, scalable approaches for complex networked dynamics.