key: cord-0010895-36an9hpe
authors: Choi, Wonjun; Lee, Deokjae; Kahng, B.
title: Critical behavior of a two-step contagion model with multiple seeds
date: 2017-06-12
journal: Phys Rev E
DOI: 10.1103/physreve.95.062115
sha: e9e58a323a26d456b4b05966068c1a60ed969c4a
doc_id: 10895
cord_uid: 36an9hpe

A two-step contagion model with a single seed serves as a cornerstone for understanding the critical behaviors and underlying mechanism of discontinuous percolation transitions induced by cascade dynamics. When the contagion spreads from a single seed, a cluster of infected and recovered nodes grows without any cluster merging process. However, when the contagion starts from multiple seeds of [Formula: see text] where [Formula: see text] is the system size, a node weakened by a seed can be infected more easily when it is in contact with another node infected by a different pathogen seed. This contagion process can be viewed as a cluster merging process in a percolation model. Here we show analytically and numerically that when the density of infectious seeds is relatively small but [Formula: see text] , the epidemic transition is hybrid, exhibiting both continuous and discontinuous behavior, whereas when it is sufficiently large and reaches a critical point, the transition becomes continuous. We determine the full set of critical exponents describing the hybrid and the continuous transitions. Their critical behaviors differ from those in the single-seed case.

A two-step contagion model with a single seed serves as a cornerstone for understanding the critical behaviors and underlying mechanism of discontinuous percolation transitions induced by cascade dynamics. When the contagion spreads from a single seed, a cluster of infected and recovered nodes grows without any cluster merging process. However, when the contagion starts from multiple seeds of O(N ) where N is the system size, a node weakened by a seed can be infected more easily when it is in contact with another node infected by a different pathogen seed. This contagion process can be viewed as a cluster merging process in a percolation model. Here we show analytically and numerically that when the density of infectious seeds is relatively small but O (1) , the epidemic transition is hybrid, exhibiting both continuous and discontinuous behavior, whereas when it is sufficiently large and reaches a critical point, the transition becomes continuous. We determine the full set of critical exponents describing the hybrid and the continuous transitions. Their critical behaviors differ from those in the single-seed case. DOI: 10.1103/PhysRevE.95.062115

Nonequilibrium dynamic transitions driven by cascade dynamics on complex networks have attracted considerable attention recently [1] [2] [3] . The spreading of epidemic disease on complex networks [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] is an instance, in which a pathogen is transmitted from an infected node (e.g., a person) to a susceptible neighbor, who then becomes infected with a certain probability. If the transmission probability is sufficiently large (small), the pathogen spreads out to a macroscopic scale (remains local). An epidemic transition occurs between these two limits. The extent of spreading also depends on the structure of an underlying network [1, 24] . When the degree distribution of a network is highly heterogeneous, diseases can spread out massively even for a small transmission probability, so that an epidemic transition point can be zero [25] . Information spreading in social media from one page to others may be modeled in a similar manner [5, 6] .

Among the several epidemic models, one of the simple contagion models is the so-called susceptible-infected-recovered (SIR) model [26, 27] , in which each node has one of three states, susceptible (denoted as S), infected (I ), or recovered (R). Initially, all the nodes are in state S except for one seed node in state I . The contagion process starts from a single node in state I . Each node in state I transmits pathogens to its neighbors in state S and infects each of them with probability κ; then it changes its state to R with a unit probability. This contact process is repeated until the system reaches an absorbing state in which no infected node is left in the system. When the probability κ is sufficiently small (large), the order parameter defined as the density of nodes in state R after the system falls into the absorbing state becomes o(N ) [O(N )]; i.e., the system falls into a subcritical (supercritical) state. In between, an epidemic transition occurs at κ c , and the system exhibits critical behavior. It is known that when the dynamics starts from a single seed on Erdős-Rényi (ER) random networks [28] , the SIR model undergoes a continuous * bkahng@snu.ac.kr percolation transition following the universal behavior of ordinary percolation.

The SIR model with multiple seeds has been considered [29] , in which two percolation transitions occur successively at κ c1 and κ c2 as κ is increased. The density of nodes in state R is finite for κ > κ c1 , whereas the density of nodes in state S disappears for κ > κ c2 . Thus, there exists a state of coexisting nodes in states R and S between κ c1 and κ c2 .

