Epidemic spreading in modular time-varying networks

We investigate the effects of modular and temporal connectivity patterns on epidemic spreading. To this end, we introduce and analytically characterise a model of time-varying networks with tunable modularity. Within this framework, we study the epidemic size of Susceptible-Infected-Recovered, SIR, models and the epidemic threshold of Susceptible-Infected-Susceptible, SIS, models. Interestingly, we find that while the presence of tightly connected clusters inhibits SIR processes, it speeds up SIS phenomena. In this case, we observe that modular structures induce a reduction of the threshold with respect to time-varying networks without communities. We confirm the theoretical results by means of extensive numerical simulations both on synthetic graphs as well as on a real modular and temporal network.


The Model
The network is defined by means of the following parameters: -the total number N of nodes in the network; -the activity distribution parameters, i.e. the lower cut-off ε and the leading exponent ν, so that F(a) ∝ a −ν for a ∈ [ε, 1]; -the lower cut-off s min , the upper limit s max and the exponent ω describing the community size distribution, i.e. P(s) ∝ s −ω for s ∈ [s min , s max ]; µ is the probability that, once active, a node will connect to a node inside the community, so that µ = 1 − µ is the probability to fire outside the community; We initialize the network extracting N activity values from the activity distribution F(a) and we then group the nodes in communities of size s drawn from a size-distribution P(s). Once we initialized the network we let it evolve following the time-varying activity driven framework. At each time step t we start with N disconnected nodes. Each node gets active with probability a i dt at each time step dt and fire to a randomly chosen node inside (outside) its own community with probability µ (µ ). At time t + dt we delete all the edges and repeat the above procedure. Each node i will then have a set of neighbors that have been contacted or have contacted the node during the network growth. The size of such a set is the integrated degree k of the node i. Of these k neighbors, k c will be inside the i's community (in-community degree) and k o = k − k c will be external to the community (out-community degree).

The Network growth
In the integrated network, each node i has a set of neighbors that have been contacted or have contacted the node during the network growth. The size of such a set is the integrated degree k. Of these k neighbors, k c are inside the i's community (i.e. the in-community degree) and k o = k − k c are external to the community (i.e. the out-community degree). Since the model is memoryless, the in-community degree k c and the out-community degree k o are decoupled and can, in fact, be treated separately. Even the activity potential a i of each node can be "split" in two components: the in-community activity µa i = (1 − µ )a i and the complementary out-community µ a i . Indeed, each node points, on average, a fraction µ of its own events toward the community, while the remaining µ are directed outside the community itself. Each node experiences a mean field of activity, µ a , coming from the community (provided that the community is large enough) and a supplementary external field µ a coming from the rest of the network.

