Epidemic spreading through direct and indirect interactions Niloy Ganguly,1,* Tyll Krueger,2,† Animesh Mukherjee,1,‡ and Sudipta Saha1,§ 1

Department of Computer Science and Engineering, Indian Institute of Technology, Kharagpur 721302, India 2 Department of Computer Science and Engineering, Wroclaw University of Technology, 50-370 Wroclaw, Poland (Received 2 June 2014; revised manuscript received 27 August 2014; published 17 September 2014) In this paper we study the susceptible-infected-susceptible epidemic dynamics, considering a specialized setting where popular places (termed passive entities) are visited by agents (termed active entities). We consider two types of spreading dynamics: direct spreading, where the active entities infect each other while visiting the passive entities, and indirect spreading, where the passive entities act as carriers and the infection is spread via them. We investigate in particular the effect of selection strategy, i.e., the way passive entities are chosen, in the spread of epidemics. We introduce a mathematical framework to study the effect of an arbitrary selection strategy and derive formulas for prevalence, extinction probabilities, and epidemic thresholds for both indirect and direct spreading. We also obtain a very simple relationship between the extinction probability and the prevalence. We pay special attention to preferential selection and derive exact formulas. The analysis reveals that an increase in the diversity in the selection process lowers the epidemic thresholds. Comparing the direct and indirect spreading, we identify regions in the parameter space where the prevalence of the indirect spreading is higher than the direct one. DOI: 10.1103/PhysRevE.90.032808

PACS number(s): 89.90.+n, 89.65.−s, 89.75.−k

I. INTRODUCTION

Spreading of contagions, e.g., information or messages in social or technological communications and pathogens in disease outbreak, has been mostly studied in the framework of contact processes taking place directly among the actors of the system [1–4] or through the carriers of pathogens, e.g., the mosquitoes in the case of malaria [5,6]. However, it is important to note that usually most of the contacts among the actors of a system take place in popular and frequently visited places [5,6]. The entity place can be interpreted in a wider context, e.g., commonly used objects, crowded spots, and popular social groups in online social network. In the context of sexually transmitted diseases, “place” can be translated to the sex workers who can act as a carrier of disease. The strategies adopted by people to select these entities one after another may crucially influence the spreading dynamics. However, the existing studies on spreading dynamics do not systematically capture the effect of these actor-place interactions and the selection strategies in general of an actor for a place, which is the prime focus of the present paper. The actors in the system, therefore, can be categorized into two classes, active and passive, where the passive entities get selected by the active entities. One possible interpretation of this scenario would be the mobile entities such as humans as the active entities and various immobile or stationary entities such as popular places or social groups as the passive ones. The dynamical interaction process between the active and passive entities can be understood through an evolving bipartite network structure where one partition corresponds to the set of active entities and the other one corresponds to the passive entities. The selection of a passive entity by an

*

[email protected] [email protected] ‡ [email protected] § Corresponding author: [email protected] †

1539-3755/2014/90(3)/032808(12)

active entity can be represented by an edge between the active and the passive entity. Thus, the contact network structure is dynamical and changes at every time step. Figure 1 describes this basic bipartite network representation. In many cases the passive entities are capable of temporary preservation of the contagion. For example, the throwboxes installed in popular places in a delay tolerant network [7,8] can store the messages temporarily. Similarly, the pathogens can survive in an infectious status for a limited period of time on the commonly used objects [9–12]. Note that this temporary storing of the contagions allows the spreading to take place in an indirect fashion. We term this kind of spreading through such indirect communications indirect spreading. On the other hand, as usual, spreading can also take place through direct communications among the actors when they come into contact in popular places, e.g., people catching air-borne diseases while visiting crowded spots and active phones exchanging photograph using Bluetooth when in close contact. We call this direct spreading. In this paper, both the direct and indirect spreading is studied primarily with a focus on understanding the role of the selection strategies adopted by the people to select the common passive entities. Based on selection strategy, there are two major categories of models. One in its basic form is a random walk [13,14], where the next destination is selected according to a probabilistic rule among the existing neighboring places. The other is based on nonuniform jump processes where places are chosen preferentially, e.g., according to some popularity values [15,16]; similarly the choice people make while joining social groups can also be modeled as a preferential selection [17]. Here the popularity values themselves are usually dynamical variables that depend on the evolutionary history of the corresponding bipartite network. To understand the effect of these selection strategies in a generalized fashion, we assume that the passive entities are preassigned with attractiveness values following some specific probability distribution. These attractiveness values characterize their popularity. This seems to be a reasonable

032808-1

©2014 American Physical Society

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

PHYSICAL REVIEW E 90, 032808 (2014)

FIG. 1. (Color online) Schematic diagram showing the spreading process running through indirect communication among active and passive entities. The status of the system at the beginning of the (a) first, (b) second, (c) third, and (d) fourth time steps in the process is shown. There are four passive entities P1 , . . . ,P4 and six active entities A1 , . . . ,A6 in the system. Here A3 is the initiator of the process. Each of the active entities comes into contact with one of the passive entities at every time step. Solid arrows denote the selection of passive entities by the active entities. The open and closed circles (squares) denote active (passive) entities without or with contagion, in susceptible or infected status, respectively, at the specific annotated time step. The dotted arrows denote the direction of transfer of the contagion. At the end of each time step, all the passive and active entities are probabilistically refreshed (by a certain value of α and β), which results, at the end of t = 3, in A1 , A2 , and P2 losing their contagion and going to the susceptible status. The changes due to an operation performed in a particular time step are reflected in the next time step.

general assumption since in most models with a priori timeor history-dependent popularity the values converge to some limiting distribution. We discuss two specific cases in detail where the selection strategies are uniform and preferential. The dynamics we consider is a special case of the susceptibleinfected-susceptible (SIS) process. We derive formulas for asymptotic prevalence, extinction probabilities, and epidemic threshold. We also obtain a simple relationship between prevalence and extinction probability. Through the analysis it turns out that the epidemic thresholds are affected by the distribution of the attractiveness values assigned to the passive entities only through its second moment. Furthermore, comparing the dynamics of the indirect and direct spreading, we find that for low densities of active entities and relatively long survival times of contagions, the indirect spreading process is more effective than the direct spreading. The rest of the paper is organized as follows. In Sec. II we provide a brief survey of the state of the art. In Sec. III we formally describe the selection strategy with a detailed illustration of two specific schemes: uniform and preferential. In Sec. IV we describe the specific metrics we focus on and in Sec. V we explain the model of the spreading process along with the associated parameters. Section VI describes in detail the derivations of all the mathematical formulas for prevalence, extinction probabilities, and epidemic thresholds, for indirect spreading. In Sec. VII we show the derivation of the same metrics for direct spreading, along with a detailed comparison between the direct and indirect scenarios. II. RELATED WORKS

