Introduction

Percolation is a prototypical model of disorder, which is often used to illustrate the emergence and the resilience of a giant cluster as links between individuals are added and deleted one by one, respectively1. A giant cluster at a transition point in the mean field limit is to good approximation a critical branching (CB) tree with unit mean number of offspring2, 3. The giant cluster of recovered nodes at a transition point of a simple epidemiological model, the so-called susceptible/infective/removed (SIR) model2, one of such percolating clusters grown in the CB processes. Percolation transition is known as a robust continuous transition4.

In a number of systems, however, the situation is more complex: Hybrid percolation transitions (HPTs) occur showing features of both second and first-order phase transitions at a transition point5,6,7,8. In these transitions, the order parameter m(z) exhibits the behaviors simultaneously as

$$m(z)=\{\begin{array}{ll}0 & {\rm{for}}\,z < {z}_{c},\\ {m}_{0}+r{(z-{z}_{c})}^{\beta } & {\rm{for}}\,z\ge {z}_{c},\end{array}$$
(1)

where m 0 and r are constants and β is the critical exponent of the order parameter, and z is a control parameter. Examples include the k-core percolation9,10,11,12, and the cascading failure (CF) model on interdependent networks13,14,15,16. In those systems, as nodes or links are removed one by one above the transition point, the order parameter, the relative size of the giant component decreases continuously, approaches a nonzero value in a critical way at the transition point, where it finally collapses to zero: A HPT occurs. Is there a universal mechanism behind this phenomenon? Can it be formulated in terms of branching processes? Even though these questions are simple and fundamental, there has been no clear answer yet.

Recently we showed on the example of the CF model that there are two kinds of critical phenomena related to the HPT17. One is carried by the behavior of the finite cascades and the other one by the order parameter (the relative size of the giant cluster). We have to distinguish between “finite” and “infinite” avalanches (the latter having the size of the giant cluster). Once an infinite avalanche occurs, the order parameter falls into an absorbing state. Therefore, occurrence of an infinite avalanche is a distinct feature of a HPT, whereas such infinite avalanche is absent for a second-order percolation transition. Thus we need to investigate what happens in the system while an infinite avalanche proceeds.

We recall the results of previous studies on k-core percolation12 and in interdependent networks18 about the temporal evolution of the giant cluster. The order parameter decreases rapidly in the early time regime, exhibits a plateau for a long time in the intermediate time regime, and decreases rapidly in the late time regime. Moreover, it was found that infinite avalanches proceed in the form of a CB process for a long time, followed by a supercritical process. There has been considerable effort to explain the mechanism leading to this scenario for specific models18,19,20. In particular, refs 19 and 20 pointed out the importance of large loops in the creation of the supercritical process for an epidemic model. However, it is still unclear whether there is a universal mechanism, which explains why, how and when such SC processes occur in the late time regime. Here we address these questions and show that there indeed exists such a universal mechanism, which governs the generally observed crossover behavior in a large class of HPT models.

In this paper, we first investigate the mechanism of the crossover behavior from the CB to SC processes using a simple epidemic model with two-step contagion processes21 that exhibits a HPT. After explaining the mechanism of the HPT on this model, we will show that the same mechanism occurs in other models. We will consider k-core percolation, the threshold model, and the CF model on interdependent networks.

Results

Two-step contagion model

We consider the epidemic model introduced in ref. 21, which is a generalization of the so-called susceptible (symbolized as S)-infected (I)-removed (R) (SIR) model by adding a weakened state (W) between susceptible and infected states. This model is referred to as the SWIR model. Various aspects of the model were studied in refs 22,23,24,25,26. Besides the usual reactions S + I → 2I and I → R of the SIR model we have the additional reactions: S + I → W + I and W + I → 2I. The reaction rate from W to I is larger than the rate from S to I. Specifically, we start the dynamics on Erdös-Rényi (ER) random graphs of N nodes with all nodes in state S but one node that is in state I. At time step n, a node in state I (denoted as I n , where subscript represents generation) is selected randomly, the states of all its neighbors are checked one by one. If the state of a neighbor is S, then this state changes either i) to I n+1 with probability κ or ii) to W with probability μ. If the state of a neighbor is W, then the state W changes to I n+1 with probability ν. We repeat the above process for all nodes in state I n and then the state I n changes to R for each associated node. Then all dynamics at time step n are completed and we move to the next time step n + 1. This dynamics continues until the system reaches an absorbing state in which no more infectious nodes remain in the system. The order parameter m(κ) is defined as the fraction of nodes in state R. Under the given reaction probabilities, a HPT occurs if the mean degree \(z > \mathrm{2/(}\sqrt{{\mu }^{2}+4\mu \nu }-\mu )\) and otherwise a continuous transition occurs. This condition is the same as that obtained in ref. 27. The transition point is κ c  = 1/z. The detailed derivations of the transition point and the condition for the HPT are presented in Supplementary Information.

