Effects of concurrency on epidemic spreading in Markovian temporal networks

The concurrency of edges, quantified by the number of edges that share a common node at a given time point, may be an important determinant of epidemic processes in temporal networks. We propose theoretically tractable Markovian temporal network models in which each edge flips between the active and inactive states in continuous time. The different models have different amounts of concurrency while we can tune the models to share the same statistics of edge activation and deactivation (and hence the fraction of time for which each edge is active) and the structure of the aggregate (i.e., static) network. We analytically calculate the amount of concurrency of edges sharing a node for each model. We then numerically study effects of concurrency on epidemic spreading in the stochastic susceptible-infectious-susceptible and susceptible-infectious-recovered dynamics on the proposed temporal network models. We find that the concurrency enhances epidemic spreading near the epidemic threshold while this effect is small in many cases. Furthermore, when the infection rate is substantially larger than the epidemic threshold, the concurrency suppresses epidemic spreading in a majority of cases. In sum, our numerical simulations suggest that the impact of concurrency on enhancing epidemic spreading within our model is consistently present near the epidemic threshold but modest. The proposed temporal network models are expected to be useful for investigating effects of concurrency on various collective dynamics on networks including both infectious and other dynamics.


Introduction
The structure of contact networks among individuals shapes the dynamics of contagion processes in a population such as the number of individuals infected, their spatial distributions, and speed of spreading [1][2][3][4].For example, a heterogeneous degree distribution, where the degree is the number of neighboring nodes that a node has, and a short average distance between nodes are two factors that usually enhance epidemic spreading on networks.In fact, empirical contact networks often vary over time on a time scale comparable to or faster than that of epidemic dynamics, probably most famously owing to mobility of human or animal individuals.This observation has naturally led to the investigation of how features of such temporal (i.e., time-varying) networks and statistical properties of contact events, such as distributions of inter-contact times, temporal correlation in inter-contact times, and appearance and disappearance of nodes and edges, affect outcomes of contagion processes [5][6][7][8][9][10].
The concurrency is a feature of temporal contact networks that has been investigated in both mathematical and field epidemiology for over two decades, while many of these studies do not specifically refer to temporal networks [11][12][13][14][15][16][17][18][19].Concurrency generally refers to multiple partnerships of an individual that overlap in time.Concurrency has been examined in particular in the context of sexually transmitted infections such as HIV, specifically regarding whether or not the concurrency was high in sub-Saharan Africa and, if so, whether or not the high concurrency enhanced HIV spreading there [11][12][13][20][21][22][23].Suppose that a partnership between individuals i and j overlaps one between i and j ′ for some time.Such concurrency may promote epidemic spreading because rapid transmission from j to j ′ through i and vice versa is possible during the concurrent partnerships (see Fig. 1(a)).In contrast, if a partnership between i and j occurs before that between i and j ′ without an overlap, i.e., without concurrency, the transmission from j to j ′ can still occur in various different ways, but the transmission from j ′ to j does not (see Fig. 1(b)).
Modeling studies are diverse in how to quantify the amount of concurrency.Many studies measure the concurrency through the mean degree, or the average contact rate for nodes [19,[24][25][26].However, in these cases, it is hard to say whether an increased epidemic size owes to a high density of edges or a large amount of concurrency [15,[27][28][29].In studies of epidemic processes on static networks, it is a stylized fact that epidemic spreading is enhanced if there are many edges with all other things being equal.Another line of theoretical and computational approach to concurrency quantifies concurrency by the heterogeneity in the degree distribution of the network [12,13,15].However, the results that a high concurrency in this sense enhances epidemic spreading is equivalent to an established result in static network epidemiology that heterogeneous degree distributions lead to a small epidemic threshold, thus promoting epidemic spreading [2][3][4]30].In some studies, the authors investigated the effect of concurrency by carefully ensuring that the degree distribution (and hence the average degree) stays the same across the comparisons while they manipulate the amount of concurrency [27][28][29].These studies suggest that, despite the different definitions of the concurrency employed in these studies, higher concurrency increases epidemic spreading (but see [27] for different results).However, mathematical models that enable us to analyze the effect of concurrency by fixing the structure of the static network are still scarce.[15] and Fig. 1 in Ref. [27].
In the present study, we focus on effects of concurrency on epidemic spreading for an arbitrary static network on top of which partnerships, infection events, and recovery events occur.In this manner, we aim to study the effect of concurrency without being affected by the effect of the network structure itself.We consider three temporal network models that have different amounts of concurrency while keeping the probability that each edge is available the same across comparisons.In terms of a concurrency measure, we analytically evaluate the amount of concurrency of edges sharing a node for the three models as well as other properties of the models.Then, we numerically study effects of concurrency on epidemic spreading using the stochastic susceptible-infectious-susceptible (SIS) and susceptible-infectious-recovered (SIR) models.
The Python codes for generating the numerical results in this article are available at Github (https: //github.com/RuodanL/concurrency).

Temporal network models
We introduce three models of undirected and unweighted temporal networks in continuous time, which we build from independent Markov processes occurring on a given static network.We use these models to compare temporal networks that have the same time-aggregated network but different amounts of concurrency.Let G be an undirected and unweighted static network with node set V = {1, ..., N } and edge set E = {e i } M i=1 , where N is the number of nodes and M is the number of edges.We assume that the network G has no self-loops.In the temporal network constructed based on undirected and unweighted static network G, by definition, each edge e i ∈ E is either active (i.e., temporarily present) or inactive (i.e., temporarily absent) at any given time t ∈ R and switches between the active and inactive states over time.We denote by G the temporal network and by G(t) the instantaneous network in G observed at time t.

Model 1
Let τ 1 and τ 2 be the duration of the active state of and the inactive state of an edge, respectively.In our model 1, we assume that τ 1 is independently drawn from probability density ψ 1 (τ 1 ) each time the edge switches from the inactive to the active state.Similarly, we draw τ 2 independently from probability density ψ 2 (τ 2 ) each time the edge switches from the active to the inactive state.The duration τ 1 or τ 2 for different edges obeys the same distributions but is drawn independently for the different edges.When ψ 1 (τ 1 ) and ψ 2 (τ 2 ) are both exponential distributions, the stochastic dynamics of the state of each edge obeys independent continuous-time Markov processes with two states, and our model 1 reduces to previously proposed models [31,32].

Model 2
In models 2 and 3, we assume that each node independently obeys a continuous-time Markov process with two states, which we refer to as the high-activity and low-activity states.Each node independently switches between the high-activity state, denoted by h, and the low-activity state, denoted by ℓ.Let a be the rate at which the state of a node changes from ℓ to h, and let b be the rate at which the state of a node changes from h to ℓ.In model 2, each edge (u, v) ∈ E, where u, v ∈ V , is active at any given time if and only if both u and v are in the h state.The dynamics of a single edge in model 2 is schematically shown in Fig. 2. The model resembles the so-called AND model in [33].An intuitive interpretation of model 2 is that two individuals chat with each other if and only if both of them want to.

Model 3
Model 3 is a variant of model 2. As in model 2, we assume that each node independently switches between the high-activity and low-activity states at rates a and b.In model 3, the edge between two nodes is active if and only if either node is in the h state (see Fig. 2 for a schematic).This model resembles the OR model proposed in [33].The intuition behind model 3 is that one person can start conversation with another person whenever either person wants to talk regardless of whether the other person wants to.

Concurrency for the three models
In this section, we first define a concurrency index for temporal networks in which partnership on each edge appears and disappears over time.Then, we calculate and compare the concurrency index and related measures for models 1, 2, and 3.

Definition of a concurrency index
Consider a temporal network G that is constructed on the underlying static network G and defined for time t ∈ [0, T ].We refer to G as the aggregate network.If two edges in G sharing a node in G are both active at t ∈ R, we say that the two edges are concurrent at time t; see Fig. 1(a) for an example.Consider a pair of edges sharing a node in the aggregate network, denoted by e i and e j , where e i , e j ∈ E. In fact, the likelihood that e i and e j are concurrent at time t depends on how likely each single edge is active at time t.We define a concurrency index that takes into account of this factor.To this end, we first define the set of time for which an edge e ∈ E is active by S(e) = {t ∈ [0, T ]; edge e is active at time t}. (3.1) We define the concurrency for the edge pair {e i , e j } ∈ S by κ(e i , e j ) = m (S(e i ) ∩ S(e j )) min{m (S(e i )) , m (S(e j ))} , where m denotes the Lebesgue measure, and S is the set of edge pairs that share a node in G.Note that where k i is the degree of the ith node in G [12,13].The numerator on the right-hand side of Eq. (3.2) is equal to the length of time for which both e i and e j are active.The denominator is a normalization constant to discount the fact that the numerator would be large if the two edges are active for long time.Because any edge e in the aggregate network G should be active sometime in [0, T ], the denominator is always positive; otherwise, we should exclude e from G.
We define the concurrency index for temporal network G by i,j such that 1≤i<j≤M and {ei,ej }∈S κ(e i , e j ).(3.4)It holds true that 0 ≤ κ(G) ≤ 1 because 0 ≤ κ(e i , e j ) ≤ 1.For empirical or numerical data, [0, T ] is the observation time window.For a stochastic temporal network model, we calculate S(e) as the expectation and in the limit of T → ∞ such that κ(G) is a deterministic quantity.Differently from κ 3 , a concurrency index proposed in seminal studies [12,13] (also see for [15] a review), κ(G) is not affected by the degree distribution of the aggregate network G.It should be noted that the calculation of κ(e i , e j ) and hence κ(G) requires the information about the aggregate network, i.e., S.

Model 1
To calculate the concurrency index for the three models of temporal networks, we first derive the probability of an arbitrary edge being active in the equilibrium for each model.In model 1, consider an arbitrary edge in the static network G.The mean duration for which the edge is active and that for which the edge is inactive are given by and respectively.Owing to the renewal reward theorem [34], the probability that the edge is active in the equilibrium, denoted by q * , is given by The concurrency index, κ(G), which we simply refer to as κ in the following text, is equal to the ratio of the time for which the two edges sharing a node are both active to the time for which an edge is active.Because the states of different edges are independent of each other, we obtain The special case of model 1 in which the edge activation and deactivation occur as Poisson processes is equivalent to previously proposed models [31,32].In this case, using ψ 1 (τ 1 ) = λ 1 e −λ1τ1 and ψ 2 (τ 2 ) = λ 2 e −λ2τ2 , we obtain (3.9)

Model 2
To analyze model 2, let us consider a pair of neighboring nodes v 1 and v 2 in the static network G. Let p * h and p * ℓ be the probability that an arbitrary node in G is in state h and ℓ in the equilibrium, respectively.Denote by p * s1s2 the probability that node v i is in state s i ∈ {h, ℓ} in the equilibrium, where i ∈ {1, 2}.Because the duration of the high-activity state and that of the low-activity state of a node obey the exponential distributions with mean 1/b and 1/a, respectively, we apply the renewal reward theorem [34] to obtain and Because the states of different nodes are independent, we obtain (3.12) and Therefore, the probability that an edge is active in the equilibrium is given by To derive the concurrency for model 2, we consider a pair of edges sharing a node in G, denoted by (v 1 , v 2 ) and (v 2 , v 3 ).Denote by p s1s2s3 the probability that node v i is in state s i ∈ {h, ℓ}, where i = 1, 2, and 3. Similar to the analysis of q * for model 2, we obtain the following stationary probabilities: ) Therefore, we obtain (3.20)

Model 3
Similarly, for model 3, we obtain and

Comparison among the three models
A large fluctuation of G(t) over time may impact concurrency [15].In this section, we compare the amount of concurrency between the three models.Proposition 1. Model 2 is more concurrent than model 1 given that q * is the same between the two models.
Proposition 2. Model 3 is more concurrent than model 1 given that q * is the same between the two models.
Proof.Because q * is the same between the two models and 0 < q * < 1, we obtain Therefore, model 3 is more concurrent than model 1.
Proposition 3. Given that q * is the same between models 2 and 3, For the sake of the present proof, let a 2 and a 3 be the rate at which the state of a node changes from ℓ to h for models 2 and 3, respectively.Likewise, let b 2 and b 3 be the rate at which the state of a node changes from h to ℓ for models 2 and 3, respectively.By imposing that q * is the same between the two models, we obtain Therefore, the difference in the concurrency between the two models, given by Eqs.(3.20) and (3.22), is given by Therefore, the concurrency of model 3 is larger than that of model 2 if and only if Using Eqs.(3.8), (3.20), and (3.22), we compare the amount of concurrency for models 1, 2, and 3 in Fig. 3.We find that, when q * > 1/2, the concurrency index for model 2 is only slightly larger than that for model 3. Figure 3: Concurrency index, κ, as a function of the stationary probability that an edge is active, q * , for models 1, 2, and 3.

Fluctuations in the node's degree
Consider the degree of individual nodes or its average over all nodes at any given time t.Let us fix their expectation and compare their statistical fluctuations, such as the standard deviation, across different models.In a model with a higher concurrency, edges sharing a node in the aggregate network are more likely to be simultaneously active or inactive, which makes the statistical fluctuation of the degree larger.Therefore, we anticipate that the concurrency of a temporal network model and the statistical fluctuation in the node's degree are interrelated.In this section, we establish such relationships for each of the three models.Let A = (A ij ) N ×N be the adjacency matrix of static network G.
To analyze the fluctuation in the node's degree in model 1, we define random variables by The adjacency matrix of the temporal network at time t, G(t), of model 1 is given by the N × N matrix be the degree of the ith node in network G(t).The average degree at time t is given by In the following text, we omit t because we discuss the fluctuations in k i (t) and ⟨k⟩(t) in the equilibrium.
Proposition 4. For model 1, it holds true that k i ∼ B(k i , q * ) and ⟨k⟩ ∼ 2 N B (M, q * ), where B(•, •) represents the binomial distribution.(We remind that k i is the degree of the ith node in G and that M is the number of edges in G.) Proof.For a proposition C, define the indicator function by For a node i, we obtain Because x ij are independent and identically distributed Bernoulli random variables and k i is the sum of k i terms, k i obeys B k i , q * .Because we have assumed that the network is undirected, we obtain Because there are M terms that comprise the summation, ⟨k⟩ obeys 2 N B (M, q * ).We denote the expectation by E and the variance by σ 2 ; the standard deviation is equal to σ.Using Proposition 4, we obtain ) To analyze models 2 and 3, we define The adjacency matrix of G(t) for model 2 is given by Using this expression, we can prove the following proposition.
Proposition 5.For model 2, we obtain (3.39) We prove Proposition 5 in Appendix A.
The adjacency matrix of G(t) for model 3 is given by Using this expression, we can prove the following proposition.
Proposition 6.For model 3, we obtain ) We prove Proposition 6 in Appendix B. Now, let us compare the variance of k i and ⟨k⟩ among the three models under the condition that the expectation of k i and ⟨k⟩ is the same across the different models.This condition is equivalent to keeping q * the same across the models for each edge.The purpose of examining the variance of the node's degree is the following.Similar to Eq. (3.3), the number of concurrent edge pairs at time t is given by Because this expression contains the second moment of the degree, we expect that the variance of the degree may be related to our concurrency measure.
In fact, we obtain the following results for the fluctuation of the degree, which are parallel to those for the concurrency index.
Proposition 7. Assume that E[k i ] is the same between models 1 and 2 for any i ∈ {1, . . ., N }.For any given q * , the variance σ Proof.We use subscripts 1 and 2 to represent the variance with respect to the probability distribution for models 1 and 2, respectively.We substitute Eq. (3.15) in Eq. (3.32) and use Eq.(3.37) to obtain Next, we compare the variance of the average degree.Because M is the number of edges of static network G and there exists i such that k i > 1 in G, we obtain for the following reason.The right-hand side of Eq. (3.45) is the number of pairs of edges.The left-hand side is the number of pairs of edges that do not share a node.These two quantities are equal to each other if and only if there is no pair of edges sharing a node, i.e., when all nodes have the degree at most 1.By substituting Eq. (3.15) in Eq. (3.34) and using Eqs.(3.39) and (3.45), we obtain is the same between models 1 and 3 for any i ∈ {1, . . ., N }.For any given q * , the variance σ 2 [k i ] is larger for model 3 than model 1 if k i > 1.Likewise, σ 2 [⟨k⟩] is larger for model 3 than model 1 if there exists i such that k i > 1.
Proof.We substitute Eq. (3.21) in Eq. (3.32) and use Eq.(3.41) to obtain (3.48) Proposition 9. Assume that E[k i ] is the same between models 2 and 3 for any i ∈ {1, . . ., N }.For any given q * , if k i > 1, we obtain ⟨k⟩] satisfies the same relationships if there exists i such that k i > 1.
We prove Proposition 9 in Appendix C.

