## Abstract

Linearized catalytic reaction equations (modelling, for example, the dynamics of genetic regulatory networks), under the constraint that expression levels, i.e. molecular concentrations of nucleic material, are positive, exhibit non-trivial dynamical properties, which depend on the average connectivity of the reaction network. In these systems, an inflation of the *edge of chaos* and multi-stability have been demonstrated to exist. The positivity constraint introduces a nonlinearity, which makes chaotic dynamics possible. Despite the simplicity of such *minimally nonlinear systems*, their basic properties allow us to understand the fundamental dynamical properties of complex biological reaction networks. We analyse the Lyapunov spectrum, determine the probability of finding stationary oscillating solutions, demonstrate the effect of the nonlinearity on the effective in- and out-degree of the *active interaction network*, and study how the frequency distributions of oscillatory modes of such a system depend on the average connectivity.

## 1. Introduction

Many complex systems in general—and living systems and cells in particular—display remarkable stability, i.e. a capacity to sustain their spatial and temporal molecular organization. Yet, their stability is dynamic, i.e. these systems—to a certain degree—are capable of adapting to changes in their physical and chemical environment. This has led several authors (Packard 1988; Langton 1990; Kauffman 1993; Mitchell *et al.* 1993) to interpret such systems as existing at the *edge of chaos*. Mathematically, the *edge of chaos* refers to regions in parameter space where the system dynamics is characterized by a maximal Lyapunov exponent (MLE), *λ*_{1}, equal to zero. In this case, small changes in parameters may cause the dynamics to switch between regular and chaotic behaviour, thereby being able to *explore* large portions of the system’s phase space. This possibility is most relevant for living systems existing in fluctuating environments. In many dynamical systems, the edge of chaos exists only for a tiny portion of parameter space, typically in sets of singular points, i.e. sets of measure zero. The dynamics of systems at the edge of chaos can become highly non-trivial, even for simple maps like the logistic map (Robledo 2005). It has been argued that living systems have evolved towards the edge of chaos by natural selection (Kauffman 1993); however, it is not clear which mechanisms allow self-organization towards these exceptional regions in parameter space.

Living systems exist in the state of quasi-stationary non-equilibrium and therefore cannot be closed systems. They require a flow of both substrate and energy to and from the system. Since long ago (Lotka 1910), rate equations for molecular dynamics have been considered. For systems to be self-sustaining, such rate equations need to be autocatalytic, i.e. some molecular species directly or indirectly catalyse their own production. For living systems, cells in particular, to be in a stationary state, the production, decay and flow rates of intercellular components effectively have to balance each other (Pross & Khodorkovsky 2004; Pross 2005). Replicating, living systems therefore in general balance between stationary states (non-replicating modes) and growth (replicating mode), limited by constraints imposed by the environment. This balance provides a natural selection criterion.

Autocatalytic systems are frequently governed by nonlinear equations for enzyme kinetics, e.g. Michaelis–Menten differential equations (Michaelis & Menten 1913), or more general replicator equations (e.g. Hofbauer & Sigmund 1998). For various reasons, linearized autocatalytic networks have been considered, for the case of abundant substrate (e.g. Jain & Krishna 1998, 2002), or for reverse engineering (Yeung *et al.* 2002; Stokić *et al.* 2009). Systems with linearized dynamics can be easily depicted in terms of directed reaction networks, where nodes represent molecular species. Two nodes, where one node directly influences (production or inhibition) the other, are connected by a directed link. The weights of such links quantify associated reaction rates; negative rates indicate inhibitory links. The weights of self-loops in the reaction network, i.e. links of a node onto itself, quantify decay rates. Recent progress in genomic and proteomic technology has begun to reveal facts about regulatory networks of organisms. There is some evidence that these directed networks show scale-free topological organization (Jeong *et al.* 2000, 2001; Maslov & Sneppen 2002; Luscombe *et al.* 2004). More recent evidence suggests topological differences between in- and out-degree distributions (Dobrin *et al.* 2004; Balázi *et al.* 2005). Basically two main approaches for modelling catalytic networks have been pursued: discrete state approaches, e.g. Boolean networks (Glass & Kauffman 1973), and continuous approaches, relying on ordinary or stochastic differential equations (Smith 1987; Mahaffy *et al.* 1992; Mestl *et al.* 1996; Chen *et al.* 1999, 2005). The relevance of noise has been experimentally demonstrated (Haitzler & Simpson 1991; Ko 1992; Guptasarma 1995; Fiering *et al.* 2000; Hasty *et al.* 2000).

