Temporal interactions facilitate endemicity in the susceptible-infected-susceptible epidemic model

Data of physical contacts and face-to-face communications suggest temporally varying networks as the media on which infections take place among humans and animals. Epidemic processes on temporal networks are complicated by complexity of both network structure and temporal dimensions. Theoretical approaches are much needed for identifying key factors that affect dynamics of epidemics. In particular, what factors make some temporal networks stronger media of infection than other temporal networks is under debate. We develop a theory to understand the susceptible-infected-susceptible epidemic model on arbitrary temporal networks, where each contact is used for a finite duration. We show that temporality of networks lessens the epidemic threshold such that infections persist more easily in temporal networks than in their static counterparts. We further show that the Lie commutator bracket of the adjacency matrices at different times is a key determinant of the epidemic threshold in temporal networks. The effect of temporality on the epidemic threshold, which depends on a data set, is approximately predicted by the magnitude of a commutator norm.


Introduction
A majority of infectious diseases, ranging from seasonal influenza to Ebola outbreaks and sexually transmitted infections, can be viewed to occur on contact networks of humans and animals, which are composed of individuals and dyadic links between them. Epidemic processes are one of the most widespread applications of network analysis. Structure of contact networks has been shown to affect, for example, the likelihood and speed of an infection penetrating into a significant part of a population, effectiveness of intervention strategies, and identification of super-spreading individuals [1][2][3].
Accumulating data evince that contact networks underlying epidemic and other processes are often highly dynamic, constituting temporal networks [4,5]. For example, links may be only occasionally used for actual physical contacts. Individuals may be socially active in some restricted periods of time. Temporality of networks may alter effects of networks on epidemic processes [6,7]. This is a practical enquiry because various instances of epidemics in human and animal populations, and also viral spreading of information in human society, seem to occur on temporally varying networks. However, our understanding of epidemic processes in temporal networks is still limited. Theory based on the branching process enables us to understand long-tail behaviour of the number of newly infected individuals [8][9][10]. Other theoretical approaches include analysis of epidemic spreading on theoretically tractable generative models of temporal networks [11][12][13][14][15][16][17][18][19][20][21]. However, these studies assume network models such that they may miss effects of the properties of temporal contact networks that are present in empirical data but not modelled. On empirical temporal networks, the most popular approach has probably been to run numerical simulations of epidemic processes on the empirical networks and their variants (e.g., [7,[22][23][24]).

Model
We consider the continuous-time SIS model on undirected temporal networks having N nodes, as schematically shown in figure 1. Each node assumes either the susceptible or infected state. An infected node infects each of its susceptible neighbours at rate β. An infected node transits to return to the susceptible state at rate γ. To model exogenous dynamics of networks, we consider an infinite sequence of adjacency matrices and sequentially apply each of them for time τ. In other words, A is applied between time t ¢ ℓ and t ¢ + ℓ ( ) 1 . We refer to each network applied for time τ as the snapshot network, or snapshot in short.
This switching network modelling of temporal networks is common in studying synchronisation processes [28][29][30]. It is also in accordance with observation of temporal network data at regular time intervals τ. Here we regard τ as a free parameter. It controls the relative time scale of the epidemic and network dynamics; multiplying τ by a constant > ( ) c 0 is equivalent to not changing τ and instead changing β and γ to b c and g c , respectively. In addition, changing (τ, β, γ) to ( t c , b c, g c) does not change the dynamics. Therefore, we set g = 1 without loss of generality.
The model simplifies over real epidemic processes. Heterogeneities are considered only through the quenched disorder of the network. A more refined model would account for heterogeneity also in the parameters, drawing β and γ separately for each node and τ for each snapshot from probability distributions. We suppress these refinements for the sake of tractability of the model.

Infections persist more easily in temporal than static networks
We denote the probability that node i (   i N 1 ) is in the infected state at time t by x i (t). Under the assumption that the states of different nodes are independent of each other, the individual-based approximation to the SIS dynamics linearised around the disease-free configuration (i.e., = ( ) x t 0 i for all i) is given by ,  denotes the transposition, and I is the identity matrix. When the network is where A is the adjacency matrix. The epidemic threshold, denoted by b c , is defined as the value of β above which infection can persist in the network. According to the individual-based approximation, b c for a static network is given by the value of β at which the leading eigenvalue of b -A I is equal to zero. The leading eigenvalue of b -A I is given by ba -1 max , where a max is the leading eigenvalue of A. Therefore, we obtain b a = 1 c m a x [3,41,42]. When the network varies in time, equation (1) yields Valdano and colleagues were the first to derive this result by analysing the SIS model in discrete time [25].
For two temporal networks, we simulated the SIS model dynamics using the quasistationary state (QS) method [43] tailored to the case of temporal networks (appendix A). We aggregated the original data over several time windows to define snapshot networks ¢ ℓ ( ) A ; the size of the aggregating window is shown in table 1. The average prevalence of infection for values of β and τ is shown in figure 2. Similarly to [25], we assumed a periodic boundary condition such that the first snapshot ensues after the last snapshot. The theoretical estimates of b c are shown by the arrows in figures 2(a) and (b), which are fairly accurate. Figure 2 suggests that the epidemic threshold, denoted by b c , decreases as τ increases in both networks. To formulate this point, we compare b c with the epidemic threshold for the SIS dynamics occurring on the aggregate, static network, denoted by b* c . We start with defining the aggregate network as the adjacency matrix . With this normalisation, the temporal and aggregate networks have the same time average of the weight of each link [30] (figure 1). The epidemic threshold for the aggregate network is given by b a = * * 1 c m a x , where a* max is the largest eigenvalue of * A . Any real-valued, continuous spectral function f that acts on the spectrum of its matrix argument and only attains finite values satisfies f f  [44]. Although a generalisation of this inequality to the case of more than two matrices is false in general, we conjecture that it remains true if the involved matrices have only non-negative entries and f is the spectral radius. We have the following evidence. First, to the best of our knowledge, all counterexamples involve matrices with negative entries [45]. Second, the theoretical results for two model temporal networks presented in section 4 are consistent with this inequality. Third, all numerical calculations performed in this paper on real data sets and synthetic networks are consistent with this inequality. By admitting this generalised inequality, and applying it to matrices bt Several remarks are in order. First, equation (4) implicitly assumes the periodic boundary condition for the temporal networks. However, because ℓ is arbitrary, equation (4) holds true for arbitrary sequences of networks. Second, in the limit t  0, equation (4) is satisfied with equality such that b b a = = * * 1 c c m a x . Third, if each snapshot is composed of a single link, numerical results suggest that the epidemic threshold increases as the temporality of the network increases (appendix B), which is opposite to the current result. In this situation, the probability that the infection is extinguished is not negligible even for a large infection rate. A theory that accounts for stochasticity, which is different from the deterministic individual-based approximation, correctly predicts the direction of the change in the epidemic threshold as τ increases (appendix B). It should be noted that a network mainly composed of isolated single links is realistic for sexually transmitted infections through monogamous relationships [11]. This situation is out of the scope of the following analysis. Our theory requires that each snapshot has a relatively large connected component such that the SIS dynamics on the snapshot are not significantly influenced by disconnected single links even if they are present. The stochasticity and the absorbing configuration (all nodes susceptible) become relevant not only for disconnected single links (appendix B) but also snapshots of increasing size as t  ¥. If τ exceeds the typical time of reaching the absorbing configuration on a snapshot, the epidemic ceases due to stochasticity.

Temporal networks with clique snapshots
To investigate the distance between the epidemic threshold for aggregate and temporal networks and its dependence on parameters, we start by calculating the epidemic threshold for two temporal network models. In the first model, snapshots consist of a disjoint union of cliques, each with + d 1 cl nodes, and isolated nodes. A clique may be a suitable model for conversation events in a small group [23,[46][47][48]. The clique size remains the same across different cliques and snapshots. The number of nodes stays constant across snapshots, but the number of cliques may depend on a snapshot.
We assume that snapshots are randomly and independently drawn from a set of possible snapshots with equal probability, which we call the random sampling with replacement. Whether the expected prevalence is positive or not can be determined by l t m = ¥ ℓ ( ) ℓ lim ln 1 max , which is equal to the maximum Lyapunov  [38] and (b) sexual contact data set [40]. See table 1 for details on the data sets. We imposed periodic boundary conditions. The theoretically obtained b c is shown by the arrows. We calculated the theoretical b c value by a bisection method on the leading eigenvalue of t ( ) T . We expanded t ( ) T up to the term of the order of t ( ) O 10 and calculated its leading eigenvalue using the power method. exponent associated with the switching linear dynamics induced by operator t ( ) T (equation (3)) as  ¥ ℓ [49][50][51]. The epidemic threshold corresponds to l = 0. It is algorithmically undecidable to determine whether l < 0 or not [52], hindering us from deriving the exact value of b c . Therefore, we evaluate the expected state vector, [ ( )] x t E , where E is the expectation. The expected state vector evolves according to Here, the summation runs over all possible snapshots, and r represents the number of possible snapshots. We denote the leading eigenvalue of t ( ) T by m ; max m t (ˆ) ln max approximates λ. It holds true that m t l (ˆ)  ln max (appendix C). In addition, we numerically verified b b »ĉ c for some networks and a range of parameter values, where the estimate of the epidemic threshold b c is obtained from m = 1 max (figure 3). As shown in appendix D, we obtain In figure 4, we test the accuracy of the theory against numerical simulations using a synthetic temporal network constructed as follows. In each snapshot, every node is independently activated with probability a i , which obeys a power-law distribution. Each activated node triggers a clique of size + d 1 cl by involving d cl other nodes drawn with equal probability. We allow multiedges in a snapshot. The degree of the aggregate network up to the leading order in terms of N is given by Figure 4(a) suggests that equation (7) (dotted line) is sufficiently close to the exact value of b c obtained through equation (6) (solid line). The small discrepancy between the exact and the approximated values are caused by the fact that cliques may overlap in the synthetic temporal networks, which the approximate theory does not assume. All these estimates accurately locate the position of the epidemic threshold obtained from direct numerical simulations ( figure 4(b)). Different temporal networks generated by the present model can have the same aggregate network. In equation (7), a* max is the same for all temporal networks sharing an aggregate network. Therefore, the epidemic threshold depends on the temporality of networks solely through d cl . Equation (7) implies that the epidemic threshold decreases as d cl increases for all values of t > 0. This observation implies that infection is more likely to be prevalent when snapshots are highly variable in the sense that some snapshots have many links and others have few links, as compared to when different snapshots have similar densities. , and r=100. The value of a i (e   a 1 i ) is distributed according to a power law with exponent equal to 3. We adjust ε such that the mean of a i equals 0.04.

Activity driven model
In the temporal networks composed of cliques (section 4.1), the degree of nodes is essentially homogeneous within each snapshot (i.e., d cl or 0). In this section, we study the case in which a snapshot consists of a disjoint union of stars, allowing snapshots to be heterogeneous in the node's degree. Each star is assumed to have one hub node connected to d hub leaves. A leaf is only adjacent to the corresponding hub. The value of d hub is assumed to be the same for different stars and snapshots. Different snapshots may contain different numbers of stars. As a special case of this model, we consider the discrete-time version of the activity driven model [17,21]. In each snapshot, every node i is activated with probability a i independently of other nodes. The variable a i plays a similar role to that in the case of clique snapshots but is distinct from it. For each activated node, d hub nodes are drawn with equal probability and connected to the activated node. Although stars in a single snapshot may overlap, we consider the case in which the overlap is rare.  (6) is shown by the solid line. The approximation, equation (7), is shown by the dotted line. The two lines severely overlap in the entire range of τ. (b) Fraction of infected nodes for the clique network model. The results for direct stochastic simulations are shown by the curves. The epidemic threshold predicted by equation (6) and equation (7) are shown by the solid and dashed arrows, respectively. (c) Epidemic threshold for the activity driven model. The exact value obtained from equation (6) is shown by the solid line. The approximation, equation (8), is shown by the dotted line. (d) Fraction of infected nodes for the activity driven model. The epidemic threshold predicted by equation (6) and equation (8) are shown by the solid and dashed arrows, respectively. We set N=2000, , r=1000, and let each a i obey the power-law distribution with the probability density function . We set h = 3 and adjusted ε to ensure á ñ = a 0.0025.
As shown in appendix E, the epidemic threshold is approximated as It should be noted that both equation (8) and a more exact estimate given by equation (39) (appendix E) converge to b a = * * 1 c m a x in the limit t  0 and to d 1 hub in the limit t  ¥. The high accuracy of equation (8) is confirmed in figure 4(c).
Equation (8) indicates that the epidemic threshold is small for a large value of d hub . This result is consistent with that for the network model with clique snapshots. In other words, if the aggregate network is the same, temporal networks that sometimes have dense snapshots and otherwise sparse snapshots would make the epidemic threshold small, as compared to temporal networks that have similar density of links across time.
According to the heterogeneous mean field approximation [17,21], the epidemic threshold for the activity driven model is equal to c (see appendix F for the derivation of a* max ). In their framework, network switching occurs sufficiently fast as compared to epidemic dynamics such that epidemic spreading is effectively occurring on the static, aggregate network. Their epidemic threshold is different from the well-known value for the configuration model [53] because the aggregate network of the activity driven network is different from the configuration model having the same degree sequence. In contrast, the present results capture how the time scale of the network dynamics as described by the activity driven model and that of the SIS dynamics interact.

Non-commuting snapshots lower the epidemic threshold
The amount of the shift in the epidemic threshold as we slow down the dynamics of the network (by increasing τ) depends on individual temporal networks. The online message network (figure 2(a)) experiences a larger shift than the sexual contact network ( figure 2(b)). In the model networks examined in section 4, the epidemic threshold is sensitive to the change in τ when d cl or d hub is large (equations (7) and (8)). In this section, we propose a quantity to predict the sensitivity of the epidemic threshold, b c , to temporality of networks, τ.
First of all, b c is independent of τ if any pair of adjacency matrices of the snapshots commutes. This is because the time evolution operator for the temporal network, t ( ) T (equation (3)), and that for the aggregate , coincide in this case. To quantify the difference between b c and b* c when matrices are non-commuting, we use Zassenhaus' formula [54] given by For instance, the first two matrices are given by By iteratively applying equation (10) with tb = s , we obtain On the basis of equation (4), this difference yields a difference between b c and b* c . Therefore, we define the degree of non-commutativity by The multiplicative constant a* 1 max in equation (16) normalises the leading eigenvalue of the aggregate network to unity. If C = 0, all pairs of the adjacency matrices of snapshots commute, and b c is independent of τ. If > C 0, at least some adjacency matrices do not commute.
We carry out numerical simulations to examine the relationship between C and the epidemic threshold. First, we generate r=50 commuting adjacency matrices of snapshots (appendix G), yielding C = 0. Then, we manipulate C by gradually swapping elements of adjacency matrices of different snapshots with the aggregate network fixed (appendix G). For each sequence of adjacency matrices thus obtained, we estimate the relative change of the epidemic threshold given by Figure 5 indicates that b D c increases roughly quadratically in C and suggests a high predictive power of C. Across several temporal network data sets summarised in table 1, the dependency of b D c on C is shown in figure 6(a). The figure also contains results for temporal networks with clique snapshots and the activity driven model. These networks were generated with different values of d cl and d hub under the condition that the aggregate network was approximately the same within the same model (appendix H). Consistently with figure 5, C is a strong determinant of the decrease in the epidemic threshold across various temporal networks. It should be noted that C and b D c are strongly correlated despite different sizes of the empirical networks.
Various types of temporal correlation in empirical data are known to affect epidemic processes [7,[22][23][24]. Therefore, the epidemic threshold may be influenced by the order of snapshots, whereas C is not. To examine this point, we calculated the epidemic threshold for each network using the order of snapshots given by random sampling with replacement. The results are shown in figure 6(b). For these decorrelated temporal networks as well, C and b D c are strongly correlated. Comparison between figures 6(a) and (b) reveals that the loss of temporal correlation somewhat decreases b D c for all data sets. However, C is clearly a stronger determinant of b D c than the order of snapshots is. Figure 5. Relationship between the epidemic threshold and the degree of non-commutativity for the synthetic temporal networks. Each point corresponds to a different sequence of snapshots manipulated from the original one such that the aggregate network is the same. The original sequence was generated with r=50 snapshots and N=10 nodes each. We used a small N value because it was computationally costly to generate snapshots with larger N. We set t = 1.

Discussion
We have provided evidence that the epidemic threshold for the SIS model on temporal networks is smaller than that for the corresponding static networks for arbitrary temporal networks. We have also shown that the degree of commutativity of the adjacency matrices of the snapshot networks predicts the impact of temporality of the network on the epidemic threshold. Our results are opposite to the previous results concluding that infection in the SIS model is less likely in temporal than static networks [11,17,20]. However, our results do not contradict theirs, for which the aggregate network obtained from the temporal network is not equal to the static network used for comparison. In contrast, we compared temporal networks with static networks such that they are the same if we ignore the temporal information in the former.
In the discrete-time SIS model on temporal networks [25], the epidemic threshold is larger for the temporal network than the corresponding aggregate network when a sequence of snapshots is randomly drawn with replacement (appendix I). This result is opposite to the current results for the continuous-time dynamics. Because a discrete-time model is a proxy to the continuous-time counterpart, which is usually more realistic, the present continuous-time framework is useful. For instance, the discrete-time SIS model implicitly assumes that, for time τ, each node is allowed to make at most one transition between the infected and susceptible states. Therefore, the propagation speed is restricted, which is not the case in the continuous-time framework. In particular, the discrete and continuous-time versions only coincide when t  0, in which limit both reduce to the SIS model on the aggregate network.
We modelled temporal networks by switching networks. In practice, we cannot manipulate the duration of each snapshot (i.e., τ) because it is specified by the data, reflecting the temporal resolution of the observation. Rather, our interpretation of τ is the relative time scale between the epidemic dynamics and network dynamics. A large τ, with which the epidemic threshold decreases, corresponds to fast epidemic dynamics relative to network dynamics. The present results imply that infection is pronounced when the network varies over time (i.e., temporal network) but only slowly. In this situation, there are times when some snapshots strongly favour infection as compared to typical snapshots, and such snapshots enhance infection more than other snapshots suppress it. It should be noted that the time scale of the epidemic dynamics does not affect the equilibrium of the SIS model in the case of static networks. A future extension of the model may relax the assumption of equal duration τ of each snapshot and explore the effect of a (possibly broad) distribution of time spans. In a further step, a dependence of τ on the nodes' state may be considered to model contact avoidance and quarantine of infected nodes.
Our commutativity result opens the way to contain infection by devising the set of snapshots without changing τ or the structure of the aggregate network. In a hospital, it may be undesirable to change aggregate interactions between doctors, nurses, and patients because the amount of the interactions may be positively correlated with service quality. We can increase the epidemic threshold (therefore, less epidemic) by designing a sequence of interactions such that the corresponding adjacency matrices commute as much as possible. The method explained in appendix G is useful in systematically generating commuting adjacency matrices. For a similar attempt in synchronisation dynamics using single-link snapshot networks, see [55]. The adjacency matrices obviously commute when the snapshots do not have any nodes in common. Therefore, designing interactions such that different sets of nodes are active at different times as much as possible may be effective at increasing the epidemic threshold. This point warrants further work. Examining the influence of τ on the epidemic threshold and the relevance of the commutativity of the adjacency matrices in more complex compartmental models of epidemic dynamics also warrants future work.

Appendix A. QS method
The QS method is used for computing the average prevalence in the SIS model in finite populations [43]. The QS method for the SIS model in static networks works as follows. We distinguish active states, which are configurations with at least one infected node, from the absorbing state, which is the configuration with no infected nodes. We keep a total of 2000 active states in memory. After one update event in the SIS model, we are either in an active or an absorbing state. If an active state is reached, it replaces a randomly chosen state in the memory with the probability proportional to the expected time to the next update, where the proportionality constant is set to 0.5. It should be noted that the probability does not exceed unity because the mean time to recovery is normalised to unity such that the expected time to the next update in the entire network is less than unity. If the absorbing state is reached, one active state is chosen uniformly at random from the memory to replace the absorbing state. After a transient of 10 3 time units, the QS is calculated as the average of the system's state over the next 10 3 time units. We modified the QS method for temporal networks. Our implementation is slightly different from that in [25]. Assume that we are currently using the nth snapshot. In other words, . The QS may depend on t --( ) t n 1 , i.e., the time since the beginning of the current snapshot. Therefore, for different values of t --( ) t n 1 and the current snapshot, we generate a memory, i.e., a list of active states to be used when the process dies. To this end, we divide the time window t [ ) 0, into those of length t ¢ . If t t < ¢ , we set t t ¢ = and the time window is not divided. For every combination of the currently used snapshot and discrete time t ¢ ¢ n , which covers t t t ¢ ¢ -< ¢ + ¢ ( )  n t n n 1 ( ¢ = ¼ n 0, 1, ), we allocate memory to store 2000 active states. When the process is in an active state at time t t -= ¢ ¢ t n n , it replaces a randomly chosen active state in the memory corresponding to t ¢ ¢ n and the current snapshot with probability 0.5. If the process reaches the absorbing state at time , the network state is replaced by an active state randomly chosen from the memory corresponding to the elapsed time t ¢ ¢ n and the current snapshot. Then, the process is restarted at time t t = + ¢ ¢ t n n . After a transient of time t r 10 3 , where r is the number of snapshots, we calculate the steady state as the fraction of the infected nodes averaged over the r 10 3 subsequent measurements conducted when snapshots switch from one to another. We used t ¢ = 0.5 and confirmed that the results remained unaffected with t ¢ = 0.1.

Appendix B. Stochastic SIS model on single-link snapshots
Consider a temporal network in which each snapshot contains disjoint single links drawn at random. The degree of every node in a snapshot is at most one. Numerical results on a temporal network generated from this model with N=2000 nodes and 500 links in each snapshot are shown in figure I1. Contrary to the results shown in the main text, the epidemic threshold increases as τ increases. This result is caused by the fact that the individual-based approximation does not work, even qualitatively, when snapshots are composed of small fragments such as disjoint single links. For such networks, stochasticity of the dynamics that the individual-based approximation does not account for plays a significant role. In short, even if the infection rate is very large, infection will die out if we apply a snapshot for long τ. To understand this situation, here we analyse the stochastic SIS model on a single link = ( ) N 2 by the master equation rather than by the individual-based approximation.
The SIS model with two distinguishable nodes has = 2 4 2 states, where each node is either susceptible or infected. We assume that the two nodes are bi-directionally coupled. We denote the time-dependent probabilities of the states by p SS (none infected), p IS (only node 1 infected), p SI (only node 2 infected), and p II (both nodes infected). It should be noted that . The probabilities evolve as To exploit the symmetry, we consider = + u p p IS SI and = r p p IS SI in place of p IS and p SI . They evolve as b = -- The equations for r, u, and p II fully describe the dynamics. The solution of equation (22), is decoupled from the dynamics of u and p II . The linear dynamics composed of equations (20) and (21) have the eigenvalues given by The corresponding left eigenvectors are given by b k b - (( ) ( ) ) 1 2 1. For initial conditions u(0) and ( ) p 0 II , the solution reads In what follows, we assume that one node is initially infected and the other is initially susceptible, yielding = ( ) u 0 1, Î -( ) { } r 0 1,1 , and = ( ) p 0 0 II . Then, the solution at time τ is given by r 0 e . 29 1 We denote the probability for nodes 1 and 2 to be infected at time τ by t ( ) p 1 and t ( ) p 2 , respectively. They are given by  The expected number of infected nodes, t t + ( ) ( ) p p 1 2 , when one node is initially infected is shown in figure I2(a). The individual-based approximation developed in the main text would predict that the number of infected nodes monotonically increases in time when β is sufficiently large. However, figure I2(a) indicates that the infection eventually dies out even for a large β value. Next, we estimated the epidemic threshold as the β value at which t t . This estimate gives a τ-dependent lower bound on the epidemic threshold in temporal networks, b c . The epidemic can only spread if the expected number of nodes after time τ is larger than the number at time 0. The relationship between the estimated epidemic threshold and τ is shown in figure I2(b). We find that the epidemic threshold increases as τ increases.

Appendix C. Relationship between the Lyapunov exponent and m max
We assume random sampling with replacement of snapshots. The maximum Lyapunov exponent is given by where ℓ is the length of a sequence of snapshots. It should be noted that, while m max is a random value for any ℓ, the maximum Lyapunov exponent is a deterministic value owing to theorem 2 in [49]. In this section, we show that m t l (ˆ)  ln max . We use theorem 1 in [51]. Suppose that any adjacency matrix of the snapshot is a matrix with non-negative entries and that, for each pair   (6) and (40), we obtain Therefore, b c converges to b a = * * 1 c m a x in the limit t  0. In the limit t  ¥, b c converges to d 1 cl , representing the fact that just one snapshot is used indefinitely long, and the epidemic threshold in a single snapshot is that of a clique.
Finally, we approximate a á ñ » * * d u , cl max for general τ values in equation (44) to obtain equation (7). It should be noted that equation (44) is independent of á ñ * d u , cl in the limit t  0 and t  ¥ such that equation (7) is exact in these two limits.   If more than one star is embedded in a snapshot, appropriately permutated versions of equation (32) are added together. By applying equation (47) to equation (6), we obtain where the summation of A 2 runs over all possible snapshots and r is the number of the possible snapshots. When = d 1 hub , we obtain å = * A r D 2 such that equation (48) is consistent with equation (41). As a special case, we consider the discrete-time version of the activity driven model [17,21]. In each snapshot, every node i is assumed to be activated with probability a i independently of the other nodes.
If node i is a hub, the probability that it connects to a node j in a snapshot is equal to d N hub . Therefore, we obtain up to the order of N 1 where we neglected the probability that both i and j are hubs and are connected to each other. If node k is a hub, the probability that it selects both nodes i and j (¹i) as leaves is equal to is the eigenvector corresponding to eigenvalue N. Because all entries of this matrix are positive, small changes in the eigenvalues will result in a matrix with positive entries. We choose eigenvalues l Î - ), e > 0 uniformly at random and calculate a new adjacency matrix as We set e = 10 in the numerical simulations. We repeat this procedure until we obtain r matrices with positive entries. The r adjacency matrices commute within themselves.
We manipulate the r commuting matrices to increase the degree of commutativity, C (equation (16)), while conserving the aggregate matrix as follows. First, we select two matrices. Then, we randomly select (i, j),   i j N 1 , and swap the (i, j) entry of the two matrices. To keep both matrices symmetric, we also swap the ( j, i) entry of the two matrices. We repeat this procedure 2000 times during which C tends to increase.
Appendix H. Generation of temporal networks with different d cl and d hub values and the same aggregate network We generate temporal networks with clique snapshots and the activity driven model with different values of d cl and d hub and the common aggregate network as follows.
For temporal networks with clique snapshots with given d cl and a, the aggregate matrix up to the order of N 1 is given by In figure 6, we chose d cl , = d 10 hub and drew a i from a power-law distribution with exponent 3 and mean á ñ = a 0.05. Figure I1. Results of stochastic simulations on a temporal network composed of disjoint links. We set N=2000 and r=1000. Each snapshot contains 500 disjoint links.
In this section, we show b b*  c,disc c , where b c,disc is the epidemic threshold for the discrete-time SIS model in which snapshots are randomly sampled with replacement from a given set.
For the SIS model in discrete time, the time evolution operator is given by equation (8) in [25] as follows: