Network structural origin of instabilities in large complex systems.

A central issue in the study of large complex network systems, such as power grids, financial networks, and ecological systems, is to understand their response to dynamical perturbations. Recent studies recognize that many real networks show nonnormality and that nonnormality can give rise to reactivity-the capacity of a linearly stable system to amplify its response to perturbations, oftentimes exciting nonlinear instabilities. Here, we identify network structural properties underlying the pervasiveness of nonnormality and reactivity in real directed networks, which we establish using the most extensive dataset of such networks studied in this context to date. The identified properties are imbalances between incoming and outgoing network links and paths at each node. On the basis of this characterization, we develop a theory that quantitatively predicts nonnormality and reactivity and explains the observed pervasiveness. We suggest that these results can be used to design, upgrade, control, and manage networks to avoid or promote network instabilities.

INTRODUCTION

The dynamical stability of large complex network systems is an intriguing problem. The basic question of what properties of such systems govern their stability has attracted much interest, which was initially sparked by a 1972 article by Robert May predicting that sufficiently large systems should be linearly unstable despite the observed stability of large real ecological systems (). While the literature on this problem and discrepancies between theory and observation has focused mostly on May’s original context [ecological networks (–), including microbiome communities (–)], the problem is relevant for large network systems in general, including financial networks (, ), power networks (), and immune system networks (). The problem acquires a new dimension when the Jacobian matrix M determining the linear stability is nonnormal () (i.e., MM ≠ M, where M denotes the transpose of M). This is because a small perturbation in such a system can cause the resulting state deviation to initially grow and become large enough to excite nonlinear instabilities, even when the system is linearly stable (–). The system’s capacity to exhibit initial growth of deviations in the linear regime is termed reactivity (), which is known to relate to a spectral property of the matrix M. Although the initial interest in nonnormality and reactivity emerged in hydrodynamics (–), these properties have recently gained attention in the study of network systems, including ecological (, ), neuronal (, ), chemical reaction (), and communication networks (), as well as in the study of pattern formation in networks (, ) and control of networks (). The literature has begun to reveal how prevalent nonnormality and reactivity are in real-world networks (), but the fundamental question of which network structural mechanisms underlie the apparent prevalence of nonnormality and reactivity has not yet been addressed.

Here, we address this question by deriving rigorous conditions for nonnormality and reactivity that can be applied to any directed network [rather than to the average over an ensemble of networks ()] and interpreted in terms of the structure of the given network. For nonnormality, the condition is that there is an imbalance between incoming and outgoing links at a node or a pair of nodes in terms of their numbers and/or weights. For reactivity, the condition is that there is an imbalance between the eigenvector centrality of a node associated with incoming network paths (including their weights) and the eigenvector centrality associated with outgoing paths. We use these conditions to show that, in a broad class of directed networks, the probability that the coupling matrices are both nonnormal and reactive approaches 1 quickly as the network size increases. We prove our results for large networks using a general network model that permits arbitrary distributions of possibly correlated in- and out-degrees (the number of incoming and outgoing links at a node, respectively) and arbitrary distributions of link weights. We also validate the prevalence of nonnormality and reactivity using a dataset of 251 real networks. This set is the largest and most diverse collection of directed networks, which also includes the largest networks (with up to nearly 8 million nodes), ever considered in this context.

These findings indicate that no additional global organization of connectivity is necessary to generically observe nonnormality and reactivity. This is important given that large-scale structures can have dynamical consequences, such as the stability promoted by trophic coherence in food-web networks () and the nonmonotonicity supported by linear chain structures in chemical reaction networks (). In addition to establishing the prevalence of nonnormality and reactivity, we develop a quantitative theory that directly relates the extent of the degree and centrality imbalances in a given network to the extent of nonnormality and reactivity, respectively. Thus, our results reveal the network structural features responsible for nonnormality and reactivity, contributing to the much needed fundamental understanding of the relationship between dynamical and structural properties of directed networks (–).

RESULTS

Nonnormality and reactivity of network systems

Given that many real network systems operate near an equilibrium, we consider the class of (nonlinear) systems whose linearization around a given reference equilibrium state is described byfor i = 1, … , n, where x is the (scalar) deviation from the reference state for the ith node. The matrix A = (A) can be regarded as the weighted adjacency matrix of the system’s directed interaction network: A ≠ 0 if node j is connected to node i, and A = 0 otherwise. The parameter α and the diagonal element A represent the node dynamics and any self-link at node i, respectively. We assume that the time scales of the dynamics in Eq. 1 are much shorter than those of the evolution of the interaction network structure, so that A can be regarded as constant. For concreteness, we also assume α = α for all i in the following unless otherwise indicated. We note, however, that our results on random networks are valid for heterogeneous α and that the presence of the heterogeneity is generally expected to increase both the nonnormality and the reactivity of the Jacobian matrix of the system (see Materials and Methods for details). For α = α, the Jacobian matrix M and the adjacency matrix A are related as M = A − αI, with I denoting the n × n identity matrix. Thus, the reference state is asymptotically stable if and only if Re λ1(M) < 0, or Re λ1(A) < α, where λ1(X) denotes the eigenvalue with the largest real part for any matrix X. Assuming A ≥ 0 for i ≠ j, as observed in many real networks, λ1(A) is guaranteed to be real by the Perron-Frobenius theorem for nonnegative matrices (), regardless of whether the network is strongly connected (i.e., any two nodes are connected by directed paths in both directions). While we have implicitly assumed one-dimensional node dynamics in Eq. 1 for clarity, we also establish an exact relation between the Jacobian and adjacency matrices for a broader class of systems. In particular, this relation shows how the nonnormality and reactivity of the adjacency matrix generically imply the same properties for the Jacobian matrix (see Supplementary Materials, section S1, for details).

