## Abstract

We analyze a collection of empirical networks in a wide spectrum of disciplines and show that strong non-normality is ubiquitous in network science. Dynamical processes evolving on non-normal networks exhibit a peculiar behavior, as initial small disturbances may undergo a transient phase and be strongly amplified in linearly stable systems. In addition, eigenvalues may become extremely sensible to noise and have a diminished physical meaning. We identify structural properties of networks that are associated with non-normality and propose simple models to generate networks with a tunable level of non-normality. We also show the potential use of a variety of metrics capturing different aspects of non-normality and propose their potential use in the context of the stability of complex ecosystems.

## INTRODUCTION

Network science (*1*–*3*) has emerged, in the past 20 years, as an essential framework to model and understand complex systems in a variety of disciplines, including physics (*1*), economics (*4*), biology (*5*), and sociology (*6*). At its core, network science views a system as a set of nodes that may be connected directly by an edge or indirectly by a succession of edges, thereby forming paths of interactions. The bridge between network structure and dynamics is generally unraveled by defining a linear dynamical model on the nodes; take, for instance, a random walk process as a simple model of diffusion or the linearization around a critical point of a nonlinear dynamical system (*7*–*10*). In each case, the process is determined by a matrix, somehow related to the adjacency matrix of the underlying network. In addition, critical aspects of the system, such as its stability and characteristic time scales, are usually described by the properties of its spectrum (*11*). Central network concepts such as the spectral gap, spectral radius, and master stability conditions all build on this interpretation. Relatedly, network spectra also appear in network algorithms, such as in community detection (*12*) or in network comparison (*13*).

The characterization of a linear system by its spectrum is canonical, but it is unreliable in situations when the linear operator is non-normal; namely, its eigenvectors do not necessarily form an orthonormal basis, and the transformation to eigenvector coordinates may involve a strong distortion of the phase space. Non-normality has a long tradition in linear algebra and dynamical systems, from early studies in hydrodynamics (*14*) to more recent works on the robustness of non-normal ecosystems (*15*) and in neuronal dynamics (*16*, *17*). Yet, these results remain focused on limited areas of science, and a systematic study of the prevalence of non-normality in real-world networks, as well as its potential impact on dynamics, is still lacking. Here, we call non-normal a network whose adjacency matrix **A** is non-normal (*18*). By definition, **A** is non-normal if it verifies **AA**^{T} ≠ **A**^{T}**A**. It is thus clear that **A** needs to be asymmetric to be non-normal or, equivalently, the network needs to be directed to be non-normal, but, as we will discuss in more detail later, asymmetry is not sufficient and certain types of network architectures are necessary to determine a strong non-normality. Given a non-normal network, other standard matrices, such as its Laplacian **L**, are also non-normal. Non-normality can hence be quantified using a standard spectral measure borrowed from matrix theory, such as Henrici’s departure from normality (*19*), , where || ⋅ ||_{F} is the Frobenius norm and λ_{i} are the eigenvalues of the adjacency matrix. A zero value is associated with a symmetric network, while the larger the values, the stronger the non-normality.

Before going further, let us briefly illustrate in more detail the influence of non-normality on the prototypical example of a linear (linearized) dynamics on a non-normal network. Without loss of generality, we consider the model , where **M** encodes the linear dynamics on the network, through the dependence on **L**, and forms a stable matrix; namely, the spectral abscissa α(**M**) = max *ℜ*σ(**M**) is not positive, with σ(**M**) being the spectrum of the matrix **M**. In the case of normal networks, the solution of the linear system would consist of a linear combination of exponentially relaxing modes, each with a characteristic time scale given by the inverse of the corresponding eigenvalue; hence, the spectral abscissa is responsible for the long-term dynamics. In situations when the network is non-normal, however, more complex patterns may emerge. Standard measures of non-normality of the matrix and their relation to dynamics are provided in Table 1. If **M** has a positive numerical abscissa, ω(**M**) = maxσ(*H*(**M**)), where *H*(**M**) = (**M** + **M**^{T})/2 is the Hermitian part of **M**, then the system can undergo a transient growth before asymptotically converging to zero, as measured by the norm of the state vector **x** (see Table 1A). This transient behavior cannot be explained by the picture provided by the spectrum of the matrix **M** (*19*) and can have a strong impact once nonlinearities are at play. In situations when the dynamics is obtained from the linearization around a critical point, this initial growth may trigger nonlinear terms, take the system far away from the equilibrium, and thus radically reshape the dynamical behavior of nonlinear systems (*18*), as shown in Table 1B.

