Dynamic stability of complex networks


 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 theory and practice. While empirical reality answers with an astounding yes, the mathematical analysis, based on linear stability theory, seems to suggest the contrary - hence, the diversity-stability paradox. Here we settle this dichotomy, by considering the interplay between topology and dynamics. We show that this interplay leads to the emergence of non-random patterns in the system's stability matrix, leading us to relinquish the prevailing random matrix-based paradigm. Instead, we offer a new matrix ensemble, which captures the dynamic stability of real-world systems. This ensemble helps us analytically identify the relevant control parameters that predict a system's stability, exposing three broad dynamic classes: In the asymptotically 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 real systems reside, in which diversity not only does not prohibit, but, in fact, enhances dynamic stability. Finally, in the sensitively stable class diversity plays no role, and hence stability is driven by the system's microscopic parameters. Together, our theory uncovers the naturally emerging rules of complex system stability, helping us reconcile the paradox that has eluded us for decades.

The study of complex social, biological and technological systems, is often directed towards dramatic events, such as cascading failures [1][2][3][4][5] or abrupt state transitions [6][7][8][9][10] .In reality, however, these represent the exception rather than the rule.In fact, the truly intriguing mathematical puzzle is how, despite enduring constant perturbations and local obstructions, many systems continue to sustain reliably stable dynamics [11][12][13][14][15] .This challenge is rooted in the diversity-stability paradox 16 , which predicts that a sufficiently large complex system, will inevitably become unstable.What, then, are the naturally emerging organizing principles, or in May's original wording -what are nature's devious strategies, to achieve this ubiquitously observed dynamic stability?Initial clues began to unveil with the mapping of empirical networks, which uncovered universally recurring topological characteristics [17][18][19] that can potentially impact the system's dynamics 20 .For example, degree heterogeneity 21 , community structure 22 and topological symmetries 23,24 The diversity-stability challenge Consider a complex system of N interacting components (nodes), whose dynamic activities x(t) = (x 1 (t), . . ., x N (t)) are driven by nonlinear pairwise interactions.The system's fixedpoints capture static states, which, unperturbed, remain independent of time.The stability of these fixed-point can be examined through their response to small perturbations δx, which, in the linear regime, can be approximated by where J represents the system's Jacobian matrix.We express J as 16 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, σ 2 ).Using I to represent the identity matrix, J's diagonal terms J ii = −C, capture the typical relaxation 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 denote this ensemble by E(P (w), C).
The system E(P (w), C) is dynamically stable if J's principle eigenvalue satisfies Re(λ) < 0. The challenge is that, according to random matrix theory, in the limit of large N we have 16 Re(λ) which is positive, and hence unstable, as long as C < √ N .Consequently, a sufficiently large system (N → ∞) becomes inevitably unstable, i.e.May's paradox.
The role of the network topology.A natural attempt to reconcile this paradox is by considering structural patterns in the interaction topology A ij , a binary network (A ij = 0, 1), which is typically sparse and, most often, highly heterogeneous 17 .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 matrix multiplication 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 diagonal weights from P (w) and C (Fig. 1a,b).
In the E(A ij , P (w), C) ensemble, the principal eigenvalue λ depends on the average nearest neighbor degree 32 κ nn , whose magnitude depends on the system's degree heterogeneity 22 .For instance, in a random network with degree distribution P (k) we have 7,32 κ nn = k 2 / k , in which the second moment k 2 increase with P (k)'s variance, and hence with A ij 's heterogeneity.Most generally, for a sufficiently broad, i.e. fat-tailed, P (k), we have 22 an asymptotic divergence with system size.As an example, consider a scale-free network, with , for which κ nn follows (4) with β = 3 − γ.Under Eq. ( 4) the principal eigenvalue in the E(A ij , P (w), C) ensemble follows 32 (Fig. 1c) in which the square-root of ( 3) is replaced by the exponent β/2.The crucial point is that, once again, we observe an asymptotic behavior that prohibits stability in the limit N → ∞.
These mathematical observations, now nearing their 50th anniversary, seem to suggest that complex networked systems must be limited in size in order to maintain stability.Reality, however, consistently confronts us with the contrary 33 .To reconcile this incongruity, we next show that under real nonlinear dynamics the structure of J is fundamentally different from the currently considered Jacobian ensembles.This allows us to fully predict a complex system's dynamic stability, 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 -indeed the paradigm -is that one can separate the role of the topology from that of the dynamics.For example, to construct J in the E(A ij , P (w), C) ensemble, one first selects the topology A ij , and then independently extracts the dynamic response terms, W ij , from P (w), and the relaxation rate C.Such construction overlooks the complex interplay between A ij and P (w), C, rooted in the nonlinear mechanisms driving the interactions between the nodes.Omitting this interplay 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 or technological.This, of course, stands in sharp contrast with the frequently observed fact that different nonlinear mechanisms, despite being run on similar networks, may express fundamentally distinct dynamic behaviors [35][36][37] .
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 each system's intrinsic dynamic mechanisms.These mechanisms are inherent to the nature of the system's interacting components: in social systems, for example, individuals interact through infection and recovery [38][39][40] , in biological networks, proteins, genes and metabolites are linked through biochemical processes [41][42][43][44] , and in ecological systems, species undergo competitive or symbiotic interactions [45][46][47] .All these processes capture built-in dynamics, representing the physical mechanisms that drive all nodes/links.Most generally, these dynamic mechanisms can be described by 35 where M 0 (x) captures the self-dynamics of all nodes, and M 1 (x), M 2 (x) describe their pairwise interactions.With the appropriate selection of these nonlinear functions Eq. ( 6) captures a broad range of frequently used models in social 38,39 , biological [41][42][43][44][45] and technological 48 systems (Fig. 1g).
The dynamic Jacobian ensemble (Fig. 1e).Extracting J from ( 6), we show in Supplementary Section 1, leads to a currently unexplored matrix ensemble, in which the weights J ij are strongly intertwined with the topology A ij , via In ( 7) and ( 8), κ 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 0 (x), M 1 (x), M 2 (x), as shown in Box I.The coefficient C, similar to the Jacobian in (2), represents the system specific rate parameters and characteristic time-scales, a parameter that can potentially be affected by external perturbation or changing environmental conditions 7 .
The result is a non-random Jacobian structure, whose entries, as opposed to the existing ensembles, cannot be selected independently from an arbitrary P (w), or set to a constant value C. Instead (7) and ( 8) define a new ensemble E(A ij , Ω), which matches, for each underlying network A ij , a predictable set of weights J ii and J ij based on the dynamic exponents Ω.Therefore, in contrast to the random matrix based treatment of E(P (w), C) or E(A ij , P (w), C), that has led the discussion since May's original analysis 16 , the newly introduced E(A ij , Ω) captures the effect of the system-specific nonlinear dynamics.Indeed, in (7) and (8), identical networks may 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 relevant empirical networks.This includes four dynamic models: Social.The susceptible-infectedsusceptible (SIS) model [38][39][40] for disease spreading; Regulatory.The Michaelis-Menten model 41 for gene regulation; Ecological.Mutualistic interactions 45 in ecology; Biochemical.Proteinprotein interactions [42][43][44] in sub-cellular networks.Applying each model to five different model and real networks, we arrive at a total of 20 combinations of networks/dynamics as a testing ground for our predicted J-ensemble (for a detailed description of all dynamic models and networks see Supplementary Sections 2 and 4.4).
Perturbing the system around its numerically obtained fixed-points, we constructed the Jacobian matrix 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 (8).For example, in our Social model we predict µ = 1, while our Regulatory dynamics are predicted to have µ = 0, both scaling relationships clearly evident in Fig. 2a and c.This means that setting all diagonal terms uniformly to −C, as in May's original formulation, misses the actual patterns that arise from the nonlinear Social/Regulatory dynamics.Similarly, the off-diagonal terms are proportional to k −1 i k 0 j in Social (Fig. 2b) and k 0 i k −2 j in Regulatory (Fig. 2d), again overruling the classic construction in which the J ij entries are extracted from an arbitrary P (w).These scaling relationships, encapsulated within our analytically predicted Ω, are independent of A ij , as indeed Fig. 2 confirms, showing recurring scaling patterns across our model (SF1, SF2, SF3) and relevant empirical (Epoch, UCIonline, PPI1 etc.) networks.Hence, Ω captures the intrinsic, and most crucially, 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 empirical social, biological and ecological networks.For each network we construct J both from the random ensemble 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 Jacobian matrices, 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 (8).These differences directly impact the network's dynamic stability.Indeed, in all cases we find that the random J has a positive λ (red), congruent with May's analysis, while the dynamic J, extracted from E(A ij , Ω) has λ < 0 (blue).Hence, these empirical networks, when coupled with relevant real-world dynamics, are, in fact stable.
Together, our analysis demonstrates that: (i) actual J matrices, extracted from real nonlinear dynamics 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 intrinsic nonlinear mechanisms (M 0 (x), M 1 (x), M 2 (x)) through Eqs. ( 7) and (8).Most crucially, while the 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 much richer stability profile, in which a broad family of empirically relevant dynamics are predictably stable.

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. ( 4).This setting allows us to observe the combined effect of dynamic nonlinearity (Ω) with one of the most ubiquitous features of real world networks, i.e. degreeheterogeneity 17 .In Supplementary Section 3 we show that extracting J from E(A ij , Ω), with a degree-heterogeneous A ij , the principal eigenvalue asymptotically follows where Q and S depend both on network structure, through β in (4), and on the nonlinear dynamics via Ω.Specifically, we show that replacing the fixed exponents, 1/2 in (3), under E(P (w), C), and β/2 in (5), under E(A ij , P (w), C), with a variable exponent that depends on the dynamics-driven Ω = (η, µ, ν, ρ).
Equations ( 9) and ( 10) represent our key result.They show that, when extracted from real nonlinear dynamics of the form (6), the asymptotic behavior of λ is different from the one predicted in May's original formulation, i.e.Eqs. ( 3) and ( 5).Most importantly, ( 9) and ( 10) have crucial implications regarding the system's dynamic stability, distinguishing between three potential stability classes (Fig. 1f): • Asymptotic instability: S > 0 (red).In case S in ( 9) is positive, we have, for sufficiently large N , Re(λ) ∼ N Q > 0. Such systems commit to the classic diversity-stability prediction, and, when large enough, become intrinsically unstable, regardless of the model parameters (C).Specifically, setting η = µ = ν = ρ = 0 we revert to the random matrix ensemble (E(A ij , P (w), C)), for which S = β > 0, thus recovering May's predicted instability.Hence, the diversity stability paradox is, indeed, a valid prediction, yet it represents a specific point in a broader class of potential dynamic behaviors, driven by the parameter S in (10).
• Asymptotic stability: S < 0 (blue).For a broad range of empirically relevant dynamics we have S < 0. Under these conditions, independently of C, the r.h.s. of ( 9) is dominated by the negative term, guaranteeing that Re(λ) < 0. Therefore, in this class, stability is not simply enabled, but rather, under N → ∞, it becomes asymptotically inevitable, even if C ≪ 1. Recalling May's original question: can a large complex system be stable?-Eq.( 9) with S < 0 provides a clear answer: it not only 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 > 1, i.e. sufficiently strong negative feedback, the system is stable, otherwise it becomes unstable.As opposed to Ω, which is intrinsic to the dynamic model, depending on the functional form of M 0 (x), M 1 (x), M 2 (x) in ( 6), C is determined by the specific model parameters.Hence, in this class stability is not an intrinsic characteristic of the system, but rather it is sensitive to Eq. ( 6)'s tunable parameters.
This classification settles the debate on complex system stability on several levels.In the context of the diversity-stability paradox, it shows that large complex systems can, and in some cases even must be stable -hence reconciling between theory and empirical observation.More broadly, the stability classifier S, identifies the relevant topological (β) and dynamic (Ω) control parameters that help analytically predict the stability class of any system within the form (6). We can therefore use S to predict a priori whether a specific combination of topology and dynamics will exhibit stable functionality or not.
To examine our stability classifier S we used our model and real networks to extract 2, 077 Jacobian matrices from the E(A ij , Ω) ensemble, with different sets of η, µ, ν and ρ.In Fig. 4a we 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 asymptotically unstable class (red, top-right) has S > 0 and consequently also λ > 0, a guaranteed instability.The asymptotically stable class (blue, bottom-left) is observed for S < 0, and has, in all cases λ < 0, 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 (∼ 4%) of our sampled J matrices were inaccurately classified by S (grey), a consequence of the approximate nature of our derivations (Supplementary Section 3).
In Fig. 4a we also present the stability matrices associated with our empirical social, biological and ecological networks (symbols) examined earlier in Figs. 2 and 3.For each of these networks we 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 ( 8), and Ω taken according to the relevant dynamic model (Fig. 1g).While all these networks are classified unstable 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 stable class (blue symbols).
Figures 2 -4, together demonstrate our complete theoretical path: First, Figs. 2 and 3 show that real J matrices exhibit the internal scaling patterns predicted in ( 7) and ( 8) -representing a steep departure from the broadly applied random matrix paradigm.Next, Fig. 4 shows that the resulting Jacobian ensemble has a rich space of potential dynamic behaviors, significantly more diverse than the currently considered ensembles.Most importantly -this space contains a broad class of asymptotically stable (blue) or sensitively stable (green) systems, that do not fall within May's paradox.
Taken together, our analysis demonstrates that nature must not rely on devious strategies in order to ensure dynamic stability.The observed stability of large complex systems emerges quite naturally thanks to the built-in nonlinear interaction mechanisms between the system's components.