Plenty of analysis on epidemic processes has been done where an underlying static contact network among the agents or actors is assumed [18–28] and the contagions are assumed to be spreading through the links of this contact network. On the other hand, various studies have also been performed to understand the dynamical effects of human mobility on spreading. These works were mostly based on the concept of metapopulation where people live in several different communities and at a certain rate people of one community visit the

other community. For example, [29,30] systematically studies the effect of interpopulation movement rate and the spread of infection. The effect of the heterogeneous connectivity structure of the metapopulation network was analyzed in [31]. With the discovery of various sophisticated mobility models, researchers have also tried to understand the effect of more realistic features of the models on the spreading dynamics, e.g., the effect of objective traveling [32], dynamical condensation [33], commuting behavior of the agents between places [34], and self-initiated behavioral changes [35,36]. Specifically, the work in [32] explicitly considers a preferential selection of the destination by the agents and its effect on the spreading process. However, in our work we do a generalized analysis where preferential selection is a special case and focus on the indirect mode of spreading. Very recently, [37] studied the impact of the heterogeneity of the amount of the time spent by the travelers in different places, while [38] studied the impact of heterogeneity in the contact patterns on the dynamics of the spreading. However, neither of these works considers the common places, paths, groups, etc., as a separate class of entity with their own speciality. Hence, they fail to capture the notion of indirect transfer of the contagions, which is the prime issue focused on in the current work. Moreover, recent investigations on human mobility patterns [15,16] from collected empirical data reveal that the selection strategy adopted by the people to visit various places is mainly preferential [39]. A systematic explicit study of the effect of any such nonuniform selection strategy on epidemic spreading is lacking. It should be noted that, even apart from the spatial contagion processes, the selection strategies adopted by people is also important in the domain of various social contagion processes [17]. In this paper we provide a generalized analytical framework that can capture the effect of a specified selection strategy on the SIS epidemic process for both direct and indirect schemes. Recently, the work in [40] developed a framework based on compressed sensing theory for reconstruction of the underlying complex network from the small amount of available real data regarding the SIS stochastic spreading dynamics. To reflect the natural diversity they considered inhomogeneous infection and recovery rates.

032808-2

EPIDEMIC SPREADING THROUGH DIRECT AND INDIRECT . . .

In contrast to this, in the present work we study the dynamical process itself, with specific attention on investigating the effect of the diversity in the selection strategy that is fundamentally a reflection of diversity in the social behaviors of humans. On the other hand, the exact relationship between survival probabilities and infection prevalence is not widely available in the literature for many of the epidemic spreading models. Recently, [41] achieved such an exact relationship for the SIS process on a static contact network based framework. It should be noted that through our analytical framework we could provide a simple closed-form equation between the infection prevalence and survival probability for indirect and direct spreading, which may be governed by any underlying arbitrary selection strategy. III. SELECTION STRATEGY

In order to have a systematic analysis of the effect of the selection strategies adopted by the active entities, we assume that each of the passive entities i is associated with an attractiveness value, denoted by xi . Then xi / N i=0 xi (N denotes the total number of passive entities in the system) defines the probability that an active entity will make its next visit to the passive entity i. Since these probabilities are invariant under scaling of the xi values by a constant factor, we can assume without loss of generality that N i=0 xi = N , hence the mean value of xi equals one. We further assume that for N → ∞ the distribution of xi converges to a limiting ∞ distribution φ(x) with 0 xφ(x)dx = 1. As we will show in the next section, the above setting captures the important special case of the preferential attachment scheme. Form of φ(x)

Under various selection strategies, φ(x) will have a different form. Uniform selection. For uniform selection of the passive entities, every passive entity will have the same probability of being selected by the active entities and hence φ(x) will be a Dirac delta function with a peak at 1, φ(x) = δ(1 − x).

(1)

Preferential selection. A selection rule is called linear preferential selection if the probability that a passive entity is selected by an active entity in the next time step is proportional to the number of active entities that have already selected the passive entity up to the current time step : di (t) Pr{di (t + 1) = [di (t) + 1]} = N , j dj (t)

(2)

where di (t) denotes the number of active entities that have selected the passive entity i up to time step t. To start the process one has to fix the initial values {di (0)}. The most common choice here is di (0) = 1 for all i. To show that the preferential attachment rule can be equivalently described by a jump process according to an a priori fixed distribution of attractiveness values of the passive entities, one has to be more precise about the timing of the updates of the di (t). Although for the infection dynamics later on we will use a discrete synchronous update, we consider in

PHYSICAL REVIEW E 90, 032808 (2014)

the following a continuous time process for the jumps of active entities to the passive entities where active entities jump with rate 1 independently of each other. This is justified because the jump dynamics does not depend on the infection state.1 Let {tk } be the sequence of event times, that is, the time steps at which a jump of some active entity to some passive entity happens. Clearly, with probability one, no multiple jumps at the same time will happen. After each jump we update di (tk ), i.e., if the active entity jumps on the passive site j at time tk , we set dj (tk ) = dj (tk−1 ) + 1. The sequence {di (tk )} can be equivalently obtained by a multicolor P´olya urn process where at time tk there are di (tk ) balls of color i in the urn. At tk+1 a ball gets drawn at random from the urn and is put back into the urn together with an extra ball of the same color. Obviously, the probability to chose color i at tk is just the probability as given by Eq. (2) for the preferential attachment rule. For each i let ui (k) be the zero-one random variable specifying whether or not an active entity jumped at tk to passive entity i. The sequence {ui (k)}1iN,k0 uniquely encodes a realization of the urn process. It is well known that {ui (k)} is an exchangeable sequence, meaning that for each i and L the probability to observe {ui (k)}Lk=0 depends only on the number L k=0 ui (k) of ones in the sequence but not on its ordering. A famous theorem by de Finetti states that if a sequence {ui (k)} is exchangeable then there exists a probability distribution {ψi } with support in [0,1]N such that the sequence {ui (k)} can be equivalently obtained by sampling an independent and identically distributed Bernoulli process with a priori probabilities σi for each color i where σi are jointly sampled from {ψi }. For the specific case of multicolor Polya urn processes the probability distribution {ψi } is known to be the Dirichlet distribution with the parameters (d1 ,d2 , . . . ,dN ), where di = di (0), i.e., the initial value of the number of balls with color i. In the special case of a uniform initial condition with di (0) = 1 for all i, the marginal distribution of σi is a β distribution with parameters (1,N − 1). See [42,43] for interesting applications in the context of evolving bipartite networks. Returning to our preferential attachment scheme of the active entities, we obtain a realization by the following procedure: First, for each passive entity i a parameter σi , sampled jointly from a Dirichlet distribution, is preassigned and (b) second, the event of a certain passive entity i to get selected by some active entity can be considered a Bernoulli process with success rate σi . The following representative result of Dirichlet distributed random variables is especially useful for largeN values since it enables us to work in the limit of N going to infinity with independent exponential distributions. Let the random variables {σ1 ,σ2 , . . . ,σN } jointly follow a Dirichlet distribution with the parameters {d1 , . . . ,dN }, where di = 1∀i. Then {σ1 ,σ2 , . . . ,σN } can be represented as {X 1 /SN ,X2 /SN , . . . ,XN /SN }, where Xi ∼ exp(1) and SN = Xi ∼ (N,1), 1 i N . 1

