text:
# Tools for network dynamicsSupported by Fundação para a Ciência e Tecnologia

## 1 Introduction

When modeling extended complex systems, the network concept appears quite often. The metabolic processes of living beings are a network with the substrates as nodes, linked together whenever they participate in the same biochemical reaction. Protein-protein as well as gene expression and regulation are also networks. Social, economic and political networks are the backbone of human society, the internet is a network, etc.\[Albert & Barabási, 2002\] \[Dorogovtsev & Mendes, 2003\]. Most studies deal with networks as statistical objects, with extensive use of the tools of statistical mechanics \[Pastor-Satorras et al., 2003\]. Much less attention has been paid to the dynamical phenomena taking place in the networks or to the behavior of the evolving networks as dynamical systems.

For several decades an intensive effort has been dedicated to the study of low-dimensional dynamical systems, leading to an extensive body of rigorous results \[Katok & Hasselblatt, 1995\]. This exploration is still proceedings at a good pace with exciting new results, for example, in the dimension theory of dynamical systems \[Pesin, 1997\] \[Barreira, 2002\] and non-uniform hyperbolicity \[Bonatti et al., 2003\]. However, the main challenge for physical applications lies on extended systems, in particular on the understanding of the dynamics leading from microscopic laws to global patterns \[Cross & Hohenberg, 1993\]. A large amount of numerical work has been done on the dynamics of these systems which led to the identification and classification of typical patterns, spatio-temporal chaos, statistical properties and multistability \[Kaneko, 1993\] \[Kaneko & Tsuda, 2000\] \[Boccaletti et al., 2001\] \[Boldrighini et al., 2000\]. Rigorous results are few, except for regular coupled map lattices \[Bunimovich & Sinai, 1988\] \[Coutinho & Fernandez, 1997a, 1997b\] \[Jiang & Pesin, 1998\] \[Afraimovich & Fernandez, 2000\] \[Gielis & MacKay, 2000\] \[Fernandez & Guiraud, 2004\].

For non-uniform coupling structures there are much less results, statistical mechanics tools being used mostly to characterize the topological features of static and evolving networks. Of course the topological structure of the network is very important for the dynamics. Form affects function and topology controls the rate at which information or diseases propagate \[Boots & Sasaki, 1999\] \[Keeling, 1999\] \[Pastor-Satorras & Vespignani, 2002\] \[Lloyd & May, 2001\], the robustness under attack and failure \[Albert et al., 2000\] or the adaptation and learning processes that take place \[Araújo & Vilela Mendes, 2000\].

The main purpose of this review is to provide a toolkit for the treatment of networks (both regular and irregular) as dynamical systems. Results from differential dynamics and ergodic theory will be presented. To deal with a network as a dynamical system, three main problems have to be addressed. First, how to characterize the dynamical behavior of the “agents” sitting on the nodes with interactions defined by the link structure of the network. Second, if the network topology is not fixed either because the links are changing in time or because the network itself is growing, to characterize the network evolution as a dynamical system. Finally to characterize the interplay between the network topology and the nature of the dynamics.

Section 2 concerns the description of network dynamics by global functions, limit cycles and multistability, with some examples illustrating both dynamics on a network and network evolution as a dynamical system. Section 3 discusses ergodic techniques (old and new) and their role in networks. Section 4 discusses the logic approach to network dynamics. On the important problem of topology versus dynamics, I will touch only briefly in the ergodic section, referring for a particular detailed example to \[Araújo et al., 2003\].

In almost all cases the reader should refer to the original papers for proofs and further developments. Occasionally, as for example in the family of ergodic parameters or in the relation between synchronization and Lyapunov spectrum, somewhat more detail is given to the exposition. This does not mean in any way that I consider this to be a more important topic than the other interesting results of many authors that are mentioned. If more detail is included it is because I believe it to be a new result not included before in any other publication.

## 2 Differential dynamics tools

### 2.1 Describing dynamics by global functions

In many networks, the node dynamics may be modelled by ordinary differential equations of the form

$$\frac{dx_i}{dt}=X_i\left(x\right)=\alpha _i+\underset{ji}{}W_{ij}f\left(x_j\right)\gamma _ix_i$$
(1)
For a neural network, the $`x_i^{}s`$ might be firing rates and the $`W_{ij}^{}s`$ synaptic intensities \[Grossberg, 1988\], for a genetic regulatory system \[Tyson & Ohtmer, 1978\] \[de Jong, 2002\] the variables $`x_i`$ would code for the concentrations of RNA, proteins or other metabolic components and $`W_{ij}`$ for the production constants (measuring the strength of $`j`$ on $`i`$), $`f()`$ being the regulation function and $`\gamma _ix_i`$ a degradation term, etc.

#### 2.1.1 Symmetric systems

Eq.(1) is a particular case of the Cohen-Grossberg form \[Cohen & Grossberg, 1983\], used by these authors to describe continuous-time neural networks,

$$\frac{dx_i}{dt}=a_i(x_i)\left\{b_i(x_i)\underset{j=1}{\overset{n}{}}W_{ij}f_j(x_j)\right\}$$
(2)
Cohen and Grossberg proved that, for the symmetric case ($`W_{ij}=W_{ji}`$), the following function

$$V\left(x_i\right)=\underset{i=1}{\overset{n}{}}^{x_i}b_i(\xi _i)f_i^{^{}}(\xi _i)𝑑\xi _i+\frac{1}{2}\underset{j,k=1}{\overset{n}{}}W_{jk}f_j(x_j)f_k(x_k)$$
(3)
is a Lyapunov function, that is

$$\frac{d}{dt}V\left(x_i\right)0$$
(4)
along the orbits if $`a_i(x_i)f_i^{^{^{}}}(x_i)>0`$. Hopfield’s \[Hopfield, 1984\] “energy” function is a particular case of this result.

The existence of a Lyapunov function is a useful device to characterize the asymptotically stable states of the network or for the synthesis of networks with a desired number of stable asymptotic solutions \[Cohen, 1992\].

In the case of symmetric connections the continuous-time result of Cohen and Grossberg has been extended to a class of discrete-time systems (\[Fogelman Soulié et al., 1989\] and references therein). For non-symmetric connections of particular form, namely

$$\mu _jW_{ij}=\mu _iW_{ji}$$
(5)
$`\mu _i>0`$, and time evolution of the connection strengths of Hebbian type

$$\frac{d}{dt}W_{ij}=\gamma _{ij}W_{ij}+f_i\left(x_i\right)f_j\left(x_j\right)$$
(6)
in Schürmann  or $`W_{ij}W_{ji}>0`$ and $`_CW_{ij}=_CW_{ji}`$ along every cycle in Fiedler & Gedeon , Lyapunov functions may also be constructed.

#### 2.1.2 General systems

The Cohen-Grossberg result has been generalized for arbitrary $`w_{ij}^{}s`$ in Vilela Mendes & Duarte , namely given

$$\begin{array}{ccc}W_{ij}& =& W_{ij}^{(S)}+W_{ij}^{(A)}\hfill \\ W_{ij}^{(S)}& =& \frac{1}{2}\left(W_{ij}+W_{ji}\right)\hfill \\ W_{ij}^{(A)}& =& \frac{1}{2}\left(W_{ij}W_{ji}\right)\hfill \\ V^{(S)}& =& _{i=1}^n^{x_i}b_i(\xi _i)f_i^{^{}}(\xi _i)𝑑\xi _i+\frac{1}{2}_{j,k=1}^nW_{jk}^{(S)}f_j(x_j)f_k(x_k)\hfill \\ H& =& _{i=1}^n^{x_i}\frac{f_i(\xi _i)}{a_i(\xi _i)}𝑑\xi _i\hfill \end{array}$$
(7)
one has the following

Theorem \[Vilela Mendes & Duarte, 1992\] If $`a_i(x_i)/f_i^{^{^{}}}(x_i)>0`$ $`x,i`$ and the matrix $`W_{ij}^{(A)}`$ has an inverse, the vector field $`\stackrel{}{x_i}`$ in Eq.(2) decomposes into one gradient and one Hamiltonian component, $`\stackrel{}{x_i}=\stackrel{(G)}{\stackrel{}{x_i}}+\stackrel{(H)}{\stackrel{}{x_i}}`$, with

$$\begin{array}{ccc}\stackrel{(G)}{\stackrel{}{x_i}}& =& \frac{a_i(x_i)}{f_i^{^{}}(x_i)}\frac{V^{(S)}}{x_i}=_jg_{ij}(x)\frac{V^{(S)}}{x_j}\hfill \\ \stackrel{(H)}{\stackrel{}{x_i}}& =& _ja_i(x_i)w_{ij}^{(A)}(x)a_j(x_j)\frac{H}{x_j}=_j\mathrm{\Gamma }_{ij}(x)\frac{H}{x_j}\hfill \end{array}$$
(8)
and

$$\begin{array}{ccc}g_{ij}(x)& =& \frac{a_i(x_i)}{f_i^{^{}}(x_i)}\delta _{ij}\hfill \\ \omega _{ij}(x)& =& a_i(x_i)^1\left(W^{(A)1}\right)_{ij}(x)a_j(x_j)^1\hfill \end{array}$$
(9)
$`\left(_j\mathrm{\Gamma }_{ij}\omega _{jk}=\delta _{ik}\right)`$. $`g_{ij}(x)`$ and $`\omega _{jk}(x)`$ are the components of the Riemannian metric and the symplectic form.

———

The conditions on $`a_i(x_i)`$, $`f_i^{^{^{}}}(x_i)`$ and $`w_{ij}^{(A)}`$ insure that $`g`$ is a well defined metric and that $`\omega `$ is non-degenerate.

The decomposition (8) is useful, for example, on the design of oscillatory networks and on the study of gated learning rules \[Howse et al., 1996\]. The nature of the dynamics in the network will depend on the relative strength of the gradient and the Hamiltonian components. Howse et al.  propose to measure this relative strength by comparing $`\frac{dV^{(S)}}{dt}`$ with $`\frac{dH}{dt}`$. However, these quantities vary in space and time and it is the compensation of the two effects that in particular regions of phase space lead to the attractors of the dynamics, for example to limit cycles (see below).

The identification, in the differential system (2), of just one gradient and one Hamiltonian component, with explicitly known potential and Hamiltonian functions, is a considerable simplification as compared to a generic dynamical system. For a general dynamical system a representation by one or two functions is possible only locally \[Vilela Mendes & Duarte, 1981\] and explicit forms for the functions are not easy to obtain \[Abarbanel & Rouhi, 1987\] \[Crehan, 1994\]. Global decomposition for general dynamical systems require one gradient and $`n1`$ Hamiltonian components \[Vilela Mendes & Duarte, 1981\], namely

$$\stackrel{}{x_i}=\underset{j}{}g_{ij}(x)\frac{V}{x_j}+\underset{k=1}{\overset{n1}{}}\underset{j=1}{\overset{n}{}}\left(\omega _{(k)}^1(x)\right)_{ij}\frac{H^{(k)}}{x_j}$$
(10)
$`\left\{\omega _{(k)}(x)\right\}`$ being a set of canonical symplectic forms adapted to each Hamiltonian component. This result is a generalization to $`n`$ dimensions of the 2-dimensional result of Roels . The first term in (10) is the dissipative component and the second one corresponds to a volume-preserving dynamical system.

The above results lead to a convenient characterization of dynamical systems of type (1) or (2). For the symmetric case the existence of a Lyapunov function guarantees global asymptotic stability of the dynamics. However not all vector fields with a Lyapunov function are differentially equivalent to a gradient field. Therefore the fact that a gradient vector is actually obtained gives additional information, namely about structural stability of the model. A necessary condition for structural stability of the gradient vector field is the non-degeneracy of the critical points of $`V^{(S)}`$, namely $`det\frac{^2V^{(S)}}{x_ix_j}0`$ at the points where $`\frac{V^{(S)}}{x_i}=0`$. In a gradient flow all orbits approach the critical points as $`t\mathrm{}`$. If the critical points are non-degenerate, the gradient flow satisfies the conditions defining a Morse-Smale field, except perhaps the transversality conditions for stable and unstable manifolds of the critical points. However because Morse-Smale fields are open and dense in the set of gradient vector fields, any gradient flow with non-degenerate critical points may always be C1-approximated by a (structurally stable) Morse-Smale gradient field. Therefore given a symmetric model of the type (2), the identification of its gradient nature provides a easy way to check its robustness as a physical model.

