# Analytical Computation of the Epidemic Threshold on Temporal Networks

Eugenio Valdano, Luca Ferreri, Chiara Poletto, Vittoria Colizza

► To cite this version:

Eugenio Valdano, Luca Ferreri, Chiara Poletto, Vittoria Colizza. Analytical Computation of the Epidemic Threshold on Temporal Networks. Physical Review X, 2015, 5 (2), pp.021005. 10.1103/Phys-RevX.5.021005. hal-01275250

HAL Id: hal-01275250

https://hal.sorbonne-universite.fr/hal-01275250

Submitted on 17 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analytical Computation of the Epidemic Threshold on Temporal Networks

Eugenio Valdano,¹,² Luca Ferreri,³ Chiara Poletto,¹,² and Vittoria Colizza¹,²,⁴,*

¹INSERM, UMR-S 1136, Institut Pierre Louis d'Epidémiologie et de Santé Publique, F-75013 Paris, France

²Sorbonne Universités, UPMC Univ Paris 06, UMR-S 1136, Institut Pierre Louis d'Epidémiologie et de Santé Publique, F-75013 Paris, France

³Dipartimento di Scienze Veterinarie, Università degli Studi di Torino, Grugliasco (TO) 10095, Italy

⁴ISI Foundation, Torino 10126, Italy

(Received 18 August 2014; published 8 April 2015)

The time variation of contacts in a networked system may fundamentally alter the properties of spreading processes and affect the condition for large-scale propagation, as encoded in the epidemic threshold. Despite the great interest in the problem for the physics, applied mathematics, computer science, and epidemiology communities, a full theoretical understanding is still missing and currently limited to the cases where the time-scale separation holds between spreading and network dynamics or to specific temporal network models. We consider a Markov chain description of the susceptible-infectious-susceptible process on an arbitrary temporal network. By adopting a multilayer perspective, we develop a general analytical derivation of the epidemic threshold in terms of the spectral radius of a matrix that encodes both network structure and disease dynamics. The accuracy of the approach is confirmed on a set of temporal models and empirical networks and against numerical results. In addition, we explore how the threshold changes when varying the overall time of observation of the temporal network, so as to provide insights on the optimal time window for data collection of empirical temporal networked systems. Our framework is of both fundamental and practical interest, as it offers novel understanding of the interplay between temporal networks and spreading dynamics.

DOI: 10.1103/PhysRevX.5.021005

Subject Areas: Complex Systems, Interdisciplinary Physics, Statistical Physics

I. INTRODUCTION

A wide range of physical, social, and biological phenomena can be expressed in terms of spreading processes on interconnected substrates. Notable examples include the spread of directly transmitted infectious diseases through host-to-host contacts [1], the spatial propagation of epidemics driven by the hosts' mobility network [1-3], the spread of cyber viruses along computer connections [4], or the diffusion of ideas mediated by social interactions [5,6]. These phenomena are the result of a complex interplay between the properties of the spreading dynamics and the network's structural and temporal features, hindering their full understanding.

A fundamental issue characterizing spreading processes is the identification of the critical condition for the wide-spreading regime, encoded in the epidemic threshold parameter. This issue is of critical importance for epidemic containment [1], as well as for control of the diffusion of information [7] and cyber viruses [4]. Extensive studies have characterized this parameter in the time-scale separation approximation, i.e., when the time scales of the spreading process and network evolution strongly differ. This includes the two limiting regimes, quenched and annealed [4,8-15]. In the first case, the network is regarded as static, as it evolves on much slower time scales than the ones characterizing the spreading process. The epidemic threshold in this case is computed from the adjacency matrix describing the network connectivity pattern [8,9]. In the second case, the underlying network evolves so rapidly with respect to the dynamical process that only its time-averaged properties are relevant to the spreading dynamics. Approaches like the heterogeneous mean field [4], the generating function [11], and percolation theory [10] provide, in this regime, estimates of the threshold in the infinite size limit.

Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*Corresponding author.
vittoria.colizza@inserm.fr
http://www.epicx-lab.com</footnote>
Recently, the extensive empirical characterization of social interactions at different scales and settings [16–22] has shown that networks often display non-Poissonian and non-Markovian temporal evolutions unfolding at time scales similar to the ones of many spreading processes of interest, stressing the need for novel theoretical tools able to overcome current limitations. Much research has focused on spreading processes occurring on time-varying networks [16, 17, 21, 23–31], modeled either as a discrete-time sequence of networks [16, 28] or as continuous-time dynamics of links [17, 26]; however, so far, only a few studies have provided an analytical calculation of the epidemic threshold in specific cases [24, 27–32]. These are all based on models for time-varying networks integrating the microscopical laws governing the network evolution, under context-specific assumptions. An analytical framework for the computation of the epidemic threshold for an arbitrary time-varying network is still missing. To fill such a gap, we present here a novel approach that, by reinterpreting the tensor formalism of multilayer networks [33, 34], extends the Markov chain approach adopted for static networks [8, 9] to their temporal counterpart. The approach is applied to discrete time-varying network models and empirical networks to highlight the role of different dynamical features on the spreading potential. The role of the observation time window is then analyzed in depth in order to provide indications on how this factor alters the estimated epidemic threshold.

