Living on the edge of chaos: minimally nonlinear models of genetic regulatory dynamics
LLiving on the edge of chaos: minimally nonlinearmodels of genetic regulatory dynamics
Rudolf Hanel , Manfred P¨ochacker and Stefan Thurner , Section for Science of Complex Systems, Medical University of Vienna, Spitalgasse 23,A-1090 Vienna Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USAE-mail: [email protected]
Abstract.
Linearized catalytic reaction equations – modeling e.g. the dynamics of geneticregulatory networks – under the constraint that expression levels, i.e. molecular concentrationsof nucleic material are positive, exhibit nontrivial dynamical properties, which depend on theaverage connectivity of the reaction network. In these systems the inflation of the edge of chaos and multi-stability have been demonstrated to exist. The positivity constraint introduces anonlinearity which makes chaotic dynamics possible. Despite the simplicity of such minimallynonlinear systems , their basic properties allow to understand fundamental dynamical propertiesof complex biological reaction networks. We analyze the Lyapunov spectrum, determine theprobability to find stationary oscillating solutions, demonstrate the effect of the nonlinearity onthe effective in- and out-degree of the active interaction network and study how the frequencydistributions of oscillatory modes of such system depend on the average connectivity.
1. Introduction
Many complex systems in general – and living systems and cells in particular – display remarkablestability, i.e. a capacity to sustain their spatial and temporal molecular organization. Yet, theirstability is dynamic, i.e. these systems – to a certain degree – are capable of adapting tochanges in their physical and chemical environment. This has led several authors [1, 2, 3, 4]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 maximalLyapunov exponent (MLE), λ , equal to zero. In this case small changes in parameters maycause the dynamics to switch between regular and chaotic behavior, thereby being able to explore large portions of the system’s phase-space. This possibility is most relevant for living systemsexisting in fluctuating environments. In many dynamical systems the edge of chaos exists onlyfor a tiny portion of parameter space, typically in sets of singular points, i.e. sets of measurezero. The dynamics of systems at the edge of chaos can become highly nontrivial, even forsimple maps like the logistic map [5]. It has been argued that living systems have evolvedtowards the edge of chaos by natural selection [2], however it is not clear which mechanismsallow self-organization around these exceptional regions in parameter space.Living systems have exist in the state of quasi-stationary nonequilibrium and therefore cannot be closed systems. They require a flow of substrate and energy to and from the system.Since long [6], rate equations for molecular dynamics have been considered. For systems to beself-sustaining, such rate equations need to be autocatalytic, i.e. some molecular species directly a r X i v : . [ c ond - m a t . d i s - nn ] M a y r indirectly catalyze their own production. For living systems, cells in particular, to be in astationary state, production, decay and flow rates of intercellular components effectively haveto balance each other, [7, 8]. Replicating, living systems therefore in general balance betweenstationary states (nonreplicating modes) and growth (replicating mode), limited by constraintsposed 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 [9], or more general replicator equations, see e.g.[10]. For various reasons linearized autocatalytic networks have been considered, for the case ofabundant substrate, see e.g. [11, 12], or for reverse engineering [13, 14]. Systems with linearizeddynamics can be easily depicted in terms of directed reaction networks, where nodes representmolecular species. Two nodes, where one node directly influences (production or inhibition)the other, are connected by a directed link. Weights of such links quantify associated reactionrates; negative rates indicate inhibitory links. Weights of self-loops in the reaction network, i.e.links of a node onto itself, quantify decay rates. Recent progress in genomic and proteomictechnology begins to reveal facts about regulatory networks of organisms. There is someevidence that these directed networks show scale-free topological organization [15, 16, 17, 18].More recent evidence suggests topological differences between in- and out-degree distributions[19, 20]. Basically two main approaches for modeling catalytic networks have been pursued:Discrete state approaches, e.g. Boolean networks [21], and continuous approaches, relying onordinary or stochastic differential equations [22, 23, 24, 25, 26]. The relevance of noise has beenexperimentally demonstrated [27, 28, 29, 30, 31].Interestingly, various models of disordered recurrent networks [21, 32] seem to share threedistinguished modes of operation: (i) stable, (ii) critical, and (iii) chaotic super-critical. Theseproperties could be generic or even universal . The importance of determining the minimumcomplexity of models exhibiting these properties has been pointed out [21] and the questionhas been raised whether these properties can already be found in linear systems. Followingthis philosophy we have recently introduced a model for genetic regulatory dynamics [33]. Thismodel is governed by sets of linear equations ddt x i = (cid:88) j A ij x j + J i + ν i , (1)where A ij is the weighted adjacency matrix of the full autocatalytic reaction network, whoseentries may be zero, positive and negative – indicating that i either has no influence on j orthe production of molecular species i is stimulated or suppressed by j , respectively. This meansthat if substrate j exists, i gets produced (or reduced) at rate A ij . x i is the concentration of themolecular species i (e.g. proteins or mRNA). J i corresponds to a flow-vector. The molecularspecies i flows into the system if J i > J i < ν i is a suitable noiseterm. Negative molecular concentration values x i do not make any sense, hence we impose the positivity condition x i ≥ , for all i . (2)In particular if x i = 0 and Eq. (1) gives ˙ x i ≤
0, then effectively ˙ x i = 0. Therefore, theconcentration x i will remain zero until Eq. (1) gives ˙ x i >
0. We refer to this as the minimallynonlinear (MNL) model.MNL models have nontrivial properties [33]: (i) They have a possibility for chaoticdynamics. (ii) MNL exhibit an inflated edge of chaos . The positivity condition causes the smallneighborhood of a singular point in parameter space (linear system without positivity condition),with MLE λ ∼
0, to form an extended region (plateau). This effect gives random strategies ofevolutionary phase-space sampling a finite chance of locating this particular region in parameterspace. This may offer an explanation for why and how complex chemical reaction systems mayave found the vicinity of the edge of chaos at all, before evolutionary self-organization couldtake over for an eventual fine tuning. (iii) MNL models show multi-stability . Perturbations (ormoderate noise levels) can push the system from one attractor of the dynamics to another.These facts raise interesting questions. (a) The existence of chaotic dynamics in MNL systemsstraight forwardly suggests to analyze the Lyapunov spectrum of the dynamics, which encodesinformation about the attractor of the dynamics. Numerical simulations indicated so far thatMNL systems exhibit several properties, which are particularly interesting for modeling livingsystems. (b) Can topological differences between in- and out-degree distributions [19, 20], beexplained by MNL dynamics? MNL dynamics can down-regulate the concentration x i of afraction of nodes i to zero. These nodes i then cease to play an active role in the dynamicsof the MNL system. The remaining nodes continue to play an active role in the dynamicsand constitute the active regulatory reaction network. The active network may have topologicalproperties that differ from the full network. (c) How probable is it to find oscillating dynamics inMNL systems, and how are fundamental frequencies of oscillatory dynamics distributed? MNLdynamics frequently shows oscillatory dynamics. This is particularly interesting, since periodicdynamics are well known in regulatory networks in the context of the cell-cycle, e.g. [34], orcircadien clocks [35]. Evidence has been presented that oscillating regulatory networks are alsoinvolved in the morphogenesis of mice [36]. Moreover, eukaryotic cells may encode informationabout extracellular environment in the frequency of stochastic intracellular events, rather thanin the concentrations of molecular species [37]. Intracellular dynamics in terms of (stochastic)rhythmic burst, may be a common mechanism of intracellular information transduction .In Section 2 we give a summary of the model [33]. In Section 3 we report results on theproperties of MNL systems. In Section 4 we conclude.
2. The stochastic MNL model
We present the MNL model as introduced in [33]. There we derived Eq. (1) by linearizing a setof nonlinear differential equations ddt y i = F i ( y ) , (3)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 , . . . x n , p , . . . p m ), where the x i , with i = 1 . . . m , are concentrations of mRNA and the p r , with r = 1 . . . m , are protein concentrations. Equation (3) can be linearized around somefixed point y = ( x , . . . x n , p , . . . p m ). The variables p r can be eliminated by the assumptionthat, around a fixed point, changes of mRNA (or protein) concentrations, δx = x − x ,translate linearly into variations of the protein (or mRNA) concentrations, δp = p − p , i.e. δp r = (cid:80) i C ri δx i , for some fixed matrix C . Assuming (thermal) fluctuations of the productionand degradation rates around average values A ij , using the law of large numbers, finally leadsto the equation ddt x i = (cid:88) j A ij ( x j − x j ) + ξ i ( t )( x i − x i ) + η i , (4)where ξ i ∈ N (0 , ¯ σ ) and η i ∈ N (0 , σ ) are independent normally distributed zero-mean randomvariables with standard deviations ¯ σ for the multiplicative noise and σ for the additive noise.Comparing with Eq. (1) the noise-term ν i can be identified, ν i = ξ i ( t )( x i − x i ) + η i , and theflow-vector J is J i = − (cid:88) j A ij x j . (5)Time series of mRNA expression levels x i ( t ) typically oscillate around fixed points (averagevalues) x = (cid:104) x i ( t ) (cid:105) t . This can be used to directly feed characteristic mRNA expression profilesnto the MNL model. Although from a purely mathematical point of view, the fixed point x = 0is a perfectly legitimate choice, it contradicts the fact that living systems are open systems andrequire non-vanishing effective flow-vectors, J (cid:54) = 0. Equation (5) immediately implies that thechoice x = 0 is incompatible with this requirement. Here we use x i = 1000, for all i .The weighted adjacency matrix A ij can be used to incorporate topological information onbiological 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 proteinsthat act together in stoichiometric complexes [38]. Since A is largely unknown experimentallywe are interested in random ensembles of matrices A , which can be parametrized with only afew parameters. We model A as a random matrix in the following way. Using terminologyfrom network theory, the out-degree k i of a node ( ≡ molecular species) i is defined as thenumber of products j that can be regulated by the molecular species i . The in-degree is thenumber of molecular species that regulate i . The ensemble of interaction networks can nowbe specified by the in- and out- degree distribution p in/out ( k ). Although in principle in- andout-degree distribution can be chosen independently, we consider identical in- and out-degreedistributions. In [33] we have compared scale-free networks with Erd¨os-R´enyi networks [39]and have noted only minor effects on the stability, i.e. the formation, of the λ ∼ N of molecular species and thenumber L of links between them, i.e. the average degree (cid:104) k (cid:105) = L/N of the network.Once the topology of a network is fixed, the actual weights A ij are sampled from anormal distribution N (0 , σ A ) with zero mean and standard deviation σ A . This assumptionis experimentally supported by e.g. [40]. We define the constant D ≡ − A ii /σ A and set σ A = 1in all numerical simulations. The time-increment used for all numerical simulations is dt = 0 . λ , measures the exponential rate with which aperturbation of a trajectory propagates over time. If λ <
0, the perturbation vanishes. If λ >
0, the perturbation grows exponentially. In [33] it was shown that for MNL systems aninterval of average network connectivities I = [ k − , k + ] exists so that λ < (cid:104) k (cid:105) < k − (AreaA), λ ∼ (cid:104) k (cid:105) ∈ I (Area B), and λ > (cid:104) k (cid:105) > k + (Area C). The two values of (cid:104) k (cid:105) , whichestimate the beginning and the end of the λ ( (cid:104) k (cid:105) ) ∼ plateau , i.e. the interval I , are given by k − = D and k + = 2 D . It also was shown, that as (cid:104) k (cid:105) gets larger than D , the number N zero of nodes i , whose concentration x i = 0 increases monotonically, due 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 thefull network, which have nonzero 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 inheritsthis 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 + dt ) = 0, if x i ( t )+ ˙ x i ( t ) dt ≤ x i ( t + dt ) = x i ( t ), if x i ( t ) + ˙ x i ( t ) dt ≤
0. The different implementation has a small effecton the plateau formation.
3. Results
We now present (i) the Lyapunov spectrum of the MNL model, (ii) the probability to findgrowing, decaying, or stable dynamics and the size and topological properties of the active catalytic network. Further we present (iv) the probability of finding oscillating time series andtheir characteristic frequencies. In all following figures the theoretical plateau interval [ k − , k + ]is marked (gray shading). igure 1. (a) Largest ten Lyapunov exponents ( λ p , p = 1 , . . . ,
10) of the Lyapunov spectrum( N = 500). The inset magnifies the plateau region. The two black dashed lines are theoreticalcurves, derived in [33], approximating λ ( (cid:104) k (cid:105) ) in the areas A and C. The intersection of thesecurves with the x-axis, λ ( (cid:104) k (cid:105) ) = 0, estimate the beginning and end of the λ ( (cid:104) k (cid:105) ) ∼ (cid:104) k (cid:105) is shown ( N = 200). (c,d,e) Lyapunovspectra, λ p ( (cid:104) k (cid:105) ) as functions of p , for areas A, B, and C. In area A spectra are all similar, inarea B the slope of λ p gets steeper with growing (cid:104) k (cid:105) , while λ ( (cid:104) k (cid:105) ) ∼
0. Comparing (d) and(e) shows a clear qualitative difference between area B and C. The Kolmogorov Sinai Entropy E KS , diagram (f), and the Kaplan York Dimension D KY , diagram (g), both divided by systemsize N , in logarithmic scale, for N = 200 and N = 500. D KY represents an upper bound for theinformation dimension of the system. E KS is calculated from Pesin’s theorem, as the sum overall positive λ p > (cid:104) k (cid:105) .In (a-g) averages are taken over 50 random realizations, time interval [200 , σ = ¯ σ = 0 . σ A = 1, D = 4, and x = 1000. .1. The Lyapunov spectrum It is straight forward to compute the full Lyapunov spectrum λ p which allows to determineproperties 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 , see e.g. [41]. While D KY gives an estimate of the dimensionof the attractor, i.e. the phase-space volume-preserving subspace of the dynamics, E KS can beinterpreted as a measure of the number of excited states in the system.In Fig. (1) we summarize numerical results. (a) shows the first ten Lyapunov exponents λ p , p = 1 . . .
10, as functions of (cid:104) k (cid:105) . Clearly, λ ∼ k − , k + ]. (b) shows how the full Lyapunov spectrum depends on (cid:104) k (cid:105) . In area A ( (cid:104) k (cid:105) < k − ) all λ p are densely arranged. In area B (the plateau) λ ∼
0, while all λ p decrease monotonically for p >
1. As a consequence, the Lyapunov spectrum gets less dense with growing (cid:104) k (cid:105) and coversan increasing range of negative values. In area C ( (cid:104) k (cid:105) > k + ) the Lyapunov spectrum gets stillless dense. Yet, this decrease in density is qualitatively different than in area C. λ is increasingin area C, while λ p , for large p , still decreases monotonically, but less pronounced than in areaB. This can be seen clearly in Fig. (1) (c-e), where Lyapunov spectra for various values of (cid:104) k (cid:105) are shown as functions of p . In (c) the values of (cid:104) k (cid:105) are chosen from area A, in (d) from area B,and in (e) from area C.Further, Fig. (1) shows the average Kolmogorov-Sinai entropy (f) and the average Kaplan-York dimension (g). Both quantities require the existence of some λ p >
0, i.e. E KS and D KY cannot be computed for area A. In area B and C averages of E KS and D KY can onlybe taken over realizations, which have at least some λ p >
0. Clearly, in area C the entropyper node E KS /N and the fraction D KY /N of the volume-preserving subspace seem to becomeindependent of system-size N , for sufficiently large (cid:104) k (cid:105) > k + . However, in the plateau region,area B, both quantities do not seem to scale with system size and finite size effects may becomerelevant. It is an interesting open question towards which limit-function E KS /N and D KY /N converge as N → ∞ . Figure 2.
Fraction of realizations which lead to exponentially growing ( λ > . λ < − .
1) and stable time series ( | λ | ≤ .
1) computed from 100 realizations, N = 500, D = 4,time interval, [200 , σ = ¯ σ = 0 . σ A = 1, and x =1000.What is the probability of finding the dynamics of MNL systems to be characterized (i)y exponential growth, (ii) exponential decay, or (iii) non-exponentially growing stationary oroscillatory dynamics? Living systems can be expected to exist close to stationary or oscillatorystates [7, 8]. Sufficiently positive and sufficiently negative MLEs, λ , in complete analogy tolinear systems, indicate exponential growth or decay. Therefore the probabilities of finding MNLsystems in one of the growth modes (i-iii) can simply be estimated by thresholding λ , for asufficiently small threshold ε >
0. Counting realizations in the MNL ensemble with (i) λ ≥ ε ,(ii) − ε ≥ λ , and (iii) ε > | λ | estimates the ratios of growth mode fractions (i-iii). Numericalresults for ε = 0 . Figure 3.
Average fractions of x i : positive ( n pos ), zero ( n zero ), or alternating ( n alt ). Averagesare taken over 1000 realizations, time interval, [500 , N = 500, D = 4, σ = ¯ σ = 0, σ A = 1, x = 1000.Recent evidence from the analysis of genetic regulatory networks [19, 20], suggests topologicaldifferences between in- and out-degree distributions. Can differences between in- and out-degreedistributions appear merely by the fact that the full interaction network A ij is different fromthe active sub-network A on ij ? Is it possible that through a symmetry-breaking mechanism theactive in- and out-degree distribution p onin ( k ) and p onout ( k ) become different from p full ( k )?We have analyzed the average properties of active sub-networks of the MNL model anddistinguish three types of nodes in MNL systems: (i) nodes i with concentrations x i > 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 andthe fraction of nodes are denoted n pos = N pos /N , n zero = N zero /N , and n alt = N alt /N . Notethat the size of the active sub-network is N on = N pos + N alt . In Fig. (3) we show these fractionsfor 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 shrinks with growing (cid:104) k (cid:105) . For (cid:104) k (cid:105) > k + the fraction n alt decreases and reaches a constant value n alt ∼ . n pos and n zero become equally large.In Fig. (4 (a) unweighted in- & out-degree shows the in- and out-degree distributions of theactive network for various (cid:104) k (cid:105) ≥ k − . The degree distribution of the full network is shown (red) igure 4. (a) Unweighted in- and out-degree distributions of the active regulatory sub-networkfor various (cid:104) k (cid:105) . Active in- and out-degree, p onin / out ( k ), are practically indistinguishable. (b)Weighted in- and out-degree distributions. In- and out-weight distributions, ρ onin / out ( φ ), of activeweights are clearly distinguishable. φ = (cid:80) A ij and the sum runs over i or j for in- and out-weightdistribution, respectively. (c) Mean, (d) standard deviation, (e) skewness and (f) kurtosis of thein/out-weight distributions. Differences between in- and out-weight distributions are found inthe standard deviation and the skewness. Averages are taken over 50 realizations, N = 500,time interval, [500 , D = 4, σ = ¯ σ = 0 . σ A = 1, x = 1000.or reference. Although in- and out-degree distribution of the active network differ substantiallyfrom the degree distribution of the full network, in- and out-degree distributions essentiallyremain identical. If we look at the weight distributions, ρ onin and ρ onout , associated with active in-and out-links in (b) weighted in- & out-degree the situation changes: differences in the in- andout-weight distributions begin to show. These differences are recognizable in (d) the standarddeviation and (e) the skewness of the weight distributions, but not in (c) the mean and (f) thekurtosis of the active weight distributions. This establishes evidence that a possible symmetrybreaking of in- and out-degree distributions of complete chemical reaction networks can arisedue to the natural nonlinearity in the dynamics of the chemical reactive systems, i.e. x i ≥ i , at all times. However, the size of this effect seems to be insufficient to explain the size oftopological differences [19, 20]. Figure 5. (a) Probabil-ity of finding oscillatingrealizations: with ex-isting fundamental fre-quency ω ∗ (blue) andboth, existing ω ∗ and ω ∗ (green). (b) Av-erage ω ∗ as a func-tion of (cid:104) k (cid:105) . (c) Stan-dard deviation of ω ∗ . N = 500, time inter-val, [1000 , D =4, σ = ¯ σ = 0, σ A =1, x = 1000. In (b)and (c) N = 500, D =6, (green circles) and N = 200, D = 4 (redsquares) are shown forcomparison.What is the fraction of MNL systems, which display oscillating dynamics and what are theirtypical frequency distributions? We find that if a particular realization of an MNL system showsoscillating dynamics, then all x i in the active network of the particular realization follow the samefundamental oscillation-pattern. The dominant frequencies ω ∗ s , s = 1 , , . . . s max , correspond toocal maxima in the power-spectra of the active x i . s max ≤ N is the maximal number ofdetectable local maxima in the power-spectrum of the MNL system dynamics. We looked forfundamental frequencies ω ∗ (if existing).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 ( ω ) ≡ (cid:104) w i P i ( ω ) (cid:105) i over all nodes of the realization. The weights, w i = P i ( ω ) − , have beenchosen inverse-proportional to the power P i ( ω ) of the lowest frequency ω in the power-spectrum. This choice turned out to be optimal for correctly detecting the dominant frequencies ω ∗ s of MNL dynamics. The first frequency ω ∗ can be found in the following way. We have searchedfor the local minimum of g ( ω ) with the smallest frequency ˆ ω . If no such local minimum existsthe time series was classified as non-oscillating. If ˆ ω exists, the fundamental frequency, ω ∗ > ˆ ω ,is determined such that g ( ω ∗ ) is the maximum of all g ( ω ), with ω > ˆ ω . Similarly, we computeda second dominant frequency by searching for the next local minimum ˆ ω > ω ∗ , and take g ( ω ∗ )to be the maximum of all g ( ω ), with ω > ˆ ω . If ˆ ω exists, ω ∗ is the second characteristicfrequency of the system.In Fig. (5) (a) shows the probability of finding a realization of the MNL model possessing afundamental frequency and a second dominant frequency. Oscillating realizations are dominantin the plateau region and the probability of finding oscillating realizations is close to certaintyfor k − < (cid:104) k (cid:105) < ( k − + k + ) /
2. For (cid:104) k (cid:105) > ( k − + k + ) / . (cid:104) k (cid:105) = k + . Furthermore, in Fig. (5) (b) the average and (c) thestandard-deviation of the fundamental frequencies ω ∗ are shown.
4. Conclusions
We presented results on properties of the MNL model. We analyzed the Lyapunov spectrumof the model and computed the Kolmogorov-Sinai Entropy and the Kaplan-York Dimension,characterizing the attractors of MNL dynamics. We analyzed stability properties by computingthe 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, which are always down-regulated 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 theplateau interval. We analyzed topological properties of the active regulatory network, consistingonly of molecular species (nodes) with nonzero concentration levels in a given time period. Wefound no symmetry-breaking in the in- and out-degree distributions of the active regulatorynetwork with respect to the full network A ij . However, we found symmetry-breaking in thein- and out-weight distributions of active networks. This indicates that in chemically reactivesystems the natural nonlinearity introduced by the positivity condition, i.e. concentrationsof molecular species can never be negative, suffices to implement a symmetry-breaking in thetopology of the system, which can actually be measured. One may speculate if the pronounceddifferences of in- and out-degree distributions as found in living organisms, have their origin insymmetry-breaking mechanisms, which later could become amplified by selective evolutionaryprocesses. Finally, we determined probabilities of finding oscillating dynamics in MNL systemsand analyzed fundamental properties of their dominant frequencies. We found that oscillatorydynamics is most likely, in fact almost certain, for average connectivities of networks chosen fromthe plateau interval. This corresponds well to the observation that regulatory networks of livingorganisms, cells in particular, frequently show sub-networks with oscillatory dynamics. Theproperties analyzed indicate that near the edge of chaos MNL system – despite the simplicityof the MNL model – display many important characteristic properties, which are expected fromliving matter. cknowledgments
Supported by Austrian Science Fund FWF project P19132.
References [1] Langton C 1990
Physica D The Origins of Order: Self-Organization and Selection in Evolution [3] Mitchell M, Hraber P and Crutchfield J 1993
Complex Systems Dynamic Patterns in Complex Systems
Europhys. News J. Phys. Chem. J. Phys. Org. Chem. Pure Appl. Chem. Biochem. Z. Evolutionary Games and Population Dynamics [11] Jain S and Krishna S 2002
Proc. Natl. Acad. Sci. Phys. Rev. Lett. Proc. Natl. Acad. Sci. BMC Bioinformatics Science
Nature
Nature
Nature
BMC Bioinformatics PNAS
J. Theor. Biol. J. Math. Biol J. Math. Biol. Physica D Pac. Symp. on Biocomputing Bioinformatics Bioessays Bioessays Proc. Natl. Acad. Sci. Cell Bioessays New J. Phys. Phys. Rev. E Nature Genetics Science
Science
Nature
PNAS Publicationes Mathematicae Pac. Symp. on Biocomputing Reviews of Modern Physics57