Mesoscale community organization governs epidemic onset and spread in metapopulations

TL;DR

The paper develops a two-scale framework for epidemic spreading in hierarchically structured metapopulations by embedding local contact networks within a diffusion-driven mobility layer. Using IBMF and a degree-based DBMF reduction, it derives an effective transmission rate tilde{β}_{μ} = β ⟨k⟩_{μ} / N and shows that internal community connectivity drives onset and spread, with stability analyzed via spectral perturbation theory. A key finding is that above-average-density metacommunities amplify transmission and skew spatial infection patterns, a result made concrete through second-order perturbation corrections and localization of Laplacian eigenvectors. The authors further demonstrate predictive power by linking the leading eigenvector of the Jacobian to the steady-state infection distribution and validating the approach against nine empirical mobility networks. Overall, the work provides a scalable, theory-grounded toolkit for vulnerability assessment and targeted interventions in structured populations, with implications for metaplex-inspired multiscale dynamics.

Abstract

Understanding how internal community structure shapes the course of epidemics remains a fundamental challenge in modeling real-world populations. Standard metapopulation models often assume uniform mixing within communities, overlooking how internal heterogeneity affects global outcomes. Here, we develop a general framework for epidemic spreading in hierarchically structured metapopulations, where individuals interact locally within dense communities and move across a broader network. We show that transmission dynamics are governed by the mesoscale organization of these communities: highly connected groups accelerate and amplify outbreaks, while less connected ones dampen spread. Through a combination of mean-field theory, spectral analysis, and stability methods, we reveal a direct link between internal connectivity and the emergence of uneven, spatially structured epidemic patterns. We further validate these predictions using real-world data, where social contact networks capture the local scale of transmission while spatial transport networks govern global connectivity, confirming the robustness of our framework across scales. These results demonstrate how community structure fundamentally governs the shape of epidemics in complex, networked populations, offering new insights into vulnerability, containment, and epidemic control.

Figures (5)

• Figure 1: Schematic of the spreading process on a hierarchically structured network. Individuals—either susceptible $S$ (blue) or infected $I$ (red)—can interact (dashed links) when co-located within the same spatial patch, represented by a metanode (light blue). Migration between metanodes (solid links) occurs via assigned positions within the destination metanode, ensuring individuals are associated with a single metanode at any given time.
• Figure 2: The mean equilibrium state, $\langle x \rangle^*$, is plotted against the recovery rate, $\gamma$ for different networks parameters realizations. In the metapopulation network, both the transport and contact networks are Erdős–Rényi graphs with wiring probabilities $q$ (between metanodes) and $p$ (within each contact network), respectively. For both panels, each scenario corresponding to a different $p$ value list is simulated 10 times using both the IBMF and DBMF$_{\mu}$ frameworks. In each simulation, a distinct distribution of $\langle k \rangle_{\mu}$ is generated. The shaded regions in the IBMF results represent the range of observed values across the 10 simulations, with the solid line within the shaded region indicating the mean outcome. The dashed line corresponds to the average result obtained from DBMF$_{\mu}$ simulations. The shaded region is constructed by setting its boundaries using the minimum and maximum values of $\langle k \rangle_{\mu}$ obtained from the simulations. The solid line within this region represents the mean value obtained from the IBMF simulations. (a) The distribution of $\langle k \rangle_{\mu}$ follows a normal distribution with a mean of 0.5 and varying variance $\sigma$. (b) The distribution of $p$ is systematically varied to include both left-skewed and right-skewed forms while maintaining a similar mean. Specifically, a gamma distribution ($\alpha = 1, \theta=2$) is utilized to generate left-skewed $p$, and the right-skewed $p$ is obtained by applying a reflection transformation to the left-skewed distribution. (c) Empirical social networks: we apply the same protocol to fixed real social graphs, plotting $\langle x \rangle^{*}$ versus $\gamma$ with IBMF (solid) and DBMF$_{\mu}$ (dashed); dataset definitions and analysis are detailed in the Supplemental Information. The set of parameters for all panels are $\Omega=10, N=50, q=0.5, \beta=0.5, D_S=0.2, D_I=0.1$.
• Figure 3: The figure illustrates the distribution of metanodes along the horizontal axis, with the vertical axis divided into two sections. Panels (a)--(c) show synthetic benchmarks where the upper plots compare the asymptotic infection state $\mathbf{I}^{*}$ obtained from IBMF simulations with the critical eigenvector $\boldsymbol{\phi}^{(\Omega+1)}_{I}$ and its successive approximations, while the lower plots compare the deviation vector $\boldsymbol{\kappa}$ with the first-order correction $\boldsymbol{\phi}^{(\Omega+1)}_{I,1}$; the imposed condition $\tilde{\beta} S^{*} - \gamma = -0.01$ ensures distinct eigenvalues, and only $\lambda^{(\Omega+1)}$ is controlled to be positive. Panels (d)--(f) display eigenvector localization patterns for the Winnipeg traffic system, the London bike-sharing network, and the Eastern Massachusetts transportation network, highlighting the emergence of strongly localized modes in empirical structures. Panel (g) reports the direct comparison between the IBMF asymptotic infection state, the critical eigenvector, and its first-order correction for the Sioux Falls traffic network, using the same parameter setting as in panels (a)--(c) to establish consistency between synthetic and empirical results. Panel (h) summarizes localization properties across the empirical datasets through a radar plot that jointly displays the Pearson correlation and normalized IPR, with methodological details provided in Appendix \ref{['sec:local_valid']}. For the synthetic networks, we use the parameters $\Omega = 20$, $N = 50$, and $q = 0.5$, with probabilities drawn from a normal distribution of mean $0.5$ and standard deviations of $0.01$, $0.05$, and $0.08$. In both the synthetic and Sioux Falls cases, the dynamical parameters are $\beta = 0.5$, $D_S = 0.2$, and $D_I = 0.1$.
• Figure 4: The graph illustrates the two types of eigenvalues, $\lambda^{(\alpha)}_0$ and $\lambda^{(\bar{\alpha})}_0$, both corresponding to $\Lambda^{(\alpha)}$ but differing in slope, associated with $D_S$ and $D_I$, respectively. Notably, only $\lambda^{(\bar{\alpha})}_0$ can become unstable or strictly stable, as illustrated by the dashed line. The inset provides a zoomed-in view, confirming that in our setting, the largest eigenvalue for $\lambda^{(\bar{\alpha})}_0$ is chosen as $\tilde{\beta} S^{*} - \gamma = -0.01$ to prevent spectral degeneracy.
• Figure 5: Localization and structural matching of Laplacian eigenvectors. Panels (a)--(d) show heatmaps of the absolute values of the Laplacian eigenvector components for networks with 10, 20, 50, and 100 nodes, respectively. Each column represents an eigenvector, ordered by increasing eigenvalue. Panels (e)--(h) display the entry-wise comparison between the eigenvector $\boldsymbol{\phi}^{(\Omega+1)}_{I,1}$ and the structural vector $\boldsymbol{\mathcal{K}}$, generated as a sawtooth function, for each corresponding network size. Panels (i) and (j) report the root-mean-square error (rMSE) and the Pearson correlation, respectively, between $\boldsymbol{\phi}^{(\Omega+1)}_{I,1}$ and $\boldsymbol{\mathcal{K}}$, averaged over 50 independent network realizations with fixed parameters. Error bars indicate standard deviations across realizations.