Various models of disordered recurrent networks (Glass & Kauffman 1973; Andrecut & Kauffman 2006) seem to share three universal modes of operation: (i) stable, (ii) critical, and (iii) chaotic super-critical. These properties could be generic or even *universal*. The importance of determining the minimum complexity of models exhibiting these properties has been pointed out (Glass & Kauffman 1973) and the question has been raised whether these properties can already be found in linear systems. Following this philosophy, we have recently introduced a model for genetic regulatory dynamics (Stokić *et al.* 2008). This model is governed by sets of linear equations,
1.1
Here *A*_{ij} is the weighted adjacency matrix of the *full* autocatalytic reaction network, whose entries may be zero, positive or negative—indicating that *i* either has no influence on *j* or the production of molecular species *i* is stimulated or suppressed by *j*, respectively. This means that if substrate *j* exists, *i* gets produced (or reduced) at rate *A*_{ij}. In (1.1) *x*_{i} is the concentration of the molecular species *i* (e.g. proteins or mRNA), and *J*_{i} corresponds to a flow vector. The molecular species *i* flows into the system if *J*_{i}>0 and out of the system if *J*_{i}<0. Finally, *ν*_{i} is a suitable noise term. Negative molecular concentration values *x*_{i} do not make sense, hence we impose the *positivity condition*
1.2
In particular, if *x*_{i}=0 and if equation (1.1) gives , then effectively the concentration *x*_{i} will remain at zero until equation (1.1) gives (which can happen through contributions from , *J*_{j} or *ν*_{i}). We refer to this as the *minimally nonlinear* (MNL) model.

MNL models have non-trivial properties (Stokić *et al.* 2008). (i) They have the possibility for chaotic dynamics. (ii) MNL exhibits an *inflated edge of chaos*. The positivity condition causes the small neighbourhood of a singular point in parameter space (linear system without positivity condition), with MLE *λ*_{1}∼0, to form an extended region (plateau). This effect gives random strategies of evolutionary phase space sampling a finite chance of locating this particular region in parameter space. This may offer an explanation for why and how complex chemical reaction systems may have found the vicinity of the edge of chaos at all, before evolutionary self-organization could take over for an eventual fine tuning. (iii) MNL models show *multi-stability*. Perturbations (or moderate noise levels) can push the system from one attractor of the dynamics to another.

These facts raise interesting questions. (i) The existence of chaotic dynamics in MNL systems straightforwardly suggests analysing the Lyapunov spectrum of the dynamics, which encodes information about the attractor of the dynamics. Numerical simulations so far have indicated that MNL systems exhibit several properties, which are particularly interesting for modelling living systems. (ii) Can topological differences between in- and out-degree distributions (Dobrin *et al.* 2004; Balázi *et al.* 2005) be explained by MNL dynamics? MNL dynamics can downregulate the concentration *x*_{i} of a fraction of nodes *i* to zero. These nodes *i* then cease to play an active role in the dynamics of the MNL system. The remaining nodes continue to play an active role in the dynamics and constitute the *active* regulatory reaction network. The active network may have topological properties that differ from those of the full network. (iii) How probable is it to find oscillating dynamics in MNL systems, and how are the fundamental frequencies of oscillatory dynamics distributed? MNL dynamics frequently shows oscillatory dynamics. This is particularly interesting, since periodic dynamics are well known in regulatory networks in the context of the cell cycle (Rustici *et al.* 2004) or circadian clocks (Alabadi *et al.* 2001). Evidence has been presented that oscillating regulatory networks are also involved in the morphogenesis of mice (Dequéant *et al.* 2006). Moreover, eukaryotic cells may encode information about extracellular environment in the frequency of stochastic intracellular events, rather than in the concentrations of molecular species (Cai *et al.* 2008). Intracellular dynamics in terms of (stochastic) rhythmic burst may be a common mechanism of intracellular information transduction.

It is clear that the proposed model is simplistic and has its limitations, the most important being that its dynamics—as given by the interaction term in equation (1.1)—is linear. Many—if not most—genetic reactions are catalytic in nature and involve ‘complexes’ (combinations of proteins are necessary before they actively stimulate genetic activity), which cannot in general be described with linear interactions. However, much of the relevance of the presented model is that it fully uncovers the essential role of the nonlinear constraint of positive concentrations. In particular, it shows which of the observed phenomena can already be explained by the constraint—and do not necessarily need to involve nonlinear dynamics in the interaction terms.

