PAPER • OPEN ACCESS Conditional random matrix ensembles and the stability of dynamical systems To cite this article: Paul Kirk et al 2015 New J. Phys. 17 083025 View the article online for updates and enhancements. Related content Random antagonistic matrices Giovanni M Cicuta and Luca Guido Molinari - On the concentration of large deviations for fat tailed distributions, with application to financial data Mario Filiasi, Giacomo Livan, Matteo Marsili et al. - Induced Ginibre ensemble of random matrices and quantum operations Jonit Fischmann, Wojciech Bruzda, Boris A Khoruzhenko et al. - Recent citations How to deal with parameters for whole-cell modelling Ann C. Babtie and Michael P. H. Stumpf - Random antagonistic matrices Giovanni M Cicuta and Luca Guido Molinari - This content was downloaded from IP address 82.8.128.158 on 16/07/2018 at 12:13 New J. Phys. 17 (2015) 083025 doi:10.1088/1367-2630/17/8/083025 PAPER Conditional randommatrix ensembles and the stability of dynamical systems PaulKirk1,2, DelphineMYRolando2, AdamLMacLean andMichael PHStumpf Centre for Integrative Systems Biology andBioinformatics, Imperial College London, London SW7 2AZ,UK 1 New address:MRCBiostatistics Unit, Cambridge Institute of PublicHealth, CambridgeCB2 0SR,UK. 2 PK andDMYR contributed equally to this work. E-mail: paul.kirk@mrc-bsu.cam.ac.uk, delphine.rolando10@imperial.ac.uk, adam.maclean09@imperial.ac.uk and m.stumpf@imperial.ac.uk Keywords: stability, dynamical systems, control Supplementarymaterial for this article is available online Abstract Randommatrix theory (RMT) has found applications throughout physics and appliedmathematics, in subject areas as diverse as communications networks, population dynamics, neuroscience, and models of the banking system.Many of these analyses exploit elegant analytical results, particularly the circular law and its extensions. In order to apply these results, assumptionsmust bemade about the distribution ofmatrix elements. Herewe demonstrate that the choice ofmatrix distribution is crucial. In particular, adopting an unrealisticmatrix distribution for the sake of analytical tractability is liable to lead tomisleading conclusions.We focus on the application of RMT to the long-standing, and at times fractious, ‘diversity-stability debate’, which is concernedwith establishingwhether large complex systems are likely to be stable. Early work (and subsequent elaborations) brought RMT to bear on the debate bymodelling the entries of a system’s Jacobianmatrix as independent and identically distributed (i.i.d.) randomvariables. These analyses were successful in yielding general results that were not tied to any specific system, but relied upon a restrictive i.i.d. assumption.Other studies took an opposing approach, seeking to elucidate general principles of stability through the analysis of specific systems.Herewe develop a statistical framework that reconciles these two contrasting approaches.Weuse a range of illustrative dynamical systems examples to demonstrate that: (i) stability probability cannot be summarily deduced from any single property of the system (e.g. its diversity); and (ii) our assessment of stability depends on adequately capturing the details of the systems analysed. Failing to condition on the structure of dynamical systemswill skew our analysis and can, even for very small systems, result in an unnecessarily pessimistic diagnosis of their stability. 1. Introduction The notion of stability of stationary solutions is central to the study of dynamical systems. Formany applied questions, it is pivotal to know that the solutionwill return to the stationary values/orbits upon perturbation [1– 4]. Examples include, but are not limited to: ecology and population dynamics; [5, 6]; (celestial)mechanics; different areas of engineering; banking systems [7]; communication networks [4]; and neuroscience [8]. Different studies have shown that very diverse complex dynamical systems can bemodelled using similar mathematical tools [9, 10]. Formal aspects of stability have been studied extensively in the analysis of such complex systems, as well as in appliedmathematics. For linear time-invariant systems, the Routh–Hurwitz criterion sets out the conditions for global stability. More generally, the local stability of an equilibrium state of a non-linear ordinary differential equation (ODE) system can be assessed by inspecting the eigenvalues of the system’s Jacobianmatrix evaluated at the equilibrium [11]. If the real parts of the eigenvalues are all negative, then the equilibrium is (locally) stable. For any non- linear systems, only local stability is implied by a negative leading eigenvalue. Given our interest in typically non- OPEN ACCESS RECEIVED 20April 2015 REVISED 30 June 2015 ACCEPTED FOR PUBLICATION 1 July 2015 PUBLISHED 13August 2015 Content from this work may be used under the terms of theCreative CommonsAttribution 3.0 licence. Any further distribution of this workmustmaintain attribution to the author(s) and the title of thework, journal citation andDOI. © 2015 IOPPublishing Ltd andDeutsche PhysikalischeGesellschaft linear systems, we here consider the spectra of the Jacobians (or the sign of the largest eigenvalue, to be precise) directly, but keep inmind that the basin of stabilitymay befinite, and potentially confined to a small region. In ecology, where the analysis has perhaps been particularly divisive, early studies suggested that complexity —or ecological diversity—was key to stability [12, 13]. First Gardner andAshby [14] (who considered the problem fromawider perspective, with implications for examples as diverse as airport traffic and the human brain), and thenMay [5, 15], investigated theway inwhich stability changes as complexity increases. In these examples, complexity was defined in terms of the number of state variables (e.g. in the caseMay considers, the number of species) and the probability of an interaction between two variables. In order to generalise the analysis of the relationship between complexity and system stability,May avoided focusing on specific examples and instead considered ensembles of Jacobians defined in terms ofmatrix probability distributions. For suitably defined randommatrix ensembles (RME) [16–18] he showed that sufficiently large or complex systems have a probability of stability close to zero. Subsequent studies considered different RMEs designed to reflect a variety of features found in real systems, and have drawn different ormore nuanced conclusions regarding stability [19–21]. Other authors, especially from the ecological sphere, have pointed out the lack of realism of this approach [22] and estimated stability probabilities for specificODE systems, either through experiments [23, 24] or by sampling values for the system’s parameters and—for each sample—identifying the equilibriumpoints (EPs) and determining their stability [25–32]. By repeatedly sampling in this way,Monte Carlo estimates of stability probabilitiesmay be obtained. The advantage of such Monte Carlo approaches is that it is possible to condition on properties such as the feasibility of EPs (i.e. whether or not they are physicallymeaningful), which can again yield different conclusions regarding stability [25]. These approaches also define RMEs, but do so implicitly andwith reference to specificODEmodels. Given the variety of conclusions that have been drawn by different authors, the choice of RME is clearly crucial in determining the stability probability [4]. Randommatrices have also, of course, a distinguished track-record in different branches of physics. FollowingWigner’s earliest work on calculating fluctuations in the eigenvalue spectra ofHamiltonians describing atomic nuclei [33], they have found use in the analysis of awhole range offluctuations in different application areas: solid state physics, chemical reactions and transition state theory, and quantum chaos [34] are only some of the areas where they, in particular in the guise of theGaussianOrthogonal Ensemble and its generalisations, have come to use. But RMEs have also been found useful in puremathematics [35, 36], and, for example, the spectral properties of theGaussianUnitary Ensemble capture the statistical properties of the zeros of Riemann’s zeta function;more recently they have also been employed in cryptography. In all these applications RMEs are used to describe fluctuationswhich are believed to be separable from the secular dynamics of the underlying system.Here our use is subtly different. Instead of considering RMEs as general descriptors of some system—this has also been the strategy ofMay and, perhaps to a lesser extent, his followers—we are trying to condition the RMEon the properties of real systems that determinewhether or not the stationary states are stable or not. This then allows us to calculate a probability for a system to become unstable upon a small butfinite perturbation. So, rather thanmaking general statements about stability, our RMEs—whichwe refer to as conditional RMEs—are explicitly geared towards being used in specific contexts. While the success of traditional RMEs in capturing universal dynamics is based on assuming symmetries and homogeneity in thematrix entries, the stability analysis of specific real-world systems requires our conditional RMEs to exhibit the same heterogeneities that characterise real-world (i.e. problem-derived) Jacobianmatrices. Wewill showbelow that this is necessary in order to understandwhen andwhy a large dynamical system can be stable, but that this fully conditioned RME should not be used to draw general conclusions, as any rule from this system-specific approach can only highlight the behaviour of that systemor systemswith similar dynamics. 2. Stability and randommatrix ensembles For any particular parametricODEmodel, the Jacobianmatrix will usually exhibit structure and dependency between its entries, andwill typically be a function of themodel parameters and the state variables. The present work addresses the question of how assessments of stability changewhen the structure and dependency present in the Jacobian is properly taken into account. For example, for the Lorenz systemofODEs (see supplementary information for details), the Jacobian is given by, σ σ σ = − − − − − J r b x y z r z x y x b ( , , , , , ) 0 1 . ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ 2 New J. Phys. 17 (2015) 083025 PKirk et al As a consequence of this structure and dependency, and regardless of howwe choose the parameters of the system, only a particular family of n × n realmatrices will be obtainable as Jacobians of the system. For example, nomatter what valueswe take for the parameters of the Lorenz system, the (1, 3)-entry of the Jacobianmatrix will always be zero, and the (1, 1)-entry will always be equal to the negative of the (1, 2)-entry. It follows that if we were interested in assessing the stability probability of one of the Lorenz system’s EPs, it would be inappropriate to calculate ∣P h(stable )using, for example, amatrix probability density function, h, that associates non-zero density withmatrices for which the (1,3)-entry is non-zero or (1, 1)≠ −(1, 2). Nevertheless,many previous analyses have failed to account for the structure and dependency present in realistic Jacobianmatrices (i.e. ones derived frommodels of real systems), instead restricting attention tomatrix probability density functions that yield analytically tractable results and assuming that the results so obtainedwere general. 2.1. An illustration To further illustrate the implications of neglecting Jacobian structure, we consider a number of examples from a family ofODEswhose Jacobians have the form − −( )ab1 1 , with ∈ a b, . In this case, the space of allmatrices can be straightforwardly represented as a 2-dimensional Cartesian coordinate plane, inwhich the abscissa describes the value taken by a and the ordinate the value taken by b (as infigure 1). More precisely, we consider systems of the form, θ θ = − + = − + x t x g y y t y g x d d ( , ), d d ( , ), (1) 1 2 Figure 1. Stability of example equilibriumpoints (EPs) ofODEs.We consider different values for a and b and show as hatched areas (labelled ‘unstable regions’) the regions of the plane for which the resultingmatrix has an eigenvalue with non-negative real part. The non-hatched area corresponds to the stable region formatrices of the form − −( )ab1 1 We illustrate regions of the planewhich correspond to the Jacobians thatmay be obtained for the variousODE systems and equilibriumpoints considered in examples 1–3 (blue shaded area, and red, green, and purple lines, as indicated).We also represent using contours the randommatrix distribution that has traditionally been considered in the literaturewhen assessing stability probabilities. 3 New J. Phys. 17 (2015) 083025 PKirk et al whose Jacobians are given by, θ θ − ∂ ∂ ∂ ∂ − y g y x g x 1 ( , ) ( , ) 1 , 1 2 ⎛ ⎝ ⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟ where θ is the vector ofmodel parameters, and x and y are the state variables. Example 1.We start by considering the following (simple, linear) choices for g1 and g2: θ θ θ θ = = g y y g x x ( , ) , ( , ) , 1 1 2 2 with θ1 and θ2 both non-zero. In this case, the Jacobian for the system is θ θ − − 1 1 , 1 2 ⎛ ⎝⎜ ⎞ ⎠⎟ which is a function only of θ1 and θ2 (and not of x or y). The EPs are given by solving the simultaneous equations: θ θ = ⇒ − + = = ⇒ − + = x t x y y t y x d d 0 0, d d 0 0. (2) 1 2 It is straightforward to show that the only solution that holds for all values of θ θ,1 2 is =x y[ , ] [0, 0]. For this system, the formof the Jacobianmeans thatwemay obtain anymatrix of the form − −( )ab1 1 , provided that each of the parameters θ1 and θ2 can take any value in .We do note, however, that in practice thismodel is likely to be of only limited interest, since it describes a system inwhich both of the interacting variables (e.g. species) will eventually become extinct. Example 2.Wenext consider a nonlinear example: θ θ θ θ = = g y y g x x ( , ) , ( , ) , 1 1 2 2 2 with θ1 and θ2 both non-zero. In this case, the Jacobian for the system is θ θ − − y1 2 1 , 1 2 ⎛ ⎝⎜ ⎞ ⎠⎟ which is a function not only of θ1 and θ2, but also y. It is straightforward to show that the only EPs that exist for all permitted values of θ1 and θ2 are: (i) EP1: =x y[ , ] [0, 0]; and (ii) EP2: = θ θ θ θx y[ , ] , 1 1 1 2 2 1 2 ⎡ ⎣⎢ ⎤ ⎦⎥. The Jacobian evaluated at EP1 is θ − − 1 0 12 ⎜ ⎟⎛⎝ ⎞ ⎠. Thus the region of the (a, b) Cartesian coordinate plane representing the possible Jacobians associatedwith EP1 is simply the line a=0. Similarly, the Jacobian evaluated at EP2 is θ θ − − 1 2 1 2 2 ⎛ ⎝⎜ ⎞ ⎠⎟, and hence the region representing the possible Jacobians associatedwith EP2 is the line =b a2 . Example 3.Weconsider a further nonlinear example: θ θ θ θ = = g y y g x x ( , ) , ( , ) , 1 1 3 2 2 4 New J. Phys. 17 (2015) 083025 PKirk et al with θ1 and θ2 both non-zero. In this case, the Jacobian for the system is θ θ − − y1 3 1 ,1 2 2 ⎛ ⎝⎜ ⎞ ⎠⎟ which is again a function of θ θ,1 2, and y. In this case, it is again straightforward to show that the only EPs that exist for all permitted values of θ1 and θ2 are: (i) EP1: =x y[ , ] [0, 0]; and (ii) EPs 2 and 3: = ± θ θ θ θ θ x y[ , ] , 11 1 3 2 3 1 2 ⎡ ⎣⎢ ⎤ ⎦⎥. The Jacobian evaluated at EP1 is again θ − − 1 0 12 ⎜ ⎟⎛⎝ ⎞ ⎠. Thus the region of the (a, b) Cartesian coordinate plane representing the possible Jacobians associatedwith EP1 is again the line a=0. The Jacobian evaluated at EP 2 or 3 is θ θ − − 1 3 1 2 2 ⎛ ⎝⎜ ⎞ ⎠⎟, and hence the region representing the possible Jacobians associatedwith these EPs is the line =b a3 . Assessing stability for these examples For anymatrix of the form = − −( )J ab1 1 , the characteristic equation λ∣ − ∣ =J I 02 may be expanded as λ λ− − − − − =ab( 1 )( 1 ) 0, i.e. λ + − =ab( 1) 0.2 The eigenvalues of J are the solutions of this equation, and are given by λ = − ± ab1 .1,2 J is in the stable region, ΛS2, provided the real parts of λ1,2 are both negative. First, we note that if ab is negative (i.e. if = −a bsgn( ) sgn( )) then the real part of both eigenvalues is−1 and hence J is in the stable region. If ab is positive, then λ = − − < −ab1 12 , so this eigenvalue is certainly negative, and it remains only to consider the sign of the other eigenvalue, λ = − + ab11 . This eigenvalue is negative if and only if ab 1.We take various choices for μ andΣ, and illustrate the resulting bivariate normal distributions using coloured contours. (A) The location of themean has an impact on stability probability: (I) represents the usual choice, μ = ⊤[0,0] ; however other choices can clearly lead to (II) lower or (III) higher stability probabilities. (B) The variances of a and bhave an impact on stability probability: e.g. forfixedmean μ = ⊤[0,0] , taking smaller or larger variances leads to, respectively, (I) higher or (II) lower stability probabilities. (C) The covariance between a and b has an impact on stability: e.g. forfixedmean μ = ⊤[0,0] , whether a and b covary negatively or positively leads to, respectively, (I) higher or (II) lower stability probabilities. 6 New J. Phys. 17 (2015) 083025 PKirk et al ∑ Λ∣ ≈ ∈θ =  ( )P N J ˆ (stable ) 1 , (4) i N S n 1 * i( ) where  X( ) is the indicator function, which is 1 ifX is true and 0 otherwise. 2.4. Conditional randommatrix ensembles The estimated stability probability defined in the previous section is the probability of a systembeing stable, conditional on a given system architecture; as discussed in section 2 and illustrated in section 2.1 such architectures do not arise without a concrete context. However, the conditions for which the circular law is believed to hold lack this connection to reality, at least formesoscopic systems. To study the effects of this context, which is encapsulated by the statistical properties of, and dependencies among, the entries in the Jacobianmatrices …θ θJ J, ,* * N(1) ( ) , we consider two further randommatrix distributions, constructed by permutation of the entries of our original RME. First, we form a newmatrix ensemble, …K K, , N(1)* ( )* , inwhich the dependency between entries is broken. For each ℓ ∈ … N{1, , } and ∈ … × …i j p p( , ) {1, , } {1, , }we set = θℓ ( )( )K Jij ij( )* * q( ) , with q drawnuniformly at random from …N{1, }. In this way, themarginal distribution of the ij-entries across the ensemble ofK*matrices is the same as themarginal distribution of ij-entries across the ensemble of θJ *matrices.Maintaining themarginal distributions ensures that the dependency between entries is the only quantity thatwe are altering: in particular, the location of zeros in thematrix and themagnitudes of interaction strengths aremaintained.We construct a further RME, …L L, , N(1)* ( )* , where for each ℓ ∈ … N{1, , } and ∈ … × …i j p p( , ) {1, , } {1, , }, we set = θℓ ( )( )L Jij rs( )* * q( ) , with q drawn uniformly at random from …N{1, }, and r and s (independently) drawn uniformly at random from …p{1, }. Now, the location of zeros in thematrix is no longer fixed; although the probability of an entry being zero is the same for the L*matrices as for the θJ *ʼs andK*ʼs.Moreover, each entry of the L*matrices is i.i.d.We henceforth refer to the θJ *matrices as the fully conditioned system (FCS) ensemble (most structure); theK*matrices as the independent ensemble (intermediate structure); and the L*matrices as the i.i.d. ensemble (least structure).We illustrate the properties of these three RMEs, and themethods for their construction, infigure 3. To further our investigationwe defined fourmore RMEs (which are presented inmore detail in the supplementary information). Thefirst onewill be referred to as the independent normal ensemble. It is constructed as follows: For each (i, j), wefit an independent normal distribution to the ij-entries of the sampled Jacobians, …θ θJ J, ,* * N(1) ( ) . That is, for each (i, j), we calculate themean, ∑μ = θ = ( )N J 1 ,i j q N ij( , ) ind 1 * q( ) and standard deviation, σ = θ = ( )Js.d. .i j ij q N ( , ) ind * 1 q( ) ⎧⎨⎩ ⎫⎬⎭ We then construct the newRME, …M M, , N(1)ind ( )ind , where for each ℓ ∈ … N{1, , } and ∈ … × …i j p p( , ) {1, , } {1, , }, we set ℓ( )M ij( ) ind to be a sample drawn from the univariate normal distribution withmean μ i j( , ) ind and standard deviation σ i j( , )ind. By construction, themean and standard deviation of the ij-entries across the ensemble ofMindmatrices are the same as themean and standard deviation of the ij-entries across the ensemble of θJ *matrices (the FCS ensemble) and across the ensemble ofK*matrices (the independent ensemble). A further ensemble is given by the independent Pearson ensemble. As in the independent normal case defined above, this newRME is defined by fitting a distribution to the ij entries of the sampled Jacobians, except that rather than using a normal distribution and just capturing themean and standard deviation, we also capture the skewness and kurtosis of the ij-entries of the θJ *matrices. That is, in addition to μ i j( , ) ind and σ i j( , )ind defined earlier, we also calculate skewness γ = θ = ( )Jskewness .i j ij q N ( , ) ind * 1 q( ) ⎧⎨⎩ ⎫⎬⎭ and kurtosis, κ = θ = ( )Jkurtosis .i j ij q N ( , ) ind * 1 q( ) ⎧⎨⎩ ⎫⎬⎭ 7 New J. Phys. 17 (2015) 083025 PKirk et al We then construct anRME, …M M, , N(1)pear ( )pear, where for each ℓ ∈ … N{1, , } and ∈ … ×i j p( , ) {1, , } … p{1, , }, we set ℓ( )M ij ( ) pear to be a sample drawn from aunivariate Pearson distributionwithmean μ i j( , ) ind, standard deviation σ i j( , )ind, skewness γ i j( , ), and kurtosis κ i j( , ). This RME thus sharesmany of the properties of the marginal distributions of ij-entries across the ensemble of θJ *matrices, but does not capture the dependencies between them. The third additional RMEwill be referred to as the i.i.d. normal ensemble. This timewewill notfit a distribution to the ij-entries of the θJ *matrices, but insteadwefit a normal distribution, using the same technique thatwe used for the independent normal ensemble defined above, to the ij-entries of the L*matrices (i.e. those from the i.i.d. ensemble ). Finally, we construct anRME that attempts to capture some of the dependencies between the entries of the θJ *matrices.We define c(M) to be the vector obtained by concatenating the columns of thematrixM (and further define −c 1be the inverse operation, so that, for example, =−c c M M( ( ))1 ). Applying c ( · ) to the matrices fromour FCSRME,we obtainN vectors of length p × p, namely: …θ θc J c J( ), , ( ).* * N(1) ( ) To these, wefit (bymaximum likelihood) a ×p p( )-variate normal distribution.We then sampleN vectors, …v v, , N1 , of length ×p p from this distribution, and form anew ensemble …M M, , N(1)mvn ( )mvn by setting = −M c v( )q q( )mvn 1 .We will call this new ensemble themultivariate normal ensemble. These new ensembles allow us to control which aspect of the structure of the FCS gives it its stability properties. For instance comparing the independent normal ensemble, the independent Pearson ensemble and the independent ensemblewe can show the impact of the differentmoments of the distribution. Themultivariate normal and FCS ensembles can be used for the same purpose in the case where dependencies are considered. More detail about the different RMEs is provided in the supplementary information. Figure 3.Randommatrix ensembles (RMEs). (A) For a givenODE system (e.g. the Lorenz equations) and equilibriumpoint, specifying a distribution for the parameters defines a random Jacobianmatrix distribution. (B) Samples from this distribution define the FCSmatrix distribution; the independent and i.i.d. distributions are obtained from this by permuting elements as illustrated. jkl m( ) is the term in row k and column l obtained in themth sample from the random Jacobianmatrix distribution. Themarginal distributions for the elements of thematrices in the three distributions. In the FCS case, these reflect the parameter distributions and the expressions for the Jacobian entries presented in (A); by construction, themarginals in the independent case are the same as for the FCS; while in the i.i.d. case, all entries have the samemarginal distribution. (D)We illustrate the joint distribution for twomatrix entries: in the FCS case, the two entries exhibit dependency, whereas in the independent and i.i.d. cases, the joint is the product of the marginals. 8 New J. Phys. 17 (2015) 083025 PKirk et al 3. Results 3.1. RME choice determines stability assessment Stability, as stressed above and in the literature, is an issue in awide variety of domains, and therefore we consider a set of systems that cover different qualitative behaviour of dynamical systems. The fourODEmodels that we consider have in common that the EPs and Jacobians can be identified analytically, whichmakes analysis straightforward; they are: (i) the Lorenz system [38]; (ii) amodel of the cell cycle due to Tyson [39]; (iii) amodel of viral dynamics due toNowak andBangham [40]; and (iv) an SEIR (susceptible-exposed-infective-recovered) population dynamicsmodel [41]. For brevity, wewill refer to these four examplemodels as: (i) ‘Lorenz’; (ii) ‘Tyson’; (iii) ‘N&B’; and (iv) ‘SEIR’. In each case, we present results for physically or biologically feasible EPs and generate 100 000matrices fromourRMEs in order to obtainMonteCarlo estimates of stability probabilities. Full details of thesemodels and their corresponding RMEs are provided in the supplementary information. Figure 4(A) shows the eigenspectra for the FCS, independent, and i.i.d. RME regimes.While the i.i.d. eigenspectra are broadly circular, we observe diverse and decidedly non-circular shapes for the other two cases, highlighting the limitations of previous analyses based upon the circular law. Figure 4(B) shows the density of eigenvalues in the complex plane for the differentmodels andRMEs. The eigenspectra distributions are typically much less dispersed for the FCS ensemble than for the other two. As shown infigure 4(C) this also leads to systematic differences in the real parts of the leading eigenvalues, which determine stability. Table 1 shows that, whatever themodel considered, the probability of stability of the i.i.d. ensemble is always very different from the stability probability of the FCS. The RMEs that includemore of the structure of the A B C Lo re nz SE IR N & B Ty so n 15 10 5 0 -5 -10 -15 -15 -10 -5 0 5 10 15 -15 -10 -10 -10 -10 -10 -5 0 0 0 0 5 10 10 1010 15 2 1.5 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 0.5 -1 -1 -1 -1 -1 -1 10-1.5 -1.5 1.5 -2 -2 2 0 0 0 0 0 10 8 6 4 2 0 -2 -4 -6 -5 50 -8 -10 10 8 6 4 2 0 -2 -4 -6 -8 -10 -10 10 -5 50-10 10 Real part FCS independent i.i.d. FCS independent i.i.d. FCSindependenti.i.d Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt Im ag in ar y pa rt 600 500 400 300 200 100 0 5000 4000 3000 2000 1000 0 -2.5 -2 -2 -2 -1.5 -0.5 0.5 0 0 0 0 -1 -1 -1 -1 1 1 1 1 1.5 2 2 2 2 2 2.5 -2 1000 800 600 400 200 0 -3 -2 -2 -2 -2 -1 0 0 0 0 1 2 2 2 2 3 600 500 400 300 200 100 0 Real part Real part Real part Maximal real part of eigenvalue Stability probabilities: Stability probabilities: Stability probabilities: Stability probabilities: 0.32109 0.03209 0.11151 1 0.56483 0.1371 1 0.57161 0.12173 0.99916 0.84059 0.10348 R el at iv e fre qu en cy R el at iv e fre qu en cy R el at iv e fre qu en cy R el at iv e fre qu en cy 0.25 0.2 0.15 0.1 0.05 0 -5 0 5 10 15 20 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0-2 -2 -1 0 1 2 3 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0 4 6 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0-5 0 5 10 Figure 4. Stability results for the four examplemodels. (A) The eigenspectra for eachmodel and randommatrix distribution, shown as scatterplots. (B) The eigenvalue distributions visualised using heatmaps (to aid visualisation, we omit pure imaginary eigenvalues). (C) The distributions ofmaximal eigenvalues together with the estimated stability probabilities. 9 New J. Phys. 17 (2015) 083025 PKirk et al system in their construction, yield stability probabilities closer to those of the FCS. Inmost cases, the univariate Pearson ensemble had a probability stability closer to that of the FCS than the univariate normal, showing that consideringmoremoments when building the RME improves the estimation of the stability of the system. In summary,Monte Carlo estimates for the stability probabilities decrease aswe decrease the amount of realism captured by theRME. Thus, failing to condition on the real-world heterogeneity and dependency present in the Jacobian can result in an unnecessarily pessimistic assessment of stability, even for these small systems. Considering RMEs inwhichwe tightly control themean, standard deviation, skewness and kurtosis of themarginal distributions of Jacobian entries, demonstrates that all of these properties also have an impact on stability. 3.2. Large dynamical systems can be stable To illustrate the effects of inadequately capturingmodel structure and parameter dependencies on the stability probabilities of larger systems, we consider extensions to the SEIRmodel (figure 5(A)).We allow formultiple subpopulations of exposed (E) individuals (in the supplementary informationwe also investigate extensions with heterogeneous infective subpopulations), enabling us to control system size.We again consider the three RMEs described above (see supplementary information for full details). Aswe increase the size of the system, the probability of stability remains 1 in the FCS case, but rapidly diminishes in the i.i.d. case (figure 5(B)). The independent case is intermediate, indicating that not only can the dependence betweenmatrix entries be important, but also their heterogeneity. Heterogeneity changes the location of the centre of thematrix p.d.f. h, and also how it stretches in different directions, whichmodifies the proportion of probabilitymass falling in the stable region, Λ ⊂ M ( )Sn n . Figure 5(C) shows how summaries of the distributions of leading eigenvalues change aswe increase the number of exposed populations, with figure 5(D) providing corresponding density estimates for a selection of these numbers. In the i.i.d. case, the distributions andmedian values shift away from stable negative values toward unstable positive values. This is in stark contrast to the FCS casewhere, regardless of the number of exposed populations, the distribution of leading eigenvalues only has support on the negative real line (and hencewe always have stability probability 1).Moreover, as we increase the number of exposed populations, the median value of the leading eigenvalue tends to becomemore negative. The independent case is again intermediate, with themedian value staying relatively constant. Figures S4 and S5 in the supplementarymaterial, wherewe consider the effect of the i.i.d. normal, the independent normal, the independent Pearson, and themultivariate normal ensembles on thesemodels, bring more evidence to our previous observations. Aswewould expect, themore of the underlying system’s structure thatwe capture using our chosenRME, the closer we get to the stability probability estimated using the FCS ensemble. The i.i.d. normalRME,which does not include anymore structure than the i.i.d. ensemble, leads to similar stability probabilities to those obtained using that RME. The independent normal, which allows the heterogeneity of variances andmeans found in the real system to be described, yields stability probabilities closer to the FCS ensemble. The independent Pearson, which also takes into account information about the skewness and kurtosis of each entry, gets even closer to the FCS ensemble, and is very similar to the independent case. Finally, themultivariate normalRME (which allows us tomodel–in a simplistic fashion–some of the dependencies between the entries of the Jacobian) results in stability probabilities that are always closer to the FCS ensemble than those obtained using the independent normalRME.We found that the stability probabilities obtained using themultivariate normalRME are not consistentlymore different fromFCS stability probabilities than those obtained using the independent PearsonRME.Thus, accounting for dependency between the Jacobian’s entriesmay, depending on the problem at hand, bemore or less important than accounting for higher order properties of themarginal distributions of the entries. Table 1.Estimated stability probabilities for the four examplemodels using our seven different RME regimes. Tyson SEIR N&B Lorenz FCS 0.32109 1 1 1 MultivariateNormal 0.1105 0.72188 0.43744 0.92382 Independent 0.11151 0.56483 0.57161 0.7225 Univariate Pearson 0 0.58252 0.93296 0.43542 UnivariateNormal 0.116 0.53502 0.26504 0.48366 I.i.dNormal 0.0377 0.14904 0.12702 0.14602 I.i.d. 0.03209 0.1371 0.12173 0.10419 10 New J. Phys. 17 (2015) 083025 PKirk et al 4.Discussion The stability of real-world systems—which is nearly ubiquitously observed [1–3]—might seemperplexing in light of classical results in randommatrix theory. By considering how randommatrices can bemade to reflect the properties of the Jacobianmatrix of real dynamical systems it becomes possible to resolve and reconcile apparent contradictions in the literature [5, 30, 31]. Figure 5.How stability changeswith system size depends on the randommatrix distribution. (A)An extension of the SEIRmodel in whichwemodel n exposed populations, …E E, , n1 . (B) Plot showing for each of the randommatrix distributions how estimated stability probability changes as we increase the number of exposed populations. Bars denote ±2 s.d.Monte Carlo error bars. (C) Plot showingmedian (filled circle) and interquartile range (bars) for the distributions of leading eigenvalues. (D)Density estimates for the distributions of leading eigenvalues. 11 New J. Phys. 17 (2015) 083025 PKirk et al In agreement with previous authors, our results demonstrate that stability probability can be affected by the mean and standard deviation of the entries in the Jacobian, as well as the dependencies between them [30, 31], and further show that properties such as the skewness and kurtosis of the entries can also have an impact. This is unsurprising: as is clear from equation (2) and illustrated infigure 2, RMEswith different properties will result in different proportions of probabilitymass fallingwithin the stable region, Λ ⊂ M ( )Sn n . Reported stability probabilities should therefore always explicitly acknowledge that they are conditioned on a particular choice of RME,which has to be carefully justified. WhileMay’smathematical study clearly shows that in some circumstances an increase in complexity can lead to instability [5], Haydon’s study highlights cases where complexity, in the shape of strong and numerous interconnections, is necessary to get higher stability [19]. Other examples where complex systems have been shown to be stable can be found in the literature, in particular Kokkoris et al show that different variances of the interaction strengthwill allow for different levels of complexity of the systemwhile keeping the same stability [30].However, none of these results can be generalised lightly. They in fact show that different systems are impacted differently by changes in complexity, and that no general prediction can bemade. Here, the FCS ensemble conditions on themodel structure, so that the RME is defined through the distribution ofmodel parameters. In our case, these distributions have been chosen by selecting plausible or interesting ranges for the parameters, and taking uniformdistributions within these ranges (in the supplementary informationwe consider alternative possibilities). In real applications, a natural choice for the distribution ofmodel parameters would be provided in the Bayesian formalismby the posterior parameter distribution (or, if no data are available, by the prior). From this, wemay obtain the posterior (or prior) predictive distributions [42] of the leading eigenvalues, fromwhichwemay derive the probability of stability. In this way, a truly realistic assessment of stabilitymay be obtained for the system, inwhichwe have conditioned on both our current understanding of the system architecture (encapsulated in ourmathematicalmodel) and our current uncertainty in themodel’s parameters. Through our analyses, we have demonstrated that identifying any single property of the RME as being the general determinant of stability ismisleading, except in some cases when the systemhas been very strictly defined [20, 27, 30, 31]. Stability probability is determined by the topology of the stable region, and howmuch probabilitymass is depositedwithin that region by theRME. This cannot be summarily deduced from any single property of the RME. At this stage it seems that system stability is system specific and that little can be gleaned fromgeneral approaches that will at best be uninformative if not entirelymisleading. In particular, we cannot assess the probability of a systembeing stable based only on its size, diversity or complexity. It is especially important to keep this inmind as the stability ofmore complicated systems is considered (see e.g. [7, 43]). This does not rule out the possibility that there are sets of rules or principles that could greatly shift the balance in favour of stability. Negative feedback, for example, is likely to lead tomore stable behaviour (even for stochastic systems). It is in principle possible to condition on such local structures (in terms of correlated entries in the Jacobian) thatmay confer or contribute to overall stability of a system [44]. To apply such knowledge to real systems, however, would require a level of certainty about the underlyingmechanisms thatwe currently lack for all but themost basic examples. Nevertheless, even in the presence of uncertainty about system structures, local negative feedback between species, for example, would tend to favour stability, whereas positive feedback (ormerely the lack of negative feedback)would typically result in amplification of initially small perturbations to a system’s behaviour [45]. Appendix.Methods Twomainmethodswere used. Thefirst used an analytical approach, while the second used a numerical approach. Thefirstmethodwas used onmodels 1 to 4, and onmodels 5 and 6 for n=1–6. In these cases the number of samples usedwas 100 000. The secondmethodwas used onmodel 5 and 6 for n=1–10 and n=20, 30,…,100, this time, for computational reasons, the sample size was 1000. A.1. The analyticalmethod For each of themodels, the EPswere found usingMatlab’s analytic equation solver, solve. The Jacobianwas also described analytically usingMatlab’s jacobian function. The different parameters were sampled from a uniformdistribution usingMatlab’s rand function. The choice of the range of the parameters will be described in detail below. 12 New J. Phys. 17 (2015) 083025 PKirk et al A.1.1. Algorithm. For each of themodel, the EPs and their stability for a different range of parameters are evaluated in the followingway. (i) Define the systemofODEs. (ii) Solve = ∀ ∈ …x i p˙ 0, {1, , }i , where p is the number of species in the system and ∈ …x i p, {1, , }i is the set of species in our system. (iii) Compute the Jacobianmatrix of our system. (iv) Sample a set of parameters for the system. (v) Evaluate the EPs for that set of parameters. (vi) If the EPs are biologically realistic, i.e. all of the species have a concentration that is positive or null, then evaluate the Jacobianmatrix at those EPs. (vii) Else, ‘reject’ that set of parameters and sample a new set. (viii) Reproduce steps 5–7 until the number of samples accepted reaches the number of samples wanted. (ix) For each Jacobian matrix obtained, compute the eigenvalues; consider their maximum real part, each system is considered as stable if and only if thatmaximum real part is strictly negative. (x) The probability of a system being stable is then the number of stable systems divided by the number of samples. In order to evaluate the stability under the independent and i.i.d. conditions we add a step between steps 8 and 9: we process the Jacobianmatrices by doing some permutations of the entries as described in section 2.4. A.2. The numericalmethod For each of themodels, the EPswere found usingMatlab’s equation solver, solve. The Jacobianwas described analytically usingMatlab’s jacobian function. The different parameters were sampled from a uniform distribution usingMatlab’s rand function. The choice of the range of the parameters will be described in detail below. Thismethodwas used only on the S(nE)IR and the SE(nI)Rmodels for computational reason (the numericalmethod becoming too expensive with big n). Thismethodworkswell in this case because for these models we know that there are 2 EPs and they are easily identified, sowe can easily segregate the cases corresponding to each of themwhen computing the probability of stability. Thefirst EP, where thewhole population is composed of recovered individuals, is not interesting because it is obviously stable, sowe only consider the other EP in each system. A.2.1. Algorithm. For each of themodel, the EPs and their stability for a different range of parameters are evaluated in the followingway. (i) Define the systemofODEs. (ii) Compute the Jacobianmatrix of our system analytically. (iii) Sample a set of parameters. (iv) Solve = ∀ ∈ …x i p˙ 0, {1, , }i , where p is the number of species in the system and ∈ …x i p, {1, , }i is the set of species in our system, using the set of parameters sampled. (v) Only keep the EP which is not a population fully composed of recovered individuals and no other type of individuals. (vi) If the EPs are biologically realistic, i.e. all of the species have a concentration that is positive or null, then evaluate the Jacobianmatrix at those EPs. (vii) Else, ‘reject’ that set of parameters and sample a new set. (viii) Reproduce steps 3–7 until the number of samples accepted reaches the number of samples wanted. (ix) For each Jacobian matrix obtained, compute the eigenvalues; consider their maximum real part, each system is considered as stable if and only if thatmaximum real part is strictly negative. 13 New J. Phys. 17 (2015) 083025 PKirk et al (x) The probability of a system being stable is then the number of stable systems divided by the number of samples. In order to evaluate the stability under the independent and i.i.d. conditions we add a step between steps 8 and 9: we process the Jacobianmatrices by doing some permutations of the entries as described in section 2.4. A.3. Parameter ranges Two criteria were used to choose the range of the parameters: first they had to be realistic; second the range had to be small enough to allow a thorough sampling of the space obtained to be computationally tractable. In the cases when the ranges were fixed arbitrarily, we verified that choosing different ranges would not impact the qualitative results. A.3.1.Model 1: the Lorenz system. To get the results obtained in themain article, we sampled the parameters from the following ranges: β ρ σ ∈ ∈ ∈ [0, 10], [1, 11], [0, 10]. ρ is considered to be bigger or equal to 1 in order to ensuremore than just the origin as a equilibriumpoint, whichmakes our studymore interesting. A.3.2.Model 2: amodel of the cell division cycle. All the parameters were taken to be uniformly distributed between 0 and 1. A.3.3.Model 3: theNowak and Banghammodel. All the parameters were taken to be uniformly distributed between 0 and 1. A.3.4.Models 4, 5 and 6: the SEIR and extended SEIRmodels. Model 4: the SEIRmodel. All the parameters were taken to be uniformly distributed between 0 and 1. Model 5: the S(nE)IRmodel. All the parameters were taken to be uniformly distributed between 0 and 1. Model 6: the SE(nI)Rmodel. All the parameters were taken to be uniformly distributed between 0 and 1. References [1] Jirsa V andDingM2004Will a large complex systemwith time delays be stable?Phys. Rev. Lett. 93 070602 [2] TimmeM,Geisel T andWolf F 2006 Speed of synchronization in complex networks of neural oscillators: analytic results based on randommatrix theoryChaos 16 015108 [3] Kriener B 2012How synaptic weights determine stability of synchrony in networks of pulse-coupled excitatory and inhibitory oscillatorsChaos 22 033143 [4] Couillet R andDebbahM2013 Signal processing in large systems: a new paradigm IEEE Signal Process.Mag. 30 24–39 [5] MayRM1972Will a large complex system be stable?Nature 238 413–4 [6] McCannKS 2000The diversity-stability debateNature 405 228–33 [7] Haldane AGandMayRM2011 Systemic risk in banking ecosystemsNature 469 351–5 [8] Feng J, Jirsa VK andDingM2006 Synchronization in networks with random interactions: theory and applicationsChaos: Interdiscip. J. Nonlinear Sci. 16 015109 [9] Boccaletti S, LatoraV,Moreno Y,ChavezMandHwangD 2006Complex networks: structure and dynamics Phys. Rep. 424 175–308 [10] Carlson JMandDoyle J 2002Complexity and robustnessProc. Natl Acad. Sci. USA 99 (Suppl. 1) 2538 [11] Strogatz SH 2001NonlinearDynamics AndChaos:With Applications To Physics, Biology, Chemistry And Engineering (Studies in Nonlinearity) (Reading,MA: Addison-Wesley) [12] EltonC 1958The Ecology of Invasions by Animals and Plants (London: Chapman andHall) [13] MacArthur R 1955 Fluctuations of animal populations and ameasure of community stability Ecology 36 533 [14] GardnerMRandAshbyWR1970Connectance of large dynamic (cybernetic) systems: critical values for stabilityNature 228 784 [15] MayR 1973 Stability andComplexity inModel Ecosystems (Princeton,NJ: PrincetonUniversity Press) [16] EdelmanA andRaoNR1999Randommatrix theoryActaNumer. 14 233–97 [17] MehtaML 2004RandomMatrices (NewYork: Academic) [18] Tao T, VuV andKrishnapurM2010Randommatrices: universality of esds and the circular lawAnn. Probab. 38 2023–65 [19] HaydonDT 2000Maximally stablemodel ecosystems can be highly connected Ecology 81 2631–6 [20] Neutel AM2002 Stability in real foodwebs: weak links in long loops Science 296 1120–3 [21] Allesina S andTang S 2012 Stability criteria for complex ecosystemsNature 483 205–8 [22] Lawlor LR 1978A comment on randomly constructedmodel ecosystemsAm.Nat. 112 445–7 [23] TilmanD andDowning J A 1994 Biodiversity and stability in grasslandsNature 367 363–5 [24] Givnish T J 1994Does diversity beget stability?Nature 371 113–4 [25] Roberts A 1974The stability of a feasible random ecosystemNature 251 607–8 [26] GilpinME1975 Stability of feasible predator-prey systemsNature 254 137–9 14 New J. Phys. 17 (2015) 083025 PKirk et al [27] King AWandPimmSL 1983Complexity, diversity, and stability—a reconciliation of theoretical and empirical resultsAm.Nat. 122 229–39 [28] PimmS1984The complexity and stability of ecosystemsNature 307 321–6 [29] McCannK,Hastings A andHuxel GR 1998Weak trophic interactions and the balance of natureNature 395 794–8 [30] Kokkoris GD, JansenVAA, LoreauMandTroumbis AY 2002Variability in interaction strength and implications for biodiversity J. Animal Ecol. 71 362–71 [31] JansenVAA andKokkoris GD2003Complexity and stability revisitedEcol. Lett. 6 498–502 [32] ChristianouMandKokkoris GD2008Complexity does not affect stability in feasiblemodel communities J. Theor. Biol. 253 162–9 [33] Wigner E P 1959 Statistical properties of real symmetricmatrices withmany dimensionsProc. 4thCanadianMathematical Congress ed MSMacPhail (Toronto: University of Toronto Press) pp 174–84 [34] Andreev AV et al 1996Quantum chaos, irreversible classical dynamics, and randommatrix theory Phys. Rev. Lett. 76 3947–50 [35] Tao T 2012Topics in RandomMatrix Theory (Providence, RI: AmericanMathematical Society) [36] SchroederMR2009Number Theory in Science andCommunication 5th edn (Berlin: Springer) [37] EvansMR et al 2013Do simplemodels lead to generality in ecology?Trends Ecol. Evol. 28 578–83 [38] Lorenz EN1963Deterministic nonperiodic flow J. Atmos. Sci. 20 130–48 [39] Tyson J J 1991Modeling the cell division cycle: cdc2 and cyclin interactions Proc. Natl Acad. Sci. USA 88 7328–32 [40] NowakMAandBanghamC1996 Population dynamics of immune responses to persistent viruses Science, NY 272 74–79 [41] Aron J L and Schwartz I B 1984 Seasonality and period-doubling bifurcations in an epidemicmodel J. Theor. Biol. 110 665–79 [42] GelmanA,Carlin J B, SternHS andRubinDB 2003BayesianData Analysis 2nd edn (London/Boca Raton, FL: Chapman and Hall/CRC) [43] Wilkinson RD2013Approximate Bayesian computation (ABC) gives exact results under the assumption ofmodel error Stat. Appl. Genet.Mol. Biol. 12 129–41 [44] Thorne TWand StumpfMPH2007Generating confidence intervals on biological networksBMCBioinform. 8 467 [45] CosentinoC andBates D 2011 Feedback Control in Systems Biology (Boca Raton, FL: CRCPress) 15 New J. Phys. 17 (2015) 083025 PKirk et al