Universal transient behavior in large dynamical systems on networks

TL;DR

This paper studies how the transient response of large, randomly coupled linear dynamical systems depends on network topology, revealing universal behaviors and deriving analytical expressions for different network types.

Contribution

It develops a formalism for analyzing transient dynamics on locally tree-like graphs and derives universal analytical results for unidirectional networks with various disorders.

Findings

Unidirectional networks show universal transient behavior depending on a single parameter.

Analytical expressions for average transient norms are derived for different disorder types.

Numerical experiments confirm the theoretical predictions on real-world networks.

Abstract

We analyze how the transient dynamics of large dynamical systems in the vicinity of a stationary point, modeled by a set of randomly coupled linear differential equations, depends on the network topology. We characterize the transient response of a system through the evolution in time of the squared norm of the state vector, which is averaged over different realizations of the initial perturbation. We develop a mathematical formalism that computes this quantity for graphs that are locally tree-like. We show that for unidirectional networks the theory simplifies and general analytical results can be derived. For example, we derive analytical expressions for the average squared norm for random directed graphs with a prescribed degree distribution. These analytical results reveal that unidirectional systems exhibit a high degree of universality in the sense that the average squared norm only depends on a single parameter encoding the average interaction strength between the individual constituents. In addition, we derive analytical expressions for the average squared norm for unidirectional systems with fixed diagonal disorder and with bimodal diagonal disorder. We illustrate these results with numerical experiments on large random graphs and on real-world networks.

Figures

Equations

$${\partial }_t{\zeta }_i=f_i(\mathbf{\zeta }),$$

$$\frac{\mathrm{d}{y}_i(t)}{\mathrm{d}t}=\sum _{k=1}^N{A}_{ik}{y}_k(t),$$

$$\mathbf{y}(t)=\mathbf{\zeta }(t)-{\mathbf{\zeta }}^{\star }$$

$$A_{ik}=\left(\frac{\partial f_{i}}{\partial\zeta_{k}}\right)\Big{|}_{\bm{\zeta}^{\star}}$$

$$S_N(t;A)=⟨|\mathbf{y}(t){|}^2⟩$$

$$S_N(t)=\overline{⟨|\mathbf{y}(t){|}^2⟩}.$$

$$S(t)=\underset{N\to \infty }{\lim }{S}_N(t).$$

$$\underset{t\to \infty }{\lim }S(t)=0(\underset{t\to \infty }{\lim }S(t)=\infty ),$$

$$t^{\star }(N)={\min }_{t\ge 0}{S}_N(t).$$

$$A_{ij}=-\mu {\delta }_{ij}+X_{ij}/\sqrt{N},$$

$$S(t)=e^{-2\mu t}{I}_0(2\rho t)\stackrel{t\to \infty }{\sim }\frac{e^{2t(\rho -\mu )}}{\sqrt{4\pi \rho t}},$$

$$A_{ij}=D_i{\delta }_{ij}+C_{ij}{J}_{ij},$$

$$p_{\mathrm{deg}}(k_{\mathrm{in}},k_{\mathrm{out}})=p_{\mathrm{deg}}(k_{\mathrm{in}})p_{\mathrm{deg}}(k_{\mathrm{out}}).$$

$$S_N(t;A)={\text{varoiint}}_{\gamma }\frac{\mathrm{d}z\mathrm{d}w{e}^{t(z+w)}}{N(2\pi \mathrm{i}{)}^2}{\mathcal{W}}_A(z,w),$$

$${\mathcal{W}}_A(z,w)=\mathrm{Tr}\left[\frac{1}{z1-A^T}\frac{1}{w1-A}\right]$$

$$S(t)={\text{varoiint}}_{\gamma }\frac{\mathrm{d}z\mathrm{d}w{e}^{t(z+w)}}{(2\pi \mathrm{i}{)}^2}\overline{\mathcal{W}}(z,w),$$

$$\overline{\mathcal{W}}(z,w)=\underset{N\to \infty }{\lim }\frac{1}{N}\overline{\mathrm{Tr}\left[(z1-A^T{)}^{-1}(w1-A{)}^{-1}\right]}$$

$$\mathbf{y}(t)=e^{At}\mathbf{y}(0),$$

$$|\mathbf{y}(t){|}^2={\mathbf{y}}^T(0)e^{A^{T}t}{e}^{At}\mathbf{y}(0).$$

$$\mathbf{y}(0)=O{\mathbf{e}}_1,$$

$${\mathbf{e}}_1=(1,0,0,\dots ,0{)}^T.$$

$${\mathbf{e}}_1^T{O}^T{e}^{A^{T}t}{e}^{At}O{\mathbf{e}}_1$$

$$⟨O_{ij}{O}_{kl}⟩=\frac{1}{N}{\delta }_{ik}{\delta }_{jl},$$

$$S_N(t;A)=\frac{1}{N}\mathrm{Tr}{e}^{A^{T}t}{e}^{At}.$$

$$f(A)=\frac{1}{(2\pi \mathrm{i})}{\oint }_{\gamma }\frac{f(z)}{z-A}\mathrm{d}z,$$

$$A=\sum _{j=1}^N{\lambda }_j|r_j⟩⟨{\ell }_j|,$$

$$S_N(t;A)=\frac{1}{N}\sum _{j,k}{\tilde{O}}_{jk}{e}^{t({\lambda }_j^{\star }+{\lambda }_k)},$$

$${\tilde{O}}_{jk}=⟨{\ell }_k|{\ell }_j⟩⟨r_j|r_k⟩$$

$$⟨|\mathbf{y}(t){|}^2⟩=\frac{1}{N}\sum _{j=1}^N{e}^{2t\mathrm{Re}[{\lambda }_j]}.$$

$${\mathcal{W}}_A(z)=\mathrm{Tr}\frac{1}{z1-A^T}$$

Full text

I Introduction

Networks of interacting constituents appear in the study of systems as diverse as ecosystems [1, 2, 3], neural networks [4, 6, 7, 5, 8], financial markets [9, 10, 11, 12], and signaling networks [13, 14, 15]; for more examples see Refs. [16, 17]. Traditionally, a strong focus has been put on whether such systems are stable at long time scales [18, 19] because stability is often associated to functionality, e.g., stable ecosystems or economies [20]. Differently, the short-time transient response of networked systems to perturbations is less understood despite being important for applications: for example, neuroscientists administer magnetic stimulations to the brain and observe distinct dynamical responses of electrical activity, which capture different connectivity states of the underlying network of neurons [21]; in ecological systems, the asymptotic dynamics does not capture the typical time scales accessible in experiments [22, 23, 24, 25, 26]; in the context of epidemics, the initial time window before vaccinations become available, makes a crucial difference in limiting the extent of the outbreak [27]. A relevant question is thus how network topology determines the early-time dynamics of large systems. For example: (i) how long does a stable system take to return to its stationary state as a function of the network topology and interaction strength among its constituents, and (ii) how long does it take to realize that a seemingly stable system is unstable after all, and disaster is looming?

We describe the state of a large dynamical system at time $t$ with $N$ real-valued variables $\zeta_{i}(t)$. For example, $\zeta_{i}(t)$ may represent the abundance of species $i$ in an ecosystem or the activity of the $i$-th neuron in the brain at time $t$. We assume that the system evolves according to a system of first-order equations

[TABLE]