Duration for the edge being inactive in model 2
In model 2, an edge is active if and only if both nodes forming the edge are in the h state, and the state of each node (i.e., h or ℓ) independently obeys a continuous-time Markov process with two states.Therefore, the duration of the edge being active obeys an exponential distribution with rate 2b.In contrast, the duration of the edge being inactive does not obey an exponential distribution, which we characterize as follows.
Proposition 10.The probability density function of the duration of the edge being inactive in model 2 is the mixture of two exponential distributions given by where ) and Proof.Consider a three-state continuous-time Markov process described as follows.Let z (with z= 0, 1, or 2) denote the number of nodes forming an edge that are in the h state.We refer to the value of z as the state of the three-state Markov chain without arising confusion with the single node's state (i.e., h or ℓ) or edge's state (i.e., active or inactive).We initialize the stochastic dynamics of nodes by setting z = 2, which corresponds to the edge being active.Consider a sequence of the state z that starts from z = 2 at time 0 and return to z = 2 for the first time.Let I n be such a sequence of the z values visiting z = 0 in total n times before returning to z = 2 for the first time, which we denote by The duration of the edge being inactive is the difference between the time of the first passage to z = 2 and the time of leaving z = 2 last time.We denote this duration by T n for the case in which z = 0 is visited n times before z = 2 is revisited for the first time.
For general n, we obtain where τ ′ 1,i is the ith duration of z = 1, which obeys the exponential distribution with rate a + b; τ ′′ 0,i is the ith duration of z = 0, which obeys the exponential distribution with rate 2a.Variables τ ′ 1,1 , . .., τ ′ 1,n+1 , τ ′′ 0,1 , . .., τ ′′ 0,n are independent of each other.Therefore, the Laplace transform of the distribution of T n is given by (3.55) The probability that sequence I n occurs is given by Let T be the duration for which the edge is inactive.Using Eqs.(3.55) and (3.56), we obtain the Laplace transform of the distribution of T as follows: 2 .
(3.57) Therefore, We verify Eq. (3.49) with numerical simulations for two parameter sets.The results are shown in Fig. 4. The dashed curves represent the exponential distributions whose mean is the same as that of the corresponding mixture of two exponential distributions.The figure suggests that the actual duration for the edge to be inactive is distributed more heterogeneously than the exponential distribution for both parameter sets.For model 3, the duration for the edge being inactive, which is equivalent to the time for which the same three-state Markov process stays in state z = 0, obeys an exponential distribution with rate 2a.The duration for the edge being active is the first passage time to z = 0 since the Markov chain has left z = 0. Therefore, the duration for the edge being active for model 3 obeys the mixture of two exponential distributions given by Eq. (3.49), but with a and b being swapped.