We first show how the nonnormality and reactivity of the system in Eq. 1 can be expressed as properties of the network structure described by A under the uniform α assumption. Noting that matrix D ≔ MM − M represents the deviation of the Jacobian matrix M from being normal, the nonnormality of the system can be quantified by the Frobenius norm (). Since D = MM − M = AA − A (and hence does not depend on α), the nonnormality of M is reduced to the nonnormality of A. To characterize the reactivity of the system in Eq. 1 mathematically, we consider the maximum exponential rate of initial growth of the state deviation vector that can result from a perturbation of the initial state. This rate is given by λ1((M + M)/2) = λ1(H) − α, where H ≔ (A + A)/2 is the symmetric part of A. Thus, the state deviation can grow if λ1(H) > α. Combining this with the stability requirement λ1(A) < α, we define reactivity as a property of the interaction network structure A (including its dependence on the reference state): A is said to be reactive if there are α values for which λ1(H) > α > λ1(A). That is, A is reactive if, within the linear regime, the system can be both stable and capable of exhibiting initial growth of state deviations. Such growth can push the system state out of the region in which the linearization in Eq. 1 is valid, potentially inducing nonlinear instabilities. The reactivity condition can be expressed asand thus we use λΔ(A) as a measure of the reactivity of the interaction network structure A [this definition can be extended to allow for negative weights by replacing λ1(A) with Re λ1(A)]. This measure is independent of α, in contrast to the reactivity for a specific α, which could be defined as λ1(H) − α following previous studies (, , , ). In the general case of heterogeneous α, both nonnormality and reactivity can be defined in the same way after absorbing the heterogeneity into the diagonal elements of the matrix A (see Materials and Methods for details).

We note that λΔ(A) ≥ 0 always holds (see Eq. 7 in Materials and Methods) and λΔ(A) > 0 implies that the initial growth rate of state vector deviation in Eq. 1 for α = 0 can be strictly larger than λ1(A) when the initial state vector x0 is chosen to be the eigenvector corresponding to λ1(H), where we use ∥·∥ to denote the 2-norm and x := (x1,…,x) to denote the state vector. This is the case because . We also note that λΔ(A) is invariant under any coordinate transformation if we concurrently apply the transformation to the observable for the system in Eq. 1, since λΔ(A) is based on the value of the observable rather than its coordinate-specific representation. For the reactivity measure defined in Eq. 2, the observable is the 2-norm of the system state vector. This choice is standard in the literature [see, e.g., (, –, , , , )], as it permits a convenient characterization through eigenvalues, but for certain network processes, different measures of deviations may be more natural and would lead to different definitions of reactivity [see, e.g., the 1-norm used in (–) and the absolute value of a one-dimensional projection used in ()].

We used the conditions and measures just defined to study nonnormality and reactivity in a large dataset of 251 real directed networks, which consists of 63 biological, 51 informational, 79 social, 39 technological, and 19 economic/game networks and avoids repetition of similar networks. The dataset used here substantially expands on an earlier study () with respect to the number of networks, the largest network size, and the diversity of network types. We verified that A is nonnormal (i.e., ∥D∥F > 0) for all 251 networks and is reactive [i.e., λΔ(A) > 0] for all but one network (see Materials and Methods for details on the network data and our findings). Quantitatively, Fig. 1A reveals a strong positive correlation between the level of nonnormality ∥D∥F and the overall level of imbalance between the nodes’ in- and out-degrees, measured by the average , where and are the in- and out-degrees of node i, respectively (i.e., the numbers of incoming and outgoing links from/to other nodes). Similarly, Fig. 1B shows a monotonic relation between the reactivity λΔ(A) and a measure of imbalance, in this case between the eigenvector centrality associated with the incoming and outgoing paths, quantified by the angle θ1 between the left and right eigenvectors corresponding to λ1(A) (to be precisely defined below). Both relations are observed for each type of networks (fig. S1) and will be theoretically derived and computationally validated in a later section.

Fig. 1.

Nonnormality and reactivity of real networks.

(A) Nonnormality versus the in- and out-degree imbalance quantified by the average 〈(din − dout)2〉 of the squared difference between and over all nodes i for the 251 real networks considered (see Materials and Methods for details on the data). The measure of nonnormality ∥D∥F is normalized by based on the scaling predicted by our theory for random networks, where is the mean squared weight over all links in the network and is the average degree; see Eq. 22. The size and color of a circle indicate the network size and link density, respectively. (B) Reactivity versus the in- and out-centrality imbalance quantified by the angle θ1 (in radian) between the left and right eigenvectors associated with the largest eigenvalue of A. The measure of reactivity λΔ(A) is normalized by λ1(A) based on our theory in Eq. 5. The plot shows only the 212 nondegenerate reactive networks [i.e., those having nondegenerate λ1(A) and satisfying λΔ(A) > 0; see Materials and Methods for details]. We observe appreciable levels of nonnormality and reactivity in a vast majority of the networks: for 248 networks (≈99%) and λΔ(A)/λ1(A) > 0.01 for 188 nondegenerate networks (≈88%). IL, Illinois.

Imbalance conditions for nonnormality and reactivity

