Component response rate variation underlies the stability of highly complex finite systems.

The stability of a complex system generally decreases with increasing system size and interconnectivity, a counterintuitive result of widespread importance across the physical, life, and social sciences. Despite recent interest in the relationship between system properties and stability, the effect of variation in response rate across system components remains unconsidered. Here I vary the component response rates (γ) of randomly generated complex systems. I use numerical simulations to show that when component response rates vary, the potential for system stability increases. These results are robust to common network structures, including small-world and scale-free networks, and cascade food webs. Variation in γ is especially important for stability in highly complex systems, in which the probability of stability would otherwise be negligible. At such extremes of simulated system complexity, the largest stable complex systems would be unstable if not for variation in γ. My results therefore reveal a previously unconsidered aspect of system stability that is likely to be pervasive across all realistic complex systems.

Introduction

In 1972, May[1] first demonstrated that randomly assembled systems of sufficient complexity are almost inevitably unstable given infinitesimally small perturbations. Complexity in this case is defined by the size of the system (i.e., the number of potentially interacting components; ), its connectance (i.e., the probability that one component will interact with another; ), and the variance of interaction strengths ()[2]. May’s finding that the probability of local stability falls to near zero given a sufficiently high threshold of is broadly relevant for understanding the dynamics and persistence of systems such as ecological[1-6], neurological[7,8], biochemical[9,10], and socio-economic[11-14] networks. As such, identifying general principles that affect stability in complex systems is of wide-ranging importance.

Randomly assembled complex systems can be represented as large square matrices () with components (e.g., networks of species[2] or banks[12]). One element of such a matrix, , defines how component affects component in the system at a point of equilibrium[2]. Off-diagonal elements () therefore define interactions between components, while diagonal elements () define component self-regulation (e.g., carrying capacity in ecological communities). Traditionally, off-diagonal elements are assigned non-zero values with a probability , which are sampled from a distribution with variance ; diagonal elements are set to −1[1,2,5]. Local system stability is assessed using eigenanalysis on , with the system being stable if the real parts of all eigenvalues (), and therefore the leading eigenvalue (), are negative ()[1,2]. In a large system (high ), eigenvalues are distributed uniformly[15] within a circle centred at (−d is the mean value of diagonal elements) and , with a radius of [1,2,5] (Fig. 1a). Local stability of randomly assembled systems therefore becomes increasingly unlikely as , , and increase.

Figure 1

Eigenvalue distributions of random complex systems. Each panel shows the real (x-axis) and imaginary (y-axis) parts of eigenvalues from random matrices. (a) A system represented by a matrix A, in which all elements are sampled from a normal distribution with and . Points are uniformly distributed within the blue circle centred at the origin with a radius of . (b) The same system as in a after including variation in the response rates of components, represented by the diagonal matrix , such that . Elements of are randomly sampled from a uniform distribution from to . Eigenvalues of are then distributed non-uniformly within the red circle centred at the origin with a radius of . (c) A different random system constructed from the same parameters as in a, except with diagonal element values of −1. (d) The same system c after including variation in component response rates, sampled from as in b.

May’s[1,2] stability criterion assumes that the expected response rates () of individual components to perturbations of the system are identical, but this is highly unlikely in any complex system. In ecological communities, for example, the rate at which population density changes following perturbation will depend on the generation time of organisms, which might vary by orders of magnitude among species. Species with short generation times will respond quickly (high ) to perturbations relative to species with long generation times (low ). Similarly, the speed at which individual banks respond to perturbations in financial networks, or individuals or institutions respond to perturbations in complex social networks, is likely to vary. The effect of such variance on stability has not been investigated in complex systems theory. Intuitively, variation in () might be expected to decrease system stability by introducing a new source of variation into the system and thereby increasing . Here I show that, despite higher , realistic complex systems (in which is high but finite) are actually more likely to be stable if their individual component response rates vary. My results are robust across commonly observed network structures, including random[1], small-world[16], scale-free[17], and cascade food web[18,19] networks.