In §2, we give a summary of the model (Stokić *et al.* 2008). In §3, we report results on the properties of MNL systems. In §4, we conclude.

## 2. The stochastic minimally nonlinear model

We present the MNL model as introduced in Stokić *et al.* (2008). There we derived equation (1.1) by linearizing a set of nonlinear differential equations
2.1
where the state vector *y* represents a collection of concentrations *y*_{i} of molecular species *i*. These molecular species include both mRNA and proteins. The state vector *y* can be written as *y*=(*x*_{1},…,*x*_{n},*p*_{1},…,*p*_{m}), where the *x*_{i}, with *i*=1,…,*n*, are concentrations of mRNA and the *p*_{r}, with *r*=1,…,*m*, are protein concentrations. Equation (2.1) can be linearized around some fixed point . The variables *p*_{r} can be eliminated by the assumption that, around a fixed point, changes of mRNA (or protein) concentrations, *δx*=*x*−*x*^{0}, translate linearly into variations of the protein (or mRNA) concentrations, *δp*=*p*−*p*^{0}, i.e. , for some fixed matrix *C*. Assuming (thermal) fluctuations of the production and degradation rates around average values *A*_{ij}, using the law of large numbers, finally leads to the equation
2.2
where and *η*_{i}∈*N*(0,*σ*) are independent normally distributed zero-mean random variables with standard deviations for the multiplicative noise and *σ* for the additive noise. Comparing with equation (1.1), the noise term *ν*_{i} can be identified, , and the flow vector *J* is
2.3
Time series of mRNA expression levels *x*_{i}(*t*) typically oscillate around fixed points (average values) *x*^{0}=〈*x*_{i}(*t*)〉_{t}. This can be used to feed characteristic mRNA expression profiles directly into the MNL model. Although, from a purely mathematical point of view, the fixed point *x*^{0}=0 is a perfectly legitimate choice, it contradicts the fact that living systems are open systems and require non-vanishing effective flow vectors, *J*≠0. Equation (2.3) immediately implies that the choice *x*^{0}=0 is incompatible with this requirement. Here, we use , for all *i*.

The weighted adjacency matrix *A*_{ij} can be used to incorporate topological information on biological networks. In our model, the decay rates, *A*_{ii}<0, have identical value, for all *i*. This assumption is wrong in general but reasonable for mRNA encoding groups of proteins that act together in stoichiometric complexes (Wang *et al.* 2002). Since *A* is largely unknown experimentally, we are interested in random ensembles of matrices *A*, which can be parametrized with only a few parameters. We model *A* as a random matrix in the following way. Using terminology from network theory, the *out-degree* *k*_{i} of a node (≡ molecular species) *i* is defined as the number of products *j* that can be regulated by the molecular species *i*. The *in-degree* is the number of molecular species that regulate *i*. The ensemble of interaction networks can now be specified by the in- and out-degree *distribution* *p*_{in/out}(*k*). Although in principle in- and out-degree distributions can be chosen independently, we consider identical in- and out-degree distributions. In Stokić *et al.* (2008), we have compared scale-free networks with Erdös–Rényi networks (Erdös & Rényi 1959) and have noted only minor effects on stability, i.e. the formation of the *λ*_{1}∼0 plateau. Here, we will only consider Erdös–Rényi networks. The associated topological ensembles for the full reaction networks are completely specified by the number *N* of molecular species and the number *L* of links between them, i.e. the *average degree* 〈*k*〉=*L*/*N* of the network.

Once the topology of a network is fixed, the actual weights *A*_{ij} are sampled from a normal distribution *N*(0,*σ*_{A}) with zero mean and standard deviation *σ*_{A}. This assumption is experimentally supported by D’Haeseler *et al.* (1999). We define the constant *D*≡−*A*_{ii}/*σ*_{A} and set *σ*_{A}=1 in all numerical simulations. The time increment used for all numerical simulations is d*t*=0.1.

The MLE, *λ*_{1}, measures the exponential rate with which a perturbation of a trajectory propagates over time. If *λ*_{1}<0, the perturbation vanishes. If *λ*_{1}>0, the perturbation grows exponentially. In Stokić *et al.* (2008) it was shown that for MNL systems an interval of average network connectivities *I*=[*k*^{−},*k*^{+}] exists so that *λ*_{1}<0 for 〈*k*〉<*k*^{−} (area A), *λ*_{1}∼0 for 〈*k*〉∈*I* (area B) and *λ*_{1}>0 for 〈*k*〉>*k*^{+} (area C). The two values of 〈*k*〉 that estimate the beginning and the end of the *λ*_{1}(〈*k*〉)∼0 *plateau*, i.e. the interval *I*, are given by *k*^{−}=*D*^{2} and *k*^{+}=2*D*^{2}. It was also shown that, as 〈*k*〉 gets larger than *D*^{2}, the number *N*_{zero} of nodes *i* whose concentration *x*_{i}=0 increases monotonically, owing to the positivity condition, until *N*_{zero}∼*N*/2. The sub-network of size *N*_{on}=*N*−*N*_{zero}, consisting of those nodes *j* of the full network that have non-zero concentration *x*_{j}, we call the *active* network. If two nodes, *i* and *j*, of the active network have a link *A*_{ij} in the full network, then the active network inherits this *active* link. We denote the adjacency matrix of the active network with *A*^{on}_{ij}.