To understand the mechanisms underlying the correlations observed above, we first examine the network structural features that are responsible for the nonnormality and reactivity of A. The measure of nonnormality introduced above can be expressed aswhere , , , and . The variables and can be regarded as generalized in- and out-degrees, respectively, since they reduce to and for unweighted networks (i.e., if A ∈ {0,1}). For distinct nodes i and j, the variable () can be interpreted as a further generalization of the in-degree (out-degree) to a pair of nodes, and it reduces to the number of common in-neighbors (out-neighbors) shared by the two nodes in the case of unweighted networks. Equation 3, despite being immediate from the definition of D and its Frobenius norm, provides an insightful decomposition of nonnormality into two types of imbalances between incoming and outgoing links: The first term is the square sum of the generalized in- and out-degree differences at individual nodes, while the second is an analogous square sum for node pairs. From Eq. 3, we see that A is normal if and only if the generalized in- and out-degrees are equal for each node and for each pair of nodes (see Fig. 2A for an illustrative example). It is thus sufficient to have just a single node whose in- and out-degrees differ in order to make A nonnormal. However, even without such a node, A can still be nonnormal if there is a node pair with an imbalance between their common weighted in-neighbors and the common weighted out-neighbors (making the second sum nonzero in Eq. 3), as illustrated in Fig. 2B. Furthermore, Eq. 3 shows that a larger total imbalance between incoming and outgoing links implies a more nonnormal A.

Fig. 2.

Topological and spectral features inducing nonnormality and reactivity.

The four-node unweighted example networks shown illustrate all three possible cases (the normal-reactive case is not possible). The eigenvector in-centrality v1 and the eigenvector out-centrality u1, normalized so that ∥v1∥ = ∥u1∥ = 1, are indicated near each node i as (v1, u1). Note that we have and as the networks are unweighted. (A) Network for which A is normal (and thus nonreactive) because the in- and out-degrees match () for all nodes i and the number of shared in- and out-neighbors match () for all pairs of nodes i and j ≠ i, as illustrated for i = 1, j = 4 (blue, shared in-neighbor; red, shared out-neighbor). (B) Network for which A is nonnormal because the second sum in Eq. 3 is nonzero due to for some i ≠ j, e.g., . However, A is nonreactive because the in- and out-centralities are balanced at each node, i.e., v1 = u1 for all i. (C) Network for which A is nonnormal and reactive. Nonnormality in this case is guaranteed not only by for some i ≠ j (e.g., for i = 1, j = 4) but also by for some i (blue node, ; red node, ). Reactivity is guaranteed by v1 ≠ u1 for all nodes (blue numbers, v1 > u1; red numbers, v1 < u1).

For reactivity, a stronger condition is needed, since nonnormal A does not have to be reactive, as illustrated by the network in Fig. 2B. In particular, for a d-regular network, defined as an unweighted network in which the in- and out-degrees are all equal to a constant integer d ≥ 0, the adjacency matrix A can be nonnormal but can never be reactive. To see this, we first note that λ1(A) is bounded between the minimum and maximum row sum of A because this eigenvalue is the spectral radius of A [see, e.g., theorem 8.1.22 in ()] and thus is equal to d for any d-regular network. Since the same argument applies to the symmetric part of A, we have λ1(H) = d, implying that λΔ(A) = d − d = 0, i.e., A is nonreactive. Nonetheless, for many d-regular networks, the second term in Eq. 3 is strictly positive, rendering A nonnormal.

Given that nonnormality is only a necessary condition for reactivity, we now present a condition guaranteeing reactivity: A is reactive if the left and right eigenspaces associated with its largest eigenvalue λ1(A) are distinct (see Materials and Methods for a proof). In the generic case in which these eigenspaces are one dimensional, the reactivity condition is equivalent to having a strictly positive angle θ1 between the right eigenvector v1 and left eigenvector u1. We always choose the acute angle so that θ1 ≤ π/2, where the quantity 1/∣cos θ1∣ is known as the eigenvalue condition number in the literature (). This condition is illustrated by the example network in Fig. 2C. While the nonorthogonality of right eigenvectors defining nonnormality has been shown to often lead to reactivity [e.g., in (, )], our condition θ1 > 0 captures reactivity more precisely, as it is equivalent to the specific nonorthogonality between the right eigenvector v1 and some other right eigenvector. The condition θ1 > 0 also translates to the existence of at least one node i for which v1 ≠ u1, where v1 is the eigenvector in-centrality, defined to be the ith component of the right eigenvector v1 (and hence associated with incoming paths to the node), and u1 is the eigenvector out-centrality, defined similarly through the left eigenvector u1 (and the outgoing paths). The inequality v1 ≠ u1 can thus be interpreted as an imbalance between the incoming and outgoing “flow” of centrality. In the example of d-regular networks above, if λ1(A) is nondegenerate (so that v1 and u1 are unique and θ1 is well defined), we find that v1 = u1 for all i and thus θ1 = 0, as expected from the result above that A is nonreactive for all d-regular networks. In general, the imbalances and v1 − u1 for individual nodes and for pairs of nodes can be either positive or negative and tend to be distributed heterogeneously across a given network, as illustrated by the three example real networks in fig. S2.

While we focus mainly on the adjacency matrix A, different interaction matrices can also be considered in Eq. 1. This includes the Laplacian matrix, which we show can be nonnormal even if A is not and vice versa (a fact that has been generally overlooked; see fig. S3 for illustrative examples).

Prevalence of nonnormality and reactivity

We now address the question of how often nonnormality and reactivity are expected to be observed by considering a model of random directed weighted networks with a given number of nodes n, a given joint probability distribution for the in- and out-degrees, and a given distribution of link weights. By allowing for arbitrary distributions for the in- and out-degrees, an arbitrary correlation between them, and an arbitrary distribution of weights, the model is capable of capturing the essential elements of real networks involved in the conditions established above for nonnormality and reactivity. This model is a generalization of the Chung-Lu-Vu model studied in () and allows us to specify a general joint distribution of in- and out-degrees (rather than the expected node degrees) and a general weight distribution. A network realization under this model is generated as follows. First, the (possibly non-integer) expected in- and out-degrees of each node are randomly drawn from the given joint distribution. Then, for each i and j, a directed link is created from node j to node i with a probability proportional to the product of the expected in-degree of node i and the expected out-degree of node j. Last, a random weight A is drawn from the given distribution for each link j → i. We note that a special case of this model produces unweighted networks, and self-links can be excluded if desired. For further details on this generalized Chung-Lu-Vu (GCLV) model, see Materials and Methods.

