Mean-field theory for double-well systems on degree-heterogeneous networks

Many complex dynamical systems in the real world, including ecological, climate, financial, and power-grid systems, often show critical transitions, or tipping points, in which the system's dynamics suddenly transit into a qualitatively different state. In mathematical models, tipping points happen as a control parameter gradually changes and crosses a certain threshold. Tipping elements in such systems may interact with each other as a network, and understanding the behavior of interacting tipping elements is a challenge because of the high dimensionality originating from the network. Here we develop a degree-based mean-field theory for a prototypical double-well system coupled on a network with the aim of understanding coupled tipping dynamics with a low-dimensional description. The method approximates both the onset of the tipping point and the position of equilibria with a reasonable accuracy. Based on the developed theory and numerical simulations, we also provide evidence for multistage tipping point transitions in networks of double-well systems.

Many tipping point phenomena in nature can be modeled by complex systems that are a collection of tipping elements interacting through a network. Examples include extinction of species in ecosystems in which different species are connected by foodweb, mutualistic, and other ecological networks [2,30], blackouts in power grids [15,23], psychiatric disorders [21,31] and the Greenland ice sheet [36,39]. Obviously, the behavior of the tipping elements is not independent of each other when they are connected as a network [9,67,68]. For instance, suppose that a tipping element, denoted by v, passes the tipping point, i.e., the value of a control parameter of the system at which a tipping point occurs. Then, it is more likely that a second tipping element directly connected with v will tip as well, resulting in a domino effect [35,64].
Accurately describing the tipping points in dynamics on networks is desirable because we may want to anticipate and prevent them as well as to know the magnitude of a tipping cascade in the network. However, high dimensionality of dynamical systems on a network and the complex structure of the network at hand generally make this task difficult. One strategy towards this goal is to reduce the dimension of the original system without loss of too much information about the original dynamics. For example, Gao, Barzel, and Barabási proposed a theory to reduce a class of dynamical system on networks, including networks of dynamical tipping elements, into one-dimensional dynamics [27]. This method, which we refer to as GBB reduction, is applicable to directed and weighted networks and uses each node's weighted in-degree (i.e., the number of incoming edges) and out-degree (i.e., the number of outgoing edges). The GBB reduction has been generalized to the methods which use the eigenvalues and eigenvectors of the adjacency matrices or related matrices of the network [46,54,72]. A related but different two-dimensional reduction locates the tipping point with a high accuracy for ecological dynamics occurring on bipartite networks representing mutualistic interactions [34].
The GBB reduction is a degree-based mean-field (DBMF) theory, also called heterogeneous mean-field theory, in the sense that the theory only uses the information about the degree of each node. In general, a DBMF theory for network dynamics assumes that nodes with the same degree are statistically the same and obey the same dynamical rule. It is expected to be accurate for large uncorrelated random networks with general degree distributions. In fact, the most common family of DBMF theory is different from the GBB reduction in that the former aims at deriving a dynamical equation that any node with a given degree obeys, one equation per degree value.
In contrast, the goal of the GBB reduction and its extensions is to derive a one-dimensional or low-dimensional reduction of the entire dynamics on the network. To avoid confusion, we save the term DBMF theory only to refer to the conventional type of the DBMF theory and do not refer to the GBB reduction as a DBMF theory below. The DBMF theory (in the conventional sense) has been successful in describing percolation [13,17,18,57,59], dynamic epidemic processes [6,10,52,56,[60][61][62][63], synchronization [20,33,41,43], random walks [5,45], voter models [3,71], and rock-scissors-paper games [53], to name a few, on degree-heterogeneous uncorrelated networks.
A similar DBMF theory for networks of tipping elements, if accurate, will be a useful tool for anticipating tipping points as a system parameter gradually changes and for understanding the impact of the degree distribution on the tipping behavior. Thus motivated, in the present study, we develop a DBMF theory for prototypical dynamical systems showing tipping behavior, i.e., the double-well system, connected as networks.

II. MODEL
We consider networks in which each node is a tipping element that obeys deterministic dynamics of a double-well system. When there is no interaction between nodes, the state of each node evolves according toẋ where r 1 , r 2 , and r 3 are constants (with r 1 < r 2 < r 3 ) and u represents an external input applied to the node. This system has a unique stable equilibrium, x * , with x * < r 1 when u < u c,ℓ , where u c,ℓ is schematically shown in Fig. 1. We call this equilibrium the lower equilibrium and denote it by x ℓ . The system has a unique stable equilibrium satisfying x * > r 3 , when u > u c,u , where u c,u is shown in Fig. 1. We call this equilibrium the upper equilibrium and denote it by x u . Both x ℓ and x u are stable when u ∈ (u c,ℓ , u c,u ).
If we initialize the system at the lower equilibrium and gradually increase u, then the system undergoes a saddle-node bifurcation at u = u c,u such that x * jumps from the lower equilibrium to the upper equilibrium. If we initialize the system with the upper equilibrium and gradually decrease u, then the system jumps to the lower equilibrium at u = u c,ℓ , implying a hysteresis. Similar models have been used for representing alternative stable states and hysteresis in ecosystems [7,11,79,81], thermohaline circulation [78], and ice sheets [48].
We consider tipping elements that obey the coupled double-well system dynamics given bẏ where N is the number of nodes, D is the coupling strength, and A = (A ij ) is the N × N adjacency matrix of the network. Equation (2) has been used for modeling ecological systems such as lake chains in which a node represents a lake [36], global interactions among climate tipping elements such as the Greenland ice sheet, the Atlantic Meridional Overturning Circulation (AMOC), and the Amazon rainforest [79], and Amazon as a network of tipping elements in which each node represents forests within a specific area [81] (also see Refs. [11,40,80] for analyses of the same model). For simplicity, we consider undirected and unweighted networks such that A ij = A ji ∈ {0, 1} for any i, j ∈ {1, . . . , N }. We exclude self-loops, i.e., we set A ii = 0 ∀i ∈ {1, . . . , N }. We denote by A ji the degree of the ith node and by p(k) the degree distribution. This model allows us to investigate tipping cascades on networks.

III. DEGREE-BASED MEAN-FIELD THEORY
To describe the dynamics given by Eq.
(2), we develop a DBMF theory, which assumes that the nodes with degree k are statistically equivalent to each other and that the states of different nodes are statistically independent of each other except their dependence on k. This assumption corresponds to a configuration model of networks, i.e., a random network with a given degree sequence or distribution. In this theory, we replace the elements of the adjacency matrix, A ij , by its ensemble average, denoted by A ij , which represents the probability that the ith and jth nodes are adjacent to each other, assuming A ij ≤ 1. For uncorrelated networks, we obtain where k is the average degree. Equation (3) implies that the probability that the jth node is a neighbor of an arbitrary node is equal to k j / (N k ). The DBMF theory based on Eq. (3) was first pioneered for percolation and the susceptible-infectious-susceptible model (see Section I for references).
(2), we obtaiṅ is the degree-weighted average of x j .
For a given Θ value, the equilibria of Eq. (4) are given by the intersection of the cubic polynomial ) and the horizontal line y = Dk i Θ + u in the x-y plane (see Fig. 1). Denote by (x (1) ,ỹ (1) ) the unique local maximum of y = (x − r 1 )(x − r 2 )(x − r 3 ) and by (x (2) ,ỹ (2) ) the unique local minimum.
For the nodes with the smallest degrees satisfying Dk i Θ + u <ỹ (2) , the curve y = (x i − r 1 )(x i − r 2 )(x i − r 3 ) and the line y = Dk i Θ + u have just one intersection at x i <x (1) such that there is a unique stable equilibrium, which we call the lower state (see the magenta line in Fig. 1). For the nodes with the largest degrees satisfying Dk i Θ + u >ỹ (1) , there is a unique stable equilibrium satisfying x i ≥ x (2) , which we call the upper state (see the red line in Fig. 1). Finally, for the nodes whose degree satisfiesỹ (2) < Dk i Θ + u <ỹ (1) , Eq. (4) allows both lower and upper states (see the blue line in Fig. 1).
We assume that all nodes are initially in their respective lower states and that we gradually increase u. A larger value of Dk i Θ + u induces the upper state for the ith node under the DBMF theory. Therefore, at a given value of u, the N q nodes with the smallest degree, i.e., the nodes whose k i is smaller than a threshold valuek, are in the lower state, where q is the fraction of nodes in the lower state. The other N (1 − q) nodes, which have k i ≥k, are in the upper state. The value of k i that satisfies Dk i Θ + u =ỹ (1) is equal tok. In other words, we obtaiñ where Θ * is the Θ value in the equilibrium.
We obtain Θ * andk from Eqs. (6) and (7) as follows. First, the right-hand side of Eq. (7), which we write Θ * 2 (k), is a continuous and monotonically decreasing function ofk, given the Θ * value with which to calculate the is also a continuous and monotonically decreasing function ofk, , which one can obtain by, for example, the bisection method. Although there may be multiple roots of Θ * 1 (k) − Θ * 2 (k) = 0 depending on the degree distribution, with a practical degree distribution, Θ * 1 (k) − Θ * 2 (k) = 0 usually has only one root. We regard that all the nodes are in the upper or lower states if the obtained solution satisfiesk ≤ k min ork > k max , respectively.