Technically, the positivity condition is implemented so that *x*_{i}(*t*+d*t*)=0, if . In Stokić *et al.* (2008), the positivity condition, equation (1.2), was implemented in a slightly different way, i.e. *x*_{i}(*t*+d*t*)=*x*_{i}(*t*), if . The different implementation has a small effect on the plateau formation.

## 3. Results

We now present (i) the Lyapunov spectrum of the MNL model, (ii) the probability to find growing, decaying or stable dynamics and (iii) the size and topological properties of the *active* catalytic network. Further, we present (iv) the probability of finding oscillating time series and their characteristic frequencies. In the following, we refer to the stable region (*k*<*k*^{−}) as area A, the plateau region [*k*^{−},*k*^{+}] as area B and the unstable region (*k*>*k*^{+}) as region C (figure 1*a*). Area B is marked by grey shading, which is also used in all the following figures to indicate the plateau region.

### (a) The Lyapunov spectrum

It is straightforward to compute the full Lyapunov spectrum *λ*_{p}, which allows us to determine properties of attractors in greater detail. In particular, we computed the *Kaplan–York dimension* *D*_{KY}, which gives an upper bound for the information dimension of the system, and the *Kolmogorov–Sinai entropy* *E*_{KS} (Eckmann & Ruelle 1985). While *D*_{KY} gives an estimate of the dimension of the attractor, i.e. the phase-space volume-preserving subspace of the dynamics, *E*_{KS} can be interpreted as a measure of the number of excited states in the system.

In figure 1, we summarize numerical results. Figure 1*a* shows the first 10 Lyapunov exponents *λ*_{p}, *p*=1,…,10, as functions of 〈*k*〉. Clearly, *λ*_{1}∼0 matches the theoretical plateau region, [*k*^{−},*k*^{+}]. Figure 1*b* shows how the full Lyapunov spectrum depends on 〈*k*〉. In area A (〈*k*〉<*k*^{−}), all *λ*_{p} are densely arranged. In area B (the plateau), *λ*_{1}∼0, while all *λ*_{p} decrease monotonically for *p*>1. As a consequence, the Lyapunov spectrum gets less dense with growing 〈*k*〉 and covers an increasing range of negative values. In area C (〈*k*〉>*k*^{+}), the Lyapunov spectrum gets still less dense. Yet, this decrease in density is qualitatively different than in area C: *λ*_{1} is increasing in area C: while *λ*_{p}, for large *p*, still decreases monotonically, but less pronounced than in area B. This can be seen clearly in figure 1*c*–*e*, where Lyapunov spectra for various values of 〈*k*〉 are shown as functions of *p*. In figure 1*c* the values of 〈*k*〉 are chosen from area A, in figure 1*d* from area B and in figure 1*e* from area C.

