Conditional random matrix ensembles and the stability of dynamical systems

To further illustrate the implications of neglecting Jacobian structure, we consider a number of examples from a family of ODEs whose Jacobians have the form $\left(\begin{array}{cc}-1 & a\\ b & -1\end{array}\right)$, with $a,b\in {\mathbb{R}}$. In this case, the space of all matrices can be straightforwardly represented as a 2-dimensional Cartesian coordinate plane, in which the abscissa describes the value taken by a and the ordinate the value taken by b (as in figure 1).

Zoom In Zoom Out Reset image size

More precisely, we consider systems of the form,

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}x}{{\rm{d}}t} & = & -x+{g}_{1}(y,\theta ),\\ \displaystyle \frac{{\rm{d}}y}{{\rm{d}}t} & = & -y+{g}_{2}(x,\theta ),\end{eqnarray} \tag{ 1 }$$

whose Jacobians are given by,

$$\begin{eqnarray*}\left(\begin{array}{cc}-1 & \displaystyle \frac{\partial }{\partial y}{g}_{1}(y,\theta )\\ \displaystyle \frac{\partial }{\partial x}{g}_{2}(x,\theta ) & -1\end{array}\right),\end{eqnarray*}$$

where θ is the vector of model parameters, and x and y are the state variables.

Example 1. We start by considering the following (simple, linear) choices for g₁ and g₂:

$$\begin{eqnarray*}{g}_{1}(y,\theta ) & = & {\theta }_{1}y,\\ {g}_{2}(x,\theta ) & = & {\theta }_{2}x,\end{eqnarray*}$$

with ${\theta }_{1}$ and ${\theta }_{2}$ both non-zero. In this case, the Jacobian for the system is

$$\begin{eqnarray*}\left(\begin{array}{cc}-1 & {\theta }_{1}\\ {\theta }_{2} & -1\end{array}\right),\end{eqnarray*}$$

which is a function only of ${\theta }_{1}$ and ${\theta }_{2}$ (and not of x or y).

The EPs are given by solving the simultaneous equations:

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}x}{{\rm{d}}t} & = & 0\Rightarrow -x+{\theta }_{1}y=0,\\ \displaystyle \frac{{\rm{d}}y}{{\rm{d}}t} & = & 0\Rightarrow -y+{\theta }_{2}x=0.\end{eqnarray} \tag{ 2 }$$

It is straightforward to show that the only solution that holds for all values of ${\theta }_{1},{\theta }_{2}$ is $[x,y]=[0,0]$. For this system, the form of the Jacobian means that we may obtain any matrix of the form $\left(\begin{array}{cc}-1 & a\\ b & -1\end{array}\right)$, provided that each of the parameters ${\theta }_{1}$ and ${\theta }_{2}$ can take any value in ${\mathbb{R}}$. We do note, however, that in practice this model is likely to be of only limited interest, since it describes a system in which both of the interacting variables (e.g. species) will eventually become extinct.

Example 2. We next consider a nonlinear example:

$$\begin{eqnarray*}{g}_{1}(y,\theta ) & = & {\theta }_{1}{y}^{2},\\ {g}_{2}(x,\theta ) & = & {\theta }_{2}x,\end{eqnarray*}$$

with ${\theta }_{1}$ and ${\theta }_{2}$ both non-zero. In this case, the Jacobian for the system is

$$\begin{eqnarray*}\left(\begin{array}{cc}-1 & 2{\theta }_{1}y\\ {\theta }_{2} & -1\end{array}\right),\end{eqnarray*}$$

which is a function not only of ${\theta }_{1}$ and ${\theta }_{2}$, but also y.

It is straightforward to show that the only EPs that exist for all permitted values of ${\theta }_{1}$ and ${\theta }_{2}$ are: (i) EP1: $[x,y]=[0,0]$; and (ii) EP2: $[x,y]=\left[\frac{1}{{\theta }_{1}{\theta }_{2}^{2}},\frac{1}{{\theta }_{1}{\theta }_{2}}\right]$.

