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.
Nonequilibrium dynamic transitions driven by cascade dynamics on complex networks have attracted considerable attention recently [1–3]. The spreading of epidemic disease on complex networks [4–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 . 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 ), infected (), or recovered (). Initially, all the nodes are in state except for one seed node in state . The contagion process starts from a single node in state . Each node in state transmits pathogens to its neighbors in state and infects each of them with probability ; then it changes its state to 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 after the system falls into the absorbing state becomes ; i.e., the system falls into a subcritical (supercritical) state. In between, an epidemic transition occurs at , 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 , the SIR model undergoes a continuous percolation transition following the universal behavior of ordinary percolation.
The SIR model with multiple seeds has been considered , in which two percolation transitions occur successively at and as is increased. The density of nodes in state is finite for , whereas the density of nodes in state disappears for . Thus, there exists a state of coexisting nodes in states and between and .
The SIR model was extended to a two-step contagion model, in which a weakened state () can exist between the and states. Accordingly, this model is called the SWIR model [7,20]. Nodes in state are involved in the reactions and , which occur in addition to the reactions and 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–22]. The order parameter, defined as the density of nodes in state 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 . The dynamic rule of the SWIR model is so simple that its underlying mechanism for the discontinuous behavior of the order parameter was disclosed . Moreover, the mechanism turned out to be universal in other models such as -core percolation [31–34], the cascading failure model on interdependent networks [35–39], and epidemic-related models [5,6,8,10,12–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.  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 often occur even in early time steps, because nodes in states and involved in that reaction can originate from different seeds (see Fig. 1). We note that the number of multiple seeds was taken as . On the contrary, in the single-seed case, such reactions rarely occur until the system reaches a characteristic dynamic step : When dynamic step is less than , the reactions and are dominant, but the number of nodes in still remains as . The contagion spreads in the form of a branching tree. When the dynamics reaches , the branching tree forms long-range loops due to finite-size effect. Once such loops form, the reaction occurs massively, in which the nodes in state were generated in early time steps. Thus, the density of nodes in state increases drastically as many as 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 .
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 nodes. Initially, nodes are selected randomly from among those nodes and assigned to state ; the other nodes are assigned to state . At each time step , the following processes are performed: (1) All the nodes in state 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 , it changes its state in one of the two ways: either to with probability or to with probability . If a neighbor is in the state , it changes to with probability , where , and are the contagion probabilities for the respective reactions. (3) All nodes in the list change their states to . 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 when is less than a critical value , and it shows a continuous transition at when for given parameter values , and , where is the mean degree of a given ER network.
In an absorbing state, each node is in one of three states: the susceptible , weakened , and recovered states. The order parameter , the density of nodes in state in an absorbing state, is written using the local tree approximation as
Similarly to , we define as the conditional probability that a node remains in state in the absorbing state, provided that it has neighbors in state and was originally in state . is defined similarly. We note that for a certain node to have neighbors in state in the absorbing state means that the node receives attempts to infect it when the recovered neighbors are in state . Thus, a node still remaining in state with neighbors in state has to be unchanged from infection attempts through the entire process. Thus, we obtain
The local tree approximation allows us to define similarly to but now at the tree level . The probability can be derived from as follows:
Equation (8) reduces to a self-consistency equation for for given epidemic parameter values in the limit . Once we obtain the solution of , we can obtain the outbreak size using Eq. (5). For ER networks, however, becomes equivalent to , 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:
Depending on the relative magnitudes of and , various solutions of the self-consistency equation 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 : The solution is stable if and only if .
The equation of state in the steady state can be obtained using . We notice that and because , as shown in Fig. 2. We examine the solutions of , which are obtained as
When there exists a range of in which has more than one solution, as shown in Fig. 2. The order parameter versus is shown in Fig. 4(a) and 4(b). In particular, when has a certain value obtained using Eq. (16) satisfies . The value at is denoted as . We also define and similarly to in Eq. (16). We note that . Depending on the magnitude of the reaction probability relative to and , the order parameter behaves differently, as follows:
When , there exists a reaction probability that satisfies the relation , and . Thus, the two solutions, and , reduce to the same one, which is denoted as . The function versus is shown in Fig. 3, and the order parameter versus is shown with the analytic solution and simulation data in Fig. 5. At , singular behavior occurs, and the order parameter behaves as on both sides. The derivation of this exponent is presented in the Appendix.
To estimate various critical exponents, we perform extensive numerical simulations on ER networks with mean degree . For simplicity, the reaction probability is set equal to , and . With these parameter values, we numerically solve Eq. (8) and determine and . Specifically, we increase the value from zero until the values and become the same. Then is determined. Moreover, at that point, the values and also become the same as . The values and will be used in numerical analysis later. For , we take 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 with in the thermodynamic limit, where .
The dashed curve in the inset of Fig. 6 are obtained from the analytical solution Eq. (8) by taking the limit . This means that the solution is relevant in the thermodynamic limit. This dashed curve follows a solid red line in the region (this point was estimated and indicated by an arrow), whereas it is deviated from the red line for . This means that the critical behavior 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 begin to deviate from the dashed curve at this . This means that the data points obtained from the system of size 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 . However, it would be impractical to perform simulations with such huge system sizes.
Following the conventional finite-size scaling theory,
The fluctuation of the order parameter diverges as . For finite systems of size , it is expected that at . From the simulation data, we obtain , as shown in Fig. 8.
With the measured values and and the analytic result , we guess and then . If we use those values, then the hyperscaling relation would hold.
At , the jump in the order parameter does not appear, and at in the thermodynamic limit. In finite systems, however, the order parameter can still exhibit a jump in some samples. Thus, the order parameter distribution accumulated over different samples exhibits two separate peaks, as shown in Fig. 9. We regard the data points of in the region (), where has the theoretical value , as those obtained from for different samples. At , in finite systems, we obtain the power-law behaviors with (Fig. 10) and with (Fig. 11).
For , the fluctuation of the order parameter behaves as with a certain scaling function . On the other hand, for , we obtain that behaves as with a certain scaling function . We numerically obtain (Fig. 12) and (Fig. 13).
On the basis of the numerically obtained values and , and the theoretical value , we estimate and . Those values are confirmed in Fig. 14 for versus , in which the data collapse well with the choices of and . The measured values of the exponents satisfy the hyperscaling relation well. Similarly, for , on the basis of the numerical values and , and the theoretical value , we obtain and . Data for for different system sizes collapse well into a single curve with the choices of and (Fig. 15). These values yield , which deviates slightly from the expected value of unity that would satisfy the hyperscaling relation. To obtain those results, we used the numerical values and . We remark that is obtained analytically.
We investigated the properties of phase transitions in the SWIR model with a finite density of initially infected seeds . A node in the state can change its state to weakened () or infected () 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 () 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 such that for , the order parameter, the density of nodes in state in the absorbing state, increases continuously with the critical exponent as is increased up to a transition point and then jumps to a finite value, followed by a continuous increase. Accordingly, the order parameter behaves as for , where is a positive constant. At , the order parameter is discontinuous by . 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 to a finite value, and thus . The fluctuation of the order parameter diverges at the transition point according to a power law with the exponent . For the correlation size exponent measured in finite systems, we find that the hyperscaling relation holds reasonably well.
As is increased, the jump shrinks and becomes zero at . For , the transition becomes continuous. We determined a complete set of critical exponents describing the phase transition at . The critical exponents are listed in Table I.
This work was supported by the National Research Foundation of Korea by Grant No. NRF-2014R1A3A2069005.
Here we introduce an analytical method to determine the critical exponents at . It was already noted in Sec. III B that for ,