for $i=1,\ldots,N$. Although the functions $f_{i}$ are arbitrary, the dynamics can be linearized close to a stationary point $\bm{\zeta}^{\star}$ for which $f_{i}(\bm{\zeta}^{\star})=0$ to yield [28, 29]

[TABLE]

where the $N$-dimensional vector

[TABLE]

encodes deviations from the stationary state and

[TABLE]

are the entries of the corresponding Jacobian matrix $A$.

The magnitude of the deviation vector $\bm{y}(t)$ as a function of time is captured by the squared norm $|\bm{y}(t)|^{2}$. In order to grasp the dynamical response to generic initial perturbations, we consider the average

[TABLE]

over initial states $\bm{y}(0)$ uniformly drawn from the sphere $|\bm{y}(0)|=\alpha$, where $\alpha$ quantifies the strength of the initial kick; without loss of generality we set $\alpha=1$. In addition, to capture the properties of a typical dynamical system, we take $A$ to be a random matrix of pairwise interactions [18, 3, 30], and we further average the squared norm over the disorder

[TABLE]

Since we will be interested in large systems, we take the limit

[TABLE]

If

[TABLE]

then we say that a system is transiently (un)stable. Note that transient stability is different from asymptotic stability, which is governed by the sign of the real part of the leading eigenvalue of $A$, see e.g. Refs. [18, 3, 31].

The observables $S_{N}(t)$ and $S(t)$ are the main objects of interest. The quantity $S(t)$ describes the transient dynamics of an infinitely large dynamical system and can be computed analytically for a large and important class of systems, as we show in this paper. On the other hand, $S_{N}(t)$ captures the dynamics of systems at finite $N$ and is the quantity that is experimentally measurable. There exists a crossover time $t^{\star}(N)$ such that for times $t<t^{\star}(N)$ theory and experiment are in correspondence, i.e. $S_{N}(t)\approx S(t)$, while for $t\gg t^{\star}(N)$ theory and experiment are in disagreement, i.e. $S_{N}(t)\gg S(t)$. The discrepancy between $S(t)$ and $S_{N}(t)$ at large $t$ is most evident for transiently stable systems, for which $\lim_{t\rightarrow\infty}S_{N}(t)=\infty$ but $\lim_{t\rightarrow\infty}S(t)=0$.

For transiently stable systems, where $S_{N}(t)$ has a non-monotonic behavior, we define the crossover time by

[TABLE]

The crossover time $t^{\star}(N)$ defined by Eq. (9) increases with $N$ and diverges for large $N$ [seemingly as a power law, see Appendix A for a detailed discussion]. For the reasons above, $S(t)$ is a good measure of the transient dynamics of large dynamical systems.

So far, the quantity $S(t)$ has only been computed for systems with fully connected topology [36, 32, 34, 33, 35], namely for

[TABLE]

where the $X_{ij}$ are independent and identically distributed (i.i.d.) entries with zero mean and finite moments. In this case,

[TABLE]

where $\rho$ is the spectral radius of the matrix $X/\sqrt{N}$ and $I_{0}(x)$ is the modified Bessel function of the first kind. Since $S(t)$ only depends on the spectral radius of the matrix $X$, it enjoys a high degree of universality.

We aim to study how $S(t)$ depends on the topology of a network, for which it is necessary to go beyond the fully connected paradigm. We assume that $A$ is the adjacency matrix of a weighted graph that is locally tree-like,

[TABLE]

where $D_{i}\in\mathbb{R}^{-}$ are the diagonal decay rates, the $J_{ij}\in\mathbb{R}$ are the coupling strengths, and the $C_{ij}\in\left\{0,1\right\}$ are the entries of the adjacency matrix of a locally tree-like, directed and simple graph. We say that a sequence of graphs is locally tree-like if in the limit $N\rightarrow\infty$ the finite neighborhood of a randomly selected node is with probability one a tree [37, 38]; loosely speaking, a graph is locally tree-like if it is large and does not contain small cycles. In this paper, we develop a mathematical method to compute $S(t)$ for the model given by Eq. (12) under the sole assumption that the graph represented by $C$ is locally tree-like.

We further illustrate the general mathematical formalism on a canonical class of random directed graphs, namely, the ensemble of adjacency matrices of weighted random graphs with a prescribed degree distribution [39, 40, 17, 41, 31]. These random graphs are used to model real-world systems, such as, the Internet [64, 39, 65], neural networks [66] and other high-dimensional systems [40, 17, 16]. In this ensemble, the matrix $A$ in Eq. (12) is defined as follows:

• •

the $D_{i}$ are i.i.d. taken from a probability density $p_{D}(x)$;

• •

the $J_{ij}$ are i.i.d. random variables with probability density $p_{J}(x)$;

• •

the $C_{ij}$ are the entries of the adjacency matrix of a random (directed) graph with a prescribed joint degree distribution $p_{\mathrm{deg}}(k_{\mathrm{in}},k_{\mathrm{out}})$ of in-degrees $k_{\mathrm{in}}$ and out-degrees $k_{\mathrm{out}}$. In addition, we assume that the degree distribution factorizes as

[TABLE]

The spectral properties of this ensemble have been studied in detail in Refs. [42, 43, 41, 31] and the asymptotic stability of dynamical systems described by this ensemble has been studied in Ref. [31]. Here, we characterize the transient response to a random perturbation for dynamical systems on random graphs with a prescribed degree distribution by deriving analytical expressions for $S(t)$.

This analytical progress is made possible because (i) our general theory simplifies for locally oriented (or unidirectional) networks (see Section III.3 for a precise statement), and (ii) directed random graphs with a prescribed degree distribution are – with high probability – locally oriented [37, 38].

Remarkably, we find that $S(t)$ is universal for directed random graphs with a prescribed degree distribution, in the sense that it only depends on the distribution $p_{J}(x)$ and the degree distribution $p_{\mathrm{deg}}(k)$ through a single parameter (see our main formula (64)). On the other hand, the dependence on $p_{D}$ is non-universal.

We compare the derived analytical results for $S(t)$ with numerical results for $S_{N}(t)$ on random graphs and on real-world graphs. We find that $S(t)$ is in excellent agreement with $S_{N}(t)$ as long as $t<t^{\star}(N)$, where $t^{\star}(N)$ is a timescale that diverges with $N$. For $t\gg t^{\star}(N)$, it holds that $S_{N}(t;A)\sim e^{2{\rm Re}[\lambda_{1}(A)]t}$, with $\lambda_{1}(A)$ the eigenvalue of $A$ with the largest real part, and as a consequence, $S(t)\ll S_{N}(t)$.

The plan of the paper is as follows. In Sec. II, we express $S(t)$ as a contour integral over the two-point correlator of a random matrix ensemble, effectively mapping a dynamical systems problem onto a random matrix theory problem. In Sec. III, we develop a mathematical formalism to compute the two-point correlator, and thus also $S(t)$, on tree or locally tree-like graphs. In Sec. IV, we consider directed random graphs with a prescribed joint degree distribution and we derive for this ensemble analytical expressions for $S(t)$. In Sec. V, we compare the obtained analytical expressions for infinitely large graphs with numerical experiments on random graphs of finite size and on real-world graphs. In Sec. VI, we discuss the obtained results and present an outlook for future research. In Appendix A, we analyze how $t^{\star}(N)$ depends on $N$. In Appendix B, we compute a contour integral that appears in our calculations. In Appendix C, we make a study of the spectra of random graphs with a prescribed degree distribution and a bimodal distribution of diagonal matrix entries, and in Appendix D, we present numerical results for random graphs with power-law degree distributions.

