Concurrency-induced transitions in epidemic dynamics on temporal networks

TL;DR

This paper investigates how the number of concurrent contacts in dynamic social networks influences epidemic thresholds, revealing that low concurrency suppresses while high concurrency enhances epidemic spread, with analytical and simulation validation.

Contribution

It provides a theoretical framework for understanding how concurrency affects epidemic thresholds on temporal networks, highlighting a phase transition in epidemic dynamics.

Findings

Low concurrency can suppress epidemics by raising the epidemic threshold.

High concurrency can promote epidemics by lowering the epidemic threshold.

The study identifies distinct phases of concurrency-induced transition in epidemic spread.

Abstract

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.

Equations

$${\rho }_1(a,t+\tau )=c_1\rho (a,t)+c_{2}m⟨\rho (t)⟩.$$

$${\rho }_2(a,a^{\prime },t+\tau )=c_3\rho (a,t)+c_4\rho (a^{\prime },t)+c_5(m-1)⟨\rho (t)⟩,$$

$$\rho (a,t+\tau )$$

$$\Theta (z,t+\tau )$$

$$\rho (a,t)=\sum _{n=1}^{\infty }{w}_n(t)a^{n-1}.$$

$$\mathbf{w}(t+\tau )=\mathbf{T}(\tau )\mathbf{w}(t).$$

$$\mathbf{T}=\left(\begin{array}{cccccc}{c}_3^{\prime }+⟨a⟩c_4^{\prime }+c_5^{\prime }& ⟨a^2⟩c_4^{\prime }+⟨a⟩c_5^{\prime }& ⟨a^3⟩c_4^{\prime }+⟨a^2⟩c_5^{\prime }& ⟨a^4⟩c_4^{\prime }+⟨a^3⟩c_5^{\prime }& ⟨a^5⟩c_4^{\prime }+⟨a^4⟩c_5^{\prime }& \cdots \\ c_1^{\prime }+c_2^{\prime }& ⟨a⟩c_2^{\prime }+c_3^{\prime }& ⟨a^2⟩c_2^{\prime }& ⟨a^3⟩c_2^{\prime }& ⟨a^4⟩c_2^{\prime }& \cdots \\ 0& c_1^{\prime }& c_3^{\prime }& 0& 0& \cdots \\ 0& 0& c_1^{\prime }& c_3^{\prime }& 0& \cdots \\ 0& 0& 0& c_1^{\prime }& c_3^{\prime }& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right).$$

$$f(\tau ,{\beta }_{\mathrm{c}})$$

$$2a{e}^{\frac{\left({\beta }_{\mathrm{c}}-1\right)\tau }{2}}\left[\cosh \left(\frac{{\kappa }_{\mathrm{c}}\tau }{2}\right)+\frac{1+3{\beta }_{\mathrm{c}}}{{\kappa }_{\mathrm{c}}}\sinh \left(\frac{-{\kappa }_{\mathrm{c}}\tau }{2}\right)\right]$$

$${\tau }_{\ast }=\ln \frac{1-(1+m)a}{1-(1+m{)}^{2}a}.$$

$$m_{\mathrm{c}}=\frac{3}{1-4a}.$$

$$p=N⟨a⟩\left[1-{\left(1-\frac{m+1}{N-1}\right)}^m\right]\approx m(m+1)⟨a⟩.$$

$$M_{\{x,y,z\},\{x,y,z\}}=-\sum _{\{x^{\prime },y^{\prime },z^{\prime }\}\ne \{x,y,z\}}{M}_{\{x^{\prime },y^{\prime },z^{\prime }\},\{x,y,z\}}.$$

$$M_{\{S,y,z\},\{I,y,z\}}=$$

$$M_{\{x,S,z\},\{x,I,z\}}=$$

$$M_{\{x,y,z-1\},\{x,y,z\}}=$$

$$M_{\{I,S,z\},\{S,S,z\}}=$$

$$M_{\{I,I,z\},\{S,I,z\}}=$$

$$M_{\{I,I,z\},\{I,S,z\}}=$$

$$M_{\{I,y,z+1\},\{I,y,z\}}=$$

$$\dot{\mathbf{p}}(t)=\mathbf{M}\mathbf{p}(t),$$

