Modeling Epidemics on Multiplex Networks: Epidemic Threshold and Basic Reproduction Number

TL;DR

This paper extends classical epidemic modeling to multiplex networks by formulating a degree-based mean-field SIR model across two layers and deriving a multiplex basic reproduction number $\mathcal{R}_0=\rho(J)$ via the Next Generation Matrix. It introduces a practical approximation $\tau$ that sums layer-wise reproduction contributions and demonstrates, through numerical and agent-based simulations, that an epidemic threshold occurs near $\mathcal{R}_0=1$ and that outbreak severity correlates with $\tau$. The results show that multiplex topology can alter transmission potential and that $\tau$ provides an interpretable, computable proxy for planning interventions. The work offers a more realistic foundation for forecasting and policy design in populations with layered social interactions.

Abstract

Accurate epidemic forecasting requires models that account for the layered and heterogeneous nature of real social interactions. The basic reproduction number $\mathcal R_0$ calculated from models that assume homogeneous mixing or single-layer contact structures have limited applicability to complex social systems. Here, we propose an expression of $\mathcal R_0$ in the context of multiplex networks, enabling the analysis of disease transmission across multiple social layers. We adapt the Degree-Based Mean-Field (DBMF) SIR model for single-layered complex networks to the multiplex setting, where each layer has its own degree distribution and infection rate. Using the Next Generation Matrix method, we derive an analytical expression for the basic reproduction number $\mathcal R_0$. Numerical integration of the multiplex DBMF equations shows that $\mathcal R_0 = 1$ marks the epidemic threshold and governs the functional dependence of key outbreak indicators. In addition to the exact result for the $\mathcal R_0$, we provide an approximation denoted as $τ$, which is easier to compute and more straightforward to interpret in terms of the parameters of the system, and shares most of the expected properties of the basic reproduction number. Stochastic agent-based simulations confirm these results, demonstrating a direct correspondence between $τ$ and the average number of secondary infections in the early epidemic phase, in line with the interpretation of $\mathcal R_0$. This research provides a robust generalization of $\mathcal R_0$ for layered contact structures, offering a more realistic basis for epidemic forecasting and the design of intervention strategies.

Figures (7)

• Figure 1: Maximum infection peak $I_\mathrm{max}= \max_t \;I_\mathrm{tot}(t)$ and final epidemic size $R_\infty= R_\mathrm{tot}(t\to\infty)$ attained in the standard SIR model. Their analytical expressions are $I_\mathrm{max} = 1-\frac{1+\ln \left(S(0) \mathcal{R}_0\right)}{\mathcal{R}_0}$ and ${R_\infty = 1 + \frac{W\left( - S(0) \mathcal{R}_0 e^{-\mathcal{R}_0}\right)}{\mathcal{R}_0}}$ respectively murray, where $W$ stands for the Lambert $W$ function.
• Figure 2: $I_\mathrm{max}$ as a function of the exact $\mathcal{R}_0$, denoted as $\rho$, given by Eq. (\ref{['R0_multiplex']}). Both layers of the network have a Poisson degree distribution. In each considered scenario, all the parameters of the system are fixed except for one, which varies in a certain range (see Table \ref{['table_poisson']} for details).
• Figure 3: $I_\mathrm{max}$ and $R_\infty$ as a function of $\tau = \frac{\beta_1}{\gamma}\frac{\langle k_1^2\rangle-\langle k_1\rangle}{\langle k_1\rangle} + \frac{\beta_2}{\gamma}\frac{\langle k_2^2\rangle-\langle k_2\rangle}{\langle k_2\rangle}$. Both layers of the network have a Poisson degree distribution, $P_i(k_i)=\frac{\lambda_i^{k_i}e^{-\lambda_i}}{k_i!}$. In each considered scenario, all the parameters of the system are fixed except for one, which varies in a certain range (see Table \ref{['table_poisson']} for details).
• Figure 4: $I_\mathrm{max}$ and $R_\infty$ as a function of $\tau$ when both layers of the network have a geometric degree distribution, $P_i(k_i)=\mu_i^{-1}(1-\mu^{-1}_i)^{k_i-1}$. See the details of each explored scenario in Table \ref{['table_geom']}.
• Figure 5: $I_\mathrm{max}$ and $R_\infty$ as a function of $\tau$. In all eight scenarios shown, the first layer of the network has a Poisson degree distribution with mean $\lambda$ while the second layer has a geometric one with mean $\mu$. See the details of each explored scenario in Table \ref{['table_mix']}.