The ingredients of dynamic stability
The parameter S in (10) 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 dynamic behavior.To understand this, consider, for example, our Social model, for which our formalism predicts Ω = (0, 1, −1, 0).This prediction is rooted in the SIS interaction mechanisms, i.e. infection vs. recovery, expressed through the functional form of M 0 (x), M 1 (x), M 2 (x) of (6).It is, therefore, insensitive to the specific parameters of the model, e.g., if the disease has a high or 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 Ω, and therefore 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 (10), β, is independent of the dynamics, determined solely by A ij , specifically by its degree distribution P (k), through (4).This parameter quantifies the fat-tailed nature of P (k), being β = 0 for a bounded distribution, and approaching unity under extreme 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 or sensitive.
The role of degree-heterogeneity.The prevalence of fat-tailed P (k) is among the defining features of real-world complex systems, from biological 49 to social 50,51 and technological 52,53 networks, with far-reaching implications on their observed dynamic behavior 7,32,35 .Our analysis indicates that this network characteristic may also play a crucial role in the context of dynamic stability.To understand this, consider a non fat-tailed P (k), such as an Erdős-Rényi, network, where the degrees follow a Poisson distribution, having β = 0.Under these conditions we have S = 0 in (10), the system has no defined asymptotic behavior, and hence it is sensitively stable -i.e. its stability depends on model parameters.Therefore, the existence of asymptotic stability/instability is a direct consequence of degree-heterogeneity, as, indeed, these forms of dynamic stability rely on β > 0.
This uncovers a new additional layer to the dynamic impact of P (k) on complex system functionality.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 7 .Such forces impact a system's functionality by perturbing its dynamic parameters, namely they affect the coefficient C. Seldom, however, do these environmental perturbations affect the system's built-in interaction mechanisms.Indeed, while the dynamic mechanisms are fixed, ingrained in the physics of the interacting components, the specific model parameters often depend on external conditions.The crucial point is, that under S < 0, a state only possible if P (k) is fat-tailed (i.e.β = 0), stability becomes asymptotically independent on C, driven solely by Ω, which is insensitive 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.This suggests 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, 1, −1, 0), and hence Eq. ( 10) predicts S = −β.Implemented on a scale-free network (β > 0), Social has S < 0, and is thus asymptotically stable.However, the same dynamics on a bounded P (k) (β = 0) becomes sensitively stable, i.e. S = 0, and therefore its stability depends on the specific value of C, which can be tuned by changing the infection/recovery rates.Consequently, while generally Social's pandemic 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 threshold vanishes under a scale-free P (k) 54 .The crucial point is, however, that while the original result, particular to the SIS model, was obtained using a dedicated model-specific analysis, our control parameter S allows us to systematically analyze the stability profile of all combinations of topology/dynamics within (6), using a single universal classification.
To observe this role of P (k), beyond the specific example of Social, we consider again E(A ij , Ω)'s principal 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, feedback, driven by the parameter C in J's diagonal (7); Fig. 4b.It is therefore, natural to enhance stability by increasing C, which, in effect, translates to strengthening each node's intrinsic negative feedback.Equation ( 9) predicts that a system becomes stable if C exceeds a critical value beyond which λ becomes negative.For asymptotically stable systems (S < 0) we have C 0 → 0, a guaranteed stability even under arbitrarily small C. In contrast, for asymptotically unstable systems (S > 0) we have C 0 → ∞, hence such systems are impossible to stabilize even under extremely large C, i.e. extreme levels of negative feedback.The point is, that if P (k) is bounded we have, always, S = 0 (sensitive stability), hence C 0 is finite, and the system's stability can be controlled by tuning the parameter C, e.g., changing the environmental conditions.
In Fig. 4c-k we extract again a set of three J matrices from E(A ij , Ω), representing systems from our three stability classes: J AS , asymptotically stable with Ω = (2, 2, 2, −1); J SS , sensitively stable with Ω = (0, −1, −2, 0); and J AU , asymptotically unstable with Ω = (1, −2, −1, 2).For each of these we plot C 0 vs. N , capturing the level of negative feedback required to ensure the system's stability.Under an Erdős-Rényi network, with bounded P (k) (β = 0) we do not observe a defined asymptotic behavior.The critical C 0 does not scale with N , indicating that sufficient perturbation to the model parameters can affect the system's stability (Fig. 4i-k, green circles).
The picture, however, fundamentally changes when we construct the same J matrices from a scale-free A ij with β = 0.6 (Fig. 4f-h), squares).Now we have S = −1.2,0 and 1.8 for J AS , J SS and J AU , respectively.As predicted in (11), under these condition, we observe a clear asymptotic behavior, in which C 0 scales with N .For J AS we have C 0 ∼ N −1.2 , while under J AU we observe C 0 ∼ N 1.8 (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 can be stabilized under finite C 0 , independently of system size N (Fig. 4g).
Hence, we observe a qualitative difference between bounded vs. scale-free P (k), in which degree heterogeneity can potentially afford the network a guaranteed stability, that is asymptotically independent of microscopic parameters.
Stabilizing vs. destabilizing dynamics.To gain qualitative insight on the role of the nonlinear dynamics, beyond the formal prediction of Eq. ( 10), let us focus on a specific test case of Ecological interactions, in which Eq. ( 6) takes the form (Fig. 1g) The self-dynamics follows logistic growth 55 and the interaction mechanism describes mutualistic dynamics 45 .The response function F (x j ) is designed to capture the positive contribution of node 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 dynamics lead to Ω = (0, 1, 1, 0), which in (10) provides S = β > 0, predicting an asymptotically unstable dynamics.The source of this instability is rooted in the strong positive feedback between 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 the limit 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, 1, 1, −2), predicting S = −β, and therefore, the same network which was unstable 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.This is expressed formally in (8) through the suppressed off-diagonal terms J ij , here -due to the negative exponent ρ = −2.This reduction in J ij , in turn, leads to a diminished λ.In that sense, the stability classifier S in (10) captures precisely this balance: if ν, ρ are large -the off-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 case S = 0 captures a balance between J ii and J ij , rendering the system sensitive to its microscopic parameters.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 system is either stable or not, we can use S to quantify the level of the system's stability.Consider again Eq. ( 11), capturing the critical value of C in (7) required to stabilize the system's dynamics.An asymptotically stable system has C 0 → 0, i.e. stability is ensured even under arbitrarily small C.
In contrast, an asymptotically unstable system has C 0 → ∞, hence cannot be stabilized under finite 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 = −β.Both are asymptotically stable, yet Eq. ( 11) indicates that Ecological is more stable: to destabilize Ecological one needs C Eco 0 The meaning is that to destabilize Ecological we much reach a more extreme reduction in C than the one required for Regulatory.In Fig. 5b we show the lines of equi-stability, comprising all dynamics with the same S value, and hence a similar level of stability/instability.These equi-stable classes may group together highly distinct dynamics, that despite their different Ω may have identical S. For example, Social and Ecological, two dynamics from fundamentally different domains, that exhibit a similar level of stability, belonging to the class S = −β (Fig. 1g).
Mixed interactions.Our analysis up to this point focused on cooperative mechanisms, in which M 1 (x) and M 2 (x) are positive, as observed, e.g., in mutualistic ecological dynamics 7,28 .To complement this, in Supplementary Section 5 we numerically analyze the case of mixed dynamics, in which a fraction f of positive interactions coexists alongside a 1 − f fraction of negative ones.We find that for a broad range of f values, our predicted asymptotic stability continues to be observed.This is attributed, once again, to a naturally emerging organizing principle, relevant in many real-world systems: the fact that the activities x i (t) in ( 6) are, by definition, non-negative.Therefore, as the system evolves in time, at any instance where a node reaches x i (t) = 0 it becomes extinct, and permanently removed from the network.This ratchet-motion prunes out nodes with many negative interacting partners, as their activities are naturally driven towards zero.As a result, the system reaches a fixed-point in which the majority of remaining links is biased towards positive interactions.Under these conditions, the emergent Jacobian matrix J Mix , while still featuring some level of negative interactions, is predominantly positive.Therefore, despite the fact that, strictly speaking, J Mix / ∈ E(A ij , Ω), it continues to approximately follow the patterns of ( 7) and ( 8), and most importantly in the present contextcontinues 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 complex system.Here, we exposed distinct patterns in the structure of J that (i) arise from the nature of the system's interaction dynamics; (ii) affect its principal eigenvalue, and thus, its stability profile.These patterns are expressed through the four dynamic exponents Ω = (η, µ, ν, ρ) in ( 7) and ( 8), linked analytically to the system's dynamics, M 0 (x), M 1 (x), M 2 (x), through the leading powers of the expansions in Eq. ( 14).This dependence on the powers (Ψ(n), Γ(n), etc.), rather than the coefficients (G n , C n , etc.) indicates that Ω is hardwired into the interaction dynamics.Indeed, the powers in (14) are determined by the dynamic model, e.g., Social or Regulatory, but not the specific model parameters, which only impact the coefficients.Therefore, our predicted Jacobian ensemble in (7) and ( 8), as well as its associated stability classifier S in (10), 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 structure into dynamic predictions [56][57][58] .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 a severe and troubling clash between spectral graph theory and empirical observation.Here we reconcile this contradiction, by showing that, in effect, we have, all this time, been looking at the wrong 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 offer to apply spectral analysis, not to the topology A ij , but rather to the derived Jacobian in ( 7) and ( 8), 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 theoretical obstacles along the path to predict, understand and control complex system behavior [59][60][61][62] .The problem is that we lack sufficient theoretical tools to treat the challenging combination of nonlinear dynamics with complex, random and often highly heterogeneous network structures.Consequently, advances are often system specific, requiring dedicated tools, and hence, lacking the breadth of a general network dynamics theory.The Jacobian ensemble presented here, transforms the nonlinear Eq. ( 6) into an effective linear system, whose weights depend on both topology and dynamics.It can, therefore, allow us to utilize the powerful tools of spectral graph theory, even in a nonlinear dynamic setting.This may offer a basis for systematic analysis of nonlinear network dynamics.