$$\mathbf{p}(t)=\exp (\mathbf{M}t)\mathbf{p}(0).$$

$$\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ c_5\end{array}\right)=\left(\begin{array}{c}{\sum }_{y,z}{\left[\exp (\mathbf{M}\tau )\right]}_{\{I,y,z\},\{I,S,0\}}\\ {\sum }_{y,z}{\left[\exp (\mathbf{M}\tau )\right]}_{\{I,y,z\},\{S,I,0\}}\\ {\sum }_{x,z}{\left[\exp (\mathbf{M}\tau )\right]}_{\{x,I,z\},\{S,I,0\}}\\ {\sum }_{x,z}{\left[\exp (\mathbf{M}\tau )\right]}_{\{x,I,z\},\{I,S,0\}}\\ {\sum }_{x,z}{\left[\exp (\mathbf{M}\tau )\right]}_{\{x,I,z\},\{S,S,1\}}\end{array}\right).$$

$$c_1=c_3=$$

$$\frac{e^{-\tau }}{2}\left[e^{-\beta \tau }+e^{-\frac{1+\beta }{2}\tau }\left(\cosh \frac{\kappa \tau }{2}+\frac{1+3\beta }{\kappa }\sinh \frac{\kappa \tau }{2}\right)\right],$$

$$c_2=c_4=$$

$$\frac{e^{-\tau }}{2}\left[-e^{-\beta \tau }+e^{-\frac{1+\beta }{2}\tau }\left(\cosh \frac{\kappa \tau }{2}+\frac{1+3\beta }{\kappa }\sinh \frac{\kappa \tau }{2}\right)\right],$$

$$p_{\{x,y,z\}}\approx P(x)P(y)P(z),$$

$${\mathbf{p}}^{\mathrm{MF}}=\left(\begin{array}{c}P(x=I)\\ P(y=I)\\ \frac{⟨z⟩}{m-1}\end{array}\right)=\left(\begin{array}{c}{\sum }_{y,z}{p}_{\{I,y,z\}}\\ {\sum }_{x,z}{p}_{\{x,I,z\}}\\ \frac{1}{m-1}{\sum }_{x,y,z}z{p}_{\{x,y,z\}}\end{array}\right).$$

$${\dot{\mathbf{p}}}^{\mathrm{MF}}$$

$${\dot{\mathbf{p}}}^{\mathrm{MF}}\approx {\mathbf{M}}^{\mathrm{MF}}{\mathbf{p}}^{\mathrm{MF}}$$

$${\mathbf{M}}^{\mathrm{MF}}=\left(\begin{array}{ccc}-1& \beta & (m-1)\beta \\ \beta & -1& 0\\ \beta & 0& -1\end{array}\right).$$

$$\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ c_5\end{array}\right)$$

$$p_3^{\mathrm{MF}}(\tau )\le \beta (1-e^{-\tau })+p_3^{\mathrm{MF}}(0)e^{-\tau }$$

Full text

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 Holme and Saramäki (2012); Holme (2015); Masuda and Lambiotte (2016). 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 Keeling and Eames (2005); Barrat et al. (2008); Pastor-Satorras et al. (2015); Porter and Gleeson (2016) 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 Bansal et al. (2010); Holme and Saramäki (2012); Masuda and Holme (2013); Holme (2015); Masuda and Lambiotte (2016).

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 Morris and Kretzschmar (1995); Kretzschmar and Morris (1996); Miller and Slim (2016). 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 Watts and May (1992); Eames and Keeling (2004); Doherty et al. (2006); Gurski and Hoffman (2016). 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 Morris and Kretzschmar (1995); Kretzschmar and Morris (1996); Bauch and Rand (2000); Leung and Kretzschmar (2015). In the present study, we use the analytically tractable activity-driven model of temporal networks Perra et al. (2012); Starnini and Pastor-Satorras (2013); Ribeiro et al. (2013); Liu et al. (2014); Zino et al. (2016) 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  Holme and Saramäki (2012); Holme (2015); Masuda and Lambiotte (2016); Masuda and Holme (2013); Zino et al. (2016); Vazquez et al. (2007); Karsai et al. (2011); Miritello et al. (2011); Stehlé et al. (2011), 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 Perra et al. (2012); Starnini and Pastor-Satorras (2013); Ribeiro et al. (2013); Liu et al. (2014); Zino et al. (2016). The number of nodes is denoted by $N$. Each node $i$ $(1\leq i\leq N)$ is assigned an activity potential $a_{i}$, drawn from a probability density $F(a)$ $(0<a\leq 1)$. Activity potential $a_{i}$ is the probability with which node $i$ is activated in a window of constant duration $\tau$. 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 $\tau$, 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 Liberzon (2003); Masuda et al. (2013); Hasler et al. (2013a); Speidel et al. (2016). A large $\tau$ implies that network dynamics are slow compared to epidemic dynamics. In the limit of $\tau\to 0$, the network blinks infinitesimally fast, enabling the dynamical process to be approximated on a time-averaged static network, as in Hasler et al. (2013a).