II Mapping onto a random matrix problem

In this section, we derive the formula

[TABLE]

which expresses the dynamical response $S_{N}(t;A)$ as a contour integral of the two-point correlator

[TABLE]

over a closed counterclockwise-oriented contour $\gamma$ that encloses all eigenvalues of $A$. The symbol $\mathbb{1}$ denotes the identity matrix of size $N$.

Analogously, we obtain

[TABLE]

where

[TABLE]

is the average two-point correlator.

Formulae (14)-(17) provide a recipe to compute $S(t)$: if we obtain an expression for the two-point correlator (17) of a random matrix ensemble, then $S(t)$ follows readily from evaluating the contour integral (14). Note that we only need to compute $\overline{\mathcal{W}}(z,w)$ for values $z,w$ that lie outside the continuous part of the spectrum.

In order to obtain (14), we first express the solution of Eq. (2) as

[TABLE]

and therefore

[TABLE]

Since $|\bm{y}(0)|^{2}=1$, there exists a matrix $O$ in the the orthogonal group $\mathrm{O}(N)$ (the group of isometries of the $N$-sphere) such that

[TABLE]

where

[TABLE]

Taking the average of (19) with respect to all initial conditions $\bm{y}(0)$ selected uniformly at random on the unit sphere is thus equivalent to taking the average of the following expression

[TABLE]

with respect to the uniform (Haar) measure on the orthogonal group. Using that [44]

[TABLE]

and Eqs. (22) and (5), we obtain

[TABLE]

We use the Dunford-Taylor formula (see [45], Eq. (5.47) on page 44) to express the right hand side of Eq. (24) as a contour integral. Let $f$ be an analytic function on the complex plane. The Dunford-Taylor formula states that

[TABLE]

where $\gamma$ is a closed counterclockwise-oriented contour that encompasses a region of the complex plane that contains all eigenvalues of $A$. Using (25) in (24), we readily obtain Eqs. (14) and (15).

The quantity $S_{N}(t;A)$ is determined by both the statistics of eigenvalues and eigenvectors of $A$. This is most clearly seen by considering the special case where $A$ is diagonalizable. In this case, $A$ can be written in its canonical form

[TABLE]

where $\lambda_{j}$ are its eigenvalues, and $|r_{j}\rangle$ and $\langle\ell_{j}|$ form a bi-orthonormal set of right and left eigenvectors. Plugging this canonical form of $A$ in Eq. (24) we obtain

[TABLE]

where

[TABLE]

encode the eigenvectors overlaps [36]. Additionally, if $A$ is a normal matrix ($[A,A^{T}]=0$), then $\langle r_{j}|r_{k}\rangle=\langle l_{k}|l_{j}\rangle=\delta_{jk}$, and

[TABLE]

Therefore, the non-orthogonality of eigenvectors is a primary source of transient behavior, since the eigenmodes can interfere constructively to deliver an initial amplification of the signal well before it eventually dies out [47, 46, 48, 49, 50, 51, 52, 53, 54].

Eq. (14) reduces the computation of $S(t)$ to a computation of the average two-point correlator $\overline{W}(z,w)$. Although the one-point correlator

[TABLE]

of sparse random graphs has been studied extensively in Ref. [42, 43, 41, 31], to our knowledge, the two-point correlator has not been considered before. In the next section, we show how to compute the two-point correlator $\mathcal{W}_{A}$ for tree matrices, and in the subsequent section we compute $\overline{\mathcal{W}}$ for the canonical model of random graphs with a prescribed degree distribution.

III Tree graphs

In this section, we present a mathematical method to compute the two-point correlator $\mathcal{W}_{A}(z,w)$, and thus also $S_{N}(t;A)$, under the sole assumption that the graph represented by $C_{ij}$ is a tree.

The mathematical method we employ is based on two ideas: the size-doubling trick presented in the first subsection, and a recursive implementation of the Schur formula, presented in the second subsection. In a third subsection we discuss how the mathematical formalism simplifies for oriented graphs. We say that a graph is oriented if

[TABLE]

for all pairs $(i,j)$.

III.1 Size doubling trick

We use the following identity

[TABLE]

which expresses the trace of

[TABLE]

in terms of the block trace of the inverse of the matrix

[TABLE]

The block trace $\mathrm{bTr}_{11}$ of a $2N\times 2N$ matrix $X$ is defined by

[TABLE]

Since $B$ is a matrix of size $2N\times 2N$ and $A$ is a matrix of size $N\times N$, we call this the size-doubling trick, which bears some similarity with the Hermitization method [55, 56] in non-Hermitian random matrix theory.

In order to proceed, we note that the block trace formula in (32) is similar to the block trace formula for the spectral density of a sparse graph, see Eq. (57) in Ref. [41]. In the next subsection we will exploit this mathematical similarity.

III.2 Recursion relations

We derive a set of recursion relations that will provide us with the diagonal elements $[B^{-1}]_{j,j}$. The recursions can be closed using that $A$ is the adjacency matrix of a (weighted) tree. In particular, we use the following property. Let $A^{(j)}$ be the matrix obtained from $A$ by removing the $j$-th row and column; $A^{(j)}$ is the adjacency matrix of the so-called cavity graph obtained by removing the $j$-th node from the original graph [57, 58, 59]. It then holds that $A^{(j)}$ is a forest of $|\partial_{j}|$ isolated trees, where we have used the notation

[TABLE]

for the neighborhood of $j$ and $|\partial_{j}|$ is the number of elements in $\partial_{j}$. We use the Schur inversion formula to derive the recursion relations. However, first we need to define the quantities that appear in the recursion, which requires a reshuffling of the elements of $B$.

III.2.1 Preliminary ordering of matrix elements

The matrix $B$ can be seen as a replicated version of the original matrix $A$, where each node of the original graph has two replicas, one with label $i$ and the other with label $i+N$ ($i=1,\ldots,N$), which are located therefore far apart in the matrix $B$. Therefore, we operate on the $2N\times 2N$ matrix $B$ to create a new matrix that preserves the same connectivity structure of the original graph (encoded in $C$). This is achieved by bundling together labels that refer to the same node.

More precisely, we permute the rows and columns of the matrix $B$. The permutation we perform defines the matrix $\tilde{B}$, whose entries are assigned according to the following operations:

[TABLE]

This permutation is performed through a similarity transformation $\tilde{B}=P^{T}BP$, where $P$ is a suitably defined permutation matrix.

We then obtain the permuted matrix $\tilde{B}$ which consists of diagonal $2\times 2$ blocks (labelled using the sans-serif font)

[TABLE]

and off-diagonal blocks of the form

[TABLE]

for $i,j=1,\ldots,N$. Note that now the matrix $\tilde{B}$ is a block matrix (formed by $2\times 2$ blocks) that inherits the same connectivity structure as the matrix $A$ (or $C$), since elements of $A$ referring to the same node - which were located far apart in the matrix $B$ - are now bundled together.

The elements of $B^{-1}$, which we need in Eq. (32), are related to the elements of $\tilde{B}^{-1}$ in the following way:

[TABLE]

since a permutation matrix is an orthogonal transformation.

Moreover, the block trace needed in Eq. (32) reads

[TABLE]

where in the first line we used Eq. (40), and in the second line the fact that the permutation matrix $P$ maps indices $j$ onto $2j-1$ if $1\leq j\leq N$.