BOX I. THE ENSEMBLE E(A ij , Ω) -HOW TO CALCULATE Ω
To construct J as appears in Eqs.(7) and (8) we must extract the exponents Ω = (η, µ, ν, ρ) from Eq. ( 6).First we use the dynamic functions M 0 (x), M 1 (x) and M 2 (x) to construct three new functions where R ′ (x) represents a derivative dR/ dx around the fixed point x.From these we extract four additional functions, which we express through a Hahn power-series expansion as where R −1 (x) and Z −1 (x) represent the inverse functions of R(x) and Z(x).The Hahn expansion 63 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 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 where the constant C > 0 captures the system's specific dynamic time-scales.Note that Ω, as opposed to C, depends only on the leading powers of ( 14), and not on the coefficients.Therefore it is unaffected by Eq. ( 6)'s microscopic parameters, capturing an intrinsic characteristic of 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 Supplementary Section 2.  C).(c) The principal eigenvalue, under this construction, diverges with the number of interacting components N , predicting that a large system cannot be stable.This non-realistic prediction captures the unresolved diversity-stability paradox.(d) Our alternative J-ensemble, E(A ij , Ω), is designed to capture the roles of both the network (A ij , β, grey) and the nonlinear dynamic mechanisms of interaction, e.g., Biochemical, Social or Regulatory (orange).These mechanisms are described by M 0 (x), M 1 (x), M 2 (x) in Eq. ( 6), 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, but rather predicted by Eqs. ( 7) and ( 8).(e) As a result, the same A ij may lead to different J matrices, as different dynamics are characterized by distinct Ω. (f) Consequently, λ follows a different form than in (c), with a variable stability classifier S, whose sign predicts the system's stability: Unstable (red) in which diversity, indeed, leads to instability; Sensitive (green), where diversity plays no role, and Stable (blue), in which diversity enhances stability.Hence, under real nonlinear dynamics a large system can (S = 0), and in some cases even must (S < 0) be stable.(g) We use four common dynamic models capturing Social, Regulatory, Ecological and Biochemical interactions to examine our theoretical analysis.For each of these dynamics we extract Ω and S, as detailed in Supplementary Section 2. * For the calculation of S in Biochemical see Supplementary Section 3.2.and Regulatory (E(A ij , Ω 2 ), 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 to the nonlinear interaction mechanisms (Biochemical vs. Regulatory).Also here we find that while the random J is unstable (λ = 6.23, 10.3 > 0), our dynamic J matrices all have λ < 0, excluded from the diversity-stability paradox.In all panels when displaying J each entry (pixel) represents an average over a 50 × 50 block of the complete Jacobian matrix.Details on all 6 empirical networks are provided in Supplementary Section 4. We extracted 2, 077 Jacobian matrices from the E(A ij , Ω) ensemble, using combinations of model/real networks with different dynamic exponents Ω.For each J we calculated the principal eigenvalue λ and the stability classifier S in (10).(a) λ vs. S for all 2, 077 J-matrices.We observe the three predicted classes: Asymptotically unstable (red) in which S > 0 and hence also λ > 0; Sensitively stable (green), where S = 0 and λ can be both positive or negative; Asymptotically stable (blue), where S < 0 and therefore λ < 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) all our empirical networks are classified as unstable (red symbols).However, the same networks under 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 accounting for the nonlinearity through Ω.(b) The value of λ emerges from the competition between J's off-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, seeking the critical C 0 , in which λ becomes negative (grey lines).(f) For the stable J AS (S = −1.2) we find that C 0 decreases with N (squares), capturing the asymptotic stability, in which for large systems (N → ∞) stability is sustained even under arbitrarily small C. The theoretical scaling predicted in Eq. ( 11) is also shown (solid line, slope −1.2).(g) For J SS we have S = 0, the critical C 0 is independent of N , hence the system's stability can be affected by finite changes to its dynamic parameters, e.g., environmental perturbations.(h) The asymptotically unstable J AU (S = 1.8) has C 0 → ∞ in the limit of large N , in perfect agreement with Eq. (11) (solid line).(i) -(k) For a bounded P (k), e.g., Erdős-Rényi, β vanishes and hence S = 0 in (10).Under these conditions, regardless of Ω, the system is always sensitively stable and therefore C 0 does not scale with N .This demonstrates the role of degree-heterogeneity in ensuring stability, in the face of changing environmental conditions.The stability classifier S (10) comprises five exponents: β characterizes the network topology and its degree heterogeneity via (4); ν, ρ determine the magnitude of J's off-diagonal terms; µ, η characterize the self-dynamic diagonal terms.This portrays stability as a balance between the positive feedback mediated by the system's interactions, strengthened with increasing ν, ρ, vs. the stabilizing effect of the diagonal terms, that are reinforced by η, µ.(b) S vs. the exponents Ω = (η, µ, ν, ρ).Increasing ν, ρ drives the system towards the unstable class (S > 0), while increasing µ, η enhances the system's stability.In addition to the discrete classification of stable dynamics (blue, S < 0), unstable dynamics (red, S > 0) and sensitive dynamics (green, S = 0), S also helps quantify the continuum level of stability, where the more negative/positive is S, the more stable/unstable is the system.Grey lines represent equi-stable systems.