For the SIS dynamics, each node takes either the susceptible or infected state. At any time, each susceptible node contracts infection at rate $\beta$ per infected neighboring node. Each infected node recovers at rate $\mu$ irrespectively of the neighbors’ states. Changing $\tau$ to $c\tau$ $(c>0)$ is equivalent to changing $\beta$ and $\mu$ to $\beta/c$ and $\mu/c$, respectively, while leaving $\tau$ unchanged. Therefore, we set $\mu=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 $\approx m\sum_{j\neq i}a_{j}(m+1)/N\propto m^{2}\langle a\rangle$, where $\langle a\rangle\equiv\int{\rm d}aF(a)a$ is the mean activity potential, we impose $m^{2}\langle a\rangle\ll 1$. (However, our method works better than the so-called individual-based approximation even when $m^{2}\langle a\rangle=0.5$, as shown in the Supplemental Material.) We denote by $\rho(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 $\langle\rho(t)\rangle\equiv\int{\rm d}aF(a)\rho(a,t)$. Let $c_{1}$ be the probability with which the hub in an isolated star graph is infected at time $t+\tau$, 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+\tau$ when only a single leaf node is infected at $t$. The probability that a hub with activity potential $a$ is infected after the duration $\tau$ of the star graph, denoted by $\rho_{1}$, is given by

[TABLE]

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 $\rho(a,t),\langle\rho(t)\rangle\ll 1$]. The probability that a leaf with activity potential $a$ that has a hub neighbor with activity potential $a^{\prime}$ is infected after time $\tau$ is analogously given by

[TABLE]

where $c_{3}$, $c_{4}$, and $c_{5}$ are the probabilities with which a leaf node with activity potential $a$ is infected after duration $\tau$ 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\leq i\leq 5)$ in the Supplemental Material. The probability that an isolated node with activity potential $a$ is infected after time $\tau$ is given by $e^{-\tau}\rho(a,t)$. By combining these contributions, we obtain

[TABLE]

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

[TABLE]

where $c_{1}^{\prime}\equiv c_{1}-e^{-\tau}$, $c_{2}^{\prime}\equiv mc_{2}$, $c_{3}^{\prime}\equiv e^{-\tau}+m\langle a\rangle(c_{3}-e^{-\tau})$, $c_{4}^{\prime}\equiv mc_{4}$, $c_{5}^{\prime}\equiv m(m-1)\langle a\rangle c_{5}$, $g(z)\equiv\int{\rm d}aF(a)z^{a}$ is the probability generating function of $a$, $\Theta(z,t)\equiv\int{\rm d}aF(a)\rho(a,t)z^{a}$, and throughout the paper the superscript $(n)$ represents the $n$th 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 $\Theta(z,t)$ and its derivative Silk et al. (2014).

We expand $\rho(a,t)$ as a Maclaurin series as follows:

[TABLE]