For weighted networks generated by the GCLV model, we show that the adjacency matrix A is almost always nonnormal in the limit of large network size n (see Supplementary Materials, section S2, for a proof for unweighted networks and its extensions to weighted and Laplacian-coupled networks). More precisely, we show that the probability of having at least one node whose generalized in-degree and out-degree are different (implying nonnormality of A, as discussed above) tends to one as n → ∞. In practice, the probability that A is nonnormal grows quickly and is very close to one even for networks with less than 10 nodes, as shown numerically for both weighted and unweighted networks in fig. S4A and table S1 for two classes of in-/out-degree distributions: (i) the gamma distribution, whose probability density function is p(x) ∼ xe−, x > 0, where a and b are parameters, and (ii) the Dirac delta distribution centered at d, which renders the model equivalent to the Erdős-Réyi (ER) networks with fixed mean degree d (and thus with n-dependent connection probability p = d/n). The same appears to hold true for other random network models, as verified in fig. S4A for the unweighted ER networks with fixed p and random d-regular networks and also verified in table S1 for the weighted versions of these random networks.

Nonnormality does not necessarily imply reactivity. For large networks, however, we can show that, if A is nonnormal, then it is also reactive in almost all cases. This is because the conditional probability ℙ(A is reactive ∣ A is nonnormal) ≥ ℙ(A is reactive), as reactivity implies nonnormality, and because we show that ℙ(A is reactive) approaches one as n → ∞ for the GCLV model (see Supplementary Materials, section S2, for a proof, including the case of Laplacian-coupled networks). More precisely, we prove that the probability of having distinct left and right eigenspaces associated with λ1(A) (sufficient for reactivity, as noted earlier) converges to one. When estimated numerically for finite n, the actual probability that A is reactive and the conditional probability that A is reactive given that it is nonnormal are again very close to one even for small n, as shown in fig. S4 (B and C), respectively (which also shows similar results for two other random network models). All these results support the observation from Fig. 1 that the reactivity of A is prevalent among real networks.

Quantitative characterization of nonnormality and reactivity

To understand the correlations observed for the real networks in Fig. 1, we now derive theoretical estimates of the nonnormality ∥D∥F and the reactivity λΔ(A) for random networks. For the nonnormality, we first consider networks generated by the unweighted GCLV model with no self-links (i.e., A ∈ {0,1} and A = 0 for all i and j) given a fixed set of in-degrees and out-degrees for all nodes. In this case, the first sum in Eq. 3 is fixed and proportional to the (constant) average single-node degree imbalance 〈(din − dout)2〉, while the second sum representing the imbalances at the node pair level is a random variable. Assuming n ≫ 1 and approximating the second sum with its expected value (see Materials and Methods for details), we havewhere 〈z〉 ≔ ∑/n denotes the average of z over nodes i and is the average degree. Thus, in the limit of large networks, nonnormality increases with the network size n as ∥D∥F ∼ n1/2, with an n-independent prefactor expressed as a simple function of the in- and out-degree imbalance 〈(din − dout)2〉, the in- and out-degree correlation 〈dindout〉, and the average degree . We also derive extensions of this formula to the weighted GCLV model (allowing self-links) and to Laplacian-coupled networks (see Materials and Methods), which have additional terms and factors involving the diagonal elements A and the statistics of link weights (see Eqs. 22, 32, and 36). For the reactivity λΔ(A), we make use of the observation that the spectral gap between the leading eigenvalue λ1(A) and the remaining eigenvalues λ2(A), …, λ(A) is often large for the adjacency matrix A of large random networks [which is the case, e.g., for the ER networks () and for the GCLV model when the mean degree is large ()]. For such networks, we expect λ1(H) to be well-approximated by λ1(H1), where is the symmetric part of the leading component and u1 and v1 are respectively the left and right (column) eigenvectors associated with the leading eigenvalue λ1(A), normalized so that ∥v1∥ = 1 and . Among the 212 real networks in Fig. 1B [which have nondegenerate λ1(A), ensuring that H1 is well defined], we indeed observe λ1(H) ≈ λ1(H1) for most, as shown in fig. S5A. For any (possibly weighted) network satisfying λ1(H) ≈ λ1(H1) and having nondegenerate λ1(A) (even when it is not strongly connected), we derive a simple expression:with a prefactor that depends monotonically on the angle θ1 between the leading left and right eigenvectors u1 and v1 (see Materials and Methods for a derivation).

We find that Eqs. 4 and 5 provide good approximations, as validated in Fig. 3 for several classes of random networks generated by the GCLV model (green, orange, and red dots) as well as for the real networks used in Fig. 1 (blue dots). In addition, they capture the general tendency observed in Fig. 1 for nonnormality and reactivity to increase with the degree imbalance and the eigenvector angle in real networks.

Fig. 3.

Validating theoretical predictions for nonnormality and reactivity.