Results

Component response rates of random complex systems

Complex systems () are built from two matrices, one modelling component interactions (), and second modelling component response rates (). Both and are square matrices. Rows in define how a given component is affected by each component in the system, including itself (where ). Off-diagonal elements of are independent and identically distributed (i.i.d.), and diagonal elements are set to as in May[1]. Diagonal elements of are positive, and off-diagonal elements are set to zero (i.e., is a diagonal matrix with positive support). The distribution of diag() over components thereby models the distribution of component response rates. The dynamics of the entire system can be defined as follows[20],

Equation 1 thereby serves as a null model to investigate how variation in component response rate () affects complex systems. In the absence of such variation (), is set to the identity matrix (diagonal elements all equal 1) and . Under these conditions, eigenvalues of are distributed uniformly[15] in a circle centred at with a radius of [1] (Fig. 1a).

Effect of on M (co)variation

The value of , and therefore system stability, can be estimated from five properties of [21]. These properties include (1) system size (), (2) mean self-regulation of components (), (3) mean interaction strength between components (), (4) the variance of between component interaction strengths (hereafter , to distinguish from and ), and (5) the correlation of interaction strengths between components, and ()[22]. Positive does not change , nor does it necessarily change or . What does change is the total variation in component interaction strengths (), and . Introducing variation in increases the total variation in the system. Variation in the off-diagonal elements of is described by the joint variation of two random variables,

Given and , Eq. 2 can be simplified,

The increase in caused by can be visualised from the eigenvalue spectra of A versus (Fig. 1). Given and , the distribution of eigenvalues of A and M lie within a circle of a radius and , respectively (Fig. 1a vs. 1b). If , positive changes the distribution of eigenvalues[23-25], potentially affecting stability (Fig. 1c vs. 1d).

Given , increases linearly with such that[26],

If , such as when models a predator-prey system in which and have opposing signs, stability increases[2]. If diagonal elements of vary independently, the magnitude of is decreased because increases the variance of without affecting the expected covariance between and (Fig. 2).

Figure 2

Complex system correlation versus stability with and without variation in component response rates. Each point represents 10000 replicate numerical simulations of a random complex system with a fixed correlation between off-diagonal elements and (, x-axis). Where real parts of eigenvalues of are negative (y-axis), is stable (black dotted line). Blue circles show systems in the absence of variation in component response rates (). Red squares show systems in which . Arrows show the range of real parts of leading eigenvalues observed. Because decreases the magnitude of , purple lines are included to link replicate simulations before (blue circles) and after (red squares) including . The range of values in which decreases the mean real part of the leading eigenvalue is indicated with grey shading. In all simulations, system size and connectance were set to and , respectively. Off-diagonal elements of A were randomly sampled from , and diagonal elements were set to −1. Elements of were sampled, .

Numerical simulations of random systems with and without

I used numerical simulations and eigenanalysis to test how variation in affects stability in random matrices with known properties, comparing the stability of A versus . Values of were sampled from a uniform distribution where and (see Supplementary Information for other distributions, which gave similar results). In all simulations, diagonal elements were standardised to ensure that −d between individual and pairs were identical (also note that ). First I focus on the effect of across values of , then for increasing system sizes () in random and structured networks. By increasing , the objective is to determine the effect of as system complexity increases toward the boundary at which stability is realistic for a finite system.

Simulation of random M across ρ

Numerical simulations revealed that results in a nonlinear relationship between and , which can sometimes increase the stability of the system. Figure 2 shows a comparison of across values for () versus M () given , , and . For (shaded region of Fig. 2), expected was lower in than . For , the lower bound of the range of values also decreased given , resulting in negative in for and . Hence, across a wide range of system correlations, variation in the response rate of system components had a stabilising effect.