Methods
To examine the effect of concurrency on epidemic spreading, we run the stochastic SIS and SIR models on three static networks.On each static network, we run stochastic dynamics of partnership according to model 1, 2, or 3 and compare the extent of epidemic spreading among the models.Each node takes either the susceptible, infectious, or recovered state, and the node's state may change over time.In both SIS and SIR models, an infectious node infects each of its susceptible neighbor independently at rate β, which we call the infection rate, and an infectious node recovers at rate µ, which we call the recovery rate, independently of the other nodes' states.The infection and recovery events occur as Poisson processes with the respective rates.Once an infectious node recovers, it turns into the susceptible state in the case of the SIS model and the recovered state in the case of the SIR model.Recovered nodes in the SIR model do not infect other nodes and are not infected by other nodes.
We consider three static networks.First, we use the Erdős-Rényi (ER) random graph with N = 200 nodes.We independently connect each pair of nodes with probability 0.05.We iterated generating a network from the ER random graph until we obtained a connected network with M = 1000 edges.The second network is a network with N = 200 nodes that has a degree distribution with a power-law tail that is proportional to k −3 , where k is the degree, generated by the Barabási-Albert (BA) model [35].We assume that each incoming node has five edges to be connected to already existing nodes according to the proportional preferential attachment rule.In other words, the probability that an existing node, denoted by v, forms a new edge with an incoming node is proportional to v's degree.The initial network is a star graph on 6 nodes.The network is connected and contains M = 975 edges.Third, we use the largest connected component of a collaboration network among researchers who had published papers in network science up to 2006 [36].The network has N = 379 nodes and M = 914 edges.The edge represents the presence of at least one paper that two authors have coauthored.We use this network as an unweighted network.
We set µ = 1 without loss of generality; multiplying the same constant to β, µ, a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 only rescales the time; a i and b i represent parameters a and b for model i.We run 5000 simulations for each epidemic process model (i.e., SIS or SIR), each static network, each partnership model (i.e., model 1, 2, or 3), and each parameter set.In each simulation of the SIS model, we started from the initial condition in which all nodes are infectious.This initial condition is not realistic.However, we use it with the aim of identifying the quasi-stationary fraction of infectious nodes without being affected by rapid extinction of infectious nodes that may happen with an initial condition having only a small number of infectious nodes even if the infection rate is large.In the SIS model, we measure the fraction of infectious nodes at time t = 10 4 and average it over 10 3 simulations at each β value.In each simulation of the SIR model, just one node selected uniformly at random is initially infectious, and all the other nodes are susceptible.For each static network, the initially infectious node is the same over the different partnership models and the different β values.For each β value, we average the fraction of recovered nodes at the end of the simulation over the 5000 simulations, which we call the final epidemic size.
We run the simulation with model 1 using the Laplace Gillespie algorithm [37,38], which is an extension of the direct method of Gillespie to the case in which inter-event times are allowed to obey a mixture of exponential distributions like Eq. (3.49).For models 2 and 3, we run the simulation using the direct method of Gillespie because these models only involve exponential distributions of inter-event times.
Apart from the variation in the β value, we consider four sets of parameter values for both SIS and SIR models.In the first set of simulations, we set a 1 = 1, b 1 = 9, a 2 = 0.5( √ 10 + 1) ≈ 2.08, b 2 = 4.5, a 3 = 0.5, and b 3 = 1.5( √ 10 + 3) ≈ 9.24.Then, we obtain q * = 0.1 for all the three partnership models.In this manner, we compare the final epidemic size for the three partnership models, which yield different amounts of concurrency, under the condition that each edge is active for the same amount of time on average.It should be noted that a large q * value will obviously lead to a larger final epidemic size with the other things being equal, and so it is necessary to compare the models with the q * value being fixed.We obtain a second parameter set by making the values of a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 for the first parameter set five times smaller, which implies that the nodes and edges flip their states five times more slowly than in the case of the first parameter set.Because multiplying these six parameters by the same constant does not change q * , we retain q * = 0.1 for the second parameter set.We refer to the first and second parameter sets as the fast and slow edge dynamics, respectively.The third parameter set is defined by 21, b 2 = 0.5, a 3 = 0.5, and b 3 = 0.5( √ 2 + 1) ≈ 1.21 such that q * = 0.5.We also consider a slower variant of edge dynamics by dividing these a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 values by five.In summary, for each of the three networks and each partnership model, we have four cases each of which consists of the combination of q * ∈ {0.1, 0.5} and either the fast or slow edge dynamics.