The network time scales
As a first insight, let us note that the in-degree time dependence can be easily approximated with a probabilistic consideration. Each node i of activity a i , within a community of size s, has s − 1 available edges. Now, for each time step of the dynamics, the edge e i j is created with probability µ(a i + a j )/(s − 1). Then, on average, each edge emanating from i is activated with probability c(a i , µ)/(s − 1) = µ(a i + a )/(s − 1), where c(a i , µ) is the number of edges intra-community. The probability P(a i , µ, s,t) for an edge pointing to i not to be activated after t time steps then reads: Since the in-degree k c ∈ [0, s − 1], we can write: In the following, we always have the dependency on a i , µ and t, so to simplify the notation we drop most of those parameters: P(a i , µ, s,t) = P (s), to avoid confusions with the community size distribution P(s). Also, k c (a i , µ, s) = k c (s) and c(a i , µ) = c. We note that Eq. 1 gives us an estimation of the characteristic time τ(s) that takes for a node of activity a i to saturate the in-degree k c → (s − 1). Indeed, we can rewrite Eq. 1 as: ( So, as expected, the saturation time (i.e. the typical time for k c to be of the same order of s) increases as the activity µa i decreases and/or the community size s grows. Generalizing the above reasoning, the characteristic time for a community to have the majority of the nodes saturated is obtained by evaluating the probability P e (s) to create (on average) an edge e i j in a community of size s in a single evolution step. In other words, it is the number of edges activated in one step divided by the total number of possible edges in the network: The probability for one edge not to be created after t time steps is then: where τ c (s) represents the typical time by which the majority of the nodes of a community has a degree k c s. Note that, in the evaluation of both τ(s) and τ c (s), we did not take into account the difference between edges pointing to a more active node and the ones pointing to a less active one. Nevertheless, this is a simple estimation that, as we will show later, correctly catches the general behaviour of the in-degree k c for any value of µ, a i and s. Besides, when computing the key features of the evolving network, we are now able to distinguish the short time range t τ c (s) (in which k c s for any activity value a i ) and the long time limit t τ c (s) (in which k c ∼ s for any activity value a i ).

The Master Equation and the P(a, k,t)
We can now write down the Master Equation (ME) for the quantities P c (s, k c ) and P o (s, k o ), that is, the probability for a node of activity a i belonging to a community of size s to have degree in (out) degree k c (k o ) at time t. In general, a i ∆t represents the probability the node i is active, where a i is the activity rate of node i. Without loss of generality we will assume ∆t = 1. To get the ME for the in-degree k c distribution, we exploit also the time-dependence for a couple of passages: where ∑ j∼i and ∑ j i are respectively two contracted notations for the sum of all the first neighbors of node i and the sum of all nodes but the neighbors of node i. The first parenthesis indicates the probability that none of the nodes in the network fire. Third (sixth) term is the probability a node i, in the instantaneous network, is active and fires to a node where, in the integrated counterpart, there is already (isn't) a link. Four (seventh) term is the probability a node j not linked to i fires to any of the nodes but i (fires to i). Fifth factor is the probability a node j already linked to i fires. After some algebra, ME can be written as: Now we pass to the continuum limit by considering t 1 and k 1. So the l.h.s becomes simply the time derivative with respect to P c (s, k c ) and, to obtain a proper convergence of the results, we can expand the probability with respect to the incommunity degree up to second order. In the regime t τ(s), we can neglect k c s and 1/s ∑ j i a j ≈ a .
where we dropped the a i index since we expect all the nodes of a given activity to behave in the same way. Now a is an activation rate and in the treatment we assume it takes small values to avoid that two nodes become active together. The solution of Eq. 8 reads: where C is a normalization constant. By following the same procedure, we recover the same results of Eq. 9 for the out-community degree k o (a,t), by substituting µ → µ and k c → k o : Since N s max , the out-degree k o N for any time t of the process, thus we assume that Eq. 10 is valid for all the time scales analyzed. Also note that, as expected, the net effect of the mixing parameter µ is just a time rescaling of the out-community and in-community activity, respectively. Then, in the t ∼ τ c (s) time range is not possible to find an analytic formula, however simulations will be run to provide, at least, a qualitative behavior. In the t τ c (s) time limit the P c (s, k c ) converges to the δ (k c − (s − 1)) distribution. In fact, all the nodes will have all their edges activated and the P c (s, k c ) time derivative goes to zero.
Let us now resume the results found in this section: where we now distinguish between the in-community degree distribution P c (s, k c ) and the out-community degree distribution P o (s, k o ). The latter however, is independent on the community size and we can then define P o (s, k o ) = P o (k o ). So far we treated the two probability functions separately when, in fact, k c and k o are bound by the relation k = k c + k o . The total degree distribution P(s, k) will then be determined by the convolution of both the P c (s, k c ) and P o (s, k − k c ): 3/9 where we integrate over all the possible arrangements of the k c edges.
In the t τ(s) limit, by substituting Eq. 11a and 12 in Eq. 13, we sum the two exponents getting: where C is, again, a normalization constant. By combining the two terms and after some algebra we get: The integration over k c gives: where Erf(x) is the error function evaluated at x. In the small time limit, if we want to evaluate the P(k) = s max s min dsP(s)P(s, k), we have to consider t min s (τ c (s)). In this way, for each community we can use Eq. 16 as the true value of the P(s, k). The integration over the different community size s is then straightforward since the terms are independent on it, giving P(k) = P(s, k). Note that this result holds for any value of µ, N and a. The computation of P(k) in the large time limit (i.e. t τ c (s) ) is more complicated and we have to assume that k c = s − 1 for each node in a community of size s, otherwise Eq. 11b put everything equal to zero. The P o (s, k o ) will still be approximated by Eq. 12. The integral now reads: The exponential can be written as: where the rise of new terms proportional to s 2 and sk makes it difficult to perform the integral. We can, however, give a solution for the simple case P(s) = δ (s −s), when all the communities have equal size. First of all, we have, for large times, P c (s, k c ) = δ (k c − (s − 1)). Then: When k < s − 1, for sure k c < s − 1 and the delta put everything equal to zero. In the other case k c = s − 1 with a certain probability P o (k − (s − 1)). Finally, equation 17 can be written as:

The average degree k(a, s,t)
We can provide a simple expression for the nodes average degree belonging in different classes. As we already showed in Eq. 2, k c (a, s,t) grows as: The average total degree k(a, s,t) for nodes of activity a belonging to communities of size s depends on the time scale we analyse the problem. For small times (i.e. t τ(s) s ): As time grows toward the regime of times comparable to the average (i.e. t ∼ τ(a, s,t) s ), we cannot approximate the in-community degree anymore but we use directly equation 2: Then, the regime of large times (i.e. t τ(s) s ) is: The above equations then predict that the average degree has a linear growth proportional to a + a for short time limit (equation 23), then a transition for t ∼ τ(s) s is followed by a second linear growth, valid for large times, proportional to µ (a + a ). These regimes correspond to the initial growth, in which both the in-community and the out-community degrees are growing linearly in time, followed by the slowing down of the in-community degree which is saturating to s − 1. Finally, the third regime is again linear and it is driven by the µ (a + a ) coefficient: it means meaning that only the out-community degree is growing.

The degree distribution ρ(k th )
Now that we have the expression of the average degree, it is straightforward to write the degree distribution. At all the time scales we found k th ∝ Ct, where C is a time-independent coefficient. Then, k th ∝ at, and it results in an equal increment in the activity and in the degree values (da = dk th ). If we use the change of variable rule, we obtain: i.e., the degree distribution has the same exponent ν as the activity distribution function.

Comparison with numerical simulations
To check the analytical predictions of Section 2 we performed numerical simulation. In particular we realized 100 representations of a network featuring: -N = 10 5 nodes with modularity µ = 0.9 evolving for 10 5 evolution steps; -activity potential distributed following the F(a) ∝ a −ν with ν = 2.1 and a ∈ [10 −3 , 1] interval; -power-law distributed community sizes P(s) ∝ s −ω with ω = 2.1 and s ∈ [10, In order to analyze the collective behavior of the nodes we group them by their activity and community size, thus defining b classes of nodes. We average over the representations of the network and for each class of nodes b we evaluate: -P c (s, k), P o (k), P(s, k) for t τ c (s) and t τ c (s); -the average degree k c (s) , k o and k(s) ; -the degree distribution ρ(k th ).
In the main discussion, we already showed that some of the above measures have a great agreement with analytical predictions. To complete the discussion, we add below Fig. 1 which proves that also the in-community probability and the in-community degree perfectly matches our expectation. In panel B) we display P c (s, k) for t ∼ τ c even if it was impossible to obtain an exact result, the comparison demonstrates that the in-community probability starts to deviate from a Gaussian distribution. Figure 1. Panels A-C) rescaled P c (k c ) probability distribution as found for a selected node class at different times (legends). Functions in panels A-B) are rescaled accordingly to the theoretical distribution given in Eq. 11a, i.e., by sending k c →k c = (k c − k c )/ k c 1/2 and plotting P c (k c ) k c 1/2 . Panel C) is rescaled accordingly to the theoretical distribution given in Eq. 11b, i.e., by sending k c → k c = k c / k c and plotting P c (k c ) k c . In panel D) we plot k c for nodes featuring different activity potential and belonging to communities of different size. The data are rescaled sending the time t →= t/τ and then plotting k c (t/τ) /s. The analytical prediction of Eq. 2 computed for a node of activity a = a and belonging to a community of size s = s . Each point is an average over 10 2 independent simulations, parameters used are: N = 10 5 , ω = 2.1, ν = 2.1, m = 1, s min = 10 and µ = 0.9 6/9