Although Lyapunov functions may in some cases be constructed for discrete-time systems \[Fogelman Soulié et al.,1989\], the natural functional representation of maps is through generating functions. This is well known for canonical maps of symplectic manifolds \[Amiet & Huguenin, 1980\] and has been generalized in \[Vilela Mendes, 1986\] for noncanonical maps.

The representation of network dynamics by global function applies to neural networks of several types \[Grossberg, 1988\], to more general networks \[Chua, 1988a, 1988b\] and, in view of an established correspondence \[Doyne Farmer, 1990\], to a large range of connectionistic systems.

### 2.2 Cycles

Existence of limit cycle oscillations in networks is an important issue \[Gouzé, 1998\] \[Plahte et al., 1995\] \[Snoussi, 1998\]. The decomposition theorems provide a tool to look for candidate orbits with limit cycle properties. Many years ago Pontryagin , studying small perturbations of Hamiltonian fields on the plane

$$\stackrel{}{x}=\frac{H}{y}+\epsilon A(x,y,\epsilon ),\stackrel{}{y}=\frac{H}{x}+\epsilon B(x,y,\epsilon )$$
(11)
introduced the notion of generating cycle $`\gamma \left(c\right)`$, lying on a level curve $`H=c`$, when the perturbed equation has a cycle that depends continuously on $`\epsilon `$, for small $`\left|\epsilon \right|`$, and tends to $`\gamma \left(c\right)`$ when $`\epsilon 0`$. Pontryagin’s result states that if $`\gamma \left(c\right)`$ is a generating cycle, then

$$I\left(c\right)=_{\gamma \left(c\right)}\left(BdxAdy\right)=0$$
(12)
the integration being along $`\gamma \left(c\right)`$ at $`\epsilon =0`$.

Further results on the existence of cycles were later proved both for weakly coupled oscillators and for more general systems with parametrized families of solutions (see \[Hoppensteadt & Izhikevich, 1997\], chapter 9 and references therein). A generalization of Pontryagin’s result to dynamical systems with constants of motion \[Duarte & Vilela Mendes, 1983\], leads to a necessary condition for the existence of a cycle using the decomposition in (10), namely

$$\left\{\left(H^{(i)}V\right)+\underset{ki}{}\omega _{(k)}\left(H^{(i)}H^{(k)}\right)\right\}𝑑\gamma _i=0$$
(13)
the integration being along a closed level curve $`\gamma _i`$ of $`H_i`$.

A similar result holds for discrete-time maps which belong to a differentiable arc with constants of motion \[Vilela Mendes & Duarte, 1982\]. A constant of motion for a map $`f`$ defined on a manifold $`M`$ is a differentiable function $`\mathrm{\Phi }:MR`$ such that for some orbit $`\gamma `$, $`\mathrm{\Phi }\gamma =`$constant. It generalizes the notion of first integral which would require this to hold for all orbits. A family of maps $`f_\epsilon `$ is called a differentiable arc with constants of motion if (i) each $`f_\epsilon `$ has a constant of motion $`\mathrm{\Phi }_\epsilon `$ for some orbit $`\gamma _\epsilon `$ ; (ii) The constant of motion $`\mathrm{\Phi }_0`$ of $`f_0`$ is a first integral in a neighborhood of $`\gamma _0`$ ; (iii) the maps $`\epsilon f_\epsilon ,\epsilon \gamma _\epsilon ,\epsilon \mathrm{\Phi }_\epsilon `$ are differentiable. Then

$$\underset{n=0}{\overset{N_01}{}}D\mathrm{\Phi }_0\left(\gamma _0\left(n+1\right)\right)|f^{^{}}\left(\gamma _0\left(n\right)\right)=0$$
(14)
$`N_0`$ being the period of the orbit $`\gamma _0`$.

Both (13) and (14) give only necessary equations for the existence of limit cycles in the composite dynamics. Nevertheless they are useful tools to identify limit cycle candidates. Sufficient conditions may also be obtained in particular low-dimensional cases \[Vilela Mendes, 1988, 2000a\].

In the same way as the Hamiltonian components of the dynamics provide a tool to look for limit cycles, the stationary points of the gradient potential provide information on the multistability of the dynamics and the nature of their basins of attraction. It is also a tool for the construction of the invariant measures of the dynamics (see below).

### 2.3 Multistability

Existence of multiple stable states with distinct basins of attraction plays a significant role in the dynamics of networks, for example in those associated to the basic processes of life. A genetic regulatory network with different stable patterns of gene activation explains the emergence of different phenotypic expressions in the absence of genetic differences \[Laurent & Kellershohn, 1999\]. Examples are also found in population dynamics \[Henson, 2000\], neural dynamics \[Skarda & Freeman, 1987\] \[Freeman, 1992\], geophysics, etc. \[Vilela Mendes, 2000a\]. Extensive numerical work has been done on classifying different patterns of multistability and their basins of attraction (see for example \[Wuensche, 2002\]). Here, I will concentrate on the dynamical mechanisms leading to the existence of multiple attractors. Although some of these results are only rigorously proven for low dimensional systems, their relevance for high dimensional systems is to be expected.

(i) In the particular case of networks with symmetric connections the attracting critical points of the potential function $`V^{(S)}`$ are stable asymptotic states of the dynamics.

(ii) Homoclinic tangencies (the Newhouse phenomenon)

Contrary to earlier conjectures that generic systems might only have finitely many attractors, Newhouse  proved that a class of diffeomorphisms in a two-dimensional manifold has infinitely many attracting periodic orbits (sinks), a result that was later extended to higher dimensions \[Palis & Viana, 1994\]. For two-dimensional manifolds the result is:

Theorem (Newhouse, Robinson ) Let $`f_\mu `$ be a $`C^3`$ map in a 2-dimensional manifold with $`C^1`$ dependence on $`\mu `$ and $`\left|det\left(T_af_{\mu _0}^n\right)\right|<1`$ and let the non-degenerate homoclinic tangency be crossed at non-zero speed at $`\mu =\mu _0`$. Then for $`\epsilon >0`$, $`(\mu _1,\mu _2)(\mu _0,\mu _0+\epsilon )`$ and a residual subset $`J(\mu _1,\mu _2)`$ such that for $`\mu J`$ , $`f_\mu `$ has infinitely many sinks.

Models of such diffeomorphisms were constructed by Gambaudo and Tresser  and Wang  proved that the Newhouse set has positive Hausdorff dimension. After these results, intense research followed on the unfolding of homoclinic tangencies and an essential question was whether, in addition to infinitely many sinks, there would also be infinitely many strange attractors near the homoclinic tangencies. The question was positively answered by Colli . The main result is:

Theorem (Colli) Let $`f_0Diff^{\mathrm{}}(M)`$ be such that $`f_0`$ has a homoclinic tangency between the stable and unstable manifolds of a dissipative hyperbolic saddle $`p_0`$. Then, there is an open set $`\mathrm{\Omega }Diff^{\mathrm{}}(M)`$ such that

(a) $`f_0\overline{\mathrm{\Omega }}`$

(b) there is a dense subset $`D\mathrm{\Omega }`$ such that for all $`fD`$ , $`f`$ exhibits infinitely many coexisting Hénon-like strange attractors.

Having established the existence of infinitely many sinks and infinitely many strange attractors near homoclinic tangencies, a question of practical importance is the stability of the phenomenon under small random perturbations of the deterministic dynamics. It turns out that the answer to this question is negative. Therefore under small random perturbations only finitely many physical measures will remain.

Theorem (Araújo ) Let $`f:MM`$ be a diffeomorphism of class $`C^r,r>1`$, of a compact connected boundaryless manifold $`M`$ of finite dimension. If $`f=f_a`$ is a member of a parametric family under parametric noise of level $`\epsilon >0`$, that satisfies the hypothesis:

There are $`KN`$ and $`\xi _0>0`$ such that, for all $`kK`$ and $`xM`$

(A) $`f^k(x,\mathrm{\Delta })B^k(x),\xi _0)`$ ;

(B) $`f^k(x,\nu ^{\mathrm{}})<<m`$ ;

then there is a finite number of probability measures $`\mu _1,\mathrm{}\mu _l`$ in $`M`$ with the properties

1. $`\mu _1,\mathrm{}\mu _l`$ are physical absolutely continuous probability measures;

2. supp$`\mu _i`$supp$`\mu _j`$ for all $`1i<jl`$ ;

3. for all $`xM`$ there are open sets $`V_1=V_1(x),\mathrm{},V_l=V_l(x)\mathrm{\Delta }`$ such that

(a) $`V_iV_j=\mathrm{}`$, $`1i<jl`$ ;

(b) $`\nu ^{\mathrm{}}\left(\mathrm{\Delta }\backslash \left(V_1\mathrm{}V_l\right)\right)=0`$ ;

(c) for all $`1il`$ and $`\nu ^{\mathrm{}}`$ a.e. $`tV_i`$ we have

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\underset{j=0}{\overset{n1}{}}\varphi \left(f^j(x,t)\right)=\varphi 𝑑\mu $$

for every $`\varphi C(M,R)`$. Moreover the sets $`V_1(x),\mathrm{},V_l(x)`$ depend continuously on $`xM`$ with respect to the distance $`d_\nu (A,B)=\nu ^{\mathrm{}}\left(A\mathrm{}B\right)`$ between $`\nu ^{\mathrm{}}`$mod $`0`$ subsets of $`\mathrm{\Delta }`$.

(iii) Small dissipative perturbations of conservative systems

Conservative systems have a large number of coexisting invariant sets, namely periodic orbits, invariant tori and cantori. By adding a small amount of dissipation to a conservative system one finds that some of the invariant sets become attractors. Of course, not all invariant sets of the conservative system will survive when the dissipation is added. However, for sufficiently small dissipation many attractors (mainly periodic orbits) have been observed in typical systems. Poon & Grebogi , Feudel et al.  and Feudel & Grebogi  have extensively studied these effects in the single and double rotor, the Hénon map and the optical cavity map. They find a large number of attractors for a small amount of dissipation, in particular in the double rotor map. The large number of coexisting stable periodic orbits has a complex interwoven basin of attraction structure, with the basin boundaries permeating most of the state space. The chaotic component of the dynamics is in the chaotic saddles embedded in the basin boundary. The systems are also found to be highly sensitive to small amounts of noise. The problem of migration between attractors and their stability in multiple-attractor systems has also been studied by other authors \[Weigel & Atlee Jackson, 1998\] \[Kaneko, 1997\] \[Dutta et al., 1999\].

Rigorous results may be obtained in particular cases using the ideas of deformation stability \[Duarte & Vilela Mendes, 1983\] \[Vilela Mendes & Duarte, 1982\] \[Vilela Mendes, 1988\] \[Lima & Vilela Mendes, 1988\]. For example, let an $`\epsilon `$family of maps be

$$\begin{array}{ccc}x^{^{}}\hfill & =\hfill & bx+y+f(x,y,\epsilon )\hfill \\ y^{^{}}\hfill & =\hfill & y+g(x,y,\epsilon )\hfill \end{array}$$
(15)
with $`f(x,y,0)=g(x,y,0)=0`$ For $`\epsilon =0`$ the map has marginally stable periodic orbits of all periods. Under perturbation some of the orbits become stable ones.

Theorem \[Vilela Mendes, 2000a\] If $`f`$ and $`g`$ are jointly $`C^2`$ in $`(x,y,\epsilon )`$ with $`f(x,y,0)=g(x,y,0)=0`$, there is an $`\overline{\epsilon }`$ such that for $`\left|\epsilon \right|<\left|\overline{\epsilon }\right|`$ an interior orbit of period $`p`$ of the unperturbed map becomes a stable orbit of the perturbed map if and only if:

(1) $`\underset{n=0}{\overset{p1}{}}_\epsilon g(x_n^{(0)},y^{(0)},0)_{\epsilon =0}=0`$