Results for the SIS model
In this section, we examine the SIS model with partnership models 1, 2, and 3. We first compare between partnership models 1 and 2. We recall that model 2 has a higher concurrency than model 1 when each edge is active with the same probability in the two models.To exclude the possible effects of the distribution of the duration for which the edge is active and that for which the edge is inactive, we use model 1 with ψ 1 (τ 1 ) = 2b 2 e −2b2τ1 and ψ 2 (τ 2 ) = C 1 λ 1 e −λ1τ2 +C 2 λ 2 e −λ2τ2 , where λ 1 , λ 2 , C 1 , and C 2 are given by Eqs.(3.50), (3.51), (3.52), and (3.53), respectively, with a = a 2 and b = b 2 .In this manner, models 1 and 2 have the identical distribution of the duration of the edge being active (i.e., ψ 1 ) and that of the edge being inactive (i.e., ψ 2 ), whereas they are different in terms of the amount of concurrency.
In Fig. 5(a)-(d), we show the relationships between the infection rate and fraction of infectious nodes at t = 10 4 for a network generated by the ER random graph.Figure 5(a) and (b) corresponds to q * = 0.1; Fig. 5(c) and (d) corresponds to q * = 0.5.Figure 5(a) and (c) corresponds to the slow edge dynamics; Fig. 5(b) and (d) corresponds to the fast edge dynamics.Figure 5(a)-(d) indicates that, in each of these combinations of a value of q * and speed of the edge dynamics, the epidemic threshold is smaller for model 2 than model 1.As a result, the fraction of infectious nodes is larger for model 2 than model 1 when the infection rate is near the epidemic threshold.These results suggest that concurrency decreases the epidemic threshold and promotes epidemic spreading when the infection rate is near the epidemic threshold, which is consistent with the main conclusions of previous studies [28,29,39].However, the epidemic threshold and the fraction of infectious nodes near the epidemic threshold are almost the same between models 1 and 2 when q * = 0.5 (see Fig. 5(c) and (d)).Furthermore, the concurrency decreases the fraction of infectious nodes when the infection rate is sufficiently larger than the epidemic threshold, with both q * = 0.1 and q * = 0.5.