The Jacobian evaluated at EP1 is $\left(\begin{array}{cc}-1 & 0\\ {\theta }_{2} & -1\end{array}\right)$. Thus the region of the (a, b) Cartesian coordinate plane representing the possible Jacobians associated with EP1 is simply the line a = 0. Similarly, the Jacobian evaluated at EP2 is $\left(\begin{array}{cc}-1 & 2/{\theta }_{2}\\ {\theta }_{2} & -1\end{array}\right)$, and hence the region representing the possible Jacobians associated with EP2 is the line $b=2/a.$

Example 3. We consider a further nonlinear example:

$$\begin{eqnarray*}{g}_{1}(y,\theta ) & = & {\theta }_{1}{y}^{3},\\ {g}_{2}(x,\theta ) & = & {\theta }_{2}x,\end{eqnarray*}$$

with ${\theta }_{1}$ and ${\theta }_{2}$ both non-zero. In this case, the Jacobian for the system is

$$\begin{eqnarray*}\left(\begin{array}{cc}-1 & 3{\theta }_{1}{y}^{2}\\ {\theta }_{2} & -1\end{array}\right),\end{eqnarray*}$$

which is again a function of ${\theta }_{1},{\theta }_{2}$, and y.

In this case, it is again straightforward to show that the only EPs that exist for all permitted values of ${\theta }_{1}$ and ${\theta }_{2}$ are: (i) EP1: $[x,y]=[0,0]$; and (ii) EPs 2 and 3: $[x,y]=\pm \left[\frac{{\theta }_{1}}{\sqrt{{\theta }_{1}^{3}{\theta }_{2}^{3}}},\frac{1}{\sqrt{{\theta }_{1}{\theta }_{2}}}\right]$.

The Jacobian evaluated at EP1 is again $\left(\begin{array}{cc}-1 & 0\\ {\theta }_{2} & -1\end{array}\right)$. Thus the region of the (a, b) Cartesian coordinate plane representing the possible Jacobians associated with EP1 is again the line a = 0. The Jacobian evaluated at EP 2 or 3 is $\left(\begin{array}{cc}-1 & 3/{\theta }_{2}\\ {\theta }_{2} & -1\end{array}\right)$, and hence the region representing the possible Jacobians associated with these EPs is the line $b=3/a.$

Assessing stability for these examples

For any system of n ODEs, the Jacobian matrix of the system evaluated at a particular EP ${{\bf{x}}}_{0}$ will be an element, J, of the set of $n\times n$ real matrices ${M}_{n}({\mathbb{R}})$. The EP ${{\bf{x}}}_{0}$ is locally stable if all of the eigenvalues of J have negative real part. An equivalent criterion is that the real part of the leading eigenvalue (i.e. the one having maximal real part) is negative.

We first consider the set of all $n\times n$ real matrices, ${M}_{n}({\mathbb{R}})$. The eigenvalues of any matrix $J\in {M}_{n}({\mathbb{R}})$ are the solutions of the characteristic equation,

$$\begin{eqnarray*}&&| J-\lambda {I}_{n}| =0,\end{eqnarray*}$$

where I_n denotes the n × n identity matrix [11]. We define ${\Lambda }_{S}^{n}\subset {M}_{n}({\mathbb{R}})$ to be the set of $n\times n$ matrices having all negative eigenvalues, and refer to this as the stable region of ${M}_{n}({\mathbb{R}})$.

The choice of a particular RME specifies a probability density function, h, on ${M}_{n}({\mathbb{R}})$. The stability probability associated with h is then

$$\begin{eqnarray}&&P({\rm{stable}}| h)={\displaystyle \int }_{J\in {\Lambda }_{S}^{n}}h(J){\rm{d}}J,\end{eqnarray} \tag{ 3 }$$

i.e. it is the total probability mass that falls within the stable region.