The stabilising effect of across increased with increasing . Figure 3 shows numerical simulations of across increasing given and ( has been lowered here to better illustrate the effect of ; note that now given , ). For relatively small systems (), never decreased the expected . But as increased, the curvilinear relationship between and decreased expected for given low magnitudes of . In turn, as increased, and systems became more complex, increased the proportion of numerical simulations that were observed to be stable (see below).

Figure 3

System correlation versus stability across different system sizes. In each panel, 10000 random complex systems are simulated for each correlation between off-diagonal elements and . Lines show the expected real part of the leading eigenvalues of M (red squares; ) versus A (blue circles; ) across , where negative values (below the dotted black line) indicate system stability. Differences between lines thereby show the effect of component response rate variation () on system stability across system correlations and sizes (). For all simulations, system connectance was . Off-diagonal elements of A were randomly sampled from , and diagonal elements were set to −1. Elements of were sampled such that , so .

Simulation of random M across S

To investigate the effect of on stability across systems of increasing complexity, I simulated random matrices at and across . One million were simulated for each , and the stability of vesus was assessed given (). For all , I found that the number of stable random systems was higher in than (Fig. 4; see Supplementary Information for full table of results), and that the difference between the probabilities of observing a stable system increased with an increase in . In other words, the potential for to affect stability increased with increasing system complexity and was most relevant for systems on the cusp of being too complex to be realistically stable. For the highest values of , nearly all systems that were stable given varying would not have been stable given .

Figure 4

Stability of large complex systems with and without variation in component response rate (γ). The log number of systems that are stable across different system sizes () given , and the proportion of systems for which variation in is critical for system stability. For each , 1 million complex systems are randomly generated. Stability of each complex system is tested given variation in by randomly sampling . Stability given is then compared to stability in an otherwise identical system in which for all components. Blue and red bars show the number of stable systems in the absence and presence of , respectively. The black line shows the proportion of systems that are stable when , but would be unstable if (i.e., the conditional probability that A is unstable given that M is stable).

I also simulated 100000 for three types of random networks that are typically interpreted as modelling three types of interspecific ecological interactions[2,27]. These interaction types are competitive, mutualist, and predator-prey, as modelled by off-diagonal elements that are constrained to be negative, positive, or paired such that if then , respectively[2] (but are otherwise identical to the purely random ). As increased, a higher number of stable relative to was observed for competitor and predator-prey, but not mutualist, systems. A higher number of stable systems was observed whenever and for competitive and predator-prey systems, respectively (note that for predator-prey systems, making stability more likely overall). The stability of mutualist systems was never affected by .

The effect of on stability did not change qualitatively across values of , , or for different distributions of (see Supporting Information).

Simulation of structured M across S

To investigate how affects the stability of commonly observed network structures, I simulated one million for small-world[16], scale-free[17], and cascade food web[18,19] networks. In all of these networks, rules determining the presence or absence of an interaction between components and constrain the overall structure of the network. In small-world networks, interactions between components are constrained so that the expected degree of separation between any two components increases in proportion to [16]. In scale-free networks, the distribution of the number of components with which a focal component interacts follows a power law; a few components have many interactions while most components have few interactions[17]. In cascade food webs, species are ranked and interactions are constrained such that a species can only feed on if the rank of .

Network structure did not strongly modulate the effect that had on stability. For comparable magnitudes of complexity, structured networks still had a higher number of stable than . For random networks, increased stability given ( and ), and therefore complexity . This threshold of complexity, above which more than were stable, was comparable for small-world networks, and slightly lower for scale-free networks (note that algorithms for generating small-world and scale-free networks necessarily led to varying ; see methods). Varying increased stability in cascade food webs for , and therefore at a relatively low complexity magnitudes compared to random predator-prey networks (). Overall, network structure did not greatly change the effect that had on increasing the upper bound of complexity within which stability might reasonably be observed.

System feasibility given