## II. DERIVATION OF THE EPIDEMIC THRESHOLD

We consider the susceptible-infected-susceptible (SIS) model [1] in discrete time, where individuals (i.e., the nodes of the network) can be in one of two mutually exclusive states—susceptible or infectious. At each time step, infectious individuals may transmit the infection to susceptible neighbors with probability $\lambda$ along each contact, and they recover spontaneously with probability $\mu$ becoming susceptible once again. We consider the temporal network forming the substrate of the spreading process to be a sequence of undirected and unweighted static networks. The generalization of the following treatment to the directed and weighted case is outlined in Ref. [35].

In order to describe the spreading dynamics on such a substrate, we extend the Markov chain approach for static networks [8, 9] to the case of temporal networks. The SIS propagation on a generic network with $N$ nodes and adjacency matrix $A$ is given by

$$p_i^{(t)} = 1 - \left[1 - (1 - \mu)p_i^{(t-1)}\right] \prod_j \left[1 - \lambda A_{ji} p_j^{(t-1)}\right], \quad (1)$$

where $p_i^{(t)}$ is the probability for the node $i$ to be in the infectious state at time $t$. The Markovian model of Eq. (1), widely adopted in different fields [12, 36], is based on the mean-field assumption of the absence of dynamical

correlations among the states of neighboring nodes [37]. For both directed and undirected networks [38, 39], the study of the asymptotic state yields the derivation of the epidemic threshold ($\lambda/\mu$) = $1/\rho(A^\dagger)$, where $\rho(A^\dagger)$ is the spectral radius of the transposed adjacency matrix $A^\dagger$ [8, 9]. This is known to be a lower bound estimate of the real epidemic threshold, approaching the real value with surprisingly high accuracy given the simplicity of the expression and its derivation [37, 40].

We extend this paradigm to a temporal network by letting the adjacency matrix in Eq. (1) depend on time:

$$p_i^{(t)} = 1 - \left[1 - (1 - \mu)p_i^{(t-1)}\right] \prod_j \left[1 - \lambda A_{ji} p_j^{(t-1)}\right]. \quad (2)$$

Here, $A^{(t)}$ is the adjacency matrix associated with the $t$th snapshot of the evolving network. In order to ensure the asymptotic solution of the SIS process in a generic temporal network, we assume periodic boundary conditions for the network dynamics. With $T$ being the total number of network time snapshots, we impose $A^{(T+1)} \equiv A^{(1)}$. This does not imply any loss of generality given that $T$ may be completely arbitrary. We notice that, as a consequence of the assumed periodic temporal dynamics of $A^{(t)}$, the asymptotic solution of Eq. (2) is, in principle, periodic, with period $T$.

We now define a more convenient representation of the coupled dynamics adopting the multilayer approach introduced in Ref. [33]. We map the temporal network to the tensor space $\mathbb{R}^N \otimes \mathbb{R}^T$, where each node is identified by the pair of indices $(i, t)$, corresponding to the node label $i$ and the time frame $t$, respectively. A multilayer representation of the temporal network can be introduced through the following rules:

(i) Each node, at time $t$, is connected to its future self at $t+1$.

(ii) If $i$ is connected to $j$ at time $t$, then we connect $i$ at time $t$ to $j$ at time $t+1$, and $j$ at time $t$ to $i$ at time $t+1$.

The second rule is termed “nondiagonal coupling” in the multilayer-network framework [34]. The first rule is consistent with the ordinal coupling in such a framework [34, 41], but unlike in that representation, no links are found connecting nodes on the same layer since layers cease to correspond to the adjacency matrices of the temporal snapshots. The so-defined network is therefore multipartite since only pairs of nodes belonging to different layers are linked together (see Fig. 1 for a schematic illustration of this transformation). While formally falling into a specific subcase of the classification introduced in Ref. [42], the proposed mapping from the network temporal sequence to a multilayer object provides a novel representation of the temporal network that preserves the information relevant for the spreading process. The tensor representation of the obtained multilayer network is the following:
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic example of the supra-adjacency matrix of the multilayer representation of the temporal network. For simplicity, we consider a network of two nodes i, j, and two time steps. The left panel represents the network as a sequence of static adjacency matrices. This is translated into a multilayer representation (right panel), where each node points to itself in the future and to the future image of its present neighbors.</figcaption></figure>

$$A_{ij}^{tt'} = \delta^{t,t'+1} [\delta_{ij} + A_{ij}^{(t)}]. \quad (3)$$

Analogously to the definition of A, we can also write, in this representation, the tensor associated with the SIS dynamics of Eq. (2), coupling together contagion and network dynamics:

$$\mathbf{M}_{ij}^{tt'} = \delta^{t,t'+1} [(1-\mu)\delta_{ij} + \lambda A_{ij}^{(t)}]. \quad (4)$$