Therefore, the objects of interest are now the elements $[\tilde{B}^{-1}]_{2j-1,2j-1}$ of the inverse matrix of $\tilde{B}$.

Defining the $2\times 2$ matrices $\mathsf{G}_{j}$ for $j=1,\ldots,N$ as

[TABLE]

the two-point correlator reads (see Eq. (41))

[TABLE]

and the one-point resolvent reads

[TABLE]

or

[TABLE]

In the next subsection we derive a set of recursion relations for $\mathsf{G}_{j}$’s using Schur inversion formula.

III.2.2 Schur formula

We employ the Schur inversion formula

[TABLE]

with $S_{D}=(A-BD^{-1}C)^{-1}$ the inverse of the Schur complement of $D$ and $S_{A}=(D-CA^{-1}B)^{-1}$ the inverse of the Schur complement of $A$.

First, we show how this works for $j=1$ and later we generalize for arbitrary $j$. In order to implement the Schur inversion formula, we represent $\tilde{B}$ with the block matrix structure of the form

[TABLE]

where $\tilde{\mathsf{B}}_{11}$ is the $2\times 2$ matrix defined in Eq. (38), $\tilde{B}_{1\star}$ and $\tilde{B}_{\star 1}$ are $2\times 2(N-1)$ and $2(N-1)\times 2$ matrices respectively, and $\tilde{B}^{(1)}$ is a $2(N-1)\times 2(N-1)$ matrix. The matrix $\tilde{B}^{(1)}$ is the $\tilde{B}$ matrix of $A^{(1)}$ obtained by removing the first column and row from the matrix $A$.

Now, we are ready to find an equation for $\mathsf{G}_{1}$, the upper-left $2\times 2$ block of $\tilde{B}^{-1}$ (see Eq. (42)), taking full advantage of the block structure in Eq. (51). The Schur formula applied to the upper-left $2\times 2$ block of $\tilde{B}$ gives the following

[TABLE]

Now, using Eq. (38) and the fact that both $\tilde{B}_{1\star}$ and $\tilde{B}_{\star 1}$ are concatenations of matrices of the form $\mathsf{A}$ (see Eq. (39)), we obtain

[TABLE]

where

[TABLE]

Note that $\mathsf{G}_{kk}^{(1)}=\mathsf{G}_{k}^{(1)}$, the matrices we defined before in (42). In the sums in Eq. (53), we omit contributions from $k,\ell\notin\partial_{1}$ because $\mathsf{A}_{1k}$ and $\mathsf{A}_{\ell 1}$ are null matrices in this case (see Eq. (39)).

Since $A$ is the adjacency matrix of a (weighted) tree graph, it holds that $\mathsf{G}_{k\ell}^{(1)}$ is null for any $k,\ell\in\partial_{1}$ with $k\neq\ell$. To show this, note the following facts:

the nodes $k$ and $\ell$ belong to distinct trees in the forest represented by $A^{(1)}$; 2. 2.

the matrix $(\tilde{B}^{(1)})^{-1}$ has the mathematical form of a resolvent matrix $(u\mathbb{1}-X)^{-1}$, where $X$ is a block matrix built out of $2\times 2$ matrices located at the edges of the graph represented by $A$; 3. 3.

the matrix element $[X^{n}]_{k,\ell}$ denotes the sum of the weights of the paths in the graph of length $n$ that connect node $k$ to node $\ell$. If there exists no path that runs form $k$ to $\ell$, then $[X^{n}]_{k,\ell}=0$ for all $n$.

Therefore, expanding $(u\mathbb{1}-X)^{-1}=\sum^{\infty}_{n=0}X^{n}/u^{n+1}$, it follows that $[(u\mathbb{1}-X)^{-1}]_{k,\ell}=0$ (and therefore also $\mathsf{G}_{k\ell}^{(1)}$) if $k$ and $\ell$ belong to distinct trees.

Hence, applying that for tree matrices the $\mathsf{G}_{k\ell}^{(1)}$ are null if $k,\ell\in\partial_{1}$ with $k\neq\ell$, we can simplify Eq. (53) to

[TABLE]

This equation can be generalized to an arbitrary node $j$ because the Schur formula works the same way upon a permutation that is equivalent to relabelling of nodes. Hence, for an arbitrary node $j=1,\ldots,N$, Eq. (55) becomes

[TABLE]

This set of equations is not closed, because in general $\mathsf{G}_{j}\neq\mathsf{G}_{j}^{(k)}$ for $k\in\partial_{j}$. To close the equations, we can repeat the same procedure on the $N$ graphs $A^{(j)}$ obtained by removing one node at a time. We then obtain

[TABLE]

where $j\in\left\{1,2,\ldots,N\right\},k\in\partial_{j}$, and $\mathsf{G}^{(k,j)}_{\ell}$ are defined analogously to Eq. (54) but on the graph where the two nodes $k,j$ have been removed.

Note that in general $\mathsf{G}^{(k,j)}_{\ell}\neq\mathsf{G}^{(k)}_{\ell}$, and the recursion has to be continued. However, for tree matrices the recursion can be closed at the second step. Indeed, since nodes $\ell$ and $j$ belong to distinct trees on the graph $A^{(k)}$, the further removal of node $j$ does not affect $\mathsf{G}^{(k)}_{\ell}$ and therefore $\mathsf{G}^{(k,j)}_{\ell}=\mathsf{G}^{(k)}_{\ell}$. We finally obtain a closed set of equations

[TABLE]

which are valid for all $k\in\partial_{j}$ on tree graphs.

Solving the Eqs. (56) and (58) on one instance of a tree graph, one would get access to the corresponding two-point correlator $\mathcal{W}_{A}(z,w)=-\sum_{j=1}^{N}[\mathsf{G}_{j}]_{1,1}$. Since the Eqs. (56) and (58) are local equations, i.e. the right-hand sides of the Eqs. (56) and (58) only depend on the local neighborhood of node $j$, we can also apply Eqs. (56) and (58) to graphs that are locally tree-like.

III.3 Exact solution for oriented trees

We show now how Eqs. (56) and (58) simplify in the case of oriented trees, i.e., trees for which all edges are unidirectional ($C_{jk}C_{kj}=0$ for all $k\neq j$).

In this case, the last terms in the denominator of the r.h.s. of Eqs. (56) and (58) have off-diagonal entries equal to zero. Therefore,

[TABLE]

where

[TABLE]

and

[TABLE]

The recursion relations given by Eqs. (60)-(61) can be solved numerically on a fixed instance and the two-point correlator is given by

[TABLE]

whereas for oriented trees the one-point resolvent reads [43, 31]

[TABLE]

In practice, we will often be interested in the dynamics on graphs that are locally tree-like and locally oriented. Eqs. (60)-(61) also apply to these graph ensembles since the right-hand side of the Eqs. (60)-(61) only depends on the local neighborhood of node $k$. In the next section, we consider an important class of random graphs that are locally tree-like and locally oriented.

IV Directed random graphs with a prescribed degree distribution

We derive analytical expressions for $S(t)$ for the case when $A$ is the adjacency matrix of a directed random graph with a prescribed joint degree distribution, as defined in the introduction. This ensemble of random graphs is locally tree-like and locally oriented, and therefore the Eqs. (60)-(61) apply. To summarize, we will obtain the general formula

[TABLE]

where

[TABLE]

is the product of the mean out-degree (or in-degree)

[TABLE]