Instead of continuous time one could also consider the discretetime scenario as we use it for the epidemic spreading but with additional microtime steps for the active entity jumps such that the update of the agent positions, and hence the update of the degrees d(i) of the passive entities, is sequential with respect to the microtime.

032808-3

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

PHYSICAL REVIEW E 90, 032808 (2014)

By the law of large numbers we have limN→∞ SN /N = 1. Furthermore, by the central limit theorem the fluctuations of √ SN around the mean are of order N, hence we can for large N well approximate the attractiveness variables xi = Nσi = Xi SNN = xi [1 + O( √1N )] by independent exponentially distributed random variables. In particular, the density of the limiting distribution of the attractiveness values is the exponential one:

The three observables studied in this paper are described below: (i) prevalence, i.e., the fraction of active and passive entities that are infected after the spreading process has continued for a long time; (ii) extinction probability, i.e., the probability that the infection dies; and (iii) epidemic thresholds, i.e., the critical values of the system parameters at which there is a transition of prevalence from zero to a nonzero value. Methodology. We study the prevalence in the associated limiting dynamical systems as nonzero fixed points. In order to analytically derive prevalence, we form recursive equations and solve them assuming the system reaches an equilibrium, that is, a fixed point. Consequently, we obtain the average prevalence. However, since we are dealing with a stochastic system, although the average prevalence is nonzero, there is always a probability that the contagion will get deleted from the system, especially when the process starts with a small number of infected entities. This is called the extinction probability. To compute the extinction probability, we use the mathematics of the branching process. To derive the epidemic threshold, we deduce the condition where the extinction probability becomes one so there is no possibility of having nonzero prevalence beyond that threshold.

(c) Infection of an active entity. We use a simple probabilistic contagion model. In indirect spreading, following the jumps, the active entities that come in contact with the already infected passive entities get infected with probability δ. However, in the case of direct spreading, we assume that an active entity gets infected with probability λ only when it comes into contact with at least one already infected active entity when they visit the same passive entity. (d) Infection of a passive entity. In the case of indirect spreading the passive entities that come into contact with at least one infected active entity get infected with probability γ . However, this step is not there for direct spreading. We describe two important issues related to buffer refreshing, which we use later in all our mathematical derivations. (i) Infection resident time. This is defined as the lifetime of a contagion until the buffer of an active (or passive) entity gets refreshed. Assuming a probabilistic deletion of the contagion at each time step, the resident time follows a geometric distribution with a mean 1−b , where b is the corresponding b buffer refresh probability. (ii) Number of active entities visiting a passive entity. We assume that in a time step each of the active entities makes a jump to one of the passive entities. According to the definition of the attractiveness values, if a passive entity has attractiveness x, the probability that it will be selected by an active entity is x . Thus, the number of active entities selecting a given passive N entity at a given time step will follow a binomial distribution with parameters Nx and M. Since we consider large N and M with constant ratio M → C, the binomial distribution can be N approximated by a Poisson distribution with mean xC. We use variants of this result in many places of our derivations where we explicitly mention the mean of the distribution. In the next two sections we derive analytical expressions for the abovementioned observables considering an arbitrary distribution φ(x) of the attractiveness values.

V. MODEL DESCRIPTION AND SYSTEM PARAMETERS

VI. INDIRECT SPREADING

In both direct and indirect spreading, we consider a system containing N passive and M active entities. We also assume that N → ∞ and M → ∞ such that M → C, where C is a N nonzero positive constant. The spreading process is initiated by one or more a priori infected active or passive entities. At each subsequent time step, a series of events is assumed to occur. The order of these events is stated below. (a) Buffer refresh. At the beginning of each time step the buffers in the passive entities and the active entities are refreshed (contagion2 is deleted) with probabilities α and β, respectively. In the case of direct spreading, contagions are refreshed with probability β from the active entities only, since passive entities do not directly participate in such spreading. (b) Passive entity selection. In both direct and indirect spreading, each active entity selects one passive entity that it visits. The selection of a passive entity is decided using a certain strategy (e.g., uniform or preferential).

In this section we provide a detailed analysis of the prevalence, extinction probabilities, and epidemic thresholds in the context of indirect spreading. We follow the model of the process as specified in Sec. V.

φ(x) = e−x .

(3)

IV. METRICS

2 We call the entity that is being spread, e.g., the pathogen in the case of disease spread or the messages or information in the case of technological or social communications, contagion.

A. Prevalence

We perceive the system as a discrete dynamical system with synchronous updates and derive recursive equations for the expected number of infected active and passive entities at time t + 1 as a function of the values at t. Each entity (active or passive) in the system may be in either of the two states infected or susceptible. Let A(t) and P (t) be the expected number of active and passive entities in the infected state at time t, respectively, with initial conditions A0 and P0 . Due to the law of large numbers, the true random variables of the number of infected active and passive entities are essentially localized around their expectation. We will therefore in the following not distinguish between them. It should also be noted that for finite N the infection will, with probability one, eventually become extinct; however, the time for this to happen is exponentially large in N . The interpretation of the asymptotic prevalence

032808-4

EPIDEMIC SPREADING THROUGH DIRECT AND INDIRECT . . .

should therefore be such that it gives the right value of infection for very large but fixed T in the limit of N → ∞. Fractions are denoted by lowercase letters, e.g., p(t) and a(t) representing PM(t) and A(t) , respectively. The probability N px (t + 1) that a given fixed passive entity with attractiveness x will be infected at time step t + 1 can be computed as the summation of the two following parts: (a) the probability that it was not infected at time t and got infected by an active entity at t + 1, which is possible if at least one infected active entity jumps to the passive entity and the passive entity gets infected from it, and (b) the probability that it was infected at time t and did not get refreshed at time t + 1. The number of infected active entities that did not refresh their buffers and could effectively transfer the infection to a passive entity is (1 − β)γ A(t). The probability that a fixed passive entity with attractiveness x is selected by an active entity is Nx . Thus, the probability that an infected active entity jumps to a passive entity and the passive entity gets infected can be calculated as 1 − (1 − Nx )(1−β)γ A(t) . Taking the limits → C as well as keeping N → ∞ and M → ∞ such that M N in mind that A(t) = CN a(t), this can be simplified as 1 − e−C(1−β)xγ a(t) . Combining parts (a) and (b), we get the following recursive equation for px (t + 1): px (t + 1) = (1 − α)px (t) + [1 − (1 − α)px (t)] × (1 − e−C(1−β)xγ a(t) ) = 1 − [1 − (1 − α)px (t)]e−C(1−β)xγ a(t) .