(2) $`\epsilon _\epsilon \underset{n=0}{\overset{p1}{}}\left\{_xg(x_n^{(0)},y^{(0)},\epsilon )+(1b)_yg(x_n^{(0)},y^{(0)},\epsilon )\right\}_{\epsilon =0}<0`$

(iv) Coupled dynamical systems near period-doubling accumulation points

An example is a system of two coupled quadratic maps \[Carvalho et al., 2001\],

$$\begin{array}{c}x_1(t+1)=1\mu _{}\left((1\epsilon )x_1(t)+\epsilon x_2(t)\right)^2\\ x_2(t+1)=1\mu _{}\left(\epsilon x_1(t)+(1\epsilon )x_2(t)\right)^2\end{array}$$
(16)
with $`x[1,1]`$, and $`\mu _{}=1.401155\mathrm{}`$ , which is the parameter value of the period doubling accumulation point. The result is that for any $`N`$ there is a sufficiently small $`\epsilon \left(N\right)`$ such that there are $`N`$ distinct stable periodic orbits.

### 2.4 Network examples

#### 2.4.1 An excitatory-inhibitory network

Excitatory-inhibitory networks exhibit a rich variety of activity patterns. They have been identified as underlying several biological processes like image segmentation, sleep rhythms, control of movement in the basal ganglia, etc. \[Terman, 2002\]. Here one considers a simple network with two populations, one composed of excitatory and the other of inhibitory cells. The equations of motion are

$$\begin{array}{ccccc}\underset{i}{\overset{}{x}}\hfill & =\hfill & \alpha _{j=N+1}^{2N}f\left(x_j\right)+\beta g\left(t\right)\hfill & & i=1,\mathrm{},N\hfill \\ \underset{i}{\overset{}{x}}\hfill & =\hfill & \alpha _{j=1}^Nf\left(x_j\right)\gamma x_i\hfill & & i=N+1,\mathrm{},2N\hfill \end{array}$$
(17)
The first population ($`i=1,\mathrm{},N`$) receives a time-dependent driving external signal $`g\left(t\right)`$ and inhibitory inputs from the second population. The second population receives excitatory inputs from the first population and has a decay rate $`\gamma `$. The activation function is

$$f\left(x\right)=1\mathrm{exp}\left(\mu x\right)$$
(18)
In terms of the global functions described before, the dynamics has a simple form

$$\underset{i}{\overset{}{x}}=\frac{V}{x_i}+\underset{j}{}\mathrm{\Gamma }_{ij}\frac{H}{x_j}$$
(19)
with

$$\begin{array}{ccc}V\hfill & =\hfill & \beta g\left(t\right)_ix_i+\frac{\gamma }{2}_ix_i^2\hfill \\ H\hfill & =\hfill & \alpha _i^{x_i}f\left(\xi \right)𝑑\xi \hfill \end{array}$$
(20)
$`\mathrm{\Gamma }`$ being the matrix111the degeneracy of the symplectic form is lifted by a particular choice of coordinates

$$\mathrm{\Gamma }=\left(\begin{array}{cccccccc}0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill & 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill \\ 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill & 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill & 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill \\ 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill \\ 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill \end{array}\right)$$
(21)
Except for the contribution of the driving signal and the decay constant, the dynamics is Hamiltonian, the symmetric connection coefficients $`w_{ij}^{(S)}`$ in (7) vanishing identically. Therefore, the reaction of the network to an external signal is simply a damped oscillation. This is illustrated in Fig.1 where, starting from a random initial condition and after a relaxation period, a square wave $`g\left(t\right)`$ with zero baseline is applied to the network. The base level of each unit depends on the initial conditions.

Notice however that, if the agents in the network are biological cells, their activation cannot be negative. Then, the dynamical system should be modified from $`\underset{i}{\overset{}{x}}=F_i\left(x\right)`$ (where $`F_i\left(x\right)`$ is the right hand side of (17)) to

$$\underset{i}{\overset{}{x}}=F_i\left(x\right)OR(sign\left(x_i\right),sign\left(F_i\right))$$
(22)
OR being the logical function

|  | $`+`$ | $``$ |
| --- | --- | --- |
| $`+`$ | 1 | 1 |
| $``$ | 1 | 0 |

(23)
In this case, because the right hand side is no longer of Cohen-Grossberg form, the decomposition (19)$``$(20) gives only qualitative information on the dynamical behavior of the system. This is illustrated in Fig.2 where the same square-wave driving force, as in Fig.1, is applied to the modified network.

#### 2.4.2 A gene regulation network

The p53 gene was one of the first tumour-suppressor genes to be identified, its protein acting as an inhibitor of uncontrolled cell growth. The p53 protein has been found not to be acting properly in most human cancers, due either to mutations in the gene or inactivation by viral proteins or inhibiting interactions with other cell products. Although apparently not required for normal growth and development, p53 is critical in the prevention of tumour development, contributing to DNA repair, inhibiting angiogenesis and growth of abnormal or stressed cells \[May & May, 1999\] \[Vogelstein et al., 2000\] \[Woods & Vousden, 2001\] \[Taylor et al., 1999\] \[Vousden, 2000\]. In addition to its beneficial anticancer activities it may also have some detrimental effects in human aging \[Sharpless & DePinho, 2002\].

The p53 gene does not act by itself, but through a very complex network of interactions \[Kohn, 1999\]. Here I will discuss a simplified network, which although not being accurate in biological detail, tends to capture the essential features of the p53 network as it is known today. In particular, several different products and biological mechanisms are lumped together into a single node when their function is identical. The network is depicted in Fig.3. The arrows and signs denote the excitatory or inhibitory action of each node on the others and the letters $`b,g,c,r,p,m,a`$ denote their intensities (or concentrations).

The p53 protein is assumed to be produced at a fixed rate ($`k_p`$) and to be degraded after ubiquitin labelling. The MDM2 protein being one of the main enzymes involved in ubiquitin labelling, the inhibitory node ($`m`$) is denoted MDM2. There is a positive feedback loop from p53 to MDM2, because the p53 protein, binding to the regulatory region of the MDM2 gene, stimulates the transcription of this gene into mRNA.

Under normal circumstances the network is “off” or operates at a low level. The main activation pathways are the detection of cell anomalies ($`a`$), like DNA damage, or aberrant growth signals, such as those resulting from the expression of several oncogenes (the p14ARF pathway, $`r`$) . They inhibit the degradation of the p53 protein, which may then reach a high level. There are several distinct activation pathways. For example, in some cases phosphorylation of the p53 protein blocks its interaction with MDM2 and in others it is a protein that binds to MDM2 and inhibits its action. However, the end result being a decrease in the MDM2 efficiency, they may both be described as inhibitory inputs to the MDM2 node.

The p53 protein controls cell growth and proliferation, either by blocking the cell division cycle, or activating apoptosis or inhibiting the blood-vessel formation ($`b`$) that is stimulated by several tumors. In our simplified p53 network all these effects are coded on the following set of equations

$$\begin{array}{ccc}\frac{dp}{dt}\hfill & =\hfill & k_pW_{pm}f_m\left(m\right)\hfill \\ \frac{dm}{dt}\hfill & =\hfill & W_{mp}f_p\left(p\right)W_{mr}f_r\left(r\right)W_{ma}a\gamma _m^{}m\hfill \\ \frac{db}{dt}\hfill & =\hfill & W_{bc}f_c\left(c\right)W_{bp}f_p\left(p\right)\hfill \\ \frac{dc}{dt}\hfill & =\hfill & W_{cg}g+W_{cb}f_b\left(b\right)W_{cp}f_p\left(p\right)\hfill \\ \frac{dr}{dt}\hfill & =\hfill & W_{rg}g\gamma _rr\hfill \end{array}$$
(24)
One should note that an increased level of cellular p53 is not by itself sufficient for it to become a transcriptional activator controlling cell growth. Conformational changes of the protein are also needed which are stimulated by the activation pathways or may be therapeutically induced. Also some viruses produce proteins that inactivate p53. All this means that in reality some of the coupling constants in Eqs.(24), (for example $`W_{cp}`$) may also be dynamical variables.

The regulation functions $`f()`$ are positive nonlinear functions with a threshold and a saturation level. By shifting variables to compensate for thresholds and rescaling the coupling constants they may be normalized by the coefficient of the linear part, that is

$$f_i\left(x_i\right)=x_i+\mathrm{}$$
(25)
With a rescaling of $`p,m,b,c,r`$ and redefinition of the constants we may consider

$$k_p=W_{mp}=W_{bc}=W_{cg}=W_{rg}=1$$
(26)
Furthermore, from the last equation in (24)

$$r\left(t\right)=\frac{1}{\gamma _r}g+\left(r\left(0\right)\frac{g}{\gamma _r}\right)e^{\gamma _rr}$$
(27)
Replacing $`r`$ by its steady state value $`g/\gamma _r`$, and rescaling $`W_{mr}`$ we are left with

$$\begin{array}{ccc}\frac{dp}{dt}\hfill & =\hfill & 1W_{pm}f_m\left(m\right)\hfill \\ \frac{dm}{dt}\hfill & =\hfill & f_p\left(p\right)W_{mr}gW_{ma}a\gamma _mm\hfill \\ \frac{db}{dt}\hfill & =\hfill & f_c\left(c\right)W_{bp}f_p\left(p\right)\hfill \\ \frac{dc}{dt}\hfill & =\hfill & g+W_{cb}f_b\left(b\right)W_{cp}f_p\left(p\right)\hfill \end{array}$$
(28)
a system of four dynamical variables and two control parameters $`g`$ and $`a`$.

Using the dynamics decomposition discussed in Section 2 one obtains

$$\begin{array}{ccc}V^{(S)}\hfill & =\hfill & ^pf_p^{^{}}\left(\xi \right)𝑑\xi ^cf_c^{^{}}\left(\xi \right)𝑑\xi +^m\left(W_{mr}g+W_{ma}a+\gamma _m\xi \right)f_m^{^{}}\left(\xi \right)𝑑\xi \hfill \\ & & +\frac{1}{2}_{x_i=p,m,b,c}W_{ij}^{(S)}f_i\left(x_i\right)f_j\left(x_j\right)\hfill \\ H\hfill & =\hfill & _{x_i=p,m,b,c}^{x_i}f_i\left(\xi \right)𝑑\xi \hfill \end{array}$$
(29)
$$g_{ij}=\frac{1}{f_i^{^{}}\left(x_i\right)}\delta _{ij}\mathrm{\Gamma }_{ij}=W_{ij}^{(A)}$$
(30)
| $`(i,j)`$ | $`W_{ij}^{(S)}`$ | $`W_{ij}^{(A)}`$ |
| --- | --- | --- |
| $`pm`$ | $`\frac{1}{2}\left(W_{pm}1\right)`$ | $`\frac{1}{2}\left(W_{pm}+1\right)`$ |
| $`pb`$ | $`\frac{1}{2}W_{bp}`$ | $`\frac{1}{2}W_{bp}`$ |
| $`pc`$ | $`\frac{1}{2}W_{cp}`$ | $`\frac{1}{2}W_{cp}`$ |
| $`mb`$ | $`0`$ | $`0`$ |
| $`mc`$ | $`0`$ | $`0`$ |
| $`bc`$ | $`\frac{1}{2}\left(W_{cb}+1\right)`$ | $`\frac{1}{2}\left(W_{cb}1\right)`$ |

(31)
The reaction of $`m`$ and $`c`$ to external stimuli ($`g`$ and $`a`$) and the production rate of $`p`$ are coded on the first three terms of the potential function $`V^{(S)}`$. For coupling constants of order unit, one sees from (31) the existence of a damped Hamiltonian oscillation for the $`pm`$ system, and a dangerous runaway behavior of $`bc`$ arising from its dominantly gradient dynamics. The action of $`p`$ on $`b`$ and $`c`$ is of mixed gradient-Hamiltonian type. Hence, from inspection of the nature of the global functions describing the dynamics, one concludes that (at least in this model) the controlling action of p53 may only be effective in particular circumstances. That is, it will depend on the initial conditions. This conclusion is now checked by a detailed study of the solutions.