Infection rate Fraction of infectious nodes
We show the fraction of infectious nodes at t = 10 4 for models 1 and 3 for the same three underlying static networks, two values of q * , and two speeds of edge dynamics in Fig. 6.The results are qualitatively similar to the comparison between models 1 and 2 shown in Fig. 5.These results with models 1 and 3 further support that concurrency decreases the epidemic threshold, whereas this effect is small for the BA model and collaboration network as well as for denser networks (i.e., q * = 0.5 compared to q * = 0.1).(a) ER, q * = 0.1, slow (b) ER, q * = 0.1, fast (c) ER, q * = 0.5, slow ER, q * = 0.5, fast BA, q * = 0.1, slow (f) BA, q * = 0.5, slow BA, q * = 0.5, fast (i) Collabo., q * = 0.1, slow (j) Collabo., q * = 0.1, fast (k) Collabo., q * = 0.5, slow (l) Collabo., q * = 0.5, fast In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.

Results for the SIR model
In this section, we examine the SIR model with the three partnership models.When comparing between models 1 and 2 and between models 1 and 3, we use the same distributions of inter-event times for model 1 as those used in section 4.1.2.
In Fig. 7(a)-(d), we show the relationships between the infection rate and final epidemic size for a network generated by the ER random graph.Figure 7(a) and (b) corresponds to q * = 0.1; Fig. 7(c) and (d) corresponds to q * = 0.5.Figure 7(a) and (c) corresponds to the slow edge dynamics; Fig. 7(b) and (d) corresponds to the fast edge dynamics.Figure 7(a)-(d) indicates that, in each of these combinations of a value of q * and speed of the edge dynamics, the final epidemic size is larger for model 2 than model 1 when the infection rate is near the epidemic threshold.This result suggests that concurrency promotes epidemic spreading near the epidemic threshold.However, the final epidemic size near the epidemic threshold is only slightly larger for model 2 than model 1 when q * = 0.5 (see Fig. 7(c) and (d)).Furthermore, the concurrency suppresses the epidemic spreading when the infection rate, and hence the final epidemic size, is larger for both q * = 0.1 and q * = 0.5.All these results are similar to those for the SIS model shown in Fig. 5(a)-(d) except that there is no noticeable difference in the epidemic threshold between models 1 and 2 in the case of the SIR model.
The results for the BA model, shown in Fig. 7(e)-(h), are similar to those for the ER random graph.The results for the collaboration network, shown in Fig. 7(i)-(l), are similar to those for the ER random graph and the BA model except that, when q * = 0.1, the final epidemic size is larger for model 2 than model 1 across the entire range of the infection rate, β, that we investigated; see Fig. 7(i) and (j).Because the final epidemic size has approximately plateaued for both models 1 and 2 at the largest β value that we investigated, it is unlikely that model 1 yields a larger final epidemic size than model 2 when β is larger.These numerical results suggest that the effect of concurrency on enhancing epidemic spreading at q * = 0.1 is stronger for the collaboration network than for the ER and BA models.In contrast, when q * = 0.5, the final epidemic size is not noticeably larger for model 2 than model 1 when β is near the epidemic threshold and smaller for model 2 than model 1 when β is larger (see Fig. 7(k) and (l)).This result is qualitatively the same as that for the ER and BA models shown in Fig. 7(c), (d), (g), and (h).

