From spherical chicken to the tipping points of a complex network
FFrom spherical chicken to the tipping points of a complex network
Gui-Yuan Shi,
1, 2, ∗ Rui-Jie Wu, Yi-Xiu Kong, and Kim Sneppen Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark Department of Physics, University of Fribourg, Fribourg, Switzerland (Dated: February 25, 2020)The outbreak of epidemics, the emergence of the financial crisis, the collapse of ecosystem, andthe explosive spreading of rumors, we face many challenges in today’s world. These real-worldproblems can be abstracted into a sequential break down of a complex system in an increasinglystressful environment. Because both the system and the environment require a large number ofparameters to describe, the break down conditions has been difficult to estimate. We use a highlysymmetric system to gauge a complex environment, which enables us to propose a scalar benchmarkto describe the environment. This allows us to prove that all the tipping points of a complex networkfall between the maximum k-core and maximum eigenvalue of the network.
We start with the story about a physicist and a spher-ical chicken. When winter comes, a farmer’s chickens allget sick and the farmer does not know what is wrong withthem. The farmer calls his neighbour, a physicist, to seeif he can figure out what is wrong. The physicist looksat the chickens and then starts scribbling in a notebook.Finally, after several gruesome calculations, he exclaims,”I have got it! But it only works for spherical chickens.”Let us consider his model seriously. The modelassumes a chicken would feel too cold because ofthe combined effects of air temperature( t a ), relativehumidity( h r ), wind speed( v w ) and sunlight intensity( I s ).These spherical chickens are put into an isotropic envi-ronment, and the radius of the spherical chickens is theonly parameter to describe them. It is not hard to imag-ine that the larger the radius, the more cold-resistant thechickens are. FIG. 1. Given an environment, the activity of a sphericalchicken increases with its radius.
Therefore, for a given environment, we find the small-est spherical chicken that can endure the coldness, itsradius can be expressed as a function of the environ-ment parameters ( t a , h r , v w , I s ). We use the critical ra- ∗ [email protected] dius R ∗ ( t a , h r , v w , I s ) to represent the effective coldnessof the given environment. If we draw contour lines for R ∗ in the 4-D parameter space ( t a , h r , v w , I s ), one can imag-ine that a real chicken will feel differently along a contourline. Because the real chicken might be more windproofthan a spherical chicken, but worse at absorbing sunlight.Thus, for a given real chicken, we want to figure outthe environmental regime that makes the chicken feel notcold. For the boundary of the regime, there exists a low-est effective coldness R l and the highest R h among thecontour lines for R ∗ ( t a , h r , v w , I s ). That is to say, if weput the chicken in an environment that effective cold-ness R ∗ that is higher than R h , for sure it will feel toocold; and if R ∗ < R l , the chicken will never feel cold. Inpractice the upper and lower bounds of effective coldnessendurance, R h and R l for any given chicken are impor-tant to know.Now we go back to the study of complex network, de-scribed by following set of differential equations dx i dt = F ( x i , N (cid:88) j =1 A ij G ( x i , x j )) , (1)here x i is the abundance of component i , the weightedadjacency matrix A ij ≥ j to i . The first term of F indicates the self-decayof i in isolation, and the second term describes the totalsupport from its direct neighbors in the network. Thatis, following previous works [1, 2] here we also focus oncooperative systems where each part of the system is con-tributing positively to existence of other parts (both F and G is increasing as function of x j , see method for adetailed analysis for properties of the equations we con-sidered).With an appropriate choice of F and G , Eq. 1 can beused to describe numerous systems include mutualisticecosystem, disease spreading, information propagation,collaboration network and so on. We can imagine thata difficult environment would cause the decay to be toofast, or the support to be inefficient, then the whole sys-tem can only stay in zero abundance state. However,increasing the connectivity or coupling strength in the a r X i v : . [ phy s i c s . s o c - ph ] F e b a b c FIG. 2. An illustration of λ ∗ and λ ∗ c . The vertical axis (cid:104) x ∗ (cid:105) shows the average abundance of components, indicating whetheror not the system collapses. a ) For homogeneous networks, λ = (cid:80) j A ij . Given the functional form F and set of environmentalparameters α , λ increases as the density of the network increases and the system is more likely to have a non-zero fixed-pointsolution. For a specific F and α , the minimum λ that allows the existence of the non-zero solution is marked as λ ∗ . b ) For E.coli network A , if we fix the network and tune the environmental parameters so that it is increasingly difficult for the system tosurvive(increasing λ ∗ ), and our theory predicts all the tipping points will emerge between λ ∗ min = k max ( A ) and λ ∗ max = ρ ( A ). c ) We condense the E.coli network by contracting 6 strongly connected subgraph of the original graph into 6 giant nodes.During the evolution the red nodes go extinct at λ ∗ = 1, orange nodes go extinct at λ ∗ = 1 .
58, and green nodes go extinctat λ ∗ = 2 .
24. Because the nodes support each other within a strongly connected subgraph, and can also contribute to thedownstream nodes but not to upstream nodes. complex network A may allow the system to work (andour chicken to lay eggs).Scientists have made efforts in finding the tippingpoints from many approaches, such as using mean-fieldapproximation[3, 4], linear stability analysis[5], second-order mean-field under quasi linear assumption [1], andlogistic approximation [2].Here we will use a framework inspired by our spheri-cal chickens. First we use a homogeneous network withconstant in-degree (cid:80) j A ij = λ (a spherical chicken withradius R ) to gauge the effective difficulty level λ ∗ (aseffective coldness R ∗ ) of an environment. This is donethrough the 1D equation (simply replace x i , x j with x and replace (cid:80) j A ij with λ ) F ( x, λG ( x, x )) = 0 , (2)obtained from Eq. 1 when all nodes have symmetricalinput. The effective difficulty level λ ∗ is the minimal λ that enable x to have positive solution (see Fig. 2 a ).Given a network A , each of its tipping points has acorresponding λ ∗ , denoted as λ ∗ c . We prove (see methods)that for any given network A , the lower and upper boundsof λ ∗ c (see Fig. 2 b ) have to be limited by properties of thenetwork A : k max ( A ) ≤ λ ∗ c ≤ ρ ( A ) . (3)Here k max ( A ) is an extension of the traditional maxi-mum k -core to the weighted networks [6, 7]. It is definedas the maximum number that allows the existence of asubgraph, in which each node receive at least k max ( A )weighted incoming edges in total within the subgraph;and ρ ( A ) is the spectral radius, i.e. the largest eigen- value, of A . Note that the largest eigenvalue ρ ( A ) isalways larger than k max ( A ) [8].Further, if the interaction term G ( x i , x j ) is a stepfunction of x j , the tipping points always follows λ ∗ = k max ( A ). In contrast, when Eq. 2’s solution from zero tononzero undergoes a continuous transition, all the tran-sition points satisfy λ ∗ = ρ ( A ).In Fig. 2 c we show the abstracted graph of the E.coli network. With this illustration we can explain the”drops” at λ ∗ = 1 and λ ∗ = 1 .
58 as shown in Fig. 2 b .These two drops are also ”tipping points” and can bestudied by our theory. But for simplicity, we only ana-lyze the tipping points that result in a total collapse ofthe network. In the following, we give a few examples toillustrate our theory. Gene regulatory networks.
First let us consider theexample of gene regulatory networks. As in the analysisof Ref [1], we wrongly assume that all inputs are positive,i.e. that all are activation regulations of type: dx i dt = − Bx i + N (cid:88) j =1 A ij x hj x hj . (4)Here the first term determines the degradation, whereasthe parameter h is the Hill coefficient that quantifies co-operativity of the gene regulation. The corresponding 1Dequation of Eq. 4 is: − Bx + λ x h x h = 0 . (5)Obviously when h < h ≥
1, the effective difficulty level λ ∗ is (note 0 = a h c * b h c * FIG. 3.
The tipping points of two gene regulatory networks
Here h = 1 , , , , ∞ . The two dashed lines represent k max (left) and ρ (right) respectively. The inset shows the changes of λ ∗ c depending on the parameter h , and the grey area showsthe predicted bounds of λ ∗ c . We recursively pruned the nodes whose in-degree or out-degree equals zero to run the simulationfaster, because these nodes do not affect the tipping points. We use β origin eff to denote the β eff of the network before pruning. a ) S. cerevisiae [9]. k max = 2, ρ = 3 , β origin eff = 3 ,
43 (4441 nodes), β eff = 3 .
30 (60 nodes). b ) E. coli [10]. k max = 2, ρ = 2 . β origin eff = 1 .
17 (1550 nodes), β eff = 2 .
21 (15 nodes). Data from https://github.com/jianxigao/NuRsE λ ∗ = Bh ( h − h − . (6) h is fixed by the details of the biological regulation,whereas B changes dependent on external stress/livingconditions. For each B , we can calculate the effectivedifficulty λ ∗ , and have the corresponding average abun-dance (cid:104) x ∗ (cid:105) through numerical simulation.(Fig. 3) At atipping point λ ∗ = λ ∗ c , the corresponding (cid:104) x (cid:105) collapse tozero.We predict that λ ∗ c will change from ρ ( A ) (when h = 1)to k max ( A ) (when h → + ∞ ). In comparison, Gao et al.suggest that λ ∗ c = β eff ≡ (cid:104) s out s in (cid:105) / (cid:104) s (cid:105) , independent of h . Here s is degree and (cid:104) s (cid:105) = (cid:104) s out (cid:105) = (cid:104) s in (cid:105) . In manysituations, β eff ≈ ρ ( A ), see Ref. [11, 12]. But in certaincases, they are quite different. Imagine two extreme case:(i) if we add a new node whose in-degree are very highbut out-degree equals zero, then β eff will decrease signifi-cantly while the largest eigenvalue and the tipping pointswill not change. Small β eff indicates a overestimation ofcollapse risk. E. coli network is an example(Fig. 3 b );(ii)for star network, β eff = N/ ρ = √ N −
1. In thiscase β eff underestimates the risk. Later we will see thelatter case can explain phenomena observed in mutualis-tic networks.Our simulations on two networks S. cerevisiae and
E.coli are shown in Fig. 3. The figure illustrates tip-ping points for different h values, with G varying fromfirst order to a step function. It is easily seen that all tip-ping points fall within our predicted region, and that thetipping point changes with onset of non-linearity, thusdiffering from the quasi linear estimates of Ref.[1].In Fig. 3, both networks collapse at λ ∗ = 1 when h → ∞ . At this limit nodes with abundance x ∗ < b , the collapses at 1 and near 1 . c ). Mutualistic ecosystems.
Morone et al. [2] studieda mutualistic ecosystem governed by: dx i dt = − x i d − x i s + γ (cid:80) Nj =1 A ij x j α + (cid:80) Nj =1 A ij x j x i (7)Here d > s > α > γ > x/ ( α + x ) equals 0when x < α and 1 when x > α , proposed the thresholdon the mutualistic benefit K γ = αs ( γ + d )( γ − d ) , and predictedthat all tipping points satisfy K γ = k max ( A ).Here we study the 1D equation − xd − x s + γλx α + λx = 0 , (8)which has a critical minimal λ ∗ that allows for positive xλ ∗ = αs ( √ γ − √ d ) . (9)This λ ∗ is coincident with K γ when d = 0. In thiscase, Morone et al [2] predicted the tipping points satisfy k max ( A ) = αs/γ . Our methodology further teaches us(see methods) that for d = 0, the transition is continu-ous around ρ ( A ) = αs/γ .Using a mutualistic network we illustrate our result inFig. 4 a . We fix s = 1 and a = 1, and plot λ ∗ c as functionof parameter d from 0 to 1. We see that λ ∗ c are between k max ( A ) and ρ ( A ), just as we predicted. As d increases, d a d b FIG. 4.
Critical λ of a real mutualistic network and a random network. We fix s = 1, a = 1. For each d , we find thecritical γ , and then calculate λ ∗ c by Eq. 9. The insets show the structure of the two networks. Red nodes represent the birds andgreen nodes are the plants. The size of a node shows the abundance of the corresponding species near tipping point. The solidblack line shows our simulation result of λ ∗ c relative to the parameter d , the grey area indicates the predicted bounds of λ ∗ c , thedashed blue line shows the prediction β eff in Ref. [1], and the dotted red line shows the prediction in Ref. [2]. a b )The random network generated by randomly swapping the links of network in a . We can find that β eff roughly equals tothe upper bound ρ ( A ), which suggests Gao et al’s [1] results are valid for random networks when transition is continuous( λ ∗ c = ρ ( A )). λ ∗ c will decrease starting from ρ ( A ). In contrast Moroneet al.’s results lead to λ ∗ c increase starting from k max ( A )to 2 k max ( A ). This means for small d , they overestimatethe risk of ecosystem collapse; and when d is large, theyunderestimate the risk. The systematic biases are alsoobserved in their simulation (Fig. 2g in Ref. [2]). Also, β eff from Ref. [1] often lies beyond the upper bound of thecritical area, indicate an underestimation of the collapserisk. This is due to the existence of a large hub thatdistorts their result(recall the star network we discussedabove).In addition, the largest eigenvalue ρ ( A ) itself is alsoa measure of nestedness [14, 15], i.e. a highly nestednetwork tends to have a large ρ ( A ). As λ ∗ c is highlycorrelated with ρ ( A ), the nested networks will have larger λ ∗ c and survive easier at higher stress. This is consistentwith the observation [16–18] that highly nested networksare more robust.Further by comparing Fig. 4 a and b , we can see thatthe real network has a larger λ ∗ c than random network,suggesting that the real networks are more robust thanrandom networks. Epidemic process.
Finally we study the SIS model in epidemic process [19], dx i dt = − x i + β (1 − x i ) N (cid:88) j =1 A ij x j (10)Here β > − x + λβ (1 − x ) x = 0, andthe effective difficulty level λ ∗ = 1 /β (for sustained en-demic). We predict the transition is continuous near β · ρ ( A ) = 1(see methods), which is in agreement withthe former studies [19]. Discussion.
In conclusion, we have proposed a uni-versal framework to study the tipping points of a complexsystem. First, we use homogeneous networks (sphericalchickens) as a benchmark to gauge the difficulty level ofan environment (determined by the parameter sets). Dueto the symmetry of the nodes in the homogeneous net-work, for any given environment we can just study the1D equation to find the sparsest homogeneous networkthat can survive in the environment. The in-degree ofthis sparsest homogeneous network λ ∗ can be regardedas the effective difficulty level of the environment, andcan usually be studied analytically.Normally a heterogeneous network will have differentresilience under different environmental parameter setcorresponding to an identical effective difficulty level. Weprove that for any given network A , the effective difficultylevel of all the tipping points λ ∗ c are between its max-imum k -core k max ( A ) and its largest eigenvalue ρ ( A ).That is to say this system can always survive in effectivedifficulty level λ ∗ = k max ( A ), and will always collapsein λ ∗ = ρ ( A ). This is particularly meaningful when inreality we need to make the network functions or mal-functions under certain conditions. With the knowledgeof the upper bounds and lower bounds, we will be able tocontrol the parameters/network that can keep the systemin a desirable state. In some special case(step function orcontinuous transition), all the tipping points of a hetero-geneous system have an identical effective difficulty level k max ( A ) or ρ ( A ). In these case, we can even know all theexplicit tipping points.In short, our new findings theoretically determine thebounds of the tipping points of a large class of dy-namical systems, and our results unveil important in-formation for controlling the networks in real world.The methodology we developed is general and may beused in other scientific disciplines, especially when theproblem can be abstracted as a complex object(realchicken)/group(heterogeneous network) evolving underan environment described by multiple parameters. Acknowledgements
This project has received fund-ing from the European Research Council (ERC) underthe European Union’s Horizon 2020 research and inno-vation program under grant agreement No 740704. [1] J. Gao, B. Barzel, and A.-L. Barab´asi, Nature , 307(2016).[2] F. Morone, G. Del Ferraro, and H. A. Makse, Naturephysics , 95 (2019).[3] R. Pastor-Satorras and A. Vespignani, Physical reviewletters , 3200 (2001).[4] A. Barrat, M. Barthelemy, and A. Vespignani, Dynami-cal processes on complex networks (Cambridge universitypress, 2008).[5] A. Lajmanovich and J. A. Yorke, Mathematical Bio-sciences , 221 (1976).[6] M. Eidsaa and E. Almaas, Physical Review E , 062819(2013).[7] Y.-X. Kong, G.-Y. Shi, R.-J. Wu, and Y.-C. Zhang,Physics Reports (2019).[8] K. Shin, T. Eliassi-Rad, and C. Faloutsos, in (IEEE, 2016) pp. 469–478.[9] S. Balaji, M. M. Babu, L. M. Iyer, N. M. Luscombe,and L. Aravind, Journal of molecular biology , 213(2006).[10] S. Gama-Castro, V. Jim´enez-Jacinto, M. Peralta-Gil, A. Santos-Zavaleta, M. I. Pe˜naloza-Spinola,B. Contreras-Moreira, J. Segura-Salazar, L. Muniz-Rascado, I. Martinez-Flores, H. Salgado, et al. , Nucleicacids research , D120 (2008).[11] F. Chung, L. Lu, and V. Vu, Proceedings of the NationalAcademy of Sciences , 6313 (2003).[12] C. Castellano and R. Pastor-Satorras, Physical ReviewX , 041024 (2017).[13] A. Sorensen, Oecologia , 242 (1981).[14] P. P. Staniczenko, J. C. Kopp, and S. Allesina, Naturecommunications , 1 (2013).[15] M. S. Mariani, Z.-M. Ren, J. Bascompte, and C. J. Tes-sone, Physics Reports (2019).[16] J. J. Lever, E. H. van Nes, M. Scheffer, and J. Bas-compte, Ecology letters , 350 (2014).[17] R. P. Rohr, S. Saavedra, and J. Bascompte, Science ,1253497 (2014).[18] S. Saavedra, R. P. Rohr, J. M. Olesen, and J. Bascompte,Ecology and evolution , 997 (2016).[19] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, andA. Vespignani, Reviews of modern physics , 925 (2015). [20] A. Berman and R. J. Plemmons, Nonnegative matricesin the mathematical sciences (SIAM, 1994) p. 27.[21] L. Collatz, Mathematische Zeitschrift , 221 (1942).[22] H. Wielandt, Mathematische Zeitschrift , 642 (1950). METHODSA. Properties
We study the systems that follow the dynamic equa-tion: dx i dt = F ( x i , N (cid:88) j =1 A ij G ( x i , x j )) , (11) x i ≥ i , weighted A ij ≥ j to i .The cooperative systems we study in this paper have thefollowing properties:(i) F ( x = ) = 0. Assume F ( x = ) >
0, then ob-viously there exist a nonzero solution. And abundancecannot be negative, so it is impossible that F ( x = ) < F ( x i > , < ∂F/∂u i ≥
0, here u i = (cid:80) Nj =1 A ij G ( x i , x j ) is the secondvariable of F ;(iv) Higher abundance produce more output: ∂G∂x j ≥ G ( x i ,
0) = 0;(v) x cannot increase to infinite, lim (cid:107) x (cid:107)→ + ∞ d (cid:107) x (cid:107) dt < B. Conditions
Here we study the following three conditions:(i) The system A follows Eq. 11 have nonzero solution;(ii) A homogeneous system with total weighted in-degree equals to k max ( A ) has positive solution. In an-other word, 1D equation F ( x, k max ( A ) G ( x, x )) = 0 haspositive solution x ∗ k ;(iii) A homogeneous system with total weighted in-degree equals to ρ ( A ) has positive solution. In anotherword, 1D equation F ( x, ρ ( A ) G ( x, x )) = 0 has positivesolution x ∗ ρ .We will show that, condition (ii) is the sufficient con-dition of condition (i), condition (iii) is the necessarycondition of condition (i); when G ( x i , x j ) is a step func-tion of x j , condition (ii) is equivalent to condition(i); andwhen F from zero solution to nonzero solution undergoesa continuous transition, condition (iii) is equivalent tocondition (i). C. Lemmas
Here we introduce two important lemmas we will uselater:(i) A ≥ B ≥
0, and B ij ≤ A ij , ∀ i, j , then ρ ( B ) ≤ ρ ( A ). [20](ii) Collatz–Wielandt theorem [21, 22]: A ≥ ∀ x > ,then min ≤ i ≤ n A ij x i x i ≤ ρ ( A ) ≤ max ≤ i ≤ n A ij x i x i (12) D. Reasoning from condition (ii) to condition (i)
By definition of k-core, we can delete the nodes anddecrease the weight of edges in system A to constructa system B that every node has total weighted in-degreeequals to k max ( A ), we denote the nodes in B by 1 , , ..., n .Obviously, x = x = ... = x ∗ k is a solution of system B .Obviously the corresponding nodes 1 , , ..., n in system A must have a solution no less than x ∗ k . E. Reasoning from condition (i) to condition (iii)
The system A has nonzero solution indicates that thereexist some effective interactions G ( x ∗ i , x ∗ j ) >
0, we removeall the nodes that have zero abundance, and all the edgesthat have zero influence G ( x ∗ i , x ∗ j ) = 0. We denote theremaining nodes as 1 , , ..., n , and the ”effective” adja-cency matrix B n ∗ n . Obviously, in system B , x ∗ i will staythe same as it in system A , and x ∗ i > ∀ ≤ i ≤ n . Wehave N (cid:88) j =1 A ij G ( x ∗ i , x ∗ j ) = n (cid:88) j =1 B ij G ( x ∗ i , x ∗ j ) (13) By Lemma (i), (ii) and ρ ( B ) = ρ ( (cid:20) B
00 0 (cid:21) ) , there existscomponent i (1 ≤ i ≤ n ) such that n (cid:88) j =1 B ij G ( x ∗ i , x ∗ j ) /G ( x ∗ i , x ∗ i ) ≤ ρ ( B ) ≤ ρ ( A ) . (14)Hence dx i dt x ∗ i = F ( x ∗ i , N (cid:88) j =1 A ij G ( x ∗ i , x ∗ j )) (15)= F ( x ∗ i , n (cid:88) j =1 B ij G ( x ∗ i , x ∗ j )) (16) ≤ F ( x ∗ i , ρ ( A ) G ( x ∗ i , x ∗ i )) (17) dx ∗ i /dt = 0, then F ( x ∗ i , ρ ( A ) G ( x ∗ i , x ∗ i )) ≥ . (18) F ( x, ρ ( A ) G ( x, x )) = 0 must have solution no less than x ∗ i . F. Reasoning from condition (i) to condition (ii)
When G ( x i , x j ) is a step function: G ( x i , x j ) = (cid:26) x j < αg ( x i ) if x j ≥ α (19)The system has nonzero solutions indicates the existenceof a group of ”symbionts” whose abundance are all largerthan α . We denote their id as 1 , , ..., n and the inter-action matrix within the group is B . The fixed-pointsolution F ( x ∗ i , N (cid:88) j =1 A ij G ( x ∗ i , x ∗ j )) = F ( x ∗ i , n (cid:88) j =1 B ij g ( x ∗ i )) = 0 (20) (cid:80) B mj is the smallest among all (cid:80) B ij , 1 ≤ i ≤ nF ( x ∗ m , n (cid:88) j =1 B mj g ( x ∗ m )) = 0 . (21)By definition, (cid:80) nj =1 B mj ≤ k max ( A ), therefore F ( x ∗ m , k max ( A ) g ( x ∗ m )) ≥ F ( x, k max ( A ) g ( x )) = 0 must have a so-lution x ∗ no less than x ∗ m . On the other hand x ∗ m ≥ α , hence F ( x ∗ , k max ( A ) g ( x ∗ )) is equivalent to F ( x ∗ , k max ( A ) G ( x ∗ , x ∗ )). G. Reasoning from condition (iii) to condition (i)
If the graph is not strongly connected, that is to say wecan divide it to different strongly connected subgraphs,and there is no feedback loop among different subgraphs.Obviously, the system will have a nonzero solution if anyof these strongly connected subgraphs in isolation has anonzero solution.On the other hand, we can always write the adjacencymatrix in a block triangular form, and easily prove thatthe largest eigenvalue of the whole graph equals to thelargest one among all the largest eigenvalues of thesestrongly connected subgraphs.The adjacency matrix A of a strongly connected sub-graph is an irreducible matrix, Perron–Frobenius Theoryguarantees that ρ ( A ) is positive, and the dominant eigen-vector such that ω T A = ρ ( A ) ω T is positive. We will usethis property in the following proof.In what follows we use Lyapunov’s second method toillustrate the local stability at x = . Note A is a irre-ducible matrix in the following proof.
1. Gene regulatory network dx i dt = − x fi + N (cid:88) j =1 (cid:101) A ij x hj x hj , (23)near x = , dx i dt = − x fi + N (cid:88) j =1 (cid:101) A ij x hj . (24)The corresponding 1D equation dxdt = − x f + ρ ( (cid:101) A ) x h . (25)Obviously, when f < h , it is locally stable at x = 0, andits nonzero solution cannot be studied near x = 0; when f > h , it is unstable at x = 0; the two results are trivial.When f = h , the 1D equation is unstable at x = 0 if andonly if ρ ( (cid:101) A ) > A when f = h .Choose Lyapunov function V ( x ) = ω · x , here ω is thedominant eigenvector of A such that ω T A = ρ ( A ) ω T .Obviously V ( x ) > x (cid:54) = 0, and dV ( x ) dt = ω · d x dt = ( ρ ( (cid:101) A ) − N (cid:88) j =1 w j x hj . (26)When ρ ( (cid:101) A ) > dV ( x ) dt > x near but x (cid:54) = .The system is unstable at x = , thus the system followsEq. 23 must have a nonzero solution. Condition (iii) canlead to condition (i) when f ≥ h . Especially, when f = h ,we have λ ∗ c = ρ ( A ).
2. Mutualistic network
Then we study the mutualistic network ruled by thefollowing dynamic equation: dx i dt = − x i d − x i s + γ (cid:80) Nj =1 A ij x j α + (cid:80) Nj =1 A ij x j x i , (27)when d = 0, dx i dt = ( − x i s + γ (cid:80) Nj =1 A ij x j α + (cid:80) Nj =1 A ij x j ) x i , (28)its nonzero solution is equivalent to the solution of thefollowing system dx i dt = − x i s + γ (cid:80) Nj =1 A ij x j α + (cid:80) Nj =1 A ij x j . (29)Near x = , dx i dt = − x i s + γα N (cid:88) j =1 A ij x j , (30)The rest part of proof is similar to the case of gene regula-tory networks. We have the conclusion that when d = 0, λ ∗ c = ρ ( A ).
3. SIS model dx i dt = − x i + (1 − x i ) N (cid:88) j =1 (cid:101) A ij x j , (31)near x = , dx i dt = − x i + N (cid:88) j =1 (cid:101) A ij x j . (32)Again, similarly to the case of gene regulatory networks,we have the conclusion that the SIS model always has λ ∗ c = ρ ( A ).
4. general