Consider first the linear approximation to the system. The solutions are, for the $`pm`$ system

$$\begin{array}{ccc}p\left(t\right)\hfill & =\hfill & \stackrel{}{p}+p^{}\left(t\right)\hfill \\ m\left(t\right)\hfill & =\hfill & \stackrel{}{m}+m^{}\left(t\right)\hfill \end{array}$$
(32)
with

$$\begin{array}{ccc}\stackrel{}{p}\hfill & =\hfill & W_{mr}g+W_{ma}a+\frac{\gamma _m}{W_{pm}}\hfill \\ \stackrel{}{m}\hfill & =\hfill & \frac{1}{W_{pm}}\hfill \end{array}$$
(33)
$$\begin{array}{cc}p^{}\left(t\right)=\hfill & e^{\gamma _m/2}\{(p\left(0\right)\stackrel{}{p})\mathrm{cos}\alpha t+\frac{1}{\alpha }(\frac{\gamma _m}{2}(p\left(0\right)\stackrel{}{p})W_{pm}(m\left(0\right)\frac{1}{W_{pm}}))\mathrm{sin}\alpha t\}\hfill \\ m^{}\left(t\right)=\hfill & e^{\gamma _m/2}\left\{\left(m\left(0\right)\frac{1}{W_{pm}}\right)\mathrm{cos}\alpha t+\frac{1}{\alpha }\left(p\left(0\right)\stackrel{}{p}\frac{\gamma _m}{2}\left(m\left(0\right)\frac{1}{W_{pm}}\right)\right)\mathrm{sin}\alpha t\right\}\hfill \end{array}$$
(34)
and $`\alpha =\sqrt{W_{pm}\gamma ^2/4}`$.

As expected, one sees a damped oscillatory behavior of the $`pm`$ system and, in the absence of stimuli ($`a=g=0`$) the $`p`$ level is small and controlled by the degradation of $`m`$.

For $`b`$ and $`c`$ one now obtains

$$\left(\begin{array}{c}b\left(t\right)\hfill \\ c\left(t\right)\hfill \end{array}\right)=\left(\begin{array}{c}W_{bp}\stackrel{}{p}\hfill \\ \frac{\left(W_{cp}\stackrel{}{p}g\right)}{W_{cb}}\hfill \end{array}\right)_0^te^{A\left(t\tau \right)}\left(\begin{array}{c}W_{bp}p^{}\left(\tau \right)\hfill \\ W_{cp}p^{}\left(\tau \right)\hfill \end{array}\right)𝑑\tau +e^{At}\left(\begin{array}{c}c\left(0\right)W_{bp}\stackrel{}{p}\hfill \\ b\left(0\right)\frac{W_{cp}\stackrel{}{p}g}{W_{cb}}\hfill \end{array}\right)$$
(35)
where $`A`$ is the matrix

$$A=\left(\begin{array}{cc}0\hfill & 1\hfill \\ W_{cb}\hfill & 0\hfill \end{array}\right)$$
(36)
This matrix has eigenvalues $`\pm \sqrt{W_{cb}}`$ implying that $`b\left(t\right)`$ and $`c\left(t\right)`$ are going to have terms proportional to $`\mathrm{exp}\left(t\sqrt{W_{cb}}\right)`$ and $`\mathrm{exp}\left(t\sqrt{W_{cb}}\right)`$. Hence $`p`$ (p53) will only have a controlling effect on cell proliferation if the coefficient of the exponentially growing terms becomes negative. Multiplying (35) on the left by the matrix $`\frac{1}{2}\left(\begin{array}{cc}\sqrt{W_{cb}}\hfill & 1\hfill \\ \sqrt{W_{cb}}\hfill & 0\hfill \end{array}\right)`$ that diagonalizes $`A`$ one obtains the coefficient of the exponentially growing term

$$\begin{array}{ccc}B\left(t\right)\hfill & =\hfill & e^{t\sqrt{W_{cb}}}\{\frac{\sqrt{W_{cb}}}{2}(c\left(0\right)W_{bp}\stackrel{}{p})+\frac{1}{2}(b\left(0\right)\frac{W_{cp}\stackrel{}{p}g}{W_{cb}})\}\hfill \\ & & _0^te^{\left(t\tau \right)\sqrt{W_{cb}}}\left(\frac{\sqrt{W_{cb}}}{2}W_{bp}+\frac{1}{2}W_{cp}\right)p^{}\left(\tau \right)𝑑\tau \hfill \end{array}$$
(37)
The conclusion is that control of cell proliferation is obtained only if $`t`$ such that $`B\left(t\right)<0`$. Therefore it depends strongly on the initial conditions. This conclusion, inferred both from the dynamical decomposition and the linear approximation is borne out by simulation of the nonlinear problem. Figs.4 and 5 show two time evolutions of the equations (24) with $`f\left(x\right)=\mathrm{tanh}\left(x\right)`$ , $`W_{pm}=W_{mr}=W_{ma}=W_{bp}=W_{cb}=W_{cp}=1`$, $`\gamma _m=0.01`$, $`g=a=1`$ and the vector field $`\underset{i}{\overset{}{x}}=F_i`$ truncated to $`\underset{i}{\overset{}{x}}=F_i`$OR$`(sign\left(x_i\right),sign\left(F_i\right))`$, because concentrations cannot become negative. The behavior depends strongly on the value of the initial conditions. In conclusion, the implication (of the model) is that unless p53 starts acting soon enough its action is useless and other means have to be used to control cell proliferation.

### 2.5 Evolving networks

In many networks found in Nature, as important as the structure of the network, is the path that the network took to reach that final state. Social or economic networks, industrial, transportation and communication networks, ecological webs, biological networks, all are examples of evolving networks. In many cases their complex structure is a simple consequence of the principles of their growth. Several network growth schemes have been studied (see Albert & Barabási , Dorogovtsev & Mendes  and Pastor-Satorras et al.  for reviews).

Network evolution occurs either by the addition or elimination of interactions between existing nodes or by the addition of new nodes. In both cases, network evolution may be looked at as a dynamical system in the space of network connections. In the case of growing networks, this dynamical point of view may also be used by considering the evolution from zero of previously vanishing connections.

This dynamical approach will be explored here. Using the global function description, discussed in Section 2.1, two types of evolving networks will be considered. The simplest situation occurs when the dynamics of the connections is derived from a potential. In this case, exact expressions for mean values and invariant measures may be obtained.

Consider

$$V_1\left(\left\{W\right\}\right)=\alpha \underset{i<j}{}W_{ij}^2\left(W_{ij}1\right)^2+\beta \underset{ijl}{}\left(W_{ij}1\right)^2W_{jl}^2$$
(38)
with the network evolving according to

$$\frac{dW_{ij}}{dt}=\frac{V_1}{W_{ij}}$$
(39)
When $`\alpha 0`$ and $`\beta =0`$, the connections evolve either to zero or to one, depending on the initial conditions. Therefore the network (with $`N`$ nodes), as a dynamical system, is a multistable system with $`2^{N\left(N1\right)/2}`$ different equilibrium points. A typical configuration, obtained from random initial conditions, is shown in Fig.6 ($`N=100`$) to which corresponds the degree distribution shown in Fig.7.

When $`\beta 0`$ the behavior is quite different, as shown in the typical configuration of Fig.8 and degree distribution Fig.9. The degree $`K_i`$ of a node $`i`$ is defined to be

$$K_i=\underset{j}{}W_{ij}$$
(40)
holding for all intermediate values of $`W_{ij}`$.

One sees that for $`\beta 0`$ some nodes are more connected than others. $`V_1\left(\left\{W\right\}\right)`$ with $`\beta 0`$ is a model for preferential attachment.

It is not be practical to obtain mean values and distributions directly from simulations. This being a multistable system many different simulations with well distributed initial conditions would be required to obtain accurate values. However, in this case, exact expressions may be obtained from the unique invariant measure for the system with small random perturbations, as discussed in Section 3.3

$$\rho ^\epsilon \mathrm{exp}\left(2\epsilon ^2V_1\left(\left\{W\right\}\right)\right)$$
(41)
As a second example consider

$$V_2\left(\left\{W\right\}\right)=\alpha \underset{i<j}{}W_{ij}^2\left(W_{ij}1\right)^2+\beta \underset{i<j}{}\underset{ki,j}{}\frac{1}{\left|ij\right|}\left(W_{ik}^2+W_{jk}^2\right)\left(\left(W_{ik}1\right)^2+\left(W_{jk}1\right)^2\right)$$
(42)
For $`\beta 0`$ a typical configuration is shown in Fig.10. The main feature is the correlation between node connections. For $`\alpha =1`$, $`\beta =0.05`$ and $`N=100`$ the sum of correlations between the node connections is around $`20`$ whereas for $`\alpha =1`$, $`\beta =0`$ it is $`4.5`$. In conclusion, $`V_2\left(\left\{W\right\}\right)`$ is a model for (approximate) node replication.

## 3 Ergodic tools

Topological and differential notions provide useful characterizations of the overall structure of phase space. However, what is more important for the applications is the dynamics in the phase space regions most frequently visited by the system. This is provided by the ergodic theory, in particular by the classification of invariant measures and their characterization by ergodic parameters.

Let a dynamical system evolve on the support of a measure $`\mu `$ which is left invariant by the dynamics. An ergodic parameter $`I_F\left(\mu \right)`$, characterizing the measure, is obtained whenever the following limit

$$I_F\left(\mu \right)=\underset{T\mathrm{}}{lim}\frac{1}{T}\underset{n=1}{\overset{T}{}}ϝ\left(f^nx_0\right)$$
(43)
exists for $`\mu `$almost every $`x_0`$. For continuous-time dynamics $`f`$ denotes the time-one map.

### 3.1 Lyapunov and conditional exponents

Lyapunov exponents are the most widely used ergodic parameters. More recently conditional exponents have also been proposed as an useful characterization of the dynamics.

Let $`f:MM`$ , with $`MR^m`$, $`\mu `$ a measure invariant under $`f`$ and $`\mathrm{\Sigma }`$ a splitting of $`M`$ induced by $`R^k\times R^{mk}`$. The conditional exponents are the eigenvalues $`\xi _i^{(k)}`$ and $`\xi _i^{(mk)}`$ of the limits

$`\underset{n\mathrm{}}{lim}\left(D_kf^n(x)D_kf^n(x)\right)^{\frac{1}{2n}}`$ (44)
$`\underset{n\mathrm{}}{lim}\left(D_{mk}f^n(x)D_kf^n(x)\right)^{\frac{1}{2n}}`$
where $`D_kf^n`$ and $`D_{mk}f^n`$ are the $`k\times k`$ and $`mk\times mk`$ diagonal blocks of the full Jacobian. For $`k=m`$ , $`\xi _i^{(m)}=\lambda _i`$ are the Lyapunov exponents.

Proposed by Pecora & Carroll  to characterize synchronization in chaotic systems, rigorous conditions for the existence of these limits have been proven in \[Vilela Mendes, 1998\]. Existence $`\mu `$-almost everywhere of both Lyapunov and conditional exponents is guaranteed by the conditions of Oseledec’s multiplicative ergodic theorem, in particular the integrability condition,

$$\mu (dx)\mathrm{log}^+T(x)<\mathrm{}$$
(45)
$`T`$ being either the Jacobian or its $`k\times k`$ and $`mk\times mk`$ diagonal blocks. The set of points where the limit is defined has full measure and

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}D_kf^n(x)u=\xi _i^{(k)}$$
(46)
with $`0uE_x^i/E_x^{i+1}`$ , $`E_x^i`$ being the subspace of $`R^k`$ spanned by eigenstates corresponding to eigenvalues $`\mathrm{exp}(\xi _i^{(k)})`$.

Based on the spectra of Lyapunov and conditional exponents, several global quantities have been defined to characterize self-organization and creation of structures in networks of multiagent systems with arbitrary connection structures. I list here the definitions and refer to \[Vilela Mendes, 2000b, 2001\] for proofs and examples.

#### 3.1.1 Structure index related to the Lyapunov spectrum