The multilayer representation and the definition of the tensor M introduce a simplified expression for Eq. (2). The tensor space can be represented in single-index notation through the isomorphism $\mathbb{R}^N \otimes \mathbb{R}^T \cong \mathbb{R}^{NT}$. In other words, similarly to the definition of the supra-adjacency matrix in Refs. [33,43,44], we can mask the tensorial origin of the space through the map $(i, t) \to \alpha = Nt + i$, with $\alpha$ running in $\{1, \dots, NT\}$, allowing us to write the network tensor M in matrix form:

$$M = \begin{pmatrix} 0 & 1 - \mu + \lambda A^{(1)} & 0 & \dots & 0 \\ 0 & 0 & 1 - \mu + \lambda A^{(2)} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 - \mu + \lambda A^{(T-1)} \\ 1 - \mu + \lambda A^{(T)} & 0 & 0 & \dots & 0 \end{pmatrix}$$

M provides a network representation of the topological and temporal dimensions underlying the dynamics of Eq. (2), which are interrelated and flattened here. Its directed nature preserves the causality of the process, while its weights account for the SIS transition probabilities. The Markov process is now described by a trajectory in $\mathbb{R}^{NT}$ where the state vector $\hat{p}_\alpha(\tau)$ represents the probability of each node to be infected at each time step t included in the interval $[\tau T, (\tau + 1)T]$. Consistently, Eq. (2) becomes

$$\hat{p}_{\alpha}(\tau) = 1 - \prod_{\beta} [1 - M_{\beta\alpha} \hat{p}_{\beta}(\tau - 1)]. \quad (5)$$

Given that vector $\hat{p}$ encodes a one-period configuration, the T-periodic asymptotic state of the SIS process is now mapped into the steady state $\hat{p}_\alpha(\tau) = \hat{p}_\alpha(\tau - 1)$. The latter can be recovered as a solution of the equilibrium equation:

$$\hat{p}_{\alpha} = 1 - \prod_{\beta} (1 - M_{\beta\alpha} \hat{p}_{\beta}), \quad (6)$$

which is formally the same as the stationary condition imposed in Eq. (1) for the static network case, and it is similar to the Markov chain approaches used to solve

contagion processes in multiplex and interconnected networks [43–45]. We can then follow Refs. [8,9] and linearize Eq. (6), recovering the necessary and sufficient condition for the asymptotically stable zero solution, $\rho(M^\dagger) < 1$ [46]. Considering that the uniform zero solution in the $\mathbb{R}^{NT}$ representation is mapped to a uniform zero solution in the original $\mathbb{R}^N$ representation, this yields the threshold condition

$$\rho(M^\dagger) = 1 \quad (7)$$

for the critical values of $\lambda$ and $\mu$ above which the transmission becomes epidemic [8,9,43–45].

The spectral radius of M can be simplified with the following relation (see Appendix A):

$$\rho(M) = \rho(P)^{1/T}, \quad (8)$$

where $P = \prod_{t=1}^T (1 - \mu + \lambda A^{(T-t)})$ represents a weighted version of the accessibility matrix [47], having connectivity weighted by $\lambda$ and waiting time weighted by $1 - \mu$. This last passage ensures a simplification of the numerical computation of the epidemic threshold, allowing an execution time scaling as $\sim TN^{5/2}$ (see Ref. [48] for an analysis of the numerical performance of our approach).
The quenched and annealed regimes can be recovered within this general framework as particular limiting solutions. In the first case, it is to be noted that the sequence of temporal snapshots naturally defines the minimum time scale of the process. In order to consider contagion dynamics that are much faster than the time-varying process of the network, we thus rely on the commonly adopted assumption regarding the temporal network as static, so $A^{(t)} \equiv A$. In this particular case, $P = (1 - \mu + \lambda A)^T$. Therefore, $\rho(M) = \rho(1 - \mu + \lambda A) = 1 - \mu + \lambda\rho(A)$. The requirement $\rho(M^\dagger) = 1$ thus recovers the expression known for the quenched case.

The study of the annealed regime is less trivial. In the assumption that $\lambda$ and $\mu$ are very small, corresponding to a very slow disease dynamics with respect to the time scale of the network evolution, it is possible to replace $P$ with its linear expansion in $\lambda/(1 - \mu)$, yielding

$$P_{\text{slow}} = (1 - \mu)^T \left[ 1 + \frac{\lambda}{1 - \mu} A \right],$$

where $A = \sum_t A^{(t)}$ is a weighted static representation of the network, formed by the sum of all the snapshots. Temporal correlations are lost, and edges count for the number of times they are active during the whole period $T$. Equation (7) for the epidemic threshold thus simplifies to $(\lambda/\mu)_{c,\text{slow}} = T/\rho(A)$ [49], and the aggregated matrix contains all the relevant information for spreading dynamics.

### III. VALIDATION AND COMPARISON WITH STOCHASTIC SIMULATIONS