(4)

In a similar way the fraction of the active entities that are infected at time t + 1 can be derived using the recursive equation ˜ a(t + 1) = (1 − β)a(t) + [1 − (1 − β)a(t)]δ(1 − α)p(t),

PHYSICAL REVIEW E 90, 032808 (2014)

Inserting expression of p˜ ∗ into Eq. (6), we get the fixed point equation for the active entities as a ∗ = (1 − β)a ∗ + [1 − (1 − β)a ∗ ]δ(1 − α) ∞ ∗ xe−C(1−β)xγ a × 1−α φ(x)dx . 1 − (1 − α)e−C(1−β)xγ a ∗ 0

For later purposes we get the following equations after rearranging the terms in Eqs. (8) and (5): (1 − β)a ∗ =1−

(1 − α)px∗ = 1 −

qx =

i

xi −→ N N→∞

0

∞

(5) (6)

=

0

∞

1−

∗ αe−C(1−β)xγ a φ(x)dx. 1 − (1 − α)e−C(1−β)xγ a ∗

∞

e−xCδl+qxCδl (1 − α)l α

l=0

=

α . 1 − (1 − α)e−xCδ(1−qa )

ρx = 1 −

px φ(x)dx

=

qak Prob {buffer lifetime is l}

Therefore, the survival probability ρx can be written as

px∗ xφ(x)dx.

p =

∞ ∞ l=0 k=0

The fixed point equation for the passive entities becomes ∗

(10)

× Prob {k active entities jump to the passive entity} ∞ ∞ 1 k l k −xCδl (xCδl) e = qa [(1 − α) α] k! k l=0

with p˜ ∗ = p∗ being the weighted total prevalence of the N xi passive entities. Since the xi are asymptotically φ distributed, for N → ∞ we get px∗i

α . 1 − (1 − α)e−C(1−β)xγ a ∗

For a given attractiveness x of a passive entity and given lifetime l after infection we have a Poisson distributed number of infected offsprings P(xCδl). Thus, we have

xi

φ(x)dx

To derive the extinction probability we formulate the following recurrence equation. Let qx be the extinction probability for a single infected passive entity with attractiveness x and qa the extinction probability of an infected active entity. Let further ξx be the random variable for the number of offsprings of a passive entity with attractiveness x with lifetime l. We have qx = E qaξx .

∗

p˜ ∗ =

xe−C(1−β)xγ a ∗ 1−(1−α)e−C(1−β)xγ a ∗

B. Extinction probabilities

˜ = where p(t) The equations for the fixed points {px∗i ,a ∗ } of the system are

a ∗ = (1 − β)a(t) + [1 − (1 − β)a ∗ ]δ(1 − α)p˜ ∗ ,

0

and

xi i N pxi (t).

αe−C(1−β)xi γ a , 1 − (1 − α)e−C(1−β)xi γ a ∗

β + (1 − β)δ(1 − α)

β ∞

(9)

px∗i = 1 −

(8)

(7)

α . 1 − (1 − α)e−xCδρa

(11)

For the active entities we denote by χ the random variable for the number of offsprings, that is, the number of infected passive entities it generates, and have

χ qa = E qxi , i

032808-5

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

PHYSICAL REVIEW E 90, 032808 (2014)

where the index i enumerates the χ offsprings of the active entity. This gives qa =

∞ ∞

Prob {buffer lifetime is k}

l=0 k=0

× Prob {l out of k passive entities get infected} ×

l

{expectation of qx (i)}

i

=

∞ ∞ k l [(1 − β)k β] × γ (1 − γ )k−l l k=0 l=0 l × E(qxi ) i

=

∞ l=0

=

∞ l=0

=

l βγ l (1 − β)l E(qxi ) (γ + β − γβ)l+1 i

βγ l (1 − β)l (γ + β − γβ)l+1

∞

l

qx xφ(x)dx 0

β ∞ . β + γ (1 − β) 1 − 0 qx xφ(x)dx

Therefore, the corresponding survival probability ρa can be written as ρa = 1 −

β + γ (1 − β)

β ∞ 0

ρx xφ(x)dx

(12)

.

Inserting Eq. (11) into Eq. (12) gives ρa = 1 −

β + γ (1 − β) 1 −

β ∞ 0

αxφ(x) 1−(1−α)e−xCδρa

dx

.

(13)

We would finally like to remark that in the case where one has ka initial infected active entities and kp initial infected passive entities with attractiveness values x1 , . . . ,xk then the extinction probability of the associated branching process is kp qxi . just given by qaka i=1 C. Relationship between survival probability and prevalence

Comparing Eqs. (10) and (11) for the passive entities and Eqs. (9) and (13) for the active entities, one sees a striking similarity. The only difference is the exchange of the positions of δ and γ in the expressions and the scaling factors 1 − β and 1 − α, respectively. Let us define ρˆa and ρˆx as the survival probabilities of a branching process where the values of γ and δ are interchanged in comparison to the branching process used for ρa and ρx . We get a simple relationship between the prevalence and survival probability: a∗ =

ρˆa , 1−β

(14)

px∗ =

ρˆx . 1−α

(15)