Crucially, the stability probability is determined by two factors: ${\Lambda }_{S}^{n}$ and h. ${\Lambda }_{S}^{n}$ is not random: for a given n it is a well-defined region of ${M}_{n}({\mathbb{R}})$. For systems of practical interest, however, this area cannot be determined analytically, and needs instead to be evaluated computationally, using e.g. Monte Carlo techniques (which are outlined below in 2.3). The results of stability analyses will therefore be completely determined by the choice of h, and how it distributes probability mass over the stable and unstable regions (see figure 2). If h is defined through the specification of an ODE model and a distribution for its parameters, then only matrices that can occur as Jacobians for that particular system will have non-zero density. Similarly, if h is defined without reference to a real system, then only some of the matrices in ${M}_{n}({\mathbb{R}})$ will be associated with non-zero density. However, for any given real ODE system, these matrices might not be obtainable as Jacobians, and—conversely—not all matrices obtainable as Jacobians will necessarily be associated with non-zero density by such a defined h, therefore limiting its relevance.

Zoom In Zoom Out Reset image size

An appropriate choice of h is thus vital. In particular, choosing h for the sake of mathematical convenience can only provide limited insight, if doing so comes at the cost of sacrificing realism. The so-called 'diversity-stability debate' [6] arose because general conclusions about stability were drawn from RMEs that were either specific to a particular system [25–32] or chosen for mathematical convenience [19, 20]—e.g. to invoke the circular law [5, 15, 21]—yielding results unlikely to be meaningful for any interesting realistic system [37].

Here we show that the dichotomy resulting from the use of different RMEs can be overcome by constructing RMEs that are appropriate for, and conditioned on, the properties of Jacobian matrices of real systems.

We consider autonomous ODE systems of the form

$$\begin{eqnarray*}&&\dot{{\bf{x}}}(t)={\bf{f}}({\bf{x}}(t);{\boldsymbol{\theta }}),\end{eqnarray*}$$

where ${\bf{x}}(t)={[{x}_{1}(t),\ldots ,{x}_{p}(t)]}^{\top }$ is the vector of state variables at time t, ${\boldsymbol{\theta }}$ is the vector of parameters, and ${\bf{f}}({\bf{x}}(t);{\boldsymbol{\theta }})={[{f}_{1}({\bf{x}}(t);{\boldsymbol{\theta }}),\ldots ,{f}_{p}({\bf{x}}(t);{\boldsymbol{\theta }})]}^{\top }.$ By definition, an EP, ${{\bf{x}}}_{{\boldsymbol{\theta }}}^{*}$, of the system has the property:

$$\begin{eqnarray*}&&{\bf{f}}({{\bf{x}}}_{{\boldsymbol{\theta }}}^{*};{\boldsymbol{\theta }})={\bf{0}}.\end{eqnarray*}$$

We include the subscript ${\boldsymbol{\theta }}$ in our notation for the EP to emphasise that its location, existence, and stability will generally depend upon the particular values taken by the parameters. We denote by ${J}_{{\boldsymbol{\theta }}}^{*}$ the Jacobian matrix of ${\boldsymbol{f}}$ evaluated at ${{\bf{x}}}_{{\boldsymbol{\theta }}}^{*}$.

We can induce an ensemble, ${\mathcal{R}}$, of Jacobian matrices by specifying a distribution, ${\mathcal{F}}$, for the parameters ${\boldsymbol{\theta }}$: a collection of N parameter vectors ${{\boldsymbol{\theta }}}^{(1)},\ldots ,{{\boldsymbol{\theta }}}^{(N)}$ sampled from ${\mathcal{F}}$ defines an ensemble of corresponding matrices ${J}_{{{\boldsymbol{\theta }}}^{(1)}}^{*},\ldots ,{J}_{{{\boldsymbol{\theta }}}^{(N)}}^{*}$. For any such RME, we may calculate a Monte Carlo estimate of the probability of stability, simply as the proportion of matrices that are stable; i.e. for which the leading eigenvalue has negative real part. That is, we obtain an estimate of the stability probability as,

$$\begin{eqnarray}&&\hat{P}({\rm{stable}}| {\mathcal{R}})\approx \displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\;{\mathbb{I}}({J}_{{{\boldsymbol{\theta }}}^{(i)}}^{*}\in {\Lambda }_{S}^{n}),\end{eqnarray} \tag{ 4 }$$

