Concurrency-induced transitions in epidemic dynamics on temporal networks

Social contact networks underlying epidemic processes in humans and animals are highly dynamic. The spreading of infections on such temporal networks can differ dramatically from spreading on static networks. We theoretically investigate the effects of concurrency, the number of neighbors that a node has at a given time point, on the epidemic threshold in the stochastic susceptible-infected-susceptible dynamics on temporal network models. We show that network dynamics can suppress epidemics (i.e., yield a higher epidemic threshold) when the nodes' concurrency is low, but can also enhance epidemics when the concurrency is high. We analytically determine different phases of this concurrency-induced transition, and confirm our results with numerical simulations.

Introduction: Social contact networks-on which infectious diseases occur in humans and animals or viral information spreads online and offline-are mostly dynamic.Switching of partners and the (usually non-Markovian) activity of individuals, for example, shape network dynamics on such temporal networks [1][2][3].A better understanding of epidemic dynamics on temporal networks is needed to help improve predictions of, and interventions in, emergent infectious diseases, to design vaccination strategies, and to identify viral marketing opportunities.This is particularly so because what we know about epidemic processes on static networks [4][5][6][7] is only valid when the time scales of the network dynamics and of the infectious processes are well separated.In fact, temporal properties of networks, such as long-tailed distributions of intercontact times, temporal and cross-edge correlation in intercontact times, and entries and exits of nodes, considerably alter how infections spread in a network [1-3, 8, 9].
In the present study, we focus on a relatively neglected component of temporal networks, i.e., the number of concurrent contacts that a node has.Even if two temporal networks are the same when aggregated over a time horizon, they may be different as temporal networks due to different levels of concurrency.Concurrency is a longstanding concept in epidemiology, in particular in the context of monogamy or polygamy affecting sexually transmitted infections [10][11][12].Modeling studies to date largely agree that a level of high concurrency (e.g., polygamy as opposed to monogamy) enhances epidemic spreading in a population.However, this finding, while intuitive, lacks theoretical underpinning.First, some models assume that the mean degree, or equivalently the average contact rate, of nodes increases as the concurrency increases [13][14][15][16].In these cases, the observed enhancement in epidemic spreading is an obvious outcome of a higher density of edges rather than a high concurrency.Second, other models that vary the level of concurrency while preserving the mean degree are numerical [10,11,17,18].In the present study, we use the analytically tractable activity-driven model of temporal networks [19][20][21][22][23] to explicitly modulate the size of the concurrently active network with the structure of the aggregate network fixed.With this machinery, we carefully treat extinction effects, derive an analytically tractable matrix equation using a probability generating function for dynamical networks, and reveal non-monotonic effects of link concurrency on spreading dynamics.We show that the dynamics of networks can either enhance or suppress infection, depending on the amount of concurrency that individual nodes have.Note that analysis of epidemic processes driven by discrete pairwise contact events, which is a popular approach [1-3, 9, 23-27], does not address the problem of concurrency because we must be able to control the number of simultaneously active links possessed by a node in order to examine the role of concurrency without confounding with other aspects.