In the following we provide a heuristic explanation of the above expressions. The basic idea is to construct a proper random graph space that for each given entity at a given time represents the potential paths of infection from the source of the infection to that entity. The probability that an entity at large time t is infected is then given by the probability that there is a connected path to the source in this random graph model, which is in essence the probability that corresponding branching process that goes backward in time is overcritical. To do so we first fix an intermediate time T0 1 such that the coverage of infected passive and active entities is roughly its asymptotic value. Let us denote the set of active and passive infected entities at time T0 by A(T0 ) and B(T0 ), respectively. To represent the infection process we first define a directed bipartite network. We generate copies for each passive entity i and denote it by the pair (i,τk (i)), where τk (i) denotes the kth lifetime of the passive entity i. Hence, {τk (i)}k0 is a sequence of independent and identically distributed random variables with values in N. The τk (i) are geometrically distributed with the parameter α and mean 1−α . Thus, each pair (i,τk (i)) corresponds to a different α vertex in the passive partition. In complete analogy we also generate copies for each active entity j and denote each by the pair (j,υl (j )), where {υl (j )}l0 is a sequence of independent and identically and geometrically distributed random variables with the parameter β and mean 1−β . Again the pairs (j,υl (j )) β correspond to nodes in the active partition. The time intervals a a Ii,l when the lth copy of active entity i exists is given by Ii,k = k k−1 [ k =0 τk (i), k =0 τk (i)] and for passive entities the analog l p p interval Ij,l is given by Ij,l = [ l−1 l =0 υl (j ), l =0 υl (j )]. Directed edges between (i,τk (i)) and (j,υl (j )) will now be generated according to the following probabilistic rule. At each p a time step t ∈ Ii,k ∩ Ij,l we first make a random allocation of the active entity i on the passive entities corresponding to the selection rules described in Sec. III. If i is allocated to j (this xi happens with probability N ) we draw at time t a directed edge from (i,τk (i)) to (j,υl (j )) with probability γ and a directed edge from (j,υl (j )) to (i,τk (i)) with probability δ. Note that the edges at different time steps are drawn independently but are not independent within a given time step due to the unique allocation of active entities to a passive entity according to the given selection strategy. Edges represent infection events with the direction chosen such that the edges always point to the source of the infection. Directed paths hence can be interpreted as going backward in time. To get the relevant network up to time T one has to truncate all nodes (i,τk (i)) for which k−1 k =0 τk (i) T and correspondingly for active partition nodes. For a passive entity i to be infected at time t > T0 , there has to be at least one directed path from (i,τk (i)) with t ∈ τk to the set of nodes in the bipartite network that represent the infected population A(T0 ) ∪ B(T0 ) at time T0 . Similarly, an active entity j is infected at time t if there is a directed path from (j,υl (j )), i.e., the representative of node j at time t if t ∈ υl (j ), to A(T0 ) ∪ B(T0 ). The probability that for a given active or passive entity i at large t there is a connecting path to A(T0 ) ∪ B(T0 ) in the limit of N → ∞ is the same as the probability that the directed, backward in time branching process describing the local growth of the network

032808-6

EPIDEMIC SPREADING THROUGH DIRECT AND INDIRECT . . .

D. Epidemic thresholds

For an arbitrary selection strategy, the epidemic threshold can now be estimated by setting the derivative with respect to ρa [on the right-hand side g(ρa ) in Eq. (13)] at zero equal to one. This gives 1 d 1 g(ρa )|ρa =0 = 1 = γ (1 − β) (1 − α)Cδ x 2 φ(x)dx dρa β α (1 − β)(1 − α) Cγ δμ2 , (16) = αβ where μ2 is the second moment of φ(x). The same equation can be also derived by setting the derivative with respect to ρx [on the right-hand side g(ρx ) of Eq. (11)] [after putting Eq. (12) into it] at zero to one. Furthermore, it is to be noted that the threshold equations can be derived also from the fixed point equation for the prevalence, i.e., Eq. (8) or (7), in a similar way. The corresponding threshold equations for preferential and uniform selection strategies can be obtained just by setting μ2 = 2 and 1, respectively, in Eq. (16). E. Discussion

In this section we discuss a few important issues regarding the process of indirect spreading. First we show the effect of a large value of C and then compare the theoretical predictions with simulation results for finite value of N . Next we try to show the importance of the second moment of the distribution of the attractiveness values in several cases through both simulation and the developed theory. Effect of C 1. It is interesting to note that for C 1 and δ > 0, from Eq. (13) one gets ρa ∼ 1 −

β , β + γ (1 − β)(1 − α)

(17)

which is independent of δ as well as any specific selection strategy. Similarly, from Eq. (11) one can obtain ρx ∼ 1 − α.

(18)

K = 1000

1

K = 100 0.8 K = 10 K=2

0.6 K=1

β

in a neighborhood of i does not become extinct. Similarly, the extinction of the corresponding branching process starting in, say, (i,τk (i)) is equivalent to the property that (i,τk (i)) is a node of a small component of the directed bipartite network that is not connected to the giant component. This explains the exchange of the infection probabilities in the correspondence between the prevalence and survival probabilities. The scaling 1 1 factors 1−β and 1−α in Eqs. (14) and (15), respectively, are due to the fact that by going backward in time we have to condition on survival for at least one time step since, in contrast to the forward in time calculations where an infected node can have zero offsprings due to the fact that it is rebuffered before the infection substep takes place, we have the role of events reverted in the backward picture. That is, the node can get infected even if it rebuffers afterward. This corresponds in the calculations just to a conditioning on survival for at least one time step. The above-constructed bipartite network, but with a forward direction of time, has been of course implicitly used in the construction of the branching process that defined the extinction probabilities in Sec. VI B.

PHYSICAL REVIEW E 90, 032808 (2014)

0.4 K = 0.5 0.2

0 0

K = 0.1

0.2

0.4

α

0.6

0.8

1

FIG. 2. (Color online) Plot of the threshold equation (16) for different values of K (K stands for Cγ δμ2 ); γ and δ are assumed to be 1.

From Fig. 2 it is shown that for large values of C, the epidemic threshold almost vanishes. From the simplified equations for the survival probabilities [Eqs. (18) and (17)] it is shown that the quantities become independent of a specific selection strategy. However, it is interesting to note that even for a very large value of C, the survival probabilities do not become 1, i.e., it still has a considerable chance of becoming extinct. It is also interesting to note that for large C, ρa does not depend on δ, whereas ρx depends only on α. Effect of small and finite N . In our analytical derivations, we assume a large value of N . Hence, to see how well the expressions can predict the observables for finite N we compare the results obtained from theory with the same obtained from simulation. Figure 3 shows this comparison for the prevalence obtained from the theory with simulation of the indirect spreading process for preferential selection for several small values of N . It is shown that after a certain increase in the value of N (approximately above 1000), the asymptotic formulas more accurately match with the simulation. However, the match is less accurate for very small values of N . This is mainly due to the fact that the approximation of the Dirichlet distributed attractiveness values associated with the passive entities (for preferential selection) by independent exp(1) random variables holds only for large N . However, it is also evident from Fig. 3 that the value of N above which one has a very good match between theory and simulation depends also on the parameter combination. For example, for higher buffer refreshing probabilities the value of N is required to be much higher for a better match, which might be due to the decrease in the effective number of passive entities participating in the spreading process. Symmetric behavior of the parameters. Equation (16) and Fig. 2 reveal a symmetric effect of the pair of parameters α and β in the epidemic thresholds. In other words, interchanging the values of α and β does not alter the threshold values of the system. The same is more apparent for the pair of parameters γ and δ from Eq. (16) itself. Since the roles of the active and the passive entities are very different in the spreading process,