where ${\mathbb{I}}(X)$ is the indicator function, which is 1 if X is true and 0 otherwise.

The estimated stability probability defined in the previous section is the probability of a system being 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 for mesoscopic systems. To study the effects of this context, which is encapsulated by the statistical properties of, and dependencies among, the entries in the Jacobian matrices ${J}_{{{\boldsymbol{\theta }}}^{(1)}}^{*},\ldots ,{J}_{{{\boldsymbol{\theta }}}^{(N)}}^{*}$, we consider two further random matrix distributions, constructed by permutation of the entries of our original RME. First, we form a new matrix ensemble, ${K}_{(1)}^{*},\ldots ,{K}_{(N)}^{*}$, in which the dependency between entries is broken. For each ${\ell }\in \{1,\ldots ,N\}$ and $(i,j)\in \{1,\ldots ,p\}\times \{1,\ldots ,p\}$ we set ${\left({K}_{({\ell })}^{*}\right)}_{{ij}}={\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{ij}}$, with q drawn uniformly at random from $\{1,\ldots N\}$. In this way, the marginal distribution of the ij-entries across the ensemble of K* matrices is the same as the marginal distribution of ij-entries across the ensemble of ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices. Maintaining the marginal distributions ensures that the dependency between entries is the only quantity that we are altering: in particular, the location of zeros in the matrix and the magnitudes of interaction strengths are maintained. We construct a further RME, ${L}_{(1)}^{*},\ldots ,{L}_{(N)}^{*}$, where for each ${\ell }\in \{1,\ldots ,N\}$ and $(i,j)\in \{1,\ldots ,p\}\times \{1,\ldots ,p\}$, we set ${\left({L}_{({\ell })}^{*}\right)}_{{ij}}={\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{rs}}$, with q drawn uniformly at random from $\{1,\ldots N\}$, and r and s (independently) drawn uniformly at random from $\{1,\ldots p\}$. Now, the location of zeros in the matrix is no longer fixed; although the probability of an entry being zero is the same for the L* matrices as for the ${J}_{{\boldsymbol{\theta }}}^{*}$'s and K*'s. Moreover, each entry of the L* matrices is i.i.d. We henceforth refer to the ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices as the fully conditioned system (FCS) ensemble (most structure); the K* 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 the methods for their construction, in figure 3.

Zoom In Zoom Out Reset image size

To further our investigation we defined four more RMEs (which are presented in more detail in the supplementary information). The first one will be referred to as the independent normal ensemble. It is constructed as follows: For each (i, j), we fit an independent normal distribution to the ij-entries of the sampled Jacobians, ${J}_{{{\boldsymbol{\theta }}}^{(1)}}^{*},\ldots ,{J}_{{{\boldsymbol{\theta }}}^{(N)}}^{*}$. That is, for each (i, j), we calculate the mean,

$$\begin{eqnarray*}&&{\mu }_{(i,j)}^{\mathrm{ind}}=\displaystyle \frac{1}{N}\displaystyle \sum _{q=1}^{N}{\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{ij}},\end{eqnarray*}$$

and standard deviation,

$$\begin{eqnarray*}&&{\sigma }_{(i,j)}^{\mathrm{ind}}={\rm{s}}.{\rm{d}}.{\left\{{\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{ij}}\right\}}_{q=1}^{N}.\end{eqnarray*}$$

We then construct the new RME, ${M}_{(1)}^{\mathrm{ind}},\ldots ,{M}_{(N)}^{\mathrm{ind}}$, where for each ${\ell }\in \{1,\ldots ,N\}$ and $(i,j)\in \{1,\ldots ,p\}\times \{1,\ldots ,p\}$, we set ${\left({M}_{({\ell })}^{\mathrm{ind}}\right)}_{{ij}}$ to be a sample drawn from the univariate normal distribution with mean ${\mu }_{(i,j)}^{\mathrm{ind}}$ and standard deviation ${\sigma }_{(i,j)}^{\mathrm{ind}}$. By construction, the mean and standard deviation of the ij-entries across the ensemble of M^ind matrices are the same as the mean and standard deviation of the ij-entries across the ensemble of ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices (the FCS ensemble) and across the ensemble of K* matrices (the independent ensemble).