Using this polynomial basis representation (the convergence is proven in the Supplemental Material), we can consider the differentiations in Eq. (4) (i.e., $\Theta^{(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 $\bm{w}(t=0)\equiv(w_{1}(0),w_{2}(0),\ldots)^{\top}=\left(p_{0},0,0,\ldots\right)^{\top}$. Epidemic dynamics near the epidemic threshold obey linear dynamics given by

[TABLE]

By substituting $\Theta(z,t)=\sum_{n=1}^{\infty}w_{n}(t)g^{(n-1)}(z)$ and $g^{(n-1)}(1)=\langle a^{n-1}\rangle$ in Eq. (4), we obtain

[TABLE]

A positive prevalence $\left<\rho(t)\right>$ (i.e., a positive fraction of infected nodes in the equilibrium state) occurs only if the largest eigenvalue of $\bm{T}(\tau)$ 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 $\beta_{\rm c}$:

[TABLE]

where $S(q)\equiv\sum_{n=0}^{\infty}(\langle a^{n+2}\rangle/\langle a\rangle^{n+2})q^{n}=(1/\langle a\rangle^{2})\left\{\left\langle(a^{2})/[1-(a/\langle a\rangle)q]\right\rangle\right\}$, $q\equiv\langle a\rangle c_{1}^{\prime}/(1-c_{3}^{\prime})$, $r\equiv\langle a\rangle c_{2}^{\prime}/(1-c_{3}^{\prime})$, $s\equiv\langle a\rangle c_{4}^{\prime}/(1-c_{3}^{\prime})$, and $u\equiv c_{5}^{\prime}/(1-c_{3}^{\prime})$ (see Supplemental Material for the derivation). Note that $f$ is a function of $\beta$ ($=\beta_{\rm c}$) through $q$, $r$, $s$, and $u$, which are functions of $\beta$. In general, we obtain $\beta_{\rm c}$ by numerically solving Eq. (8), but some special cases can be determined analytically.

In the limit $\tau\to 0$, Eq. (8) gives $\beta_{\rm c}=\left[m\left(\langle a\rangle+\sqrt{\langle a^{2}\rangle}\right)\right]^{-1}$, which coincides with the epidemic threshold for the activity-driven model derived in the previous studies Perra et al. (2012); Liu et al. (2014). In fact, this $\beta_{\rm c}$ value is the epidemic threshold for the aggregate (and hence static) network, whose adjacency matrix is given by $A_{ij}^{*}\approx m(a_{i}+a_{j})/N$  Speidel et al. (2016); Masuda and Lambiotte (2016), as demonstrated in Fig. S1.

For general $\tau$, if all nodes have the same activity potential $a$, and if $m=1$, we obtain $\beta_{\rm c}$ as the solution of the following implicit equation:

[TABLE]

where $\kappa_{\rm c}=\sqrt{\beta_{\rm c}^{2}+6\beta_{\rm 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 $\tau$ and $\beta$ values shown in different colors. Equations (8) and (9) describe the numerical results fairly well. When $m=1$, the epidemic threshold increases with $\tau$ and diverges at $\tau\approx 0.1$ [Fig. 2(a)]. Furthermore, slower network dynamics (i.e., larger values of $\tau$) reduce the prevalence for all values of $\beta$. In contrast, when $m=10$, the epidemic threshold decreases and then increases as $\tau$ increases [Fig. 2(b)]. The network dynamics (i.e., finite $\tau$) 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 $\beta_{\rm c}$ by the individual-based approximation (Speidel et al. (2016), see Supplemental Material for the derivation), which may be justified when $m\gg 1$, is consistent with the numerical results and our theoretical results only at small $\tau$ [a dashed line in Fig. 2(b)]. Qualitatively similar results are found, when the activity potential $a$ is power-law distributed [Figs. 2(c) and 2(d)].

To illuminate the qualitatively different behaviors of the epidemic threshold as $\tau$ 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 $\tau=\tau_{*}$. In the limit $\beta\to\infty$, infection starting from a single infected node in a star graph immediately spreads to the entire star graph, leading to $c_{i}\to 1$ $(1\leq i\leq 5)$. By substituting $c_{i}\to 1$ in Eq. (8), we obtain $f(\tau_{*},\beta_{\rm c}\to\infty)=0$, where

[TABLE]

When $\tau>\tau_{*}$, 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 Keeling and Ross (2008); Simon et al. (2011); Hindes et al. (2016). Second, although $\beta_{\rm c}$ eventually diverges as $\tau$ increases, there may exist $\tau_{\rm c}$ such that $\beta_{\rm c}$ at $\tau<\tau_{\rm c}$ is smaller than the $\beta_{\rm c}$ value at $\tau=0$. Motivated by the comparison between the behavior of $\beta_{\rm c}$ at $m=1$ and $m=10$ (Fig. 2), we postulate that $\tau_{\rm c}$ ($>0$) exists only for $m>m_{\rm c}$. Then, we obtain ${\rm d}\beta_{\rm c}/{\rm d}\tau=0$ at $(\tau,m)=(0,m_{\rm c})$. The derivative of Eq. (8) gives $\partial f/\partial\tau+(\partial f/\partial\beta_{\rm c})({\rm d}\beta_{\rm c}/{\rm d}\tau)=0$. Because ${\rm d}\beta_{\rm c}/{\rm d}\tau=0$ at $(\tau,m)=(0,m_{\rm c})$, we obtain $\partial f/\partial\tau=0$, which leads to

[TABLE]

When $m<m_{\rm c}$, network dynamics (i.e., finite $\tau$) always reduce the prevalence for any $\tau$ [Figs. 2(a) and 2(c)]. When $m>m_{\rm c}$, a small $\tau$ raises the prevalence as compared to $\tau=0$ (i.e., static network) but a larger $\tau$ reduces the prevalence [Figs. 2(b) and 2(d)].

The phase diagram based on Eqs. (S49) and (11) is shown in Fig. 3(a). The $\beta_{\rm 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 $\beta_{\rm c}$ has the same value for all $m$ at $\tau=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 event as some temporal network models do Tantipathananandh et al. (2007); Stehlé et al. (2010); Zhao et al. (2011).

Discussion: Our analytical method shows that the presence of network dynamics boosts the prevalence (and decreases the epidemic threshold $\beta_{\rm c}$) when the concurrency $m$ is large and suppresses the prevalence (and increases $\beta_{\rm c}$) when $m$ is small, for a range of values of the network dynamic time scale $\tau$. This result lends theoretical support to previous claims that concurrency boosts epidemic spreading Watts and May (1992); Morris and Kretzschmar (1995); Kretzschmar and Morris (1996); Bauch and Rand (2000); Eames and Keeling (2004); Doherty et al. (2006); Eaton et al. (2011); Perra et al. (2012); Leung and Kretzschmar (2015); Gurski and Hoffman (2016). 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 $\tau$, regardless of the level of concurrency. As shown in Figs. 3 and S6, the transition to the “die out” phase occurs at values of $\tau$ 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 formulation-dissolution models Kretzschmar and Morris (1996); Bauch and Rand (2000); Eames and Keeling (2004); Leung and Kretzschmar (2015); Gurski and Hoffman (2016) and individual-based approximations Valdano et al. (2015); Speidel et al. (2016); Rocha and Masuda (2016). 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 thank Leo Speidel for discussion. We thank the SocioPatterns collaboration (http:// www.sociopatterns.org) for providing the data set. T.O. acknowledges the support provided through JSPS Research Fellowship for Young Scientists. J.G. acknowledges the support provided through Science Foundation Ireland (Grants No. 15/SPP/E3125 and No. 11/PI/1026). N.M. acknowledges the support provided through JST, CREST, and JST, ERATO, Kawarabayashi Large Graph Project.

I Prevalence on the aggregate network

II When the low-activity assumption is violated

Here we consider the situation in which the low-activity assumption $m^{2}\langle a\rangle\ll 1$ is violated. When $m\ll N$, the expected number of star graphs that a star graph overlaps with is given by

[TABLE]

If $p\ll 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\geq 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. If $p\geq 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)].

III Derivation of $c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, and $c_{5}$

We consider SIS dynamics on a star graph with $m$ leaves and derive $c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, and $c_{5}$. Let us denote the state of the star graph by $\{x,y,z\}~{}(x,y\in\{S,I\},0\leq z\leq m-1)$, where $x$ and $y$ are the states of the hub and a specific leaf node, respectively, and $z$ is the number of infected nodes in the other $m-1$ leaf nodes. Although a general network with $m+1$ nodes allows $2^{m+1}$ states, using this notation, we can describe SIS dynamics on a star graph by a continuous-time Markov process with $4m$ states Simon et al. (2011).

We denote the transition rate matrix of the Markov process by $\bm{M}$. Its element $M_{\{x^{\prime},y^{\prime},z^{\prime}\},\{x,y,z\}}$ is equal to the rate of transition from $\{x,y,z\}$ to $\{x^{\prime},y^{\prime},z^{\prime}\}$. The diagonal elements are given by

[TABLE]

The rates of the recovery events are given by

[TABLE]

The rates of the infection events are given by

[TABLE]

The other elements of $\bm{M}$ are equal to [math]. Let $p_{\{x,y,z\}}(t)$ be the probability for a star graph to be in state $\{x,y,z\}$ at time $t$. Because

[TABLE]

where $\bm{p}(t)$ is the $4m$-dimensional column vector whose elements are $p_{\{x,y,z\}}(t)$, we obtain

[TABLE]

Note that $c_{1}$ and $c_{2}$ are the probabilities with which $x=I$ at time $\tau$, 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 $\tau$, when the initial state is $\{S,I,0\}$, $\{I,S,0\}$, and $\{S,S,1\}$, respectively. Therefore, we obtain

[TABLE]

When $m=1$, Eq. (S12) yields

[TABLE]

where $\kappa=\sqrt{\beta^{2}+6\beta+1}$, and $c_{5}$ is not defined.

When $m\gg 1$, we can apply an individual-based approximation Speidel et al. (2016); Valdano et al. (2015); Pastor-Satorras et al. (2015). We assume that the state of each node is statistically independent of each other, i.e.,

[TABLE]

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^{\rm MF}_{1}$ and $p^{\rm MF}_{2}$, respectively, and $z$ obeys a binomial distribution with parameters $m-1$ and $p^{\rm MF}_{3}$, where $\bm{p}^{\rm MF}\equiv(p_{1}^{\rm MF},p_{2}^{\rm MF},p_{3}^{\rm MF})^{\top}$ is given by

[TABLE]

By substituting Eq. (S10) in the time derivative of Eq. (S16), we obtain

[TABLE]

If $p^{\rm MF}_{3}\ll 1$, $\bm{p}^{\rm MF}$ obeys linear dynamics given by

[TABLE]

where

[TABLE]

In a similar fashion to the derivation of Eq. (S12), we obtain

[TABLE]

We estimate the extent to which Eq. (S20) is valid as follows. First, we need $m\gg 1$, because the initial condition $p^{\rm MF}_{3}=1/(m-1)$ should satisfy $p^{\rm MF}_{3}\ll 1$. Second, $p^{\rm MF}_{3}$ must satisfy

[TABLE]

because $p^{\rm MF}_{1}\leq 1$ in Eq. (S17). To satisfy $p^{\rm MF}_{3}\ll 1$, we need $\tau<1/\beta$. This condition remains unchanged by re-scaling $(\tau,\beta)$ to $(c\tau,\beta/c)$. These two conditions are sufficient for this approximation to be valid. If $m\gg 1$ is violated, the individual-based approximation significantly underestimates the epidemic threshold for any finite $\tau$ because it ignores the effect of stochastic dying-out. If $\tau<1/\beta$ 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 $\tau\to\infty$, whereas analytical [Eq. (8)] and numerical (Fig. 2) results diverge at a finite $\tau$.

IV Derivation of Eq. (8)

At the epidemic threshold, the largest eigenvalue of $\bm{T}$ is equal to unity. Let $\bm{v}=(v_{1},v_{2},\ldots)^{\top}$ be the corresponding eigenvector of $\bm{T}$. We normalize $\bm{v}$ such that $\sum_{j=1}^{\infty}v_{j}=1$. By substituting Eq. (7) in $\bm{v}=\bm{T}\bm{v}$, we obtain the system of equations

[TABLE]

Equation (S24) gives

[TABLE]

where

[TABLE]

By combining Eqs. (S23) and (S25), we obtain

[TABLE]

where

[TABLE]

Because $\bm{v}$ is normalized, we obtain

[TABLE]

Equation (S22) leads to

[TABLE]

where,

[TABLE]

By substituting Eq. (S30) in Eq. (S31), we obtain

[TABLE]

which is Eq. (8) in the main text. If all nodes have the same activity potential $a$, Eq. (S34) is reduced to

[TABLE]

V Convergence of the Maclaurin series

We derive the condition under which the Maclaurin series in Eq. (5) converges for any $t$ when $\beta\leq\beta_{\rm c}$. First, at $t=0$, the series converges because $\bm{w}(0)=(p_{0},0,0,...)^{\top}$.

Second, consider a finite $t$. It should be noted that the series is only defined at $t$ that is a multiple of $\tau$. Because $T_{ij}=0$ $(i\geq j+2)$ in Eq. (7), we obtain

[TABLE]

Therefore, the series converges.

Third, we consider the limit $t\to\infty$. If $\beta<\beta_{\rm c}$, because

[TABLE]

we obtain

[TABLE]

Therefore, the series converges. For $\beta=\beta_{\rm c}$, we consider the convergence of the series when

[TABLE]

where $\bm{v}$ is the eigenvector of $\bm{T}$ given by Eq. (S30), and $b$ is a constant. Because Eq. (S30) yields

[TABLE]

the radius of convergence is equal to $\langle a\rangle/q$. To ensure convergence, we require that

[TABLE]

Because $c_{i}$ $(1\leq i\leq 5)$ are probabilities, we obtain

[TABLE]

By substituting Eqs. (S42) and (S43) in the definitions of $c_{1}^{\prime}$ and $c_{2}^{\prime}$, we obtain

[TABLE]

By substituting Eqs. (S44) and (S45) in Eq. (S26), we obtain

[TABLE]

Inequalities (S42)–(S46) hold with equality in the limit $\beta\to\infty$. Hence, a sufficient condition for convergence is given by

[TABLE]

Equation (S47) holds true in practical situations because the assumption $m^{2}\langle a\rangle\ll 1$ guarantees that $m\langle a\rangle\ll 1$ and $a_{i}$ is a probability.

VI Epidemic threshold under the individual-based approximation

When $m\gg 1$, the epidemic threshold can be obtained by the individual-based approximation Speidel et al. (2016); Valdano et al. (2015); Pastor-Satorras et al. (2015). We assume that all nodes have the same activity potential $a$. By substituting Eq. (S20) in Eq. (S35), we obtain

[TABLE]

Equation (S48) agrees with the value derived in Speidel et al. (2016). Note that this approximation is valid only for small $\tau$ ($\tau<1/\beta_{\rm c}$).

VII Derivation of $\tau_{*}$ for general activity distributions

In the limit $\beta\to\infty$, we obtain $c_{i}\to 1$ ($1\leq i\leq 5$). For general activity distributions, $f(\tau_{*},\beta_{\rm c}\to\infty)=0$ leads to

[TABLE]

where

[TABLE]

VIII Derivation of $m_{\rm c}$ for general activity distributions

At $m=m_{\rm c}$, an infinitesimal increase in $\tau$ from [math] to $\Delta\tau$ does not change the $\beta_{\rm c}$ value. For general activity distributions, by setting $\partial f/\partial\tau=0$ for $f$ given by Eq. (S34), we obtain

[TABLE]

IX Activity-driven model with a reinforcement process

We carried out numerical simulations for an extended activity-driven model in which link dynamics are driven by a reinforcement process Karsai et al. (2015). The original activity-driven model is memoryless Perra et al. (2012). 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.

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\leq i\leq N)$ in (a) and (c) and $a_{i}=0.005$ $(1\leq i\leq N)$ in (b) and (d). Figure S3 indicates that the extended model only slightly changes the epidemic threshold.

X Stochastic $m$

We consider the case in which the strength of concurrency, $m$, is not constant. To analyze this case, we change the definitions of $c_{1}^{\prime}$, $c_{2}^{\prime}$, $c_{3}^{\prime}$, $c_{4}^{\prime}$, and $c_{5}^{\prime}$ to

[TABLE]

where ${\mathrm{E}}[\cdot]$ is the expectation with respect to the distribution of $m$. The mean degree is given by $\langle k\rangle=2a{\mathrm{E}}[m]$. Using Eqs. (S53)–(S57) instead of $c_{i}^{\prime}$ $(1\leq i\leq 5)$, we derived the epidemic threshold in the same manner as the derivation of Eq. (8). The phase diagrams of the epidemic threshold when $m$ obeys a truncated Poisson distribution and a power-law distribution are shown in Figs. S4(a) and S4(b), respectively. We obtain $\beta_{\rm c}=1/\langle k\rangle$ at $\tau=0$. We set the activity potential of all nodes $a=\langle k\rangle/(2{\mathrm{E}}[m])$ such that the epidemic threshold is the same for all ${\mathrm{E}}[m]$ at $\tau=0$. We numerically calculated $m_{\rm c}$ at which $\tau_{\rm c}=0$. For the power-law distribution of $m$, we cannot make ${\mathrm{E}}[m]$ smaller than $m_{\rm c}$ because the distribution does not have a probability mass at $m=0$ by definition. However, the phase diagrams in the case of both the truncated Poisson and power-law distributions of $m$ are qualitatively similar to the case of constant $m$.

To gain analytical insights, we calculated the phase diagrams when $m$ is equal to $m_{1}$ and $m_{2}$ with probabilities $\tilde{p}$ and $1-\tilde{p}$, respectively. We varied $\tilde{p}$ between [math] and $1$. Here again, we set the activity potential of all nodes $a=\langle k\rangle/(2{\mathrm{E}}[m])$ such that the epidemic threshold is the same for all ${\mathrm{E}}[m]$ at $\tau=0$. The phase diagram [Fig. S4(c)] is again qualitatively similar to that found in the case of constant $m$.

XI Heterogeneous activity distributions

We analyzed the phase diagram for different distributions of activity potentials to confirm the robustness of the results shown in the main text. We consider an exponential distribution and a power-law distribution with exponent $2.5$. We numerically calculate the epidemic threshold by solving Eq. (8) and derive $\tau_{*}$ and $m_{\rm c}$ from Eqs. (S49) and (S52), respectively. The phase diagrams for the exponential and power-law distributions are shown in Figs. S5(a) and S5(b), respectively. These results are qualitatively similar to those found when all nodes have the same activity potential value.

XII Temporal networks composed of cliques

We consider the case in which an activated node creates a clique (a fully-connected subgraph) with $m$ randomly chosen nodes instead of a star graph. This situation models a group conversation among $m+1$ people. We only consider the case in which all nodes have the same activity potential $a$. The mean degree for a network in a single time window is given by $\langle k\rangle=m(m+1)a$. The aggregate network is the complete graph. We impose $m^{2}a\ll 1$ so that cliques in the same time window do not overlap.

As in the case of the activity-driven model, we denote the state of a clique by $\{x,y,z\}~{}(x,y\in\{S,I\},0\leq z\leq m-1)$, where $x$ and $y$ are the states of the activated node and another specific node, respectively, and $z$ is the number of infected nodes in the other $m-1$ nodes. The transition rate matrix of the SIS dynamics on this temporal network model is given as follows. The rates of the recovery events are given by Eqs. (S3), (S4), and (S5). The rates of the infection events are given by

[TABLE]

We obtain $c_{i}$ $(1\leq i\leq 5)$ from $\bm{M}$ in the same fashion as in the case of the activity-driven model. Because of the symmetry inherent in a clique, we obtain $c_{1}=c_{3}$ and $c_{2}=c_{4}=c_{5}$. Therefore, Eq. (S35) is reduced to

[TABLE]

Calculations similar to the case of the activity-driven model lead to

[TABLE]

The phase diagram shown in Fig. S6 is qualitatively the same as those for the activity-driven model (Fig. 3). Note that, in Fig. S6, we selected the activity potential value $a$ to force $\beta_{\rm c}$ to be independent of $m$ at $\tau=0$, i.e.,

[TABLE]

Although Eq. (S67) coincides with the expression of $\tau_{*}$ for the activity-driven model [Eq. (10)], $\tau_{*}$ as a function of $m$ is different between the activity-driven model [a solid line in Fig. 3(a)] and the present clique network model (a solid line in Fig. S6). This is because the values of $a$ are different between the two cases when $m\geq 2$.

XIII Empirical activity distributions

The epidemic threshold and prevalence when $F(a)$ is constructed from empirical contact data at a workplace, obtained from the SocioPatterns project Genois et al. (2015), are shown in Figs. S7(a) and S7(b) for $m=1$ and $m=10$, respectively. The results for $F(a)$ constructed from email communication data at a research institution, obtained from the Stanford Network Analysis Platform Paranjape et al. (2015), are shown in Figs. S7(c) and S7(d) for $m=1$ and $m=10$, respectively. These results are qualitatively similar to those shown in Fig. 2.