032808-7

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

1

α=0.1, β=0.1, γ=0.4, δ=0.8

0.95 Active entity prevalence

PHYSICAL REVIEW E 90, 032808 (2014)

α=0.4, β=0.1, γ=0.3, δ=0.3

0.6

0.9 0.4

0.85 0.8 0.75

α=0.1, β=0.1,

Passive entity prevalence

0.2

(a) γ=0.8, δ=0.4

0.7 1 10 0.8

Theory Simulation Theory Simulation

α=0.1, β=0.4, γ=0.3, δ=0.3

(b) 2

0 3 x 103 101 0.4

3

10

10

α=0.1, β=0.1, γ=0.4, δ=0.8

0.3

2

3

3

10

3 x 10

10

α=0.4, β=0.1, γ=0.3, δ=0.3

0.7 0.2 0.6 0.5

(c)

0.4 1 10

α=0.1, β=0.1, γ=0.8, δ=0.4

0.1

2

0 1 3 x 103 10

3

10

α=0.1, β=0.4, γ=0.3, δ=0.3

(d)

10

2

3

10

N (= M)

10

3 x 103

FIG. 3. (Color online) Comparison between the the active (and passive) entity prevalence obtained from theory and the same obtained from the simulation of the indirect spreading process under preferential selection for different finite values of N and under a few different combinations of the parameters. The theoretical results are calculated from Eqs. (8) and (7) by replacing φ(x) using Eq. (3). In the simulation, for each different parameter combination and for each different value of N , the indirect spreading processes were run for 1000 time steps and the results are computed by averaging the achieved prevalence values from time steps 500 to 1000. The results are plotted on a linear-logarithmic scale.

this symmetry is not so obvious. However, this symmetry is true for neither the prevalence (see Fig. 3) nor the survival probability [see Eqs. (11) and (13) and even their simplified forms (17) and (18)]. Effect of second moment of φ(x). The second moment μ2 of φ(x) plays a very crucial role in the dynamics of spreading. Plots of the epidemic threshold (16) in Fig. 2 show that for given values of C, γ , and δ, the subcritical region (area above a specific curve in Fig. 2) decreases with μ2 . Thus, higher diversity in the selection of the passive entities implies a higher chance of getting infected. Figure 4 illustrates this fact through

the plot of the prevalence, survival probability, and epidemic threshold when φ(x) is a power law (see details in the figure caption). Due to diverging μ2 , the subcritical region almost vanishes for the power-law exponent ζ = 2.1 [Fig. 4(a)], while for ζ = 3.1, a considerable region is subcritical [Fig. 4(b)]. Similarly, Fig. 4(a) shows that, as ζ goes below 3.0, the survival probability rapidly increases even under high α and β and low γ and δ. A comparison of prevalence under preferential (μ2 = 2) and uniform (μ2 = 1) selection in Fig. 5 shows a similar effect of μ2 on prevalence (see the figure caption for details).

0.05 0.8

0.04

0.6

0.6

0.03

0.4

0.4

α

(a)

0.2

ρa

0.8

10

−5

0.02

Theory

α= 0.5, β=0.5 α= 0.6, β=0.6 α= 0.7, β=0.7 α= 0.8, β=0.8 α= 0.9, β=0.9

0

10

(c)

(b) 0.2

0.01 2

0

0 0

0.2 0.4 0.6 0.8

β

0

0.2 0.4 0.6 0.8

0 2

2.4

3 2.8

τ

3.2

4 3.6

4

FIG. 4. (Color online) (a) and (b) Active entity prevalence obtained through simulation of the indirect spreading with C = 1 (N = 200 and M = 200) in the α − β plane where the attractiveness parameters associated with the passive entities follow a power-law distribution −1 x −ζ ( xmin ) , where the exponent ζ is equal to (a) 3.1 and (b) 2.1. For different ζ , we adjust the value of xmin in order to make the mean φ(x) = xζ min μ1 of the distribution equal to 1. The xmin values are (a) 0.54 and (b) 0.88. The empirically measured values of the second moment μ2 are (a) 5.71 and (b) 24.5. The dark black points indicate zero prevalence, while the gray and white points indicate nonzero prevalence. (c) Plot of the active entity survival probabilities derived from Eq. (13) for different values of the power-law exponent ζ , under a few different combinations of the parameters α and β (the values of transmission probabilities γ and δ are both 0.1 and the value of C is 1). The inset in (c) shows the same plot on a log-log scale. 032808-8

EPIDEMIC SPREADING THROUGH DIRECT AND INDIRECT . . . Uniform

PHYSICAL REVIEW E 90, 032808 (2014) Preferential 1

1

1

0.8

0.8

0.8

0.8

0.6

0.6

0.6

Active entity prevalence

α = 0.12

0.6

1

α = 0.2

0.4

0.4

(a)

0.2 0 0

0.2 0.4 0.6 0.8

0.2 1

0 0

α = 0.32

α = 0.44

0.4

0.4

(d)

(c)

(b)

0.2

0.2

0.2 0.4 0.6 0.8

1

0 0

0.2 0.4 0.6 0.8

1

0 0

0.2 0.4 0.6 0.8

1

β

FIG. 5. (Color online) Comparison between the active entity prevalence achieved in preferential selection and uniform selection of the passive entities in an indirect spreading process with C = 1 (N = 100 and M = 100) for a few different values of α and β (γ = δ = 1). VII. DIRECT SPREADING

In the following we derive the fixed point equation for the prevalence, the extinction probabilities, and the threshold equations in the context of a direct spreading process. We follow the model of the process as specified in Sec. V. A. Prevalence

In the case of direct spreading, the fraction of active entities that are infected at time step t + 1 can be calculated by combining the following two parts. (a) The first is the fraction of active entities that were not infected at time t or refreshed their buffer at t, but visited a certain passive entity where at least one infected active entity also visited. The number of infected active entities that did not refresh their buffers and could effectively transmit the infection to the others is A(t)λ(1 − β). Thus, the number of active entities that can effectively transmit the infection to the others by selecting a passive entity with attractiveness x to visit is distributed as P[xCa(t)λ(1 − β)]. Hence, the probability that at least one such active entity visits a passive entity i with attractiveness value xi can be calculated as 1 − e−xi Cat λ(1−β) . Taking the average over all possible attractiveness values, we get N xi (1 − e−xi Cat λ(1−β) ). N i=1

spreading can be obtained as a ∗ = 1 − [1 − a ∗ (1 − β)]

∞

xφ(x)e−xλCa