and the second moment $\overline{J^{2}}=\int\mathrm{d}x~{}p_{J}(x)x^{2}$ of the bond disorder. Eq. (64) is one of the main results of this paper.

From Eq. (64), we observe that $S(t)$ is universal in the sense that it only depends on $p_{J}(x)$ and $p_{\rm deg}(k)$ through the single parameter $r$. On the other hand, the dependence on $p_{D}$ is explicit.

The present section is organized as follows. In the first subsection, we discuss the spectral properties of adjacency matrices of directed random graphs with a prescribed joint degree distribution. In the second subsection, we derive the formula (64). In the third and fourth subsections, we derive explicit expressions for $S(t)$ by computing the contour integral in Eq. (64) for fixed diagonal disorder

[TABLE]

and for bimodal diagonal disorder,

[TABLE]

respectively.

IV.1 Spectral properties and asymptotic stability

The spectra of random graphs with a prescribed degree distribution have been studied in detail in the limit $N\rightarrow\infty$ in Refs. [43, 41, 31]. If $c>1$, then the graph contains a giant strongly connected component [39], which contributes a deterministic continuous part to the spectrum and may also contribute deterministic outliers.

The boundary of the continuous part of the spectrum consists of all $\lambda_{\rm b}$ solving

[TABLE]

The deterministic outliers $\lambda_{\mathrm{isol}}$ solve the relation

[TABLE]

Note that (i) outliers are always real, (ii) if $\overline{J}=0$, there are no outliers, (iii) the boundary of the continuous part and the location of outlier(s) are universal in the sense that they depend on $p_{J}(x)$ and $p_{\rm deg}(k)$ only through a single parameter (either $r$ or $c\overline{J}$).

IV.2 General diagonal disorder

Since directed random graphs with a prescribed degree distribution are locally tree-like and locally oriented, the Eqs. (60)-(61) apply. We obtain an analytical expression for the average

[TABLE]

by taking the ensemble average of Eq. (60).

Using that (i) the right-hand-side of Eq. (60) is identical to the right-hand-side of Eq. (61) for any $j,k$ with $C_{jk}=1$; (ii) all variables on the r.h.s. of Eq. (60) are statistically independent, as degree-degree correlations are absent; and (iii) all nodes in the ensemble are statistical equivalent, we obtain the equation

[TABLE]

Solving (72) for $\overline{\alpha}$ and using

[TABLE]

we obtain the general formula Eq. (64) for $S(t)$.

It remains to perform the contour integral in (64) for some specific choice of the diagonal disorder, which is the subject of the next subsections.

IV.3 Fixed diagonal disorder

We compute the contour integral in Eq. (64) in the case of fixed disorder given by Eq. (67). The integral (64) can be performed using residues. Changing variables $z^{\prime}=z+\mu$ and $w^{\prime}=w+\mu$, we obtain

[TABLE]

Quite remarkably, Eq. (74) for sparse oriented graphs and Eq. (11) for fully connected structures share the same functional form, and therefore the two models fall into the same universality class. Indeed, if $\overline{J}=0$, then $r$ is the spectral radius $\rho$ of $A$ (see (69)).

For large $t$,

[TABLE]

where the exponent describes the asymptotic growth or decay at a rate $r-\mu$. The rate $r-\mu$ contains a positive contribution $r$, which is the amplification of the initial perturbation when it spreads throughout the network, and a negative contribution $-\mu$, which is the local decay rate.

For $t\approx 0$, it holds that

[TABLE]

independent of the network structure.

While the initial response is independent of the network structure and leads to a decrease in $S(t)$, the response at larger $t$ will depend on how the initial perturbation spreads through the network. In particular, if $r-\mu>0$, then the response $S(t)$ is non-monotonic.

It is insightful to interpret the result given by Eq. (75) in the light of the spectral properties of the ensemble, as discussed in Sec. IV.1. The boundary of the continuous part of the spectrum is according to Eq. (69) a circle of radius $r$ centered around $-\mu$. As a consequence, the eigenvalue with the maximum real part, which belongs to the continuous spectrum, is given by

[TABLE]

This implies that the asymptotic dynamics of $S(t)$, given by Eq. (75) for large $t$, reads

[TABLE]

and the asymptotic rate is governed by the boundary of the continuous part of the spectrum. The outlier

[TABLE]

plays no role in the behavior of $S(t)$, even when $\lambda_{\rm isol}>\lambda^{\prime}_{\rm b}$. This is because the initial perturbation is random and therefore – in the limit $N\rightarrow\infty$ – it stands orthogonal to the one-dimensional eigenspace spanned by the outlier.

IV.4 Bimodal diagonal disorder

In the case of bimodal diagonal disorder with probability density $p_{D}$ given in Eq. (68), we evaluate the integral in Eq. (64) using residues. After some calculations, presented in Appendix B, we obtain

[TABLE]

where

[TABLE]

is a series involving confluent hypergeometric functions ${}_{1}F_{1}$ [79]. The confluent hypergeometric function is defined by the series

[TABLE]

where

[TABLE]

is the falling factorial and $(a)^{\underline{0}}=1$. Although not immediately evident, the formula (80) is symmetric in $\mu_{1}$ and $\mu_{2}$ upon the exchange $q\to 1-q$ as it should, due to a transformation formula of the confluent hypergeometric function upon change of sign of the main argument.

The formula (80) is rather complicated, but it has an appealing dynamical interpretation. In this example, the dynamical system consists of two sub-populations with decay rates $\mu_{1}$ and $\mu_{2}$, respectively. Neglecting interactions between the subpopulations, each of these would evolve in isolation according to Eq. (74), albeit with reduced connectivities $qc$ and $(1-q)c$, respectively. The first two terms in (80) describe precisely the dynamics of the populations in isolation. The third term instead describes the dynamical interference between the two sub-populations. We expect that the interference will be important in the limit of large $t$.

For $t\approx 0$, it holds that

[TABLE]

which is again independent of the network structure.

The analysis of the $t\rightarrow\infty$ limit of Eq. (80) is more complicated because of the nontrivial form of the interference term. Nevertheless, based on the analysis in the previous subsection for fixed diagonal disorder, we expect that for large $t$

[TABLE]

where

[TABLE]

is the value of $\lambda_{b}$, located at the boundary of the continuous part of the spectrum, with the largest real part. In the present example with $p_{D}$ the sum of two Dirac distributions, we obtain (see Appendix C)

[TABLE]

where

[TABLE]

Hence, for bimodal diagonal disorder even the asymptotic rate $\lambda^{\prime}_{\rm b}$ does not admit a simple expression, which clarifies why it is not a simple task to analyze the Eq. (80) in the limit of large $t$.

V Numerical experiments on random matrices and on empirical networks

We compare the obtained analytical expression for $S(t)=\lim_{N\rightarrow\infty}S_{N}(t)$, given by Eq. (74) and Eqs. (80)-(81), with numerical results for $S_{N}(t)$ on random matrices and on empirical networks.

V.1 Random graphs

We analyze the dynamics described by Eq. (2) with $A$ being the adjacency matrix of a weighted random graph with a prescribed degree distribution (see the Introduction for the definition of this ensemble). Here, we use a Poissonian distribution for the in-degrees (and out-degrees), namely,

[TABLE]

in Eq. (13). In addition, in Appendix D, we show results for power-law degree distributions, which are relevant to describe real-world networks, see Refs. [16, 60, 61, 62, 63].