For complex systems in which individual system components represent the density of some tangible quantity, it is relevant to consider the feasibility of the system. Feasibility assumes that values of all components are positive at equilibrium[6,28,29]. This is of particular interest for ecological communities because population density () cannot take negative values, meaning that ecological systems need to be feasible for stability to be biologically realistic[28]. While my results are intended to be general to all complex systems, and not restricted to species networks, I have also performed a feasibility analysis on all matrices tested for stability. I emphasise that is not interpreted as population density in this analysis, but instead as a fundamental property of species life history such as expected generation time. Feasibility was unaffected by and instead occurred with a fixed probability of , consistent with a recent proof by Serván et al.[30] (see Supplementary Information). Hence, for pure interacting species networks, variation in component response rate (i.e., species generation time) does not affect stability at biologically realistic species densities.

Targeted manipulation of γ

To further investigate the potential of to be stabilising, I used a genetic algorithm. Genetic algorithms are heuristic tools that mimic evolution by natural selection, and are useful when the space of potential solutions (in this case, possible combinations of values leading to stability in a complex system) is too large to search exhaustively[31]. Generations of selection on value combinations to minimise demonstrated the potential for to increase system stability. Across , sets of values were found that resulted in stable systems with probabilities that were up to four orders of magnitude higher than when (see Supplementary Information), meaning that stability could often be achieved by manipulating values rather than elements (i.e., by manipulating component response rates rather than interactions between components).

Discussion

I have shown that the stability of complex systems might often be contigent upon variation in the response rates of their individual components, meaning that factors such as rate of trait evolution (in biological networks), transaction speed (in economic networks), or communication speed (in social networks) need to be considered when investigating the stability of complex systems. Variation in component response rate is more likely to be critical for stability in systems that are especially complex, and it can ultimately increase the probability that system stability is observed above that predicted by May’s[1] classically derived criterion. The logic outlined here is general, and potentially applies to any complex system in which individual system components can vary in their reaction rates to system perturbation.

It is important to recognise that variation in component response rate is not stabilising per se; that is, adding variation in component response rates to a particular system does not increase the probability that the system will be stable. Rather, highly complex systems that are observed to be stable are more likely to have varying component response rates, and for this variation to be critical to their stability (Fig. 4). This is caused by the shift to a non-uniform distribution of eigenvalues that occurs by introducing variation in (Fig. 1), which can sometimes cause all of the real components of the eigenvalues of the system matrix to become negative, but might also increase the real components of eigenvalues.

My focus here is distinct from Gibbs et al.[24], who applied the same mathematical framework to investigate how a diagonal matrix (equivalent to in my model) affects the stability of a community matrix given an interaction matrix within a generalised Lotka-Volterra model, where . Gibbs et al.[24] analytically demonstrated that the effect of on system stability decreases exponentially as system size becomes arbitrarily large () for a given magnitude of complexity . My numerical results do not contradict this prediction because I did not scale , but instead fixed and increased to thereby increase total system complexity (see Supplemental Information for results simulated across and ). Overall, I show that component response rate variation increases the upper bound of complexity at which stability can be realistically observed, meaning that highly complex systems are more likely than not to vary in their component response rates, and for this variation to be critical for system stability.

Interestingly, while complex systems were more likely to be stable given variation in component response rate, they were not more likely to be feasible, meaning that stability was not increased when component values were also restricted to being positive at equilibrium. Feasibility is important to consider, particularly for the study of ecological networks of species[6,25,28,30] because population densities cannot realistically be negative. My results therefore suggest that variation in the rate of population responses to perturbation (e.g., due to differences in generation time among species) is unlikely to be critical to the stability of purely multi-species interaction networks (see also Supplementary Information). Nevertheless, ecological interactions do not exist in isolation in empirical systems[20], but instead interact with evolutionary, abiotic, or social-economic systems. The relevance of component response rate for complex system stability should therefore not be ignored in the broader context of ecological communities.