A structure (in a collective system) is a phenomenon with a characteristic scale very different from the scale of the elementary units in the system. In a multi-agent system, a structure in space is a feature at a length scale larger than the characteristic size of the agents and a structure in time is a phenomenon with a time scale larger than the cycle time of the individual agent dynamics. A (temporal) structure index may then be defined by

$$S=\frac{1}{N}\underset{i=1}{\overset{N_s}{}}\frac{T_iT}{T}$$
(47)
where $`N`$ is the total number of components (agents) in the coupled system, $`N_s`$ is the number of structures, $`T_i`$ is the characteristic time of the structure $`i`$ and $`T`$ is the cycle time of the isolated agents (or, alternatively the characteristic time of the fastest structure). A similar definition applies for a spatial structure index, by replacing characteristic times by characteristic lengths.

Structures are collective motions of the system. Therefore their characteristic times are the characteristic times of the separation dynamics, that is, the inverse of the positive Lyapunov exponents. Hence, for the temporal structure index, one may write

$$S=\frac{1}{N}\underset{i=1}{\overset{N_+}{}}\left(\frac{\lambda _0}{\lambda _i}1\right)$$
(48)
the sum being over the positive Lyapunov exponents $`\lambda _i`$. $`\lambda _0`$ is the largest Lyapunov exponent of an isolated component or some other reference value.

The temporal structure index diverges whenever a Lyapunov exponent approaches zero from above. Therefore the structure index diverges at the points where long time correlations develop. Also, when in a multiagent network the coupling between the agents increases, the positive part of the Lyapunov spectrum contracts leading to an effective dimension reduction and to partial synchronization effects \[Vilela Mendes, 1999\].

#### 3.1.2 Exponent entropies and dynamical selforganization

Self-organization in a system concerns the dynamical relation of the whole to its parts. The conditional Lyapunov exponents, being quantities that separate the intrinsic dynamics of each component from the influence of the other parts in the system, provide a measure of dynamical selforganization $`I_\mathrm{\Sigma }(\mu )`$

$$I_\mathrm{\Sigma }(\mu )=\underset{k=1}{\overset{N}{}}\left\{h_k(\mu )+h_{mk}(\mu )h(\mu )\right\}$$
(49)
the sum being over all relevant partitions $`\mathrm{\Sigma }_k=R^k\times R^{mk}`$ and

$$h_k(\mu )=\underset{\xi _i^{(k)}>0}{}\xi _i^{(k)};h_{mk}(\mu )=\underset{\xi _i^{(mk)}>0}{}\xi _i^{(mk)};h(\mu )=\underset{\lambda _i>0}{}\lambda _i$$

are the exponent entropies, that is, the sums over positive conditional and Lyapunov exponents.

$`I_\mathrm{\Sigma }(\mu )`$ may also be given the following dynamical interpretation: Lyapunov exponents measure the rate of information production or, equivalently, they define the dynamical freedom of the system, in the sense that they control the amount of change that is needed today to have an effect on the future. In this sense the larger a Lyapunov exponent is, the freer the system is in that particular direction, because a very small change in the present state will induce a large change in the future. The conditional exponents have a similar interpretation concerning the dynamics as seen from the point of view of each agent and his neighborhood \[Vilela Mendes, 2000b\]. However the actual information production rate is given by the sum of the positive Lyapunov exponents, not by the sum of the conditional exponents. Therefore, $`I_\mathrm{\Sigma }(\mu )`$ is a measure of apparent dynamical freedom (or apparent rate of information production).

Being constructed as functions of well defined ergodic limits, both $`I_\mathrm{\Sigma }(\mu )`$ and $`S`$ are also well defined ergodic parameters. They characterize the dynamics of multiagent networks and, in addition, also provide some insight on the relation between dynamics and the topology of the network \[Araújo et al., 2003\].

### 3.2 The problem of pattern discovery: Computational mechanics

The ultimate practical goal in the study of dynamical systems is the construction of models by which these systems might be predicted and (or) controlled to some useful purpose. In extended systems with many degrees of freedom, unless exact solutions are known, even a knowledge of the (microscopic) equations of motion might not be very useful to predict the collective features and patterns that the system generates. Crutchfield and collaborators \[Crutchfield & Young, 1989\] \[Crutchfield, 1994\] \[Shalizi & Crutchfield, 2001\] have developed a program of pattern discovery and construction of minimal models inferred directly from the data generated by the dynamical systems. Central to this approach is the notion of causal state. Given the knowledge of the infinite past of a system, a causal state is an equivalence class of pasts that have the same conditional distribution of futures. Denoting by $`\stackrel{}{s}`$ and $`\stackrel{}{s}`$ the semi-infinite past and future time sequences of coded states of the system, two past sequences $`\stackrel{}{s_1}`$ and $`\stackrel{}{s_2}`$ belong to the same causal state if

$$P\left(\stackrel{}{s}|\stackrel{}{s_2}\right)=P\left(\stackrel{}{s}|\stackrel{}{s_2}\right)$$

for all $`\stackrel{}{s}`$(except perhaps in a zero measure set). The dynamics of the system is then characterized by the set of causal states and the transition probabilities between them. That is, the system is mapped into a non-deterministic automaton called an $`\epsilon `$machine. Minimality and uniqueness of the $`\epsilon `$machines has been proved. Although more general, this approach bears some relation to the reconstruction of hidden Markov processes and to grammatical inference.

As a tool for network dynamics this approach might be useful whenever analytical equations are intractable or unknown. Reconstruction algorithms for $`\epsilon `$machines were developed in some cases \[Hanson & Crutchfield, 1997\] \[Crutchfield & Feldman, 1997\]. For extended system with many degrees of freedom and irregular connections, one problem might be the large size of the causal state alphabet. Nevertheless this is a very interesting general approach that might be useful to map network dynamics onto probabilistically equivalent automata.

### 3.3 Construction of invariant measures

In general a deterministic system has a multitude of invariant measures. However, some of them have little practical interest, because they are not stable for small random perturbations. Because systems in Nature are subjected to perturbations, only the stable measures are physical measures. In some cases it is possible to use the properties of the deterministic system to identify the physical measures. For example, in Axiom A systems a unique physical measure may be identified with the Sinai-Bowen-Ruelle (SBR) measure, a measure absolutely continuous along unstable manifolds. However in most cases, for example in the multistable systems so frequent in natural networks, the SBR characterization is useless. Instead, it is better to study the stochastic differential equation that is obtained from (1) by addition of a small noise term

$$dx_i=X_i\left(x\right)dt+\epsilon \sigma \left(X\right)dW_t$$
(50)
$`W_t`$ being a Wiener process and $`\sigma \left(X\right)`$ a $`X`$dependent diffusion coefficient. A great deal of information on the invariant measure for this process may be obtained using the theory of small random perturbations of dynamical systems \[Freidlin & Wentzell, 1984, 1994\] \[Kifer, 1974\].

If, in the decomposition (8), $`X\left(x\right)`$ has only a gradient component, an explicit form for the invariant measure may be obtained. If

$$X\left(x\right)=_{(g)}V\left(x\right)$$
(51)
$`_{(g)}`$ being the gradient in the metric

$$ds^2=a_{ij}\left(x\right)dx_idx_j$$
(52)
with $`\sigma \left(x\right)`$ in (50) chosen such that

$$a_{ij}\left(x\right)=\left(\sigma \left(x\right)\sigma ^{}\left(x\right)\right)_{ij}^1=g_{ij}\left(x\right)$$
(53)
then, the density of the invariant measure is

$$\rho ^\epsilon \left(x\right)=C_\epsilon \mathrm{exp}\left(2\epsilon ^2V\left(x\right)\right)$$
(54)
as may be easily checked from the forward Kolmogorov equation. In this case, finding the stable minima and level sets of $`V\left(x\right)`$ one characterizes the multistability of the network, their basins of attraction and, from the values of $`V\left(x\right)`$ in these sets, the relative occurrence probability of each attractor.

For general $`X\left(x\right)`$, small $`\epsilon `$ estimates of the invariant measure for (50) are also possible. Here the crucial role is played by the functional

$$S_{oT}\left(\phi \right)=\frac{1}{2}_0^T\underset{ij}{}a_{ij}\left(\phi _t\right)(\stackrel{}{\phi _t^i}X^i\left(\phi _t\right))(\stackrel{}{\phi _t^j}X^j\left(\phi _t\right))dt$$
(55)
and the infimum

$$U(x,y)=inf\{S_{0T}\left(\phi \right):\phi _0=x,\phi _T=y,t[0,T]\}$$
(56)
taken over intervals $`[0,T]`$ of arbitrary length.

An equivalence relation is established between points in the domain by $`xy`$ if $`U(x,y)=U(y,x)=0`$. Let the domain be partitioned into a number of compacta $`\left\{K_i\right\}`$ with each $`\omega `$limit set of the deterministic dynamics contained entirely in one compactum and $`xy`$ inside each compactum. Then, the (small $`\epsilon `$) asymptotics of the invariant measure is obtained from the invariant measure of the Markov chain of transitions between the compacta. For sufficiently small $`\epsilon `$ the measure of each compactum is approximated by

$$\mathrm{exp}\left\{\epsilon ^2\left(W\left(K_i\right)\underset{i}{\mathrm{min}}W\left(K_i\right)\right)\right\}$$
(57)
where

$$W\left(K_i\right)=\underset{gG\left(i\right)}{\mathrm{min}}\underset{\left(mn\right)g}{}V(K_m,K_n)$$
(58)
$`V(K_m,K_n)`$ is the minimum of the function (56) between points in compacta $`K_m`$ and $`K_n`$ and the sum runs over graphs that have exactly one closed cycle and this cycle contains the compactum $`K_i`$. For proofs I refer to \[Freidlin & Wentzell, 1984\].

### 3.4 A family of ergodic parameters

Ergodic parameters like the Lyapunov and the conditional exponents, are global functions of the invariant measure. However, the invariant measure itself contains more information. Ergodic parameters being defined by infinite-time limits, these quantities will fluctuate and, in general, fluctuations will not be Gaussian. The quantity describing the fluctuations is again an ergodic parameter and the same reasoning applies in turn to its fluctuations, etc. \[Ruelle, 1987\]. Therefore, to characterize the measure, a larger set of parameters is needed. To construct this larger set from the fluctuations is not very practical and a different approach will be followed here, namely a variational approach.

In a restricted sense, a variational principle states that the equations of motion may be written in the form $`\delta S=0`$, where $`S`$ is a functional of the dynamical variables and $`\delta `$ is the Gateaux derivative. Only a limited set of dynamical systems may be described by a variational principle in this restricted sense. However, if one only requires that $`\delta S=0`$ and the equations of motion possess the same set of solutions, essentially all differential equation problems admit a variational formulation \[Tonti, 1992\]. Let

$$\underset{i}{\overset{}{x}}=X_i\left(x\right)$$
(59)
be a differentiable continuous-time dynamical system and $`S`$ be the functional

$$S=_0^Tdtd\tau \underset{i}{}\{\underset{i}{\overset{}{x}}\left(t\right)X_i\left(x\left(t\right)\right)\}g(t,\tau )\{\underset{i}{\overset{}{x}}\left(\tau \right)X_i\left(x\left(\tau \right)\right)\}$$
(60)
where $`g(t,\tau )`$ is a symmetric kernel ($`g(t,\tau )=g(\tau ,t)`$). Let us compute the Gateaux derivative for variations in space restricted by the boundary conditions

$$u\left(0\right)=u\left(T\right)=0$$
(61)
From

$$\delta _uS=_0^T𝑑t𝑑\tau \underset{i,k}{}u_k(t)\left\{\delta _{ki}\frac{dg(t,\tau )}{dt}_kX_i\left(x\left(t\right)\right)g(t,\tau )\right\}\left\{\underset{i}{\overset{}{x}}\left(\tau \right)X_i\left(x\left(\tau \right)\right)\right\}$$
(62)
we have

Lemma: The equations of motion (59) and the critical points of the functional ($`\delta _uS=0`$) have the same set of solutions if