In Fig. 1, we compare the obtained expressions for $S(t)$ with estimates of $S_{N}(t)$ obtained from numerical experiments. We consider three qualitatively different scenarios: (a) the system is both transiently and asymptotically stable ($\lambda^{\prime}_{{\rm b}}<0$ and $\lambda_{\mathrm{isol}}<0$), (b) the system is transiently stable but asymptotically unstable ($\lambda^{\prime}_{{\rm b}}<0$ and $\lambda_{\mathrm{isol}}>0$), and (c) the system is unstable ($\lambda^{\prime}_{{\rm b}}>0$). Since the theory in Eq. (64) is obtained by taking the limit $N\to\infty$ at fixed time, while in the simulations we work with a fixed system size $N$ and look at its time evolution, there will be a crossover time $t^{\star}(N)$ after which the theory and simulations are expected to diverge (see Appendix A for a detailed discussion). However, for $t<t^{\star}(N)$, we observe perfect agreement between theory and simulations, proving that the precise connectivity structure and the type of bond disorder do not matter in the transient dynamics of large dynamical systems. Only the combination $r^{2}=c\overline{J^{2}}$ plays a role, as predicted by the theory. In case (a), theory and simulations are in very good agreement, and $1/(2(r-\mu))$ provides the typical relaxation time to stationarity.

We clearly see the effect of the crossover time $t^{\star}(N)$ in cases (b) and (c) of Fig. 1 [but this would also become apparent in case (a) when $t$ is large enough and if a different scale were used]. For $t\gg t^{\star}(N)$, the contribution from the eigenmode with largest – and positive – real part $S_{N}(t;A)\sim a_{N}e^{2t\mathrm{Re}[\lambda_{1}]}$ (where $a_{N}$ is independent of $t$ and vanishes for large $N$) will govern the large-$t$ behavior of every single realization of the $N\times N$ numerical experiment, and as a consequence also of the average $S_{N}(t)=\overline{S_{N}(t;A)}$.

In order to capture both the transient and asymptotic behaviors, we can modify $S(t)$ by including the effect of the largest eigenmode as

[TABLE]

where $a$ and $b$ are two (albeit non-universal) fitting parameters. This is illustrated in Fig. 1 (b) and (c). We observe empirically that universality is broken after $t^{\star}(N)$: this is because the exponential growth of $S_{N}(t;A)$ at large times – for every individual realization of the matrix $A$ – is very sensitive to the precise value of $\lambda_{1}(A)$, which of course fluctuates from one realization to another. The average $S_{N}(t)=\overline{S_{N}(t;A)}$ is therefore disproportionally affected by the largest $\lambda_{1}(A)$ across the sample used for the simulations.

Note also that in the panel (b) the fitted exponent $b\approx\lambda_{\mathrm{isol}}$, whereas in panel (c) the exponent $b$ is different from $\max\mathrm{Re}[\lambda_{\rm b}]$.

V.2 Empirical networks

We test how well Eq. (74) describes the transient dynamics of a dynamical system defined on a real-world graph. To this aim, we numerically solve the dynamics in Eq. (2) for a matrix of the form Eq. (12), with $C_{ij}$ the entries of the adjacency matrix of a real world network. We consider two examples of real-world networks, namely, a food web and a signaling network.

For the food web example, we have chosen the food web of Otago Harbour, an intertidal mudflat ecosystem in New Zealand, whose adjacency matrix is determined in Ref. [69]. This adjacency matrix contains 180 nodes and 1924 edges. Nodes represent species, which include among others, plants, annelids, birds, and fish. Links represent trophic interactions between species, including among others, predation, parasitic castration, and macroparasitism interactions. The food web of Otago Harbour is an (almost) oriented network: $97\%$ of the links are oriented. The mean out-degree is $c=1924/180\approx 10.7$.

For the example of a cellular signaling network, we have considered a network of molecular interactions between signaling proteins in human cells. We have collected data from the NCI-Nature Pathway Interaction database [70], which is now available in the Network Data Exchange, NDEx [71, 72, 73]. We have extracted the adjacency matrix formed from the nodes in the 2-step neighborhood of the protein kinases MAPK1, AKT1, JAK1 and the protein APC. The resultant adjacency matrix contains 2288 nodes and 23207 edges. The mean out-degree is $c\approx 13.3$ and $70\%$ of the edges are oriented.

In Fig. 2 we present simulated dynamics on these networks. Numerical results show that the theory describes very well the transient dynamics on the food web (for all $t<t^{\star}(N)$), in spite of its small size $N=180$, while the dynamics of the signaling network is less well captured. Plotting the full spectra (see insets of Fig. 2), we see that the largest eigenvalue of the food web is almost equal to $\sqrt{c}$ (predicted by Eq. (69)), whereas for the signaling network it is much larger than $\sqrt{c}$, which clarifies why the theory works better for the food web.

In order to further demonstrate the relevance of Eq. (74) in modeling dynamics on real networks, we compare $S_{N}(t)$ with a naive theory based on either an exponential decay $e^{-2t\mu}$ or on exponential growth $e^{2t(r-\mu)}$. We observe that the naive theories do not capture well the transient response of the dynamical system. This is because the response of the real system consists of an initial decay at rate $\mu$ and asymptotic growth at rate $r-\mu$, which govern how the initial random perturbation propagates through the networks.

VI Conclusions and Outlook

We have developed a mathematical formalism that allows one to study how the transient response to an initial perturbation of a large dynamical system, captured by the observable $S_{N}(t)$, depends on the topology of the underlying network of interactions. The developed method allows one to compute the limiting value $S(t)=\lim_{N\rightarrow\infty}S_{N}(t)$ for graphs that have a locally tree-like structure. As an example, we have studied the transient dynamics on directed random graphs with prescribed degree distribution, which are canonical models for real-worlds systems, such as, the Internet [64, 39, 65], neural networks [66] and other high-dimensional systems [40, 17, 16]. We have found that the transient response of large systems is universal, in the sense that $S(t)$ only depends on the network topology through the single parameter $r^{2}=c\overline{J^{2}}$, which is the product of the mean degree and the second moment of the distribution of interactions strengths. On the other hand, the dependence on the statistics of diagonal elements, given by the distribution $p_{D}$, is non-universal.

The developed method is fairly general, and therefore various other examples can be considered. For example, it should be straightforward to extend the presented results to random graphs with correlations between in-degrees and out-degrees [31] or to non-oriented systems defined on regular graphs [67]. More challenging, but still feasible, is to extend the theory to graphs that contain many small cycles [76, 77, 78].

From a random matrix theory point of view, in the present paper we have developed the mathematical theory for the 2-point correlator $\mathcal{W}_{A}(w,z)$ of sparse random graphs. We can extend the present formalism to account for higher order $2n$-correlators, which provide information on higher-order moments of $S_{N}$ [75].

Since the theoretical response functions, given by Eqs. (74) and (80), only depend on a few parameters and describe well the transient response of dynamical systems defined on real-world networks, we believe that Eqs. (74) and (80) can be used to infer properties of networks from time-series data.

Acknowledgements.

WT is supported by ETIUDA scholarship 2018/28/T/ST1/00470 from National Science Center and the Diamond Grant 0225/DIA/2015/44 from the Polish Ministry of Science and Higher Education. WT is grateful to King’s College London for hospitality, where this work was done. PV and IN acknowledge funding by the Engineering and Physical Sciences Research Council (EPSRC) through the Centre for Doctoral Training in Cross Disciplinary Approaches to Non-Equilibrium Systems (CANES, Grant Nr. EP/L015854/1). We thank Ahmed Zaid for making us aware of a sign typo in some equations of the paper.