In the following, we validate the analytical method and compare its predictions with the behavior of a simulated SIS spreading process. For this purpose, we consider six networks, three of which are built from models for time-varying networks; the other three are obtained from empirical measures. The first network, ER, is formed by a sequence of random Erdős-Rényi graphs [50] with a given number of nodes and edges. It represents a simple and completely uncorrelated example of a temporal network. The second network, ACTIVITY, is a realization of the activity-driven model [28] where each node is assigned an activity potential, representing the probability of being active in a certain snapshot. Once activated, the node establishes a fixed number of connections that are renewed at each snapshot. We consider a heterogeneous activity distribution so that the obtained networks are characterized by a temporally uncorrelated sequence of snapshots with a heterogeneous aggregated degree distribution. BURSTY is built from the model introduced in Ref. [23] and accounts for a heterogeneous activation pattern describing a sequence of homogeneous networks where the intercontact time is power-law distributed. Size and period are chosen arbitrarily for all these networks since the choice of these parameters does not impact the method validation, as

discussed in more detail in the following section. As examples of real time-varying networks, we choose data sets describing different kinds of human contact: HT09 is the network of face-to-face proximity during a 2.5-day scientific conference [16]; SEX is a 1-year network of sexual contacts between prostitutes and their clients [18]; SCHOOL is a contact network describing one day in a high school [51]. Size, period, and topological properties are constrained by the measurements and are very diverse. Further information about the six networks can be found in Appendix B and in Table I.

To verify the validity of the proposed analytical expression, we numerically solve the Markov equation (2). For given $\lambda$ and $\mu$, we iterate the equation until the periodic state is reached and compute the average prevalence over a period $\langle i_{\text{MC}} \rangle = \sum_{i,t} p_i^{(t)}/(TN)$. Predictions are then also compared with the threshold behavior obtained from numerical simulations of the stochastic and microscopic SIS dynamics on the evolving networks. We use the quasistationary state method [52] (see Appendix C) to measure the average prevalence $\langle i_{\text{sim}} \rangle$ over the time series for different values of $\lambda$, after an initial transient time is discarded.

Figure 2 shows $\langle i_{\text{MC}} \rangle$ and $\langle i_{\text{sim}} \rangle$ as functions of $\lambda$ for two different values of $\mu$ ($\mu = 0.2$ and $\mu = 0.5$) for all networks under study. The average prevalence displays the expected transition behavior. The solution of the Markov chain equation $\langle i_{\text{MC}} \rangle$ is equal to zero for small values of $\lambda$ until the critical value of $\lambda$ is reached, after which a rapid growth is observed signaling an epidemic affecting a finite fraction of the network. The transition is well predicted by the analytical expression of Eq. (7). The threshold behavior obtained from numerical simulations is also very similar to the mean-field prediction. The two curves of $\langle i_{\text{MC}} \rangle$ and $\langle i_{\text{sim}} \rangle$ are nearly superimposed, showing that the mean-field approximation in Eq. (2) is valid in all conditions of

<table><caption>TABLE I. Temporal networks considered for the validation. The first three networks are single realizations obtained from synthetic models for time-varying networks; the other three are empirical networks. The ER model is a sequence of random graphs with 500 nodes and 750 edges, so $\langle k \rangle = 3$. The ACTIVITY model is a sequence of snapshots built with parameter values: $\Delta t = 1$, $m = 2$, $\eta = 10$, $\gamma = 2.8$, $\epsilon = 3 \times 10^{-2}$, in the notation of Ref. [28]. The BURSTY network is built with a power-law distributed interactivation time, with exponent $-2$, and cutoff equal to the period of the network. For the real networks, the collection time is the total time considered in the data set.</caption><thead><tr><th>Network</th><th>Number of nodes</th><th>Period T</th><th>Aggregating window</th></tr></thead><tbody><tr><td>ER [50]</td><td>500</td><td>13</td><td>...</td></tr><tr><td>ACTIVITY [28]</td><td>1000</td><td>20</td><td>...</td></tr><tr><td>BURSTY [23]</td><td>500</td><td>50</td><td>...</td></tr><tr><td>HT09 [16]</td><td>113</td><td>30</td><td>1 hour</td></tr><tr><td>SCHOOL [51]</td><td>787</td><td>42</td><td>10 mins</td></tr><tr><td>SEX [18]</td><td>6866</td><td>13</td><td>28 days</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Validation of the analytical method and comparison with microscopic numerical simulations. Top panel: Network models. Bottom panel: Empirical networks. Cross symbols represent $\langle i_{\text{MC}} \rangle$ as a function of $\lambda$ for two different values of $\mu$ ($\mu = 0.2$ in blue and $\mu = 0.5$ in red), i.e., the average prevalence obtained from the numerical solution of the Markov chain, Eq. (2). Circles represent $\langle i_{\text{sim}} \rangle$, i.e., the average prevalence obtained from stochastic microscopic numerical simulations of the SIS process, for the same values of $\mu$. The arrows indicate the analytical predictions of the threshold from Eq. (8) (one single realization of the network models is considered in each panel).</figcaption></figure>

network size and average connectivity here considered. The presence of correlations shows its effects in proximity to the transition, which is smoother for $\langle i_{\text{sim}} \rangle$ with respect to $\langle i_{\text{MC}} \rangle$. This is particularly evident for the network HT09 and is a consequence of its small size ($N = 113$). In Ref. [53], we report the analysis of the dynamical correlations.

