Universal transient behavior in large dynamical systems on networks

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.


I. INTRODUCTION
Networks of interacting constituents appear in the study of systems as diverse as ecosystems [1][2][3], neural networks [4][5][6][7][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 ζ i (t).For example, ζ 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 for i = 1, . . ., N .Although the functions f i are arbitrary, the dynamics can be linearized close to a stationary point ζ ⋆ for which f i (ζ ⋆ ) = 0 to yield [28,29] where the N -dimensional vector encodes deviations from the stationary state and are the entries of the corresponding Jacobian matrix A.
The magnitude of the deviation vector y(t) as a function of time is captured by the squared norm |y(t)| 2 .In order to grasp the dynamical response to generic initial perturbations, we consider the average over initial states y(0) uniformly drawn from the sphere |y(0)| = α, where α quantifies the strength of the initial kick; without loss of generality we set α = 1.In addition, to capture the properties of a typical dynamical system, we take A to be a random matrix of pairwise interactions [3,18,30], and we further average the squared norm over the disorder Since we will be interested in large systems, we take the limit If then we say that a system is transiently (un)stable.Note that transient stability is different from asymptotic stability, which arXiv:1906.10634v5 [nlin.AO] 14 Jan 2024 is governed by the sign of the real part of the leading eigenvalue of A, see e.g.Refs.[3,18,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 ⋆ (N ) such that for times t < t ⋆ (N ) theory and experiment are in correspondence, i.e. S N (t) ≈ S(t), while for t ≫ t ⋆ (N ) theory and experiment are in disagreement, i.e. S N (t) ≫ 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→∞ S N (t) = ∞ but lim t→∞ S(t) = 0.
For transiently stable systems, where S N (t) has a nonmonotonic behavior, we define the crossover time by The crossover time t ⋆ (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 [32][33][34][35][36], namely for where the X ij are independent and identically distributed (i.i.d.) entries with zero mean and finite moments.In this case, where ρ is the spectral radius of the matrix X/ √ 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, where D i ∈ R − are the diagonal decay rates, the J ij ∈ R are the coupling strengths, and the C ij ∈ {0, 1} 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 → ∞ 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 [17,31,[39][40][41].These random graphs are used to model real-world systems, such as, the Internet [39,64,65], neural networks [66] and other highdimensional systems [16,17,40].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 deg (k in , k out ) of in-degrees k in and outdegrees k out .In addition, we assume that the degree distribution factorizes as The spectral properties of this ensemble have been studied in detail in Refs.[31,[41][42][43] 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 C 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 deg (k) through a single parameter (see our main formula (60)).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 ⋆ (N ), where t ⋆ (N ) is a timescale that diverges with N .For t ≫ t ⋆ (N ), it holds that S N (t; A) ∼ e 2Re[λ1(A)]t , with λ 1 (A) the eigenvalue of A with the largest real part, and as a consequence, S(t) ≪ 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 ⋆ (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 which expresses the dynamical response S N (t; A) as a contour integral of the two-point correlator over a closed counterclockwise-oriented contour γ that encloses all eigenvalues of A. The symbol 1 denotes the identity matrix of size N .Analogously, we obtain where 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 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 and therefore Since |y(0)| 2 = 1, there exists a matrix O in the the orthogonal group O(N ) (the group of isometries of the N -sphere) such that where Taking the average of (19) with respect to all initial conditions y(0) selected uniformly at random on the unit sphere is thus equivalent to taking the average of the following expression e T 1 O T e A T t e At Oe 1 (22) with respect to the uniform (Haar) measure on the orthogonal group.Using that [44] and Eqs. ( 22) and ( 5), we obtain 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 where γ 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 where λ j are its eigenvalues, and |r j ⟩ and ⟨ℓ j | form a biorthonormal set of right and left eigenvectors.Plugging this canonical form of A in Eq. ( 24) we obtain where encode the eigenvectors overlaps [36].Additionally, if A is a normal matrix ([A, A T ] = 0), then ⟨r j |r k ⟩ = ⟨l k |l j ⟩ = δ jk , and 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 [46][47][48][49][50][51][52][53][54].
Eq. ( 14) reduces the computation of S(t) to a computation of the average two-point correlator W (z, w).Although the one-point correlator of sparse random graphs has been studied extensively in Ref. [31,[41][42][43], 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 W A for tree matrices, and in the subsequent section we compute 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 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 for all pairs (i, j).

A. Size doubling trick
We use the following identity which expresses the trace of in terms of the block trace of the inverse of the matrix The block trace bTr 11 of a 2N × 2N matrix X is defined by Since B is a matrix of size 2N × 2N and A is a matrix of size N × 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.

B. 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 |∂ j | isolated trees, where we have used the notation for the neighborhood of j and |∂ j | is the number of elements in ∂ 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.

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, . . ., N ), which are located therefore far apart in the matrix B. Therefore, we operate on the 2N × 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 B, whose entries are assigned according to the following operations: (37) This permutation is performed through a similarity transformation B = P T BP , where P is a suitably defined permutation matrix.
We then obtain the permuted matrix B which consists of diagonal 2 × 2 blocks (labelled using the sans-serif font) and off-diagonal blocks of the form for i, j = 1, . . ., N .Note that now the matrix B is a block matrix (formed by 2 × 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 B−1 in the following way: since a permutation matrix is an orthogonal transformation.Moreover, the block trace needed in Eq. ( 32) reads where in the first line we used Eq. ( 40), and in the second line the fact that the permutation matrix P maps indices Therefore, the objects of interest are now the elements Defining the 2 × 2 matrices G j for j = 1, . . ., N as the two-point correlator reads (see Eq. ( 41)) and the one-point resolvent reads or In the next subsection we derive a set of recursion relations for G j 's using Schur inversion formula.

Schur formula
We employ the Schur inversion formula with 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 B with the block matrix structure of the form where B11 is the 2 × 2 matrix defined in Eq. ( 38), B1⋆ and B⋆1 are 2 × 2(N − 1) and 2(N − 1) × 2 matrices respectively, and is the 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 G 1 , the upper-left 2 × 2 block of B−1 (see Eq. ( 42)), taking full advantage of the block structure in Eq. (47).The Schur formula applied to the upper-left 2 × 2 block of B gives the following Now, using Eq. ( 38) and the fact that both B1⋆ and B⋆1 are concatenations of matrices of the form A (see Eq. ( 39)), we obtain where Note that G (1) k , the matrices we defined before in (42).In the sums in Eq. ( 49), we omit contributions from k, ℓ / ∈ ∂ 1 because A 1k and A ℓ1 are null matrices in this case (see Eq. ( 39)).
Since A is the adjacency matrix of a (weighted) tree graph, it holds that G (1) kℓ is null for any k, ℓ ∈ ∂ 1 with k ̸ = ℓ.To show this, note the following facts: 1. the nodes k and ℓ belong to distinct trees in the forest represented by A (1) ; 2. the matrix ( B(1) ) −1 has the mathematical form of a resolvent matrix (u1 − X) −1 , where X is a block matrix built out of 2 × 2 matrices located at the edges of the graph represented by A; 3. the matrix element [X n ] k,ℓ denotes the sum of the weights of the paths in the graph of length n that connect node k to node ℓ.If there exists no path that runs form k to ℓ, then [X n ] k,ℓ = 0 for all n.
Therefore, expanding (u1 Hence, applying that for tree matrices the G kℓ are null if k, ℓ ∈ ∂ 1 with k ̸ = ℓ, we can simplify Eq. ( 49) to 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, . . ., N , Eq. ( 51) becomes This set of equations is not closed, because in general for k ∈ ∂ 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 where j ∈ {1, 2, . . ., N } , k ∈ ∂ j , and G (k,j) ℓ are defined analogously to Eq. ( 50) but on the graph where the two nodes k, j have been removed.
Note that in general ℓ , and the recursion has to be continued.However, for tree matrices the recursion can be closed at the second step.Indeed, since nodes ℓ and j belong to distinct trees on the graph A (k) , the further removal of node j does not affect We finally obtain a closed set of equations which are valid for all k ∈ ∂ j on tree graphs.Solving the Eqs.( 52) and ( 54) on one instance of a tree graph, one would get access to the corresponding two-point correlator Since the Eqs.( 52) and (54) are local equations, i.e. the right-hand sides of the Eqs.( 52) and ( 54) only depend on the local neighborhood of node j, we can also apply Eqs. ( 52) and (54) to graphs that are locally tree-like.

C. Exact solution for oriented trees
We show now how Eqs. ( 52) and ( 54) simplify in the case of oriented trees, i.e., trees for which all edges are unidirectional (C jk C kj = 0 for all k ̸ = j).
In this case, the last terms in the denominator of the r.h.s. of Eqs. ( 52) and ( 54) have off-diagonal entries equal to zero.Therefore, and The recursion relations given by Eqs. ( 56)-( 57) can be solved numerically on a fixed instance and the two-point correlator is given by whereas for oriented trees the one-point resolvent reads [31,43] In practice, we will often be interested in the dynamics on graphs that are locally tree-like and locally oriented.Eqs. ( 56)-( 57) also apply to these graph ensembles since the righthand side of the Eqs.( 56)-( 57) 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.( 56)-( 57) apply.To summarize, we will obtain the general formula where is the product of the mean out-degree (or in-degree) and the second moment J 2 = ´dx p J (x)x 2 of the bond disorder.Eq. ( 60) is one of the main results of this paper.From Eq. ( 60), we observe that S(t) is universal in the sense that it only depends on p J (x) and p 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 (60).In the third and fourth subsections, we derive explicit expressions for S(t) by computing the contour integral in Eq. ( 60) for fixed diagonal disorder and for bimodal diagonal disorder, A. Spectral properties and asymptotic stability The spectra of random graphs with a prescribed degree distribution have been studied in detail in the limit N → ∞ in Refs.[31,41,43].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 λ b solving The deterministic outliers λ isol solve the relation Note that (i) outliers are always real, (ii) if 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 deg (k) only through a single parameter (either r or cJ).

B. General diagonal disorder
Since directed random graphs with a prescribed degree distribution are locally tree-like and locally oriented, the Eqs.( 56)-( 57) apply.We obtain an analytical expression for the average by taking the ensemble average of Eq. (56).
Using that (i) the right-hand-side of Eq. ( 56) is identical to the right-hand-side of Eq. ( 57) for any j, k with C jk = 1; (ii) all variables on the r.h.s. of Eq. ( 56) are statistically independent, as degree-degree correlations are absent; and (iii) all nodes in the ensemble are statistical equivalent, we obtain the equation Solving ( 68) for α and using we obtain the general formula Eq. ( 60) for S(t).
It remains to perform the contour integral in (60) for some specific choice of the diagonal disorder, which is the subject of the next subsections.

C. Fixed diagonal disorder
We compute the contour integral in Eq. ( 60) in the case of fixed disorder given by Eq. ( 63).The integral (60) can be performed using residues.Changing variables z ′ = z + µ and w ′ = w + µ, we obtain Quite remarkably, Eq. ( 70) 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 J = 0, then r is the spectral radius ρ of A (see (65)).For large t, where the exponent describes the asymptotic growth or decay at a rate r − µ.The rate r − µ contains a positive contribution r, which is the amplification of the initial perturbation when it spreads throughout the network, and a negative contribution −µ, which is the local decay rate.For t ≈ 0, it holds that 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 −µ > 0, then the response S(t) is non-monotonic.
It is insightful to interpret the result given by Eq. ( 71) in the light of the spectral properties of the ensemble, as discussed in Sec.IV A. The boundary of the continuous part of the spectrum is according to Eq. ( 65) a circle of radius r centered around −µ.As a consequence, the eigenvalue with the maximum real part, which belongs to the continuous spectrum, is given by This implies that the asymptotic dynamics of S(t), given by Eq. ( 71) for large t, reads and the asymptotic rate is governed by the boundary of the continuous part of the spectrum.The outlier plays no role in the behavior of S(t), even when λ isol > λ ′ b .This is because the initial perturbation is random and therefore -in the limit N → ∞ -it stands orthogonal to the onedimensional eigenspace spanned by the outlier.

D. Bimodal diagonal disorder
In the case of bimodal diagonal disorder with probability density p D given in Eq. ( 64), we evaluate the integral in Eq. ( 60) using residues.After some calculations, presented in Appendix B, we obtain where is a series involving confluent hypergeometric functions 1 F 1 [79].The confluent hypergeometric function is defined by the series where is the falling factorial and (a) 0 = 1.Although not immediately evident, the formula ( 76) is symmetric in µ 1 and µ 2 upon the exchange q → 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 ( 76) is rather complicated, but it has an appealing dynamical interpretation.In this example, the dynamical system consists of two sub-populations with decay rates µ 1 and µ 2 , respectively.Neglecting interactions between the subpopulations, each of these would evolve in isolation according to Eq. ( 70), albeit with reduced connectivities qc and (1−q)c, respectively.The first two terms in (76) 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 ≈ 0, it holds that which is again independent of the network structure.The analysis of the t → ∞ limit of Eq. ( 76) 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 S(t) ∼ e λ ′ b t (81) where is the value of λ 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) where Hence, for bimodal diagonal disorder even the asymptotic rate λ ′ b does not admit a simple expression, which clarifies why it is not a simple task to analyze the Eq. ( 76) 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 →∞ S N (t), given by Eq. ( 70) and Eqs. ( 76)-( 77), with numerical results for S N (t) on random matrices and on empirical networks.

A. 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, 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 (λ ′ b < 0 and λ isol < 0), (b) the system is transiently stable but asymptotically unstable (λ ′ b < 0 and λ isol > 0), and (c) the system is unstable (λ ′ b > 0).Since the theory in Eq. ( 60) is obtained by taking the limit N → ∞ 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 ⋆ (N ) after which the theory and simulations are expected to diverge (see Appendix A for a detailed discussion).However, for t < t ⋆ (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 = cJ 2 plays a role, as predicted by the theory.In case (a), theory and simulations are in very good agreement, and 1/(2(r − µ)) provides the typical relaxation time to stationarity.
We clearly see the effect of the crossover time t ⋆ (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 ≫ t ⋆ (N ), the contribution from the eigenmode with largest -and positive -real part S N (t; A) ∼ a N e 2tRe [λ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 × N numerical experiment, and as a consequence also of the average S N (t) = 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 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 ⋆ (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 λ 1 (A), which of course fluctuates from one realization to another.The average S N (t) = S N (t; A) is therefore disproportionally affected by the largest λ 1 (A) across the sample used for the simulations.Note also that in the panel (b) the fitted exponent b ≈ λ isol , whereas in panel (c) the exponent b is different from max Re[λ b ].

B. Empirical networks
We test how well Eq.( 70) 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 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 ≈ 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 ≈ 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 ⋆ (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 √ c (predicted by Eq. ( 65)), whereas for the signaling network it is much larger than √ c, which clarifies why the theory works better for the food web.
In order to further demonstrate the relevance of Eq. ( 70) in modeling dynamics on real networks, we compare S N (t) with a naive theory based on either an exponential decay e −2tµ or on exponential growth e 2t(r−µ) .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 3) are compared with the theoretical prediction S(t) (blue and red lines), given by Eq. ( 70), as well as with, e 2t(r−µ) and e −2µt (black and magenta lines).We have weighted the networks with couplings Jij that are i.i.d.random variables drawn from the distribution pJ (x) = (1/2)δ(x − 1) + (1/2)δ(x + 1) and we have used a single realization of the matrix J for both the food web and the signaling network.The diagonal elements are fixed to a constant −µ = −r + 0.27, such that the system is asymptotically unstable.Insets (a) and (b) show the spectra of the adjacency matrices A for the food web and the signaling network considered with random couplings.The red circle has radius r = √ c and is the predicted boundary of the spectrum according to (65).We have estimated SN (t) = ⟨|y(t)| 2 ⟩ by simulating the dynamics on the generated networks for 25 realizations of the initial condition y(0).system consists of an initial decay at rate µ and asymptotic growth at rate r − µ, 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 →∞ 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 [39,64,65], neural networks [66] and other high-dimensional systems [16,17,40].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 = cJ 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 nonuniversal.
The developed method is fairly general, and therefore vari-ous 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 2point correlator W A (w, z) of sparse random graphs.We can extend the present formalism to account for higher order 2ncorrelators, which provide information on higher-order moments of S N [75].
Since the theoretical response functions, given by Eqs. ( 70) and ( 76), only depend on a few parameters and describe well the transient response of dynamical systems defined on realworld networks, we believe that Eqs. ( 70) and ( 76) can be used to infer properties of networks from time-series data.
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) ≈ S N (t) for small enough t (transient regime) while S(t) ≪ S N (t) when t is large enough (asymptotic regime).Moreover, we observe that the crossover time t ⋆ (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 ⋆ (N ) depends on N .For this we use its definition given by Eq. ( 9).We consider four random matrix ensembles: (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 µ = −1; (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 µ = 2; (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 µ = 5; (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 indegrees and out-degrees given by 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 µ = 5.
In cases (i), (ii) and (iv) the numerical results for t ⋆ are well fitted by the power law t ⋆ = αN β .In the case (iii) a pure power-law curve does not describe the data as well as the function t ⋆ = αN β 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 ⋆ with N in these cases.
In conclusion, numerical results indicate that t ⋆ 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.Dashed black line is the curve S(t) = e −2µt I0(2ρt) (see Eq. ( 11)).Markers are averages over 2000 matrix samples with 25 realizations of the initial conditions for each sample.

Analytical treatment of t ⋆ (N ) for the normal Ginibre ensemble
We analyze the N -dependence of t ⋆ (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 ) a. Definition of the normal Ginibre ensemble We consider dynamical systems of the type given by Eq. (2) with FIG.4: Numerical results for the crossover time t ⋆ as a function of N for the four considered models: (i) Ginibre ensemble, (ii) Gaussian Orthogonal Ensemble, (iii) Ginibre ensemble with an outlier, (iv) sparse random graph with Poissonian connectivity, mean degree c = 2 and Gaussian bond disorder.Solid lines denote fitted functions to the empirical data.In cases (i) and (ii) we have set µ such that the system is at the edge of stability (in the N → ∞ limit, the average of the leading eigenvalue λ1 = 0), and therefore S(t) decays asymptotically as t −1/2 .The parameters in (iii) and (iv) are set such that the system is transiently stable but asymptotically unstable, as in the case (b) of the top panel of Fig. 1.In cases (i)-(iii), markers are averages over 2000 matrix samples with 25 realizations of the initial conditions each.In case (iv), we have used 20 matrix samples with 25 realizations each.In Fig. 6, we plot the spectra for a few matrix instances together with the theoretical boundary given by Eq. (C1).If r 2 is small enough, then the continuous part of the spectrum consists of two disconnected sets, which are centered around −µ 1 and −µ 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 λ 1 , i.e., the eigenvalue with the largest real part, we determine the real solutions of Eq. (C1).If where D = r 2 µ 2 1 q − 2r 2 µ 1 µ 2 q + r 2 µ 2 2 q + r 4 q 2 , (C4) then Eq. (C1) admits two real solutions, namely, 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 Solving Eq. (C2), we find that there may exist two outlier eigenvalues, namely, where Hence, the leading eigenvalue 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. ( 70) 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 326 < 0 and there is no outlier).We observe that the theoretical expression Eq. ( 70) 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 ⟨k γ ⟩ = ∞ for γ ≥ 2 as long as t < t ⋆ (N ), which is consistent with results on the leading eigenvalue of directed random graphs with power-law degree distributions obtained in [31].

FIG. 2 :
FIG.2: SN (t) = ⟨|y(t)| 2 ⟩ for the simulated dynamics of Eq.(2) on two examples of real-world graphs (circles for a food web with N = 180 and c = 10.7, and stars for a signaling network with N = 2288 and c = 13.3) are compared with the theoretical prediction S(t) (blue and red lines), given by Eq. (70), as well as with, e 2t(r−µ) and e −2µt (black and magenta lines).We have weighted the networks with couplings Jij that are i.i.d.random variables drawn from the distribution pJ (x) = (1/2)δ(x − 1) + (1/2)δ(x + 1) and we have used a single realization of the matrix J for both the food web and the signaling network.The diagonal elements are fixed to a constant −µ = −r + 0.27, such that the system is asymptotically unstable.Insets (a) and (b) show the spectra of the adjacency matrices A for the food web and the signaling network considered with random couplings.The red circle has radius r = √ c and is the predicted boundary of the spectrum according to(65).We have estimated SN (t) = ⟨|y(t)| 2 ⟩ by simulating the dynamics on the generated networks for 25 realizations of the initial condition y(0).

FIG. 3 :
FIG.3: Markers denote numerical results for SN (t) = ⟨|y(t)| 2 ⟩ for model(2) with matrix (A1) as a function of t and for different values of N .The Xij are i.i.d.random variables taken from a Gaussian distribution with zero mean and variance 1/N and we have set µ = −1.Dashed black line is the curve S(t) = e −2µt I0(2ρt) (see Eq. (11)).Markers are averages over 2000 matrix samples with 25 realizations of the initial conditions for each sample.

= 4 , r 2 = 64 FIG. 6 :
FIG.6: Spectra of adjacency matrices of weighted random graphs with the bimodal distribution (C1) on the diagonal, a Gaussian distribution pJ (x), and an in-degree (or out-degree distribution) p deg (k) = e −c c k /k!.We used the parameters c = 2, µ1 = 5, µ2 = 14, q = 0.5 and the values for J and r 2 = cJ 2 are provided below the figures.The red line denotes the theoretical value of the boundary of the absolute continuous spectrum, according to Eq. (C1) while the markers denote the eigenvalues of one single matrix instance of size N = 4000.