$$K(t,\tau )=\delta _{ki}\frac{dg(t,\tau )}{dt}_kX_i\left(x\left(t\right)\right)g(t,\tau )$$
(63)
is invertible.

Remarks:

a) If $`K(t,\tau )`$ is not invertible, the solutions of the equations of motion are still critical points of the functional, but this one might have other solutions.

b) A variational principle, with only $`u\left(0\right)=0`$ being required, may also be obtained by choosing a kernel such that $`g(t,T)=0`$.

The critical points of the $`S`$ functional contain the same information as the equations of motion. Therefore the dynamics may be characterized by the properties of the critical points, in particular by their Hessian matrix. Computing the second Gateaux derivative on the orbits one obtains

$$\delta _{u,v}^2S|_{\delta S=0}=_0^T𝑑t𝑑\tau \underset{i,j}{}u_i\left(t\right)v_j\left(\tau \right)H_{ij}(t,\tau )$$
(64)
with

$$H_{ij}(t,\tau )=\underset{k}{}\left\{\begin{array}{c}\frac{dg(t,\tau )}{dt}_jX_i\left(x\left(\tau \right)\right)_iX_k\left(x\left(t\right)\right)g(t,\tau )_jX_k\left(x\left(\tau \right)\right)\\ +\frac{d^2g(t,\tau )}{dtd\tau }\delta _{ij}+_iX_j\left(x\left(t\right)\right)\frac{dg(t,\tau )}{d\tau }\end{array}\right\}$$
(65)
Now assume that the symmetric kernel $`g(t,\tau )`$ is a function of finite support of $`t\tau `$

$$g(t,\tau )=g\left(t\tau \right)=0for\left|t\tau \right|>r$$
(66)
Define

$$J_{0,T}^{(n)}=_0^TTr\left(H^n(t,\tau )\right)𝑑t𝑑\tau $$
(67)
as well as

$$J_{0,T}^{(n)}=J_{0,T}^{(n)}+_{Tr}^{T+r}\left|Tr\left(H^n(t,\tau )\right)\right|𝑑t𝑑\tau +_r^r\left|Tr\left(H^n(t,\tau )\right)\right|𝑑t𝑑\tau $$

Then

$$J_{0,T_1+T_2}^{(n)}J_{0,T_1}^{(n)}+J_{T_1,T_2}^{(n)}$$

and we are in the conditions of Kingman’s sub-additive ergodic theorem. Taking limits, if both $`X_i`$ and $`_iX_k`$ are bounded $`J_{0,T}^{(n)}`$and $`J_{0,T}^{(n)}`$ differ only by a finite quantity and one concludes:

Theorem: If $`\mu `$ is an invariant measure of the dynamics in (59), $`X_i`$ and $`_iX_k`$ are bounded and there is $`M0`$ such that $`J_{0,T}^{(n)}M`$ for sufficiently large $`T`$, then the limit

$$I_n\left(\mu \right)=\underset{T\mathrm{}}{lim}\frac{1}{T}J_{0,T}^{(n)}$$
(68)
exists and

$$\underset{T\mathrm{}}{lim}\frac{1}{T}J_{0,T}^{(n)}d\mu =\underset{T\mathrm{}}{lim}\frac{1}{T}J_{0,T}^{(n)}𝑑\mu $$

$`I_n\left(\mu \right)`$ for $`n=1,2,\mathrm{}`$ is a family of ergodic parameters for the $`\mu `$measure preserving dynamics.

A similar construction for discrete-time maps may be found in \[Vilela Mendes, 1984\] \[Carreira et al., 1991\].

### 3.5 Synchronization, mode-locking and dynamical correlations

The onset of correlated motions in coupled many-agent systems is a phenomenon of widespread occurrence in many scientific fields. The most dramatic effect is the synchronization of assemblies of coupled dynamical systems which, when in isolation, may have quite different rhythms \[Pikovski et al., 2001\]. Examples are biological rhythms \[Winfree, 1967\] like the pacemaker cells in the heart \[Peskin, 1975\], neural systems \[Golomb & Hansel, 2000\], synchronous metabolism \[Aldridge & Pye, 1976\], flashing fireflies \[Buck, 1988\], laser arrays \[Jiang & McCall, 1993\], even fads and social trends may be interpreted as synchronization of distinct agent dynamics. The study of the correlated behavior of many-agent dynamics is also closely related to the problem of control in extended dynamical systems.

I will consider both the coupled behavior of non-chaotic systems (oscillators with distinct individual frequencies) and of systems with isolated chaotic dynamics. In both cases one may distinguish between globally coupled systems and systems where each agent has a limited range or number of interacting partners.

For systems of oscillators the canonical example is the Kuramoto model \[Kuramoto, 1984, 1991\],

$$\frac{d\theta _i}{dt}=\omega _i+\frac{K}{N1}\underset{j=1}{\overset{N}{}}\mathrm{sin}\left(\theta _j\theta _i\right)$$
(69)
with $`K>0`$ and the frequencies $`\omega _i`$ randomly distributed around a central value $`\omega _0`$ with the shifted Cauchy distribution

$$p\left(\omega \right)=\frac{\gamma }{\pi \left[\gamma ^2+\left(\omega \omega _0\right)^2\right]}$$
(70)
A great deal of work has been done on this model (for a review see Strogatz ). The existence of a synchronized cluster is characterized by the order parameter

$$r\left(t\right)=\left|\frac{1}{N}\underset{j=1}{\overset{N}{}}e^{i\theta _j\left(t\right)}\right|$$
(71)
It is found that in the $`N\mathrm{}`$ and $`t\mathrm{}`$ limit, $`r=0`$ for $`K<2\gamma `$ and $`r=\sqrt{1\left(2\gamma /K\right)}`$ for $`K2\gamma `$. That is, there is a coupling threshold above which part of the oscillators starts to synchronize. Figs.11 and 12 show the nonsynchronized (at $`K=\gamma `$) and the synchronized (at $`K=5\gamma `$) behavior for $`100`$ oscillators. The upper plot displays the color-coded values of the oscillator variables at the end of each unit time interval. The lower plots show the time evolution of the order parameter. Fig.13 compares the numerically computed Lyapunov spectrum in the synchronized and non-synchronized situations. One sees that even below the synchronization threshold ($`K=2\gamma `$), part of the Lyapunov exponents becomes negative, meaning that there are many contracting directions, implying an effective dimension-reduction of the asymptotic behavior of the system. This clearly suggests that synchronization is not the whole story and that even before synchronization strong correlations must develop between the dynamics of the individual oscillators.

A type of correlation, of which synchronization is a limiting case is mode-locking. Mode-locking is the entrainment of some integer combination of the frequencies to zero. It also plays an important role in the dynamics of coupled oscillators \[MacKay, 1994\]. However even if all the effective frequencies are incommensurable, the existence of negative Lyapunov directions, implies the existence of dynamical correlations between the oscillators. What is important is the dimension of the invariant measure and the correlations may be characterized by the eigenvectors of the Lyapunov spectrum. These notions are better clarified in a simple model with exactly computable Lyapunov spectrum. Let the dynamics of an assembly of discrete-time oscillators be

$$x_i\left(t+1\right)=x_i\left(t\right)+\omega _i+\frac{k}{N1}\underset{j=1}{\overset{N}{}}f_\alpha \left(x_jx_i\right)$$
(72)
with $`x_i[0,1)`$ and $`f_\alpha \left(x_jx_i\right)=\alpha \left(x_jx_i\right)`$ $`\left(mod1\right)`$ and the $`\omega _i^{}s`$ distributed according to $`p\left(\omega \right)`$, as above.

The Lyapunov spectrum is composed of one isolated zero and $`\mathrm{log}\left(1\frac{N}{N1}\alpha k\right)`$ ($`N1`$)-times. However, although (for all $`k>0`$) $`N1`$ contracting directions are always present, it is only for sufficiently large $`k`$ that synchronization effects emerge as shown in Figs.14 and 15. Nevertheless dynamical correlations do exist for all $`k`$, no matter how small and the Lyapunov dimension is always one. In this case, the eigenvectors of the Lyapunov spectrum may be exactly computed and the correlations explicitly identified. This is illustrated in Fig.16.

So far I have dealt with coupled oscillators, that is, with systems which have individual nonchaotic dynamics. Another important field with many practical applications refers to the case where the individual node dynamics is chaotic. Synchronization of chaotic systems has been extensively studied (for a review see Pecora et al. ) and is still a field of current research \[Wei et al., 2002\]. However, as in the oscillators, for networks of chaotic elements the interesting phenomena go beyond synchronization effects. Correlations and self-organization effects may be characterized by ergodic parameters. I refer to \[VilelaMendes, 1999, 2000b\] for illustrative examples of networks of chaotic elements both globally connected and with a limited range of interactions. The Lyapunov spectrum and the entropies associated to the conditional exponents provide a characterization of the emergent phenomena. It should be noticed that dynamical correlations play an important role in the organization of the dynamics even when there is no reduction of the Lyapunov dimension \[Vilela Mendes, 1999\]. As before these correlations are associated to the eigenvectors of the Lyapunov spectrum.

### 3.6 Dynamics and network topology

The topology of the network connections has a determining effects on the dynamical phenomena taking place in the network \[Watts, 1999\]. Local clustering and far-reaching connections in the network influence the spread of infectious diseases \[Boots & Sasaki, 1999\] \[Keeling, 1999\] \[Pastor-Satorras & Vespignani, 2002\] \[Lloyd & May, 2001\] or social fads. The nature and range of the couplings influences the speed of signal propagation and even the computational abilities of the network. On the other hand the topological structure of the networks by influencing the dynamics, may have a feedback effect on network growing and therefore also evolutionary significance.

In particular the small world topology \[Watts & Strogatz, 1998\] \[Watts, 1999\] \[Strogatz, 2001\] (with both small path length and large clustering) has been found to modify or enhance coherent behavior effects \[Lago-Fernandez et al., 2000\] \[Gade & Hu, 2000\] and in general influence the dynamics in the network \[Yang, 2001\] \[Kulkarni et al., 1999\]. An attempt has also been made to relate the ergodic parameters of a dynamical system, living on the network, with the changes of topological structure associated to the small-world features. It turns out that whereas the emergence of short path lengths is associated to a modification of the Lyapunov spectrum, the transition from the small world to the random regimen is characterized by the conditional exponent entropies (for details see Araújo et al. ).

## 4 The logic approach to network dynamics

Thomas and collaborators \[Thomas & D’Ari, 1990\] \[Thomas, 1991\] \[Thieffry et al., 1993\] \[Snoussi & Thomas, 1993\] \[Thomas et al., 1995\] have developed logical tools to analyze biological regulatory networks. This treatment seems particularly appropriate to deal with regulatory networks where most interactions are characterized by a threshold and a saturation plateau. The regulator is usually inefficient below a threshold concentration and its effect rapidly levels off above the threshold.

The elements of the network and their interactions are represented by discrete variables, functions and parameters. Because some variables have several distinct actions, one often needs to consider more than two logical levels. One may also consider the threshold values $`s^{(i)}`$ separating the logical values. Then, the most general logical variable takes values in the set

$$\{0,s^{(1)},1,s^{(2)},2,\mathrm{}\}$$
(73)
The state of the system is described by the state vector $`xyz\mathrm{}`$ containing the logical values of the variables and the evolution of the system by a vector of functions $`XYZ\mathrm{}`$ representing the excitatory or inhibitory interactions in the network. For example, for the graph of interactions

$$\underset{\underset{-}{}}{\underset{|}{x}\stackrel{+}{}y\stackrel{}{}\underset{|}{z}}$$
(74)
the logical functions are

$$X=\stackrel{}{z},Y=x,Z=\stackrel{}{y}$$
(75)
$`X=0`$ means that the product $`x`$ is not being produced and $`X=1`$ means that the product $`x`$ is being produced, not that $`x=1`$ immediately. For example, if the initial state is $`000`$, after a certain time the state might become $`100`$ or $`001`$, depending on whether the time delay $`\left(t_x\right)`$ for production of $`x`$ is smaller or larger than the time delay $`\left(t_z\right)`$ for the production of $`z`$. Notice that although the image $`XYZ`$ of the state $`000`$ is $`101`$, the next state is either $`100`$ or $`001`$ because it is highly unlikely that $`t_x=t_z`$. The state $`001`$ is a stable state because it equals its image, whereas $`100`$ further evolves to $`110`$ or $`101`$, of which $`110`$ is stable but $`101`$ is not, etc.