(1−β)

dx. (21)

0

B. Extinction probabilities

Let X be a random variable denoting the number of infected active offsprings produced by an infected active entity during its lifetime. The probability that the lifetime of the buffer of the active entity is l is given by (1 − β)l β. Now, in these l time steps the active entity visits l different passive entities and at each of them it meets with a certain number of other people who it infects. The number of other active entities that get infected by this active entity at a passive entity with attractiveness x is distributed as P(xCλ). Thus, the distribution of the number of active entities that get infected in the whole period of l time steps will be given by P( lm=1 x (m) Cλ). Here x (m) is just the random variable denoting the attractiveness of those passive entities that are visited by an infected active entity at the mth time step of its lifetime. Thus, the extinction probabilities can be derived in a way similar to that in Sec. VI B as qa = E qaX qak Prob {buffer lifetime is l} = k

× Prob{k active entities get infected} k ∞ ∞ l l 1 k (m) (1 − β) β qa · · · x Cλ = k! k l=0 m=1 l l (m) × exp − x Cλx x (m) φ(x (m) )dx (m)

For large N , the expression can be approximated by the following integral: ∞ xφ(x)(1 − e−xi Cat λ(1−β) )dx. (19) 0

m=1

(b) The second part is the fraction of active entities that were infected at time t but did not refresh their buffers. Combining parts (a) and (b) and using Eq. (19) we get the following recursive equation for a(t + 1): a(t + 1) = a(t)(1 − β) + [1 − a(t) + βa(t)] ∞ −xi Cat λ(1−β) × xφ(x)(1 − e )dx .

∗

l=0

= (20)

0

m=1

and for N → ∞ one obtains l ∞ (1 − β)l β qa = xe−xCλ(1−qa ) φ(x)dx

1 − (1 − β)

β . xe−xCλ(1−qa ) φ(x)dx

For the survival probabilities we have

Thus, similar to Sec. VI A, the fixed point equation for the asymptotic prevalence of the infected active entities for direct 032808-9

ρa = 1 −

1 − (1 − β)

β . xe−xCλρa φ(x)dx

(22)

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

PHYSICAL REVIEW E 90, 032808 (2014)

FIG. 6. (Color online) (a)–(c) Plot of the difference between the active entity prevalence values in indirect and direct spreading through preferential selection using Eqs. (8) and (21) after replacing φ(x) by Eq. (3). (d)–(f) Same as (a)–(c) but for the survival probabilities [obtained from Eqs. (13) and (22)]. The plots are shown in the γ -β plane for a few different values of C. The thick and thin black lines are the plots of the epidemic thresholds for indirect and direct spreading [Eqs. (16) and (24)], respectively. Subcritical (right) and overcritical (left) regions are separated out by these threshold lines for each different kind of spreading. In (b) and (e) indirect and direct threshold lines coincide with each other. In (a), (c), (d), and (f) the regions inside the two lines indicate the cases where one is overcritical while the other is subcritical. The shaded parts denote the regions where the prevalence and the survival probability are higher in indirect spreading than in direct spreading. The white parts in the overcritical regions indicate the cases where prevalence and survival probability in direct spreading are higher than or equal to the same in indirect spreading.

C. Relationship between survival probability and prevalence

E. Discussion

Comparing Eqs. (21) and (22) and after few steps of simplifications we get the following simple relationship between the survival probability and the prevalence for direct spreading: ρa . (23) a∗ = 1−β

The analytical expressions for survival probabilities, prevalence, and epidemic thresholds derived for direct spreading show that the effect of the ratio C, transmission probability γ , the second moment of φ(x), i.e., μ2 , and the buffer refresh probability β are completely analogous to the indirect spreading case. One of the principal motivations behind the analysis of the direct spreading is to compare it with the case of indirect spreading, which is done next.

The explanation in terms of a backward in time branching process is analogous to the explanation in the indirect spreading case with the important difference that the corresponding directed graph is not bipartite. Furthermore, since there is only one transmission probability involved, the backward and forward associated branching processes are equivalent up to a scaling by 1 − β. D. Epidemic thresholds

Using the same technique used in Sec. VI D, the epidemic threshold condition can now be obtained from Eq. (22) as follows: d −xCγρa 1−β 1− (1− β) φ(x)dx = 1, xe dρa ρa =0 which gives 1=

1−β Cλ β

x 2 φ(x)dx = (1 − β)

Cλμ2 . β

(24)

The same equation for the epidemic threshold can also be derived from Eq. (21) by differentiating with respect to a ∗ , setting a ∗ = 0, and equating with 1. For (1 − β) Cγβμ2 > 1 there exists a unique nonzero solution that gives the survival probability for the associated branching process starting with one initially infected active entity.

F. Comparison of direct and indirect spreading

It is understood that indirect spreading in general would be more effective than direct spreading as the contagion lingers in the place even after the agent has left it. However, a detailed comparison shows that this is not universally true. To compare direct and indirect spreading under the same settings, e.g., for the same selection strategy, we consider α = β and γ = δ = λ. Hence, in the following we only use the parameters β, γ , and C. The comparison of the threshold equation for direct and indirect spreading [Eqs. (16) and (24)] yields an interesting observation. We find two different regimes for the value of Cμ2 : For Cμ2 < 1, the plot of the threshold equation for direct spreading is strictly below the same for the indirect spreading, while for Cμ2 > 1, the indirect spreading is strictly below the same for the direct spreading. This implies that, for Cμ2 < 1, there are parameter configurations where the indirect spreading is overcritical but the direct spreading is still subcritical and vice versa. Figures 6 and 7 show such regions for preferential selection (μ2 = 2) in the γ -β and γ -C planes, respectively. They also highlight the regions where the active entity survival

032808-10

EPIDEMIC SPREADING THROUGH DIRECT AND INDIRECT . . .

PHYSICAL REVIEW E 90, 032808 (2014)

FIG. 7. (Color online) Same as in Fig. 6 but in the γ -C plane for a few different values of β.

probabilities and the prevalence values are higher in indirect spreading than the corresponding direct spreading. Observations. It is shown in Fig. 6 that for preferential selection, except the region where β is small and γ is high, for values of C less than or equal to 0.5, direct spreading is more efficient than indirect spreading. The same is also depicted in Fig. 7. VIII. CONCLUSION

Understanding the effect of nonuniform selection rules followed by the active entities was the prime target of the analysis. A mathematical framework has been introduced that attempts to study this effect through the assignment of attractiveness values, following an arbitrary distribution, to the passive entities. One of the important findings in this direction is the fact that the diversity in the popularity in the frequently visited places as well as commonly used objects implies an enlargement of the overcritical region in the parameter space as well as wider spread of the contagion. For the special but important case of the preferential attachment rule we could derive explicit asymptotic results. Moreover, we obtain a simple relationship between survival probabilities and the prevalence, which is interesting in the context of a deeper understanding