SIR and SIS processes on modular activity driven networks
We present together all the results about SIR and SIS models obtained in synthetic networks. We start our simulations by setting 1% of randomly selected nodes as initial infected seeds, the other parameters are fixed as γ = 0.01, m = 3, ν = 2.1, N = 10 5 and Y = 0.5. Panels A) and C) set the exponent of the distribution of community sizes ω = 1.5, while panels B) and D) have ω = 2.1. The qualitative picture is unchanged due to selecting a different value of ω. The modular structure becomes irrelevant for µ → 0, whilst it significantly modifies the spread of the disease when µ 0. Modularity slows down contagion processes in SIR models, it favors the disease outbreak in SIS models. Moreover, quantitatively, in SIR models the presence of larger communities lead to higher epidemic size, in SIS models it lower the epidemic threshold. Moreover, in SIS models and when µ → 1, the epidemic threshold is likely to increase due to major limitations in reaching an endemic size. This phenomenon is particularly visible in Panel C), solid blue curve. Figure 2. Panels A-B) take R max of each R ∞ curve and plot it as a function of µ. Panels C-D) ξ SIS , that is β /γ in correspondence of L max , and plot it as a function of the modularity. In red curves we set s min = 100, in blue curves s min = 10. Solid curves are obtained by drawing community sizes from a power law distribution, 95% confidence interval is in gray. Dashed curves have a constant community size equal to the average value of the power law. Each point is an average of 10 2 independent simulations.

