## Abstract

We study an evolutionary model of a complex system that evolves under catalytic dynamics and Darwinian selection and exhibits spontaneous growth, stasis and then a collapse of its structure. We find that the typical lifetime of the system increases sharply with the diversity of its components or species. We also find that the prime reason for crashes is a naturally occurring internal fragility of the system. This fragility is captured in the network organizational character and is related to a reduced multiplicity of pathways or feedback loops between its components. These results apply to several generalizations of the model as well. This work suggests new parameters for understanding the robustness of evolving molecular networks, ecosystems, societies and markets.

## 1. Introduction

Several systems involving complex networks of interacting components exhibit dramatic crashes during the course of their evolution. Examples include mass extinctions in the biosphere as evidenced in the palaeontological record (Sepkowski 1984; Raup 1991; Erwin 2006), collapses of ecosystems (Paine 1969; May 1973; Pimm 1991; Solé & Bascompte 2006) and civilizations (Tainter 1990; Diamond 2006) and crashes of economic systems (Schumpeter 1939; Sornette 2004). While some of these catastrophic events are caused by large external perturbations such as meteorite impacts (Raup 1991; Hallam & Wignall 1997), earthquakes (Nur & Burgess 2008), famines, wars and infections (Diamond 2005), for the vast majority of them no single dramatic cause can be traced. An alternative reason for crashes, explored in this paper, is a fragility in the internal organization of these systems that naturally develops in the course of their evolution, making them vulnerable to small perturbations.

Empirical data characterizing the ‘internal fragility’ or ‘robustness’ of evolving complex systems are scarce. One of the chief problems in collecting data is a lack of clarity about what to look for; we do not know what system parameters can characterize its poisedness for a crash. Hence, a key step in identifying possible signatures of fragility is to construct theoretical and mathematical models of systems that exhibit repeated catastrophes in the course of their time evolution, whose analysis can reveal structural and dynamical features that make them vulnerable to such events. An important aspect of a complex system's organizational structure is the underlying interaction network of its components (Watts & Strogatz 1998; Bhalla & Iyengar 1999; Albert & Barabasi 2002; Dorogovtsev & Mendes 2003), hence we need in particular to study examples of evolving networks that exhibit crashes and recoveries.

Our model system (Jain & Krishna 1998, 2001, 2002*a*,*b*) exhibits crashes and recoveries for an evolving network of interacting populations, with Darwinian selection and dynamic feedback loops playing an important role in system evolution. In this paper, we discuss how the average system lifetime between crashes depends upon the diversity of system components, which has not been explored earlier. This is made possible by going to much larger system sizes (higher diversity) than considered earlier. We also explore for the first time in this model the dependence of the average number of feedback loops in the network on diversity, and find that this number is strongly correlated with average system lifetime. This shows that the number of feedback loops is a network structural signature that is relevant for crashes.

## 2. The model

The system consists of *s* nodes, whose network of interactions is specified completely by its adjacency matrix *C*≡(*c*_{ij}), *i*, *j*=1, …, *s*. A node may represent a molecular species in a ‘prebiotic pond’. The model is motivated by the origin of life problem (Dyson 1985; Bagley *et al*. 1991; Kauffman 1993), but may be more generally valid. The element *c*_{ij}=1 if species *j* ‘catalyses’ the growth of species *i* and zero otherwise. Also, *c*_{ii}=0 for all *i*, corresponding to the exclusion of self-catalysing species. Relaxing the above restrictions by allowing links with different weights and negative signs does not change the qualitative behaviour of the model.

Using the adjacency matrix, we write an equation for the population dynamics of the species given by
2.1
Here, is the rate of change of the population of species *i*. The first term on the right takes into account the positive effect of all the species that catalyse species *i*, each one having an effect proportional to its population. The second term is a constant mortality term.

The rate of growth of the relative populations, *x*_{i}, can be computed from equation (2.1). Applying that equation to the numerator and denominator of one gets an equation independent of *ϕ*:
2.2
For any given graph *C*, the dynamics described by equation (2.2) flows to a steady state (a fixed point of equation (2.2)) in which all *x*_{i} become time independent constants.