The good agreement between the computed epidemic threshold, the solution of the Markov chain equation, and the numerical simulation results is thus maintained under a range of different temporal network properties (presence or absence of temporal correlations, heterogeneous vs homogeneous distributions characterizing temporal and structural observables, and the possible presence of community structure as in the case of school) and sizes (from approximately $10^2$ nodes to $10^4$). It is important to mention that periodic boundary conditions in the case of real networks may, in principle, induce nonexisting phenomena (such as, for example, temporal paths [54]) that could alter the threshold estimation by influencing the spreading process. We analyze the effect of our technical assumption of adopting periodic constraints in the following section, also in relation to data availability and collection.

## IV. OPTIMAL DATA COLLECTION TIME

Available data sets characterizing empirical networks only account for a portion of the real contact process, and

the extent of the recording time window may affect the prediction of the epidemic threshold. One may expect that, when the data-collection period is long enough, the data would represent an approximately complete reconstruction of the temporal network properties, thus leading to an accurate estimate of the epidemic threshold. Given the resources needed for the setup of data-collection deployments, here we explore the role of the period $T$ aimed at identifying a minimum length of observation of the contact process that is optimal in providing a reliable characterization of the spreading potential.

We thus compute the epidemic threshold from Eq. (8) for increasingly larger values of the period $T$ up to the entire data-collection time window, for the three empirical networks under study. Figure 3 shows a saturation behavior for $\lambda_c$, indicating that the data-collection period is long enough to characterize the epidemic dynamics. Such behavior and its associated relaxation time strongly depend on the network's typical time scale and on the temporal variability of its structure. In more detail, a simple structural measure—the variation of the average degree along the period—is shown to strongly impact the predicted $\lambda_c$ (Fig. 3 and Ref. [55]). This is particularly evident in the SCHOOL network, where the daily activity of students determines considerable variations in the average degree and induces marked oscillations in the resulting threshold. In addition to the empirical networks, we also consider the BURSTY network model as it includes nontrivial temporal correlations.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Epidemic threshold estimated from different period lengths. $\lambda_c(T)$ is the epidemic threshold computed by considering only the first $T$ snapshots of the network. For each panel, the blue (red) curve corresponds to $\mu = 0.2$ ($\mu = 0.5$), and its scale is reported on the left (right) side of the plot. The gray bar chart shows the mean network degree associated with the snapshot at time $T$. Bar charts are in linear scale, and min and max values are placed on the corresponding bars. For the empirical networks, the measure of the real time is also reported.</figcaption></figure>

In this case, the critical transmission probability $\lambda_c$ rapidly saturates to a constant value (Fig. 3), and an even more rapid saturation is observed for the other two network models (see Ref. [56]). The average degree is indeed relatively stable, so small temporal windows are enough to fully characterize $\lambda_c$.

Different values of the recovery probability $\mu$ lead, in general, to similar behaviors of $\lambda_c$ towards saturation, differing essentially by a scaling factor. The effect of $\mu$ on saturation time is instead visible for the BURSTY network. In this case, when the period length is smaller than the average duration of the infection, the truncation in the intercontact time distribution clearly alters the estimation of the threshold.

These results indicate that it is possible to identify a minimum length of the observation window of a real system for contact data collection, highlighting the presence of well-defined properties and patterns characterizing the system that can be captured in a finite time.

## V. CONCLUSION

Being able to provide a reliable and accurate estimation of the epidemic threshold for a spreading process taking place on a given networked system is of the utmost importance, as it allows predictions of the likelihood of a wide-spreading event and identification of containment measures (crucial for infectious disease epidemics) or strategies for enhancing the propagation (desired in the

case of information diffusion). While analytical approaches have so far targeted only specific contexts, our framework allows the analytical computation of the epidemic threshold on an arbitrary temporal network, requiring no assumption about the network topology or time variation. The proposed approach is based on the spectral decomposition of the flattened matrix representation of the topological and temporal structure of the network, extending the Markov chain model introduced for the static network case to its temporal counterpart. The predicted epidemic threshold, validated against the numerical solution of the model, also reproduces the behavior observed in stochastic microscopic numerical simulations of the spreading process with high accuracy.

The technical requirement of periodic conditions does not limit the general applicability of our approach, as the method is valid for an arbitrary period length. Moreover, this feature allows us to inform data-collection endeavors on the time period of observation of the system required to fully characterize its spreading properties. Our focus on the discrete time formulation of the process is prompted by the study of several empirical networks for which time is naturally discrete, the time step being dictated by the resolution of the data-collection procedure. Extensions to the continuous-time case would be needed when the continuous-time description is more appropriate, as for example with some modeling approaches.
Our framework thus introduces a multilayer formulation of spreading phenomena on time-varying networks that opens the path to new theoretical understandings of the complex interplay between the two temporal processes, disentangling the role of the network's dynamical features, such as activation rate, temporal correlations, and temporal resolution.

## ACKNOWLEDGMENTS