Figure 1 :
Figure1: Solving the diversity-stability paradox.(a) A complex system is captured by an interaction network A ij , whose diversity can be quantified by β in(4).In May's random matrix framework, to assess dynamic stability, we construct the Jacobian J with weights extracted independently from the distribution P (w), and diagonal term set to −C (constant).(b) The resulting stability matrix J ij , belonging to the ensemble E(A ij , P (w), C).(c) The principal eigenvalue, under this construction, diverges with the number of interacting components N , predicting that a large system cannot be stable.This non-realistic prediction captures the unresolved diversity-stability paradox.(d) Our alternative J-ensemble, E(A ij , Ω), is designed to capture the roles of both the network (A ij , β, grey) and the nonlinear dynamic mechanisms of interaction, e.g., Biochemical, Social or Regulatory (orange).These mechanisms are described by M 0 (x), M 1 (x), M 2 (x) in Eq. (6), 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, but rather predicted by Eqs.(7) and (8).(e) As a result, the same A ij may lead to different J matrices, as different dynamics are characterized by distinct Ω. (f) Consequently, λ follows a different form than in (c), with a variable stability classifier S, whose sign predicts the system's stability: Unstable (red) in which diversity, indeed, leads to instability; Sensitive (green), where diversity plays no role, and Stable (blue), in which diversity enhances stability.Hence, under real nonlinear dynamics a large system can (S = 0), and in some cases even must (S < 0) be stable.(g) We use four common dynamic models capturing Social, Regulatory, Ecological and Biochemical interactions to examine our theoretical analysis.For each of these dynamics we extract Ω and S, as detailed in Supplementary Section 2. * For the calculation of S in Biochemical see Supplementary Section 3.2.