A further ensemble is given by the independent Pearson ensemble. As in the independent normal case defined above, this new RME 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 the mean and standard deviation, we also capture the skewness and kurtosis of the ij-entries of the ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices. That is, in addition to ${\mu }_{(i,j)}^{\mathrm{ind}}$ and ${\sigma }_{(i,j)}^{\mathrm{ind}}$ defined earlier, we also calculate skewness

$$\begin{eqnarray*}&&{\gamma }_{(i,j)}^{\mathrm{ind}}=\mathrm{skewness}{\left\{{\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{ij}}\right\}}_{q=1}^{N}.\end{eqnarray*}$$

and kurtosis,

$$\begin{eqnarray*}&&{\kappa }_{(i,j)}^{\mathrm{ind}}=\mathrm{kurtosis}{\left\{{\left({J}_{{{\boldsymbol{\theta }}}^{(q)}}^{*}\right)}_{{ij}}\right\}}_{q=1}^{N}.\end{eqnarray*}$$

We then construct an RME, ${M}_{(1)}^{\mathrm{pear}},\ldots ,{M}_{(N)}^{\mathrm{pear}}$, where for each ${\ell }\in \{1,\ldots ,N\}$ and $(i,j)\in \{1,\ldots ,p\}\;\times \{1,\ldots ,p\}$, we set ${\left({M}_{({\ell })}\right)}_{{ij}}^{\mathrm{pear}}$ to be a sample drawn from a univariate Pearson distribution with mean ${\mu }_{(i,j)}^{\mathrm{ind}}$, standard deviation ${\sigma }_{(i,j)}^{\mathrm{ind}}$, skewness ${\gamma }_{(i,j)}$, and kurtosis ${\kappa }_{(i,j)}$. This RME thus shares many of the properties of the marginal distributions of ij-entries across the ensemble of ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices, but does not capture the dependencies between them.

The third additional RME will be referred to as the i.i.d. normal ensemble. This time we will not fit a distribution to the ij-entries of the ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices, but instead we fit a normal distribution, using the same technique that we 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 an RME that attempts to capture some of the dependencies between the entries of the ${J}_{{\boldsymbol{\theta }}}^{*}$ matrices. We define c(M) to be the vector obtained by concatenating the columns of the matrix M (and further define ${c}^{-1}$ be the inverse operation, so that, for example, ${c}^{-1}(c(M))=M$). Applying $c(\cdot )$ to the matrices from our FCS RME, we obtain N vectors of length p × p, namely: $c({J}_{{{\boldsymbol{\theta }}}^{(1)}}^{*}),\ldots ,c({J}_{{{\boldsymbol{\theta }}}^{(N)}}^{*}).$ To these, we fit (by maximum likelihood) a $(p\times p)$-variate normal distribution. We then sample N vectors, ${v}_{1},\ldots ,{v}_{N}$, of length $p\times p$ from this distribution, and form a new ensemble ${M}_{(1)}^{\mathrm{mvn}},\ldots ,{M}_{(N)}^{\mathrm{mvn}}$ by setting ${M}_{(q)}^{\mathrm{mvn}}={c}^{-1}({v}_{q})$. We will call this new ensemble the multivariate 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 ensemble we can show the impact of the different moments of the distribution. The multivariate 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.

Stability, as stressed above and in the literature, is an issue in a wide variety of domains, and therefore we consider a set of systems that cover different qualitative behaviour of dynamical systems. The four ODE models that we consider have in common that the EPs and Jacobians can be identified analytically, which makes analysis straightforward; they are: (i) the Lorenz system [38]; (ii) a model of the cell cycle due to Tyson [39]; (iii) a model of viral dynamics due to Nowak and Bangham [40]; and (iv) an SEIR (susceptible-exposed-infective-recovered) population dynamics model [41]. For brevity, we will refer to these four example models 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 000 matrices from our RMEs in order to obtain Monte Carlo estimates of stability probabilities. Full details of these models 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 different models and RMEs. The eigenspectra distributions are typically much less dispersed for the FCS ensemble than for the other two. As shown in figure 4(C) this also leads to systematic differences in the real parts of the leading eigenvalues, which determine stability.