The potential importance of component response rate variation was most evident from the results of simulations in which the genetic algorithm was used in attempt to maximise the probability of system stability. The probability that some combination of component response rates could be found to stabilise the system was shown to be up to four orders of magnitude higher than the background probabilities of stability in the absence of any component response rate variation. Instead of manipulating the interactions between system components, it might therefore be possible to manipulate only the response rates of individual system components to achieve stability. Hence, managing the response rates of system components in a targeted way could potentially facilitate the stabilisation of complex systems through a reduction in dimensionality.

A general mathematical framework encompassing shifts in eigenvalue distributions caused by a diagonal matrix γ has been investigated[23] and recently applied to questions concerning species density and feasibility[24,25], but γ has not been interpreted as rates of response of individual system components to perturbation. My model focuses on component response rates for systems of a finite size, in which complexity is high but not yet high enough to make the probability of stability unrealistically low for actual empirical systems. For this upper range of system size, randomly assembled complex systems are more likely to be stable if their component response rates vary (e.g., for parameter values in Fig. 4). Variation in component response rate might therefore be critical for maintaining stability in many highly complex empirical systems. These results are broadly applicable for understanding the stability of complex networks across the physical, life, and social sciences.

Methods

Component response rate (γ) variation

In a synthesis of eco-evolutionary feedbacks on community stability, Patel et al.[20] model a system that includes a vector of potentially changing species densities (n) and a vector of potentially evolving traits (x). For any species or trait , change in species density () or trait value () with time () is a function of the vectors n and x,

In the above, and are functions that define the effects of all species densities and trait values on the density of a species and the value of trait , respectively. Patel et al.[20] were interested in stability when the evolution of traits was relatively slow or fast in comparison with the change in species densities, and this is modulated in the above by the scalar . The value of thereby determines the timescale separation between ecology and evolution, with high modelling relatively fast evolution and low modelling relatively slow evolution[20].

I use the same principle that Patel et al.[20] use to modulate the relative rate of evolution to modulate rates of component responses for components. Following May[1,32], the value of a component at time () is affected by the value of () and ’s marginal effect on (), and by ’s response rate (),

In matrix notation[32],

In the above, is a diagonal matrix in which elements correspond to individual component response rates. Therefore, defines the change in values of system components and can be analysed using the techniques of May[1,23,32]. In these analyses, row means of are expected to be identical, but variation around this expectation will naturally arise due to random sampling of off-diagonal elements and finite . In simulations, the total variation in row means that is attributable to is small relative to that attributable to , especially at high . Variation in specifically isolates the effects of differing component response rates, hence causing differences in expected row means.

Construction of random and structured networks

I used the R programming language for all numerical simulations and analyses[33]. Purely random networks were generated by sampling off-diagonal elements from an i.i.d. with a probability (unsampled elements were set to zero). Diagonal elements were set to −1. Elements of were simulated i.i.d. from a distribution with positive support (typically ). Random matrices with correlated elements and were built using Cholesky decomposition. Competitor networks in which off-diagonal elements were constructed by first building a random , then flipping the sign of any elements in which . Similarly, mutualist networks were constructed by building a random , then flipping the sign of elements where . Predator-prey networks were constructed by first building a random , then flipping the sign of either or if .

Small-world networks were constructed using the method of Watts and Strogatz[16]. First, a regular network[16] was created such that components were arranged in a circle. Each component was initially set to interact with its closest neighbouring components on each side, where was an even natural number (e.g., for the regular network forms a ring in which each component interacts with its two adjacent neighbours; see Supplemental Material for examples). Each interaction between a focal component and its neighbour was then removed and replaced with with a probability of . In replacement, a new component was randomly selected to interact with the focal component; selection was done with equal probability among all but the focal component. The resulting small-world network was represented by a square binary matrix in which 1s represented interactions between components and 0s represented the absence of an interaction. A new random matrix was then generated with elements sampled i.i.d. from . To build the interaction matrix , I used element-wise multiplication , then set . I set and simulated small-world networks across all combinations of and .