At the transition point κ c , a single infected node can trigger an infinite avalanche of size O(N) with a certain probability P ∞. With the remaining probability, finite avalanches occur and their sizes are o(N). When an infinite avalanche occurs, as shown in Fig. 1(a), the order parameter remains almost zero (o(N)) for long time up to the characteristic time n c (N), beyond which it increases rapidly and reaches its final, O(N) value in a short time period. To see how an infinite avalanche proceeds at a microscopic level, we trace an infection dynamics in the view of branching processes as shown schematically in Fig. 2. At the step n = 0, a single infectious seed is present. At each time step, infected and weakened nodes are generated following the aforementioned rule. Because the probability to generate an infectious node per each edge is 1/z and the mean number of edges outgoing from the infected parent node is z, a single infected node can be generated on average. There is some probability that a weakened node is created. Thus a CB tree of recovered nodes is generated. We notice that, although during the CB process many W nodes are created, there are very few nodes produced from them in state I (consequently nodes in state R) through the reaction W + I → I + I as shown in Fig. 1(b). However, as the dynamics proceeds and approaches n c (N), the reaction W + I → I + I occurs more frequently and the branching ratio to generate a node in state I through this reaction becomes non-negligible.

Figure 1
figure 1

Evolution of the fractions of nodes in each state and of each reaction type. (a) Plot of the fraction of nodes in states R (blue, solid curve), W (green, dashed curve) and I (red, dotted curve) as a function of generation n. Inset: Plot of the maximum slope of the curve R(n) vs N (•) (left vertical axis). The maximum slopes are independent of N. Plot of the characteristic time n c (N) vs N (Δ) (right vertical axis). The fitted straight line has slope 0.35. (b) Plot of the branching ratios as a function of generation n for several types of reactions. Here A (B) represents the mean number of offspring that change their state from S to I (W) by the their parents in state I each reaction S + I → 2I (S + I → W + I). C represents the mean number of offspring that change their state from W to I by the reaction W + I → 2I. For both (a,b), data are obtained from a single realization of an infinite avalanche on an ER network with mean degree z = 8 of system size N = 5.12 × 106 using the coefficients κ = 1/8, μ = 1/16 and ν = 0.9.

Full size image
Figure 2
figure 2

Schematic picture of the epidemic spreading process in the SWIR model. (a) The process begins from an infectious node. (b) It can infect a susceptible node among its neighbors and change the state thereof from S to I and can also change the state of another neighbor from S to W. This type of process persists for a long time and the critical branching tree is constructed. After a long O(N 1/3) time passes, an infectious node can contact a node in state W that was created much earlier and change its state from W to I in (d). In addition, (d) the I-node in (c) infects a susceptible neighbor and changes its state to I. Thus, a SC process occurs, leading to the jump in the order parameter.

Full size image