Infection rate
Final epidemic size In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.The inset of (i) shows the magnification of the main panel because the final epidemic size is small in the entire range of the infection rate in this case.
In Fig. 8, we compare the final epidemic size in the SIR model between models 1 and 3 for the same three underlying static networks, two values of q * , and two speeds of edge dynamics.The results are largely similar to the comparison between models 1 and 2. A notable difference from the comparison between models 1 and 2 is that, for the ER random graph and BA model combined with q * = 0.1 and the slow edge dynamics, model 3 yields a larger final epidemic size than model 1 across the entire range of the infection rate values investigated (see Fig. 8 BA, q * = 0.5, slow (h) BA, q * = 0.5, fast (i) Collabo., q * = 0.1, slow (j) Collabo., q * = 0.1, fast (k) Collabo., q * = 0.5, slow (l) Collabo., q * = 0.5, fast

Infection rate
Final epidemic size In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.

Theoretical results
We analytically derived an expression of the epidemic threshold for the SIS model on arbitrary underlying static networks combined with partnership models 1, 2, and 3. We assume that each network composed of partnership edges lasts for time h before switching to another partnership network and that the static network G(t) in each time window of length h is generated as an i.i.d.In fact, our partnership models 1, 2, and 3 imply that G(t) at different times are correlated with each other.Therefore, the i.i.d.assumption is for facilitating theoretical analysis.
We found that, to the first order of h, the epidemic threshold is the same for the three models.We show the detail in Appendix D. This result is apparently not consistent with the numerical results for the SIS model shown in Figs. 5 and 6, which support that the epidemic threshold decreases as the concurrency increases.In fact, the i.i.d.assumption made for our theoretical analysis corresponds to the limit of infinitely fast switching of edges.This is because a faster edge dynamics for a fixed value of h implies that G(t) and G(t ′ ), where t ′ ̸ = t, are less correlated with each other for given t and t ′ .Therefore, we carried out additional numerical simulations by making the values of a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 five times larger than those for the fast edge dynamics, i.e., 25 times larger than those for the slow edge dynamics.Then, we found that the epidemic threshold for the SIS model in the case of high concurrency (i.e., models 2 and 3) is much closer to the epidemic threshold in the case of low concurrency (i.e., model 1), as compared to when the edge dynamics are slower (i.e., Figs. 5 and 6).See Fig. 9 for numerical results.We conclude that the theoretical results in this section explain our numerical results when the edge dynamics are sufficiently faster than the epidemic dynamics.

Discussion
We investigated three models of edge dynamics.They allow us to generate cases with different amounts of concurrency while keeping the probability that each edge is active (i.e., q * ) and the structure of the aggregate network the same.We theoretically evaluated the concurrency of each model and compared it among the models.Then, we numerically observed using the stochastic SIS and SIR models that quasi-equilibrium fractions of infectious nodes and the final epidemic size, respectively, were larger for the edge dynamics models with the higher concurrency (i.e., models 2 or 3 compared to model 1) if the infection rate was near the epidemic threshold.In the SIS model, the epidemic threshold was also smaller when the concurrency was higher.These effects of the concurrency surrounding the epidemic threshold were larger when the edge dynamics is slower, i.e., with q * = 0.1 than with q * = 0.5.At the same time, we found that the effect of concurrency on epidemic spreading was modest up to our numerical efforts.It was particularly the case for the BA model and the collaboration network, which are more realistic than the ER random graph in terms of the heterogeneity in the node's degree.Furthermore, high concurrency suppressed epidemic spreading at large infection rates.
A previous study numerically showed that concurrency has a greater effect on increasing an epidemic potential in sparser networks [11].They measured the epidemic potential defined by the proportion of ordered pairs of nodes that are reachable in the sense that it is possible to travel from one node to the other node along a time-respecting path in the given temporal network composed of time-stamped edges.
Our results are consistent with theirs because the effect of concurrency on enhancing epidemic spreading is stronger when q * is smaller in all the cases investigated (see Figs. 5, 6, 7, and 8).The larger effect of concurrency on epidemic spreading for sparser networks may be because networks with large q * , which leads to a large number of edges, have more time-respecting paths from an infectious node to susceptible nodes transmitting infection even in the absence of concurrency [11].
Another previous study concluded that the impact of concurrency on the epidemic size is fairly limited except in an early stage of epidemic dynamics [27].Despite the dependence of epidemic spreading on the density of contacts, which we discussed in the last paragraph, our findings are at large consistent with theirs.In other words, we found that the presence of concurrency little affected or even decreased the quasiequilibrium fraction of infectious nodes in the SIS model and the final epidemic size in the SIR model in many cases, in particular for the BA model and collaboration networks.The reason why high concurrency decreases epidemic spreading at large infection rates is unclear.This phenomenon may be because, at high infection rates, infection can easily spread from nodes to nodes even without concurrency, such that many contact events on different edges, which would occur at close times in the case of high concurrency, are wasted.Note that our modeling framework guarantees that the total number of contact events on each edge is the same between low and high concurrency cases on average.Our analytical result showing that the epidemic threshold is independent of the concurrency (i.e., see section 4.2) also supports that the impact of concurrency on enhancing epidemic spreading may not be large.
In contrast, our previous study analytically showed that the concurrency enhances epidemic spreading in terms of the epidemic threshold value for the SIS model [28,39].The diversity of these results may be due to different strategies for modeling concurrency.For example, the authors of Ref. [27] designed their dynamic network models by ensuring that the number of partners that an individual has in a long term is the same across the individuals, rendering the aggregate network the complete graph.A different study that also carefully controlled the amount of interaction for each edge to be the same across comparisons made the same homogeneity assumption [29].In contrast, Refs.[28,39] and the present study have assumed that the number of partners that an individual has in a long term may be heterogeneously distributed.Examining similarities and differences between these existing models of concurrency, including the model proposed in the present study, is an outstanding question.This being said, an overall conclusion based on the present numerical results is that the concurrency only modestly or little influences extents of epidemic spreading in many cases.This result indicates that other factors such as static network structure [1][2][3][4] and burstiness [5][6][7]40] may be a larger contributor to epidemic dynamics than concurrency.Further comparing the impact of concurrency on epidemic spreading and that of other network factors warrants future work.
Concurrency is not only relevant to sexually transmitted infections, as most studies of concurrency have focused on [15,[20][21][22][23]41], or even to epidemic spreading in general.Other network dynamics such as opinion formation, synchronization, and information spreading and cascading on social networks may also be affected by concurrency.Our assumption that the partnership between any given pair of individuals can reoccur after it has disappeared is unrealistic in most cases of sexual partnerships.However, this assumption is considered to be natural for describing other types of dynamic social contacts such as face-to-face encounters and online communications, which underlie social dynamics apart from sexually transmitted infections.Investigating impacts of concurrency on these different social dynamics using the present models also warrants future work.i ̸ = j, are independent of each other, we obtain We also obtain (A.4)