We thank Ciro Cattuto, Mario Giacobini, Sandro Meloni, Yamir Moreno, Nicola Perra, Mason Porter, and Alessandro Vespignani for fruitful discussions. This work has been partially funded by the EC-Health Contract No. 278433 (PREDEMICS) to V. C. and C. P.; the ANR Contract No. ANR-12-MONU-0018 (HARMSFLU) to V. C.; the PHC Programme Galilee Contract No. 28144NH, the Italian Ministry of Health Contract No. IZS AM 04/11 RC, the EC-ANIHWA Contract No. ANR-13-ANWA-0007-03 (LIVEepi) to E. V., C. P., V. C.; the "Pierre Louis" School of Public Health of UPMC, Paris, France to E. V.; and the Lagrange Project on Complex Systems-CRT and ISI Foundation to L. F.

## APPENDIX A: PROOF OF EQ. (8)

Computing the eigenvalues of $M^\dagger$ means solving the equation $\det(x - M^\dagger) = 0$, where the determinant is computed on the $\mathbb{R}^{NT}$ space ($\det_{NT}$). Given that $x - M^\dagger$ is composed of $T^2$ blocks of size $N \times N$, we can use the findings in Ref. [57] to reduce the dimensionality of the problem, i.e., $\det_{NT} \rightarrow \det_N$. Moreover, given that several blocks of $x - M^\dagger$ are zero, the general result in Ref. [57] simplifies to $\det_{NT}(x - M^\dagger) = (-1)^{NT}\det_N(x^T - P)$. Equation (8) immediately follows.

## APPENDIX B: NETWORKS CONSIDERED FOR VALIDATION

In this section, we provide the details of the six networks used for the validation of the threshold expression.

**ER:** The network is formed by a sequence of random Erdős-Rényi graphs [50] with 500 nodes and 750 edges, so $\langle k \rangle = 3$.

**ACTIVITY:** In the activity-driven model by Perra et al. [28], nodes that are active in a certain snapshot establish a fixed number of connections (in our case, 2) with other nodes picked at random (both active and inactive). All links are renewed after every snapshot. The activity potential is assigned by sampling numbers $x \in [\epsilon, 1]$ from a power-law distribution (in our case, with exponent 2.8 and $\epsilon = 3 \times 10^{-2}$) and then converting them to activity potentials $a = 1 - e^{-\eta x}$. $\eta$ is a free parameter used to tune the average degree; here, $\eta = 10$.

**BURSTY:** The network is built from the model introduced by Rocha et al. [23], where a node becomes active at time $t$ with a probability that depends on the time it was last

