key: cord-024742-hc443akd
authors: Liu, Quan-Hui; Xiong, Xinyue; Zhang, Qian; Perra, Nicola
title: Epidemic spreading on time-varying multiplex networks
date: 2018-12-03
journal: nan
DOI: 10.1103/physreve.98.062303
sha: 
doc_id: 24742
cord_uid: hc443akd

Social interactions are stratified in multiple contexts and are subject to complex temporal dynamics. The systematic study of these two features of social systems has started only very recently, mainly thanks to the development of multiplex and time-varying networks. However, these two advancements have progressed almost in parallel with very little overlap. Thus, the interplay between multiplexity and the temporal nature of connectivity patterns is poorly understood. Here, we aim to tackle this limitation by introducing a time-varying model of multiplex networks. We are interested in characterizing how these two properties affect contagion processes. To this end, we study susceptible-infected-susceptible epidemic models unfolding at comparable timescale with respect to the evolution of the multiplex network. We study both analytically and numerically the epidemic threshold as a function of the multiplexity and the features of each layer. We found that higher values of multiplexity significantly reduce the epidemic threshold especially when the temporal activation patterns of nodes present on multiple layers are positively correlated. Furthermore, when the average connectivity across layers is very different, the contagion dynamics is driven by the features of the more densely connected layer. Here, the epidemic threshold is equivalent to that of a single layered graph and the impact of the disease, in the layer driving the contagion, is independent of the multiplexity. However, this is not the case in the other layers where the spreading dynamics is sharply influenced by it. The results presented provide another step towards the characterization of the properties of real networks and their effects on contagion phenomena.

Social interactions take place in different contexts and modes of communication. On a daily basis we interact at work, in the family, and across a wide range of online platforms or tools, e.g., Facebook, Twitter, emails, mobile phones, etc. In the language of modern network science, social networks can be conveniently modeled and described as multilayer networks [1] [2] [3] [4] [5] . This is not a new idea. Indeed, the intuition that social interactions are stratified in different layers dates back several decades [6] [7] [8] . However, the digitalization of our communications and the miniaturization of devices has just recently provided the data necessary to observe, at scale, and characterize the multilayer nature of social interactions.

As in the study of single layered networks, the research on multilayer graphs is divided in two interconnected areas. The first deals with the characterization of the structural properties of such entities [1, 4] . One of the central observations is that the complex topology describing each type of interaction (i.e., each layer) might be different. Indeed, the set and intensity of interactions in different contexts (e.g., work, family, etc.) or platforms (e.g., Facebook, Twitter, etc.) is not the same. Nevertheless, layers are coupled by individuals active across two or more of them. The presence of such coupling as well as its * n.perra@greenwich.ac.uk degree is often referred to as multiplexity. Another interesting feature of multilayer graphs is that the connectivity patterns in different layers might be topologically and temporally correlated [9] [10] [11] . The second area of research instead considers the function, such as sustaining diffusion or contagion processes, of multilayer networks [1, 12, 13] . A large fraction of this research aims at characterizing how the complex structural properties of multilayer graphs affect dynamical processes unfolding on their fabric. The first important observation is that disentangling connections in different layers gives rise to complex and highly nontrivial dynamics function of the interplay between interlayer and intralayer connections [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] . A complete summary of the main results in the literature is beyond the scope of the paper. We refer the interested reader to these recent resources for details [1] [2] [3] [4] [5] .