[1] E. Cator and P. Van Mieghem, Phys. Rev. E 87, 012811 (2013). [2] H. K. Lee, P.-S. Shim, and J. D. Noh, Phys. Rev. E 87, 062812 (2013). [3] S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, Phys. Rev. E 86, 041125 (2012). [4] M. Bogu˜na´ , C. Castellano, and R. Pastor-Satorras, Phys. Rev. Lett. 111, 068701 (2013). [5] R. Anderson and R. M. May, Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, 1992). [6] F. Brauer and C. Castillo-Chvez, Mathematical Models for Communicable Diseases, CBMS-NSF Regional Conference Series in Applied Mathematics (SIAM, Philadelphia, 2012). [7] N. Banerjee, M. D. Corner, and B. N. Levine, INFOCOM 2007, Proceedings of the 26th International Conference on Computer Communications (IEEE, Piscataway, NJ, 2007).

of the relationship between the SIS spreading dynamics and the percolation dynamics of susceptible-infectious-removed models. The comparison between the direct and indirect spreading shows a surprising phenomenon: For lesser values of Cμ2 , there are regions in the parameter space where indirect spreading has higher survival probability, prevalence, and a larger overcritical region than the direct spreading. Interpreted in epidemic language, the finding implies the following for a given selection strategy: If the number of active entities is sufficiently small compared to the number of passive entities, indirect spreading becomes a dominating mechanism in disease spreading. An important remaining open issue is the stochastic analysis of the system for small values of N . Here one sees a clear mismatch with the asymptotic theory. ACKNOWLEDGMENTS

S.S. thanks Tata Consultancy Services and Samsung Pvt. Ltd. for financial assistance. T.K. thanks IIT Kharagpur and Wroclaw University of Technology for financial support for his research stay in India. Specially, T.K. thanks IIT Kharagpur for hospitality and providing excellent working environment.

[8] W. Zhao et al., in 2006 IEEE International Conference of Mobile Adhoc and Sensor Systems (MASS) (IEEE, Piscataway, NJ, 2006), pp. 31–40. [9] K. M. Hewitt, C. P. Gerba, S. L. Maxwell, and S. T. Kelley, PLoS ONE 7, e37849 (2012). [10] S. A. Boone and C. P. Gerba, Appl. Environ. Microbiol. 73, 1687 (2007). [11] A. Kramer, I. Schwebke, and G. Kampf, BMC Infect. Dis. 6, 130 (2006). [12] S. A. Sirsat, J.-K. K. Choi, B. A. Almanza, and J. A. Neal, J. Environ. Health 75, 8 (2013). [13] T. Camp, J. Boleng, and V. Davies, Trends Appl. 2, 483 (2002). [14] C. Bettstetter, SIGMOBILE Mob. Comput. Commun. Rev. 5, 55 (2001).

032808-11

GANGULY, KRUEGER, MUKHERJEE, AND SAHA

PHYSICAL REVIEW E 90, 032808 (2014)

[15] A. Mei and J. Stefa, in INFOCOM 2009, IEEE (IEEE, Piscataway, NJ, 2009), pp. 2106–2113. [16] S. Lim, C. Yu, and C. Das, in Proceedings 2006 31st IEEE Conference onLocal Computer Networks (IEEE, Piscataway, NJ, 2006), pp. 231–238. [17] A. Ghosh, S. Ha, E. Crabbe, and J. Rexford, IEEE J. Sel. Areas Commun. 31, 2673 (2013). [18] M. J. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals (Princeton University Press, Princeton, NJ, 2008). [19] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Phys. Rep. 424, 175 (2006). [20] P. Csermely, Weak Links: Stabilizers of Complex Systems from Proteins to Social Networks (Springer, Berlin, 2006). [21] A. Vespignani, Nat. Phys. 8, 32 (2012). [22] M. Perc, J. G´omez-Garde˜nes, A. Szolnoki, L. M. Flor´ıa, and Y. Moreno, J. R. Soc. Interface 10, 20120997 (2013). [23] L. A. Meyers, B. Pourbohloul, M. E. Newman, D. M. Skowronski, and R. C. Brunham, J. Theor. Biol. 232, 71 (2005). [24] T. Zhou, J.-G. Liu, W.-J. Bai, G. Chen, and B.-H. Wang, Phys. Rev. E 74, 056109 (2006). [25] J. G´omez-Garde˜nes, V. Latora, Y. Moreno, and E. Profumo, Proc. Natl. Acad. Sci. USA 105, 1399 (2008). [26] S. Meloni, A. Arenas, and Y. Moreno, Proc. Natl. Acad. Sci. USA 106, 16897 (2009).

[27] R. Parshani, S. Carmi, and S. Havlin, Phys. Rev. Lett. 104, 258701 (2010). [28] C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 105, 218701 (2010). [29] L. A. Rvachev and I. M. Longini, Jr., Math. Biosci. 75, 3 (1985). [30] M. J. Keeling and P. Rohani, Ecol. Lett. 5, 20 (2002). [31] V. Colizza and A. Vespignani, J. Theor. Biol. 251, 450 (2008). [32] M. Tang, Z. Liu, and B. Li, Europhys. Lett. 87, 18005 (2009). [33] M. Tang, L. Liu, and Z. Liu, Phys. Rev. E 79, 016108 (2009). [34] V. Belik, T. Geisel, and D. Brockmann, Phys. Rev. X 1, 011001 (2011). [35] S. Meloni et al., Sci. Rep. 1, 62 (2011). [36] B. Wang, L. Cao, H. Suzuki, and K. Aihara, Sci. Rep. 2, 887 (2012). [37] C. Poletto, M. Tizzoni, and V. Colizza, J. Theor. Biol. 338, 41 (2013). [38] L. Wang, Z. Wang, Y. Zhang, and X. Li, Sci. Rep. 3, 1468 (2013). [39] A.-L. Barab´asi and R. Albert, Science 286, 509 (1999). [40] Z. Shen, W.-X. Wang, Y. Fan, Z. Di, and Y.-C. Lai, Nat. Commun. 5, 4323 (2014). [41] R. R. Wilkinson and K. J. Sharkey, PLoS ONE 8, e69028 (2013). [42] N. Ganguly, S. Ghosh, T. Krueger, and A. Srivastava, Theor. Comput. Sci. 466, 20 (2012). [43] S. Saha, N. Ganguly, A. Mukherjee, and T. Krueger, Phys. Rev. E 89, 042812 (2014).

032808-12