Appendix A $N$-dependence of the transition time $t^{\star}(N)$

We study here how the crossover time $t^{\star}(N)$ depends on $N$. The crossover time determines when the system changes from its transient regime, where $S_{N}(t)=\overline{\langle|\mathbf{y}(t)|^{2}\rangle}\approx S(t)=\lim_{N\rightarrow\infty}S_{N}(t)$, to the asymptotic regime, where $S(t)\ll S_{N}(t)$. For transiently stable systems, we define the crossover time $t^{\star}(N)$ by Eq. (9).

In the first subsection, we determine the $N$-dependence of $t^{\star}(N)$ with a brute-force numerical analysis on four random matrix ensembles. In the second subsection, we analyze $t^{\star}(N)$ in the case of the normal Ginibre ensemble, a case for which we have the rare luxury of being able to carry out the analytical treatment in full. We obtain an explicit analytical expression for $S_{N}(t)$, which can be used to evaluate numerically its minimum value. We find that $t^{\star}(N)\approx 1.15\sqrt{N}$.

A.1 Numerical results for four ensembles

We study the first-order dynamics in Eq. (2) for

[TABLE]

In Fig. 3 we present numerical results for $S_{N}(t)$ as a function of $t$ for different values of $N$ in the case of $X_{ij}$ being i.i.d.  random variables taken from a Gaussian distribution. Fig. 3 shows that $S(t)\approx S_{N}(t)$ for small enough $t$ (transient regime) while $S(t)\ll S_{N}(t)$ when $t$ is large enough (asymptotic regime). Moreover, we observe that the crossover time $t^{\star}(N)$ from the transient regime to the asymptotic regime grows as a function of $N$.

In Fig. 4 we analyse how the crossover time $t^{\star}(N)$ depends on $N$. For this we use its definition given by Eq. (9). We consider four random matrix ensembles:

1. (i)

Ginibre matrices with zero mean [74]: the $X_{ij}$ are i.i.d. random variables taken from a Gaussian distribution with zero mean and variance $1/N$, and we set $\mu=-1$; 2. (ii)

Gaussian Orthogonal Ensemble: the $X_{ij}=X_{ji}$ are i.i.d. real-valued random variables taken from a Gaussian distribution with zero mean and variance $1/N$, and we set $\mu=2$; 3. (iii)

Ginibre matrices with nonzero mean: the $X_{ij}$ are the i.i.d.random variables taken from a Gaussian distribution with mean $6/N$ and variance $20/N$, and we set $\mu=5$; 4. (iv)

Adjacency matrices of sparse random graphs: we set $X_{ij}=J_{ij}C_{ij}$ with $C_{ij}$ the adjacency matrix of a random graph with a prescribed joint distribution of in-degrees and out-degrees given by

[TABLE]

with a mean in-degree (out-degree) $c=2$. In addition, we weigh the edges with couplings $J_{ij}$ that are i.i.d. random variables taken from a Gaussian distribution with mean $3$ and variance $1$. We set $\mu=5$.

In cases (i), (ii) and (iv) the numerical results for $t^{\star}$ are well fitted by the power law $t^{\star}=\alpha N^{\beta}$. In the case (iii) a pure power-law curve does not describe the data as well as the function $t^{\star}=\alpha N^{\beta}\log N$, i.e., a power-law with a logarithmic correction.

Interestingly, for the Ginibre ensemble and GOE (cases (i) and (ii)) the fitted exponents are approximately $1/2$ and $2/3$, respectively, which is the exponent governing the scale of typical fluctuations at the edge of the spectrum [80, 81]. This pattern, however, does not seem to carry over to cases (iii) and (iv). We are not aware of a theory that can shed more light on the scaling of $t^{\star}$ with $N$ in these cases.

In conclusion, numerical results indicate that $t^{\star}$ diverges with $N$ with a law that might be related to the scale of typical fluctuations at the edge of the spectrum. In the next section, we analyze in more detail the normal Ginibre ensemble, which allows for a more complete analytical treatment.

A.2 Analytical treatment of $t^{\star}(N)$ for the normal Ginibre ensemble

We analyze the $N$-dependence of $t^{\star}(N)$ for a rare example where a full analytical treatment is possible, namely the normal Ginibre ensemble. For this ensemble, we derive an explicit analytical expression for $S_{N}(t)$. Evaluating numerically its minimum, we obtain that

[TABLE]

A.2.1 Definition of the normal Ginibre ensemble

We consider dynamical systems of the type given by Eq. (2) with

[TABLE]

where $X^{\prime}$ is a matrix drawn from the normal Ginibre ensemble.

Normal Ginibre matrices are obtained by retaining the normal term in the Schur decomposition of a complex Ginibre matrix. The entries of the complex Ginibre ensemble [74] are defined as follows,

[TABLE]

where $x_{jk}$ and $y_{jk}$ are i.i.d. random variables from a Gaussian distribution with zero mean and variance $1/2$. The matrices in the ensemble are non-symmetric. The joint pdf of matrix entries reads

[TABLE]

where $C_{N}$ is the normalization constant.

In order to obtain a normal Ginibre ensemble, we implement the Schur decomposition

[TABLE]

with $U$ a unitary matrix, $\Lambda$ the diagonal matrix of eigenvalues, and $T$ a strictly upper triangular matrix. From Eq. (97) it is apparent that the upper-triangular matrix $T$ is the source of non-normality. Therefore, we drop the upper-triangular matrix to obtain

[TABLE]

which defines a matrix sample drawn from the normal Ginibre ensemble.

Implementing the Schur decomposition in Eq. (96), we find that $P_{X}$ factorizes into

[TABLE]

where

[TABLE]

is the distribution of complex eigenvalues in the normal Ginibre ensemble.

A.2.2 Calculation of $S(t)$

Because of Eqs. (5), (6), and (24), we can write

[TABLE]

Taking the limit $N\rightarrow\infty$ and using that the eigenvectors of normal (Ginibre) matrices are orthonormal, we obtain that

[TABLE]

where the integration is over the entire complex plane and $\rho(z,\bar{z})$ is the spectral density of the normal Ginibre ensemble in the limit $N\rightarrow\infty$. Since in this limit, the spectral density is uniform in the disk of unit radius [82],

[TABLE]

with $\mathbf{1}_{T}$ the indicator function (equal to $1$ when $T$ is true and [math] when $T$ is false), we obtain the expression

[TABLE]

where $I_{0}$ and $I_{1}$ are modified Bessel functions. This result is the analog of Eqs. (11) and Eqs. (64) for dense and normal matrices.

To compute $t^{\star}$ as a function of $N$, though, the analysis above is not sufficient, and we need to resort to a finite-$N$ calculation.

A.2.3 Calculation of $S_{N}(t)$

For finite $N$, Eq. (102) is modified as

[TABLE]

where $\rho_{N}(z,\bar{z})$ is now the spectral density of the normal Ginibre ensemble at finite matrix sizes $N$. It holds that [83]

[TABLE]

Therefore, we obtain

[TABLE]

It is convenient to rescale $r\to\frac{s}{\sqrt{N}}$ to get

[TABLE]

The integral can be performed with some work to yield eventually the exact result