7/9 5 Real network: APS dataset
In the data each author of an article is described as a node. An undirected link between two different authors is drawn if they collaborated in the same article. We used the dataset from Ref 1 which spans a period between 1893 and 2006. To have the average degree in each instantaneous network as comparable as possible (see Fig. 3), we select a period of ten years, from January 1997 to December 2006. In this time window we register 96940 scholars who create 692667 connections. When we simulate SIR and SIS models on top of APS temporal network, we use periodic boundary conditions to let the disease dynamics evolve without late-time constrains.

Figure 3.
Average degree for each month in the selected subset of APS dataset (red circles), compared with the global average of all the ten years considered (solid black line). We label each month with increasing integer numbers from 1 to 120, where 1 represents the beginning of our sample, January 1997, and 120 the end, December 2006. We observe an increasing number of collaborations through years.
Using OSLOM, in Fig.4A we show that the integrated APS network is modular. Then, we apply the following degreepreserving randomization technique to destroy the network's community structure. We choose randomly a source node S 1 and, among its neighbors, we select randomly a target node T 1 . We do the same for other two nodes S 2 and T 2 . If the two pairs are equivalent or if a multi-edge will be created, we start back by selecting S 1 again and repeating the instructions. Otherwise, we swap T 1 and T 2 to have the new undirected links: S 1 − T 2 and S 2 − T 1 . The edges are chosen within the same temporal network and the number of swaps is equal to the number of edges in that instantaneous network. So, the above procedure is applied for each instantaneous network. At the end, we integrate the randomized temporal network and use OSLOM to detect the community structure. In Fig.4B, we qualitatively prove that our degree preserving randomization destroys the network's community structure. The options used in OSLOM are: -uw (to study undirected networks); -cp0.99 (to have communities as large as possible); -hr (to avoid to consider hierarchies); -r100 (to repeat 100 times the community detection). Since OSLOM finds communities in a non-deterministic way, last option is useful to get rid of stochastic fluctuations and have a more reliable community structure. Finally, we evaluate quantitatively the modularity of the two networks. For the original APS network, Q = 0.6685, and for its randomized counterpart Q = 0.0937.