In order to determine the crossover point between CB and SC we recall that the size of the largest cluster at the critical average degree z c of ordinary percolation of the ER graph is O(N 2/3)28, 29. The giant cluster at criticality has the topology of a CB tree such that the branching process persists up to the steps O(N 1/3) beyond which finite-size effects appear in the form of short-range and long-range loops28 (see also Supplementary Information). In the epidemic models on ER networks the average degree is above the percolation threshold (z > z c , otherwise global spreading would be trivially impossible), however, the reaction probability κ c  = 1/z assures just the critical branching probability by which the infection proceeds. Thus the growing cluster of R nodes can be considered as if a critical ER cluster would develop on the ER supercritical graph. In accordance with this picture the probability distribution of the generation at which a loop is formed in CB processes shows a peak at a characteristic generation n c (N) ~ O(N 1/3) (Fig. 3). This means that long-range loops begin to form mostly when a CB tree is grown up to n c (N). Based on this, we conclude that before n c the I-state nodes are almost entirely generated through the CB tree and the W-state nodes accumulate to an extent of O(N 2/3) because the number of W-state nodes is proportional to that of I-state nodes. Around n c (N) the loops become important, and due to the long range links the reaction W + I → I + I occurs over the entire system, with W-state nodes having been generated at all times to O(N 1/3) (see Supplementary Information). The accumulated population of W-state nodes and the possibility of long-range loop formation lead to an increase of the number of infected offspring above the critical value resulting in the SC process and eventually in the jump of the order parameter. We remark that up to the characteristic generation n c (N), the population of recovered nodes in state R is less than or equal to O(N 2/3), sublinear to the system size O(N), whereas beyond n c (N), the population suddenly increases to O(N). Thus, we regard the characteristic generation n c (N) as the so-called golden time, during which one needs to take some actions to control and prevent a pandemic outbreak. We also remark that the formation of long-range loops needed for a discontinuous percolation transition was first observed and conjectured in a model of two interacting epidemics19, 20. However, the connection between the length scale of long loops and finite-size scaling of the ordinary percolation was missing, so that the scale of golden time could not be predicted.

Figure 3
figure 3

Distributions of loop lengths for different system sizes. Scaling plot of the probability \({P}_{{n}_{I}}\) of the generation n I at which a loop is formed in critical branching processes on ER networks. Data for different system sizes are well collapsed onto a single curve with the scaling form of \({P}_{{n}_{I}}{N}^{\mathrm{1/3}}\) as a function of n I /N 1/3. Data are obtained from ER network with mean degree z = 8 far away from the transition point z c  = 1.

Full size image

k-core percolation

k-core percolation is known as a prototypical model that exhibits a HPT. The k-core subgraph is identified on a graph (here the ER graph with mean degree z) as follows. One starts with removing all nodes that have degree less than k. These removals may decrease the degrees of remaining nodes. If degrees of some nodes become less than k, then those nodes are removed as well. This process is repeated until no more node is removed. For z > z c , a k-core subgraph remains after the pruning process and its size is O(N). The order parameter is the relative size of the k-core subgraph.

Here we remove a randomly chosen node from the k-core subgraph and repeat the pruning process once again. Near z c , this process can remove all nodes (infinite avalanche of size O(N)) or a fraction of nodes (finite avalanche of size o(N)) from the k-core subgraph, each of which contributes to a discontinuous or continuous change of the order parameter in the thermodynamic limit, leading to a first-order or second-order transition, respectively30. As it was shown earlier in ref. 30 the critical exponents of the k-core percolation model are of two kinds: those associated with the order parameter and those with finite avalanches. This is the typical behavior at hybrid percolation transitions induced by cascade dynamics.

We focus on the infinite avalanches at z c from the perspective of branching processes. Let us consider a k-core subgraph configuration at z c in which each node has at least k degree and the deletion of node i leads to the collapse of the entire system. The node i is regarded as an infectious seed node (I). We check the degrees of neighbors of the node i. If a neighbor of i has degree k, it is regarded as a susceptible node (S), and changes its state to I because it will be deleted after node i gets deleted. If a neighbor has degree \(\ell > k\), then it is regarded as a generalized weakened node and denoted as \({W}_{\ell -k}\). Now its state changes to W l−k−1. The subscript \(\ell -k\) refers to the threshold and decreases as the neighbors of that node are deleted. When it becomes zero during an avalanche, the state \({W}_{\ell -k}\) becomes W and the node has the same role as weakened nodes in the SWIR model. This node gets infected when it contacts an I-state node once more. In analogy with the process in the SWIR model the infective state in the n-th generation or branching step is denoted by I n . Once the dynamics in n-th step is completed, the nodes in state I n are deleted. A schematic illustration for a specific example of the avalanche dynamics is presented in Supplementary Information.