IV. NUMERICAL RESULTS
In this section, we mainly investigate the accuracy of the DBMF theory developed in Section III against direct numerical simulations in locating the tipping point and approximating {x 1 , . . . , x N } at equilibrium.

A. Numerical methods and networks
We set r 1 = 1, r 2 = 2, and r 3 = 5. We either fix D and vary u (Section IV B) or fix u and vary D (Section IV C). For the given u and D values, we set the initial conditions to x i = 0.01 (with i = 1, . . . , N ) and run simulations for 30000 time units to ensure that the equilibrium has been reached. At a relatively small value of u or D, all the nodes are in their lower state at equilibrium.
Then, we gradually increase u or D.
We use the following networks. First, we use networks generated by the Barabási-Albert (BA) model [4] that produces scale-free networks with p(k) ∝ k −3 , where ∝ represents "proportional to".
We refer to these networks as BA networks. To generate a BA network, we start with a network with two nodes connected by an edge and connect each of the two edges emanating from each new node to an existing node according to the linear preferential attachment rule. We obtain k ≈ 4, where ≈ represents "approximately equal to".
Second, we use the Holme-Kim network model [32], which is a modification of the BA model for high clustering (i.e., a large number of triangles). We construct networks with average degree k ≈ 10 by setting the number of edges that each new node brings in to five. We set the probability of constructing a triangle for the each new edge to 0.5. Third, we use the largest connected component (LCC) of a coauthorship network of researchers in network science, which has N = 379 nodes and 914 undirected edges [58]. Each node in this network represents a researcher publishing a paper in network science up to year 2006. An edge exists between two nodes if the two researchers coauthored at least one paper. Finally, we use the LCC of the hamsterster social network, which has N = 1788 nodes and 12476 undirected edges [44]. A node in this network represents a user of hamsterster.com. Two users are adjacent if they have a friendship relation on the website.
We compare the performance of the DBMF theory with that of the GBB reduction. The GBB reduction for the coupled double-well system (see Eq. (2)) is given bẏ [27]. When the GBB reduction is accurate, the dynamics and equilibria of x obtained from Eq. (8) closely approximate those of the observable x eff , called the effective state, given by where we obtain {x 1 , . . . , x N } from direct numerical simulations of the original N -dimensional dynamical system, Eq.
(2). To compare the accuracy of the DBMF theory with the GBB, we measure the effective state calculated from the DBMF theory given by where