Zoom In Zoom Out Reset image size

Table 1 shows that, whatever the model 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 include more of the structure of the system in their construction, yield stability probabilities closer to those of the FCS. In most cases, the univariate Pearson ensemble had a probability stability closer to that of the FCS than the univariate normal, showing that considering more moments when building the RME improves the estimation of the stability of the system.

In summary, Monte Carlo estimates for the stability probabilities decrease as we decrease the amount of realism captured by the RME. 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 in which we tightly control the mean, standard deviation, skewness and kurtosis of the marginal distributions of Jacobian entries, demonstrates that all of these properties also have an impact on stability.

To illustrate the effects of inadequately capturing model structure and parameter dependencies on the stability probabilities of larger systems, we consider extensions to the SEIR model (figure 5(A)). We allow for multiple subpopulations of exposed (E) individuals (in the supplementary information we 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). As we 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 between matrix entries be important, but also their heterogeneity. Heterogeneity changes the location of the centre of the matrix p.d.f. h, and also how it stretches in different directions, which modifies the proportion of probability mass falling in the stable region, ${\Lambda }_{S}^{n}\subset {M}_{n}({\mathbb{R}})$.

Zoom In Zoom Out Reset image size

Figure 5(C) shows how summaries of the distributions of leading eigenvalues change as we 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 and median values shift away from stable negative values toward unstable positive values. This is in stark contrast to the FCS case where, regardless of the number of exposed populations, the distribution of leading eigenvalues only has support on the negative real line (and hence we always have stability probability 1). Moreover, as we increase the number of exposed populations, the median value of the leading eigenvalue tends to become more negative. The independent case is again intermediate, with the median value staying relatively constant.

Figures S4 and S5 in the supplementary material, where we consider the effect of the i.i.d. normal, the independent normal, the independent Pearson, and the multivariate normal ensembles on these models, bring more evidence to our previous observations. As we would expect, the more of the underlying system's structure that we capture using our chosen RME, the closer we get to the stability probability estimated using the FCS ensemble. The i.i.d. normal RME, which does not include any more 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 and means 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, the multivariate normal RME (which allows us to model–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 normal RME. We found that the stability probabilities obtained using the multivariate normal RME are not consistently more different from FCS stability probabilities than those obtained using the independent Pearson RME. Thus, accounting for dependency between the Jacobian's entries may, depending on the problem at hand, be more or less important than accounting for higher order properties of the marginal distributions of the entries.

For each of the models, the EPs were found using Matlab's analytic equation solver, solve. The Jacobian was also described analytically using Matlab's jacobian function.

The different parameters were sampled from a uniform distribution using Matlab's rand function. The choice of the range of the parameters will be described in detail below.

A.1.1. Algorithm

For each of the models, the EPs were found using Matlab's equation solver, solve. The Jacobian was described analytically using Matlab's jacobian function. The different parameters were sampled from a uniform distribution using Matlab's rand function. The choice of the range of the parameters will be described in detail below.

This method was used only on the S(nE)IR and the SE(nI)R models for computational reason (the numerical method becoming too expensive with big n). This method works well in this case because for these models we know that there are 2 EPs and they are easily identified, so we can easily segregate the cases corresponding to each of them when computing the probability of stability. The first EP, where the whole population is composed of recovered individuals, is not interesting because it is obviously stable, so we only consider the other EP in each system.

A.2.1. Algorithm

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 the main article, we sampled the parameters from the following ranges:

$$\begin{eqnarray*}\beta & \in [0,10],\\ \rho & \in [1,11],\\ \sigma & \in [0,10].\end{eqnarray*}$$

ρ is considered to be bigger or equal to 1 in order to ensure more than just the origin as a equilibrium point, which makes our study more interesting.

A.3.2. Model 2: a model of the cell division cycle

A.3.3. Model 3: the Nowak and Bangham model

A.3.4. Models 4, 5 and 6: the SEIR and extended SEIR models