Figure 4(a) shows the branching ratio as a function of branching step n for an infinite avalanche of k-core percolation. We find again that the CB process continues up to the characteristic step n c (N) when it changes to the SC process. By the crossover time n c large number of nodes get their degrees reduced to k so that they become W-nodes. The SC process is again driven by the meeting of an old W node with a new I node, W + I → I + I. Such a reaction sets up the rapid SC process and the entire collapse of the k-core subgraph.

Figure 4
figure 4

Evolution of each reaction type for several models. (a) For k-core percolation with k = 3, plot of the branching ratios as a function of generation n for each type of reactions during an infinite avalanche. A represents the ratio of removed nodes with original degree z = 3 (z > 3), which corresponds to the reaction S + I → 2I for the SWIR model. B does the ratio of the reaction W k−l>0 + I → W + I, which corresponds to the ratio of generating weakened nodes denoted by W. C does the ratio of nodes changing their degrees to z = k = 3 (W + I → 2I). A + C represents the total branching ratio of I. Data are obtained from a network with N = 5.12 × 106 at a transition point. (b) For the threshold model, a similar plot. A and B represent the mean number of I and W offsprings generated by the reactions S + I → 2I and \({W}_{\lfloor {k}_{i}{q}_{i}\rfloor -{m}_{i}}+I\to W+I\), respectively. C does the ratio of nodes of I by the reaction corresponding to W + I → 2I. The sum of the mean I offspring from A and C represents the total branching ratio of I. Data are obtained from a network with N = 5.12 × 106 at a transition point. (c) For the CF model, a similar plot. A represents the ratio for S + I → 2I. B is the mean number of new W-state nodes (k A or k B becomes unity for the first time). C is the mean number of I-state nodes transformed from W nodes (W + I → 2I). A + C represents the total branching ratio of I. Data are obtained from a network with N = 5.12 × 106 at a transition point.

Full size image

The threshold model

The threshold model was introduced in ref. 31 for understanding the spread of fads, cultural traits, the diffusion of norms, and innovations, on social networks. In this model, each node i is assigned its threshold value q i and exists in one of two states, either active or inactive state. An inactive node i with m i active neighbors among k i total neighbors (degree) becomes active when its fraction of active neighbors, m i /k i exceeds its threshold value q i . This threshold model is known to exhibit a hybrid phase transition when mean degree z is sufficiently large. Here we show that the mechanism underlying this hybrid phase transition is the same as we observed in the previous instances.

To illustrate how the universal mechanism works in the threshold model, we reconsider the rule of the threshold model in the perspective of the SWIR model in the following way: We match up active nodes in the threshold model with either infectious I or recovered R nodes in the SWIR model. Among the active nodes, an I-state node is the node that becomes active at the preceding step. The other active nodes are regarded as R-state nodes. Inactive nodes are matched up with either susceptible S or weakened W nodes in the SWIR model: (i) A node satisfying k i q i  < 1 from the beginning is regarded as susceptible node. (ii) A node satisfying k i q i  > 1 is regarded as a generalized weakened node and denoted as \({W}_{\lfloor {k}_{i}{q}_{i}\rfloor -{m}_{i}}\) similarly to the k-core percolation case. Then the dynamics proceeds following the same way as in the k-core percolation.

We performed simulations with a single threshold value m i  = 0.16 for all nodes on ER networks with mean degree z c  = 7.47707 which is the transition point for the given threshold value. At this point, the cascade dynamics becomes critical, so that the avalanche size distribution follows a power law.

We obtain the branching ratios as a function of dynamic step (generation) n for several types of reactions for the threshold model, which is shown in Fig. 4(b). Here we also find a crossover from a CB to SC process similar to that of the SWIR model. Again the accumulation of a sufficient number of weakened nodes during the CB process and their activations through long-range loops are the underlying mechanism of the SC behavior.

The CF model on interdependent networks