Figure 3 :
Figure3: The dynamic J ensemble in empirical networks.(a) Using the real social network, UCIonline, we extracted the Jacobian J from the random ensemble E(A ij , P (w), Ω) (left).The resulting J matrix has λ = 25.4 > 0, and is, therefore unstable, as per May's prediction.Implementing our Social dynamics (SIS) we constructed J again, this time using our predicted E(A ij , Ω) (right).The resulting J is visibly different, owing to the emergent patterns of (7) and(8).Most importantly, it has λ = −0.80< 0, and hence despite the random matrix prediction, the system is, in fact, stable.(b) Similar results are obtained from Epoch, with λ = 49.6 > 0 under E(A ij , P (w), Ω) (left), vs. λ = −0.94under 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 be stable (λ = −16.0 and −289).(e)-(f) The biological networks PPI1 and PPI2 are matched with three potential J matrices: random (E(A ij , P (w), C), left), Biochemical (E(A ij , Ω 1 )), center) and Regulatory (E(A ij , Ω 2 ), 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 to the nonlinear interaction mechanisms (Biochemical vs. Regulatory).Also here we find that while the random J is unstable (λ = 6.23, 10.3 > 0), our dynamic J matrices all have λ < 0, excluded from the diversity-stability paradox.In all panels when displaying J each entry (pixel) represents an average over a 50 × 50 block of the complete Jacobian matrix.Details on all 6 empirical networks are provided in Supplementary Section 4.