The SIR model was extended to a two-step contagion model, in which a weakened state (W ) can exist between the S and I states. Accordingly, this model is called the SWIR model [7, 20] . Nodes in state W are involved in the reactions S + I → W + I and W + I → 2I , which occur in addition to the reactions S + I → 2I and I → R in the SIR model. The properties of the epidemic transition in the SWIR model were extensively investigated for the single-seed case [7, 10, 15, [20] [21] [22] . The order parameter, defined as the density of nodes in state R after an absorbing state is reached, displays a discontinuous transition, whereas other physical quantities such as the outbreak size distribution exhibit critical behaviors. Thus, the phase transition occurring in the SWIR model with a single seed is regarded as a mixed-order phase transition [22] . The dynamic rule of the SWIR model is so simple that its underlying mechanism for the discontinuous behavior of the order parameter was disclosed [30] . Moreover, the mechanism turned out to be universal in other models such as k-core percolation [31] [32] [33] [34] , the cascading failure model on interdependent networks [35] [36] [37] [38] [39] , and epidemic-related models [5, 6, 8, 10, [12] [13] [14] 16, 17] .

Here we investigate the phase transitions of the SWIR model with multiple seeds. The model with multiple seeds has been investigated in Refs. [15, 20, 40] : The authors of Refs. [15, 40] used the mean-field approach and performed numerical simulations, obtaining the phase diagram as a function of the reaction rates. The order parameter exhibits either a discontinuous or continuous transition depending on the density of the infectious seeds and mean degree of a given network [15, 40] . In Ref. [20] the discontinuous transition is regarded as a spinodal transition, because there is no coexistence phase in the system while the order parameter jumps. Even though such results were obtained, the properties of the phase transitions and critical behaviors were not deeply investigated yet.

Here we reveal that the spread of contagion in the SWIR model with multiple seeds proceeds differently from that in the SWIR model with a single seed: in the multiple-seed case, the reactions W + I → 2I often occur even in early time steps, because nodes in states W and I involved in that reaction can originate from different seeds (see Fig. 1 ). We note that the number of multiple seeds was taken as O(N). On the contrary, in the single-seed case, such reactions rarely occur until the system reaches a characteristic dynamic step n c (N ) ∼ N 1/3 : When dynamic step n is less than n c (N), the reactions S + I → 2I and I → R are dominant, but the number of nodes in R still remains as o(N ). The contagion spreads in the form of a branching tree. When the dynamics reaches n c (N ), the branching tree forms long-range loops due to finite-size effect. Once such loops form, the reaction W + I → 2I occurs massively, in which the nodes in state W were generated in early time steps. Thus, the density of nodes in state R increases drastically as many as O(N ) in short time steps. Due to these different contagion mechanisms, the properties of epidemic transitions in the multiple seed case become different from those in the single-seed case. We will determine the full set of critical exponents describing the phase transitions in the multiple seed case, and compare them with those obtained in the single-seed case [22] . This paper is organized as follows: In Sec. II we present the rules of the SWIR model in detail. In Sec. III we set up the self-consistency equation to derive the mean-field solution using the local tree approximation of the order parameter for the epidemic transition on the ER networks. We show that, depending on the initial density of infectious nodes, different types of phase transition can occur. In Sec. IV we report numerical results for the epidemic transitions. In the final section, a summary and discussion are presented.

The SWIR model with multiple seeds is simulated on ER networks with N nodes. Initially, Nρ 0 nodes are selected randomly from among those N nodes and assigned to state I ; the other nodes are assigned to state S. At each time step n, the following processes are performed: (1) All the nodes in state I are listed in random order. (2) The states of the neighbors of each node in the list are updated sequentially as follows: If a neighbor is in state S, it changes its state in one of the two ways: either to I with probability κ or to W with probability μ. If a neighbor is in the state W , it changes to I with probability η, where κ, μ, and η are the contagion probabilities for the respective reactions. (3) All nodes in the list change their states to R. This completes the time step, and we repeat the above processes until the system reaches an absorbing state in which no infectious node is left in the system. The reactions are summarized as follows:

The order parameter exhibits a discontinuous transition at a transition point κ c when ρ 0 is less than a critical value ρ c , and it shows a continuous transition at κ c when ρ 0 = ρ c for given parameter values z, μ, and η, where z is the mean degree of a given ER network.

In an absorbing state, each node is in one of three states: the susceptible S, weakened W , and recovered R states. The order parameter m(κ), the density of nodes in state R in an absorbing state, is written using the local tree approximation as

(5) The first term in Eq. (5), ρ 0 , is the initial density of infected nodes. In the second term, the factor (1 − ρ 0 ) represents the probability that a node is originally in state S. P d (k) is the probability that a randomly selected node has degree k; q is the probability that an arbitrarily chosen edge leads to a node that is in state R but not infected through the chosen edge in the absorbing state. Thus, P d (k) k q (1 − q) k− is the probability that a node has degree k and of them are in state R in the absorbing state. P R ( ) is the conditional probability that a node is finally in state R, provided that it was originally in state S and its neighbors are in state R in the absorbing state.

Similarly to P R ( ), we define P S ( ) as the conditional probability that a node remains in state S in the absorbing state, provided that it has neighbors in state R and was originally in state S. P W ( ) is defined similarly. We note that for a certain node to have neighbors in state R in the absorbing state means that the node receives attempts to infect it when the recovered neighbors are in state I . Thus, a node still remaining in state S with neighbors in state R has to be unchanged from infection attempts through the entire process. Thus, we obtain

Next, the probability P W ( ) is given as

where j denotes the number of attacks that a node sustains before it changes to state W . Using the relation P S ( ) + P W ( ) + P R ( ) = 1, one can determine P R ( ) in terms of P S and P W . The local tree approximation allows us to define q n similarly to q but now at the tree level n. The probability q n+1 can be derived from q n as follows:

where the factor kP d (k)/z is the probability that a node connected to a randomly chosen edge has degree k. As a particular case, when the network is an ER network with

Equation (8) reduces to a self-consistency equation for q for given epidemic parameter values in the limit n → ∞. Once we obtain the solution of q, we can obtain the outbreak size m(κ) using Eq. (5) . For ER networks, however, m(κ) becomes equivalent to q, thus the solution of the self-consistency equation (8) yields the order parameter. We remark that the methodology we used here is similar to those used in previous studies of epidemic spreading on complex networks [6, 10, 15, 20, 21] . Hereafter, we set μ = κ for convenience and define a function:

Using formula (9), we approximate F (m,ρ 0 ) in the limit m → 0 as

where

For convenience, we neglect the higher-order terms and redefine F (m,ρ 0 ) as

Depending on the relative magnitudes of a and b, various solutions of the self-consistency equation F (m,ρ 0 ) = 0 can be obtained. However, we need to check whether these solutions are indeed physically relevant in the steady state when we start epidemic dynamics from a certain initial condition. The stability criterion was established in a previous work [22] : The solution F (m 0 ,ρ 0 ) = 0 is stable if and only if ∂F (m,ρ 0 )/∂m| m=m 0 < 0.

The equation of state in the steady state can be obtained using F (m,ρ 0 ) = 0. We notice that F (m = 0,ρ 0 ) = ρ 0 /(1 − ρ 0 ) > 0 and F (m = ∞,ρ 0 ) = −∞ because c < 0, as shown in Fig. 2 . We examine the solutions of dF (m,ρ 0 )/dm = 0, which are obtained as where a, b, and c are given in formulas (12)- (14) . Note that a depends on ρ 0 . At these extreme points m ± , F (m,ρ 0 ) has either a local maximum or a local minimum. For a given ρ 0 , z, and η, both m ± values exist, and they are positive in the range of κ d < κ < κ a , where b 2 /9c 2 − a/3c = 0 at κ = κ d , and a = 0 at κ = κ a . For a given z and η, diverse types of phase transitions occur depending on ρ 0 . When ρ 0 is less than a certain value ρ c , the order parameter jumps at a transition point. On the other hand, when ρ 0 ρ c , the order parameter increases continuously with κ. At ρ 0 = ρ c , m + = m − = m 0 and F (m 0 ,ρ c ) = 0 at κ = κ d = κ c , as schematically shown in the blue (lower) curve in Fig. 3 .

When ρ 0 < ρ c , there exists a range of κ in which F (m,ρ 0 ) = 0 has more than one solution, as shown in Fig. 2 . The order parameter m versus κ is shown in Fig. 4 (a) and 4(b). In particular, when κ has a certain value κ − c , m − obtained using Eq. (16) satisfies F (m − ,ρ 0 ) = 0. The m − value at κ − c is denoted as m − 0 . We also define κ + c and m + 0 similarly to m + in Eq. (16) . We note that κ + c < κ − c . Depending on the magnitude of the reaction probability κ relative to κ + c and κ − c , the order parameter behaves differently, as follows:

(1) For κ < κ + c , there exists one stable solution m = m d (κ), which increases slowly with κ. It is obtained that to the same one, which is denoted as m 0 . The function F (m,ρ c ) versus m is shown in Fig. 3 , and the order parameter m versus κ is shown with the analytic solution and simulation data in Fig. 5 . At κ c , singular behavior occurs, and the order parameter m behaves as m − m 0 ∼ |κ − κ c | 1/3 on both sides. The derivation of this exponent 1/3 is presented in the Appendix.

To estimate various critical exponents, we perform extensive numerical simulations on ER networks with mean degree z = 8. For simplicity, the reaction probability μ is set equal to κ, and η = 1/2. With these parameter values, we numerically solve Eq. (8) and determine ρ c ≈ 0.00747762 and κ c ≈ 0.108021. Specifically, we increase the value ρ 0 from zero until the values m − 0 and m + 0 become the same. Then ρ c is determined. Moreover, at that point, the values κ − c and κ + c also become the same as κ c ≈ 0.108021. The values ρ c and κ c will be used in numerical analysis later. For ρ 0 < ρ c , we take ρ 0 = 2 × 10 −3 in the simulations. We take the average over 50 different dynamics samples for each of 1600-4000 network configurations. Thus, 80 000-200 000 configuration averages were taken to obtain each data point.

Analytically we obtained that the order parameter behaves as m − 0 − m d (κ) ∼ ( κ) β with β = 1/2 in the thermodynamic limit, where κ ≡ κ − c − κ. The dashed curve in the inset of Fig. 6 are obtained from the analytical solution Eq. (8) by taking the limit n → ∞. This means that the solution is relevant in the thermodynamic limit. This dashed curve follows a solid red line in the region κ < κ * ≈ 10 −4.2 (this point was estimated and indicated by an arrow), whereas it is deviated from the red line for κ > κ * . This means that the critical behavior m − 0 − m d ∼ ( κ) 1/2 appears in the region κ < κ * . On the other hand, the data points in the main panel were obtained by simulations for different system sizes. The dashed curve was the same drawn in the inset. A particular point on the dashed curve at κ * is indicated by an arrow, which locates at the same κ * in the inset. As we may see, the data points obtained from N = 1.6384 × 10 8 begin to deviate from the dashed curve at this κ * . This means that the data points obtained from the system of size N = 1.6384 × 10 8 for κ < κ * belong to the critical region, whereas other data points do to the noncritical region. This means that finite-size scaling behavior of the order parameter that appears in the form of plateaus in Fig. 6 has to be checked with the data points for the systems of sizes N (10 8 ). However, it would be impractical to perform simulations with such huge system sizes.

Following the conventional finite-size scaling theory,

at κ − c . We check this relation in Fig. 7 . For small system sizes N , β/ν seems to be about 0.2, whereas it is estimated to be ≈ 0.24 for large N . Again the crossover occurs between the system sizes N = 10 7 and 10 8 . We could obtain a more reliable value for the exponent ratio β/ν for somewhat larger system sizes, but that is impractical.

The fluctuation of the order parameter

For finite systems of size N , it is expected that χ ∼ N γ /ν at κ = κ − c . From the simulation data, we obtain γ /ν ≈ 0.5, as shown in Fig. 8 .

With the measured values β/ν ≈ 0.24 and γ /ν ≈ 0.5 and the analytic result β = 1/2, we guessν = 2 and then γ = 1. If we use those values, then the hyperscaling relation 2β + γ =ν would hold.

c in a double logarithmic scale. Data are obtained for ER networks with mean degree z = 8. ρ 0 = 2 × 10 −3 is used. The slope of the data point corresponds to β/ν. As the system size is increased, crossover behavior appears in the slope from −0.2 to −0.24, which indicates that β/ν ≈ 0.24 in the limit N → ∞.

At ρ 0 = ρ c , the jump in the order parameter does not appear, and m + 0 = m − 0 at κ = κ c in the thermodynamic limit. In finite systems, however, the order parameter can still exhibit a jump in some samples. Thus, the order parameter distribution p(m) accumulated over different samples exhibits two separate peaks, as shown in Fig. 9 . We regard the data points of We numerically obtain γ /ν ≈ 0.69 (Fig. 12) and γ /ν ≈ 0.6 ( Fig. 13) .

On the basis of the numerically obtained values β/ν ≈ 0.53 and γ /ν ≈ 0.69, and the theoretical value β = 1/3, we estimateν ≈ 2.179 and γ ≈ 1.5. Those values are confirmed in Fig. 14 for χN −γ /ν versus (κ c − κ)N 1/ν , in which the data collapse well with the choices ofν ≈ 2.13 ± 0.1 and γ /ν ≈ 0.69. The measured values of the exponents satisfy the hyperscaling relation (2β + γ )/ν ≈ 0.996 well. Similarly, for κ > κ c , on the basis of the numerical values β /ν ≈ 0.164 and γ /ν ≈ 0.6, and the theoretical value β = 1/3, we obtain ν ≈ 2.03 and γ ≈ 1.218. Data for χ for different system sizes collapse well into a single curve with the choices of ν = 2.13 ± 0.1 and γ /ν = 0.6 (Fig. 15) . These values yield (2β + γ )/ν ≈ 0.91 − 0.93, which deviates slightly from the expected value of unity that would satisfy the hyperscaling relation. To obtain those results, we used the numerical values ρ c ≈ 0.00747762 and κ c ≈ 0.108021. We remark that β = β = 1/3 is obtained analytically. FIG. 13 . Plot of the susceptibility χ , the fluctuation of the order parameter m u , as a function of the system size N at κ c ≈ 0.108021. γ /ν is estimated to be ≈ 0.6 ± 0.005. Data are obtained for ER networks with mean degree z = 8. ρ 0 = ρ c is taken as ≈ 0.00747762.

FIG. 14. Scaling plot of the susceptibility χ for κ < κ c in the form χN −γ /ν versus (κ − κ c )N 1/ν , whereν ≈ 2.13 and γ ≈ 1.47 are used. Data are obtained for ER networks with mean degree z = 8. ρ 0 = ρ c is taken as ≈ 0.00747762.

We investigated the properties of phase transitions in the SWIR model with a finite density ρ 0 of initially infected seeds [7] . A node in the state S can change its state to weakened (W ) or infected (I ) when it comes in contact with an infected node from the same or a different root. A weakened node can also change its state to infected (I ) when it contacts an infected node from the same or a different root. The reaction probabilities κ and μ in Eqs. (1) and (2), respectively, serve as control parameters. For convenience, we take κ = μ. We found that for a given network, there exists a critical density of seeds ρ c such that for ρ 0 < ρ c , the order parameter, the density of nodes in state R in the absorbing state, increases continuously with the critical exponent β = 1/2 as κ is increased up to a transition point κ − c and then jumps to a finite value, followed by a continuous increase. Accordingly, the order parameter behaves as m(κ) = m − 0 − b(κ − c − κ) 1/2 for κ < κ − c , where b is a positive constant. At κ − c , the order parameter is discontinuous FIG. 15 . Scaling plot of the susceptibility χ for κ > κ c in the form χ N −γ /ν versus (κ − κ c )N 1/ν , whereν ≈ 2.13 and γ ≈ 1.28 are used. Data are obtained for ER networks with mean degree z = 8. ρ 0 = ρ c is taken as ≈ 0.00747762. by m = m u (κ − c ) − m − 0 . Thus, the order parameter itself exhibits a hybrid phase transition. This pattern is different from that for the single-seed case, in which the order parameter jumps from m = 0 to a finite value, and thus β = 0. The fluctuation of the order parameter diverges at the transition point κ − c according to a power law with the exponent γ . For the correlation size exponentν measured in finite systems, we find that the hyperscaling relation 2β + γ =ν holds reasonably well.

As ρ 0 is increased, the jump shrinks and becomes zero at ρ c . For ρ 0 = ρ c , the transition becomes continuous. We determined a complete set of critical exponents describing the phase transition at κ c . The critical exponents are listed in Table I . where only nonzero terms are considered. Since δκ and (δm) 3 are two lowest terms in Eq. (A2) and their coefficients have the opposite sign to each other, δm ∼ (δκ) 1/3 when δκ 1. Thus for both cases of κ < (>)κ c , the critical exponents β = β = 1/3.

Proc. Natl. Acad. Sci. USA