Appendix B Proof of Proposition 6
For model 3, we obtain 2)a 3 , which is equivalent to 1 2 < q * ≤ 1.By applying Eq. (3.24) to Eqs. (3.39) and (3.43) and using Eq.(3.45), we obtain for each ℓ ∈ {1, . . ., M }.Then, we obtain where p(t) = [p 1 (t), . . ., p N (t)] ⊤ , p i (t) = E[x i (t)] is the probability that node v i is infectious at time t, ⊤ represents the transposition, q(t) = [q i1j1 (t), . . ., is the joint probability that v i is infectious and node v j is susceptible at time t, and the (N + M ) × (N + M ) matrix Ȃ is defined by where C + = max(C, 0) and C − = max(−C, 0) denote the positive and negative parts of the incidence matrix C, respectively.
Proof.We adapt Eq. (3.20) in Ref. [47] to the case in which the network is defined by G, the infection rate is independent of edges, and the recovery rate is independent of nodes.Specifically, by replacing B ′ , D, D ′ 1 , D ′ 2 in Eq. (3.20) in Ref. [47] by β Ξ, µI, µI, µI, respectively, we obtain where ϵ(t) is entry-wise nonnegative for every t ≥ 0. Equation (D.8) implies where matrix A is given by

D.2 A lower bound on the decay rate for temporal networks
Let G1 = (V, Ẽ1 ), . .., GL = (V, ẼL ) be directed and unweighted networks having the common node set V = {v 1 , . . ., v N }.Let {σ k } ∞ k=0 be independent and identically distributed random variables following a probability distribution on the set {1, . . ., L}. Variable σ k indexes the kth network to be used.Let h > 0 be arbitrary.For a real number x, let ⌊x⌋ denote the maximum integer that does not exceed x.Define the stochastic temporal network G by G(t) = Gσ ⌊t/h⌋ (D.13) for all t ≥ 0. Equation (D.13) implies that G(kh + τ ) = Gσ k for all k ≥ 0 and τ ∈ [0, h).In the temporal network models 1, 2, and 3 used in the main text, the states of individual edges or nodes independently repeat flipping in continuous time in a Poissonian manner, which induces network switching.Therefore, the duration for a single time-independent network, h, is not a constant.Furthermore, the time-independent networks before and after a single flipping of an edge's or node's state are not independent of each other.For these two reasons, the present G is not the same as models 1, 2, or 3.However, one can enforce the present G and any of models 1, 2, and 3 to have the same probability that each type of time-independent network, Gℓ , occurs by setting an appropriate probability distribution for { G1 , . . ., GL } for the present G.We assume such a probability distribution in the following text.We consider the stochastic SIS model taking place in G.
Definition 12.The decay rate of the SIS model over G is defined by where all nodes are assumed to be infected at t = 0.
Definition 12 states that N i=1 p i (t), which is equal to the expected number of infected nodes at time t, roughly decays exponentially in time in proportion to e −γt .Because the number of infected nodes always becomes zero in finite time, the SIS model on networks always has a positive decay rate.In fact, exact computation of the decay rate is computationally demanding because the decay rate is equal to the modulus of the largest real part of the non-zero eigenvalues of a 2 N ×2 N matrix representing the infinitesimal generator of the Markov chain [43].Therefore, bounds of the decay rate that only require computation of much smaller matrices are available [43,[47][48][49][50][51].Here we develop such a bound for temporal network G.
Let us label the set of time-aggregated edges, E = L ℓ=1 Ẽℓ , as E = {e 1 , . . ., e M }.Define the timeaggregated graph G = (V, E).Let C ∈ R N ×M denote the incidence matrix of G.Although the notations G, E, and C are common to those used in section D.1, this should not arise confusions in the following text.For each t ≥ 0, let E(t) denote the set of edges of the network at time t, i.e., G(t).In other words, E(t) = Ẽσ ⌊t/h⌋ .We also define the matrix H(t) ∈ R M ×M by H(t) ℓm = 1 if e ℓ ∈ E(t), j ℓ = i m , and j m ̸ = i ℓ , 0 otherwise, (D. 15) where ℓ, m ∈ {1, . . ., M }, and the diagonal matrix  Then, the decay rate of the SIS model over the temporal network G is greater than or equal to −h −1 log ρ(F).

D.3 Epidemic threshold for small h
In this section, we consider the case in which the switching interval, h, is sufficiently small.The first-order expansion of matrix F with respect to h yields Therefore, h −1 log ρ(F) = λ max = O(h).The combination of this asymptotic evaluation and Proposition 13 suggests that the epidemic threshold is the value of β/µ at which λ max = 0. We rewrite the epidemic threshold in terms of the spectral radius of the relevant matrices.As in the proof of Corollary 5 in Ref. [47], we rewrite the random matrix A(0) as   We keep E[Ξ(0)] constant across the different models to inspect the effect of different amounts of concurrency on epidemic spreading under the condition that the aggregate network is the same across the comparisons.Therefore, with this analysis, we do not find a difference in the epidemic threshold across our different models.

Figure 1 :
Figure 1: Schematic of part of temporal networks with different amounts of concurrency.We depict two edges sharing a node in each case.(a) Concurrent partnerships.(b) Non-concurrent partnerships.A shaded box represents a duration for which a partnership is present on the edge.The thick lines are examples of time-respecting paths transmitting infection from j to j ′ (shown in orange) or from j ′ to j (shown in magenta).This figure is inspired by Fig. 1 in Ref. [15] and Fig. 1 in Ref. [27].

Figure 2 :
Figure 2: Schematic illustration of models 2 and 3.In model 2, edge (1, 2) is active if and only if both node 1 and node 2 are in the h state.Otherwise, the edge is inactive.In model 3, edge (1, 2) is active if and only if either node 1 or node 2 is in the h state.

Figure 4 :
Figure 4: Distributions of the duration of the edge being inactive in model 2. The shaded bars represent numerically obtained distributions calculated on the basis of 5 × 10 5 samples.The solid lines represent the mixture of two exponential distributions, i.e., Eq. (3.49).The dashed lines represent the exponential distribution whose mean is the same as that for Eq.(3.49).(a) a = 2.0, b = 1.0.(b) a = 1.0, b = 2.0.

Figure 5 :
Figure 5: Comparison of the quasi-equilibrium fraction of infectious nodes in the SIS model between models 1 and 2. (a)-(d) ER random graph.(e)-(h) BA model.(i)-(l) Collaboration network.In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics."Collabo." is a shorthand for the collaboration (i.e., co-authorship) network.

Figure 6 :
Figure 6: Comparison of the quasi-equilibrium fraction of infectious nodes in the SIS model between models 1 and 3. (a)-(d) ER random graph.(e)-(h) BA model.(i)-(l) Collaboration network.In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.

Figure 7 :
Figure 7: Comparison of the final epidemic size in the SIR model between models 1 and 2. (a)-(d) ER random graph.(e)-(h) BA model.(i)-(l) Collaboration network.In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.The inset of (i) shows the magnification of the main panel because the final epidemic size is small in the entire range of the infection rate in this case.

Figure 8 :
Figure 8: Comparison of the final epidemic size in the SIR model between models 1 and 3. (a)-(d) ER random graph.(e)-(h) BA model.(i)-(l) Collaboration network.In panels (a), (e), and (i), we set q * = 0.1 and use the slow edge dynamics.In panels (b), (f), and (j), we set q * = 0.1 and use the fast edge dynamics.In panels (c), (g), and (k), we set q * = 0.5 and use the slow edge dynamics.In panels (d), (h), and (l), we set q * = 0.5 and use the fast edge dynamics.

Figure 9 :
Figure 9: Comparison of the quasi-stationary fraction of infectious nodes in the SIS model among models 1, 2, and 3 under faster edge dynamics.We use the edge dynamics that are five times faster than the fast edge dynamics used in Figs. 5 and 6. (a)-(d) ER random graph.(e)-(h) BA model.(i)-(l) Collaboration network.In panels (a), (e), and(i), we set q * = 0.1 and compare models 1 and 2. In panels (b), (f), and (j), we set q * = 0.5 and compare models 1 and 2. In panels (c), (g), and (k), we set q * = 0.1 and compare models 1 and 3.In panels (d), (h), and (l), we set q * = 0.5 and compare models 1 and 3. We remind that we cannot compare models 1, 2, and 3 in a single figure panel because we need to use different variants of model 1 in the comparison with models 2 and 3 to make the comparison fair.
2, edge (1, 2) is active if and only if both node 1 and node 2 are in the h state.Otherwise, the edge is inactive.In model 3, edge (1, 2) is active if and only if either node 1 or node 2 is in the h state.
(ξ 1 (t), . .., ξ M (t)), ∈ {1, . .., M }.Note that H(t) is similar to but different from the non-backtracking matrix of network G(t) because the definition of H(t) does not require e m ∈ E(t).The following proposition gives an upper bound on the decay rate of the SIS model over the temporal network G.
. From Corollary 11, we obtain d dt By combining Eq. (D.20) and the definition of the temporal network given by (D.13), we obtain Because the random variables {σ k } ∞ k=0 are independently and identically distributed, we take the expectation with respect to σ in Eq. (D.22) to obtain