A central role in this formalism is played by the oriented circuits (closed series of interactions where each element appears only once), called feedback loops. Feedback loops are positive or negative according to whether they have an even or odd number of negative interactions. Positive loops generate multistability and negative loops generate homeostasis, that is, the variables in the loop tend to middle range values, with or without oscillations.

The role of the parameters in the logical approach is to allow for the distinction between weak and strong interactions. Therefore the parameters are actually real values, not logical variables. For example, in the following interaction graph

$$x_{}^+y^+$$
(76)
$`x`$ has two values $`(0,1)`$ because it has one action only and $`y`$ has three values $`(0,1,2)`$ because it has two actions and thus two thresholds. The corresponding logical functions are

$$X=\stackrel{}{y^1},Y=x^1+y^2$$
(77)
in which all the variables are Boolean. The function $`Y`$ means that if $`x`$ is above its threshold or $`y`$ above the second threshold ($`x^1=1`$ or $`y^2=1`$) then $`y`$ is going to be produced. However we might give different weights to the variables by writing

$$X=d_x\left(K_1\stackrel{}{y^1}\right),Y=d_y(K_2x^1+K_3y^2)$$
(78)
where the $`K_i^{^{}}s`$ are real variables and $`d_x`$ and $`d_y`$ operators that discretize the value in brackets. Here the $`+`$ sign is the real algebraic sum not the logical one.

States with variable values $`xyz\mathrm{}`$ that involve only logical values are called regular states and those which involve some threshold values are called singular states. It is found that each feedback loop can be characterized by a singular logical state located on the thresholds. For appropriate parameter values, this state is stationary and the corresponding loop functional. In this context, functional means that if the loop is positive it actually produces multistability and if negative it generates homeostasis. A technique for the analysis of the network consists in dissociating it into its feedback loops and checking the dynamics of each loop.

Other logical approaches have been developed to analyze networks: the Boolean network models \[Kaufman, 1993\] \[Somogyi & Sniegoski, 1996\], hybrid models using logical and continuous variables \[Glass, 1975\] \[Lewis & Glass, 1991\] and rule-based formalisms \[Brutlag et al., 1991\] \[Meyers & Friedland, 1984\] \[Shimada et al., 1995\] \[Hofestädt & Meineke, 1995\].

A general problem in the logical approaches to network dynamics is to establish the correspondence of the logical description to the corresponding set of differential equations. For the multilevel approach (with parameters) of Thomas and collaborators, if the thresholds are represented by Hill functions, the correspondence becomes perfect when the Hill functions become very steep. The problem of establishing the correspondence between logical descriptions and smooth dynamical systems has also been addressed in other contexts. For example, in \[Martins & Vilela Mendes, 2001\] a correspondence is established between a type of neural networks and a logical programming language and in \[Martins et al., 2001\] a correspondence between controlled dynamical systems and context-dependent languages. Similar techniques may be used for general networks \[Aguirre et al., 2004\].

References

Abarbanel, H. D. I. & Rouhi, A.  “Hamiltonian structures for smooth vector fields”, Phys. Lett. A 124, 281-286.

Afraimovich, V. & Fernandez, B.  “Topological properties of linearly coupled expanding map lattices”, Nonlinearity 13, 973-993.

Aguirre, C., Martins, J. & Vilela Mendes, R.  “From differential equations to multilevel logic by grammatical inference”, in preparation.

Albert, R., Jeong, H. & Barabási, A.  “Error and attack tolerance of complex networks”, Nature 406, 378-382.

Albert, R. & Barabási, A.-L.  “Statistical mechanics of complex networks”, Rev. Mod. Phys. 74, 47-97.

Aldridge, J. & Pye, E. K.  “Cell density dependence of oscillatory mechanisms”, Nature 259, 670-671.

Amiet, J.-P. & Huguenin, P.  “Generating functions of canonical maps”, Helvetica Physica Acta 53, 377-397.

Araújo, T. & Vilela Mendes, R.  “Function and form in networks of interacting agents, Complex Systems 12, 357-378.

Araújo, T., Vilela Mendes, R. & Seixas, J. ; “A dynamical characterization of the small-world phase”, Phys. Lett. A 319, 285-289.

Araújo, V.  “Attractors and time averages for random maps”, Ann. Inst. H. Poincaré Anal. Non-linéaire 17, 307-369.

Barreira, L.  “Hyperbolicity and recurrence in dynamical systems: A survey of recent results”, math.DS/0210267, Resenhas 5, 171-230.

Boccaletti, S., Burguete, J., González-Viñas, W., Mancini, H. & Valladares, D. (Eds.)  “Space-time chaos: Characterization, control and synchronization” (World Scientific, Singapore)

Boldrighini, C., Bunimovich, L. A., Cosimi, G., Frigio, S. & Pellegrinotti, A.  “Ising-type and other transitions in one-dimensional coupled map lattices with sign symmetry”, J. Stat. Phys. 102, 1271-1283.

Bonatti, C., Diaz, L. J. & Viana, M.  “Dynamics beyond uniform hyperbolicity. A global geometric and probabilistic approach”, IMPA preprint A-221.

Boots, M. & Sasaki, A.  “Small worlds and the evolution of virulence: infection occurs locally and at a distance”, Proc. Roy. Soc. London B 266, 1933-1938.

Brutlag, D. L., Galper, A. R. & Mills, D. H.  “Knowledge-based simulation of DNA metabolism: Prediction of enzime action”, Computer Applications in the Biosciences 7, 9-19.

Buck, J.  “Synchronous rhythmic flashing of fireflies”, Quart. Rev. Bio. 63, 265-289.

Bunimovich, L. A. & Sinai, Ya.  “Spacetime chaos in coupled map lattices”, Nonlinearity 1, 491-516.

Carreira, A., Hongler, M. O. & Vilela Mendes, R.  “Variational formulation and ergodic invariants”, Phys. Lett. A 155, 388-396.

Carvalho, R., Fernandez, B. & Vilela Mendes, R.  “From synchronization to multistability in two coupled quadratic maps”, Phys. Lett. A 285, 327-338.

Chua, L. O. & Yang, L. \[1988a\] “Cellular neural networks: Theory”, IEEE Trans. Circuits Syst. 35, 1257-1272.

Chua, L. O. & Yang, L. \[1988b\] “Cellular neural networks: Applications”, IEEE Trans. Circuits Syst. 35, 1273-1290.

Cohen, M. A. & Grossberg, S.  “Absolute stability of global pattern formation and parallel memory storage by competitive neural networks”, IEEE Transactions on Syst., Man and Cybern. 13, 815-826.

Cohen, M. A.  “The construction of arbitrary stable dynamics in nonlinear neural networks”, Neural Networks 5, 83-103.

Colli, E.  “Infinitely many coexisting strange attractors”, Ann. Inst. H. Poincaré Anal. Non-linéaire 15, 539-579.

Coutinho, R. & Fernandez, B. \[1997a\] “Extended symbolic dynamics in bistable coupled map lattices: Existence and stability of fronts”, Physica D 108, 60-80.

Coutinho, R. & Fernandez, B. \[1997b\] “On the global orbits in a bistable coupled map lattice”, Chaos 7, 301-310.

Crehan, P.  “On the local Hamiltonian structure of vector fields”, Mod. Phys. Lett. A 9, 1399-1405.

Cross, M. C. & Hohenberg, P. C.  “Pattern formation outside of equilibrium”, Rev. Mod. Phys. 65, 851-1112.

Crutchfield, J. P.  “The calculi of emergence: Computation, dynamics and induction”, Physica D 75, 11-54.

Crutchfield, J. P. & Young, K.  “Inferring statistical complexity”, Phys. Rev. Lett. 63, 105-108.

Crutchfield, J. P. & Feldman, D. P.  “Statistical complexity of simple 1D spin systems”, Phys. Rev. E 55, 1239R-1243R.

de Jong, H.  “Modeling and simulations of genetic regulatory systems: A literature review”, J. Comp. Biology 9, 67-103.

Dorogovtsev, S. N. & Mendes, J. F.  “Evolution of networks: From biological nets to the internet and WWW” (Oxford Univ. Press, Oxford)

Doyne Farmer, J.  “A Rosetta stone for connectionism”, Physica D 42, 153-187.

Duarte, J. T. & Vilela Mendes, R.  “Deformation of Hamiltonian dynamics and constants of motion in dissipative systems”, J. Math. Phys. 24, 1772-1778.

Dutta, M., Nusse, H. E., Ott, E., Yorke, J. A. & Yuan, G.  “Multiple attractor bifurcations: A source of unpredictability in piecewise smooth systems”, Phys. Rev. Lett. 83, 4281-4284.

Fernandez, B. & Guiraud, P.  “Route to chaotic synchronization in coupled map lattices: Rigorous results”, Discrete Contin. Dyn. Systems \- Ser. B 4, 435-455.

Feudel, U., Grebogi, C., Hunt, B. R. & Yorke, J. A.  “Maps with more than 100 coexisting low-period periodic attractors”, Phys. Rev. E 54, 71-81.

Feudel, U. & Grebogi, C.  “Multistability and the control of complexity”, Chaos 7, 597-604.

Fiedler, B. & Gedeon, T.  “A class of convergent neural network dynamics”, Physica D 111, 288-294.

Fogelman Soulié, F., Mejia, C., Goles, E. & Martinez, S.  “Energy functions in neural networks with continuous local functions”, Complex Systems 3, 269-293.

Freeman, W. J.  “Tutorial in Neurobiology: From single neurons to brain chaos”, Int. J. of Bifurcation and Chaos 2, 451-482.

Freidlin, M. I. & Wentzell, A. D.  “Random perturbations of dynamical systems” (Springer-Verlag, New York)

Freidlin, M. I. & Wentzell, A. D.  “Random perturbations of Hamiltonian systems”, Memoirs of the American Mathematical Society no. 523.

Gade, P. M. & Hu, C.-K. ; “Synchronous chaos in coupled map lattices with small-world interactions”, Phys. Rev. E 62, 6409-6413.

Gambaudo, J. M. & Tresser, C.  “Simple models for bifurcations creating horseshoes”, J. Stat. Phys. 32, 455-476.

Gielis, G. & MacKay, R. S.  “Coupled map lattices with phase transition” Nonlinearity 13, 867-888.

Glass, L.  “Classification of biological networks by their qualitative dynamics”, J. Theor. Biology 54, 85-107.

Golomb, D. & Hansel, D.  “The number of synaptic inputs and the synchrony of large sparse neuronal networks”, Neural Comput. 12, 1095-1139.

Gouzé, J.-L.  “Positive and negative circuits in dynamical systems”, J. of Biological Systems 6, 11-15.

Grossberg, S.  “Nonlinear neural networks: Principles, mechanisms and architectures”, Neural Networks 1, 17-61.

Hanson, J. E. & Crutchfield, J. P.  “Computational mechanics of cellular automata: An example”, Physica D 103, 169-189.

Henson, S. M.  “Multiple attractors and resonance in periodically forced population models”, Physica D 140, 33-49.

Hofestädt, R. & Meineke, F.  “Interactive modeling and simulation of biochemical networks”, Computers in Biology and Medicine 25, 321-334.

Hopfield, J. J.  “Neurons with graded response have collective properties like those of two-state neurons”, Proc. Nat. Acad. Sciences USA 81, 3088-3092.

Hoppensteadt, F. C. & Izhikevich, E. M.  “Weakly coupled neural networks” (Springer, NY)

Howse, J. W., Abdallah, C. T. & Heileman, G. L.  “Gradient and Hamiltonian dynamics applied to learning in neural networks”, in Advances in Neural Information Processing Systems 8 (MIT Press)

Jiang, Z. & McCall, M.  “Numerical simulation of a large number of coupled lasers”, J. Opt. Soc. Am. 10, 155.