Scale-free networks were constructed using the method of Albert and Barabási[17]. First, a saturated network (all components interact with each other) of size was created. New components were then added sequentially to the network; each newly added component was set to interact with randomly selected existing components. When the system size reached , the distribution of the number of total interactions that components had followed a power-law tail[17]. The resulting network was represented by an binary matrix , where 1s and 0s represent the presence and absence of an interaction, respectively. As with small-world networks, a random matrix was generated, and . Diagonal elements were set to −1. I simulated scale-free networks across all combinations of and .

Cascade food webs were constructed following Solow and Beet[18]. First, a random matrix was generated with off-diagonal elements sampled i.i.d. so that . Each component in the system was ranked from to . If component had a higher rank than component and , then was multiplied by −1. If had a lower rank than and , then was multiplied by −1. In practice, this resulted in a matrix with negative and positive values in the lower and upper triangles, respectively. Diagonal elements of were set to −1 and . I simulated cascade food webs for .

System feasibility

Dougoud et al.[28] identify the following feasibility criteria for ecological systems characterised by interacting species with varying densities in a generalised Lotka-Volterra model,

In the above, n* is the vector of species densities at equilibrium. Feasibility is satisfied if all elements in n* are positive. The matrix is the identity matrix, and the value is the strength of intraspecific competition (diagonal elements). Diagonal values are set to −1, so . The variable is a normalisation parameter that modulates the strength of interactions () for . Implicitly, here underlying strong interactions. Hence, , so in the above, a diagonal matrix of −1s () is added to , which has a diagonal of all zeros and an off-diagonal affecting species interactions (i.e., the expression relates to May’s[1] stability criterion[28] by , and hence for my purposes ). Given , the above criteria is therefore reduced to the below (see also Serván et al.[30]),

To check the feasibility criteria for , I therefore evaluated ( elements were sampled i.i.d. from ). Feasibility is satisfied if all of the elements of the resulting vector are positive.

Genetic algorithm

Ideally, to investigate the potential of for increasing the proportion of stable complex systems, the search space of all possible diag() vectors would be evaluated for each unique . This is technically impossible because can take any real value between 0–2, but even rounding to reasonable values would result in a search space too large to practically explore. Under these conditions, genetic algorithms are highly useful tools for finding practical solutions by mimicking the process of biological evolution[31]. In this case, the practical solution is finding vectors of diag() that decrease the most positive real eigenvalue of . The genetic algorithm used achieves this by initialising a large population of 1000 different potential diag() vectors and allowing this population to evolve through a process of mutation, crossover (swaping values between vectors), selection, and reproduction until either a diag() vector is found where all or some “giving up” critiera is met.

For each , the genetic algorithm was run for 100000 random (, ). The genetic algorithm was initialised with a population of 1000 different diag() vectors with elements sampled i.i.d. from . Eigenanalysis was performed on the resulting from each , and the 20 diag() vectors resulting in with the lowest each produced 50 clonal offspring with subsequent random mutation and crossover between the resulting new generation of 1000 diag() vectors. Mutation of each in a diag() vector occurred with a probability of 0.2, resulting in a mutation effect of size being added to generate the newly mutated (any values that mutated below zero were multiplied by −1, and any values that mutated above 2 were set to 2). Crossover occurred between two sets of 100 diag() vectors paired in each generation; vectors were randomly sampled with replacement among but not within sets. Vector pairs selected for crossover swapped all elements between and including two randomly selected with replacement (this allowed for reversal of vector element positions during crossover; e.g., ). The genetic algorithm terminated if a stable was found, 20 generations occurred, or if the mean fitness increase between generations was less than 0.01 (where fitness was defined as for M).

Supplementary Information.