Epidemic mean-field thresholds of SIS dynamics on temporal multiplex networks with activity

Abstract We investigate the susceptible–infected–susceptible (SIS) process on a temporal multiplex network where a persistent static layer is overlaid with an activity-driven layer of one-step contacts. Simulations on three random network families (Erdős–Rényi, Watts–Strogatz, Barabási–Albert; \(N=10{,}000\) nodes each) show a clear extinction–persistence transition that shifts systematically as the number of temporal contacts per activation increases. To explain this behavior, we linearize the dynamics near the disease-free equilibrium and replace the stochastic temporal contacts by their mean effect, yielding a compact mean-field threshold \(\tau _c = 1/(\lambda _1 + 2am)\), where \(\lambda _1\) is the spectral radius of the static adjacency matrix, a is the activation probability, and m is the number of contacts per activation. This formula accurately tracks the empirical phase boundary for Erdős–Rényi and Watts–Strogatz backbones. For the Barabási–Albert (scale-free) network, however, a finite-N correction is required: the spectral formula underestimates the empirical threshold by approximately 42% at \(N=10{,}000\), a discrepancy attributable to degree heterogeneity and finite-size effects on scale-free networks. We introduce an empirically calibrated correction factor and provide a corrected threshold formula that substantially improves agreement with simulations. We also report surviving-runs prevalence estimates, extinction probabilities, and early-time growth rates as complementary diagnostics for locating the phase boundary. Acknowledgements This work has been supported by the USDA National Institute of Food and Agriculture, grant number 2022-67015-38059 through the NSF/NIH/USDA/BBSRC/BSF/NSFC Ecology and Evolution of Infectious Diseases Program. Author information Authors and Affiliations Corresponding author Ethics declarations Competing interest The authors declare that they have no competing interest. Use of AI and AI-assisted technologies During preparation of this work, the authors used AI-assisted tools to improve sentence clarity in selected passages and to assist with code documentation, cleanup, and debugging in the simulation software. The authors reviewed and edited all such content and take full responsibility for the content of this publication. Additional information Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Appendix A Derivation of the threshold on a static+activity-driven multiplex Appendix A Derivation of the threshold on a static+activity-driven multiplex A.1 Model recap and assumptions We consider a multiplex with a fixed (static) layer encoded by the symmetric adjacency matrix \(A\in \mathbb {R}^{N\times N}\) (zero diagonal), and a temporal activity-driven layer: at each discrete step \(t\rightarrow t+1\), every node activates independently with probability \(a\in (0,1)\) and, if active, creates \(m\in \mathbb {N}\) undirected links to uniformly chosen distinct targets; these links persist for one step and then disappear. SIS dynamics (per step): an infected node recovers with probability \(\mu \in (0,1)\); a susceptible node meeting an infected neighbor becomes infected with probability \(\beta \in (0,1)\) (independently per infectious contact within the step). We define \(\tau =\beta /\mu \). Throughout the derivation we use standard approximations for threshold analysis: (i) Pair-independence for the static layer: we write the probability that node i is susceptible while static neighbor j is infected as \((1-p_i)\,p_j\), i.e. we treat the states of i and j as statistically independent. On a static network the true joint probability \(P(S_i,I_j)\) deviates from this product because neighbor states are dynamically correlated; as a consequence \(\tau _c=1/\lambda _1\) is an approximation rather than an exact identity (Pastor-Satorras et al. 2015). (ii) Linearization around the DFE: we drop \(\mathcal {O}(p_i p_j)\) terms, valid when \(p_i\ll 1\) near the onset. (iii) Mean-field averaging of the temporal contacts: we replace the stochastic one-step temporal adjacency by its expectation. (iv) Homogeneous temporal mixing: the temporal exposure in the activity-driven layer is uniform across nodes. Assumptions (ii)–(iv) are classical in spectral/mean-field threshold analyses and are appropriate near the onset. A.2 Node-level update and linearization Let \(p_i(t)\in [0,1]\) denote the infection probability of node i at step t and let \(\bar{p}(t)=\tfrac{1}{N} \sum _{j=1}^N p_j(t)\) denote the population average. The discrete-time update can be written as where \(\Pi _i(t)\) is the probability that at least one infectious contact reaches node i during step t (across both layers). Close to the DFE, \(p_i(t)\ll 1\), so the product \((1-p_i(t))\,\Pi _i(t)\) can be linearized as \(\Pi _i(t)\) up to \(\mathcal {O}(p_i\Pi _i)=\mathcal {O}(p^2)\) terms. Subtracting \(p_i(t)\) from Eq. (11) then yields Interpreting one step as unit time (or taking the small-step limit) gives the continuous-time linearized form A.3 Expected infectious pressure \(\Pi _i(t)\) The probability \(\Pi _i(t)\) can be decomposed into contributions from the static and temporal layers. Under the linearization, where \(x_i^{\text {temp}}(t)\) is the expected number of infectious temporal contacts that reach node i during step t. Lemma A.1 (Expected temporal exposure) Under homogeneous activity a and m one-step contacts per activation, the mean temporal infectious exposure to node i during step t is where the approximation holds for large N or near the DFE (\(p_i(t)\ll 1\)). Proof Decompose \(X_i^{\textrm{temp}}(t)\) into outbound and inbound infected temporal contacts: \(\square \) Outbound (when i activates): with probability a, node i selects m distinct targets uniformly from \(N-1\) candidates. By linearity of expectation, the expected number of infected targets among those m is Inbound (when an infected node chooses i): for each \(j\ne i\), with probability \(p_j(t)\) node j is infected; conditional on infection, j activates with probability a and includes i among its m choices with probability \(m/(N-1)\). Summing indicators over \(j\ne i\), Adding the two contributions gives Eq. (14). A.4 Vector form Let \(\textbf{p}(t)=[p_1(t),\dots ,p_N(t)]^\top \) and let \(J:=\tfrac{1}{N}\textbf{1}\textbf{1}^\top \) denote the rank-one averaging operator so that \(J\,\textbf{p}=\bar{p}\,\textbf{1}\). Stacking Eq. (13) for all i yields Thus, the temporal layer contributes a rank-one all-to-all term \(2\beta a m\,J\). A.5 Dominant-eigenvector reduction Because A is nonnegative and (assumed) irreducible, by Perron–Frobenius there exists a strictly positive unit eigenvector \(\textbf{v}_1\) with \(A\textbf{v}_1=\lambda _1\textbf{v}_1\), where \(\lambda _1\) is the spectral radius of A. Near the DFE, the growth is dominated by this mode, so we project the dynamics along \(\textbf{v}_1\) and by substituting Eq. (16) into Eq. (15) and left-multiplying by \(\textbf{v}_1^\top \) (with \(\textbf{v}_1^\top \textbf{v}_1=1\)) gives By Cauchy–Schwarz, \(0< s=\textbf{v}_1^\top \textbf{1}\le \Vert \textbf{v}_1\Vert \,\Vert \textbf{1}\Vert =\sqrt{N}\), hence \(0< \tfrac{s^2}{N}\le 1\), with equality when \(\textbf{v}_1=\tfrac{1}{\sqrt{N}}\textbf{1}\) (homogeneous degree). Under the homogeneous temporal mixing assumption (empirically accurate for ER/WS backbones), we take \(\tfrac{s^2}{N} \approx 1\), where Therefore the DFE is unstable (and the infection grows initially) iff \(\tau (\lambda _1+2am)>1\). A.6 Threshold and critical integer number of temporal contacts Equation (18) implies the epidemic threshold or, equivalently, \(\beta _c=\mu /(\lambda _1+2am)\). For fixed \((\beta ,\mu ,a)\) and varying m, the smallest integer number of temporal contacts per activation that crosses the threshold is A.7 Remarks on accuracy and bounds The rank-one form in Eq. (15) implies that the effective linear operator is \(\beta (A+2am\,J)-\mu I\). Since J has eigenvalues \(\{1,0,\dots ,0\}\) with eigenvector \(\textbf{1}\), standard eigenvalue perturbation (Weyl inequalities) gives with equality on the right when \(\textbf{v}_1\) is exactly collinear with \(\textbf{1}\). Thus Eq. (19) is exact for regular/homogeneous graphs and provides an accurate approximation for many random-like networks (e.g., ER/WS), as confirmed by the simulations; for strongly heterogeneous backbones (e.g., BA), small deviations can appear, but the additive law remains an excellent predictor of the onset. Rights and permissions Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. About this article Cite this article Marjani, M., Chatterjee, S., Bhattacharyya, S. et al. Epidemic mean-field thresholds of SIS dynamics on temporal multiplex networks with activity-driven layers. Appl Netw Sci (2026). https://doi.org/10.1007/s41109-026-00824-0 Received: Accepted: Published: DOI: https://doi.org/10.1007/s41109-026-00824-0