(A) Computed nonnormality versus the corresponding prediction from Eq. 4 for the real networks in Fig. 1A (blue). Plots are also shown for random networks generated by the GCLV model: unweighted networks with the gamma or power-law distribution for the in- and out-degrees (green) and weighted networks with gamma-distributed degree distributions and exponentially distributed weights (orange; mean μ = 0.7, 1.0, or 1.3) or Gaussian-distributed weights (red; mean one and standard deviation (SD) σ = 0.1 or 0.2), keeping only positive weights. For each class of networks, we used 20 equally spaced values of n between 103 and 105 on the logarithmic scale. For each n, we show 20 realizations of the model, with the parameters of the distribution drawn randomly for each network realization (see Materials and Methods for details). The continuous lines correspond to Eq. 22. The plots for the random networks and the corresponding lines are shifted vertically to avoid overlapping. (B) Reactivity versus the eigenvector angle θ1. The continuous curves correspond to Eq. 5. The plot shows the same random networks as in (A) and the real networks in Fig. 1B, to which Eq. 5 is applicable.

Our results show that A is more nonnormal and reactive when the network is weighted. In particular, Eqs. 3 and 22 indicate that allowing for weights in random networks can only increase the probability of observing imbalances that induce nonnormality (i.e., the probability that ∥D∥F ≠ 0) and that a larger variance for the weights leads to a larger extent of nonnormality. Likewise, the eigenvector condition and Eq. 5 indicate that the presence of weights, whose randomness would be reflected in the eigenvectors, is expected to increase in-/out-centrality imbalances, leading to a higher probability that A is reactive and to a larger extent of reactivity.

The tendency for A to be dominated by the leading eigenvalue λ1(A) and the corresponding eigen-component A1 also explains why large directed networks are almost always both nonnormal and reactive. To see this, we first note that A1 satisfies the following exact relation linking its nonnormality and reactivity to each other and to the eigenvector angle:where we denote λ1 = λ1(A1) = λ1(A) for brevity. Here, we note that A1 and A have the same eigenvector angle θ1 associated with λ1, and we introduce to quantify the nonnormality of A1. From Eq. 6, it immediately follows that, for the leading component A1, nonnormality (∥D1∥F > 0), reactivity [λΔ(A1) > 0], and having a strictly positive eigenvector angle (θ1 > 0) are all mathematically equivalent to each other. Moreover, Eq. 6 shows that an increase in any one of these three measures implies an increase in all the other measures. Building on these observations, we find that random networks whose A is dominated by A1 generally satisfy Eq. 6 approximately with ∥D1∥F and λΔ(A1) replaced by ∥D∥F and λΔ(A), respectively, which implies that nonnormality and reactivity are approximately equivalent for such networks. This finite-n approximate equivalence complements the rigorous results we established above in the limit of large network size. For the real networks in Fig. 3B, Eq. 6 approximately holds true when ∥D1∥F and λΔ(A1) are replaced by ∥D∥F and λΔ(A), respectively (fig. S5, B to D).

For networks satisfying Eqs. 4 and 6 approximately, we see that, if λ1 = λ1(A) is bounded as n increases, nonnormality and reactivity would scale with n as ∥D∥ ∼ n1/2 and λΔ(A) ∼ n1/4, respectively. For the real networks, λΔ(A) tends to increase with n (even when normalized by 〈w2〉, the root mean square of the link weights), which is a trend observed even more strongly when the networks are randomized while holding the in- and out-degrees fixed (Fig. 4A). For random networks with a power-law degree distribution p(x) ∼ x−β and Gaussian-distributed link weights, the reactivity λΔ(A) scales with n with an exponent that depends on the power-law parameter β (Fig. 4B). The increase of λΔ(A) with n indicates that, as the system becomes larger, there will be a wider range of α for which the system can simultaneously exhibit a more pronounced transient response to a small perturbation and stronger linear stability. Specifically, Fig. 4C establishes that, if the node stability parameter α in Eq. 1 has an n-dependence α(n) ∼ nℓ with a constant ℓ, then there is a region in the β versus ℓ parameter space (shaded red) for which the maximum initial growth rate of perturbations given by λ1(H) − α(n) increases with n even though the first-order dynamics become more stable [i.e., λ1(A) − α(n) decreases]. Below this region the growth rate increases while linear stability decreases with n, which is doubly destabilizing, and above this region the opposite is observed.

Fig. 4.

Increase of reactivity with network size.

(A) Reactivity versus network size n for the real networks in Fig. 3B (blue circles). We normalize λΔ(A) by to facilitate comparison between networks with different scales for link weights. The black line indicates the least-squares fit (on the logarithmic scale), reflecting an increasing trend for the reactivity. Red dots indicate randomized versions of the real networks (see Materials and Methods for details), which more closely follow the trend line. (B) Scaling of λΔ(A) with n for the GCLV model with the same discrete power-law degree distribution p(x) ∼ x−β used in Fig. 3, but for a given value of β and uncorrelated and , with link weights drawn from the positive domain of the Gaussian distribution with mean one and SD 0.1. For each β, the plot shows the average over 200 network realizations. We observe that the scaling exponent for λΔ(A) depends on the parameter β of the degree distribution. (C) Reactivity−stability phase diagram for the system in Eq. 1 with an n-dependent parameter α = α(n) ∼ nℓ, where ℓ is a constant, and for the same power-law distribution with parameter β as in (B). Shaded in red is the region in which the maximum growth rate of state deviations λ1(H) − α(n) increases with n despite the concurrent increase of linear stability [i.e., decrease of λ1(A) − α(n)]. See Materials and Methods for details on the procedure to identify the boundaries of this region.

DISCUSSION