B. Tipping under gradual increases in u
We set D = 0.001 and vary u in this section. We show the dependence of x eff on u for a BA network with N = 10 2 nodes, a BA network with N = 10 3 nodes, the coauthorship network, and the hamsterster network in Figs. 2(a), 2(b), 2(c), and 2(d), respectively. We find that, in numerical simulations, the tipping point does not occur at the same value of u for the different nodes, resulting in gradual transitions from the lower state to the upper state as one increases u; we will examine this phenomenon in Section IV D . Therefore, we need to define which tipping points we want to anticipate by theory. Because we are usually interested in critical points at which a macroscopic number of nodes (i.e., O(N )) starts to experience a major transition, we operationally define the tipping point as the value of u at which x eff exceeds 1.5, denoted byũ 1.5 (see the dashed lines in Fig. 2), for the first time as we gradually increase u. Roughly speaking, 5-10% of the nodes have switched from the lower state to the upper state at u =ũ 1. 5 . Figures 2(a)-(d) suggest that the DBMF theory better approximatesũ 1.5 than the GBB reduction does. We quantitatively confirm this claim in Table I, which compares the error in the tipping point defined by |u GBB −ũ 1.5 | for GBB reduction and |u DBMF −ũ 1.5 | for our DBMF theory, where u GBB and u DBMF are the tipping point estimated by the GBB reduction and the DBMF theory, respectively. We note that we do not add a generalization of the GBB reduction [46,72] to the present comparison because the DBMF theory and the GBB reduction only use the degree of each node, whereas the generalizations of the GBB reduction use the full adjacency matrix of the network [54].
The accuracy of the different approximations depends on the definition of the tipping point. To show this, consider a tipping point defined as the value of u at which x eff exceeds 2.5 for the first time as one increases u, denoted byũ 2.5 (see the dotted lines in Fig. 2). With this definition, the DBMF theory is better at locating the tipping point than the GBB reduction for two networks and vice versa for the other two networks (see Table I).
We also confirmed that, for BA networks of different sizes, the error in locating the tipping point is systematically smaller for the DBMF theory than the GBB reduction when we useũ 1.5 as the tipping point (see Fig. 2(e)), whereas the opposite is the case withũ 2.5 (see Fig. 2(f)). With bothũ 1.5 andũ 2.5 , the error increases as N increases.
We show the approximation error for the GBB reduction and the DBMF theory for Holme-Kim networks in Figs. 2(g) and 2(h) forũ 1.5 andũ 2.5 , respectively. The results are similar to the case of BA networks (see Figs. 2(e) and 2(f)) despite the presence of high clustering in the Holme-Kim networks. Specifically, the DBMF theory is better approximating the tipping point than the GBB reduction for the tipping pointũ 1.5 (see Fig. 2(g)), whereas the opposite holds true except for small networks in the case ofũ 2.5 (see Fig. 2(h)). TABLE I. Approximation error for the DBMF theory and GBB reduction. We set D = 0.001 and gradually increase u to determine tipping pointsũ 1.5 andũ 2.5 . We set u = 0 and gradually increase D to determine tipping pointsD 1.5 andD 2.5 .ũ We next compared the performance of the DBMF theory and the GBB reduction in approximating the the effective state, x eff , off the tipping point, which has been a main goal of recent papers on dimension reduction techniques for dynamical systems on networks [27,42,46,[72][73][74].
To quantify the accuracy at approximating x eff , we measure the relative error between the results obtained from direct numerical simulations and the theoretical estimate, which we denote by ǫ.

C. Tipping under gradual increases in D
Next, we examine the accuracy of the DBMF theory and GBB reduction when we set u = 0 and vary D. In Fig. 4(a), (b), (c), and (d), we show the relationship between x eff and D for a BA network with N = 10 2 nodes, a BA network with N = 10 3 nodes, the coauthorship network, and the hamsterster network, respectively. For approximating the location of the tipping point, we find that our DBMF theory is more accurate than the GBB reduction when we define the tipping point byD 1.5 , i.e, the value of D at which x eff exceeds 1.5 for the first time when we increase the value of D (see the dashed lines in Fig. 4). If we define the tipping point byD 2.5 (see the dotted lines in Fig. 4), the DBMF theory is more accurate on the BA network with N = 10 2 nodes ( Fig. 4(a)) and the hamsterster network (Fig. 4(d)) and vice versa on the BA network with N = 10 3 nodes ( Fig. 4(b)) and the coauthourship network (Fig. 4(c)); also see Table I.
With BA networks with different numbers of nodes, the DBMF theory locatesD 1.5 more accurately than the GBB reduction (see Fig. 4(e)). The opposite is the case when we define the tipping point byD 2.5 (see Fig. 4(f)) except for small networks.
We also show the error in approximating the tipping point for Holme-Kim networks with different numbers of nodes in Figs. 4(g) and 4(h). The DBMF theory locates the tipping point with a higher accuracy than the GBB reduction with the tipping point being defined byD 1.5 ( Fig. 4(g)), whereas the opposite is the case for the combination ofD 2.5 and large N (Fig. 4(h)). All these results are similar to those when we varied u instead of D (see Section IV B).

D. Multistage transitions
For the given D value, the DBMF theory predicts a range of parameter value u ∈ (u c,ℓ , u c,u ) over which the different nodes tip from the lower to the upper state. Using Eq. (4) and Fig. 1, we obtain u c,ℓ =ỹ (1) − Dk max Θ * ,ℓ and where and We obtain k max ≫ k min in most networks. In contrast, the values of Θ * ,ℓ and Θ * ,u are comparable although Θ * ,ℓ < Θ * ,u . Therefore, except when k max and k min are close to each other, including the case of regular networks (i.e., those in which all nodes have the same degree), one obtains u c,ℓ < u c,u such that there is a nonvanishing range of u in which some nodes are in the lower state and the other nodes are in the upper state. Therefore, the entire N -dimensional dynamical system on the network does not show a single tipping point but undergoes multiple tipping points in uncorrelated degreeheterogeneous networks even when N → ∞. We refer to this phenomenon as multistage transition.
To verify this prediction, we numerically simulated double-well dynamics on several networks.
We selected 17 networks from the archive of empirical networks in the "networkdata" R package [69] that were a mix of human and animal social networks with between 34 and 379 nodes (and mean = 116.2). We simplified the networks by making them undirected and unweighted, retaining only the largest connected component, and allowing no self-loops. When several networks drawn from a published work were available, we used a network with close to 100 nodes. We kept and gradually increase D. A circle represents a network. We gathered networks from [69] which originated in published research [1, 12, 16, 19, 22, 24-26, 28, 29, 51, 58, 65, 70, 76, 77, 83].
increasing the bifurcation parameter (i.e., u or D) until at least 90% of nodes were in the upper state at equilibrium. Then, we calculated the range of the bifurcation parameter over which node transitions occurred as the difference between the first value when at least 10% of nodes were in the upper state and the first value when at least 90% of nodes were in the upper state.
In Fig. 5(a), we show the relationship between the observed range of u and the range of node degrees (i.e., k max −k min ) for each network when we fix D = 0.001 and vary u. Each circle represents a network. Equations (11) and (12) imply that the size of the range of u is roughly proportional to k max − k min , which the results shown in Fig. 5(a) support. The Pearson correlation coefficient between the size of the range of u and k max − k min is 0.87. In Fig. 5(b), we show the relationship between the range of D over which multiple transitions occur and 1/k min − 1/k max when we fix u = 0 and vary D. Note that D ∈ (ỹ (1) − u)/(k max Θ * ,ℓ ), (ỹ (1) − u)/(k min Θ * ,u ) implies that the size of the range of D is roughly proportional to 1/k min − 1/k max . Figure 5(b) is also supportive of the DBMF theory; the Pearson correlation coefficient between the two quantities is 0.59. In both cases, the heterogeneity in the degree distribution in terms of k min and k max predicts the presence of multistage transitions fairly well.

V. DISCUSSION
We have developed a DBMF theory for a standard double-well system on networks. The DBMF theory has turned out to be more accurate than the GBB reduction in locating the tipping point if we define it as the bifurcation parameter value at which a small fraction of nodes has tipped. We have also shown that the DBMF theory is less accurate than the GBB reduction at approximating the effective state substantially above the tipping point.
The original spectral method is a generalization of the GBB reduction and uses the leading eigenvector of the adjacency matrix to linearly combine {x 1 , . . . , x N } and realize a low-dimensional reduction of a family of dynamical systems on networks [46,72]. Similar to the GBB reduction, the original spectral method is accurate at approximating the effective state, whereas it is not accurate at locating the tipping point [42]. We recently extended this spectral method to the case of the nonleading eigenvector that minimizes a mathematically derived error [54]. By design, this spectral method, which we call the modified spectral method, is more accurate than the original spectral method, and presumably than the GBB reduction and DBMF theory, at approximating the effective state except near the tipping point. However, the modified spectral method is far less accurate at locating the tipping point than the original spectral method, and also apparently than the GBB reduction and DBMF theory [54]. Therefore, it seems that the accuracy at locating the tipping point and that at approximating the effective state sufficiently far from the tipping point are in a trade-off relationship. Specifically, among the methods discussed here, the DBMF theory is the most accurate at the former task, the modified spectral method is the most accurate at the latter task, and the GBB and the original spectral method are in between. Looking more closely into this trade-off, including searching a theory that more accurately locates the tipping point than the present DBMF theory, warrants future work.
There exists a critical value of the parameter at which the lower and upper equilibrium states have the same potential energy and thus have the same resilience [8,75]. In this case, there exists the tipping point at which half nodes belong to the lower equilibrium and the other half to the upper equilibrium. With this definition of the tipping point, i.e., u =ũ 2.5 or D =D 2.5 in our case, we have shown that our DBMF theory is not accurate. In the present study, our primary definition of the tipping point is the value of the parameter at which a small fraction of nodes (5-10%) has transitioned from the lower to upper state (i.e., u =ũ 1.5 or D =D 1.5 ). We chose such a small fraction because, in practice, we are often interested in the situation in which a small but noticeable fraction of nodes, rather than half the nodes, has tipped to enter an undesirable state. For example, the tipping point for the flip to non-forest ecosystems in eastern, southern and central Amazonia is estimated to be at 20-25% deforestation [50]. We have shown that our DBMF theory is more accurate than the GBB reduction with this latter definition of tipping.
In the present study, we have assumed that the influence of the jth node on the ith node is proportional to A ij x j , where A ij is the (i, j) element of the adjacency matrix. In contrast, doublewell systems on networks in which the nodes are diffusively coupled have also been studied. Such a model represents, for example, chemical reaction diffusion [37,38,82]. A DBMF theory was previously developed for the case of diffusive coupling [37]. However, adapting their theory to be able to self-consistently determine the mean field (i.e., Θ * in our formulation) and the threshold degree (i.e.,k), and thus further pursuing how nodes with different degrees behave and comparing the results with the case of coupling by the adjacency matrix warrants future work.
We have provided theoretical and numerical evidence of multistage transitions in the present double-well system on networks. In particular, our DBMF theory suggests that the multistage transition is not a finite-size effect because the DBMF theory is more accurate in large uncorrelated networks in general. In fact, our DBMF theory depends on the number of nodes, N , through Eq. (7). However, an increase in N does not mitigate the validity of Eqs. (11) and (12), which suggests that multistage transitions are likely to occur in networks of different sizes. Investigating multistage transitions for a wider variety of networks and dynamical systems, as well as their practical implications, also warrants future work.