Jiang, M. & Pesin, Ya.  “Equilibrium measures for coupled map lattices: Existence, uniqueness and finite-dimensional approximations”, Comm. Math. Phys. 193, 675-711.

Kaneko, K. (Ed.)  “Theory and applications of coupled map lattices”  (Wiley, London)

Kaneko, K.  “Dominance of Milnor attractors and noise-induced selection in a multiattractor system”, Phys. Rev. Lett. 78, 2736-2739.

Kaneko, K. & Tsuda, I.  “Chaos and beyond” (Springer, Berlin)

Katok, A. & Hasselblatt, B.  “Introduction to the modern theory of dynamical systems” (Encyclopedia of Mathematics and its Applications 54, Cambridge Univ. Press)

Kauffman, S. A.  “The origins of order: Self-Organization and selection in evolution” (Oxford Univ. Press, New York)

Keeling, M. J.  “The effects of local spatial structure on epidemiological invansions”, Proc. Roy. Soc. London B 266, 859-867.

Kifer, Yu. I.  “On small random perturbations of some smooth dynamical systems”, Mat. USSR-Izv. 8, 1083-1107.

Kohn, K. W.  “Molecular interaction map of the mammalian cell cycle control and DNA repair systems”, Molecular Biol. of the Cell 10, 2703-2734.

Kulkarni, R. V., Almaas, E. & Stroud, D.  “Evolutionary dynamics in the Bak-Sneppen model on small-world networks”, cond-mat/9905066.

Kuramoto, Y.  “Chemical Oscillations, Waves and Turbulence” (Springer, Berlin)

Kuramoto, Y.  “Collective synchronization of pulse-coupled oscillators and excitable units”, Physica D 50, 15-30.

Lago-Fernandez, L. F., Huerta, R., Corbacho, F. & Siguenza, J.  “Fast response and temporal coherent oscillations in small world networks”, Phys. Rev. Lett. 84, 2758-2761.

Laurent, M. & Kellershohn, N.  “Multistability: a major means of differentiation and evolution in biological systems”, Trends in Biochem. Sci. 24, 418-422.

Lewis, J. E. & Glass, L.  “Steady states, limit cycles and chaos in models of complex biological networks”, Int. J. Bifurcation and Chaos 1, 477-483.

Lima, R. & Vilela Mendes, R.  “Stability of invariant circles in a class of dissipative maps”, Nonlinear Analysis 12, 1061-1067.

Lloyd, A. L. & May, R. M.  “How virus spread among computers and people”, Science 292, 1316-1317.

MacKay, R. S.  “Mode-locking and rotational chaos in networks of oscillators: A mathematical framework”, J. Nonlinear Sci. 4, 301-314.

Martins, J. F., Dente, J. A., Pires, A. J. & Vilela Mendes, R.  “Language identification of controlled systems: Modeling, control and anomaly detection”, IEEE Trans. on Syst. Man and Cybernetics 31, 234-242.

Martins, J.F. & Vilela Mendes, R.  “Neural networks and logical reasoning systems: A translation table”, Int. Journal of Neural Systems 11, 179-186.

May, P. & May, E.  “Twenty years of p53 research: Structural and functional aspects of the p53 protein”, Oncogene 18, 7621-7636.

Meyers, S. & Friedland, P.  “Knowledge-based simulation of genetic regulation in bacteriophage lambda”, Nucleic Acids Research 12, 1-9.

Newhouse, S.  “Nondensity of Axiom A(a) on $`S^2`$ ”, Proc. A. M. S. Symp. Pure Mat. 14, 191-202.

Newhouse, S.  “Diffeomorphisms with infinitely many sinks”, Topology 13, 9-18.

Newhouse, S.  “The abundance of wild hyperbolic sets and non-smooth stable sets for diffeomorphisms”, Publ. Math. I.H.E.S. 50, 101-151.

Palis, J. & Viana, M.  “High dimension diffeomorphisms displaying infinitely many periodic attractors”, Annals of Math. 140, 207-250.

Pastor-Satorras, R. & Vespignani, A.  “Immunisation of complex networks”, Phys. Rev. E 65, 036104.

Pastor-Satorras, R., Rubi, M. & Diaz-Guilera, A.  (Eds.); “Statistical mechanics of complex networks” (Springer, Berlin)

Pecora, L. M. & Carroll, T. L.  “Synchronization in chaotic systems”, Phys. Rev. Lett. 64 (1990) 821-824; “Driving systems with chaotic signals”, Phys. Rev. A 44 (1991) 2374-2383.

Pecora, L. M., Carrol, T. L. & Heagy, J. F.  “Synchronization in chaotic systems, Concepts and applications” in Handbook of Chaos Control, ch. 10, H. G. Schuster (Ed.), (Wiley-VCH, Weinheim)

Pesin, Ya.  “Dimension theory in dynamical systems: Contemporary views and applications” (Chicago Lectures in Mathematics, Chicago Univ. Press)

Peskin, C. S.  “Mathematical Aspects of Heart Physiology”, Courant Inst. Math. Sciences, New York.

Pikovsky, A., Rosenblum, M. & Kurths, J. ; “Synchronization. A universal concept in nonlinear sciences” (Cambridge Univ. Press)

Plahte, E., Mestl, T. & Omholt, S. W.  “Stationary states in food web models with threshold relationships”, J. of Biological Systems 3, 569-577.

Pontryagin, L. S.  “On dynamical systems close to Hamiltonian systems” (Russian), Zh. Eksp. Teor. Fiz. 4, 234-238.

Poon, L. & Grebogi, C.  “Controlling complexity”, Phys. Rev. Lett. 75, 4023-4026.

Robinson, C.  “Bifurcation to infinitely many sinks”, Comm. Math Phys. 90, 433-459.

Roels, J.  “On the local decomposition of a vector field on a symplectic surface as the sum of a gradient and a Hamiltonian field” Comptes Rendus Acad. Sci. Paris 278 SerA., 29-31.

Ruelle, D.  “Theory and experiment in the ergodic study of chaos and strange attractors” in Proc. VIII Int. Congress on Mathematical Physics, pgs. 273-282, M. Mebkhout and R. Sénéor (Eds.) (World Scientific, Singapore)

Schürmann, B.  “Stability and adaptation in artificial neural systems”, Phys. Rev. A 40, 2681-2688.

Shalizi, C. R. & Crutchfield, J. P.  “Computational mechanics: Pattern and prediction, structure and simplicity”, J. Stat. Phys. 104, 817-879.

Sharpless, N. E. & DePinho, R. A.  “p53: Good cop/Bad cop”, Cell 110, 9-12.

Shimada, T., Hagiya, M., Arita, M., Nishizaki, S. & Tan, C. L.  “Knowledge-based simulation of regulatory action in lambda phage”, Int. J. of Artificial Intelligence Tools 4, 511-524.

Skarda, C. A. & Freeman, W. J.  “How brains make chaos in order to make sense of the world”, Behavioral and Brain Sciences 10, 161-195.

Snoussi, E. H. & Thomas, R.  “Logical identification of all steady states: The concept of feedback loop-characteristic state”, Bull. Math. Biology 55, 973-991.

Snoussi, E. H.  “Necessary conditions for multistationarity and stable periodicity”, J. of Biological Systems 6, 3-9.

Somogyi, R. & Sniegoski, C. A.  “Modeling the complexity of genetic networks: Understanding multigenic and pleiotropic regulation”, Complexity 1, 45-63.

Strogatz, S. H.  “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators”, Physica D 143, 1-20.

Strogatz, S. H.  “Exploring complex networks, Nature 410, 268-276.

Taylor, W. R. et al.  “Mechanisms of G2 arrest in response to overexpression of p53”, Mol. Biol. of the Cell 10, 3607-3622.

Terman, D.  “Dynamics of singularly perturbed neuronal networks” in An Introduction to Mathematical Modeling in Physiology, Cell Biology and Immunology, Proc. Symp. in Applied Math. vol. 59, pages 1-32, Am. Math. Soc., Providence 2002.

Thieffry, D., Colet, M. & Thomas, R. ; “Formalization of regulatory networks: A logical method and its automatization”, Math. Modeling Sci. Computing 2, 144-151.

Thomas, R. & D’Ari  “Biological Feedback” (CRC Press, Boca Raton)

Thomas, R.  “Regulatory networks seen as asynchronous automata: A logical description”, J. Theor. Biology 153, 1-23.

Thomas, R., Thieffry, D. & Kaufman, M.  “Dynamical behavior of biological regulatory networks - I. Biological role of feedback loops and practical use of the concept of loop-characteristic state”, Bull of Math. Biology 57, 247-276.

Tonti, E.  “Variational formulation for every nonlinear equation”, Trieste ICTP report SMR/92-20 and references therein.

Tyson, J. J. & Ohtmer, H. G.  “The dynamics of feedback control circuits in biochemical pathways”, Progress in Theor. Biology 5, 1-62.

Vilela Mendes, R. & Duarte, J. T.  “Decomposition of vector fields and mixed dynamics”, J. Math. Phys. 22, 1420-1422.

Vilela Mendes, R. & Duarte, J. T.  “Arcs of discrete dynamics and constants of motion”, Lett. Math. Phys. 6, 249-252.

Vilela Mendes, R.  “A variational formulation for dissipative maps”, Phys. Lett. A 104, 391-395.

Vilela Mendes, R.  “Generating functions for noncanonical maps”, Lett. Math. Phys. 11, 289-292.

Vilela Mendes, R.  “Deformation stability of periodic and quasiperiodic motion in dissipative systems” in Deformation theory of algebras and structures and applications, M. Hazewinkel and M. Gerstenhaber (Eds.) (Kluwer Acad. Pub., Dordrecht)

Vilela Mendes, R. & Duarte, J. T.  “Vector fields and neural networks”, Complex Systems 6, 21-30.

Vilela Mendes, R.  “Conditional exponents, entropies and a measure of dynamical self-organization”, Phys. Lett. A 248, 167 - 171.

Vilela Mendes, R.  “Clustering and synchronization with positive Lyapunov exponents”, Phys. Lett. A 257, 132-138.

Vilela Mendes, R. \[2000a\] “Multistability in dynamical systems” in ”Dynamical systems: From crystal to chaos”, J.-M.Gambaudo, P. Hubert, P. Tisseur and S. Vaienti (Eds.), pp. 105-113 (World Scientific, Singapore)

Vilela Mendes, R. \[2000b\] “Characterizing self-organization and coevolution by ergodic invariants”, Physica A 276, 550-571.

Vilela Mendes, R.  “Structure generating mechanisms in agent-based models”, Physica A 295, 537-561.

Vogelstein, B., Lane, D. & Levine, A. J.  “Surfing the p53 network”, Nature 408, 307-310.

Vousden, K. H.  “p53: Death star”, Cell 103, 691-694.

Wang, X.-J.  “The Newhouse set has a positive Hausdorff dimension”, Comm. Math. Phys. 131, 317-332.

Watts, D. J. & Strogatz, S. H.  “Collective dynamics of small-world networks”, Nature 393, 440-442.

Watts, D. J.  “Small Worlds: The dynamics of networks between order and randomness” (Priceton Univ. Press, Princeton NJ)

Wei, G. W., Zhan, M. & Lai, C.-H.  “Tayloring wavelets for chaos control”, Phys. Rev. Lett. 89, 284103.

Weigel, R. & Atlee Jackson, E.  “Experimental implementation of migrations in multiple attractor systems”, Int. J. of Bifurcation and Chaos 8, 173-178.

Winfree, A. T.  “Biological rhytms and the behavior of populations of coupled oscillators”, J. Theor. Bio. 16, 15-42.

Woods, D. B. & Vousden, K. H.  “Regulation of p53 function”, Exp. Cell Reserach 264, 56-66.

Wuensche, A.  “Basins of attraction in network dynamics: A conceptual framework for biomolecular networks”, SFI Working paper 02-02-004 to appear in Modularity in Development and Evolution, G. Schlosser and G. P. Wagner (Eds.) (Chicago Univ. Press.)

Yang, X.-S.  “Chaos in small-world networks”, Phys. Rev. E 63, 046206.