Our demonstration that nonnormality and reactivity are stronger and more prevalent for larger networks implies that conventional modal stability analysis alone is not sufficiently informative: The system can be less stable against perturbations than the modal analysis indicates, and this gap can grow with the system size. Even if such a system is linearly stable and the perturbation is well within the linear regime, the transient response may bring the system sufficiently far from the equilibrium that nonlinear instabilities can be induced, possibly in the form of a cascade (). Our findings thus suggest that, following a perturbation, there is a parameter region for which larger networks can exhibit larger oscillations in the linearized system and thus the possibility of a rare but substantial instability in the nonlinear system. Does such an instability cause a permanent transition to a different state? If so, to which state does the system transition? The answers to these questions depend on the size and direction of the perturbation, the extent of reactivity in the system, and the global structure of the state space. The latter, in particular, requires specific knowledge of the nonlinearity of the system under consideration.

Given our validation with extensive network data and general mathematical analysis, these conclusions are applicable to a wide range of real network systems and are thus not limited to the networks previously considered in addressing the complexity–stability problem. In particular, knowing that nonnormality and reactivity can be caused by the imbalance of incoming/outgoing connections or of the node’s in-/out-centrality will likely be useful in designing large complex technological (and possibly synthetic biological) systems. Importantly, our results show that nonuniform link weights and self-links, which are aspects of complexity often neglected in the systematic study of networks, tend to increase nonnormality. While here we focused on reactivity as a widely used measure of transient response, we anticipate similar results for other measures, such as pseudospectra (), singular values (), and the perturbation-averaged evolution of the state vector (). We also suggest that the structural characterization of nonnormality and reactivity can be relevant for the study of generalized networks of current interest, such as those accounting for temporal, multilayer, and higher-order (nonpairwise) interactions. Ultimately, the demonstration that nonnormality and reactivity tend to be more prevalent and more consequential for stability as the system grows in size opens up new vistas of network dynamics in complex systems.

MATERIALS AND METHODS

Impact of αi heterogeneity on nonnormality and reactivity

Any heterogeneity in α can be subtracted from the first term and absorbed into the second in Eq. 1 by replacing α and A with and A + (α − α), respectively. In Eq. 31, this would contribute to the heterogeneity of β, extending the validity of the nonnormality approximation in that equation to the general case of heterogeneous α. It then follows from Eqs. 23 and 31 that the expected nonnormality ∥D∥F for the GCLV model generally increases with the heterogeneity of α. Combining this with the approximate relation between nonnormality and reactivity in Eq. 6, we see that the reactivity also tends to increase with the heterogeneity of α in random networks. Moreover, the modification of α and A in Eq. 1 used here can also be applied to Eq. 5, which would extend the reactivity approximation formula to heterogeneous α. The conclusion above on the impact of α heterogeneity in random networks holds true in particular when α and A are uncorrelated, while specific correlations can, in principle, lead to a reduction in the heterogeneity of β and thus in .

After absorbing α heterogeneity into A in Eq. 1 as in the previous paragraph, the nonnormality and reactivity of the system can be defined in the same fashion, i.e., as ∥D∥F = ∥AA − A∥ and λΔ(A) = λ1(H) − λ1(A), respectively, but using the modified A. While this definition allows us to derive results for random networks (as described in the previous paragraph), the real networks were analyzed under the uniform α assumption since the values of α were not available in the dataset.

Sufficient condition for reactivity

To establish that having distinct (real) left and right eigenspaces associated with λ1(A) implies reactivity [i.e., λΔ(A) > 0], we will prove the contrapositive: λΔ(A) = 0 implies that these eigenspaces coincide. Suppose that λΔ(A) = 0, which implies λ1(H) = λ1(A). We will seek to show that any right eigenvector is also a left eigenvector and vice versa. Let v be a (real) right eigenvector of A associated with the eigenvalue λ1(A). Without loss of generality, we can assume v to be normalized so that v = 1. Then, we havewhere we note that v = v(A + A)v/2 = λ1(A)v = λ1(A). Since λ1(H) = λ1(A), the inequality in Eq. 7 becomes an equality, implying that v is a solution of the maximization problem and thus a (right) eigenvector of H associated with λ1(H), i.e., Hv = λ1(H)v. This, together with 2H = A + A and Av = λ1(A)v, yields v = (2λ1(H) − λ1(A))v = λ1(A)v. We thus conclude that v is a left eigenvector of A associated with λ1(A), in addition to being a right eigenvector. To show the opposite direction, we now let v be any left eigenvector of A associated with the eigenvalue λ1(A). This again implies v = λ1(A) and turns Eq. 7 into an equality, implying that v is a (left) eigenvector of H associated with λ1(H), i.e., v = v1(H). From this, an argument similar to the one above shows that v is also a right eigenvector of A associated with the eigenvalue λ1(A). Therefore, we conclude that the left and right eigenspaces of A corresponding to λ1(A) must coincide, which completes the proof.

Real network data

We retrieved raw data for 534 distinct directed networks available from Koblenz Network Collection (KONECT) (, ), the Netzschleuder Network Catalogue and Repository (), the Colorado Index of Complex Networks (ICON) (), and the Matrix Market Repository (). From KONECT, Netzschleuder, and ICON, we obtained all distinct directed non-bipartite networks that were available in the form of adjacency lists. From the Matrix Market Repository, we retrieved a total of 17 directed nonbipartite networks among the largest available in each technological or economic application domain to compensate for the relatively few networks of those types available from the other three data sources.