Initially, the matrix *C* is sparse and random with on average *m* links per node, with *m*<1 (i.e. each off-diagonal entry of *C* is independently chosen to be unity with a probability *p*=*m*/(*s*−1), and zero with probability 1−*p*). To introduce evolution into the model, we note that the pond can be washed by nearby tides, floods or storms that can flush out some of the contents of the pond. We introduce ‘selection’ by imposing that the species with the lowest relative population in the steady state gets removed from the system (Bak & Sneppen 1993); we eliminate the corresponding node and all its links from the graph. (If there is more than one such species, we choose one at random.) Furthermore, such a fluctuation can bring new species into the pond; we assume for simplicity that a single new node gets added to the graph whose links with the existing ones are made randomly with the same average connectivity *m*. After each such fluctuation, the populations evolve according to equation (2.2) with a fixed *C* to reach a new fixed point, whereafter the above update sequence is repeated.

The model implicitly assumes the existence of two time scales in the evolutionary process. One is a short time scale characterizing the population dynamics of the species in which the network remains constant. The second time scale, on which the graph is updated, characterizes the evolution of the network (arrival and death of species) and is taken to be longer than the time in which populations reach their attractor. The sharp sequential separation between these two kinds of processes is a model artefact for mathematical simplicity. Nevertheless, these two processes—internal system dynamics and evolution of the network—do often occur on different time scales in several complex systems. Note that, since one node is eliminated and introduced at each time step in the model, the average lifetime of a node is *s* update time steps. Thus, evolutionary time, in units of the typical lifetime of a species, is marked in the model by *n*/*s*, where *n* is the number of graph updates. In this paper, our main concern is another, much longer time scale that corresponds to the lifetime of the system as a whole. Unlike the two smaller time scales mentioned earlier, this new time scale is not put in by hand in the model but arises dynamically. It corresponds to the time scale over which system level fragilities build up spontaneously in the network of interacting species and cause large crashes.

As an example of a complex system that is far removed from the prebiotic chemical network, but to which the above model might still apply, we mention ecology of firms and other agents in an economic system. Here, a non-zero value of *c*_{ij} would represent a catalytic link from the economic output of firm *j* to that of firm *i*. For example, firm *j* may provide sustenance to firm *i* by means of finance, raw materials or inputs used in goods manufactured by *i*, trained manpower, services, etc. The ‘population’ variable *y*_{i} would in this context represent a variable describing the well being or economic performance of firm *i*, e.g. its profit, assets, turnover, etc. The first term in equation (2.1) captures the fact that a catalyst *j* can only contribute to *i*'s well being if it is performing well itself (in proportion to its own performance *y*_{j} to leading order). The second term would denote a drain on *i*'s resources in terms of interest payments on finance, raw material costs, wages, other expenses, etc. that would be in proportion to its own profits, assets, turnover, etc. The graph update rules would signify the removal of under-performing firms in a competitive scenario and the birth of new firms with new relationships to the existing firms. The present model, albeit highly stylized, thus captures the evolution of an economic network of interacting firms, in a manner different from other models of economic and financial systems (Palmer *et al*. 1994; Bouchaud & Cont 1998; Farmer 2002; Marsili *et al*. 2004; Challet *et al*. 2005).

## 3. Results

### 3.1 Lifetime of the system

At each graph update the system suffers a structural perturbation that modifies *C*. The perturbation is small in that only one species is updated, affecting the links of only ∼*m*∼*O*(1) number of species. Since the update of *C* depends on populations, the long-time dynamics of the populations is highly nonlinear in spite of the simplicity of equations (2.1) and (2.2). The typical dynamics is shown in figure 1 where the number, *s*_{1}, of populated species (whose steady state *x*_{i}>0) is plotted against time (*n*/*s*) for three values of *s*=100, 300, 500 and fixed *m*=0.25. Initially, the graph is sparse and the number *s*_{1} is small. After a certain time, *s*_{1} begins to grow and soon reaches its maximum value *s*. Thereafter, the system exhibits a stasis for a certain time in which *s*_{1} fluctuates between *s* and *s*−1. In this latter state that we call the ‘organized state’, all species except possibly the one being picked for replacement have *x*_{i}>0. Thereafter, the system experiences a collapse in which *s*_{1} drops to a fraction of *s*. This is followed by a recovery and a repetition of the same kind of dynamics.