Figure 4 :
Figure4: Three classes of dynamic stability in real and model networks.We extracted 2, 077 Jacobian matrices from the E(A ij , Ω) ensemble, using combinations of model/real networks with different dynamic exponents Ω.For each J we calculated the principal eigenvalue λ and the stability classifier S in(10).(a) λ vs. S for all 2, 077 J-matrices.We observe the three predicted classes: Asymptotically unstable (red) in which S > 0 and hence also λ > 0; Sensitively stable (green), where S = 0 and λ can be both positive or negative; Asymptotically stable (blue), where S < 0 and therefore λ < 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) all our empirical networks are classified as unstable (red symbols).However, the same networks under 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 accounting for the nonlinearity through Ω.(b) The value of λ emerges from the competition between J's off-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, seeking the critical C 0 , in which λ becomes negative (grey lines).(f) For the stable J AS (S = −1.2) we find that C 0 decreases with N (squares), capturing the asymptotic stability, in which for large systems (N → ∞) stability is sustained even under arbitrarily small C. The theoretical scaling predicted in Eq. (11) is also shown (solid line, slope −1.2).(g) For J SS we have S = 0, the critical C 0 is independent of N , hence the system's stability can be affected by finite changes to its dynamic parameters, e.g., environmental perturbations.(h) The asymptotically unstable J AU (S = 1.8) has C 0 → ∞ in the limit of large N , in perfect agreement with Eq. (11) (solid line).(i) -(k) For a bounded P (k), e.g., Erdős-Rényi, β vanishes and hence S = 0 in(10).Under these conditions, regardless of Ω, the system is always sensitively stable and therefore C 0 does not scale with N .This demonstrates the role of degree-heterogeneity in ensuring stability, in the face of changing environmental conditions.

Figure 5 :
Figure 5: The ingredients of dynamic stability.(a)The stability classifier S (10) comprises five exponents: β characterizes the network topology and its degree heterogeneity via (4); ν, ρ determine the magnitude of J's off-diagonal terms; µ, η characterize the self-dynamic diagonal terms.This portrays stability as a balance between the positive feedback mediated by the system's interactions, strengthened with increasing ν, ρ, vs. the stabilizing effect of the diagonal terms, that are reinforced by η, µ.(b) S vs. the exponents Ω = (η, µ, ν, ρ).Increasing ν, ρ drives the system towards the unstable class (S > 0), while increasing µ, η enhances the system's stability.In addition to the discrete classification of stable dynamics (blue, S < 0), unstable dynamics (red, S > 0) and sensitive dynamics (green, S = 0), S also helps quantify the continuum level of stability, where the more negative/positive is S, the more stable/unstable is the system.Grey lines represent equi-stable systems.

Figures
Figures

Figure 3 The
Figure 3

Figure 4 Three
Figure 4