Further, figure 1*f*,*g* shows the average Kolmogorov–Sinai entropy and the average Kaplan–York dimension. Both quantities require the existence of some *λ*_{p}>0, i.e. *E*_{KS} and *D*_{KY} cannot be computed for area A. In areas B and C, averages of *E*_{KS} and *D*_{KY} can only be taken over realizations that have at least some *λ*_{p}>0. Clearly, in area C, the entropy per node *E*_{KS}/*N* and the fraction *D*_{KY}/*N* of the volume-preserving subspace seem to become independent of system size *N*, for sufficiently large 〈*k*〉>*k*^{+}. However, in the plateau region, area B, the two quantities do not seem to scale with system size and finite-size effects may become relevant. It is an interesting open question as to which limit function *E*_{KS}/*N* and *D*_{KY}/*N* converge towards as .

### (b) Stability of minimally nonlinear systems

What is the probability of finding the dynamics of MNL systems to be characterized by (i) exponential growth, (ii) exponential decay, or (iii) non-exponentially growing stationary or oscillatory dynamics? Living systems can be expected to exist close to stationary or oscillatory states (Pross & Khodorkovsky 2004; Pross 2005). Sufficiently positive and sufficiently negative MLEs, *λ*_{1}, in complete analogy to linear systems, indicate exponential growth or decay. Therefore, the probabilities of finding MNL systems in one of the growth modes (i)–(iii) can simply be estimated by thresholding *λ*_{1}, for a sufficiently small threshold *ε*>0. Counting realizations in the MNL ensemble with (i) *λ*_{1}≥*ε*, (ii) −*ε*≥*λ*_{1} and (iii) *ε*>|*λ*_{1}| estimates the ratios of growth mode fractions (i)–(iii). Numerical results for *ε*=0.1 are given in figure 2. Clearly, dynamics of type (iii) is favoured in the plateau region. Yet, to a much lesser extent, stationary and oscillatory dynamics also can be found in areas A and C.

### (c) The active genetic regulatory sub-network

Recent evidence from the analysis of genetic regulatory networks (Dobrin *et al.* 2004; Balázi *et al.* 2005) suggests topological differences between in- and out-degree distributions. Can differences between in- and out-degree distributions appear merely by the fact that the full interaction network *A*_{ij} is different from the *active* sub-network *A*^{on}_{ij}? Is it possible that through a symmetry-breaking mechanism the active in- and out-degree distributions *p*^{on}_{in}(*k*) and *p*^{on}_{out}(*k*) become different from *p*^{full}(*k*)?

We have analysed the average properties of active sub-networks of the MNL model and distinguish three types of nodes in MNL systems: (i) nodes *i* with concentrations *x*_{i}>0 for all times, (ii) nodes *j* with *x*_{j}=0 for all times, and (iii) nodes that alternate between on and off. The associated numbers of nodes are *N*_{pos}, *N*_{zero} and *N*_{alt}, where *N*=*N*_{pos}+*N*_{zero}+*N*_{alt}, and the fraction of nodes is denoted *n*_{pos}=*N*_{pos}/*N*, *n*_{zero}=*N*_{zero}/*N* and *n*_{alt}=*N*_{alt}/*N*. Note that the size of the active sub-network is *N*_{on}=*N*_{pos}+*N*_{alt}. In figure 3, we show these fractions for a network with *N*=500. Alternating nodes are most important in the plateau region, where *N*_{zero} starts to grow, i.e. the active networks shrink with growing 〈*k*〉. For 〈*k*〉>*k*^{+}, the fraction *n*_{alt} decreases and reaches a constant value *n*_{alt}∼0.05; *n*_{pos} and *n*_{zero} become equally large.

Figure 4*a* shows the in- and out-degree distributions of the active network for various 〈*k*〉≥*k*^{−}. The degree distribution of the full network is shown (black) for reference. Although in- and out-degree distributions of the active network differ substantially from the degree distribution of the full network, in- and out-degree distributions essentially remain identical. If we look at the weight distributions, *ρ*^{on}_{in} and *ρ*^{on}_{out}, associated with active in- and out-links in figure 4*b*, the situation changes: differences in the in- and out-weight distributions begin to show. These differences are recognizable in the standard deviation (figure 4*d*) and the skewness (figure 4*e*) of the weight distributions, but not in the mean (figure 4*c*) and the kurtosis (figure 4*f*) of the active-weight distributions. This establishes evidence that a possible symmetry breaking of in- and out-degree distributions of complete chemical reaction networks can arise owing to the natural nonlinearity in the dynamics of the chemical reactive systems, i.e. *x*_{i}≥0 for all *i*, at all times. However, the size of this effect seems to be insufficient to explain the size of topological differences (Dobrin *et al.* 2004; Balázi *et al.* 2005).

### (d) Oscillating modes in minimally nonlinear systems

What is the fraction of MNL systems that displays oscillating dynamics and what are their typical frequency distributions? We find that if a particular realization of an MNL system shows oscillating dynamics, then all *x*_{i} in the active network of the particular realization follow the same fundamental oscillation pattern. The dominant frequencies , , correspond to local maxima in the power spectra of the active *x*_{i}. Here, is the maximal number of detectable local maxima in the power spectrum of the MNL system dynamics. We looked for fundamental frequencies (if they exist).

Technically, we identified and in the following way. We computed the power spectrum *P*_{i}(*ω*) for each time series *x*_{i}(*t*) in a particular realization. We took the weighted average *g*(*ω*)≡〈*w*_{i}*P*_{i}(*ω*)〉_{i} over all nodes of the realization. The weights, *w*_{i}=*P*_{i}(*ω*_{1})^{−1}, have been chosen to be inverse proportional to the power *P*_{i}(*ω*_{1}) of the lowest frequency *ω*_{1} in the power spectrum. This choice turned out to be optimal for correctly detecting the dominant frequencies of the MNL dynamics. The first frequency can be found in the following way. We have searched for the local minimum of *g*(*ω*) with the smallest frequency . If no such local minimum exists, the time series was classified as non-oscillating. If exists, the fundamental frequency, , is determined such that is the maximum of all *g*(*ω*), with . Similarly, we computed a second dominant frequency by searching for the next local minimum , and take to be the maximum of all *g*(*ω*), with . If exists, is the second characteristic frequency of the system.

Figure 5*a* shows the probability of finding a realization of the MNL model possessing a fundamental frequency and a second dominant frequency. Oscillating realizations are dominant in the plateau region and the probability of finding oscillating realizations is close to certainty for *k*^{−}<〈*k*〉<(*k*^{−}+*k*^{+})/2. For 〈*k*〉>(*k*^{−}+*k*^{+})/2 this high probability decreases, but still has a value of about 0.6 for 〈*k*〉=*k*^{+}. Furthermore, the average (figure 5*b*) and the standard deviation (figure 5*c*) of the fundamental frequencies are shown.

## 4. Conclusions

We presented results on properties of the MNL model. We analysed the Lyapunov spectrum of the model and computed the Kolmogorov–Sinai entropy and the Kaplan–York dimension, characterizing the attractors of MNL dynamics. We analysed stability properties by computing the probabilities for finding exponentially growing, decaying and non-exponentially growing (stable) dynamics and found that stable dynamics plays a dominant role in the plateau interval [*k*^{−},*k*^{+}]. We determined characteristic fractions of concentration levels that are always downregulated to zero, are always positive or are alternating, i.e. oscillating between zero and non-zero concentration levels. Nodes with alternating concentration levels are dominating in the plateau interval. We analysed topological properties of the active regulatory network, consisting only of molecular species (nodes) with non-zero concentration levels in a given time period. We found no symmetry breaking in the in- and out-degree distributions of the active regulatory network with respect to the full network *A*_{ij}. However, we found symmetry breaking in the in- and out-weight distributions of active networks. This indicates that, in chemically reactive systems, the natural nonlinearity introduced by the positivity condition, i.e. concentrations of molecular species can never be negative, suffices to implement a symmetry breaking in the topology of the system, which can actually be measured. One may speculate if the pronounced differences of in- and out-degree distributions as found in living organisms have their origin in symmetry-breaking mechanisms, which later could become amplified by selective evolutionary processes. Finally, we determined probabilities of finding oscillating dynamics in MNL systems and analysed fundamental properties of their dominant frequencies. We found that oscillatory dynamics is most probable, in fact almost certain, for average connectivities of networks chosen from the plateau interval. This corresponds well to the observation that regulatory networks of living organisms, cells in particular, frequently show sub-networks with oscillatory dynamics. The properties analysed indicate that, near the *edge of chaos*, an MNL system—despite the simplicity of the MNL model—displays many important characteristic properties, which are expected from living matter.

## Acknowledgements

This work was supported by Austrian Science Fund FWF project P19132.

## Footnotes

One contribution of 13 to a Theme Issue ‘Complex dynamics of life at different scales: from genomic to global environmental issues’.

- This journal is © 2010 The Royal Society