To eliminate redundancy among these networks, we identified and excluded duplicates by comparing the description, number of nodes, number of links, and other network statistics. Since the resulting set still contained a disproportionately large number of wiki-talk and wiki-link networks (28 and 167, respectively), we removed this bias by excluding all but the largest network from each of these two groups. This led to our final selection of 251 networks, consisting of 62 biological networks, 51 informational networks, 79 social networks, 40 technological networks, and 19 economic/game networks. The network size ranges from n = 10 to n = 7.9 × 106, and the link density ranges from 4.0 × 10−7 to 1.0. The adjacency matrix A of each network was constructed as follows. For each link, the corresponding A was set to the weight from the data if available and set to one otherwise. Multiple links between the same pair of nodes were combined into a single link whose weight equals the sum of the weights of the original links. For a self-link, the corresponding A was set based on the information from the data if available and A = 0 otherwise. The constructed adjacency matrices were then used to compute the nonnormality measure ∥D∥F, reactivity measure λΔ(A), degree imbalance 〈(din − dout)2〉, eigenvector angle θ1, and other properties of the networks. Of these 251 networks, 38 had degenerate or nearly degenerate λ1(A), which we numerically identified using the criterion θ1 > 1.57 [recalling that θ1 ≤ π/2 by definition and that θ1 = π/2 when λ1(A) is degenerate].

The minimum level of nonnormality observed among the 251 networks in Fig. 1A was for a 8192-node network of interactions in a square dielectric waveguide, with . The only nonreactive network in the dataset was one representing child-parent relationships from an online genealogical website called WikiTree, for which λΔ(A) and θ1 are both estimated to be zero to machine precision. Among the reactive networks that are nondegenerate, shown in Fig. 1B, the minimum value of λΔ(A)/λ1(A) ≈ 3.1 × 10−6 was observed for the 3140-node network of intercounty migration in the United States.

Directed random network model

To sample a network realization of size n from the GCLV model, we first draw the expected in-degree and out-degree of different nodes i = 1, …, n independently from a common joint probability distribution that does not depend on the network size n and has finite second moments. For each i, the random variables and are not restricted to be integers and need not to be independent (and thus can be correlated). Following (), the actual in- and out-degrees and of the network are then determined as a result of randomly creating a directed link from node j to node i with the probabilityfor each i and j. We then independently draw random weights A > 0 for the resulting links from a given distribution, while we set A = 0 if node j is not connected to node i. Note that this process can create self-links with random weights, but we could choose to prohibit them by setting A = 0 for all i. The unweighted version of this model is obtained if we instead set A = 1 for all links.

For the (marginal) distributions of and , we assume that ρ ≤ c for some constant and that the extreme values and asymptotically follow the so-called generalized extreme value (GEV) distributions () after appropriate normalization (see below for more details). We note that the first assumption is only a slight addition to the second, as the second implies that ρ → 0 in probability as n → ∞ (see Supplementary Materials, section S2.1, for a proof). We also assume that , so that we have almost surely (which can be shown using the strong law of large numbers). Denoting the mean of the expected in- or out-degree by (which does not depend on i), we assume . This holds true because unless the joint distribution of and is singular with all probability density concentrated at . These assumptions guarantee that, with probability one, the means of the actual in-degrees and out-degrees over the Bernoulli distribution of A indeed match and , respectively, in the limit of n → ∞. Thus, the given joint distribution of and can be interpreted as the expected degree distribution for this random network model.

For , the GEV assumption mentioned above is more precisely described as follows: there exist sequences of normalization constants a > 0 and b such thatand . Here, the GEV cumulative distribution function Gγ(x) is given byThe parameter γ is called the extreme value index and is determined solely by the distribution of . For , we assume the same but possibly with different a, b, and γ. The GEV distributions are the only ones that can arise as the limit distribution of the maximum of independent and identically distributed random variables if the limit in Eq. 9 exists, and they include the well-known Gumbel, Fréchet, and Weibull distribution as special cases.

The GEV assumption is mild and is satisfied for almost all commonly encountered distributions. It can be verified using the explicit conditions for Eq. 9 given in (). For example, one of these conditions (theorem 1.1.8) shows that it is satisfied by any power-law distribution with exponent β > 2, minimum expected degree , and the density function given by p(x) = Cx−β if , and p(x) = 0 otherwise, where is the normalization constant. The extreme value index in this case is γ = 1/(β − 1), and the sequences of constants in Eq. 9 are and . Since β > 2 (which is also the condition for this distribution to have a finite mean), we have γ < 1, implying and .

Nonnormality approximation for unweighted networks

Our derivation of Eq. 4 is for large random networks generated by the unweighted version of the GCLV model. Since Eq. 4 expresses as a function of the in-degrees and out-degrees , the derivation is based on calculating the conditional expected value , which we denote simply by in this section (in which we also use similar simplified notations for other expected values, variances, and probabilities). Taking the (conditional) expected value of Eq. 3, we obtain where we define (recalling the assumption A = 0, which will be relaxed in the next section). We note that, while and are held fixed, and are random.

We now seek to estimate in Eq. 11 by calculating and Var(Y), which can be broken down to individual terms and . We first note that . To compute this probability, we defineto represent the probability that the in-degree of node i is d and the probability that the out-degree of node j is d, respectively. Assuming the random variable fixed for the moment, we havewhere ∑α denotes the summation over all α = {α1, …, α} ∈ {0,1} in which exactly d elements are equal to one, and k: α = 0, 1 denotes all k such that α = 0, 1, respectively. Equation 13 also used the following approximations for large network size n (with d fixed): 1 − ρ ≈ e−ρ valid when ρ is small (which is the case because ρ → 0 in probability as n → ∞; see Supplementary Materials, section S2.1); and (based on the strong law of large numbers); and [a special case of the theorem proved in ()]. Analogously, holding fixed temporarily, we haveand similar arguments also show that

We can now calculate and from the definition of the conditional expected values asWe note that these expressions do not depend on nor , and hence, the estimates are valid even when these variables are allowed to be random. We thus havewhere we used that the high-order terms and are negligible for large n. Summing these over k (including the terms k = i, j, which are small compared to other terms and hence do not affect the sums for large n), we have