Despite the incredible growth of this area of network science over the last years, one particular aspect of multilayer networks is still largely unexplored: the interplay between multiplexity and the temporal nature of the connectivity patterns especially when dynamical processes unfolding on their fabric are concerned [13] . This should not come as a surprise. Indeed, the systematic study of the temporal dynamics even in single layered graphs is very recent. In fact, the literature has been mostly focused on time-integrated properties of networks [27, 28] . As result, complex temporal dynamics acting at shorter timescales have been traditionally discarded. However, the recent technological advances in data storing and collection are providing unprecedented means to probe also the temporal dimension of real systems. The access to this feature is allowing to discover properties of social acts invisible in time aggregated data sets, and is helping characterize the microscopic mechanisms driving their dynamics at all timescales [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] . The advances in this arena are allowing to investigate the effects such temporal dynamics have on dynamical processes unfolding on time-varying networks. The study of the propagation of infectious diseases, ideas, rumors, or memes, etc., on temporal graphs shows a rich and nontrivial phenomenology radically different than what is observed on their static or annealed counterparts [29, .

Before going any further, it is important to notice how in their more general form, multilayer networks might be characterized by different types of nodes in each layer. For example, modern transportation systems in cities can be characterized as a multilayer network in which each layer captures a different transportation mode (tube, bus, public bikes, etc.) and the links between layers connect stations (nodes) where people can switch mode of transport [4, 12] . A particular version of multilayer networks, called multiplex, is typically used in social networks. Here, the entities in each layer are of the same type (i.e., people). The interlayer links are drawn only to connect the same person in different layers.

In this context, we introduce a model of time-varying multiplex networks. We aim to characterize the effects of temporal connectivity patterns and multiplexity on contagion processes. We model the intralayer evolution of connections using the activity-driven framework [29] . In this model of time-varying networks, nodes are assigned with an activity describing their propensity to engage in social interactions per unit time [29] . Once active, a node selects a partner to interact. Several selection mechanisms have been proposed, capturing different features of real social networks [30, 31, [62] [63] [64] . The simplest, that will be used here, is memoryless and random [29] . In these settings, active nodes do not have a preference towards individuals previously connected. Despite such simplification, the mechanism allows capturing the heterogeneity of time-integrated connectivity patterns observed in real networks while guaranteeing mathematical tractability [29] . The multiplexity or coupling between layers is modulated by a probability p. If p = 1 all nodes are present in all layers. If p = 0, the multiplex is formed by M disconnected graphs. We consider p as a parameter and explore different regime of coupling between layers. Furthermore, each layer is characterized by an activity distribution. We consider different scenarios in which the activity of coupling nodes, which are present in different layers (regulated by p), is uncorrelated as well as others in which it is instead positively or negatively correlated. In these settings, we study the unfolding of susceptible-infected-susceptible (SIS) epidemic processes [65] [66] [67] . We derive analytically the epidemic threshold for two layers for any p and any distributions of activities. In the limit of p = 1 we find analytically the epidemic threshold for any number of layers. Interestingly, the threshold is a not trivial generalization of the correspondent quantity in the monoplex (single layer network). In the general case 0 < p < 1 we found that the threshold is a decreasing function of p. Positive correlations of coupling nodes push the threshold to smaller values with respect to the uncorrelated and negatively correlated cases. Furthermore, when the average connectivity of two layers is very different, the critical behavior of the system is driven by the more densely connected layer. In such a scenario the epidemic threshold is not affected by the multiplexity, its value is equivalent to the case of a monoplex, and the coupling affects only the layer featuring the smaller average connectivity.

The paper is organized as follows. In Sec. II we introduce the multiplex model. In Sec. III we study first both analytically and numerically the spreading of SIS processes. Finally, in Sec. IV we discuss our conclusions.

We first introduce the multiplex model. For simplicity of discussion, we consider the case in which the system is characterized by M = 2 layers A and B. However, the same approach can be used to create a multiplex with any number of layers. Let us define N as the number of nodes in each layer. In general, we have three different categories of nodes: N A , N B , and N o . They describe, respectively, the number of nodes that are present only in layer A, B, or in both. The last category is defined by a parameter p: the coupling between layers (multiplexity). Thus, on average, we have N A = N B = (1 − p)N and N o = pN. As mentioned in the introduction, the temporal dynamics in each layer is defined by the activity-driven framework [29] . Thus, each noncoupling node is characterized by an activity extracted from a distribution f A (a) or f B (a) which captures its propensity to be engaged in a social interaction per unit time. Observations in real networks show that the activity typically follows a heavy-tailed distribution [29] [30] [31] 41, 62, 68] . Here, we assume that activities follow power laws, thus, f x (a) = c x a −γ x with x = [A, B] and a 1 to avoid divergences. Coupling nodes instead are characterized by a joint activity distribution h(a A , a B ). As mentioned in the introduction, real multiplex networks are characterized by correlations across layers. In particular, the study of a wide range of real systems shows a complex and case dependent phenomenology in which the topological features (i.e., static connectivity patterns) of coupling nodes can be either positively or negatively correlated [9] . Furthermore, researchers found evidence of positive temporal correlations between the activation patterns across layers [10, 11] . To account for such observations and explore their effects on spreading processes, we consider three simple prototypical cases in which the activities of coupling nodes in the two layers are (i) uncorrelated, or (ii) positively and (iii) negatively correlated. To simplify the formulation and to avoid adding other parameters, in case of positive and negative correlations we adopt the following steps. We first extract N o activities from the two distributions f x (a). Then, we order them. In the case of positive correlation, a node that has the rth activity in A will be assigned to the correspondent activity in B. In other words, the first node will be assigned to the highest activity extracted from f A (a) and the highest value extracted from f B (a). The second will be assigned to the second highest activity extracted from both distributions, etc. In the case of negative correlations instead, a node that has the rth activity in A will be assigned the (pN − r + 1)th in B. In other words, the first node will be assigned to the highest activity in A and to the lowest activity in B. The second node will be assigned to the second highest activity in A and to the second lowest in B, etc.

In these settings, the temporal evolution of the multiplex network is defined as follows. For each realization, we randomly select pN nodes as coupling nodes between layers. At each time step t, note the following:

(i) Each node is active with a probability defined by its activity.

(ii) Each active node creates m x links with randomly selected nodes. Multiple connections with the same node in the same layer within each time step are not allowed.

(iii) Coupling nodes can be active and create connections in both layers.

(iv) At time step t + t all connections are deleted and the process restarts from the first point.

All connections have the same duration of t. In the following, we set, without lack of generality, t = 1. At each time the topology within each layer is characterized, mostly, by a set of disconnected stars of size m x + 1. Thus, at the minimal temporal resolution each network looks very different than the static or annealed graphs we are used to seeing in the literature [69] . However, it is possible to show that, integrating links over T time steps in the limit in which T N, the resulting network has a degree distribution that follows the activity [29, 31, 70] . In other words, the heterogeneities in the activity distribution translate in heterogenous timeaggregated connectivity patterns typically of real networks. Thus, as observed in real temporal networks the topological features at different timescales are very different than the late (or time-integrated) characteristics [38] .

At each time step the average degree in each layer can be computed as

where E x t is the number of links generated in each layer at each time step. Furthermore, a x = da f x (a)a and a x o = da a da B h(a A , a B )a x are the average activity of noncoupling and coupling nodes in each layer, respectively. Similarly, the total average degree (often called overlapping degree [71] ), at each time step, is

Thus, the average connectivity, at each time step, is determined by the number of links created in each layer, and by the interplay between the average activity of coupling and noncoupling nodes. As shown in Fig. 1 (a), Eq. (2) describes quite well the behavior of the average overlapping degree which is an increasing function of the multiplexity p. Indeed, the larger the fraction of coupling nodes, the larger the connectivity of such nodes across layers. As we will see in the next section, this feature affects significantly the unfolding on contagion processes. In Fig. 1 (b) we show the integrated degree distribution of the overlapping degree for different p. The plot clearly shows how the functional form is defined by the activity distributions of the two layers which in this case are equal. An increase in the fraction of coupling nodes does not change the distribution of the overlapping degree; it introduces a vertical shift which, however, is more visible for certain values of k.

In order to understand how the interplay between multiplexity and temporal connectivity patterns affects dynamical processes, we consider SIS contagion phenomena spreading on the multiplex model introduced in the previous section. In this prototypical epidemic model each node can be in one of two compartments. Healthy nodes are susceptible to the disease and thus in the compartment S. Infectious nodes instead join the compartment I . The natural history of the disease is defined as follows. A susceptible, in contact with an infected node, might get sick and infectious with probability λ. Each infected node spontaneously recovers with rate μ, thus staying infectious for μ −1 time steps, on average. One crucial feature of epidemic models is the threshold which determines the conditions above which a disease is able to affect a macroscopic fraction of the population [65] [66] [67] . In case of SIS models, below the threshold the disease dies out, reaching the so called disease-free equilibrium. Above threshold instead, the epidemic reaches an endemic stationary state. This can be captured running the simulations for longer times and thus estimating the fraction of infected nodes for t → ∞: i ∞ . In general, in a multiplex network, such fraction might be different across layers. Thus, we can define i x ∞ .

To characterize the threshold we could study the behavior of such fraction(s) as function of λ/μ. Indeed, the final number of infected nodes acts as order parameter in a second order phase transition, thus defining the critical conditions for the spreading [66] . However, due to the stochastic nature of the process, the numerical estimation of the endemic state, especially in proximity of the threshold, is not easy. Thus, we adopt another method measuring the lifetime of the process L [72] . This quantity is defined as the average time the disease needs either to die out or to infect a macroscopic fraction Y of the population. The lifetime acts as the susceptibility in phase transitions thus allows a more precise numerical estimation [72] .

In the case of single layer activity-driven networks, in which partners of interactions are chosen at random and without memory of past connections, the threshold can be written as (see Ref. [29] for details)

Thus, the conditions necessary for the spread of the disease are set by the interplay between the features of the disease (left side) and the dynamical properties of the time-varying networks where the contagion unfolds (right side). The latter are regulated by first and second moments of the activity distribution and by the number of connections created by each active node (i.e., m). It is important to notice that Eq. (3) considers the case in which the timescale describing the evolution of the connectivity patterns and the epidemic process is comparable. The contagion process is unfolding on a time-varying network. In the case when links are integrated over time and the SIS process spreads on a static or annealed version of the graph, the epidemic threshold will be much smaller [29, 73, 74] . This is due to the concurrency of connections which favors the spreading. In fact, by aggregating the connections over time, the degree of each node increases thus facilitating the unfolding of the disease. In this limit of timescale separation between the dynamics of and on networks, the evolution of the connectivity patterns is considered either much slower (static case) or much faster (annealed case) with respect to the epidemic process. In the following, we will only consider the case of comparable timescales. This is the regime of time-varying networks which is extremely relevant for a variety of spreading processes ranging from sexual transmitted diseases and influenzalike illnesses to rumors and information propagation [28] .

What is the threshold in the case of our multiplex and timevarying network model? In the limit p = 0 the number of coupling nodes is zero. The two layers are disconnected, thus, the system is characterized by two independent thresholds regulated by the activity distributions of the two layers. The most interesting question is then what happens for p > 0.

To find an answer to this conundrum, let us define I , respectively, as the number of noncoupling nodes of activity a and as the number of coupling nodes of activity a A and a B . The implicit assumption we are making by dividing nodes according to their activities is that of statistical equivalence within activity classes [66, 75] . In these settings, we can write the variation of the number of infected noncoupling nodes as function of time as

where we omitted the dependence of time. The first term on the right-hand side considers nodes recovering, thus leaving the infectious compartment. The second and third terms account for the activation of susceptibles in activity class a (S x a = N x a − I x a ) that select infected nodes (noncoupling and coupling) as partners and get infected. The last two terms instead consider the opposite: infected nodes activate, select as partners noncoupling and coupling nodes in the activity class a infecting them as a result.

Similarly, we can write the expression for the variation of coupling nodes of activity classes a A and a B as

The general structure of the equation is similar to the one we wrote above. The main difference is, however, that coupling nodes can be infected and can infect in both layers. The first term in the right-hand side accounts for the recovery process. The next four (two for each element in the sum in y) consider the activation of susceptible nodes that select as partners both noncoupling and coupling infected nodes and get infected. The last four terms account for the reverse process. In order to compute the epidemic threshold we need to define four auxiliary functions, thus defining a closed system of differential equations. In particular, we define x = da I x a a and o x = da A da B I o a A ,a B a x . For simplicity, we will skip the detailed derivation here (see the Appendix for the details). By manipulating the previous three differential equations we can obtain four more, one for each auxiliary function. The condition for the spreading of the disease can be obtained by studying the spectral properties of the Jacobian matrix of such system of seven differential equations. In particular, if the largest eigenvalue of the Jacobian matrix is larger than zero, the system of equations will not be stable and, consequently, the number of infected nodes will increase. Thus, the epidemic threshold can be obtained by studying the conditions for which this holds.

As sanity check, let us consider first the limit p = 0. In this case, each layer acts independently and we expect the threshold of each to follow Eq. (3). This is exactly what we find. In particular, the system of equations can be de-coupled in two different subsets (one per layer) which are governed by two Jacobian matrices whose largest eigenvalues are

where a n x = da f x (a)a n . Thus, the spreading process will be able to affect a finite fraction of the total population in case either of these two eigenvalues is larger than zero, which implies λ μ > (m x a x + m x a 2 x ) −1 as expected. It is important to notice that in case of a multiplex network the disease might be able to spread in one layer but not in the other. However, in case the condition for the spreading is respected in both layers, they will experience the disease.

Let us consider the opposite limit: p = 1. As described in details in the Appendix, the condition for the spreading of the disease reads as

where a n

Interestingly, the threshold is a function of the first and second moments of the activity distributions of the coupling nodes which are modulated by the number of links each active node creates, plus a term which encodes the correlation of the activities of such nodes in the two layers.

Before showing the numerical simulations to validate the mathematical formulation, an important observation is in order. In this limit, effectively, we could think the multiplex as a multigraph: a single layer network with two types of edges. In case the joint probability distribution of activity is h(a A , a B ) = f (a A )δ(a B − a A ), thus two activities are exactly the same, and m A = m B the threshold reduces to Eq. (3) (valid for a single layer network) in which the number of links created by active nodes is 2m. However, for a general form of the joint distribution and in case of different number of links created by each active node in different layers, this correspondence breaks down.

In all the following simulations, we set N = 10 5 , = 10 −3 , μ = 0.015, Y = 0.3, start the epidemic process from a 1% of nodes selected randomly as initial seeds, and show the averages of 10 2 independent simulations. In Fig. 2 we show the first results considering a simple scenario in which m A = m B = 1 and the exponents for the distributions of activities are the same γ x = 2.1. The first observation is that in all three cases the analytical solutions (vertical dotted lines) agree with the results from simulations. The second observation is that in case of positive correlation between the activities of nodes in two layers, the threshold is significantly smaller than in the other two cases. This is not surprising as the nodes sustaining the spreading in both layers are the same. Thus, effectively, active nodes are capable to infect the double number of other nodes. As we mentioned above, many real multiplexes are characterized by different types of positive correlations. When thinking at real outbreaks, the effect of such feature on the spreading process suggests quite a worrying scenario. However, it is important to remember that real multiplexes are sparse, thus characterized by values of multiplexity which are far from the limit p = 1 [9] . As we will see below, this aspect plays a crucial role in the more realistic cases of 0 < p < 1. The thresholds of the uncorrelated and negatively correlated cases are very similar. In fact, due to the heterogeneous nature of the activity distributions, except for few nodes in the tails, the effective difference between the activities matched in reverse or random order is not large, for the majority of nodes. In Figs. 2(b) and 2(c) we show the behavior of the threshold as function of the activity exponents and the number of links created by active nodes in the two layers. For a given distribution of activity in a layer, increasing the exponent in the other (thus reducing the heterogeneity in the activity distribution) results in an increase of the threshold. This is due to the change of the first and second moments which decrease as a result of the reduced heterogeneity. In the settings considered here, if both exponents of activity distributions are larger than 2.6, the critical value of λ becomes larger than 1, as shown in Fig. 2(b) . Thus, in such region of parameters, the disease will not be able to spread. For a given number of links created in a layer by each active node, increasing the links created in the other layer results in a quite rapid reduction of the threshold. This is due to the increase of the connectivity and thus the spreading potential of active nodes.

In the limit p = 1, we are able to obtain an expression for the threshold of an SIS process unfolding on M layers. It is important to stress that this scenario is rather unrealistic. Indeed, in real multiplexes the majority of nodes is present only in one or two layers [9] . Nevertheless, we can argue that understanding the behavior of the threshold also in this case is of theoretical interest. With this observation in mind, the analytical condition for the spreading of the disease can be written as (see the Appendix for details) 

where x = [A, B, . . . , Z] and z > y implies an alphabetical ordering. The first observation is that in case h(a y , a z ) = f (a y )δ(a z − a y ) ∀ y, z ∈ x, thus the activity is the same for each node across each layer, Eq. (8) reduces to

which is the threshold for a single layer activity-driven network in which m → Mm. This is the generalization of the correspondence between the two thresholds we discussed above for two layers. The second observation is that, in general, increasing the number of layers decreases the epidemic threshold. Indeed, each new layer increases the connectivity potential of each node and thus the fragility of the system to the contagion process. Figure 3(a) shows the analytical behavior of the epidemic threshold up to M = 10 for the simplest case of uncorrelated and positively correlated activities between layers confirming this result. In Fig. 3(b) we show the comparison between the analytical results and the numerical simulations. The plot shows a perfect match between the two. Furthermore, the two plots confirm the effects of positive correlations which facilitate the spreading of the disease.

We now turn the attention to the most interesting and realistic cases which are different from the two limits of null and total coupling of nodes considered above. For a general value of p, we could not find a general closed expression for the epidemic threshold. However, the condition for the spreading can be obtained by investigating, numerically, the spectral properties of the Jacobian (see the Appendix for details). In Fig. 4 we show the lifetime of SIS spreading processes unfolding on a multiplex network for three different values of p. Figure 4(a) shows the uncorrelated case and the dashed vertical lines describe the analytical predictions. The first observation is that the larger the multiplexity between two layers, the smaller the threshold. This should not come as a surprise. In fact, as shown previously in Fig. 1 , the average connectivity in the system increases as a function of p. Thus, increasing the fraction of nodes active in both layers increases the spreading power of such nodes when they get infected. The second observation is that the analytical predictions match remarkably well the simulations. The bottom panel shows instead the case of positive correlation between the activities of coupling nodes in the two layers. Also in this case, the larger the multiplexity the smaller the epidemic threshold. The comparison between the two panels highlights It is important to notice that also here our analytical predictions match remarkably well the numerical simulations.

In Fig. 5 we show the behavior of the (analytical) epidemic threshold as function of p for three types of correlations.

The results confirm what is discussed above. The larger the multiplexity, the smaller the threshold. Negative and null correlations of coupling nodes exhibit very similar thresholds. Instead, positive correlations push the critical value to smaller values. Furthermore, the smaller the multiplexity, the smaller the effect of positive correlations as the difference between the thresholds increases as function of p. This is a reassuring result. Indeed, as mentioned above, among the properties of real multiplex systems we find the presence of positive topological and temporal correlations as well as low values of multiplexity. The first feature favors the spreading of diseases. Luckily instead, the second property reduces the advantage of positive correlations pushing the threshold to higher values which are closer to the case of negative and null correlations. It is also important to notice how the threshold of a multiplex network (p > 0) is always smaller than the threshold of a monoplex (p = 0) with the same features. Indeed, the presence of coupling nodes effectively increases the spreading potential of the disease, thus reducing the threshold. However, the presence of few coupling nodes (p ∼ 0) does not significantly change the threshold; this result and the effect of multiplexity on the spreading power of diseases is in line with what was already discussed in the literature for static multiplexes [1, 22] .

In Fig. 6 , we show how the epidemic threshold varies when the average connectivity of the two layers is progressively different and asymmetric. In other words, we investigate what happens when one layer has a much larger average connectivity than the other. This situation simulates individuals engaged in two different social contexts, one characterized by fewer interactions (e.g., close family interactions) and one instead by many more connections (e.g., work environment). In the figure, we consider a multiplex network in which the layer A is characterized by m A = 1. We then let m B vary from 1 to 10 and measure the impact of this variation on the epidemic threshold for different values of p. For simplicity, we considered the case of uncorrelated activities in the two layers, but the results qualitatively hold also for the other types of correlations. Few observations are in order. As expected, the case p = 0 is the upper bound of the epidemic threshold. However, the larger the asymmetry between the two layers, thus, the larger the average connectivity in the layer B, the smaller the effect of the multiplexity on the threshold. Indeed, while systems characterized by m B = 1 and higher multiplexity feature a significantly smaller threshold with respect to the monoplex, for m B 3 such differences become progressively negligible and the effects of multiplexity vanish. In this regime, the layer with the largest average connectivity drives the spreading of the disease. The connectivity of layer B effectively determines the dynamics of the contagion, and thus the critical behavior is not influenced by coupling nodes. Interestingly, this result generalizes to the case of time-varying systems which have been found in the case of static multiplexes. Indeed, Cozzo et al. [26] showed how the threshold of contact based processes is driven by the dominant layer which roughly corresponds to the layer featuring higher connectivity.

In order to get a deeper understanding on this phenomenon, we show the asymptotic number of infected nodes in each layer for m A = 1 and m B = 10 in Fig. 7 . For any λ above the threshold the fraction of infected nodes in layer B [ Fig. 7(b) ] is larger than in layer A [ Fig. 7(a) ] and is independent of the fraction of coupling nodes. As discussed above, in these settings the layer B is driving the contagion process and the imbalance between the connectivity patterns is large enough to behave as a monoplex. However, for layer A the contagion process is still highly influenced by p. Indeed, as the fraction of coupling nodes increases, layer A is more and more influenced by the contagion process unfolding in B. Overall, these results are qualitatively similar to the literature of spreading phenomena in static multilayer networks [14] .

We presented a time-varying model of multiplex networks. The intralayer temporal dynamics follows the activity-driven framework which was developed for single layered networks (i.e., monoplexes). Thus, nodes are endowed with an activity that describes their propensity, per unit time, to initiate social interactions. We define the multiplexity as a free parameter p regulating the fraction of coupling nodes between layers. The activities of such nodes are considered, in general, to be different but potentially correlated. In these settings, we studied how multiplexity and temporal connectivity patterns affect dynamical processes unfolding on such systems. To this end, we considered a prototypical model of infectious diseases: the SIS model. We derived analytically the epidemic threshold of the process as function of p. In the limit p = 0, the system is constituted by disconnected networks that behave as monoplexes. In the opposite limit instead (i.e., p = 1) the epidemic threshold is a function of the first and second moments of the activity distributions as well as by their correlations across layers. We found that systems characterized by positive correlations are much more fragile to the spreading of the contagion process with respect to negative and null correlations. As several real multiplex systems feature positive topological and temporal correlations [9] [10] [11] , this result depicts a worrying scenario. Luckily, real multiplexes are also sparse, thus characterized by multiplexity values far from limit p = 1 [9] . The threshold also varies as a function of the number of layers M. Indeed, with perfect coupling each node is present and potentially active in each layer. Thus, the larger M, the smaller the epidemic threshold as the spreading potential of each node increases. In the general and more realistic case 0 < p < 1, we could not find a closed expression for the epidemic threshold. However, the critical conditions for the spreading can be calculated from the theory by investigating numerically the spectral properties of the Jacobian matrix describing the contagion dynamics. Also in this case, positive correlations of activities across layers help the spreading by lowering the epidemic threshold, while negative and null correlations result in very similar thresholds.

Moreover, the lower the multiplexity, the larger the epidemic threshold. Indeed, the case of disconnected monoplexes (i.e., p = 0) is the upper bound for the threshold. Furthermore, the difference between the thresholds in the case of positive and the other two types of correlations decreases by lowering the multiplexity. Considering the features of real multiplexes, this is a rather reassuring result. In fact, on one side the spreading is favored by positive correlations. On the other, the effect of such correlations is far less important for small values of multiplexity, which are typical of real systems. Interestingly, the role of the multiplexity is drastically reduced in case the average connectivity in one layer is much larger than the other. In this scenario, which mimics the possible asymmetry in the contact patterns typical of different social contexts (e.g., family or work environment), one layer drives the contagion dynamics and the epidemic threshold is indistinguishable from the monoplex. However, the multiplexity is still significantly important in the other layer as the fraction of nodes present in both layers largely determines the spreading dynamics.

Some of these results are qualitatively in line with the literature of contagion processes unfolding on static and annealed multiplexes. However, as known in the case of single layered graphs, time-varying dynamics induces large quantitative differences [27, 28] . Indeed, the concurrency and order of connections are crucial ingredients for the spreading and neglecting them, in favor of static and annealed representations, generally results in smaller thresholds. While the limits of timescale separation might be relevant to describe certain types of processes, they might lead to large overestimation of the spreading potential of contagion phenomena.

The model presented here comes with several limitations. In fact, we considered the simplest version of the activitydriven framework in which, at each time step, links are created randomly. Future work could explore the role of more realistic connectivity patterns in which nodes activate more likely a subset of (strong) ties and/or nodes are part of communities of tightly linked individuals. Furthermore, we assumed that the activation process is Poissonian and the activity of each node is not a function of time. Future work could explore more realistic dynamics considering bursty activation and aging processes. All these features of real time-varying networks have been studied at length in the literature of single layered networks but their interplay with multiplexity when dynamical processes are concerned is still unexplored. Thus, the results presented here are a step towards the understanding of the temporal properties of multiplex networks and their impact on contagion processes unfolding on their fabric. Integrating over all activity spectrum of Eq. (4), it obtains the following equation:

Thus, Eq. (A1) can be further simplified as 

and

Integrating over all activity spectrum of Eq. (5), it obtains the following equation:

Multiplying both sides of Eq. (4) by a x , and integrating over all activity spectrum, we get the following equation:

Replacing x with A and B in Eq. (A6), respectively, we have

and

In the same way, multiplying both sides of Eq. (5) by a x and integrating over all activity spectrum, it obtains the following two equations:

when x is replaced with A and B, respectively. When the system enters the steady state, we have . Thus, the critical condition is determined by the following Jacobian matrix: 

If the largest eigenvalue of J is larger than zero, the epidemic will outbreak. Otherwise, the epidemic will die out. Specifically, if p = 0, two layers are independent, thus regulated by two independent Jacobians, and we can get the following two eigenvalues: For 0 < p < 1, we could not find a general analytical expression for the eigenvalues of J . However, the critical transmission rate can be determined by finding the value of λ leading the largest eigenvalue of J to zero. In other words, rather than solving explicitly the characteristic polynomial |J − I | = 0 and defining the condition for the spreading max > 0 as done above, we can determine the critical value of λ as the value corresponding to the largest eigenvalue to be zero [76, 77] . · · · · · · · · · · · ·

Multilayer Networks. Structure and Function

Multilayer networks

Mathematical Formulation of Multilayer Networks

The structure and dynamics of multilayer networks

Multiplex Networks: Basic Formalism and Structural Properties

Social Network Analysis: Methods and Applications

Multiplexity in adult friendships

Social capital in the creation of human capital

Measuring and modeling correlations in multiplex networks

Quantifying dynamical spillover in co-evolving multiplex networks

Effects of temporal correlations in social multiplex networks

The physics of spreading processes in multilayer networks

Spreading processes in multilayer networks

Diffusion Dynamics on Multiplex Networks

Global stability for epidemic models on multiplex networks

Epidemic spreading and bond percolation on multilayer networks

Optimal percolation on multiplex networks

Resource control of epidemic spreading through a multilayer network

Dynamical Interplay Between Awareness and Epidemic Spreading in Multiplex Networks

Epidemic spreading and risk perception in multiplex networks: a self-organized percolation method

Multiple routes transmitted epidemics on multiplex networks

Epidemics in partially overlapped multiplex networks

Immunization of epidemics in multiplex networks

Asymmetrically interacting spreading dynamics on complex layered networks

Conditions for Viral Influence Spreading Through Multiplex Correlated Social Networks

Contactbased social contagion in multiplex networks

Temporal networks

Modern temporal network theory: a colloquium

Activity driven modeling of time-varying networks

Time varying networks and the weakness of strong ties

Asymptotic theory of time-varying social networks with heterogeneous activity and tie allocation

From calls to communities: a model for time-varying social networks

The dynamical strength of social ties in information spreading

Persistence and periodicity in a dynamic proximity network

What's in a crowd? analysis of face-to-face behavioral networks

From seconds to months: an overview of multi-scale dynamics of mobile telephone calls

Persistence of social signatures in human communication

Fundamental structures of dynamic social networks

Face-to-face interactions

Random Walks and Search in Time Varying Networks

Quantifying the effect of temporal resolution on time-varying networks

Controlling Contagion Processes in Activity Driven Networks

Contagion dynamics in time-varying metapopulation networks

Epidemic spreading in time-varying community networks

Immunization strategies for epidemic processes in time-varying contact networks

Random walks on temporal networks

Analytical Computation of the Epidemic Threshold on Temporal Networks

Causality driven slow-down and speed-up of diffusion in non-markovian temporal networks

Spatio-temporal networks: reachability, centrality and robustness

Random walk centrality for temporal networks

Importance of individual events in temporal networks

Bursts of vertex activation and epidemics in evolving networks

Attractiveness and activity in internet communities

Contrasting effects of strong ties on SIR and SIS processes in temporal networks

Committed activists and the reshaping of status-quo social consensus

Betweenness Preference: quantifying Correlations in the Topological Dynamics of Temporal Networks

Bursty communication patterns facilitate spreading in a threshold-based epidemic dynamics

Birth and death of links control disease spreading in empirical contact networks

The basic reproduction number as a predictor for epidemic outbreaks in temporal networks

Statistical physics of vaccination

Social phenomena: From Data Analysis to Models

Burstiness and tie reinforcement in time-varying social networks

Random walks on activity-driven networks with attractiveness

Epidemic spreading in modular time-varying networks

Modeling Infectious Disease in Humans and Animals

Modeling dynamical processes in complex socio-technical systems

Epidemic processes in complex networks

The role of endogenous and exogenous mechanisms in the formation of r&d networks

Networks. An Introduction

Topological properties of a time-integrated activity-driven network

Structural measures for multiplex networks

Nature of the Epidemic Threshold for the Susceptible-Infected-Susceptible Dynamics in Networks

Telling tails explain the discrepancy in sexual partner reports

Concurrent partnerships and transmission dynamics in networks

Dynamical Processes on Complex Networks

A unified approach to percolation processes on multiplex networks

Clustering determines the dynamics of complex contagions in multiplex networks

Further, the maximum eigenvalue of matrix J M can be calculated asThus, the critical transmission rate isFurther, if the activities of the same node in each layer are the same, the above equation can be simplified as follows:(A20)