Although a system can initially be close to an asymptotically stable equilibrium, it can leave this state even when a moderate external perturbation occurs due to non-normality (*19*). This effect is even more notable once one includes stochastic forces to the model, e.g., exogenous or demographic perturbations due to the surrounding environment (*20*), as they may push even a linear, stable system out of equilibrium when it is non-normal. Again, this behavior cannot be properly captured by the spectrum of the linear model, nor can it be described by the numerical abscissa, which only determines the short-term behavior of the dynamics. To describe the long-term consequences of these perturbations, another important tool is the pseudospectrum σ_{ϵ}(**M**) = {σ(**M** + **E**):∥**E**∥≤ϵ}, for ϵ > 0 (*19*), from which one can compute the ϵ–pseudospectral abscissa, replacing somehow the role of the spectral abscissa and eventually the Kreiss constant K, which provides a direct measure of the size of the transient amplification (see Table 1). By definition, a complex number *z* is an eigenvalue of **M** if a bounded inverse of *z***I** − **M** does not exist. The pseudospectrum is based on a less strict definition and defines regions of the complex plane where ∥(*z***I** − **M**)^{−1}∥ is larger than a prescribed positive number ϵ^{−1}. By its very first definition, the pseudospectrum defines regions of the complex plane where eigenvalues of a matrix can be found because of a small perturbation, **M** + Δ**M**, with ∥Δ**M**∥ < ϵ. These perturbations lead to small variations of the spectrum in the case of normal matrices, but they can become much more important in the case of non-normal matrices. In particular, even small perturbations can make a linearly stable system unstable. Note that this effect may have important practical consequences for networks, as the precise value of edge weights is often unknown (*21*), and empirical measurements of networks are prone to missing edges (*22*).

As we have discussed, non-normality may strongly affect linear and nonlinear dynamical systems on networks and, more generally, their behavior. The contributions of this work are manifold. First, we show that a strong non-normality is widespread in complex networks empirically observed in a variety of domains. As a second step, we reveal the organization behind non-normality and show that non-normality is associated with a combination of absence of cycles (*23*), low reciprocity (*24*), and hierarchical organization (*25*). We also propose a simple model for growing networks based on preferential attachment reproducing our observations. Last, we consider in detail a Lotka-Volterra model applied to a real-world network and show that the use of network metrics for non-normality helps to understand the dynamics of the system.

## RESULTS

### Non-normal networks: Empirical data and the shape of non-normality

As a first step, we have considered a large set of directed, real-world networks from different disciplines, including biology, sociology, communication, transport, and many more. Results reported in Table 2 (see also the more complete table presented in the Supplementary Materials) show values of standard measures of non-normality, including the numerical abscissa, the ϵ–pseudospectral abscissa, and the normalized Henrici’s departure from normality, , all revealing that the networks present a strong non-normality.

