Emergent stability in complex network dynamics
Chandrakala Meena, Chittaranjan Hens, Suman Acharyya, Simcha Haber, Stefano Boccaletti, Baruch Barzel
DDynamic stability of complex networks
Chandrakala Meena , , Chittaranjan Hens , Simcha Haber , Stefano Boccaletti − & Baruch Barzel , , ∗
1. Department of Mathematics, Bar-Ilan University, Ramat-Gan, Israel2. Gonda Multidisciplinary Brain Research Center, Bar-Ilan University, Ramat-Gan, Israel3. Physics and Applied Mathematics Unit, Indian Statistical Institute, Kolkata, India4. CNR - Institute of Complex Systems, Via Madonna del Piano 10, Sesto Fiorentino I-50019, Italy5. Unmanned Systems Research Institute, Northwestern Polytechnical University, Xi’an 710072, China6. Moscow Institute of Physics and Technology (National Research University), 9 Institutskiy per.,Dolgoprudny, Moscow Region 141701, Russian Federation * Correspondence : [email protected] In memory of Prof. Robert May
Will a large complex system be stable? This question, first posed by May in 1972,captures a long standing challenge, fueled by a seeming contradiction between the-ory and practice. While empirical reality answers with an astounding yes, the math-ematical analysis, based on linear stability theory, seems to suggest the contrary -hence, the diversity-stability paradox. Here we settle this dichotomy, by consideringthe interplay between topology and dynamics. We show that this interplay leadsto the emergence of non-random patterns in the system’s stability matrix, leadingus to relinquish the prevailing random matrix-based paradigm. Instead, we offer anew matrix ensemble, which captures the dynamic stability of real-world systems.This ensemble helps us analytically identify the relevant control parameters thatpredict a system’s stability, exposing three broad dynamic classes: In the asymp-totically unstable class, diversity, indeed, leads to instability `a la
May’s paradox.However, we also expose an asymptotically stable class, the class in which most realsystems reside, in which diversity not only does not prohibit, but, in fact, enhancesdynamic stability. Finally, in the sensitively stable class diversity plays no role, andhence stability is driven by the system’s microscopic parameters. Together, ourtheory uncovers the naturally emerging rules of complex system stability, helpingus reconcile the paradox that has eluded us for decades.
The study of complex social, biological and technological systems, is often directed towardsdramatic events, such as cascading failures or abrupt state transitions . In reality,however, these represent the exception rather than the rule. In fact, the truly intriguingmathematical puzzle is how, despite enduring constant perturbations and local obstructions,many systems continue to sustain reliably stable dynamics . This challenge is rooted inthe diversity-stability paradox , which predicts that a sufficiently large complex system, willinevitably become unstable. What, then, are the naturally emerging organizing principles, orin May’s original wording - what are nature’s devious strategies , to achieve this ubiquitouslyobserved dynamic stability?Initial clues began to unveil with the mapping of empirical networks, which uncovered universallyrecurring topological characteristics that can potentially impact the system’s dynamics .1 a r X i v : . [ n li n . AO ] A ug or example, degree heterogeneity , community structure and topological symmetries - allnaturally emerging features of real-world networks - indeed, impact the system’s dynamics .However, on their own, these topological features are insufficient to ensure stability , leavingMay’s mathematical challenge unresolved.Here, we show that the desired devious strategies arise not just from the network topology , butmainly from its interplay with the system’s intrinsic, often nonlinear, interaction dynamics .Constrained by the physical mechanisms that drive the system’s interacting components, thesedynamics provide the desired, naturally emerging, design principles that determine the system’sstability. Our analysis shows that, under real-world dynamics, stability can be fully predictedfrom a small set of relevant parameters that capture the combined contribution of both topology and dynamics. Most crucially, we identify a broad class of frequently encountered mechanisms,for which a large complex system can, and, at times, even must be stable, therefore settling thelong overdue diversity-stability debate . The diversity-stability challenge
Consider a complex system of N interacting components (nodes), whose dynamic activities x ( t ) = ( x ( t ) , . . . , x N ( t )) are driven by nonlinear pairwise interactions. The system’s fixed-points capture static states, which, unperturbed, remain independent of time. The stability ofthese fixed-point can be examined through their response to small perturbations δ x , which, inthe linear regime, can be approximated byd δ x d t = J δ x + O ( δ x ) , (1)where J represents the system’s Jacobian matrix. We express J as J = W − CI. (2)The off-diagonal elements W ij , extracted from an arbitrary distribution P ( w ), capture the direct(linear) impact of node j on node i ; often P ( w ) is taken to be a normal distribution N (0 , σ ).Using I to represent the identity matrix, J ’s diagonal terms J ii = − C , capture the typicalrelaxation times of all the nodes, which are dictated by the system’s specific rate parameters.Together, P ( w ) and C define the random matrix ensemble from which J is extracted. We denotethis ensemble by E ( P ( w ) , C ).The system E ( P ( w ) , C ) is dynamically stable if J ’s principle eigenvalue satisfies Re( λ ) <
0. Thechallenge is that, according to random matrix theory, in the limit of large N we have Re( λ ) ∼ √ N (cid:18) − C √ N (cid:19) , (3)which is positive, and hence unstable, as long as C < √ N . Consequently, a sufficiently largesystem ( N → ∞ ) becomes inevitably unstable, i.e . May’s paradox. The role of the network topology . A natural attempt to reconcile this paradox is by consid-ering structural patterns in the interaction topology A ij , a binary network ( A ij = 0 , . As J is designed to capture the direct response between i and j , its terms J ij vanish in case there is no direct i, j link. We, therefore,write the Jacobian in (2) as J = A ⊗ W − CI , where the Hadamard product ⊗ represents matrix2ultiplication element by element. Hence J ij = 0 in case there is no link between i and j ,and is assigned a weight from P ( w ) if i and j interact. This provides the Jacobian ensemble E ( A ij , P ( w ) , C ), in which one first selects the topology A ij , then assign off-diagonal and diagonalweights from P ( w ) and C (Fig. 1a,b).In the E ( A ij , P ( w ) , C ) ensemble, the principal eigenvalue λ depends on the average nearestneighbor degree κ nn , whose magnitude depends on the system’s degree heterogeneity . Forinstance, in a random network with degree distribution P ( k ) we have κ nn = (cid:104) k (cid:105) / (cid:104) k (cid:105) , in whichthe second moment (cid:104) k (cid:105) increase with P ( k )’s variance, and hence with A ij ’s heterogeneity. Mostgenerally, for a sufficiently broad, i.e . fat-tailed, P ( k ), we have κ nn ∼ N β , (4)an asymptotic divergence with system size. As an example, consider a scale-free network, with P ( k ) ∼ k − γ (2 < γ < κ nn follows (3.2) with β = 3 − γ . Under Eq. (3.2) theprincipal eigenvalue in the E ( A ij , P ( w ) , C ) ensemble follows (Fig. 1c)Re( λ ) ∼ N β (cid:32) − CN β (cid:33) , (5)in which the square-root of (3) is replaced by the exponent β/
2. The crucial point is that, onceagain, we observe an asymptotic behavior that prohibits stability in the limit N → ∞ . These mathematical observations, now nearing their th anniversary, seem to suggest that com-plex networked systems must be limited in size in order to maintain stability. Reality, however,consistently confronts us with the contrary . To reconcile this incongruity, we next show thatunder real nonlinear dynamics the structure of J is fundamentally different from the currentlyconsidered Jacobian ensembles. This allows us to fully predict a complex system’s dynamicstability, and, down the line, disentangle May’s theoretical gridlock. The role of dynamics
The common thread that binds the present treatments of the diversity stability paradox - indeedthe paradigm - is that one can separate the role of the topology from that of the dynamics. Forexample, to construct J in the E ( A ij , P ( w ) , C ) ensemble, one first selects the topology A ij , andthen independently extracts the dynamic response terms, W ij , from P ( w ), and the relaxationrate C . Such construction overlooks the complex interplay between A ij and P ( w ) , C , rootedin the nonlinear mechanisms driving the interactions between the nodes. Omitting this inter-play is tantamount to ignoring the role of the system’s actual nonlinear dynamics. Indeed,in E ( A ij , P ( w ) , C ), if two networks share a similar topology, they will have an equally similar J (up to statistical variations), regardless of their interaction dynamics - social, biological ortechnological. This, of course, stands in sharp contrast with the frequently observed fact that dif-ferent nonlinear mechanisms, despite being run on similar networks, may express fundamentallydistinct dynamic behaviors .Hence, to truly predict dynamic stability, we derive the mathematical link between A ij and P ( w ) , C , seeking the relevant ensemble from which to extract real-world Jacobian matrices. Dynamic mechanisms (Fig. 1d). To construct realistic J matrices we must consider eachsystem’s intrinsic dynamic mechanisms. These mechanisms are inherent to the nature of the3ystem’s interacting components: in social systems, for example, individuals interact throughinfection and recovery , in biological networks, proteins, genes and metabolites are linkedthrough biochemical processes , and in ecological systems, species undergo competitiveor symbiotic interactions . All these processes capture built-in dynamics, representing the physical mechanisms that drive all nodes/links. Most generally, these dynamic mechanisms canbe described by d x i d t = M (cid:0) x i ( t ) (cid:1) + N (cid:88) j =1 A ij M (cid:0) x i ( t ) (cid:1) M (cid:0) x j ( t ) (cid:1) , (6)where M ( x ) captures the self-dynamics of all nodes, and M ( x ) , M ( x ) describe their pairwiseinteractions. With the appropriate selection of these nonlinear functions Eq. (1.1) captures abroad range of frequently used models in social , biological and technological systems (Fig. 1g). The dynamic Jacobian ensemble (Fig. 1e). Extracting J from (1.1), we show in Supplemen-tary Section 1, leads to a currently unexplored matrix ensemble, in which the weights J ij arestrongly intertwined with the topology A ij , via J ii ∼ − Cκ η nn k µi (7) J ij ∼ A ij k νi k ρj . (8)In (7) and (1.4), κ nn and k i , k j are extracted from the topology A ij and its degree distribution P ( k ). The four exponents, Ω = ( η, µ, ν, ρ ) are, on the other hand, fully determined by the dynamics , M ( x ) , M ( x ) , M ( x ), as shown in Box I . The coefficient C , similar to the Jacobianin (2), represents the system specific rate parameters and characteristic time-scales, a parameterthat can potentially be affected by external perturbation or changing environmental conditions .The result is a non-random Jacobian structure, whose entries, as opposed to the existing en-sembles, cannot be selected independently from an arbitrary P ( w ), or set to a constant value C . Instead (7) and (1.4) define a new ensemble E ( A ij , Ω), which matches, for each underlyingnetwork A ij , a predictable set of weights J ii and J ij based on the dynamic exponents Ω. There-fore, in contrast to the random matrix based treatment of E ( P ( w ) , C ) or E ( A ij , P ( w ) , C ), thathas led the discussion since May’s original analysis , the newly introduced E ( A ij , Ω) capturesthe effect of the system-specific nonlinear dynamics. Indeed, in (7) and (1.4), identical networksmay yield highly distinctive Jacobian matrices, depending on whether the interactions are, e.g. ,social, biological or ecological, as each dynamics is characterized by its unique set of exponentsΩ (Fig. 1g).To examine this prediction we implemented the dynamics of Fig. 1g on a set of model and rel-evant empirical networks. This includes four dynamic models:
Social . The susceptible-infected-susceptible (SIS) model for disease spreading;
Regulatory . The Michaelis-Menten model forgene regulation; Ecological . Mutualistic interactions in ecology; Biochemical . Protein-proteininteractions in sub-cellular networks. Applying each model to five different model andreal networks, we arrive at a total of 20 combinations of networks/dynamics as a testing groundfor our predicted J -ensemble (for a detailed description of all dynamic models and networks seeSupplementary Sections 2 and 4.4).Perturbing the system around its numerically obtained fixed-points, we constructed the Jacobian4atrix J for each of our 20 systems. In Fig. 2 we find that, indeed, the diagonal and off-diagonal terms of J follow the predicted scaling of (7) and (1.4). For example, in our Socialmodel we predict µ = 1, while our Regulatory dynamics are predicted to have µ = 0, bothscaling relationships clearly evident in Fig. 2a and c. This means that setting all diagonal termsuniformly to − C , as in May’s original formulation, misses the actual patterns that arise fromthe nonlinear Social/Regulatory dynamics. Similarly, the off-diagonal terms are proportionalto k − i k j in Social (Fig. 2b) and k i k − j in Regulatory (Fig. 2d), again overruling the classicconstruction in which the J ij entries are extracted from an arbitrary P ( w ). These scalingrelationships, encapsulated within our analytically predicted Ω, are independent of A ij , as indeedFig. 2 confirms, showing recurring scaling patterns across our model (SF1, SF2, SF3) and relevantempirical (Epoch, UCIonline, PPI1 etc.) networks. Hence, Ω captures the intrinsic, and mostcrucially, hitherto overlooked, contribution of the nonlinear dynamics to the structure of J .To further observe J in an empirical setting, we examine in Fig. 3 the stability of our empiricalsocial, biological and ecological networks. For each network we construct J both from the randomensemble E ( A ij , P ( w ) , C ) (left) and from the dynamic ensemble E ( A ij , Ω) (right), where we takeΩ in accordance with the relevant dynamic mechanisms (Fig. 1g). We find that the two Jacobianmatrices, random vs. dynamical, despite belonging to the same network, are visibly different,as the dynamic Jacobian incorporates the emergent non-random patterns predicted in (7) and(1.4). These differences directly impact the network’s dynamic stability. Indeed, in all cases wefind that the random J has a positive λ (red), congruent with May’s analysis, while the dynamic J , extracted from E ( A ij , Ω) has λ <
Hence, these empirical networks, when coupledwith relevant real-world dynamics, are, in fact stable .Together, our analysis demonstrates that: (i) actual J matrices, extracted from real nonlineardynamics are fundamentally distinct from random matrices; (ii) contrary to the random J ensembles, real Jacobians feature scaling patterns in which topology ( P ( k )) and dynamics (Ω)are deeply intertwined; (iii) these patterns can be analytically traced to the system’s intrinsicnonlinear mechanisms ( M ( x ) , M ( x ) , M ( x )) through Eqs. (7) and (1.4). Most crucially, whilethe random ensembles E ( P ( w ) , C ) and E ( A ij , P ( w ) , C ) lead to the instability paradox of Eqs.(3) and (5), as we demonstrate in Fig. 3, and show analytically below E ( A ij , Ω) has a muchricher stability profile, in which a broad family of empirically relevant dynamics are predictablystable.
Topology, dynamics and the emergence of stability
To address May’s paradox systematically we analyze the E ( A ij , Ω) ensemble under a fat-tailed P ( k ), as captured by Eq. (3.2). This setting allows us to observe the combined effect of dynamicnonlinearity (Ω) with one of the most ubiquitous features of real world networks, i.e . degree-heterogeneity . In Supplementary Section 3 we show that extracting J from E ( A ij , Ω), with adegree-heterogeneous A ij , the principal eigenvalue asymptotically followsRe( λ ) ∼ N Q (cid:18) − CN S (cid:19) , (9)where Q and S depend both on network structure, through β in (3.2), and on the nonlineardynamics via Ω. Specifically, we show that S = β (1 + ν + ρ − µ − η ) , (10)5eplacing the fixed exponents, 1 / E ( P ( w ) , C ), and β/ E ( A ij , P ( w ) , C ),with a variable exponent that depends on the dynamics-driven Ω = ( η, µ, ν, ρ ).Equations (9) and (3.31) represent our key result. They show that, when extracted from realnonlinear dynamics of the form (1.1), the asymptotic behavior of λ is different from the onepredicted in May’s original formulation, i.e . Eqs. (3) and (5). Most importantly, (9) and (3.31)have crucial implications regarding the system’s dynamic stability, distinguishing between threepotential stability classes (Fig. 1f): • Asymptotic instability:
S > S in (9) is positive, we have, for suffi-ciently large N , Re( λ ) ∼ N Q >
0. Such systems commit to the classic diversity-stabilityprediction, and, when large enough, become intrinsically unstable, regardless of the modelparameters ( C ). Specifically, setting η = µ = ν = ρ = 0 we revert to the random matrixensemble ( E ( A ij , P ( w ) , C )), for which S = β >
0, thus recovering May’s predicted insta-bility. Hence, the diversity stability paradox is, indeed, a valid prediction, yet it representsa specific point in a broader class of potential dynamic behaviors, driven by the parameter S in (3.31). • Asymptotic stability:
S <
S <
0. Under these conditions, independently of C , the r.h.s. of (9) is dominatedby the negative term, guaranteeing that Re( λ ) <
0. Therefore, in this class, stability isnot simply enabled , but rather, under N → ∞ , it becomes asymptotically inevitable, evenif C (cid:28)
1. Recalling May’s original question: can a large complex system be stable ? - Eq.(9) with
S < can , but, in fact, must be stable. • Sensitive stability: S = 0 (green). Under S = 0 the system lacks an asymptotic behavior,and therefore, its stability depends on C in (9). If C > i.e . sufficiently strong negativefeedback, the system is stable, otherwise it becomes unstable. As opposed to Ω, which isintrinsic to the dynamic model , depending on the functional form of M ( x ) , M ( x ) , M ( x )in (1.1), C is determined by the specific model parameters . Hence, in this class stability isnot an intrinsic characteristic of the system, but rather it is sensitive to Eq. (1.1)’s tunableparameters.This classification settles the debate on complex system stability on several levels. In the con-text of the diversity-stability paradox, it shows that large complex systems can , and in somecases even must be stable - hence reconciling between theory and empirical observation. Morebroadly, the stability classifier S , identifies the relevant topological ( β ) and dynamic (Ω) con-trol parameters that help analytically predict the stability class of any system within the form(1.1). We can therefore use S to predict a priori whether a specific combination of topologyand dynamics will exhibit stable functionality or not.To examine our stability classifier S we used our model and real networks to extract 2 , E ( A ij , Ω) ensemble, with different sets of η, µ, ν and ρ . In Fig. 4awe show the principal eigenvalue λ vs. S for the entire 2 ,
077 Jacobian sample. As predicted,we find that the parameter S sharply splits the sample into three classes. The asymptoticallyunstable class (red, top-right) has S > λ >
0, a guaranteed instability.The asymptotically stable class (blue, bottom-left) is observed for
S <
0, and has, in all cases λ < i.e . stable dynamics, as predicted. Finally, for S = 0 we observe sensitive stability, with λ having no asymptotic positive/negative assignment (green). A small fraction ( ∼ J matrices were inaccurately classified by S (grey), a consequence of the approximatenature of our derivations (Supplementary Section 3).6n Fig. 4a we also present the stability matrices associated with our empirical social, biologicaland ecological networks (symbols) examined earlier in Figs. 2 and 3. For each of these networkswe constructed J via the E ( A ij , P ( w ) , C ) ensemble, assigning random weights along the links,as well as via the E ( A ij , Ω) ensemble, with weights determined by (7) and (1.4), and Ω takenaccording to the relevant dynamic model (Fig. 1g). While all these networks are classifiedunstable under the random ensemble (red symbols), as per May’s diversity-stability principle,once fit with J from our dynamic ensemble they all transition into the asymptotically stableclass (blue symbols).Figures 2 - 4, together demonstrate our complete theoretical path: First, Figs. 2 and 3 show thatreal J matrices exhibit the internal scaling patterns predicted in (7) and (1.4) - representing asteep departure from the broadly applied random matrix paradigm. Next, Fig. 4 shows thatthe resulting Jacobian ensemble has a rich space of potential dynamic behaviors, significantlymore diverse than the currently considered ensembles. Most importantly - this space contains abroad class of asymptotically stable (blue) or sensitively stable (green) systems, that do not fallwithin May’s paradox. Taken together, our analysis demonstrates that nature must not rely on devious strategies inorder to ensure dynamic stability. The observed stability of large complex systems emergesquite naturally thanks to the built-in nonlinear interaction mechanisms between the system’scomponents.
The ingredients of dynamic stability
The parameter S in (3.31) reduces a network’s dynamic stability into five relevant exponents.The first four Ω = ( η, µ, ν, ρ ) are determined by the nonlinear dynamics, Social, Regulatory,Ecological etc., and are therefore intrinsic to the system’s inherent interaction mechanisms.These exponents are independent of the topology A ij or of the microscopic model parameters,all of which are encapsulated within C . Therefore they are hardwired into the system’s dynamicbehavior. To understand this, consider, for example, our Social model, for which our formalismpredicts Ω = (0 , , − , i.e .infection vs. recovery, expressed through the functional form of M ( x ) , M ( x ) , M ( x ) of (1.1).It is, therefore, insensitive to the specific parameters of the model, e.g. , if the disease has a highor low infection rate, or if it is transmitted via physical contact or aerosols. Such distinctions,expressed through the model parameters , may impact the coefficient C , but do not affect Ω, andtherefore have no bearing on the stability classifier S . Hence, S is characteristic of, e.g. , the SIS model , not of its specific parameters .The remaining exponent in (3.31), β , is independent of the dynamics, determined solely by A ij , specifically by its degree distribution P ( k ), through (3.2). This parameter quantifies thefat-tailed nature of P ( k ), being β = 0 for a bounded distribution, and approaching unity underextreme levels of degree-heterogeneity. Therefore, together, S captures the roles of both topology and dynamics, whose interplay determines the system’s stability class - stable, unstable orsensitive. The role of degree-heterogeneity . The prevalence of fat-tailed P ( k ) is among the definingfeatures of real-world complex systems, from biological to social and technological networks, with far-reaching implications on their observed dynamic behavior . Our analysisindicates that this network characteristic may also play a crucial role in the context of dynamicstability. To understand this, consider a non fat-tailed P ( k ), such as an Erd˝os-R´enyi, network,7here the degrees follow a Poisson distribution, having β = 0. Under these conditions we have S = 0 in (3.31), the system has no defined asymptotic behavior, and hence it is sensitivelystable - i.e . its stability depends on model parameters. Therefore, the existence of asymptoticstability/instability is a direct consequence of degree-heterogeneity, as, indeed, these forms ofdynamic stability rely on β > P ( k ) on complex system func-tionality. Consider, for example, the factors that drive a system towards the loss of stability.Most often such events result from external stress or changes in environmental conditions . Suchforces impact a system’s functionality by perturbing its dynamic parameters, namely they affectthe coefficient C . Seldom, however, do these environmental perturbations affect the system’sbuilt-in interaction mechanisms. Indeed, while the dynamic mechanisms are fixed, ingrained inthe physics of the interacting components, the specific model parameters often depend on exter-nal conditions. The crucial point is, that under S <
0, a state only possible if P ( k ) is fat-tailed( i.e . β (cid:54) = 0), stability becomes asymptotically independent on C , driven solely by Ω, which isinsensitive to specific model parameters. Hence, for a sufficiently large system ( N → ∞ ), if P ( k )is fat-tailed, stability becomes asymptotically robust against any external perturbation. Thissuggests that degree-heterogeneity serves as a dynamically stabilizing topological characteristic,in the face of a persistently fluctuating environment .To illustrate this, consider again the SIS model (Social), for which Ω = (0 , , − , S = − β . Implemented on a scale-free network ( β > S <
0, andis thus asymptotically stable. However, the same dynamics on a bounded P ( k ) ( β = 0) becomessensitively stable, i.e . S = 0, and therefore its stability depends on the specific value of C , whichcan be tuned by changing the infection/recovery rates. Consequently, while generally Social’spandemic state is only stable if C in (9) is greater than unity, in a scale-free environment it is always stable, even when C →
0. This observation, predicted here directly from our formalism,recovers an already well-established result: the fact that in the SIS model, the epidemic thresholdvanishes under a scale-free P ( k ) . The crucial point is, however, that while the original result,particular to the SIS model, was obtained using a dedicated model-specific analysis, our controlparameter S allows us to systematically analyze the stability profile of all combinations oftopology/dynamics within (1.1), using a single universal classification.To observe this role of P ( k ), beyond the specific example of Social, we consider again E ( A ij , Ω)’sprincipal eigenvalue λ in (9). Its structure portrays stability as a balance between the positive, i.e . destabilizing, effect mediated by the network interactions, vs. the negative, stabilizing, feed-back, driven by the parameter C in J ’s diagonal (7); Fig. 4b. It is therefore, natural to enhancestability by increasing C , which, in effect, translates to strengthening each node’s intrinsic neg-ative feedback. Equation (9) predicts that a system becomes stable if C exceeds a critical value C ∼ N S , (11)beyond which λ becomes negative. For asymptotically stable systems ( S <
0) we have C → C . In contrast, for asymptotically unstablesystems ( S >
0) we have C → ∞ , hence such systems are impossible to stabilize even underextremely large C , i.e . extreme levels of negative feedback. The point is, that if P ( k ) is boundedwe have, always, S = 0 (sensitive stability), hence C is finite, and the system’s stability can becontrolled by tuning the parameter C , e.g ., changing the environmental conditions.8n Fig. 4c-k we extract again a set of three J matrices from E ( A ij , Ω), representing systems fromour three stability classes: J AS , asymptotically stable with Ω = (2 , , , − J SS , sensitivelystable with Ω = (0 , − , − , J AU , asymptotically unstable with Ω = (1 , − , − , C vs. N , capturing the level of negative feedback required to ensurethe system’s stability. Under an Erd˝os-R´enyi network, with bounded P ( k ) ( β = 0) we do notobserve a defined asymptotic behavior. The critical C does not scale with N , indicating thatsufficient perturbation to the model parameters can affect the system’s stability (Fig. 4i-k, greencircles).The picture, however, fundamentally changes when we construct the same J matrices from ascale-free A ij with β = 0 . S = − . , . J AS , J SS and J AU , respectively. As predicted in (11), under these condition, we observe a clear asymptoticbehavior, in which C scales with N . For J AS we have C ∼ N − . , while under J AU we observe C ∼ N . (solid lines), capturing the intrinsic stability/instability of these two systems. Finally,in J SS , having S = 0, the system does not feature any asymptotic behavior, and therefore canbe stabilized under finite C , independently of system size N (Fig. 4g). Hence, we observe a qualitative difference between bounded vs. scale-free P ( k ) , in which degreeheterogeneity can potentially afford the network a guaranteed stability, that is asymptoticallyindependent of microscopic parameters . Stabilizing vs. destabilizing dynamics . To gain qualitative insight on the role of the non-linear dynamics, beyond the formal prediction of Eq. (3.31), let us focus on a specific test caseof Ecological interactions, in which Eq. (1.1) takes the form (Fig. 1g)d x i d t = Bx i ( t ) (cid:18) − x i ( t ) K (cid:19) + N (cid:88) j =1 A ij x i ( t ) F (cid:0) x j ( t ) (cid:1) . (12)The self-dynamics follows logistic growth and the interaction mechanism describes mutualisticdynamics . The response function F ( x j ) is designed to capture the positive contribution ofnode j ’s activity to that of node i . Often this contribution is taken to follow F ( x j ) = x j ,assuming a linear contribution of a species’ abundance to its symbiotic partners. Such dynamicslead to Ω = (0 , , , S = β >
0, predicting an asymptoticallyunstable dynamics. The source of this instability is rooted in the strong positive feedbackbetween interacting nodes, in which an increase in x j is linearly transferred to x i ’s growth rate.An alternative formulation, which we used in our simulations in Fig. 2e,f, takes F ( x j ) = x j / (1 + x j ). This represents a saturating function, in which j ’s impact on i is bounded, as, indeed, in thelimit x j → ∞ we have F ( x j ) →
1. Under this response function, the exponent ρ changes from ρ = 0, observed in the linear F ( x j ), to ρ = −
2, under the saturating alternative. Consequently,we now have Ω = (0 , , , − S = − β , and therefore, the same network which wasunstable under F ( x j ) = x j is now asymptotically stable thanks to F ( x j )’s saturation.From an intuitive perspective, as F ( x j ) becomes bounded, the interaction is weakened . Thisis expressed formally in (1.4) through the suppressed off-diagonal terms J ij , here - due to thenegative exponent ρ = −
2. This reduction in J ij , in turn, leads to a diminished λ . In thatsense, the stability classifier S in (3.31) captures precisely this balance: if ν, ρ are large - theoff-diagonal terms J ij , which represent positive feedback, dominate the system’s dynamics, S becomes positive and the system is asymptotically unstable. Conversely, if µ, η are increased,the diagonal terms in (7) carry sufficient weight to asymptotically stabilize the system. The case9 = 0 captures a balance between J ii and J ij , rendering the system sensitive to its microscopicparameters. Hence, despite their formal mathematical nature ( Box I ), the four exponents η, µ, ν and ρ provide qualitative insight into the dynamic ingredients that govern the network’s stability(Fig. 5). Levels of dynamic stability . While stability is, in its essence, a binary classification - a systemis either stable or not, we can use S to quantify the level of the system’s stability. Consider againEq. (11), capturing the critical value of C in (7) required to stabilize the system’s dynamics. Anasymptotically stable system has C → i.e . stability is ensured even under arbitrarily small C .In contrast, an asymptotically unstable system has C → ∞ , hence cannot be stabilized underfinite parameters. These effects are pronounced more strongly as S becomes more negative(stable) or more positive (unstable). As an example, let us compare Regulatory dynamics with h = 1 , a = 2, under which we have S = − β/
2, vs. Ecological dynamics, where S = − β . Bothare asymptotically stable, yet Eq. (11) indicates that Ecological is more stable: to destabilizeEcological one needs C Eco0 ∼ N − β , while for Regulatory it is sufficient to have C Reg0 ∼ N − β/ ,hence C Eco0 (cid:28) C Reg0 . The meaning is that to destabilize Ecological we much reach a moreextreme reduction in C than the one required for Regulatory. In Fig. 5b we show the linesof equi-stability, comprising all dynamics with the same S value, and hence a similar level ofstability/instability. These equi-stable classes may group together highly distinct dynamics, thatdespite their different Ω may have identical S . For example, Social and Ecological, two dynamicsfrom fundamentally different domains, that exhibit a similar level of stability, belonging to theclass S = − β (Fig. 1g). Mixed interactions . Our analysis up to this point focused on cooperative mechanisms, inwhich M ( x ) and M ( x ) are positive, as observed, e.g ., in mutualistic ecological dynamics .To complement this, in Supplementary Section 5 we numerically analyze the case of mixeddynamics, in which a fraction f of positive interactions coexists alongside a 1 − f fraction ofnegative ones. We find that for a broad range of f values, our predicted asymptotic stabilitycontinues to be observed. This is attributed, once again, to a naturally emerging organizingprinciple, relevant in many real-world systems: the fact that the activities x i ( t ) in (1.1) are,by definition, non-negative. Therefore, as the system evolves in time, at any instance where anode reaches x i ( t ) = 0 it becomes extinct, and permanently removed from the network. Thisratchet-motion prunes out nodes with many negative interacting partners, as their activities arenaturally driven towards zero. As a result, the system reaches a fixed-point in which the majorityof remaining links is biased towards positive interactions. Under these conditions, the emergentJacobian matrix J Mix , while still featuring some level of negative interactions, is predominantlypositive. Therefore, despite the fact that, strictly speaking, J Mix / ∈ E ( A ij , Ω), it continues to approximately follow the patterns of (7) and (1.4), and most importantly in the present context- continues to allow a broad class of asymptotically stable systems.
Discussion and outlook
The linear stability matrix J carries crucial information on the dynamic behavior of a complexsystem. Here, we exposed distinct patterns in the structure of J that (i) arise from the natureof the system’s interaction dynamics; (ii) affect its principal eigenvalue, and thus, its stabilityprofile. These patterns are expressed through the four dynamic exponents Ω = ( η, µ, ν, ρ ) in(7) and (1.4), linked analytically to the system’s dynamics, M ( x ) , M ( x ) , M ( x ), through theleading powers of the expansions in Eq. (1.65). This dependence on the powers (Ψ( n ) , Γ( n ),etc.), rather than the coefficients ( G n , C n , etc.) indicates that Ω is hardwired into the interaction10ynamics. Indeed, the powers in (1.65) are determined by the dynamic model, e.g ., Social orRegulatory, but not the specific model parameters, which only impact the coefficients. Therefore,our predicted Jacobian ensemble in (7) and (1.4), as well as its associated stability classifier S in (3.31), both capture highly robust and distinctive characteristics of the system’s dynamics,that cannot be perturbed or otherwise affected by shifting environmental conditions.Graph spectral analysis represents a central mathematical tool to translate network structureinto dynamic predictions . A network’s spectrum, i.e . its set of eigenvalues and eigenvectors,captures information on its dynamic time-scales, potential states, and - in the present context- its dynamic stability. In that sense, the long-standing diversity-stability paradox exposed asevere and troubling clash between spectral graph theory and empirical observation. Here wereconcile this contradiction, by showing that, in effect, we have, all this time, been looking at thewrong spectrum . Indeed, in its current application, spectral analysis is applied to the network topology , namely we seek the graph’s eigenvalues. Focusing on the graph topology, however,loses information on the nonlinear dynamics that occur on that graph. Therefore, here, we offerto apply spectral analysis, not to the topology A ij , but rather to the derived Jacobian in (7)and (1.4), which, thanks to Ω preserves the information of both topology and dynamics.In a broader perspective, the interplay between topology and dynamics is one of the major theo-retical obstacles along the path to predict, understand and control complex system behavior .The problem is that we lack sufficient theoretical tools to treat the challenging combination ofnonlinear dynamics with complex, random and often highly heterogeneous network structures.Consequently, advances are often system specific, requiring dedicated tools, and hence, lack-ing the breadth of a general network dynamics theory. The Jacobian ensemble presented here,transforms the nonlinear Eq. (1.1) into an effective linear system, whose weights depend on bothtopology and dynamics. It can, therefore, allow us to utilize the powerful tools of spectral graphtheory, even in a nonlinear dynamic setting. This may offer a basis for systematic analysis ofnonlinear network dynamics. Data availability . Upon acceptance, all codes/data to reproduce our analysis will be madeavailable online. 11
OX I. THE ENSEMBLE E ( A ij , Ω) - HOW TO CALCULATE ΩTo construct J as appears in Eqs. (7) and (1.4) we must extract the exponents Ω = ( η, µ, ν, ρ )from Eq. (1.1). First we use the dynamic functions M ( x ) , M ( x ) and M ( x ) to constructthree new functions R ( x ) = − M ( x ) M ( x ) , Y ( x ) = M ( x ) R (cid:48) ( x ) , Z ( x ) = R ( x ) M ( x ) , (13)where R (cid:48) ( x ) represents a derivative d R/ d x around the fixed point x . From these we extractfour additional functions, which we express through a Hahn power-series expansion as M (cid:0) Z − ( x ) (cid:1) = ∞ (cid:88) n =0 G n x Ψ( n ) , Y (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 C n x Γ( n ) ,M (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 K n x Π( n ) , M (cid:48) (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 L n x Θ( n ) , (14)where R − ( x ) and Z − ( x ) represent the inverse functions of R ( x ) and Z ( x ). The Hahnexpansion is a generalization of the Taylor expansion to include negative and real powers.Hence Ψ( n ) , Γ( n ) , Π( n ) and Θ( n ) represent series of real powers in ascending order, the leading( i.e . smallest) powers assigned the index n = 0. These leading powers directly provide Ω via µ = 2 − Γ(0) , ν = − Π(0) , ρ = − Θ(0) , η = − Ψ(0)( µ − ν − ρ ) . (15)Hence, to construct J ij ∈ E ( A ij , Ω) we first generate the network A ij , then extract the degrees k i of all nodes and the nearest neighbor degree κ nn . The resulting J ij satisfies J ii ∼ − Cκ η nn k µi (16) J ij ∼ A ij k νi k ρj , (17)where the constant C > C , depends only on the leading powers of (1.65), and not on the coefficients . There-fore it is unaffected by Eq. (1.1)’s microscopic parameters, capturing an intrinsic characteristicof each dynamic model. The detailed derivation of Ω appears in Supplementary Section 1,followed by a step by step application on all our dynamic models (Fig. 1g) in SupplementarySection 2. 12 𝒊𝒋 , 𝜷 𝑷(𝒘)
The random ensemble
𝔼(𝑨 𝒊𝒋 , 𝑷 𝒘 , 𝑪) 𝑨 𝒊𝒋 , 𝜷 𝛀 = (𝜼, 𝝁, 𝝂, 𝝆) 𝐝𝒙 𝒊 𝐝𝒕 = 𝑴 𝒙 𝒊 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝑴 𝒙 𝒊 𝑴 (𝒙 𝒋 ) The dynamic ensemble
𝔼(𝑨 𝒊𝒋 , 𝛀) 𝝀 ∼ 𝑵 𝑸 𝑺 𝝀 ∼ 𝑵 𝜷𝟐 𝜷𝟐 Social 𝐝𝒙 𝒊 𝐝𝒕 = −𝒇𝒙 𝒊 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒈 𝟏 − 𝒙 𝒊 𝒙 𝒋 𝛀 = (𝟎,𝟏,−𝟏,𝟎) 𝑺 = −𝜷
Regulatory 𝐝𝒙 𝒊 𝐝𝒕 = −𝑩𝒙 𝒊𝒂 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒋𝒉 𝒋𝒉 𝛀 = 𝟎,𝒂 − 𝟏𝒂 ,𝟎,−𝒉 + 𝟏𝒂 𝑺 = −𝒉𝒂 𝜷
Ecological 𝐝𝒙 𝒊 𝐝𝒕 = 𝑩𝒙 𝒊 𝒊 𝑲 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝑭(𝒙 𝒋 ) 𝛀 = (𝟎,𝟏,𝟏,−𝟐) 𝑺 = −𝜷 Biochemical * 𝐝𝒙 𝒊 𝐝𝒕 = 𝑭 − 𝑾𝒙 𝒊 − ෩𝑩 𝒋=𝟏 𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝒙 𝒋 𝛀 = −𝟏,𝟏,−𝟏,𝟎 𝑺 = −𝜷
𝑺 > 𝟎 →
Unstable
𝑺 = 𝟎 →
Sensitive
𝑺 < 𝟎 →
StableUnstable (a) (b)(d) (e) (c) (f) (g) 𝑱 𝒊𝒋 𝑱 𝒊𝒋 Figure 1:
Solving the diversity-stability paradox . (a) A complex system is captured by aninteraction network A ij , whose diversity can be quantified by β in (3.2). In May’s random matrixframework, to assess dynamic stability, we construct the Jacobian J with weights extractedindependently from the distribution P ( w ), and diagonal term set to − C (constant). (b) Theresulting stability matrix J ij , belonging to the ensemble E ( A ij , P ( w ) , C ). (c) The principaleigenvalue, under this construction, diverges with the number of interacting components N ,predicting that a large system cannot be stable. This non-realistic prediction captures theunresolved diversity-stability paradox. (d) Our alternative J -ensemble, E ( A ij , Ω), is designed tocapture the roles of both the network ( A ij , β , grey) and the nonlinear dynamic mechanisms ofinteraction, e.g. , Biochemical, Social or Regulatory (orange). These mechanisms are describedby M ( x ) , M ( x ) , M ( x ) in Eq. (1.1), from which we extract the four exponents of Ω. Here,the diagonal and off-diagonal weights of J are not selected independently from P ( w ) , C , butrather predicted by Eqs. (7) and (1.4). (e) As a result, the same A ij may lead to different J matrices, as different dynamics are characterized by distinct Ω. (f) Consequently, λ follows adifferent form than in (c), with a variable stability classifier S , whose sign predicts the system’sstability: Unstable (red) in which diversity, indeed, leads to instability; Sensitive (green), wherediversity plays no role, and Stable (blue), in which diversity enhances stability. Hence, underreal nonlinear dynamics a large system can ( S = 0), and in some cases even must ( S < S , as detailed in Supplementary Section 2. ∗ For the calculation of S in Biochemicalsee Supplementary Section 3.2. 13 oc i a l R e g u l a t o r y E co l o g i c a l B i oc h e m i c a l (a) (b)(c) (d)(e) (f)(g) (h) Figure 2:
Emerging patterns in the dynamic ensemble E ( A ij , Ω). We implemented Social,Regulatory, Ecological and Biochemical dynamics (Fig. 1g) on a set of relevant model (darkblue) and empirical (light blue) networks. Perturbing all nodes we numerically constructedeach system’s Jacobian, J , capturing real Jacobian matrices, as obtained from actual nonlineardynamic models. (a) The diagonal terms J ii vs. degree k i as obtained from Social dynamics(symbols). We observe the scaling relationship of Eq. (7) with µ = 1 (solid line) confirmingour predicted scaling for this dynamics. (b) The off diagonal terms J ij vs. prediction (1.4) with ν = − , ρ = 0 (symbols). Once again, we observe a perfect agreement with our theoreticalprediction (solid line). We also include results obtained from two relevant empirical networks,Epoch and UCIonline , capturing online social dynamics. (c) - (d) Similar results are observedunder Regulatory dynamics ( µ = ν = 0 , ρ = −
2) on both model and real networks (PPI1 andPPI2 ); (e) - (f) Ecological dynamics ( µ = ν = 1 , ρ = −
2, empirical networks: ECO1 andECO2 ); (g) - (h) Biochemical dynamics ( µ = 1 , ν = − , ρ = 0, empirical networks: PPI1 andPPI2). As predicted, we find that real Jacobian matrices exhibit the non-random patterns of Eqs.(7) and (1.4), and therefore cannot be modeled via the random ensemble E ( A ij , P ( w ) , C ), butrather through our dynamic Jacobian family E ( A ij , Ω). Data in all panels are logarithmicallybinned . Details on the numerical calculation of J , log-binning and the construction of allmodel/real networks appear in Supplementary Section 4.14 a) (b)(c) (d)(e) (f) Figure 3:
The dynamic J ensemble in empirical networks . (a) Using the real socialnetwork, UCIonline, we extracted the Jacobian J from the random ensemble E ( A ij , P ( w ) , Ω)(left). The resulting J matrix has λ = 25 . >
0, and is, therefore unstable, as per May’sprediction. Implementing our Social dynamics (SIS) we constructed J again, this time using ourpredicted E ( A ij , Ω) (right). The resulting J is visibly different, owing to the emergent patternsof (7) and (1.4). Most importantly, it has λ = − . <
0, and hence despite the randommatrix prediction, the system is, in fact, stable. (b) Similar results are obtained from Epoch,with λ = 49 . > E ( A ij , P ( w ) , Ω) (left), vs. λ = − .
94 under E ( A ij , Ω) (right). (c)-(d)Ecological networks ECO1 and ECO2, analyzed via E ( A ij , P ( w ) , Ω) (left) vs. E ( A ij , Ω) (right),taking Ω from our Ecological dynamics. Once again we find that while the random J is unstable( λ = 186 and 227) once the dynamics is properly treated through Ω the system is found to bestable ( λ = − . − J matrices: random ( E ( A ij , P ( w ) , C ), left), Biochemical ( E ( A ij , Ω )), center)and Regulatory ( E ( A ij , Ω ), right). This exemplifies the fact that under the dynamics ensemble E ( A ij , Ω) the same network can lead to different dynamics behaviors, i.e . distinct J , owing tothe nonlinear interaction mechanisms (Biochemical vs. Regulatory). Also here we find thatwhile the random J is unstable ( λ = 6 . , . > J matrices all have λ < J each entry (pixel)represents an average over a 50 ×
50 block of the complete Jacobian matrix. Details on all 6empirical networks are provided in Supplementary Section 4.15 𝝀 Asymptotically unstable
𝑺 > 𝟎
Asymptotically stable
𝑺 < 𝟎
SocialRegulatoryEcologicalBiochemical 𝑱 𝑱 𝒊𝒋 − 𝑪 𝑱 𝒊𝒊 Positive feedback Negative feedback (a) (b)
𝑺 = 𝟎 𝝀 𝝀
𝑪 𝑪 𝑪 (c) (d)
𝑵 𝑵 𝑵 (f) (g)(i) (j) S e n s i t i v e 𝑺 = 𝑪 𝑺 = 𝟎 𝑺 = 𝟎 𝑺 = 𝟎𝑺 = −𝟏.𝟐 𝑺 = 𝟏.𝟖 𝝀 = 𝟎𝝀 = 𝟎𝝀 = 𝟎 𝑪 𝑪 𝑪 𝑱 𝐀𝐒 𝑱 𝐒𝐒 𝑱 𝐀𝐔 𝑵 𝑵 𝑵 𝑪 Scale-free Scale-free Scale-freeErdős-Rényi Erdős-Rényi Erdős-Rényi (e)(h)(k) 𝝀 𝑪 𝑪 𝑪 𝑪 Figure 4:
Three classes of dynamic stability in real and model networks . We extracted2 ,
077 Jacobian matrices from the E ( A ij , Ω) ensemble, using combinations of model/real networkswith different dynamic exponents Ω. For each J we calculated the principal eigenvalue λ and thestability classifier S in (3.31). (a) λ vs. S for all 2 , J -matrices. We observe the three predictedclasses: Asymptotically unstable (red) in which S > λ >
0; Sensitively stable(green), where S = 0 and λ can be both positive or negative; Asymptotically stable (blue),where S < λ <
0. Our classification was inaccurate in ∼
4% of the cases(grey). Empirical networks (solid symbols): Under the random ensemble E ( A ij , P ( w ) , C ) allour empirical networks are classified as unstable (red symbols). However, the same networksunder E ( A ij , Ω), each with its relevant Ω, become dynamically stable (blue symbols). Hence,networks rendered unstable via the existing paradigm, are, in fact, stable, when accountingfor the nonlinearity through Ω. (b) The value of λ emerges from the competition between J ’soff-diagonal terms, representing positive feedback, and the strength of the diagonal terms J ii (negative feedback). Therefore one can force a system towards stability ( λ <
0) by increasing C in (7). (c) - (e) Taking three specific J matrices from each class, we plot λ vs. C , seekingthe critical C , in which λ becomes negative (grey lines). (f) For the stable J AS ( S = − .
2) wefind that C decreases with N (squares), capturing the asymptotic stability, in which for largesystems ( N → ∞ ) stability is sustained even under arbitrarily small C . The theoretical scalingpredicted in Eq. (11) is also shown (solid line, slope − . J SS we have S = 0, thecritical C is independent of N , hence the system’s stability can be affected by finite changesto its dynamic parameters, e.g ., environmental perturbations. (h) The asymptotically unstable J AU ( S = 1 .
8) has C → ∞ in the limit of large N , in perfect agreement with Eq. (11) (solidline). (i) - (k) For a bounded P ( k ), e.g ., Erd˝os-R´enyi, β vanishes and hence S = 0 in (3.31).Under these conditions, regardless of Ω, the system is always sensitively stable and therefore C does not scale with N . This demonstrates the role of degree-heterogeneity in ensuring stability,in the face of changing environmental conditions.16 𝒊𝒋 ∼ 𝒌 𝒊𝝂 𝒌 𝒋𝝆 − 𝑪 𝑱 𝒊𝒊 ∼ −𝑪𝜿 𝐧𝐧 𝜼 𝒌 𝝁 Positive feedback Negative feedback
𝑺 = 𝜷 𝟏 + 𝝂 + 𝝆 − 𝝁 − 𝜼𝑨 𝒊𝒋 Degree- heterogeneity (a)(b)
Asymptotically unstableAsymptotically stable 𝑺 −𝟒𝜷 Figure 5:
The ingredients of dynamic stability . (a) The stability classifier S (3.31) com-prises five exponents: β characterizes the network topology and its degree heterogeneity via(3.2); ν, ρ determine the magnitude of J ’s off-diagonal terms; µ, η characterize the self-dynamicdiagonal terms. This portrays stability as a balance between the positive feedback mediatedby the system’s interactions, strengthened with increasing ν, ρ , vs. the stabilizing effect of thediagonal terms, that are reinforced by η, µ . (b) S vs. the exponents Ω = ( η, µ, ν, ρ ). Increasing ν, ρ drives the system towards the unstable class ( S > µ, η enhances thesystem’s stability. In addition to the discrete classification of stable dynamics (blue,
S <
S >
0) and sensitive dynamics (green, S = 0), S also helps quantify thecontinuum level of stability, where the more negative/positive is S , the more stable/unstable isthe system. Grey lines represent equi-stable systems.17 eferences [1] S.V. Buldyrev, R. Parshani, G. Paul, H.E. Stanley and S. Havlin. Catastrophic cascade of failures ininterdependent networks. Nature , 464:1025–1028, 2010.[2] I. Dobson, B.A. Carreras, V.E. Lynch and D.E. Newman. Complex systems analysis of series of blackouts:Cascading failure, critical points, and self-organization.
Chaos , 17:026103, 2007.[3] D. Duan, C. Lv, S. Si, Z. Wang, D. Li, J. Gao, S. Havlin, H.E. Stanley and S. Boccaletti. Universal behaviorof cascading failures in interdependent networks.
Proc. Natl. Acad. Sci. USA , 116:22452–57, 2019.[4] A.E. Motter and Y.-C. Lai. Cascade-based attacks on complex networks.
Physical Review E , 66:065102,2002.[5] P. Crucitti, V. Latora and M. Marchiori. Model for cascading failures in complex networks.
Phys. Rev. E ,69:045104–7, 2004.[6] D. Achlioptas, R.M. D’Souza and J. Spencer. Explosive percolation in random networks.
Science , 323:1453–1455, 2009.[7] J. Gao, B. Barzel and A.-L. Barab´asi. Universal resilience patterns in complex networks.
Nature , 530:307–312,2016.[8] S. Boccaletti, J.A. Almendral, S. Guana, I. Leyvad, Z. Liua, I.S. Nadal, Z. Wang and Y. Zou. Explosivetransitions in complex networks structure and dynamics: Percolation and synchronization.
Physics Reports ,660:1–94, 2016.[9] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. Hwang. Complex networks: Structure anddynamics.
Physics Reports , 424:175–308, 2006.[10] Michael M Danziger, Ivan Bonamassa, Stefano Boccaletti, and Shlomo Havlin. Dynamic interdependenceand competition in multilayer networks.
Nature Physics , 15(2):178–185, 2019.[11] K.Z. Coyte1, J. Schluter1 and K.R. Foster1. The ecology of the microbiome: Networks, competition, andstability.
Science , 350:663–666, 2015.[12] D. Tilman, P.B. Reich and J.M.H. Knops. Biodiversity and ecosystem stability in a decade-long grasslandexperiment.
Nature , 441:629632, 2006.[13] R.V. Sol´e and J.M. Montoya. Complexity and fragility in ecological networks.
Proc. R. Soc. London Ser. B ,268:2039–2045, 2001.[14] H.I. Schreier, Y. Soen and N. Brenner. Exploratory adaptation in large random networks.
Nature Commu-nications , 8:14826, 2017.[15] L.M. Pecora and T.L. Carroll. Master stability functions for synchronized coupled systems.
Phys. Rev. Lett. ,80:2109, 1998.[16] R.M. May. Will a large complex system be stable?
Nature , 238:413 – 414, 1972.[17] G. Caldarelli.
Scale-free networks: complex webs in nature and technology . Oxfrod University Press, NewYork, 2007.[18] A.-L. Barab´asi, H. Jeong, E. Ravasz, Z. N´eda, A. Schubert and T. Vicsek.
Physica A , 311:590, 2002.[19] D.J. Watts and S.H. Strogatz. Collective dynamics of ’small-world’ networks.
Nature , 393:440–442, 1998.[20] B. Barzel and O. Biham. Quantifying the connectivity of a network: The network correlation functionmethod.
Phys. Rev. E , 80:046104–15, 2009.[21] A.-L. Barab´asi and R. Albert. Emergence of scaling in random networks.
Science , 286:509–512, 1999.[22] M.E.J. Newman.
Networks - an introduction . Oxford University Press, New York, 2010.[23] S. Allesina and S. Tang. Stability criteria for complex ecosystems.
Nature , 483:205208, 2012.[24] W. Tarnowski, I. Neri and P. Vivo. Universal transient behavior in large dynamical systems on networks.
Phys. Rev. Research , 2:023333, 2020.[25] S. Suweis, F. Simini, J.R. Banavar and A. Maritan. Emergence of structural and dynamical properties ofecological mutualistic networks.
Nature , 500(7463):449–452, 2013.[26] W. Feng and K. Takemoto. Heterogeneity in ecological mutualistic networks dominantly determines com-munity stability.
Scientific reports , 4:5912, 2014.
27] A. Mougi and M. Kondoh. Diversity of interaction types and ecological community stability.
Science ,337(6092):349–351, 2012.[28] T. Okuyama and J.N. Holland. Network structural properties mediate the stability of mutualistic commu-nities.
Ecology Letters , 11(3):208–216, 2008.[29] E. Th´ebault and C. Fontaine. Stability of ecological communities and the architecture of mutualistic andtrophic networks.
Science , 329(5993):853–856, 2010.[30] D. Gravel, F. Massol and M.A. Leibold. Stability and complexity in model meta–ecosystems.
NatureCommunications , 7(12457):1–8, 2016.[31] Sitabhra Sinha and Sudeshna Sinha. Evidence of universality for the may-wigner stability theorem for randomnetworks with local dynamics.
Physical Review E , 71(2):020902, 2005.[32] G. Yan, N.D. Martinez and Y.-Y. Liu. Degree heterogeneity and stability of ecological networks.
Journal ofThe Royal Society Interface , 14:131, 2017.[33] C. Jacquet, C. Moritz, L. Morissette, P. Legagneux, F. Massol, P. Archambault and D. Gravel. No complex-itystability relationship in empirical ecosystems.
Nature Communications , 7:12573, 2016.[34] K.S. McCann. The diversity–stability debate.
Nature , 405(6783):228–233, 2000.[35] B. Barzel and A.-L. Barab´asi. Universality in network dynamics.
Nature Physics , 9:673 – 681, 2013.[36] U. Harush and B. Barzel. Dynamic patterns of information flow in complex networks.
Nature Communica-tions , 8:2181, 2017.[37] C. Hens, U. Harush, R. Cohen, S. Haber and B. Barzel . Spatiotemporal propagation of signals in complexnetworks.
Nature Physics , 15:403, 2019.[38] R. Pastor-Satorras, C. Castellano, P. Van Mieghem and A. Vespignani. Epidemic processes in complexnetworks.
Rev. Mod. Phys. , 87:925–958, 2015.[39] P.S. Dodds, R. Muhamad and D.J. Watts. An experimental study of search in global social networks.
Science ,301:827–829, 2003.[40] D. Brockmann, V. David and A.M. Gallardo. Human mobility and spatial disease dynamics.
Reviews ofNonlinear Dynamics and Complexity , 2:1, 2009.[41] G. Karlebach and R. Shamir. Modelling and analysis of gene regulatory networks.
Nature Reviews , 9:770–780,2008.[42] J.D. Murray.
Mathematical Biology . Springer, Berlin, 1989.[43] B. Barzel and O. Biham. Binomial moment equations for stochastic reaction systems.
Phys. Rev. Lett. ,106:150602–5, 2011.[44] B. Barzel and O. Biham. Stochastic analysis of complex reaction networks using binomial moment equations.
Phys. Rev. E , 86:031126, 2012.[45] C.S. Holling. Some characteristics of simple types of predation and parasitism.
The Canadian Entomologist ,91:385–398, 1959.[46] J.N. Holland, D.L. DeAngelis and J.L. Bronstein. Population dynamics and mutualism: functional responsesof benefits and costs.
American Naturalist , 159:231–244, 2002.[47] E.L. Berlow, J.A. Dunne, N.D. Martinez, P.B. Stark, R.J. Williams and U. Brose. Simple prediction ofinteraction strengths in complex food webs.
Proceedings of the National Academy of Sciences , 106:187–191,2009.[48] J.F. Hayes and T.V.J. Ganesh Babu.
Modeling and Analysis of Telecommunications Networks . John Wiley& Sons, Inc., Hoboken, NJ, USA, 2004.[49] A.-L. Barab´asi and Z.N. Oltvai. Network biology: understanding the cell’s functional organization.
Nat.Rev. Gen. , 5:101, 2004.[50] T. Opsahl and P. Panzarasa. Clustering in weighted networks.
Social Networks , 31:155–163, 2009.[51] J.-P. Eckmann, E. Moses and D. Sergi. Entropy of dialogues creates coherent structures in e-mail traffic.
Proc. Natl. Acad. Sci. USA , 101:14333–7, 2004.[52] M. Faloutsos, P. Faloutsos and C. Faloutsos. On power-law relationships of the Internet topology.
Compute.Commun. Rev. , 29:251–262, 1999.
53] G. Caldarelli and M. Catanzaro.
Networks: A very short introduction . Oxfrod University Press, New York,2012.[54] M. Bogu˜n´a, R. Pastor-Satorras and A. Vespignani. Absence of epidemic threshold in scale-free networks withdegree correlations.
Phys. Rev. Lett. , 90:028701–4, 2003.[55] R.M. May.
Theoretical Ecology: Principles and Application , 1:78–104, 1976.[56] J.A. Almendral and A. D´ıaz-Guilera. Dynamical and spectral properties of complex networks.
New Journalof Physics , 9:187–187, 2007.[57] P. Van Mieghem. Epidemic phase transition of the SIS type in network.
Europhysics Letters , 97(4):48004,2012.[58] P. Van Mieghem.
Graph Spectra for Complex Networks . Cambridge University Press, Cambridge, UK, 2010.[59] V.S. Afraimovich and L.A. Bunimovich. Dynamical networks: interplay of topology, interactions and localdynamics.
Nonlinearity , 20:1761–1771, 2007.[60] C. Kirst, M. Timme and D.Battaglia. Dynamic information routing in complex networks.
Nature Commu-nications , 7:11061, 2016.[61] A. Barrat, M. Barth´elemy and A. Vespignani.
Dynamical Processes on Complex Networks . CambridgeUniversity Press, Cambridge, 2008.[62] G. Yan, G. Tsekenis, B. Barzel, J.-J. Slotine, Y.-Y. Liu and A.-L. Barab´asi. Spectrum of controlling andobserving complex networks.
Nature Physics , 11:779786, 2015.[63] L. Schmetterer and K. Sigmund (Eds.).
Hans Hahn Gesammelte Abhandlungen Band 1/Hans Hahn CollectedWorks Volume 1 . Springer, Vienna, Austria, 1995.[64] H. Yu et al.
High-quality binary protein interaction map of the yeast interactome network.
Science , 322:104–110, 2008.[65] J.F. Rual et al . Towards a proteome-scale map of the human protein-protein interaction network.
Nature
Journalof the American Society for Information Science and Technology , 61:2417–2425, 2010. 𝒊𝒋 , 𝜷 𝑷(𝒘)
The random ensemble
𝔼(𝑨 𝒊𝒋 , 𝑷 𝒘 , 𝑪) 𝑨 𝒊𝒋 , 𝜷 𝛀 = (𝜼, 𝝁, 𝝂, 𝝆) 𝐝𝒙 𝒊 𝐝𝒕 = 𝑴 𝒙 𝒊 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝑴 𝒙 𝒊 𝑴 (𝒙 𝒋 ) The dynamic ensemble
𝔼(𝑨 𝒊𝒋 , 𝛀) 𝝀 ∼ 𝑵 𝑸 𝑺 𝝀 ∼ 𝑵 𝜷𝟐 𝜷𝟐 Social 𝐝𝒙 𝒊 𝐝𝒕 = −𝒇𝒙 𝒊 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒈 𝟏 − 𝒙 𝒊 𝒙 𝒋 𝛀 = (𝟎, 𝟏, −𝟏, 𝟎) 𝑺 = −𝜷
Regulatory 𝐝𝒙 𝒊 𝐝𝒕 = −𝑩𝒙 𝒊𝒂 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒋𝒉 𝒋𝒉 𝛀 = 𝟎, 𝒂 − 𝟏𝒂 , 𝟎, − 𝒉 + 𝟏𝒂 𝑺 = − 𝒉𝒂 𝜷
Ecological 𝐝𝒙 𝒊 𝐝𝒕 = 𝑩𝒙 𝒊 𝒊 𝑲 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝑭(𝒙 𝒋 ) 𝛀 = (𝟎, 𝟏, 𝟏, −𝟐) 𝑺 = −𝜷 Biochemical * 𝐝𝒙 𝒊 𝐝𝒕 = 𝑭 − 𝑾𝒙 𝒊 − ෩𝑩 𝒋=𝟏 𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝒙 𝒋 𝛀 = −𝟏, 𝟏, −𝟏, 𝟎 𝑺 = −𝜷
𝑺 > 𝟎 →
Unstable
𝑺 = 𝟎 →
Sensitive
𝑺 < 𝟎 →
StableUnstable (a) (b)(d) (e) (c) (f) (g) 𝑱 𝒊𝒋 𝑱 𝒊𝒋 Figure 1 - full size .21 oc i a l R e g u l a t o r y E co l o g i c a l B i oc h e m i c a l (a) (b)(c) (d)(e) (f)(g) (h) Figure 2 - full size .22 a) (b)(c) (d)(e) (f)
Figure 3 - full size .23 𝝀 Asymptotically unstable
𝑺 > 𝟎
Asymptotically stable
𝑺 < 𝟎
SocialRegulatoryEcologicalBiochemical 𝑱 𝑱 𝒊𝒋 − 𝑪 𝑱 𝒊𝒊 Positive feedback Negative feedback (a) (b)
𝑺 = 𝟎 𝝀 𝝀
𝑪 𝑪 𝑪 (c) (d)
𝑵 𝑵 𝑵 (f) (g)(i) (j) S e n s i t i v e 𝑺 = 𝑪 𝑺 = 𝟎 𝑺 = 𝟎 𝑺 = 𝟎𝑺 = −𝟏. 𝟐 𝑺 = 𝟏. 𝟖 𝝀 = 𝟎𝝀 = 𝟎𝝀 = 𝟎 𝑪 𝑪 𝑪 𝑱 𝐀𝐒 𝑱 𝐒𝐒 𝑱 𝐀𝐔 𝑵 𝑵 𝑵 𝑪 Scale-free Scale-free Scale-freeErdős-Rényi Erdős-Rényi Erdős-Rényi (e)(h)(k) 𝝀 𝑪 𝑪 𝑪 𝑪 Figure 4 - full size .24 𝒊𝒋 ∼ 𝒌 𝒊𝝂 𝒌 𝒋𝝆 − 𝑪 𝑱 𝒊𝒊 ∼ −𝑪𝜿 𝐧𝐧 𝜼 𝒌 𝝁 Positive feedback Negative feedback
𝑺 = 𝜷 𝟏 + 𝝂 + 𝝆 − 𝝁 − 𝜼𝑨 𝒊𝒋 Degree- heterogeneity (a)(b)
Asymptotically unstableAsymptotically stable 𝑺 −𝟒𝜷 Figure 5 - full size .25 ynamic stability of complex networksSupplementary information
August 18, 2020 ontents E ( A ij , Ω) x ∗ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Evaluating (cid:104) M ( x ) (cid:105) (cid:12) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Diagonal terms D ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 The off-diagonal terms W ij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 Piecing together the J ij puzzle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.4.1 Impact of P ( k ) and κ nn . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 J ij ∈ E ( A ij , Ω) S under negative interactions . . . . . . . . . . . . . . . . . . . 24 J ij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.3 Logarithmic binning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.4 Model and empirical networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 The Jacobian ensemble E ( A ij , Ω) To analyze complex system stability we seek the real
Jacobian matrix J ij , as extracted from thesystem’s specific nonlinear interaction mechanisms. We first rewrite Eq. (6) of the main text asd x i d t = F i (cid:0) x ( t ) (cid:1) , (1.1)where F i (cid:0) x ( t ) (cid:1) = M (cid:0) x i ( t ) (cid:1) + N (cid:88) j =1 A ij M (cid:0) x i ( t ) (cid:1) M (cid:0) x j ( t ) (cid:1) . (1.2)We denote (1.1)’s fixed-point(s) by x ∗ = ( x , . . . x N ) (cid:62) , where we omit the the term t to capturetheir independence on time. These fixed-points are obtained by solving the equilibrium equation F i ( x ∗ ) = 0 . (1.3)To assess the dynamic stability of each of (1.3)’s solutions we track their response to smallperturbations, via the Jacobian J ij = ∂F i ( x ) ∂x j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ , (1.4)whose structure we obtain below. Writing J = A ⊗ W − I ⊗ D (1.5)we treat separately the diagonal terms J ii = D ii and the off-diagonal terms J ij = A ij W ij ( I isthe identity matrix and ⊗ is the Hadamard product). ∗ We consider systems of the form (1.1) that exhibit at least one fully positive fixed-point x ∗ .Using Eq. (1.3) we write M ( x i ) + N (cid:88) j =1 A ij M ( x i ) M ( x j ) = 0 , (1.6)which can be expressed as R ( x i ) N (cid:88) j =1 A ij M ( x j ) = 1 , (1.7)where 1 ( x i ) = − M ( x i ) M ( x i ) . (1.8)Denoting i ’s degree by k i = (cid:80) Nj =1 A ij we write (1.7) as R ( x i ) = 1 k i (cid:104) M ( x ) (cid:105) i, (cid:12) , (1.9)in which (cid:104) M ( x ) (cid:105) i, (cid:12) = 1 k i N (cid:88) j =1 A ij M ( x j ) (1.10)represents an average of the function M ( x j ) over all of node i ’s direct network neighborhood ,which we captured by the (cid:12) symbol. Assuming the function R ( x ) is invertible around x i weobtain i ’s fixed-point activity as x i = R − ( q i ) , (1.11)where R − ( x ) is the inverse function of R ( x ) and q i = 1 k i (cid:104) M ( x ) (cid:105) i, (cid:12) (1.12)is i ’s inverse degree, which approaches zero in the limit of large k i .To understand the dependence of activity on degree, we seek the average steady-state activityover all nodes with degree k as x ( k ) = 1 | G ( k ) | (cid:88) i ∈ G ( k ) x i , (1.13)where G ( k ) = { i ∈ [1 , N ] | k i = k } is the group of all nodes with degree k , and | G ( k ) | is thenumber of nodes in that group. To approximate x ( k ) we first evaluate the average neighborhoodactivity of all nodes in G ( k ) via (cid:104) M ( x ) (cid:105) k, (cid:12) = 1 | G ( k ) | (cid:88) i ∈ G ( k ) (cid:104) M ( x ) (cid:105) i, (cid:12) , (1.14)capturing the average value of M ( x ) surrounding nodes with degree k . Most generally, theactivities x j may depend, to some extent, on the degree k i of their neighbor i , however thisdependence is often unknown, and, most importantly, relatively weak . Indeed, i is but one of j ’s k j neighbors, who are typically a random selection, representative of the network’s degreedistribution . Therefore, i ’s degree k i is only weakly correlated with j ’s activity x j . To capturethis dependence we approximate (1.14) as (cid:104) M ( x ) (cid:105) k, (cid:12) ≈ f ( k ) (cid:104) M ( x ) (cid:105) (cid:12) , (1.15)where 2 M ( x ) (cid:105) (cid:12) = 1 N N (cid:88) i =1 k i N (cid:88) j =1 A ij M ( x j ) (1.16)represents the average of the function M ( x j ) over all network neighborhoods, independent of i or k . In (1.15) the term (cid:104) M ( x ) (cid:105) (cid:12) captures a characteristic aggregated function of the system’sfixed-point, capturing the value of M ( x j ) averaged over all nearest neighbor nodes, but notspecifically over i ’s neighbors - hence the approximation sign ≈ , rather than =. The function f ( k ) corrects this average to account for the fact that in in the l.h.s. of (1.15) the average iscarried out selectively over nodes whose neighbor has degree k . For example, f ( k ) = 1 representsthe perfect case of no degree-correlations, under which x j is fully independent of its neighbor’sdegree k . More generally, we assume f ( k ) to be a weak function, i.e . sub-polynomial in k , forexample f ( k ) ∼ log n ( k ) . (1.17)We can now approximate x ( k ) in (1.13), taking x i from (1.11), and (cid:104) M ( x ) (cid:105) i, (cid:12) from (1.15) toobtain x ( k ) = R − ( q ) , (1.18)where q = 1 f ( k ) (cid:104) M ( x ) (cid:105) (cid:12) k . (1.19)In the asymptotic limit of large k , the sub-polynomial contribution of f ( k ) becomes negligible,and hence we have q ∼ (cid:104) M ( x ) (cid:105) − (cid:12) k − . (1.20) (cid:104) M ( x ) (cid:105) (cid:12) To link (cid:104) M ( x ) (cid:105) (cid:12) to the network A ij we use a mean-field approximation. We first define thenearest neighbor activity as x nn = (cid:104) x (cid:105) (cid:12) = 1 N N (cid:88) i =1 k i N (cid:88) j =1 A ij x j , (1.21)and its corresponding nearest neighbor degree as κ nn = (cid:104) k (cid:105) (cid:12) = 1 N N (cid:88) i =1 k i N (cid:88) j =1 A ij k j . (1.22)Hence, the average nearest neighbor node is characterized by activity x nn and degree κ nn . While x nn depends on the system’s dynamics (1.1), κ nn is fully determined by the topology A ij , linkedto the network’s degree distribution P ( k ). In the absence of degree correlations we have nn = (cid:104) k (cid:105)(cid:104) k (cid:105) , (1.23)with (cid:104) k n (cid:105) capturing the n th moment of P ( k ). For a homogeneous network in which P ( k ) isbounded we have κ nn ≈ (cid:104) k (cid:105) , however, if the network is highly heterogeneous, i.e . P ( k ) is fat-tailed - as indeed, many real networks are - we find that κ nn (cid:29) (cid:104) k (cid:105) . More generally, in casethe network has degree correlations we write κ nn = ∞ (cid:88) k (cid:48) =1 ∞ (cid:88) k =1 kP ( k | k (cid:48) ) P ( k (cid:48) ) (1.24)in which the conditional probability P ( k | k (cid:48) ) captures the dependence between neighboring nodes.In both cases, (1.23) and the more general (1.24), under extreme degree-heterogeneity, κ maydiverge with system size as κ nn ∼ N β . (1.25)We can now link (cid:104) M ( x ) (cid:105) (cid:12) in (1.16) to κ nn and x nn , as, indeed, (cid:104) M ( x ) (cid:105) (cid:12) is designed to capturethe mean value of M ( x ) over all nearest neighbor nodes, namely nodes with typical degree κ nn and activity x nn . We, therefore, approximate (1.16) as (cid:104) M ( x ) (cid:105) (cid:12) ≈ M ( x nn ) , (1.26)in which we swap the order of operations from (cid:104) M ( x ) (cid:105) to M ( (cid:104) x (cid:105) ). Such approximation becomesexact in the limit where M ( x ) is linear or when x j are narrowly distributed. For example, in ourBiochemical model we have M ( x j ) = x j , a linear function, and in our Social model the activities x j are bounded, hence narrowly distributed between zero and unity. Approximation (1.26) isalso justified when M ( x ) is sub-linear, e.g ., a saturating function M ( x → ∞ ) →
1. Suchsaturation is present, for instance, in our Regulatory and Ecological dynamics. Under theseconditions, even if x j are broadly distributed, M ( x j ) are uniform, and hence their averageis (roughly) concentrated around M ( (cid:104) x (cid:105) ). For a super-linear M ( x ), combined with broadlydistributed x j , approximation (1.26) becomes inaccurate.To evaluate M ( x nn )’s dependence on κ nn we use (1.9) to write R ( x nn ) = 1 κ nn f ( κ nn ) M ( x nn ) , (1.27)where we substitute (cid:104) M ( x ) (cid:105) i, (cid:12) by f ( κ nn ) M ( x nn ) following approximations (1.15) and (1.26),and use the nearest neighbor degree κ nn in place of k . We therefore arrive at Z ( x nn ) = 1 f ( κ nn ) κ nn ≡ q (cid:12) , (1.28)where Z ( x ) = R ( x ) M ( x ) is a dynamic function, fully determined by M ( x ) , M ( x ) and M ( x )in (1.1). In (1.29) we use the notation q (cid:12) to signify the inverse degree ( q ) associated with thetypical neighborhood ( (cid:12) ). By inversion we obtain4 nn = Z − ( q (cid:12) ) , (1.29)and hence M ( x nn ) = M (cid:0) Z − ( q (cid:12) ) (cid:1) . (1.30)To obtain the asymptotic scaling of M ( x nn ) on κ nn we use the Hahn expansion to express Z ( x ) in the form of a power series around q (cid:12) → i.e . κ nn → ∞ . Hence we write M (cid:0) Z − ( q (cid:12) ) (cid:1) = ∞ (cid:88) n =0 G n q Ψ( n ) (cid:12) , (1.31)where Ψ( n ) is a set of real powers in ascending order with n . In the limit of large κ nn (small q (cid:12) ) the expansion in (1.31) is dominated by the leading power Ψ(0), predicting that M ( x nn ) ∼ κ ξ nn , (1.32)where ξ = − Ψ(0) . (1.33)The exponent ξ is fully determined by the dynamic functions M ( x ) , M ( x ) and M ( x ) throughthe Hahn expansion in (1.31). It is independent of A ij , as well as of other microscopic parameters,such as the specific rate-constants used in (1.1), which are encapsulated within the coefficients G n , but have no impact on the powers Ψ( n ).Using approximation (1.26) we substitute M ( x nn ) in (1.32), to replace the (cid:104) M ( x ) (cid:105) (cid:12) term inthe denominator of Eq. (1.19), obtaining q ∼ κ − ξ nn k − . (1.34)Next we use these results, specifically (1.18) and (1.34), to obtain the terms of the Jacobian, J ij and J ii , and their dependence on κ nn and k . D ii Using (1.2) and (1.4) we write the diagonal Jacobian terms as D ii = ∂F i ( x ) ∂x i (cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ = M (cid:48) ( x i ) + M (cid:48) ( x i ) N (cid:88) j =1 A ij M ( x j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ , (1.35)where M (cid:48) ( x ) = ∂M /∂x and M (cid:48) ( x ) = ∂M /∂x . Next we use (1.10) to express the sum on ther.h.s. obtaining D ii = M (cid:48) ( x i ) (cid:12)(cid:12)(cid:12) x = x ∗ + M (cid:48) ( x i ) k i (cid:104) M ( x ) (cid:105) i, (cid:12) (cid:12)(cid:12)(cid:12) x = x ∗ , (1.36)5n which we substitute the summation by the neighbor average (cid:104) M ( x ) (cid:105) i, (cid:12) . Finally, expressingthe fixed-point x ∗ via (1.11) and averaging over all nodes with degree k as done in (1.13) - (1.19),we arrive at D ( k ) = M (cid:48) (cid:0) R − ( q ) (cid:1) + kM (cid:48) (cid:0) R − ( q ) (cid:1) f ( k ) (cid:104) M ( x ) (cid:105) (cid:12) , (1.37)capturing the average diagonal term D ii , associated with a node i ∈ G ( k ).Next we use (1.8) to write M ( x ) = − M ( x ) /R ( x ), which provides M (cid:48) ( x ) = − M (cid:48) ( x ) R ( x ) + M ( x ) R (cid:48) ( x ) R ( x ) . (1.38)Substituting (1.38) and (1.19) into (1.37) we arrive at D ( k ) = − M (cid:48) (cid:0) R − ( q ) (cid:1) q + M (cid:0) R − ( q ) (cid:1) R (cid:48) (cid:0) R − ( q ) (cid:1) q + M (cid:48) (cid:0) R − ( q ) (cid:1) q , (1.39)where we used the fact that R ( R − ( q )) = q . Collecting all terms we write D ( k ) = 1 q Y (cid:0) R − ( q ) (cid:1) , (1.40)where Y ( x ) = M ( x ) R (cid:48) ( x ), similar to Z ( x ) above, is a dynamic function fully determined by M ( x ) and M ( x ), independent of A ij .To obtain the asymptotic scaling of D ( k ) with k we, once again, use the Hahn expansion toexpress Y ( x ) in the form of a power series around q → i.e . k → ∞ . Writing Y (cid:0) R − ( q ) (cid:1) = ∞ (cid:88) n =0 C n q Γ( n ) , (1.41)we obtain the leading power in the expansion of D ( k ) as D ( k ) ∼ q − , exact up to higherpowers in q , which vanish under q →
0. Using the scaling of Eq. (1.34), this predicts D ( k ) ∼ κ ξµ nn k µ , (1.42)where µ = 2 − Γ(0) , (1.43)and ξ is taken from (1.33).Equation (1.42), valid in the limit of large κ nn and k , captures the typical magnitude of thediagonal terms D ii , and their dependence on each node i ’s degree k i . The scaling exponents µ and ξ are fully determined by Γ( n ) in (1.41) and Ψ( n ) in (1.31), namely they are affected by the powers of the dynamic functions M ( x ) and M ( x ). At the same time µ and ξ are independent ofthe coefficients C n , G n , hence they are not sensitive to specific parameters, such as rate constantsin (1.1). While Eq. (1.42) provides the scaling with k , D ( k )’s specific magnitude may, generally,depend on the coefficients, however for sufficiently large k , such dependencies have little impact6n the spectrum of J ij . Finally, the sign of D ( k ), positive or negative, is determined by Y ( x ),which, in turn depends on M ( x ) and M ( x ). Stability can only be ensured in case D ( k ) < i.e . negative self-feedback, and hence below we write D ( k ) ∼ − Cκ ξµ nn k µ , (1.44)where C > W ij The off-diagonal Jacobian terms are given by W ij = ∂F i ( x ) ∂x j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ = ∂∂x j M ( x i ) + M ( x i ) N (cid:88) j =1 A ij M ( x j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ . (1.45)Collecting only the terms that explicitly depend on x j we obtain W ij = M ( x i ) A ij M (cid:48) ( x j ) (cid:12)(cid:12)(cid:12) x = x ∗ , (1.46)which, taking the fixed-point activities from (1.11) provides W ij = M (cid:0) R − ( q i ) (cid:1) A ij M (cid:48) (cid:0) R − ( q j ) (cid:1) . (1.47)We seek the average dependence of W ij on the degrees k i and k j , hence we consider the group G ( k , k ) = { ( i, j ) | k i = k ; k j = 2; A ij = 1 } , comprising all links ( A ij = 1) linking nodes withdegrees k , k . We then average over the relevant Jacobian term as W ( k , k ) = 1 | G ( k , k ) | (cid:88) ( i,j ) ∈ G ( k ,k ) W ij , (1.48)allowing us to evaluate the magnitude of an average (non zero) term whose row index is associatedwith a node of degree k , and whose column index is associated with a node of degree k . Using(1.46) we write W ( k , k ) = 1 | G ( k , k ) | (cid:88) ( i,j ) ∈ G ( k ,k ) M ( x i ) A ij M (cid:48) ( x j ) = (cid:104) M ( x i ) M (cid:48) ( x j ) (cid:105) G ( k ,k ) , (1.49)a product average over all node pairs in G ( k , k ). In case the two activities x i and x j areuncorrelated, we can approximate the r.h.s. of (1.49) by a product over the single averages, i.e . (cid:104) M ( x i ) (cid:105) G ( k ,k ) (cid:104) M (cid:48) ( x j ) (cid:105) G ( k ,k ) . Each of these averages can be then evaluated by (cid:104) M ( x i ) (cid:105) G ( k ,k ) ≈ M (cid:0) R − ( q ) (cid:1) (1.50) (cid:104) M (cid:48) ( x j ) (cid:105) G ( k ,k ) ≈ M (cid:48) (cid:0) R − ( q ) (cid:1) , (1.51)where we use the fact that in G ( k , k ) the activity x i is on average R − ( q ), i.e. x ( k ) (1.13),7nd the activity x j is R − ( q ), namely x ( k ). In reality, however, the two activities are, to someextent, correlated by virtue of being nearest neighbors ( A ij = 1). We therefore follow a similartrack as in (1.15) and write (1.49) as W ( k , k ) ≈ f ( k , k ) M (cid:0) R − ( q ) (cid:1) M (cid:48) (cid:0) R − ( q ) (cid:1) , (1.52)capturing the k , k correlations via the sub-polynomial function f ( k , k ).We can now obtain the scaling of W ( k , k ) by examining the behavior of M (cid:0) R − ( q ) (cid:1) and M (cid:48) (cid:0) R − ( q ) (cid:1) in the limit of large k and k , or subsequently small q , q . We, therefore expressthese functions via the Hahn expansion as M (cid:0) R − ( q ) (cid:1) = ∞ (cid:88) n =0 K n q Π( n ) (1.53) M (cid:48) (cid:0) R − ( q ) (cid:1) = ∞ (cid:88) n =0 L n q Θ( n ) , (1.54)and take only the leading term in the limit q →
0, providing W ( k , k ) = f ( k , k ) K q Π(0)1 L q Θ(0)2 + · · · . (1.55)As we are only interested in the asymptotic scaling with k , k we neglect all terms that do notcontribute the scaling, namely f ( k , k ) , K and L . Finally, substituting q with κ − ξ nn k − as perEq. (1.34), we obtain W ( k , k ) ∼ κ ξ ( ν + ρ )nn k ν k ρ (1.56)where ν = − Π(0) , ρ = − Θ(0) , (1.57)and ξ is taken from (1.33). J ij puzzle We have now obtained the interplay between structure and dynamics, expressed through themagnitude of the diagonal/off-diagonal Jacobian entries, and their dependence on the network’sdegree distribution via k i , k j and κ nn . This dependence is given in terms of asymptotic scalingrelationships, as captured by the model dependent exponents µ, ν, ρ, ξ . This leaves a degree offreedom, determined by the specific parameters of each model to add an arbitrary coefficient asa pre-factor to the actual J ij terms. For example, the expression D ( k ) ∼ − Cκ ξµ nn k µ in Eq. (1.42)should be taken to mean D ( k ) = − C ( κ nn , k ) κ ξµ nn k µ + · · · (1.58)where C ( κ nn , k ) is a sub-polynomial function, which is roughly constant, i.e . C ( κ nn , k ) ≈ C , and8 · · represents lower powers of κ nn and k , that are negligible in the asymptotic limit. In thecontext of stability the diagonal terms of J ij are assumed to be negative, capturing a stabilizingnegative feedback, hence the − C term in (1.58), which is negative if C > J ij = A ij W ij + δ ij D ii = − δ ij Cκ ξµ nn k µi + A ij κ ξ ( ν + ρ )nn k νi k ρj , (1.59)where the first term on the r.h.s. represents the (negative) diagonal entries, and second termcaptures the off-diagonal entries, which are non-zero only if A ij = 1; δ ij is the Kronecker deltafunction. As we are only interested in the sign of the principle eigenvalue, but not in its specificmagnitude, we have the degree of freedom to multiply J ij by an arbitrary constant. We thereforenormalize J ij by κ − ξ ( ν + ρ )nn , providing J ii ∼ − Cκ η nn k µi (1.60) J ij ∼ A ij k νi k ρj (1.61)for the diagonal and off-diagonal terms respectively, with η = ξ ( µ − ν − ρ ) - recovering Eqs. (7)and (8) of the main text.Equations (1.60) and (1.61) describe the asymptotic structure of the diagonal and off-diagonalJacobian terms, as extracted from the general dynamics of Eq. (1.1). The resulting J ij ischaracterized by several distinct structural and dynamic inputs: A ij , the network topology,which determines the non-vanishing off-diagonal elements. It also determines the degrees k i , k j and the average neighbor degree κ nn . The exponentsΩ = ( η, µ, ν, ρ ) (1.62)are fully independent of the network topology, extracted from the dynamic functions M ( x ) , M ( x )and M ( x ), i.e . the interaction mechanisms driving Eq. (1.1). These exponents are universalin the sense that they do not depend of the specific model parameters , but rather on the model itself, e.g. , Social, Regulatory etc. The coefficient C , in contrast, is non-universal, andits value is determined by the specific rate-constants and time-scales driving Eq. (1.1); herewe do not attempt to predict the magnitude of this coefficient. In the limit of sufficientlylarge k i and k j , and, where applicable - in the limit of large κ nn , the specific finite value of C has negligible impact on the principle eigenvalue of J ij as we explicitly show in Sec. 3. Hencestability is asymptotically determined by the exponents Ω, independent of C . The meaning isthat the model can be asymptotically stable or unstable, regardless of its specific parameters . P ( k ) and κ nn In this derivation we considered only the scaling of J ij on the degrees k i , k j , and on the nearestneighbor degree κ nn . We disregarded coefficients, such as C n or G n in (1.41) and (1.31), or sub-polynomial functions like f ( k ) in (1.17). The rationale is that in the limit of large degrees theJacobian spectrum is mainly impacted by these asymptotic scaling relationships, and has littledependence on finite size coefficients that do not scale with k, κ nn , and hence also do not grow9ignificantly as a function of system size N . Indeed, May’s paradox is rooted precisely in suchscaling, in which the principle eigenvalue grows asymptotically as ∼ N c (in May’s work c = 1 / .Therefore, Eqs. (1.60) and (1.61) are especially relevant when P ( k ) is fat-tailed, e.g. , scale-free,where the hub-degrees and the nearest neighbor degree can potentially scale with N . The role of κ nn . While k i is a node specific attribute, that captures a specific dependencybetween the i th diagonal term and i ’s degree, the pre-factor κ η nn represents a network aggregatedparameter, indeed - a constant , whose impact is often negligible in the asymptotic limit of large k . We include it in our analysis, however, because under extreme degree-heterogeneity, we mayobserve that κ nn diverges as κ nn ∼ N β (1.25), and therefore can potentially impact the systemstability under N → ∞ . For example, in a scale-free network where P ( k ) ∼ k − γ , we can use κ nn = (cid:104) k (cid:105) / (cid:104) k (cid:105) in (1.23) to obtain κ nn ∼ N γ < N − γ ≤ γ < N ) γ ≥ , (1.63)which scales with N as long as γ <
3. There are, however, broad conditions, that arise quitenaturally in many real systems, in which the κ nn term in (1.60) can be neglected, significantlysimplifying the stability analysis: • Finite κ nn . In case γ ≥ κ nn no longer scales with N , it behaves as a constant and hasno impact in the asymptotic limit. Under these conditions it suffices to write Eq. (1.60)as J ii ∼ − Ck µi . • Bounded activities . In some models the activities x i are bounded. For example, inspreading processes, from epidemics to cascading failures, the activities satisfy 0 ≤ x i ≤ x ( k → ∞ ) ∼
1. The result is that the nearest neighbor activity x nn , associatedwith degree κ nn is itself bounded, and even if κ nn → ∞ `a la Eq. (1.63), we still have x nn ∼
1. Consequently M ( x nn ) in (1.30) also approaches a constant value, and thereforethe leading power in the expansion (1.31) is Ψ(0) = 0. This provides, based on Eq. (1.33), ξ = 0, which in turn leads to η = 0 in (1.60), again resulting in J ii ∼ − Ck µi , independentof κ nn . • Saturating M ( x ). Another common feature in many relevant models is that M ( x j →∞ ) →
1. This represents the saturating impact of node j on its nearest neighbor i , asfrequently observed in regulatory processes or in population dynamics. Once again, wehave M ( x nn ) ∼
1, providing ξ = 0, and consequently η = 0 in (1.60).10 he ensemble E ( A ij , Ω) summary Staring from Eq. (1.1) we use M ( x ) , M ( x ) and M ( x ) to construct the functions R ( x ) = − M ( x ) M ( x ) , Y ( x ) = M ( x ) R (cid:48) ( x ) , Z ( x ) = R ( x ) M ( x ) . (1.64)From these we extract the four relevant power-series expansions M (cid:0) Z − ( x ) (cid:1) = ∞ (cid:88) n =0 G n x Ψ( n ) , Y (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 C n x Γ( n ) ,M (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 K n x Π( n ) , M (cid:48) (cid:0) R − ( x ) (cid:1) = ∞ (cid:88) n =0 L n x Θ( n ) (1.65)whose leading powers determine the dynamic exponents Ω = ( η, µ, ν, ρ ) as µ = 2 − Γ(0) , ν = − Π(0) , ρ = − Θ(0) , η = − Ψ(0)( µ − ν − ρ ) . (1.66)To construct J ij ∈ E ( A ij , Ω) we first generate the network A ij , then extract the degrees k i ofall nodes and the nearest neighbor degree κ nn from Eq. (1.22). The resulting J ij satisfies J ii ∼ − Cκ η nn k µi (1.67) J ij ∼ A ij k νi k ρj , (1.68)where the constant C >
Analyzing the dynamic models
To demonstrate our formalism we examined several commonly used dynamic models, as listedin Fig. 1g of the main text. Below we extract the exponents Ω for each of these models.
As an example of Social dynamics we used epidemic spreading via the Susceptible-Infected-Susceptible (SIS) model, in which x i ( t ) captures the probability of infection of node i . Denotingthe susceptible state by S and the infected state by I , the model includes the following transitions I f −→ S (2.1) I + S g −→ I, (2.2)capturing the processes of recovery at a rate f and infection at rate g . This gives rise to thedynamic equation d x i d t = − f x i + N (cid:88) j =1 A ij g (1 − x i ) x j , (2.3)characterized by the dynamic functions M ( x ) = − f x, M ( x ) = 1 − x and M ( x ) = gx .To obtain the relevant J ij ensemble, we seek the exponents Ω = ( η, µ, ν, ρ ). We begin bytranslating the given functions M ( x ) , M ( x ) , M ( x ) in to the relevant functions shown in (1.64),providing R ( x ) = − M ( x ) M ( x ) = 1 − xf x (2.4) Y ( x ) = M ( x ) R (cid:48) ( x ) = x − f x (2.5) Z ( x ) = R ( x ) M ( x ) = gf − gf x. (2.6)Inverting R ( x ) and Z ( x ), we write R − ( x ) = 1 f x + 1 , (2.7) Z − ( x ) = 1 − fg x, (2.8)allowing us to construct the Hahn expansions of (1.65) as12 (cid:0) Z − ( x ) (cid:1) = gZ − ( x ) = g − f x (2.9) Y (cid:0) R − ( x ) (cid:1) = R − ( x ) − f × (cid:0) R − ( x ) (cid:1) = − x − f x (2.10) M (cid:0) R − ( x ) (cid:1) = 1 − R − ( x ) = f xf x + 1 = f x − f x + f x − · · · (2.11) M (cid:48) (cid:0) R − ( x ) (cid:1) = g. (2.12)Each of these functions is expressed as a Hahn power-series, in some cases a finite polynomial, e.g ., (2.9) or (2.10), and in others an infinite series, where we only write the leading terms around x →
0. Here, coincidentally, the last function (2.12) is a constant, with the only power being x . We can list the relevant powers in these Hahn expansions as Ψ(0) = 0 , Γ(0) = 1 , Π(0) = 1and Θ(0) = 0, which, using (1.66) provides µ = 2 − Γ(0) = 1 (2.13) ν = − Π(0) = − ρ = − Θ(0) = 0 (2.15) η = − Ψ(0)( µ − ν − ρ ) = 0 . (2.16) We used the Michaelis-Menten model to capture gene Regulatory dynamics . Here, Eq. (1.1)tracks the level of gene expression x i ( t ), as regulated via the sub-cellular network A ij , providingd x i d t = − Bx ai + N (cid:88) j =1 A ij x hj x hj . (2.17)The self-dynamic term M ( x ) = − Bx a captures biochemical processes , such as degradation( a = 1) or dimerization ( a = 2). The interaction terms M ( x ) = 1 , M ( x ) = x h / x h describegenetic activation, a switch-like dynamics, which ranges from M (0) = 0 to M ( x → ∞ ) = 1 tocapture activation of gene i by gene j . The Hill coefficient h governs the rate of saturation of M ( x ).First, we construct the three functions summarized in (1.64): R ( x ) = − M ( x ) M ( x ) = 1 Bx a (2.18) Y ( x ) = M ( x ) R (cid:48) ( x ) = − aBx a +1 (2.19) Z ( x ) = R ( x ) M ( x ) = x h B ( x a + x a + h ) . (2.20)Inverting R ( x ), we write 13 − ( x ) = B − a x − a , (2.21)allowing us to construct the Hahn expansions (1.65) Y (cid:0) R − ( x ) (cid:1) = − aB (cid:18) Bx (cid:19) − a +1 a = aB a x a +1 a (2.22) M (cid:0) R − ( x ) (cid:1) = 1 = x (2.23) M (cid:48) (cid:0) R − ( x ) (cid:1) = hB − h − a x − h − a (cid:16) B − ha x − ha (cid:17) = hB h +1 a x h +1 a − hB h +1 a x h +1 a + · · · . (2.24)In each of these expansions we write the leading terms in the x → i.e . a pure monomial, and in (2.24) we show the firsttwo terms of the relevant Hahn expansion. To obtain Ω we extract only the leading powers Γ(0) = ( a + 1) /a, Π(0) = 0 , Θ(0) = ( h + 1) /a , ignoring the coefficients , e.g ., aB /a or hB ( h +1) /a .We can now use (1.66) to extract the dynamic exponents as µ = 2 − Γ(0) = a − a (2.25) ν = − Π(0) = 0 (2.26) ρ = − Θ(0) = − h + 1 a . (2.27)To obtain η we must calculate M ( Z − ( x )), requiring us to invert the function Z ( x ) in (2.20),which becomes analytically prohibitive for general a and h . Fortunately, in Regulatory dynamicswe can calculate this term directly. We first recall, from Eq. (1.30) that M ( Z − ( q (cid:12) )) = M ( x nn ),where x nn is the average neighbor activity and q (cid:12) ∼ /κ nn is the average neighbor degree. Fora general node with degree k we can use (1.18) together with (2.21) to write x ( k ) = B − a q − a ∼ k a , (2.28)a positive scaling with k , predicting that x ( k → ∞ ) → ∞ , namely that a node’s activity increaseswith its degree. Consequently, under a sufficiently large κ nn , the nearest neighbor activity x nn also diverges as κ /a nn , or alternatively, as q − /a (cid:12) . Hence we write M ( x nn ) = x h nn x h nn ∼ q − ha (cid:12) q − ha (cid:12) = 11 + q ha (cid:12) = 1 − q ha (cid:12) + · · · , (2.29)observing directly that the leading power of M ( x nn ), and therefore also of M ( Z − ( x )) isΨ(0) = 0. As a result, we obtain, using (1.66) our final exponent η = − Ψ(0)( µ − ν − ρ ) = 0 . (2.30)This is, in fact, precisely the case of Bounded M ( x ), discussed in Sec. 1.4.1, for which we apriori predicted η = 0. 14 .3 Ecological dynamics We consider mutualistic eco-systems, such as plant-pollinator networks, in which the interactingspecies exhibit symbiotic relationships. The species populations follow the dynamic equationd x i d t = Bx i ( t ) (cid:18) − x i ( t ) K (cid:19) + N (cid:88) j =1 A ij x i ( t ) F (cid:0) x j ( t ) (cid:1) . (2.31)The self dynamics M ( x ) = Bx (cid:18) − xK (cid:19) (2.32)captures logistic growth: when the population is small, the species reproduce at a rate B , yet,as x i approaches the carrying capacity of the system K , growth is hindered by competition overlimited resources . The mutualistic inter-species interactions are captured by M ( x ) = x and M ( x ) = F ( x ), where F ( x ) represents the functional response , describing the positive impactthat species j has on species i . This functional response can take one of several forms : Type slowromancapi@: linear impact F ( x ) = x. (2.33) Type slowromancapii@: saturating impact F ( x ) = x x . (2.34) Type slowromancapiii@: a generalization of Type slowromancapii@, where F ( x ) = x h x h . (2.35)In our simulations we used Type slowromancapii@ mutualistic interactions, therefore M ( x ) = x x . (2.36)For simplicity, we set the coefficients to B = K = 1.For R ( x ) , Y ( x ) and Z ( x ) we have R ( x ) = − M ( x ) M ( x ) = 1 x − Y ( x ) = M ( x ) R (cid:48) ( x ) = − x ( x − (2.38) Z ( x ) = R ( x ) M ( x ) = xx − . (2.39)Inverting R ( x ) and Z ( x ), we write 15 ocial 𝐝𝒙 𝒊 𝐝𝒕 = −𝒇𝒙 𝒊 + 𝒋=𝟏 𝑵 𝑨 𝒊𝒋 𝒈 𝟏 − 𝒙 𝒊 𝒙 𝒋 𝛀 = (𝟎, 𝟏, −𝟏, 𝟎) 𝑺 = −𝜷
Regulatory 𝐝𝒙 𝒊 𝐝𝒕 = −𝑩𝒙 𝒊𝒂 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒋𝒉 𝒋 𝒉 𝛀 = 𝟎, 𝒂 − 𝟏𝒂 , 𝟎, − 𝒉 + 𝟏𝒂
𝑺 = − 𝒉𝒂 𝜷 Ecological 𝐝𝒙 𝒊 𝐝𝒕 = 𝑩𝒙 𝒊 𝒙 𝒊 𝑲 + 𝒋=𝟏𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝒙 𝒋 𝒋 𝛀 = (𝟎, 𝟏, 𝟏, −𝟐)
𝑺 = −𝜷
Biochemical 𝐝𝒙 𝒊 𝐝𝒕 = 𝑭 − 𝑾𝒙 𝒊 − ෩𝑩 𝒋=𝟏 𝑵 𝑨 𝒊𝒋 𝒙 𝒊 𝒙 𝒋 𝛀 = −𝟏, 𝟏, −𝟏, 𝟎 𝑺
𝐍𝐞𝐠 = −𝜷
Figure 6:
Dynamic models - summary . Our four dynamic models and their associatedexponents Ω. For each model we also show the stability classifier S . For Biochemical we use S Neg , see Sec. 3.2 R − ( x ) = x + 1 x , (2.40) Z − ( x ) = 12 x (cid:16) (cid:112) x (cid:17) , (2.41)where in Z − ( x ) we choose the positive solution, corresponding to the non-zero fixed-point of(2.31). Next, we construct the Hahn expansions of (1.65) as M (cid:0) Z − ( x ) (cid:1) = Z − ( x )1 + Z − ( x ) = 1 + √ x x + 1 + √ x = 1 − x + · · · (2.42) Y (cid:0) R − ( x ) (cid:1) = − R − ( x ) (cid:0) R − ( x ) − (cid:1) = x + x (2.43) M (cid:0) R − ( x ) (cid:1) = x + 1 x = x − + 1 (2.44) M (cid:48) (cid:0) R − ( x ) (cid:1) = 1 (cid:0) R − ( x ) + 1 (cid:1) = x x + 4 x + 1 = x − x + · · · , (2.45)whose leading powers are Ψ(0) = 0 , Γ(0) = 1 , Π(0) = − µ = 2 − Γ(0) = 1 (2.46) ν = − Π(0) = 1 (2.47) ρ = − Θ(0) = − η = − Ψ(0)( µ − ν − ρ ) = 0 , (2.49)once again, predicting η = 0, this time thanks both to the Bounded activities (0 < x i < K ) andto the Saturating nature of the interaction function M ( x ) (Sec. 1.4.1).16 .4 Biochemical dynamics As a Biochemical model we consider protein-protein interactions (PPI), which are driven bythree processes: φ → P i , describing the synthesis of the i th protein P i at a rate F ; P i → φ ,describing protein degradation at rate W ; P i + P j (cid:10) P i P j describing the binding and unbindingof a pair of interacting proteins at rates B and U respectively. The hetero-dimer P i P j undergoesdegradation P i P j → φ at rate Q . The dynamical equations for this system are d x i d t = F − W x i ( t ) + N (cid:88) j =1 U x ij ( t ) − B N (cid:88) j =1 A ij x i ( t ) x j ( t ) (2.50)d x ij d t = BA ij x i ( t ) x j ( t ) − ( U + Q ) x ij ( t ) , (2.51)where x i ( t ) is the concentration of P i and x ij ( t ) is the concentration of the hetero-dimer P i P j .Under time-scale separation we assume that the hetero-dimer concentration is at steady-state,setting d x ij / d t = 0 in (2.51). This provides us withd x i d t = F − W x i ( t ) − ˜ B N (cid:88) j =1 A ij x i ( t ) x j ( t ) , (2.52)where the effective binding rate is ˜ B = QB/ ( U + Q ). Hence, in the effective equation (2.52) wehave M ( x ) = F − W x , M ( x ) = − ˜ Bx and M ( x ) = x , providing the dynamic functions (1.64) R ( x ) = − M ( x ) M ( x ) = ˜ BxF − W x (2.53) Y ( x ) = M ( x ) R (cid:48) ( x ) = − F ˜ B x ( F − W x ) (2.54) Z ( x ) = R ( x ) M ( x ) = ˜ Bx F − W x . (2.55)The inverse functions are, therefore R − ( x ) = F x ˜ B + W x (2.56) Z − ( x ) = W B x − (cid:115) BFW x , (2.57)where in Z − ( x ) we choose the positive solution, corresponding to the positive fixed-point of(2.52). We can now compose the functions in (1.65), obtaining17 (cid:0) Z − ( x ) (cid:1) = Z − ( x ) = (cid:114) F ˜ B x + W B F x + · · · (2.58) Y (cid:0) R − ( x ) (cid:1) = − F ˜ B R − ( x ) (cid:0) F − W R − ( x ) (cid:1) = F ˜ B x + F ˜ BW x (2.59) M (cid:0) R − ( x ) (cid:1) = − ˜ BF x ˜ B + W x = − F x + F W ˜ B x + · · · (2.60) M (cid:48) (cid:0) R − ( x ) (cid:1) = 1 , (2.61)allowing us to extract the leading powers as Ψ(0) = 1 / , Γ(0) = 1 , Π(0) = 1 and Θ(0) = 0.These powers provide the dynamic exponents via (1.66), providing, for Biochemical µ = 2 − Γ(0) = 1 (2.62) ν = − Π(0) = − ρ = − Θ(0) = 0 (2.64) η = − Ψ(0)( µ − ν − ρ ) = − . (2.65) J ij ∈ E ( A ij , Ω) The stability of Eq. (1.1) can be characterized by the principal eigenvalue λ of J ij in (1.67) -(1.68). As the entries depend on k i and k j , the value of λ is driven by the limit of large k -specifically, in a degree-heterogeneous network, it is determined by the nearest neighbor degree κ nn , which can potentially diverge with the system size N as κ nn ∼ N β . (3.1)As opposed to the dynamic exponents Ω, the exponent β is a topological exponent, determinedby the network’s degree distribution P ( k ). In a scale-free network, for example β is extractedfrom (1.63).To approximate this hub-periphery structure, we use a star network, in which a single hub islinked to k = κ nn ∼ N β (3.2)peripheral nodes. This captures the environment of a typical hub in, e.g ., a scale-free network.Indeed, the hubs in a scale-free network can be viewed as a collection of weakly coupled stars ,that are only sparsely linked to each other . Under these conditions we have the ( k +1) × ( k +1)18etwork A ij = . . .
11 0 0 . . . . . . , (3.3)which using (1.67) and (1.68), provides the corresponding Jacobian as J ij = − C ˜ κ η nn k µ . . .
00 1 . . . . . . + k ν . . . k ν k ρ . . . k ρ . . . , (3.4)where C >
0. In (3.4) we used ˜ κ nn to express the nearest neighbor degree in our star construction,which is potentially distinct from κ nn of the originally approximated scale-free network. Lackingan a priori estimate for this parameter we express is as˜ κ nn ∼ N α , (3.5)leaving us a degree of freedom to later tune α such that the prediction from our star constructionbest captures the observed results from a real scale-free network. Hence, β , characterizing thestar-hub ( k ), is extracted from A ij via (3.1) and α is a tunable parameter, which we select forthe star-model to best fit the complete scale-free network results.We emphasize that the star approximation in (3.3) by no means captures the complete behaviorof a real scale-free network, as indeed it represents but a crude representation of an isolatedhub environment. However, our goal here is to examine stability vs. instability - a feature thatonly depends on the sign of the principal eigenvalue, not on its specific value, and is, therefore,insensitive to the detailed structure of A ij . As we show in Fig. 4 of the main text, the starapproximation, while highly stylized, is, indeed, sufficient to capture this J ij characteristic.To obtain the principal eigenvalue we solve the linear equation J v = λ v . (3.6)Using the symmetry of (3.4), we seek a solution of the form v = ( a, b, b, . . . , b ) (cid:62) , allowing us toreduce (3.6) into − C ˜ κ η nn k µ a + kk ν b = λa (3.7) k ρ a − C ˜ κ η nn b = λb. (3.8)Note that in v the specific values of a, b have no significance, only the ratio a/b , as we are onlyinterested in the direction of the eigenvector, not its magnitude. We therefore arbitrarily set a = 1, allowing us to solve (3.7) - (3.8) and obtain19 = 12 (cid:18) − C ˜ κ η nn ( k µ + 1) + (cid:113) C ˜ κ η nn ( k µ + 1) − C ˜ κ η nn k µ + 4 k ν + ρ (cid:19) , (3.9)where, of the two solutions, we selected only the one in which the square-root is added (ratherthan subtracted), as we seek the largest eigenvalue. Seeking the limitlim N →∞ (cid:0) λ (cid:1) , (3.10)we use (3.2) and (3.5) to rewrite (3.9) as λ ∼ (cid:16) − CN αη (cid:0) N βµ + 1 (cid:1) + (cid:113) C N αη ( N βµ + 1) + 4 (cid:2) − C N αη + βµ + N β (1+ ν + ρ ) (cid:3)(cid:19) , (3.11)in which we replace ˜ κ nn and k by N α and N β , respectively. In (3.11) we ignore the pre-factorsof the N -scaling, focusing only on the powers ( α, β ), hence substituting the equality sign = withthe asymptotic scaling operator ∼ . The case where µ >
0. First we analyze λ under µ >
0. Here, we write N βµ (cid:29)
1, simplifying(3.11) in the limit N → ∞ into λ ∼ (cid:18) − CN αη + βµ + (cid:113) C N αη + βµ ) + 4 (cid:2) − C N αη + βµ + N β (1+ ν + ρ ) (cid:3)(cid:19) . (3.12)Extracting the common terms out of the product we rewrite (3.12) as λ ∼ CN αη + βµ − (cid:115) (cid:18) − N − βµ + 1 C N σ (cid:19) , (3.13)where σ = β (cid:18) ν + ρ − µ − αβ η (cid:19) . (3.14)Note that N − βµ ≤ µ >
0. Therefore, in case σ ≥
0, the N σ term dominates ther.h.s. of the equation, providing, under N → ∞ , λ ∼ C √ C N σ + αη + βµ , (3.15)which, since C > i.e . the system is unstable .Next we consider the case σ <
0. Under these conditions N σ (cid:28)
1, allowing us to expand thesquare-root in (3.13) to first order, providing λ ∼ CN αη + βµ (cid:32) − (cid:18) − N − βµ + 1 C N σ (cid:19)(cid:33) = − CN αη + 1 C N σ + αη + βµ . (3.16)20e can rewrite this in the form of Eq. (9) of the main text, obtaining λ ∼ N Q (cid:18) − CN S (cid:19) , (3.17)where Q = σ + αη + βµ (3.18) S = σ + βµ. (3.19)As N Q is guaranteed to be positive, the sign of λ depends on S : in case S > CN − S (cid:28)
1, and hence we have λ ∼ N Q >
0, an unstable dynamics. Ifhowever
S <
0, we have − CN − S → −∞ , predicting λ <
0, regardless of C , an asymptoticallystable system. Consequently the system is stable as long as S <
0, which using (3.19), andtaking σ from (3.14), predicts the stability condition β (cid:18) ν + ρ − µ − αβ η (cid:19) < . (3.20)This condition reduces stability into a small set of relevant parameters: β , characterizing thenetwork topology A ij , and Ω = ( η, µ, ν, ρ ), associated with the dynamics M ( x ) , M ( x ) , M ( x ).The non-universal parameter C becomes irrelevant in the limit of sufficiently large N . Finally, α can be selected to tune the star approximation optimally to the case of a real scale-free A ij , as we do below. Note that condition (3.20) can also be expressed as σ + βµ <
0. Thisimplies that if σ > βµ the system is unstable. This instability condition already contains thepreviously obtained condition of σ >
0, that led to Eq. (3.15). Therefore, Eq. (3.20) is sufficientto characterize the system’s stability.Finally, the case S = 0 in (3.17) represents sensitive stability, in which λ ’s value is not asymp-totically defined, and rather it depends on the coefficient C . In this class, stability in no longera characteristic of the dynamic model , but rather of its specific rate constants, as encapsulatedwithin C . A trivial example is when β = 0, i.e . a non fat-tailed P ( k ), for example - Erd˝os-R´enyi.Indeed, as we discuss in the main text, in such networks, stability can be tuned by the modelparameters, lacking a defined asymptotic behavior. The case where µ ≤
0. Here, we have N βµ ≤
1, and hence Eq. (3.11) can be approximated by λ ∼ (cid:18) − CN αη + (cid:113) C N αη + 4 (cid:2) − C N αη + βµ + N β (1+ ν + ρ ) (cid:3)(cid:19) , (3.21)leading to λ ∼ CN αη − (cid:115) (cid:18) − N βµ + 1 C N ω (cid:19) , (3.22)where ω = β (cid:18) ν + ρ − αβ η (cid:19) . (3.23)21nce again, we have N βµ ≤
1, this time since µ <
0. Therefore, as before, if ω > λ becomesdominated by the N ω term, following λ ∼ C √ C N ω + αη , (3.24)which is always positive, i.e. unstable . If, however, ω <
0, we use a linear approximation towrite λ ∼ CN αη (cid:32) − (cid:18) − N βµ + 1 C N ω (cid:19)(cid:33) = − CN αη + βµ + 1 C N ω + αη , (3.25)once again - a competition between a positive vs. a negative term. Collecting the powers we,again, rewrite (3.25) in the form λ ∼ N Q (cid:18) − CN S (cid:19) , (3.26)where this time Q = ω + αη (3.27) S = ω − βµ. (3.28)Stability is ensured if, for N → ∞ , the negative term dominates, namely if S <
0. Using (3.23)to express ω this provides β (cid:18) ν + ρ − µ − αβ η (cid:19) < , (3.29)recovering precisely the condition in (3.20). Tuning α . The parameter α represents a degree of freedom, rooted in the fact that the starapproximation captures an inaccurate depiction of a real scale-free network. This approximationcan be optimized if we select α such that the star best captures the complete scale-free networkenvironment. We find that setting α = β S = β (1 + ν + ρ − µ − η ) , (3.31)as appears in Eq. (10) of the main text. The stability condition (3.31) is driven by five distinct exponents. The first four exponentsΩ = ( η, µ, ν, ρ ) are determined by the dynamic model - Social, Regulatory, Ecological etc. -22
𝑺 𝑺𝑺 (a) (b)(c) (d) 𝜶𝜷 = 𝟎
Error :
𝟏𝟏. 𝟓% 𝜶𝜷 = 𝟏𝟒
Error : 𝜷 = 𝟏 Error :
𝟒. 𝟑% 𝜶 𝜷 = 𝟏 Error :
𝟓. 𝟓%
Figure 7:
Tuning the parameter α . In Eqs. (3.20) and (3.29) we have the degree of freedom toset α/β for the star approximation to best predict the stability of complete scale-free networks.We used our set of 2 , J matrices (Fig. 4a of main text) to predict stability using differentvalues of α/β . Quantifying the Error by the fraction of mis-classified J matrices (grey dots),we find that setting α/β = 1 / ∼
96% correct classification(Error= 4 . A ij or of microscopic model parameters, encapsulated within C , and are therefore hardwired into the system’s dynamic behavior. For example, our formalism predicts that the SISmodel (Social) has Ω = (0 , , − ,
0) (Sec. 2.1). This prediction is characteristic of the SIS model,namely it is an intrinsic feature of the SIS interaction mechanisms of infection and recovery. Itis, therefore, insensitive to the specific rates of the model - hence Ω remains unchanged if, e.g. ,a disease has a high or low infection rate, or if it is transmitted via physical contact or aerosols.These will impact the constant C in (1.67), but will have no impact on the universal scaling.Similarly, Ω is unaffected by A ij . Therefore, regardless of whether the disease spreads along thestandard social network (Flu) or via sexual transmission (AIDS), as long as the mechanism isSIS (or any other mechanism for that matter) Ω remains the same.The remaining exponent in (3.31), β , on the other hand, is independent of the dynamics, anddetermined solely by A ij , specifically by its degree distribution P ( k ), capturing the divergenceof κ nn in the limit of large N . For example, under a scale-free topology, where P ( k ) ∼ k − γ , wehave κ nn following Eq. (1.63), predicting β = γ < − γ ≤ γ < γ ≥ . (3.32)Therefore under γ ≥
3, we have β = 0 and stability becomes sensitive to parameters ( C ). On theother hand, under a scale-free topology ( γ <
3) we observe the asymptotic stability/instabilitydriven by S (cid:54) = 0. S under negative interactions Our formalism is designed to treat cooperative interactions, in which the interaction term A ij M ( x i ) M ( x j ) >
0. Strictly speaking, our analysis is not designed to treat adversarial inter-actions, in which the mean-field derivation of Sec. 1 may fail to capture the complex patternsof all activities x i ( t ). Specifically, in the presence of negative interactions the steady-states x i cannot always be expressed as a single variable function of k i , as appears in Eq. (1.18). Rather,it also depends on i ’s neighborhood, e.g ., if the neighborhood has high activities, the adversarialinteractions will drive x i to be small, and vice versa. In certain cases, such as in our Biochemicalmodel, this effect is marginal, all neighborhoods have (roughly) the same activity, and henceour analysis continues to apply. The star approximation, however, cannot distinguish betweenpositive and negative interactions - as the derivation yields an identical outcome in both cases.Indeed, this approximation overlooks the potentially stabilizing effect of the negative off-diagonalterms.Therefore, for simplicity, to evaluate S under negative interactions we consider a fully connectednetwork of k ∼ N β (3.33)nodes, `a la May’s original formulation. This simplification may not allow us to fully evaluatethe contribution of A ij , as, indeed, we have replaced the network by a single large clique,however, we are only interested in the stability of the dynamics , which remains independent of24his simplifying assumption. With all nodes now having degree k , the resulting Jacobian takesthe form J ij = − Ck η k µ . . . k µ . . . . . . k µ − k ν + ρ . . . k ν + ρ k ν + ρ . . . k ν + ρ ... . . . ... k ν + ρ k ν + ρ . . . , (3.34)in which both the diagonal and the off-diagonal terms are preceded by a minus sign.For a matrix of the form (3.34) the principal eigenvector is v = (1 , − , , . . . , (cid:62) , whose associ-ated eigenvalue is λ = − Ck η + µ + k ν + ρ . (3.35)Using (3.33) we rewrite (3.35) in the form λ = N Q Neg (cid:18) − CN S Neg (cid:19) , (3.36)where now Q Neg = β ( ν + ρ ) (3.37) S Neg = β ( ν + ρ − µ − η ) . (3.38)Hence under negative dynamics the stability classifier is different from (3.29), being S Neg = S − β .For example, in Biochemical, where Ω = ( − , , − ,
0) we have S Neg = − β <
0, predicting astable dynamics.
To numerically test our predictions we constructed Eq. (1.1) for each of the systems in Sec. 2,using the appropriate A ij (Scale-free, Erd˝os-R´enyi, empirical, etc.). We then used a fourth-orderRunge-Kutta stepper (Matlab’s ode45) to numerically solve the resulting equations. Startingfrom an arbitrary initial condition x i ( t = 0), i = 1 , . . . , N we allowed the system to reach steady-state by waiting for ˙ x i →
0. To numerically realize this limit we implemented the terminationcondition N max i =1 (cid:12)(cid:12)(cid:12)(cid:12) x i ( t n ) − x i ( t n − ) x i ( t n )∆ t n (cid:12)(cid:12)(cid:12)(cid:12) < ε, (4.1)25here t n is the time stamp of the n th Runge-Kutta step and ∆ t n = t n − t n − . As the systemapproaches a steady-state, the activities x i ( t n ) become almost independent of time, and thenumerical derivative ˙ x i = x i ( t n ) − x i ( t n − ) / ∆ t n becomes small compared to x i ( t n ). The con-dition (4.1) guarantees that the maximum of ˙ x i /x i over all activities x i ( t n ) is smaller than thepre-defined termination variable ε . In our simulations, across the different dynamics we tested,we set ε ≤ − , a rather strict condition, to ensure that our system is sufficiently close to the true steady-state. J ij Once the steady state x = ( x , . . . , x N ) (cid:62) is reached we construct the numerical Jacobian bysubstituting the numerically obtained states x i into (1.35) and (1.46). This represents thesystem’s actual stability matrix, incorporating its specific dynamics on the relevant network A ij . For example, consider our Social model in Eq. (2.3), where M ( x ) = − f x, M ( x ) = 1 − x and M ( x ) = gx , and hence M (cid:48) ( x ) = − f, M (cid:48) = − M (cid:48) = g . Once we obtain x we introduceall numerically calculated x i into D ii and W ij , which for Social take the form D ii = − f + − N (cid:88) j =1 A ij gx j (4.2)and W ij = A ij (1 − x i ) g. (4.3)In Fig. 2 of the main text we compare the scaling of these numerically estimated D ii and W ij vs. the theoretically predicted ensemble E ( A ij , Ω), shown in (1.67) and (1.68).
Our main theoretical prediction focuses on scaling relationships, such as D ( k ) ∼ k µ , which weobserve by their linear slope in a log-log plot. To construct such plots we employed logarithmicbinning . First we divide all nodes into B bins B ( b ) = (cid:8) i = 1 , . . . , N (cid:12)(cid:12) c b − < k i ≤ c b (cid:9) , (4.4)where b = 1 , ..., B and c is a constant. In (4.4) the b th bin includes all nodes i whose degrees k i are between c b − and c b . The parameter c is selected such that the unity of all bins ∪ Bb =1 B ( b )includes all nodes, hence we set c B = max Ni =1 k i . We then plot the average degree of the nodesin each bin k b = (cid:104) k i (cid:105) i ∈ B ( b ) = 1 | B ( b ) | (cid:88) i ∈ B ( b ) k i (4.5)versus the average D ( k ) term of nodes in that bin26 ( k b ) = (cid:104) D ii (cid:105) i ∈ B ( b ) = 1 | B ( b ) | (cid:88) i ∈ B ( b ) D ii . (4.6)In a similar fashion we plot W ij ∼ k νi k ρj , this time applying the binning on k νi k ρj , instead of k i .Therefore, the bins are defined as B ( b ) = (cid:8) ( i, j ) (cid:12)(cid:12) A ij = 1 , c b − < k νi k ρj ≤ b w (cid:9) , (4.7)such that in each bin we include all node pairs whose product k νi k ρj is within a given range. Asabove, we then plot the average W ij in each bin, providing the real Jacobian terms J Real b = (cid:104) W ij (cid:105) ( i,j ) ∈ B ( b ) = 1 | B ( b ) | (cid:88) ( i,j ) ∈ B ( b ) W ij (4.8)vs. the average value of k νi k ρj in that bin J Th b = (cid:104) k νi k ρj (cid:105) ( i,j ) ∈ B ( b ) = 1 | B ( b ) | (cid:88) ( i,j ) ∈ B ( b ) k νi k ρj . (4.9) To test our predictions we used model and real networks, as summarized below: ER . An Erd˝os-R´enyi random network with N = 6 ,
000 nodes and an average degree of (cid:104) k (cid:105) = 6. SF1,2,3 . A set of binary scale-free networks, constructed via the configuration model , with N = 6 ,
000 nodes, (cid:104) k (cid:105) = 6 and degree distribution following P ( k ) ∼ k − γ with γ = 2 , . UCIonline (Social). An instant messaging network from the University of California Irvine ,capturing 61 ,
040 transactions between 1 ,
893 users during a T = 218 day period. Connectingall individuals who exchanged messages throughout the period, we obtain a network of 1 , ,
670 links, exhibiting a fat-tailed degree distribution.
Email Epoch (Social). This dataset monitors ∼ × emails exchanged between 3 , ∼ , giving rise to a scale-free social network with31 ,
885 binary links.
PPI1 (Regulatory, Biochemical). The yeast scale-free protein-protein interaction network, con-sisting of 1 ,
647 nodes (proteins) and 5 ,
036 undirected links, representing chemical interactionsbetween proteins . PPI2 (Regulatory, Biochemical). The human protein-protein interaction network, a scale-freenetwork, consisting of N = 2 ,
035 nodes (protein) and L = 13 ,
806 protein-protein interactionlinks . ECO1,2 (Ecological). To construct the mutualistic ecological networks we collected data onsymbiotic interactions of plants and pollinators in Carlinville Illinois from . The resulting 456 × ,
429 network M ik is a bipartite graph linking the 456 plants with their 1 ,
429 pollinators. Whena pair of plants is visited by the same pollinator they mutually benefit each other indirectly, byincreasing the pollinator populations. Similarly pollinators sharing the same plants also posses27n indirect mutualistic interaction. Hence we can collapse M ik to construct two mutualisticnetworks: The 1 , × ,
429 pollinator network ECO1 and the 456 ×
456 plant network ECO2.The resulting networks are B kl = (cid:88) i =1 M ik M il (cid:80) ns =1 M is , (4.10)for the pollinator network (ECO1), and A ij = , (cid:88) k =1 M ik M jk (cid:80) ns =1 M sk , (4.11)for the plant network (ECO2). In both networks the numerator equals to the number of mutualplants ( B kl ) or pollinators ( A ij ). For each mutual plant i (pollinator k ) we divide by the overallnumber of plants (pollinators) that share i ( k ). Hence, the weight of the mutualistic interactionin, e.g. , A ij is determined by the density of mutual symbiotic relationships between all plants,where: (i) the more mutual pollinators k that plants i and j share the stronger the mutualisticinteraction between them; (ii) on the other hand the more plants pollinated by k the smaller isits contribution to each plant. A similar logic applies also for the pollinator network B ij . Thisprocess potentially allows us to have isolated components, e.g. , single disconnected nodes. Thestate of these isolated nodes is decoupled from the state of the rest of the network, and hencein our analysis we only focused on the giant connected component of A ij and B ij , comprisingall 456 plants, rendering A ij to be a fully connected component, but only 1 ,
044 pollinators,eliminating 385 isolated pollinators in B ij . Our analytical derivation revolves around cooperative interactions, in which A ij ∈ { , } , i.e .non-negative, and the functions M ( x ) , M ( x ) in (1.1) are also positive for all x . The onlyexception is our Biochemical model, in which the hetero-dimer interaction ( − A ij x i x j ) capturesdepletion of i and j , and hence all interactions are negative. Some applications, however, involvea mixture of cooperative and adversarial interactions, requiring us to analyze signed networksin which A ij ∈ {− , , } . A classic example is in population dynamics, which, in many cases,represent a mixture of mutualistic interactions, coexisting alongside competitive and predatorydynamics. Strictly speaking, under these conditions, the Jacobian J ij cannot be be extractedfrom the E ( A ij , Ω) ensemble. Yet, as we next show, even then, our theory can still help usunderstand the emergence of stability.Consider the Ecological modeld x i d t = Bx i ( t ) (cid:18) − x i ( t ) K (cid:19) + N (cid:88) j =1 A ij x i ( t ) x j x j , (5.1)28 b) (c) (d) (g) (h) (i) (a) (e) (f) Figure 8:
Emerging stability under mixed interactions . (a) Starting with a scale-free A ij with 10% negative links, the system reaches a fixed-point in which nodes with many negativelinks become extinct. As a result the surviving 5 ,
981 node-network is slightly biased towardspositive links (9 .
8% vs. the initial 10% negative). (b) The Jacobian matrix continues to followthe analytically predicted scaling J ii ∼ k µi , selectively for the (majority) of positive k i . (c)Similarly the off-diagonal terms also follow the predicted J ij ∼ k νi k ρj (for k i , k j >
0) . (d)-(f)Similar analysis for an initial 25% negative links. (g)-(i) Even under 50% negative links thesurviving network Jacobian (4 ,
635 nodes) can still be approximated by E ( A ij , Ω).29n which the the existing links in A ij include a random fraction of 0 ≤ f ≤ A ij = −
1) and a remaining 1 − f fraction of links that are positive ( A ij = 1). As the systemevolves in time, the activities x i ( t ) are positively impacted by the cooperative interactions, butsuppressed by the adversarial ones. Therefore, nodes with many negatively interacting partnerswill have their activities driven towards x i ( t ) →
0. The crucial point is that once a node reacheszero activity, it becomes extinct, and effectively removed from the network, together with its setof links. As a consequence the system will reach a steady state, in which the surviving nodes arebiased towards positive k i . This represents a quite general, and most importantly, a naturallyemerging organizing principle that leads real world A ij to have a majority positive set of links .To exemplify this, we constructed a scale free A ij with N = 6 ,
000 nodes and L = 18 ,
000 links,of which a fraction f = 0 . x i ( t ) ≤ (cid:15) we set it permanently tozero (extinction), effectively removing it from the network; we used (cid:15) = 10 − . At its finalstate the network is reduced to N = 5 ,
981 surviving nodes, featuring a (slight) decrease in thefraction of negative links, from 10% to 9 .
8% (Fig. 8a). This is a minor effect, of course, however,as we increase the negative link fraction, f , this positive-link bias becomes more pronounces.Indeed, under an initial condition where f = 0 .
25, the system reached a final state in whichthe negative links are reduced to 22% (Fig. 8d), and, finally, starting with a balanced positive-negative network ( f = 0 . J matrices include bothpositive and negative weights, and hence, strictly speaking J / ∈ E ( A ij , Ω). However, thanks tothe pruning effect, the actual J matrices remain predominantly positive, and the majority ofsurviving nodes continue to have k i >
0. Hence, while imperfect, we still observe that J follows(approximately) the scaling patterns of the E ( A ij , Ω) ensemble. Indeed, Fig. 8b,c,e,f,h,i confirmsthat, despite the inclusion of negative links, J ii and J ij continue to follow (1.67) and (1.68) with µ = 1 , ν = 1 and ρ = −
2, as predicted for Ecological.Taken together, this exemplifies that both under fully positive and under mixed interactions,real-world Jacobian matrices are, indeed, non-random, exhibiting (exactly or approximately) theunique interplay between topology and dynamics expressed in E ( A ij , Ω). These patterns requireno devious strategies , but rather arise from naturally emerging organizing principles: for positivedynamics - the scaling Ω, and for mixed interactions that, in combination with the pruning ofnegative links via extinction. 30 eferences [1] B. Barzel and A.-L. Barab´asi. Universality in network dynamics.
Nature Physics , 9:673 –681, 2013.[2] M.E.J. Newman.
Networks - an introduction . Oxford University Press, New York, 2010.[3] J. Gao, B. Barzel and A.-L. Barab´asi. Universal resilience patterns in complex networks.
Nature , 530:307–312, 2016.[4] G. Caldarelli.
Scale-free networks: complex webs in nature and technology . Oxfrod Univer-sity Press, New York, 2007.[5] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. Hwang. Complex networks:Structure and dynamics.
Physics Reports , 424:175–308, 2006.[6] L. Schmetterer and K. Sigmund (Eds.).
Hans Hahn Gesammelte Abhandlungen Band1/Hans Hahn Collected Works Volume 1 . Springer, Vienna, Austria, 1995.[7] R.M. May. Will a large complex system be stable?
Nature , 238:413 – 414, 1972.[8] P.S. Dodds and D.J. Watts. A generalized model of social and biological contagion.
Journalof Theoretical Biology , 232:587–604, 2005.[9] G. Karlebach and R. Shamir. Modelling and analysis of gene regulatory networks.
NatureReviews , 9:770–780, 2008.[10] B. Barzel and O. Biham. Binomial moment equations for stochastic reaction systems.
Phys.Rev. Lett. , 106:150602–5, 2011.[11] R.M. May. Simple mathematical models with very complicated dynamics.
Nature , 261:459–467, 1976.[12] C.S. Holling. Some characteristics of simple types of predation and parasitism.
The Cana-dian Entomologist , 91:385–398, 1959.[13] E.O. Voit.
Computational Analysis of Biochemical Systems . Cambridge University Press,New York, NY, 2000.[14] H.I. Schreier, Y. Soen and N. Brenner. Exploratory adaptation in large random networks.
Nature Communications , 8:14826, 2017.[15] S. Milojevi´c. Power-law distributions in information science: making the case for logarithmicbinning.
Journal of the American Society for Information Science and Technology , 61:2417–2425, 2010.[16] T. Opsahl and P. Panzarasa. Clustering in weighted networks.
Social Networks , 31:155–163,2009.[17] J.-P. Eckmann, E. Moses and D. Sergi. Entropy of dialogues creates coherent structures ine-mail traffic.
Proc. Natl. Acad. Sci. USA , 101:14333–7, 2004.[18] H. Yu et al.
High-quality binary protein interaction map of the yeast interactome network.
Science , 322:104–110, 2008. 3119] J.F. Rual et al . Towards a proteome-scale map of the human protein-protein interactionnetwork.