The term in Eq. 11 can now be estimated for large n assince is negligible compared to Var(Y), and we have(retaining the i = j term, which is small compared to the sum of the other terms). Substituting this into Eq. 11 yields Eq. 4.

Nonnormality approximation for weighted networks

Here, we consider the weighted version of the GCLV model, assuming that the distribution of link weights squared has a finite mean and a finite variance Var(w2). Assuming first that there are no self-links, i.e., A = 0, we will show that the nonnormality ∥D∥F can be approximated for large n asThis extends Eq. 4 to weighted networks. When the link weight variance Var(w2) is small, this equation yields the scaling . As in the case of unweighted networks, the derivation is based on taking the conditional expectation given fixed values of and . The weighted version of Eq. 11 reads:where Y = ∑(A − A) and thus . Following the derivation of Eq. 17 while accounting for the random weights, we haveThis leads to the weighted version of Eq. 21:To approximate the first term on the right side of Eq. 23, we note thatand that there are exactly and nonzero terms in and , respectively. Thus, given fixed values of and , the expectation and variance of can be estimated asand hence we haveCombining Eqs. 23, 26, and 30 and approximating by its conditional expectation yield Eq. 22.

The above results can be further extended to allow for any given assignment of self-links. Assume that with diagonal elements β, i = 1,2, ⋯, n. Then, Y = ∑(A − A) − (A − A)(β − β) and Y remains unchanged. Note thatTherefore, we have

Extension of the nonnormality approximation to Laplacian-coupled networks

Here, we consider the nonnormality of the Laplacian matrices of networks generated by the weighted GCLV model under the same conditions as in the derivation of Eq. 22. The Laplacian matrix of a network with adjacency matrix A is defined by L ≔ K − A, where K denotes the diagonal matrix with K = ∑. A straightforward calculation yieldswhere and, in particular, . Following Eq. 30, we haveIn addition, following Eqs. 20 and 26, we obtainTherefore, we have

Reactivity approximation for networks with dominant largest eigenvalue

Here, we derive Eq. 5 assuming λ1(H) ≈ λ1(H1) and the nondegeneracy of λ1(A). We normalize the left and right eigenvectors associated with λ1(A), so that we have and , and thus, the angle θ1 between the left and right eigenvectors is given by . Noting that the symmetric part of A1 can be written as and approximating the eigenvector associated with its largest eigenvalue by a linear combination αu1 + βv1, we havewhere we definedThe last maximum in Eq. 37 can be computed as the largest generalized eigenvalue μ of the matrix pencil (Φ(θ1), Ψ(θ1)) and satisfies the equationThis equation can be explicitly solved to yield . Substituting this into Eq. 37 and using the definition of λΔ, we obtain Eq. 5 as

Random networks used in Fig. 3

To generate these networks, we used the GCLV model with a given correlated identical distributions of the expected in-degree and the expected out-degree for each node i. To realize such a joint distribution, we first drew and (independently for each node i) from a common distribution, which was either (i) the gamma distribution, with probability density function p(x) ∼ (x − dmin)e−, x > dmin, where dmin is the minimum degree, a > 0 is the shape parameter, and b > 0 is the rate parameter; or (ii) a discrete power-law distribution with the probability mass function p(x) ∼ x−β for integers x = dmin, dmin + 1, …, dmax, where β > 0 is the scaling exponent, dmin is the minimum degree, and is the maximum degree imposed to ensure that the model assumption ρ ≤ c is satisfied (and here, we set c = 0.99). For the gamma distribution, the mean degree d and the parameter 1/b were drawn randomly from the intervals [10,50] and [2,10], respectively, and the parameter a was then set to be a = b(d − dmin) (to ensure that the mean degree equals d). For the power-law distribution, the scaling exponent β was drawn randomly from the interval [2,4]. For both distributions, we used the minimum degree dmin = 10.

After generating and , correlation was added between them using a parameter ρ randomly chosen from the interval [−1,1] to increase the range of ∥D∥F and λΔ(A) observed. The correlation was created by sorting and for a randomly chosen subset of nodes (with mean fraction ∣ρ∣) in the same order for and if ρ > 0 (leading to positive correlation) and in the opposite order if ρ < 0 (leading to negative correlation).

For the weighted networks, the random weights were drawn either from the exponential distribution with mean μ and variance μ2 (for μ = 0.7, 1.0, or 1.3) or from the positive domain of the Gaussian distribution with mean one and variance σ2 (for σ = 0.1 or 0.2). For this figure, we allowed self-links in the GCLV model (with random weights drawn from the same distribution as the other links).

Randomization of real networks in Fig. 4A

For a given network from the dataset, we generated its randomization using the GCLV model. The in- and out-degrees of node i in the network were used as the expected degrees and , respectively, to generate a random network topology. Note that, for nodes i and j with the right side of Eq. 8 exceeding one, we set the connection probability ρ = 1. The link weights of the network were resampled (with replacement) to generate random weights.

Determination of the shaded region in Fig. 4C

First, for each β, we computed the averages of λ1(H) and λ1(A) over 10 realizations of the GCLV model for 15 values of n, equally spaced on the logarithmic scale between n = 102 and n = 105. We then performed a least-squares fit of these averages as functions of n on the logarithmic scale for n ≥ 300 and used the resulting slopes as the scaling exponents for λ1(H) and λ1(A) for the given β (shown in the background by blue and green dots, respectively). The boundary curves were obtained by fitting these scaling exponents with seventh-order polynomials. The shaded area between these two curves thus represents the parameter region in which the scaling exponent ℓ for α(n) is smaller than that for λ1(H) but larger than that for λ1(A), implying that, as n increases, λ1(H) − α(n) increases while λ1(A) − α(n) decreases.