We consider here ER interdependent networks in the single layer representation of ref. 32. In this picture we have a single ER graph but with two types of links (A and B), for each having the average degree z. The order parameter is the relative size of the giant mutually connected cluster (GMCC), in which every pair of nodes are connected following each type of links. The CF model exhibits a HPT at the transition point z c 15, 17.

As with k-core percolation, the removal of a node from the GMCC can induce further removal of nodes from the GMCC. This avalanche can be infinite or finite, each of which contributes to the discontinuity of the order parameter or the critical behavior of the HPT, respectively17. Here we focus on the infinite avalanches at z c which leads to the collapse of the entire GMCC.

We consider the avalanche process in the view of a branching process of removed nodes15. To describe the avalanche process in terms of the SWIR model, we determine the effective degrees k A (j) and k B (j) of each node j for each type of links. The effective degree k A (j) (k B (j)) is defined as the number of A-type (B-type) of links of the node j following which one can reach O(N) nodes. Each node in the GMCC has k A  ≥ 1 and k B  ≥ 1. We explain how to determine the effective degrees of each node in simulations in Supplementary Information.

The cascading dynamics proceeds in the following way: An avalanche is initiated by removing a node chosen randomly from the GMCC. During an avalanche, we identify removed nodes at each time step, then the effective degrees of the neighbors may decrease. As a result one or both type of the effective degrees of some neighbors can become zero. Then, they are removed from the GMCC at the next time step, i.e., they are infected and removed at the next time step. Such avalanche process propagates to all neighbors of those infected nodes recursively until no more node is removed.

If a node is removed at a time step, it is regarded as an I-state node and it becomes R-state node at the next time step. If one or both effective degrees of a node is unity from the beginning, the node is regarded as a S-state node because it can be infected (i.e., removed) by contacting an infected node (i.e., losing the unit effective degree). If one or both effective degrees of a node become unity during an avalanche, we identify the state of that node as W because the node became vulnerable as a result of contacting infected nodes. This view enables us to understand the correspondence between the cascading dynamics of the CF model and the dynamics of the SWIR model. A specific example of the avalanche dynamics is presented in Supplementary Information.

Figure 4(c) shows the branching ratio as a function of branching steps n for an infinite avalanche of the CF model. We find that a CB process persists and the generating ratio of the weakened nodes is constant with some fluctuations. The number of infected offspring from weakened nodes is negligible up to the characteristic step n c (N), beyond which it increases rapidly. Thus the generation ratio of infected offspring exceeds unity beyond n c (N): a collapse of the giant MCC takes place.

Discussion

We disclosed the universal mechanism of the HPT induced by cascade dynamics on ER networks. We have shown that during the CB processes, the order parameter sustains up to the time step of O(N 1/3), and W-state nodes accumulate to an extent of O(N 2/3). In the SWIR model the nodes are in W state if they have already got into contact with infected node but only partial infection (or weakening of the immune system) took place. The corresponding (generalized) weakened state in k-core percolation is that of a node, which has already lost some neighbors but the number of alive neighbors is still above k. In the CF model the nodes with reduced effective degrees (but larger than zero) correspond to the generalized weakened nodes. Finally, for the threshold model the nodes having infected (active) neighbors but less than required by the threshold criterion are the weakened ones. Those nodes may be thought of as powder keg in explosive percolation. After the CB processes, those W-state nodes change their state to state I in the way of a SC process. Such reactions are achieved along long-range loops of length O(N 1/3) presented in finite systems. As a consequence, infected nodes are generated abundantly in a short time, leading to a discontinuity of the order parameter. This explains that the SWIR model and k-core percolation do not exhibit discontinuous transitions in low dimensional Euclidean space because of the absence of long-range connections in such space. Moreover, we showed that the mechanism is universal for diverse systems such as the multi-stage contagion models including the SWIR model and the threshold model, k-core percolation and the CF model on the interdependent networks. We expect more models also to belong to this category. Finally, we regarded the characteristic generation n c (N) ~ O(N 1/3) as the golden time during which one can control a pandemic outbreak of macroscopic disaster, using for instance the explosive percolation idea33. Because, for n < n c (N), the number of damaged nodes is sublinear as O(N 2/3) to the system O(N), while for n > n c (N), it is linear as O(N).