active, $t'$. If active, it then forms a link with another active node. All links are removed before proceeding to the following snapshots. To enforce a BURSTY interevent time distribution, the probability of becoming active is sampled with the distribution $(t - t')^{-\alpha_1}e^{-\alpha_2(t-t')}$, with $\alpha_1 = 2$ and $\alpha_2 = 5 \times 10^{-4}$.

**HT09:** The data set was collected by the Sociopatterns group [16], and it records the interactions among participants at a scientific conference. Links represent face-to-face proximity recorded by wearable radio-frequency identification (RFID) tags. Time resolution of the signal is 20 sec. Each tag emits wave packets that may be recorded by other tags, thus signaling proximity. Tags were embedded in conference badges; their signal intensity was set to be detectable at a maximum distance of 1 m and completely shielded by the human body. With these settings in effect, only close proximity in a face-to-face position resulted in interaction [58].

**SEX:** This is a network of sexual contacts between female prostitutes and their male clients as inferred from posts on a Brazilian online escort forum where customers could rate their experience with a certain sex worker. The date of the post was taken as a proxy for the time of the sexual intercourse. The (anonymized) data set can be found in Ref. [18].

**SCHOOL:** This network represents the face-to-face proximity interactions during a day in a high school [51]. Students and staff were given wearable RFID sensors, and proximity was recorded in a similar fashion as for HT09 [59].

## APPENDIX C: ESTIMATION OF THE EPIDEMIC THRESHOLD FROM NUMERICAL SIMULATIONS

The computation of $\langle i_{sim} \rangle$ in proximity to the transition is difficult because surviving configurations are rare and a very large number of realizations of the process is needed to collect substantial statistics. We use the quasistationary state (QS) method [52,60] to overcome this difficulty and increase our computational efficiency. The QS method is based on the idea of constraining the system in an active state. Every time the absorbing state (i.e., no infected) is reached by the system, it is substituted with an active configuration that is randomly taken from the history of the simulation. In particular, 50 active configurations for each network snapshot are kept in memory. Whenever an active configuration is reached, it replaces one of the 50 with probability 0.2. When the absorbing state is reached, an active configuration is chosen among these 50 of that particular snapshot. For each simulation, after a relaxation time of $3 \times 10^3$ time steps, statistics are collected during $10^5$ time steps. The method produces a time series that is long enough to accurately compute the observables $\langle i_{sim} \rangle$.
[1] M. J. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals (Princeton University Press, Princeton, 2008).

[2] V. Colizza, A. Barrat, M. Barthelemy, A.-J. Valleron, and A. Vespignani, Modeling the Worldwide Spread of Pandemic Influenza: Baseline Case and Containment Interventions, PLoS Med. 4, e13 (2007).

[3] A. J. Tatem and D. L. Smith, International Population Movements and Regional Plasmodium Falciparum Malaria Elimination Strategies, Proc. Natl. Acad. Sci. U.S.A. 107, 12222 (2010).

[4] R. Pastor-Satorras and A. Vespignani, Epidemic Spreading in Scale-Free Networks, Phys. Rev. Lett. 86, 3200 (2001).

[5] W. Goffman and V. A. Newill, Generalization of Epidemic Theory: An Application to the Transmission of Ideas, Nature (London) 204, 225 (1964).

[6] D. J. Daley and D. G. Kendall, Epidemics and Rumours, Nature (London) 204, 1118 (1964).

[7] S. Bikhchandani, D. Hirshleifer, and I. Welch, A Theory of Fads, Fashion, Custom, and Cultural Change in Informational Cascades, J. Polit. Econ. 100, 992 (1992).

[8] Y. Wang, D. Chakrabarti, C. Wang, and C. Faloutsos, Epidemic Spreading in Real Networks: An Eigenvalue Viewpoint, in Proceedings of the 22nd International Symposium on Reliable Distributed Systems (IEEE, Los Alamitos, CA, 2003), pp. 25-34.

[9] S. Gómez, A. Arenas, J. Borge-Holthoefer, S. Meloni, and Y. Moreno, Discrete Time Markov Chain Approach to Contact-Based Disease Spreading in Complex Networks, Europhys. Lett. 89, 38009 (2010).

[10] R. Cohen, K. Erez, D. ben Avraham, and S. Havlin, Resilience of the Internet to Random Breakdowns, Phys. Rev. Lett. 85, 4626 (2000).

[11] M. E. J. Newman, Spread of Epidemic Disease on Networks, Phys. Rev. E 66, 016128 (2002).

[12] M. Boguñá, C. Castellano, and R. Pastor-Satorras, Nature of the Epidemic Threshold for the Susceptible-Infected-Susceptible Dynamics in Networks, Phys. Rev. Lett. 111, 068701 (2013).

[13] C. Castellano and R. Pastor-Satorras, Thresholds for Epidemic Spreading in Networks, Phys. Rev. Lett. 105, 218701 (2010).

[14] A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Localization and Spreading of Diseases in Complex Networks, Phys. Rev. Lett. 109, 128702 (2012).

[15] A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge, England, 2008).

[16] L. Isella, J. Stehlé, A. Barrat, C. Cattuto, J.-F. Pinton, and W. Van den Broeck, What's in a Crowd? Analysis of Face-to-Face Behavioral Networks, J. Theor. Biol. 271, 166 (2011).

[17] G. Miritello, E. Moro, and R. Lara, Dynamical Strength of Social Ties in Information Spreading, Phys. Rev. E 83, 045102 (2011).

[18] L. E. C. Rocha, F. Liljeros, and P. Holme, Information Dynamics Shape the Sexual Networks of Internet-Mediated Prostitution, Proc. Natl. Acad. Sci. U.S.A. 107, 5706 (2010).

[19] P. Bajardi, A. Barrat, F. Natale, L. Savini, and V. Colizza, Dynamical Patterns of Cattle Trade Movements, PLoS One 6, e19869 (2011).

[20] C. T. Butts, Revisiting the Foundations of Network Analysis, Science 325, 414 (2009).

[21] M. Karsai, M. Kivelä, R. K. Pan, K. Kaski, J. Kertész, A.-L. Barabási, and J. Saramäki, Small but Slow World: How Network Topology and Burstiness Slow Down Spreading, Phys. Rev. E 83, 025102 (2011).

[22] P. Holme and J. Saramäki, Temporal Networks, Phys. Rep. 519, 97 (2012).

[23] L. E. C. Rocha and V. D. Blondel, Bursts of Vertex Activation and Epidemics in Evolving Networks, PLoS Comput. Biol. 9, e1002974 (2013).

[24] L. Ferreri, P. Bajardi, M. Giacobini, S. Perazzo, and E. Venturino, Interplay of Network Dynamics and Heterogeneity of Ties on Spreading Dynamics, Phys. Rev. E 90, 012812 (2014).

[25] J. Luis Iribarren and E. Moro, Impact of Human Activity Patterns on the Dynamics of Information Diffusion, Phys. Rev. Lett. 103, 038702 (2009).

[26] A. Vazquez, B. Rácz, A. Lukács, and A.-L. Barabási, Impact of Non-Poissonian Activity Patterns on Spreading Processes, Phys. Rev. Lett. 98, 158702 (2007).

[27] K. T. D. Eames and M. J. Keeling, Monogamous Networks and the Spread of Sexually Transmitted Diseases, Math. Biosci. 189, 115 (2004).

[28] N. Perra, B. Gonçalves, R. Pastor-Satorras, and A. Vespignani, Activity Driven Modeling of Time Varying Networks, Sci. Rep. 2, 469 (2012).

[29] E. Volz and L. Ancel Meyers, Epidemic Thresholds in Dynamic Contact Networks, J. R. Soc. Interface 6, 233 (2009).

[30] T. Gross, C. J. Dommar D'Lima, and B. Blasius, Epidemic Dynamics on an Adaptive Network, Phys. Rev. Lett. 96, 208701 (2006).

[31] M. Taylor, T. J. Taylor, and I. Z. Kiss, Epidemic Threshold and Control in a Dynamic Network, Phys. Rev. E 85, 016103 (2012).

[32] S. Liu, N. Perra, M. Karsai, and A. Vespignani, Controlling Contagion Processes in Activity Driven Networks, Phys. Rev. Lett. 112, 118702 (2014).

[33] M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, Mathematical Formulation of Multilayer Networks, Phys. Rev. X 3, 041022 (2013).

[34] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, Multilayer Networks, J. Comp. Networks 2(3), 203 (2014).

[35] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevX.5.021005 for an extension to the weighted and directed cases.

[36] P. Van Mieghem, The N-Intertwined SIS Epidemic Network Model, Computing 93, 147 (2011).

[37] J. P. Gleeson, S. Melnik, J. A. Ward, M. A. Porter, and P. J. Mucha, Accuracy of Mean-Field Theory for Dynamics on Real-World Networks, Phys. Rev. E 85, 026106 (2012).

[38] P. Van Mieghem and R. van de Bovenkamp, Non-Markovian Infection Spread Dramatically Alters the
*Susceptible-Infected-Susceptible Epidemic Threshold in Networks*, Phys. Rev. Lett. **110**, 108701 (2013).

[39] C. Li, H. Wang, and P. Van Mieghem, *Epidemic Threshold in Directed Networks*, Phys. Rev. E **88**, 062802 (2013).

[40] E. Cator and P. Van Mieghem, *Second-Order Mean-Field Susceptible-Infected-Susceptible Epidemic Threshold*, Phys. Rev. E **85**, 056111 (2012).

[41] P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, *Community Structure in Time-Dependent, Multiscale, and Multiplex Networks*, Science **329**, 277 (2010).

[42] K. Wehmuth, A. Ziviani, and E. Fleury, *A Unifying Model for Representing Time-Varying Graphs*, arXiv: 1402.3488.

[43] E. Cozzo, R. A. Baños, S. Meloni, and Y. Moreno, *Contact-Based Social Contagion in Multiplex Networks*, Phys. Rev. E **88**, 050801 (2013).

[44] H. Wang, Q. Li, G. D'Agostino, S. Havlin, H. Eugene Stanley, and P. Van Mieghem, *Effect of the Interconnected Network Structure on the Epidemic Threshold*, Phys. Rev. E **88**, 022801 (2013).

[45] C. Granell, S. Gómez, and A. Arenas, *Dynamical Interplay between Awareness and Epidemic Spreading in Multiplex Networks*, Phys. Rev. Lett. **111**, 128701 (2013).

[46] S. Elaydi, *An Introduction to Difference Equations*, 3rd ed. (Springer, New York, NY, 2005).

[47] H. H. K. Lentz, T. Selhorst, and I. M. Sokolov, *Unfolding Accessibility Provides a Macroscopic Approach to Temporal Networks*, Phys. Rev. Lett. **110**, 118701 (2013).

[48] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.5.021005 for performance of the algorithm used to compute the threshold.

[49] N. Masuda, K. Klemm, and V. M. Eguiluz, *Temporal Networks: Slowing Down Diffusion by Long Lasting Interactions*, Phys. Rev. Lett. **111**, 188701 (2013).

[50] P. Erdős and A. Rényi, *On the Evolution of Random Graphs*, Publ. Math. Inst. Hungar. Acad. Sci. **5**, 17 (1960).

[51] M. Salathé, M. Kazandjieva, J. W. Lee, P. Levis, M. W. Feldman, and J. H. Jones, *A High-Resolution Human Contact Network for Infectious Disease Transmission*, Proc. Natl. Acad. Sci. U.S.A. **107**, 22020 (2010).

[52] S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, *Epidemic Thresholds of the Susceptible-Infected-Susceptible Model on Networks: A Comparison of Numerical and Theoretical Results.*, Phys. Rev. E **86**, 041125 (2012).

[53] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.5.021005 for an analysis of dynamical correlations.

[54] R. K. Pan and J. Saramäki, *Path Lengths, Correlations, and Centrality in Temporal Networks*, Phys. Rev. E **84**, 016105 (2011).

[55] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.5.021005 for further plots showing the impact of average degree variations on threshold fluctuations.

[56] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.5.021005 for the saturation of the threshold in the other network models.

[57] P. D. Powell, *Calculating Determinants of Block Matrices*, arXiv:1112.4379.

[58] Data sets and additional information are available on www.sociopatterns.org.

[59] Data can be found on www.salathegroup.com.

[60] K. Binder and D. W. Heermann, *Monte Carlo Simulation in Statistical Physics, An Introduction* (Springer, Berlin, 2010).