Model:
We consider the following continuous-time susceptible-infected-susceptible (SIS) model on a discrete-time variant of activity-driven networks, which is a generative model of temporal networks [19][20][21][22][23].The number of nodes is denoted by N .Each node i (1 ≤ i ≤ N ) is assigned an activity potential a i , drawn from a probability density F (a) (0 < a ≤ 1).Activity potential a i is the probability with which node i is activated in a window of constant duration τ .If activated, node i creates m undirected links each of which connects to a randomly selected node (Fig. 1).If two nodes are activated and send edges to each other, we only create one edge between them.However, for large N and relatively small a i , such events rarely occur.After a fixed time τ , all edges are discarded.Then, in the next time window, each node is again activated with probability a i , independently of the activity in the previous time window, and connects to randomly selected nodes by m undirected links.We repeat this procedure.Therefore, the network changes from one time window to another and is an example of a switching network [28][29][30][31].A large τ implies that network dynamics are slow compared to epidemic dynamics.In the limit of τ → 0, the network blinks infinitesimally fast, enabling the dynamical process to be approximated on a time-averaged static network, as in [30].
For the SIS dynamics, each node takes either the susceptible or infected state.At any time, each susceptible node contracts infection at rate β per infected neighboring node.Each infected node recovers at rate µ irrespectively of the neighbors' states.Changing τ to cτ (c > 0) is equivalent to changing β and µ to β/c and µ/c, respectively, while leaving τ unchanged.Therefore, we set µ = 1 without loss of generality.Analysis: We calculate the epidemic threshold as follows.First, we formulate SIS dynamics near the epidemic threshold on a static star graph, which is the building block of the activity-driven model, while explicitly considering extinction effects.Second, we convert the obtained set of linear difference equations into a tractable mathematical form with the use of a probability generating function of an activity distribution.Third, the epidemic threshold is obtained from an implicit function.For the sake of the analysis, we assume that star graphs generated by an activated node, which we call the hub, are disjoint from each other.Because a star graph with hub node i overlaps with another star graph with probability ≈ m j =i a j (m + 1)/N ∝ m 2 a , where a ≡ daF (a)a is the mean activity potential, we impose m 2 a ≪ 1. (However, our method works better than the so-called individual-based approximation even when m 2 a = 0.5, as shown in the Supplemental Material.)We denote by ρ(a, t) the probability that a node with activity a is infected at time t.The fraction of infected nodes in the entire network at time t is given by ρ(t) ≡ daF (a)ρ(a, t).Let c 1 be the probability with which the hub in an isolated star graph is infected at time t+τ , when the hub is the only infected node at time t and the network has switched to a new configuration right at time t.Let c 2 be the probability with which the hub is infected at t + τ when only a single leaf node is infected at t.The probability that a hub with activity potential a is infected after the duration τ of the star graph, denoted by ρ 1 , is given by In deriving Eq. ( 1), we considered the situation near the epidemic threshold such that at most one node is infected in the star graph at time t [and hence ρ(a, t), ρ(t) ≪ 1].
The probability that a leaf with activity potential a that has a hub neighbor with activity potential a ′ is infected after time τ is analogously given by where c 3 , c 4 , and c 5 are the probabilities with which a leaf node with activity potential a is infected after duration τ when only that leaf node, the hub, and a different leaf node is infected at time t, respectively.We derive formulas for c i (1 ≤ i ≤ 5) in the Supplemental Material.The probability that an isolated node with activity potential a is infected after time τ is given by e −τ ρ(a, t).By combining these contributions, we obtain To analyze Eq. ( 3) further, we take a generating function approach.With this approach, one trades a probability distribution for a probability generating function whose derivatives provide us with useful information about the distribution such as its moments.Furthermore, it often makes analysis easier, in particular linear analysis.By multiplying Eq. ( 3) by z a and averaging over a, we obtain where ≡ daF (a)z a is the probability generating function of a, Θ(z, t) ≡ daF (a)ρ(a, t)z a , and throughout the paper the superscript (n) represents the nth derivative with respect to ln z.Although Eq. ( 3) is an infinite dimensional system of linear difference equations, Eq. ( 4) is a single difference equation of Θ(z, t) and its derivative [32].
We expand ρ(a, t) as a Maclaurin series as follows: Using this polynomial basis representation (the convergence is proven in the Supplemental Material), we can consider the differentiations in Eq. (4) (i.e., Θ (1) (z, t) and g (1) (z)) as an exchange of bases and convert Eq. ( 4) into a tractable matrix form.Let p 0 be the fraction of initially infected nodes, which are selected uniformly at random, independently of a.We represent the initial condition as w(t = 0) ≡ (w 1 (0), w 2 (0), . ..) ⊤ = (p 0 , 0, 0, . ..) ⊤ .Epidemic dynamics near the epidemic threshold obey linear dynamics given by By substituting Θ(z, t) = ∞ n=1 w n (t)g (n−1) (z) and g (n−1) (1) = a n−1 in Eq. ( 4), we obtain A positive prevalence ρ(t) (i.e., a positive fraction of infected nodes in the equilibrium state) occurs only if the largest eigenvalue of T (τ ) exceeds 1, because in this situation the probability of being infected grows in time, at least in the linear regime.Therefore, we get the following implicit function for the epidemic threshold, denoted by where , s, and u, which are functions of β.In general, we obtain β c by numerically solving Eq. ( 8), but some special cases can be determined analytically.
In the limit τ → 0, Eq. ( 8) gives , which coincides with the epidemic threshold for the activity-driven model derived in the previous studies [19,22].In fact, this β c value is the epidemic threshold for the aggregate (and hence static) network, whose adjacency matrix is given by A * ij ≈ m(a i + a j )/N [3,31], as demonstrated in Fig. S1.
For general τ , if all nodes have the same activity potential a, and if m = 1, we obtain β c as the solution of the following implicit equation: 2ae where κ c = β 2 c + 6β c + 1.The theoretical estimate of the epidemic threshold [Eq.( 8); we use Eq. ( 9) in the case of m = 1] is shown by the solid lines in Figs.2(a) and 2(b).It is compared with numerically calculated prevalence values for various τ and β values shown in different colors.Equations ( 8) and ( 9) describe the numerical results fairly well.When m = 1, the epidemic threshold increases with τ and diverges at τ ≈ 0.1 [Fig.2(a)].Furthermore, slower net-work dynamics (i.e., larger values of τ ) reduce the prevalence for all values of β.In contrast, when m = 10, the epidemic threshold decreases and then increases as τ increases [Fig.2(b)].The network dynamics (i.e., finite τ ) impact epidemic dynamics in a qualitatively different manner depending on m, i.e., the number of concurrent neighbors that a hub has.Note that the estimate of β c by the individual-based approximation ( [31], see Supplemental Material for the derivation), which may be justified when m ≫ 1, is consistent with the numerical results and our theoretical results only at small τ [a dashed line in Fig. 2 To illuminate the qualitatively different behaviors of the epidemic threshold as τ increases, we determine a phase diagram for the epidemic threshold.We focus our analysis on the case in which all nodes share the activity potential value a, noting that qualitatively similar results are also found for power-law distributed activity potentials [Fig.3(b)].We calculate the two boundaries partitioning different phases as follows.First, we observe that the epidemic threshold diverges at τ = τ * .In the limit β → ∞, infection starting from a single infected node in a star graph immediately spreads to the entire star graph, leading to c i → 1 (1 ≤ i ≤ 5).By substituting c i → 1 in Eq. ( 8), we obtain f (τ * , β c → ∞) = 0, where When τ > τ * , infection always dies out even if the infection rate is infinitely large.This is because, in a finite network, infection always dies out after a sufficiently long time due to stochasticity [35][36][37].Second, although β c eventually diverges as τ increases, there may exist τ c such that β c at τ < τ c is smaller than the β c value at τ = 0. Motivated by the comparison between the behavior of β c at m = 1 and m = 10 (Fig. 2), we postulate that τ c (> 0) exists only for m > m c .Then, we obtain dβ c /dτ = 0 at (τ, m) = (0, m c ).The derivative of Eq. ( 8 We simulate the stochastic SIS dynamics using the quasistationary state method [33], as in [31], and calculate the prevalence averaged over 100 realizations after discarding the first 15 000 time steps.We set the step size ∆t = 0.002.Qualitatively similar results are obtained for the variant of the activity-driven model with a reinforcement mechanism of link creation [34] (Fig. S3).
which leads to When m < m c , network dynamics (i.e., finite τ ) always reduce the prevalence for any τ [Figs.10) and ( 11) is shown in Fig. 3(a).The β c values numerically calculated by solving Eq. ( 8) are also shown in the figure.It should be noted that the parameter values are normalized such that β c has the same value for all m at τ = 0. We find that the dynamics of the network may either increase or decrease the prevalence, depending on the number of connections that a node can simultaneously have, extending the results shown in Fig. 2.
These results are not specific to the activity-driven model.The phase diagram is qualitatively similar for randomly distributed m (Fig. S4), for different distributions of activity potentials (Fig. S5), and for a different model in which an activated node induces a clique instead of a star (Fig. S6), modeling a group conversation obeys a power-law distribution with exponent 3 (ǫ ≤ ai ≤ 0.9).We set k = 0.1 at m = 1 and adjust the value of a and ǫ such that βc takes the same value for all m at τ = 0.In the "die out" phase, infection eventually dies out for any finite β.In the "suppressed" phase, βc is larger than the βc value at τ = 0.In the "enhanced" phase, βc is smaller than the βc value at τ = 0.The solid and dashed lines represent τ * [Eq.( 10)] and τc, respectively.The color bar indicates the βc values.In the gray regions, βc > 100.
event as some temporal network models do [38][39][40].Discussion: Our analytical method shows that the presence of network dynamics boosts the prevalence (and decreases the epidemic threshold β c ) when the concurrency m is large and suppresses the prevalence (and increases β c ) when m is small, for a range of values of the network dynamic time scale τ .This result lends theoretical support to previous claims that concurrency boosts epidemic spreading [10,11,[13][14][15][16][17][18][19]41].The result may sound unsurprising because a large m value implies that there exists a large connected component at any given time.However, our finding is not trivial because a large component consumes many edges such that other parts of the network at the same time or the network at other times would be more sparsely connected as compared to the case of a small m.We confirmed that qualitatively similar results are found when the activity potentials were constructed from two empirical social contact networks (Fig. S7).Our results confirm that a monogamous sexual relationship or a small group of people chatting face to face, as opposed to polygamous relationships or large groups of conversations, hinders epidemic spreading, where we compare like with like by constraining the aggregate (static) network to be the same in all cases.
For general temporal networks, immunization strategies that decrease concurrency (e.g., discouraging polygamy) may be efficient.Restricting the size of the concurrent connected component (e.g., size of a conversation group) may also be a practical strategy.
Another important contribution of the present study is the observation that infection dies out for a sufficiently large τ , regardless of the level of concurrency.As shown in Figs. 3 and S6, the transition to the "die out" phase occurs at values of τ that correspond to network dynamics and epidemic dynamics having comparable time scales.This is a stochastic effect and cannot be captured by existing approaches to epidemic processes on temporal networks that neglect stochastic dying out, such as differential equation systems for pair formulationdissolution models [11,[15][16][17][18] and individual-based approximations [31,42,43].Our analysis methods explicitly consider such stochastic effects, and are therefore expected to be useful beyond the activity-driven model (or the clique-based temporal networks analyzed in the Supplemental Material) and the SIS model. .We set m = 5 and a = 0.01.

WHEN THE LOW-ACTIVITY ASSUMPTION IS VIOLATED
Here we consider the situation in which the lowactivity assumption m 2 a ≪ 1 is violated.When m ≪ N , the expected number of star graphs that a star graph overlaps with is given by If p ≪ 1 is violated, a star graph would overlap with others such that the actual concurrency is larger than m.In the extreme case of p ≥ 1, almost all star graphs overlap with each other such that the concurrency is not sensitive to m.In this situation, our results overestimate the epidemic threshold because our analysis does not take into account infections across different star graphs.The rates of the infection events are given by M {I,S,z},{S,S,z} = zβ, (S6) M {I,I,z},{S,I,z} = (z + 1)β, (S7) M {I,I,z},{I,S,z} = β, (S8) The other elements of M are equal to 0. Let p {x,y,z} (t) be the probability for a star graph to be in state {x, y, z} at time t.Because where p(t) is the 4m-dimensional column vector whose elements are p {x,y,z} (t), we obtain Note that c 1 and c 2 are the probabilities with which x = I at time τ , when the initial state is {I, S, 0} and {S, I, 0}, respectively, and that c 3 , c 4 , and c 5 are the probabilities that y = I at time τ , when the initial state is {S, I, 0}, {I, S, 0}, and {S, S, 1}, respectively.Therefore, we obtain When m = 1, Eq. (S12) yields where κ = β 2 + 6β + 1, and c 5 is not defined.When m ≫ 1, we can apply an individual-based approximation [1,4,5].We assume that the state of each node is statistically independent of each other, i.e., p {x,y,z} ≈ P (x)P (y)P (z), (S15) where P (x), for example, is the probability that the hub takes state x.We have suppressed t in Eq. (S15).Under the individual-based approximation, x and y obey Bernoulli distributions with parameters p MF where In a similar fashion to the derivation of Eq. (S12), we obtain We estimate the extent to which Eq. (S20) is valid as follows.First, we need m ≫ 1, because the initial condition p MF ≪ 1, we need τ < 1/β.This condition remains unchanged by rescaling (τ, β) to (cτ, β/c).These two conditions are sufficient for this approximation to be valid.If m ≫ 1 is violated, the individual-based approximation significantly underestimates the epidemic threshold for any finite τ because it ignores the effect of stochastic dying-out.If τ < 1/β is violated, the approximation [dashed lines in Fig. 2(b) and (d)] underestimates the epidemic threshold because dynamics on the star graph deviate from the linear regime.In particular, the epidemic threshold obtained from the approximation [Eq.(S48)] remains finite even in the limit τ → ∞, whereas analytical [Eq.( 8)] and numerical (Fig. 2) results diverge at a finite τ .

DERIVATION OF EQ. (8)
At the epidemic threshold, the largest eigenvalue of T is equal to unity.Let v = (v 1 , v 2 , . ..) ⊤ be the corresponding eigenvector of T .We normalize v such that ∞ j=1 v j = 1.By substituting Eq. ( 7) in v = T v, we obtain the system of equations Equation (S24) gives where By combining Eqs.(S23) and (S25), we obtain where where, By substituting Eq. (S30) in Eq. (S31), we obtain − qr − qs + qrs − q 2 u − rs = 0, (S34) which is Eq. ( 8) in the main text.If all nodes have the same activity potential a, Eq. (S34) is reduced to
Second, consider a finite t.It should be noted that the series is only defined at t that is a multiple of τ .Because T ij = 0 (i ≥ j + 2) in Eq. ( 7), we obtain Therefore, the series converges.Third, we consider the limit t → ∞.where v is the eigenvector of T given by Eq. (S30), and b is a constant.Because Eq. (S30) yields the radius of convergence is equal to a /q.To ensure convergence, we require that max Equation (S47) holds true in practical situations because the assumption m 2 a ≪ 1 guarantees that m a ≪ 1 and a i is a probability.

EPIDEMIC THRESHOLD UNDER THE INDIVIDUAL-BASED APPROXIMATION
When m ≫ 1, the epidemic threshold can be obtained by the individual-based approximation [1,4,5].We assume that all nodes have the same activity potential a.By substituting Eq. (S20) in Eq. (S35), we obtain Equation (S48) agrees with the value derived in [1].Note that this approximation is valid only for small τ (τ < 1/β c ).

DERIVATION OF mc FOR GENERAL ACTIVITY DISTRIBUTIONS
At m = m c , an infinitesimal increase in τ from 0 to ∆τ does not change the β c value.For general activity distributions, by setting ∂f /∂τ = 0 for f given by Eq. (S34), We carried out numerical simulations for an extended activity-driven model in which link dynamics are driven by a reinforcement process [6].The original activitydriven model is memoryless [7].In the extended model, an activated node i connects to a node j that i has already contacted with probability 1/(n i +c) and to a node j that i has not contacted with probability c/(n i + c), where n i denotes the number of nodes that node i has already contacted.
(b)].Qualitatively similar results are found, when the activity potential a is power-law distributed [Figs.2(c) and 2(d)].

FIG. 2 .
FIG. 2. Epidemic threshold and the numerically-simulated prevalence when m = 1 (a),(c) and m = 10 (b),(d).In (a) and (b), all nodes have the same activity potential value a.The solid lines represent the analytical estimate of the epidemic threshold [Eq.(8); we plot Eq. (9) instead in (a)].The dashed lines represent the epidemic threshold obtained from the individual-based approximation (Supplemental Material).The color indicates the prevalence.In (c) and (d), the activity potential (ǫ ≤ ai ≤ 0.9, 1 ≤ i ≤ N ) obeys a power-law distribution with exponent 3.In (a)-(d), we set N = 2000 and adjust the values of a and ǫ such that the mean degree is the same ( k = 0.1) in the four cases.We simulate the stochastic SIS dynamics using the quasistationary state method[33], as in[31], and calculate the prevalence averaged over 100 realizations after discarding the first 15 000 time steps.We set the step size ∆t = 0.002.Qualitatively similar results are obtained for the variant of the activity-driven model with a reinforcement mechanism of link creation[34] (Fig.S3).

FIG. 3 .
FIG.3.Phase diagrams for the epidemic threshold, βc, when the activity potential is (a) equal to a for all nodes, or (b) obeys a power-law distribution with exponent 3 (ǫ ≤ ai ≤ 0.9).We set k = 0.1 at m = 1 and adjust the value of a and ǫ such that βc takes the same value for all m at τ = 0.In the "die out" phase, infection eventually dies out for any finite β.In the "suppressed" phase, βc is larger than the βc value at τ = 0.In the "enhanced" phase, βc is smaller than the βc value at τ = 0.The solid and dashed lines represent τ * [Eq.(10)] and τc, respectively.The color bar indicates the βc values.In the gray regions, βc > 100.

Supplemental 2 − 1
FIG. S1.Prevalence on the aggregate (hence static) network whose adjacency matrix is given (in the limit N → ∞) by A * ij = m(ai + aj)/N [1, 2].The lines represent the numerical results for the delta function (i.e., all nodes have same activity potential) and power-law activity distributions.The arrows indicate βc = m a + a 2 −1

If p ≥ 1 ,
the individual-based approximation describes the numerical results more accurately than our method does [Figs.S2(c) and S2(d)].However, even at a moderately large value of p (= 0.5), our method is more accurate than the individual-based approximation [Figs.S2 (a) and S2(b)].
FIG. S2.Epidemic threshold and numerically calculated prevalence when the low-activity assumption is violated.We set m = 1 in (a) and (c), m = 10 in (b) and (d), p = 0.5 in (a) and (b), and p = 1.5 in (c) and (d).The solid and dashed lines represent the epidemic threshold obtained from Eq. (8) and that obtained from the individual-based approximation, respectively.All nodes are assumed to have the same activity potential a = 0.25 in (a), a = 0.0045 in (b), a = 0.75 in (c), and a = 0.0136 in (d).We calculated the prevalence averaged over 100 simulations after discarding the first 15000 time steps of each simulation.We set N = 1000 and ∆t = 0.002.
The numerically calculated prevalence is compared between the original model [Figs.S3(a) and S3(b)] and the extended model with c = 1 [Figs.S3(c) and S3(d)].We replicate Figs.2(a) and 2(b) in the main text as Figs.S3(a) and S3(b) as reference.All nodes are assumed to have the same activity potential a i = 0.05 (1 ≤ i ≤ N ) in (a) and (c) and a i = 0.005 (1 ≤ i ≤ N ) in (b) and (d).
FigureS3indicates that the extended model only slightly changes the epidemic threshold.

2 FIG
FIG. S7. Results for activity potentials derived from empirical data.The epidemic threshold and numerically simulated prevalence are shown for m = 1 (a),(c) and m = 10 (b),(d).In (a) and (b), the activity potential is constructed from contact data obtained from the SocioPatterns project[8].This data set contains contacts between pairs of N = 92 individuals measured every 20 seconds.In (c) and (d), the activity potential is constructed from email communication data at a research institution, obtained from the Stanford Network Analysis Platform[9].Although the original edges are directed, we treat them as undirected.We assume that each email exchange event corresponds to a one-minute contact.We calculate the degree of each node per minute averaged over time, denoted by ki , and define the activity potential as ai = [ ki − k /2] /m.In (c) and (d), we used N = 439 individuals satisfying ai > 0 (some individuals exchanged few emails such that ai < 0).We set ∆t = 0.001.