[TABLE]

where ${}_{1}F_{1}(a;b;z)$ is the Kummer confluent hypergeometric function, defined in Eq. (82).

A.2.4 Computing the crossover time $t^{\star}(N)$

The crossover time $t^{\star}$ is defined as the minimum of Eq. (109). By equating the derivative of Eq. (109) to zero, we obtain that $t^{\star}$ satisfies

[TABLE]

We have solved Eq. (LABEL:trans) numerically for $N<400$ and obtained that $t^{\star}$ is well fitted by $t^{\star}(N)\approx\alpha N^{\beta}$ with a fitted exponent $\beta\approx 1/2$, see Fig. 5. A careful asymptotic analysis of Eq. (LABEL:trans) is beyond the scope of this paper and is left for future investigation.

Appendix B Derivation of Eqs. (80) and (81) for $S(t)$ in the case of bimodal disorder

We explicitly solve the contour integral in Eq. (64) for the case of a bimodal distribution

[TABLE]

with $\mu_{1},\mu_{2}\geq 0$ and $q\in[0,1]$, and derive the expression for $S(t)$ given by Eqs. (80) and (81).

Without loss of generality we can consider a shifted distribution

[TABLE]

since

[TABLE]

which follows from the variable transformation $z\to z-\mu_{1}$ and $w\to w-\mu_{1}$ in Eq. (64).

Therefore, we evaluate the quantity

[TABLE]

Using

[TABLE]

in Eq. (64), we find that

[TABLE]

where $r^{2}=c\overline{J^{2}}$ is as defined in Eq. (65).

B.1 Series expression for $\hat{S}(t)$

In order to compute the contour integral in Eq. (116), we express $\hat{S}(t)$ as an infinite series, which we can integrate term by term. To this aim, we use the geometric series

[TABLE]

in Eq. (116), and find

[TABLE]

Using the binomial theorem, we obtain the double sum

[TABLE]

where we have denoted

[TABLE]

B.2 Performing the contour integral

To compute the contour integrals in $\varphi(n,m;q,-\mu;t)$, we need to evaluate the residues at $z=0$ and $z=-\mu$, and sum them up,

[TABLE]

The two residues are given by

[TABLE]

and

[TABLE]

respectively, where

[TABLE]

is the Tricomi hypergeometric function (see Eq. (13.2.7) in Ref. [79]), $\Gamma(x)=\int^{\infty}_{0}y^{x-1}e^{-y}{\rm d}y$ is the Gamma function, and ${}_{1}F_{1}(a;b;z)$ is the Kummer confluent hypergeometric function as defined in Eq. (82). Note that although the definition given by Eq. (124) does not apply for integer $b$, the Tricomi hypergeometric function exists also for integer $b$ by continuity of the rhs of Eq. (124).

We illustrate the derivation of Eq. (122). Using the Leibniz formula for the derivative of a product, we obtain

[TABLE]

where $(a)^{\underline{s}}$ is the falling factorial as defined in Eq. (83). Taking the limit $z\rightarrow 0$ in Eq. (LABEL:Tric2) and rearranging terms with the help of the Kummer transformation [79, Eq. (13.2.40)]

[TABLE]

we arrive at Eq. (122).

Equation (LABEL:Tric2) is obtained using an analogous reasoning.

Finally, thanks to the identity [79, Eq. (13.2.41)] the sum of (122) and (LABEL:Tric2) can be simplified as

[TABLE]

in terms of the Kummer confluent hypergeometric function defined in Eq. (82).

B.3 Meaning of the different terms

Plugging Eqs. (127) and (120) into the double series Eq. (119), provides us with an explicit expression for $\hat{S}(t)$. Unfortunately, this mathematical expression is not meaningful yet. However, we will see that the $n=0$, $n\in\left\{1,2,\ldots,m\right\}$, and $n=m+1$ terms in the double series can be given an appealing interpretation, making sense of the formula.

Isolating the $n=0$ and $n=m+1$ terms in (119), we obtain

[TABLE]

Using Eqs. (113) and (128), using the explicit forms

[TABLE]

and using the series expressions

[TABLE]

for the modified Bessel function $I_{0}$, we obtain the final result given by Eqs. (80) and (81).

Appendix C Spectra of random graphs with a prescribed degree distribution and bimodal diagonal disorder

We analyze the spectral properties of adjacency matrices $A$ of weighted random graphs with a prescribed degree distribution — as defined in the introduction of this paper – for the case where the distribution of diagonal elements is given by Eq. (68).

If the mean in-degree $c>1$, then the spectrum consists of a continuous part with boundary given by Eq. (69) and with outliers solving Eq. (70). Substitution of Eq. (68) in Eq. (69) leads to the relation

[TABLE]

for the boundary of the continuous part of the spectrum. Analogously, substituting Eq. (68) in Eq. (70) we obtain that outliers solve

[TABLE]

In Fig. 6, we plot the spectra for a few matrix instances together with the theoretical boundary given by Eq. (133). If $r^{2}$ is small enough, then the continuous part of the spectrum consists of two disconnected sets, which are centered around $-\mu_{1}$ and $-\mu_{2}$, respectively. On the other hand, for large $r$ the continuous spectrum is a simply connected set in the complex plane.

In order to obtain the leading eigenvalue $\lambda_{1}$, i.e., the eigenvalue with the largest real part, we determine the real solutions of Eq. (133). If

[TABLE]

where

[TABLE]

then Eq. (133) admits two real solutions, namely,

[TABLE]

and

[TABLE]

In this scenario, the continuous part of the spectrum is simply connected. This is illustrated in panels (b) and (c) of Fig. 6. On the other hand, if

[TABLE]

then Eq. (133) admits four real solutions, namely, $\lambda_{\rm b,1}$, $\lambda_{\rm b,2}$,

[TABLE]

and

[TABLE]

This the scenario illustrated in panel (a) of Fig. 6.

It follows that

[TABLE]

Solving Eq. (134), we find that there may exist two outlier eigenvalues, namely,

[TABLE]

where

[TABLE]

Hence, the leading eigenvalue

[TABLE]

Appendix D Power-law random graphs

Degree distributions of real-world graphs are often well described by power laws, see for instance Refs. [16, 60, 61, 62, 63]. Therefore, we compare in this appendix the theoretical result Eq. (74) for $S(t)$ with numerical results of $S_{N}(t)$ on random graphs with a power law degree distribution.

In Figure 7 we present results for random graphs with the prescribed degree distribution

[TABLE]

where $\zeta(x)$ is the Riemann-zeta function. We have set $\overline{J}=0$ and $\mu$ large enough such that we are in the stable case, as in sub-panel (a) of Fig. 1 in the main text ( all eigenvalues are negative, as ${\rm max_{\lambda_{{\rm b}}}}{\rm Re}[\lambda_{\rm b}]=\sqrt{\zeta(3)/\zeta(2)}-1.326<0$ and there is no outlier). We observe that the theoretical expression Eq. (74) for $S(t)$ captures very well the empirical results at small enough $t$ and finite $N$. On the other hand, for larger values of $t$ we observe significant sample-to-sample fluctuations.

We can thus conclude that the theoretical expression $S(t)$ describes well the dynamical response on random graphs with $\langle k^{\gamma}\rangle=\infty$ for $\gamma\geq 2$ as long as $t<t^{\star}(N)$, which is consistent with results on the leading eigenvalue of directed random graphs with power-law degree distributions obtained in [31].