As a next step, we investigate the type of network organization associated with non-normality. The directedness and low reciprocity of a network are necessary conditions for non-normality, but they are by no means sufficient. For instance, a *k*-regular directed ring, whose adjacency matrix is circulant, is normal because of its rotational symmetry (*18*). The condition **AA**^{T} ≠ **A**^{T}**A** is instead satisfied when the network is hierarchical, that is, when nodes have a rank and edges with a strong weight tend to flow from nodes with a small rank to nodes with a high rank (or vice versa). These organizations are known to be prevalent in different types of networks (*25*–*28*), for instance, through the concepts of dominance hierarchies in social ecology, trophic levels in food webs, and social status in social networks. The inequality becomes maximum when the network is a directed acyclic graph (DAG), such that the matrix takes an upper triangular form after proper relabeling of the nodes. Again, DAGs find several applications, for instance, in the case of citation or causal networks. On the basis of these intuitions, we estimate the level of hierarchy of a real-world network as follows. Given an adjacency matrix, we search for the best nodes ordering such that the total weight in the upper triangular part of the matrix is maximal (see Methods), thus allowing one to identify the DAG of maximum weight embedded in the network. A simple measure of asymmetry is then the unbalance Δ between the number of entries in the upper and lower triangular part of the relabeled adjacency matrix **Ã**, that is, Δ := |*K*^{<} − *K*^{>}|/*K*, where *K*^{<} = ∑_{i <j}*Ã*_{ij}, *K*^{>} = ∑_{i >j}*Ã*_{ij}, and *K* = *K*^{<} + *K*^{>}. Δ thus measures the asymmetry of the adjacency matrix after relabeling, and it provides a structural indicator of non-normality that can be compared with standard spectral measures, such as the normalized Henrici’s departure from normality. In Fig. 1B, we show that the normalized Henrici’s index strongly correlates with Δ across a variety of real-world networks, hence reinforcing the connection between structural hierarchy and dynamical non-normality (see also Methods). Last, note that non-normality and a hierarchical structure are global properties of a network. As an illustration, we have considered in the Supplementary Materials the case of networks built from different combinations of the same constituting blocks or motifs (*25*, *29*–*31*), and we observe that different levels of non-normality can emerge from the combination of two directed motifs, from strong non-normality and a DAG structure for the whole graph to weakly non-normal patterns, where the presence of a cycle prevents the dynamical flow to accumulate on a small number of nodes.

### Mechanistic model

Here, we propose a simple mechanistic model, denoted by the nSF network, leading to the formation of networks with tunable levels of non-normality. The model builds on the seminal ideas of de Solla Price (*32*), later leading to the family of preferential attachment models (*33*). As discussed before, critical ingredients of the model should be the low reciprocity of the directed edges and the presence of a hierarchical structure. We thus consider a growing network where, at each step, a new node *j* draws a directed edge to a previously existing node *i*, with a probability proportional to its in-degree. Note that this type of process is expected to lead to the formation of power-law in-degree distributions, but this is not our concern here. Node *i* also has the possibility of reciprocating and creating a link directed at *j*, as in (*34*). Asymmetry can be included in different ways, either by endowing edges with a weight and imposing that 0 ≤ *w*_{ji} ≪ *w*_{ij} or by considering unweighted edges and assuming that the reciprocal edge is created with a probability *p*_{i → j} ≪ 1, see Fig. 2. Hierarchy is then induced by the ordering of the nodes in terms of their arrival time. As expected, the stronger the inequalities, the stronger the non-normality of the resulting networks. We have investigated the non-normality of the resulting networks and found a notable similarity with the relation between Δ and Henrici’s departure from normality observed in real-world networks, as shown in Fig. 1. Note that we have also considered variants of other classical network models, such as Erdös-Rényi (ER) (*35*) and Watts-Strogatz (WS) (*36*) (see the Supplementary Materials), but their lack of hierarchical structure prevents the formation of strong non-normality, as can be seen in their pseudospectral properties (see Fig. 3 and the Supplementary Materials).

### Application to the stability of complex ecosystems

The hierarchical structure of non-normal networks allows for the introduction of interesting connections with dynamical systems. Here, we focus on stability, a central concept to understand the emergence of collective phenomena (*10*). The importance of network structure for stability is well established since the seminal works of May (*7*) and Allesina and Tang (*8*) in the context of ecology. For instance, choosing the interaction strengths from a normal distribution , May proved that an ecosystem loses its stability above a critical size, as a consequence of the circular law (*37*). To understand the interplay of non-normality and dynamics, we analyze the master stability function (*3*, *10*), which is a general tool that allows one to infer about the (in)stability of a networked dynamical system; it often relies on the use of the spectrum of some suitable linearization, while hereby conditions for (in)stability are determined through the pseudospectrum of the linearized system (*1*). The latter represents the generalized Lotka-Volterra (GLV) model, which is popular for understanding competition and mutualism among interacting species (*7*–*9*). The set of equations governing the dynamics of trophic interactions is given by (*9*)(1)

Here, *r*_{i} are the intrinsic rates of (i) birth if *r*_{i} > 0, meaning that species *i* can reproduce itself in absence of other species and in abundance of resources; (ii) death if *r*_{i} < 0 in the sense that the population of species *i* will decline in absence of other species (e.g., preys). The positive constants *s*_{i} represent the finite carrying capacity of the ecosystem (limited resources) and prevent the species *i* from growing indefinitely. An important role is played by the community matrix **M**, whose entries *M*_{ij} (respectively *M*_{ji}) represent the influence of species *j* on *i* (respectively *i* on *j*). We also assume that *M*_{ii} = 0 , ∀ *i*; namely, the community matrix describes only interspecies interactions, and intraspecies interactions have been cast into *s*_{i}. In the following, we adopt the method of Chen and Cohen (*38*), as already explained in the literature (*8*, *9*). More precisely, we hypothesize the existence of a positive equilibrium solution **x*** that, without loss of generality, can be assumed to be of the form , for all *i*, after a suitable choice of the growth/death rates *r*_{i}. At this point, the master stability function of the GLV model depends solely on the spectrum (pseudospectrum) of the matrix **M** − diag(s), namely, the community matrix from which we remove the matrix whose diagonal contains the interspecies strengths *s*_{i}. The problem is hence mapped to a framework where stability directly depends on species interactions.

The vulnerability of the system is visible in Fig. 4, where the spectra (black dots) shift from the left to the right of the imaginary axis once mutualism increases. Although the system remains stable for a strong competitive setting (Fig. 4, A and B), we can observe (colored curves) that the ϵ–pseudospectral abscissa [see (*19*) and Table 1] is positive and larger for structured systems than for random ones (*8*). To represent the former systems, we used the nSF networks obtained with the generation model previously presented, the latter exhibiting features very similar to the real trophic relations. This implies that the system can easily be destabilized by (relatively) small fluctuations due to demographic, thermal, or endogenous noise that are always present in the surrounding environment and are amplified because of the non-normality (see the Supplementary Materials). This remark can have important consequences in the understanding of the problem of coexistence of multiple species in a harsh competitive environment, e.g., in the case of the paradox of the plankton (*39*), for which field observations are at odds with the principle of competitive exclusion (*8*).

## DISCUSSION AND CONCLUSIONS

We have shown that a large number of real-life networks are strongly non-normal and that a characterization of their properties solely by spectral methods may be misleading. Non-normality induces a strong dependency on fluctuations and needs to be considered with care when performing a linear stability analysis of nonlinear systems. Despite the fact that the non-normality is well studied and that its importance has been recognized in a variety of domains, a systematic analysis of its importance and effect in large-scale networks was still lacking. Our first contribution is thus not only the identification of what appears to be a ubiquitous property of directed networks but also the introduction of new methods in the toolbox of network science to generate non-normal networks and capture the effect of non-normality on their dynamics. Potential applications have recently been explored, for instance, in pattern formation on networks (*18*) and in epidemic spreading in metapopulation models (*18*). Overall, these findings emphasize that non-normality is a critical component of complex systems and that specific tools are necessary to complement standard methods based on eigenspectra, which are prevalent in network science. More specifically, this new perspective may shed light on how to explain the diversity of species in ecosystems (*40*), the origin of cascade failures in power grids (*41*), or the spread of epidemics in mobility networks (*42*), just to mention a few possible applications.

## METHODS

### Measures of non-normality

A real matrix **M** is said to be non-normal if it is not diagonalizable by a unitary matrix; namely, its eigenvectors are not orthogonal to each other (*19*). The numerical abscissa has been introduced in population dynamics (*15*) with the term of reactivity, and it is defined by ω(**M**) = supσ(*H*(**M**)), where *H*(**M**) = (**M** + **M**^{T})/2 is the Hermitian part of **M**. This is a very natural concept; however, it does not allow the computation of the maximum amplification of the initial conditions exhibited by linear stable non-normal systems. For this reason, one has to resort to the pseudospectrum (*19*), σ_{ϵ}(**M**), which is defined for all ϵ > 0 as the spectrum of the perturbed matrix **M** + **E**, for any matrix ∥**E**∥ ≤ ϵ. From the ϵ–pseudospectral abscissa, α_{ϵ}(**M**) = supℜσ_{ϵ}(**M**), we can obtain the Kreiss constant, = sup_{ϵ>0}α_{ϵ}(**M**)/ϵ, and eventually the lower bound on the orbit size(2)

Let us observe that the latter provides a straightforward bound on the amplification envelope defined in (*15*). Moreover, is more informative that reactivity, in fact a stable system, can exhibit a small amplification even if ω(**M**) > 0 is very large.

Henrici’s index is based on the observation that the Frobenius norm of a normal matrix is given by , where λ_{i} are the eigenvalues of the matrix; one can thus define Henrici’s departure from normality (*19*) for a non-normal matrix **M** by(3)

It attains its minimum once the matrix is normal and then increases as long as the matrix deviates from normality. To compare systems with different sizes, we define the normalized index .

For a generic matrix with binary (respectively positive) entries, one can define the imbalance between the number (respectively the total sum) of entries in the upper and lower triangular part, using the language of networks(4)where are the entries of the final relabeled matrix.

While it can be relatively easy to determine a DAG and compute Δ, once we have a drawing of a (small enough) network, this task becomes hard starting from the adjacency matrix or a large network. We observed that the simple operation of relabeling the nodes can change the value of Δ and that the latter increases the larger the number of entries in the adjacency matrix in the upper triangular part, namely, links *i* → *j*, where *j* > *i*. Having in mind these observations, we designed an algorithm aiming at maximizing Δ once couples of nodes are relabeled; i.e., rows and columns of the adjacency matrix are swapped. To overcome the combinatorial difficulty of the problem, we resorted to a simulated annealing method (*43*) to get an accurate solution in a relatively short time. A pseudocode is presented in the Supplementary Materials, and the generic convergence behavior of the maximization process is shown in fig. S2.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/12/eaau9403/DC1

Section S1. Non-normal matrices and their pseudospectra

Section S2. Global structure of non-normal networks

Section S3. Models for generation of non-normal networks

Section S4. An intuitive meaning of the pseudospectrum

Section S5. Pseudospectra of real non-normal networks

Section S6. Extended table of real non-normal networks

Fig. S1. Time evolution of the norm of the solution of the linear ordinary differential equation .

Fig. S2. Convergence of the maximization algorithm.

Fig. S3. Henrici’s departure from normality versus Δ.

Fig. S4. Motifs organization and non-normality.

Fig. S5. Spectra and pseudospectra of non-normal networks.

Fig. S6. Pseudospectra for real networks I.

Fig. S7. Pseudospectra for real networks II.

Table S1. Some figures for real webs.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**Gratitude goes to D. Fanelli and P. K. Maini for useful discussions.

**Funding:**The work of M.A. and T.C. presents research results of the Belgian Network DYSCO, funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The work of M.A. was also supported by an FRS-FNRS Postdoctoral Fellowship.

**Author contributions:**M.A. and T.C. designed the study, conducted the formal analysis, and analyzed the data. All authors organized and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).