In this paper, we focus on crashes whereby, in a single update step, the number of populated species *s*_{1} goes from *s* or *s*−1 to a fraction less than or equal to *h* of *s*. We present results for *h*=0.50 and 0.75. While the absolute number of crashes depends upon *h*, the qualitative results are not very sensitive to its value. As shown in figure 1, for fixed *m*, the frequency of crashes comes markedly down with increasing *s*. Similarly, if we increase *m* for fixed *s*, the number of crashes again decreases markedly.

For given values of *m* and *s*, there is a typical lifetime before the network collapses. We define this time *τ* as the number of update steps spent in the organized state in a given run (typically 10^{6} steps long), divided by the number of crashes observed during that run. Each run is parametrized by *s* and *m* and crashes are defined with respect to the parameter *h*. Hence *τ* depends upon *s*, *m* and *h*. The dependence of *τ* on *s* and *m* is shown in figure 2, for *h*=0.75 (again we present the rescaled lifetime *τ*/*s*). For fixed *m*, *τ*/*s* grows exponentially with *s*. This behaviour is consistent with the empirical relation
3.1
The coefficient *α*(*m*,*h*) is an increasing function of *m* (and a weak function of *h*) whose quantitative behaviour is discussed later.

These results show that the system is more stable against crashes as its diversity, *s*, increases for fixed connectivity, *m*, and also as its connectivity, *m*, increases for a fixed diversity. We emphasize that even for low connectivity the system can be stabilized against collapse by increasing its diversity. It turns out that in the organized state, the average connectivity of the species is close to ∼1+*m*; hence, for the values of *m* given above the average connectivity is only slightly above one. Even such sparsely connected systems are stabilized in this model by a sufficient amount of diversity.

### 3.2 The number of feedback loops in the network is correlated with system lifetime

We now attempt to understand this behaviour in terms of the structure of the graph near and far from a crash. The organized state has the structure of an *autocatalytic set* (ACS). An ACS is a subgraph, each of whose nodes has at least one incoming link from a node belonging to the same subgraph (Kauffman 1993). In the organized state, all the species except possibly the one being picked for replacement are part of the ACS (Jain & Krishna 1998). The ACS consists of a *core* and a *periphery*. The core comprises the set of nodes (along with their mutual links) from which there is a directed path to every other node in the ACS. All other nodes and links in the ACS constitute the periphery. Examples of the graph (ACS with core and periphery) observed in the organized state are shown in figure 3. By definition, there is no directed path from a periphery node to any core node. The core, by virtue of closed paths inside it, is a ‘self-sustaining’ structure in the sense that all the core nodes would be populated even if the only links present in the graph are those in the core. By contrast, the periphery nodes would become depopulated if the links from the core to the periphery were to be removed. In this sense, the periphery nodes are ‘parasites’ that are sustained by the core.

While there is always by definition at least one path from every core node to every other core node, the number of such paths is significantly different between a normal organized state and a state poised for a crash. In the typical organized state, there are several paths from each core node to another (figures 3*a*,*c*). In such configurations, no single node addition or deletion can cause a crash. However, the number of paths between core nodes drops to a much lower value just before a crash (figures 3*b*,*d*). Then, a single node change can disrupt the core and cause most species in the network to be depopulated.

We computed the number of distinct, non-intersecting closed paths of all lengths in the core of the graph at different time steps using a standard graph theoretic algorithm (see the review by Deo (1974) and references therein). In figure 4, we plot the frequency distribution of the number of these paths in the organized state (filled circles). The distribution shows a peak whose position, *N*_{p}, is dependent upon *s* and *m*. A plot of ln(*N*_{p}) against *s* for various values of *m* is shown in figure 5. This is consistent with the empirical formula
3.2
We note that loops in other graph ensembles have also been counted (Bianconi & Marsili 2005).

In figure 4, the open symbols show the distribution of closed paths in the core just before crashes. Its peak occurs at a much smaller value than *N*_{p} (note that the *x*-axis scale is logarithmic). This is also evident from figure 3 (the cores in (*b*) and (*d*) have much fewer closed paths than in (*a*) and (*c*)).

We find a strong correlation between the coefficients *α*(*m*) and *β*(*m*). This is shown in figure 6 where *α*(*m*,*h*) for two values of *h* and *β*(*m*) are plotted against *m*. It is seen that the dependence of *α* on *h* is weak, as mentioned before, and that *α* and *β* have a similar dependence upon *m*. Thus *N*_{p} and *τ*/*s* have a similar dependence on *m* and *s*. This close correspondence between a structural property such as the number of loops in the graph in the organized phase and a dynamical property such as the lifetime of that phase is one of the surprising results we have found.

This suggests an explanation of why a higher diversity and density of links enhances stability against crashes in this model. Diversity increases the number of closed paths in the core and thus provides a buffer against crashes by ensuring alternative routes of sustenance in the event of loss of core nodes. Crashes occur typically when the core has thinned out, and such fragile states take longer to be realized when there is a larger number of paths in the core to begin with.

We remark that unlike the models considered by May (1972, 1973), in the present model crashes are not due to the instability of the population dynamics with respect to fluctuations of *x*_{i} about its fixed point. For any graph, the attractor of equation (2.2), denoted ** X**, is an eigenvector of

*C*corresponding to its largest eigenvalue. For a generic non-negative matrix

*C*that arises in the organized state, the largest eigenvalue is non-degenerate and the corresponding eigenvector

**is therefore a unique, global attractor (independent of initial conditions), stable against perturbations of the**

*X**x*

_{i}(all eigenvalues of the appropriate Jacobian arising from equation (2.2), evaluated at

**=**

*x***, are negative). Crashes occur because the graph changes structurally in a graph update, and sometimes (especially where the core has thinned out) a small change in the graph at just one node can cause the population attractor of the new graph to be quite different from the old one, with many of the species populated in the former attractor unpopulated in the new one. To make this explicit, let**

*X**n*be a graph update time step at which a crash occurs in the organized state, i.e.

*s*

_{1}(

*n*−1) equals

*s*or

*s*−1, and

*s*

_{1}(

*n*)≤

*hs*. The graphs at the two consecutive time steps are denoted

*C*(

*n*−1) and

*C*(

*n*), and the corresponding attractors of equation (2.2) (steady state relative population vectors) as

**(**

*X**n*−1) and

**(**

*X**n*). Before the graph update at

*n*, the populations are in the configuration

**(**

*X**n*−1) that is a stable fixed point of equation (2.2) with

*C*=

*C*(

*n*−1) (the appropriate Jacobian has all eigenvalues negative, independent of parameter values

*s*and

*m*). After the graph update

**(**

*X**n*−1) is no longer the attractor of equation (2.2) since

*C*has changed from

*C*(

*n*−1) to

*C*(

*n*). The populations move to the configuration

**(**

*X**n*) that is a stable fixed point of equation (2.2) with

*C*=

*C*(

*n*). The crash occurs because a small structural change in

*C*(note that the matrices

*C*(

*n*−1) and

*C*(

*n*) differ in just

*O*(1) number of entries in the row and column corresponding to the replaced node) has caused a large change in the attractor of equation (2.2) (the number of non-zero components of

**(**

*X**n*−1) is

*s*or

*s*−1 and of

**(**

*X**n*) is ≤

*hs*). Our results imply that as diversity and hence the number of feedback loops increases, structural changes that arise in the course of evolution are less likely to cause such serious changes in the population attractors.

## 4. Discussion

As in real evolutionary systems, the model generates several dynamical time scales. The model has only two parameters: system size or diversity, *s*, and the average connectivity of a new node, *m* (the latter being typically *O*(1)). In spite of its extreme simplicity, the time scales that dynamically appear have a wide range of dependence on *s*, including logarithmic, power law and exponential. The time scale for the appearance of an ACS is independent of *s* and the time scale of its growth is ∼ln *s* (Jain & Krishna 1998, 2001; at constant *m*, in scaled units of time as used in figure 1). Once a crash sets in it occurs fast—on a time scale ∼1/*s* in the present version of the model. The fast collapse with a relatively slower recovery seen in the model is an observed feature in the fossil record as well as stock markets. The lifetime of the system between its growth and collapse has turned out to be the time scale that is the most sensitive to its diversity, namely, ∼e^{αs}, as shown here. Such a dependence means that there is a threshold scale of diversity set by 1/*α*, such that if diversity is well above this scale, the system is robust to crashes, but if it is close to or lower it is vulnerable.

The dynamics of growth and collapse in our model is different from other existing models, including various models of extinction studied in the literature (see the review by Newman & Palmer (2003) and references therein). The seed for the growth of complexity in this model is a small feedback loop (usually a two-cycle) that arises in the network by chance. The cooperativity implicit in this autocatalytic structure causes its nodes to have much higher populations than other nodes. Under a selection dynamics that preferentially preserves nodes with higher population, such a structure is stable and grows in complexity until it spans the whole system. Then, the same selection dynamics causes its components, erstwhile cooperators, to become competitors. This happens because the eliminated node now belongs to the ACS, whereas before spanning it was outside the ACS. This leads to internal organizational restructuring of the ACS, and, on a certain time scale, when its internal feedback loops become sparse, to fragility. Thus, we have here an example of how the very success and domination of a certain organizational structure changes the effective rules of the game leading to the collapse of the structure (for another such example see the model of Cohen *et al*. (2001)). This is reminiscent of how certain civilizations and organizations collapse (Tainter 1990; Diamond 2006). The role of feedback loops in a network structure that evolves under both selection and stochastic forces is also characteristic of several real evolutionary systems.

It is significant that in several generalizations of this model introduced earlier (Jain & Krishna 2002*b*; Krishna 2004), we find that the lifetime of the evolutionary network increases exponentially with the diversity. This includes models where the network allows self loops, both positive and negative links (that inhibit species production), and links with varying strengths. This also includes a model in which the number of species *s* is itself a dynamical variable. Here, instead of eliminating the least populated node, all species whose relative population falls below a specified threshold are eliminated at every graph update time step, and only one new node is brought into the network. Our main result that system lifetime grows exponentially with *s* and that fragility to crashes arises when the density of feedback pathways thins out holds for all these models. Details will be presented elsewhere (Mehrotra *et al*. in preparation).

It would be interesting to consider more general classes of network models showing crashes and recoveries (including those that capture the dynamics of organisms, ecosystems, societies, etc. more realistically) and to explore the extent to which they share the behaviour of our simple idealized model. The linear terms in equation (2.1), motivated by catalytic chemical production, would need to be replaced by more complex functional responses used in population dynamics. The construction and analysis of such models showing non-trivial crashes and recoveries of network structure remains as an important open problem.

Mathematical models of ecosystems suggest that several factors determine ecosystem stability under various types of perturbations (see the reviews of McCann (2000), McKane & Drossel (2006), Montoya *et al*. (2006), Pascual *et al*. (2006) and references therein). The suggested importance of the multiplicity of sustenance pathways of species (MacArthur 1955) is analogous to the result we have found above. Note that as in the core of our graphs, so in ecosystems at the most basic level there exists several feedback loops between plants and microbial communities that feed on detritus and restore soil nutrients and other components of the atmosphere that are recycled in nature. Disruption of these feedback pathways would, beyond a certain point, be catastrophic for the ecosystem as a whole. Most ecosystem models concerned with stability typically take into account only the plants and higher trophic levels, and exclude microbes (see, however, the work of Neutel *et al*. (2007)), detritus and inorganic molecules that provide essential feedback loops. Our work suggests that newer and perhaps clearer patterns may emerge when models and field data are considered that include these components along with other trophic levels. It may also be of interest to explore connections with networks of mutualistic interactions among species that contain nested feedback loops and are regarded as being important for biodiversity and evolution of ecosystems (Bascompte *et al*. 2003; Thompson 2005).

To summarize, we have considered the effect of structural perturbations of species (node/link) removal and introduction that arise in the natural course of evolution in a highly self-organized network. This work shows how an increase in diversity and link density can contribute to long-term system stability against crashes caused by such perturbations by increasing the cooperative routes of sustenance in the network.

## Acknowledgments

We thank Sandeep Krishna for collaboration during the initial phase of this work and Areejit Samal for help with graph visualization. S.J. acknowledges support from the Robustness programme of the Santa Fe Institute and a grant from the Department of Biotechnology, Government of India.

## Footnotes

- Received September 20, 2008.
- Accepted October 31, 2008.

- © 2008 The Royal Society