- Email: [email protected]
Full characterization of a strange attractor
Physica D 48 (1991) 65-90 North-Holland
FULL CHARACTERIZATION
OF A STRANGE A’ITRACTOR
CHAOTIC DYNAMICS IN LOW-DIMENSIONAL REPLICATOR SYSTEMS
Wolfgang SCHNABL, Peter F. STADLER, Christian FORST and Peter SCHUSTER Institut fiir Theoretische Chemie der Uniuersitiit Wien, WiihringerstralJe 17, A-1090 Wien, Austria Received 17 January 1990 Revised manuscript received 8 September 1990 Accepted 8 September 1990 Communicated by A.T. Winfree
Two chaotic attractors observed in Lotka-Volterra equations of dimension n = 3 are shown to represent two different cross-sections of one and the same chaotic regime. The strange attractor is studied in the equivalent four-dimensional catalytic replicator network. Analytical expressions are derived for the Lyapunov exponents of the flow. In the centre of the chaotic regime the strange attractor was characterized by the Lyapunov dimension (D, = 2.06 f 0.02) and the Rtnyi fractal dimensions. The set corresponding to the attractor represents a multifractul. Its singularity spectrum is computed. One route in parameter space leading into the chaotic regime and crossing it was studied in detail. It leads to a Feigenbaum-type cascade of bifurcations. Before the chaos is fully developed the dynamic system passes an internal crisis through a period-two tangent bifurcation. Then a second sequence of period-doubling bifurcations leading to another chaotic regime is observed, which eventually ends in a crisis and then after a rather complicated periodic regime the attractor finally disappears. Trajectories in the intermediate periodic regime show transient chaotic bursts of intermittency type. A series of one-dimensional maps is derived from a properly chosen Poincare cross-section which illustrates structural changes in the attractor. Mutations are included in the catalytic replicator network and the changes in the dynamics observed are compared with the predictions of an approach based on perturbation theory. The most striking result is the gradual disappearance of complex dynamics with increasing mutation rates. The parameter space spanned by the mutation rate and selected parameters boundaries
of the replication network is split in regions of complex periodic behaviour near the chaotic regime. of these regions are likely to occur but the data available do not allow definitive conclusions yet.
Fractal
1. Introduction The search for chaotic dynamics became an issue in many scientific disciplines from physics to biology. The relevance of strange attractors in modeling real systems and the problem how to trace chaos in the experimental data were important questions already five years ago [l] and remained in the centre of interest since then [2-41. Numerical evidence for strange attractors was found in many model equations but unlike simple attractors they cannot be characterized by a few coordinates or coordinates and frequencies. In the pioneering works of several groups efficient techniques to analyse chaotic dynamics and fractal sets were developed. (For a recent review of the available repertoire of methods see refs. [5, 61.) In order to provide a sufficiently detailed basis for the comparison of model calculations and experimental data full characterization of strange attractors is required. We present such a study on a chaotic attractor which appears in the solutions of a differential equation which is of interest in biophysics and in theoretical ecology as well as in other disciplines. Autocatalytic reaction networks which involve two classes of catalytic processes - autocatalytic instruction for reproduction as well as material- and process-specific catalysis- were postulated as models for 0167-2789/91/$03.50
0 1991-
Elsevier
Science
Publishers
B.V. (North-Holland)
66
W. Schnabl et al. /Full characterization of a strange attractor
studies on prebiotic and early biological evolution scenarios [7, 81. Now they serve also as simple models showing the characteristic features of nonlinear dynamics in very different fields like molecular biology, population genetics, theoretical ecology or dynamical game theory [9]. Qualitative analysis and systematic numerical investigations revealed a very rich dynamics of the corresponding kinetic differential equations [lo]. The name replicator equations was coined for these systems because they reveal several unique - in the standardized form to be discussed in section 2 - fairly easy to prove features [ll]. The four-species replicator system shows chaotic dynamics for certain parameter values [121. According to Hofbauer [13] this autocatalytic network with four species is - apart from a transformation of the time scale - equivalent to a three-species Lotka-Volterra equation. Possible existence of very complex dynamical behaviour in Lotka-Volterra models was predicted by Smale [14]. Indeed, two different strange attractors were reported for the Lotka-Volterra system [15-181, which describes communities of one predator species living on two prey species. In this contribution we shall show that the two strange attractors lie on different cross-sections through one and the same chaotic regime. The number of degrees of freedom is minima1 for the occurrence of strange attractors and hence the chaotic regime in question houses the simplest strange attractors possible in Lotka-Volterra and replicator equations - replicator systems have normalized variables and hence the number of independent degrees of freedom is the number of species but one.
2. The replicator model Replicator equations are based on the simplest possible molecular mechanism of catalysed replication and mutation. The basic reactions are of the form (A) +Ij+Ii%
I,+I,+Ij,
i,j,k=l,...,
n.
(1)
Herein Ij is the template which is replicated and Ii the catalyst, (A) represents an overall notation for substrates needed for replication. It is assumed that the substrates are present in excess- their concentrations are buffered and hence do not change during the course of reactions. There are n3 possible reactions producing the n different species I, (k = 1,. . . , n) of the autocatalytic reaction network. Clearly it is very hard-if not impossible -to handle such a large number of parameters. In order to be able to analyze the reaction network simplifying assumptions are inevitable. It is straightforward to treat catalysed replication and mutation as independent superpositions: mutation frequencies depend on template (Ii) and target (I,), but not on the catalyst (Ii). The remaining 2n2 rate parameters are properly understood then as elements of two n x n matrices: a replication rate matrix A k (AjJ and a mutation matrix Q A (Qkj). The rate constants of the individual reactions are products of a replication and a mutation factor. Qkj .Aji, expressing the following sequence of processes: Ii catalyses the replication of I j which yields I k as an error copy. Error-free replication of species Ii- described by eq. (1) with I, = Ij - is assumed to occur with frequency Q,,. A mutation from Ij to I, (k #j) occurs with frequency Qkj. Every copy has to be either correct or a mutant and thus Q is a stochastic matrix (Xi =, Qkj = 1 Vj = 1,. . . , n). Concentrations of individual species are denoted by ck = [Ik]. As variables we use relative-or normalized - concentrations, xk = c&j”= icj, which obviously fulfil X; = ixk = 1.
U! Schnabl et al. /Full
characterization
Mass action kinetics applied to the autocatalytic reaction-mutation form dxk -
dt
61
of a strange attractor
network (1) yields an ODE of the
n =f,=
QkiAijxixj-x,@,
c
n,
k = l,...,
i,j=l
where an automatically adjusted dilution flux
k Aijxixj
CD=
i,j=l
is introduced in order to keep the sum of relative concentrations meaningful domain of eq. (2) is the simplex
s, A
(x,,...,
normalized
to unity. The physically
txi=l .
xn)l xi20,
i
I
i=l
If AijrO, i,j= l,..., n then S, is a compact, foreward invariant set for the ODE (2). An intensity of mutation can be measured by the mean mutation rate E: 1 E = n(n - 1)
5 i~i
Qij*
(3)
i,j=l
Throughout
this paper we use a particularly simple type of the mutation matrix
Qij= 1 - (n - l)e,
if i=j
and
Qij=e,
if
i#j.
It represents a special case of the house of cards model of mutation [19, 201, frequently used in theoretical population genetics. The important special case of error-free replication is represented by a unit mutation matrix, Q A IQ,, = S,,) with 6,, being the Kronecker symbol. Then eq. (1) is turned into the well known replicator equation :
~ AkiXi i=l
~
AijXiXj
(4)
i,j=l
The linear replicator equation (4) for n is flow-equivalent to the n - 1 species Lotka-Volterra
,
k=l,...,
n-l.
The variables are related by means of the transformation yk=xk/x,,
fork=l,...,n-1.
equation:
(5) introduced by Hofbauer [13]:
(6)
68
W. Schnabl et al. /Full
3. The chaotic regime in parameter
characterization
of a strange attractor
space
Two distinct chaotic attractors were found in Lotka-Volterra equations for three species: Vance [17] discovered a “quasi-cyclic” trajectory in a one-predator two-prey model and Arneodo, Coullet and Tresser [15] found a one-parameter family of strange attractors (ACT attractor). We used Hofhauer’s transform to construct the replicator equations which are equivalent to the models mentioned above. Then a barycentric transform [211 was applied to Vance’s model in order to shift the interior fixed point into the centre of the simplex S, - the point c = f
e
bkizi
-
i=l
k
bijziz/
i,j=l
has an interior fixed point i = (i,, i,, . . . , 2,). It can be shown [221 that there exists a di’omorphism which relates this replicator equation to another one in the notation of eq. (4) with ~~~~~
and
f=i(l,l,...,
1).
I
Both the ACT family and the Vance attractor (V) correspond to nonrobust phase portraits because of zero or almost zero (see e.g. A,, or A,, in A,) off-diagonal elements in the reaction matrices A. In order to make the numerical computations easy to reproduce all parameter values were rounded to three digits. Then the Vance model has a replication matrix 0
A,=
0.537 -0.535 0.536
0.063 0 0.38 - 0.032
0 - 0.036 0 - 0.004
0.437 - 0.001 0.655 0
I
(7)
and the ACT attractors are found at
J 0 A ACT=
1.1 -0.5 \ 1.7+/L
0.5 0 1 -1-p
-0.1 - 0.6 0 -0.2
0.1 '
0 0. 0,
(8)
Note that we use a parameter p which corresponds to p - 1.5 in ref. [15]. In order to find out whether or not these two attractors are related in the twelve-dimensional parameter space we performed a systematic numerical survey in the two-dimensional subspace (p, VI which connects the two models:
!
0
A(P, v) =
1 .l - 0.563~ 0.035v - 0.5 1.7 +/_L- 1.164~
0.5 - 0.437v 01 - 0.62~ - 1 - /L + 0.968~
-0.1 + O.lu 0.6 0 + 0.564~ - 0.2 + 0.196~
(9)
W. Schnabl et al. /Full
69
characterization of a strange attractor
1.0
.a .6 .4
.2
0.
-.2
-.4
-.6
-.a
t -1.0
1
! -1.5
’
’
’
’
! -1.0
’
’
’
’
! -.5
’
’
’
’
!
’
’
’
0
’ .5
v
Fig. 1. A two-dimensional cross-section through the chaotic regime of eq. (4). The cross-section is defined by the two parameters w and Y in the replication matrix A (9). The density of hatching refers to the dynamics of the attractor. In the high-density zone chaotic dynamics is observed. The lowest density indicates simple limit cycles, which undergo sequences of period doublings in the intermediate region. Numbers 1 to 14 indicate regions of qualitatively different dynamics on the boundary of the simplex S,. The corresponding phase portraits are shown in fig. 2.
We note that all three matrices of replication constants fulfil E;= ,Aij = 0.5 - a property which will be used together with the existence of a fixed point in the centre of the simplex later on. In the (IL, v) plane the ACT model lies on the line - 0.2 I p 5 0.05 and v = 0. The Vance attractor is close to the point (0,l). Both models are embedded in a Z-shaped two-dimensional manifold of strange attractors (fig. 1). In scanning through the manifold we could not find any sudden changes in the shape of the attractor. The mean revolution time also varies in apparently continuous manner all over the chaotic regime. A common feature of all phase portraits exhibiting chaotic behaviour (attractive or only transient) is a saddle focus Pi,, on the [l, 2,3] plane, which represents the stable manifold. The unstable manifold thus points into the interior of the simplex. We found that Silnikov’s condition [23] is always satisfied: the real positive eigenvalue is larger in magnitude than the real parts of the two complex conjugate eigenvalues. Although there is no homoclinic orbit possible for Plz3, a close relation to !%lnikov’s condition is evident, as was pointed out by Arneodo et al. [161: there is a heteroclinic orbit involving three saddles instead of
70
W. Schnabl et al. /Full
characterization
of a strange attractor
one and this orbit replaces the homoclinic one. For generalized Lotka-Volterra models with additional cubic terms Samardzija and Greller [24] found a different type of chaos arising from a fractal torus. For small values of p the fixed point in the centre of S, (c) is stable. Chaos is reached via a Hopf bifurcation and a subsequent cascade of period-doubling bifurcations. At first the limit cycle, and then the chaotic attractor, both grow steadily with increasing p values. At a parameter value around p = 0.21 the chaotic attractor collides with the stable manifold of a saddle on the [l, 2,41 plane and disappears suddenly. This behaviour is known as limit crisis [251: at the bifurcation point we find a heteroclinic orbit connecting the [l, 41 saddle, a saddle on the [l, 21 edge, and the saddle focus in the [l, 2,3] plane (see
A
a
Fig. 2. A sketch of fourteen flows on the boundaries of the simplex S, which can coexist with the chaotic attractor The numbers 1 to 14 are used as markers in fig. 1 to indicate where the phase portraits occur in parameter space.
in the interior.
71
W. Schnabl et al. /Full characterization of a strange attractor
fig. 7 below). The evidence for the orbit is derived from numerical integration. We cannot provide an analytical proof for the saddle connection between the saddle focus and Pi,. The saddle connections on the boundary follow from the invariance and the subsimplices and the classification of the flow of three-component replicator equations [26, 271. Beyond the limit crisis we find a field of transient chaos: the trajectory seems to draw an image of the chaotic attractor which in parameter space lies on the other side of the limit crisis before it ends up in a stable fixed point. For small values of v the frequency of the oscillations increases with decreasing parameter values-revolution times become shorter and shorter. Eventually the attractor collides with the stable manifold of the saddle on the [1,2] edge. The invariant hyperplane Ci,ixk = 1 is a repellor in this area and therefore integration is extremely sensitive to numerical errors. The chaotic attractor in the interior of the simplex S, coexists with a rather large variety of robust flows on the surface. Fig. 2 shows fourteen boundary flows consistent with chaos on the (II, v) plane. Fig. 3 sketches the phase portraits on the ACT and the Vance model together with a typical chaotic trajectory for each system. Direct evidence from the inspection of computed trajectories was supplemented by Fourier transformation [28]. Some typical examples of Fourier spectra in periodic and chaotic regimes of our attractor are shown in fig. 4. They are used as a sensitive indicator of the occurrence of chaos. In the four spectra chosen we observe a complete transition from a periodic window into the chaotic regime caused by a small change in the parameters CAP < 0.02).
4. Measures of the chaotic attractors The most common measure for a chaotic attractor is the spectrum of Lyapunov exponents [291. Let J(t) = dF(x(x(O), t)) be the Jacobian of the vector field 9 along a trajectory with starting point x(O).
x3
x2
A Fig, 3. Phase portraits and chaotic ACT model (p = 0, v = 0).
trajectories in the Vance
x3
x2
0 (7) and the ACT model (8). (A) Vance
model (CL= 0, v = 1) and (B)
W. Schnabl et al. /Full
characterization
of a strange attractor
I
u=-0.0536
8
6
I
4
2
0
-2 ~ 01
.02
.03
.04
.05
.06
.a7
.oa
.09
.I0
Frequency
A log[Al u=-0.0529
.Ol
.02
.03
.04
.05 6
.06
.07
.oa
.09
Frequency
Fig. 4. Fourier spectra of selected trajectories of the replicator model (9) with the parameter values v = 0 and (A) p = -0.0536, period-3 window, (B) -0.0529 period-3 window after period doubling, (C) -0.05 at the borderline of chaos and (D) -0.04, well developed chaos.
73
W. Schnabl et al. /Full characterization of a strange attractor
logCAl
logCAl a
.01
.02
.03
.04
.05
D Fig. 4.
Continued
.06
.07
.OB
.09
Frequency
74
W. Schnabl et al. /Full
characterization
of a strange attractor
Now we define an operator -Y(t)
A’? exp/oj(r)dT).
The time ordering operator ? has to be introduced since Jacobian matrices evaluated at different times - J(T) and J(T’) - usually do not commute. Let pi(f) be the eigenvalues of 1(t), then the Lyapunov exponents are given by *_= 1
lnh-%(t)ll
lim
f--to=
t
(11)
*
In most cases it is not possible to calculate Lyapunov exponents for the flows of differential equations analytically. An appropriate numerical algorithm has been given by Bennetin et al. 1301, a well tested FORTRAN program is available in the literature [31]. In the four-species replicator system there are three physical dimensions imbedded in a fourdimensional state space. The hyperplane
is invariant. Accordingly d(t) has three eigenvectors gi lying in this plane and their components fulfil C,$) = 0. The fourth eigenvector is pointing out of the hyperplane. Thus we expect to find three physically meaningful Lyapunov exponents. If all three of them are negative the attractor is a sink. For a periodic, quasiperiodic or chaotic attractor at least one exponent has to be zero: this exponent corresponds to an eigenvector which points in the direction of the trajectory. In particular, we have for periodic orbits A, I A, < 0, A, = 0, for quasiperiodic orbits A, < 0, A, = A, = 0 and finally for a chaotic attractor one Lyapunov exponent is negative, one is zero and one is positive, A, < 0, A, = 0, A, > 0. The sum of all Lyapunov exponents fulfils the equation
(12) For error-free divF=
replication - if the mutation matrix Q is the unit matrix - the divergence div F reduces to e
Akjxj - (n + 2)
k,j=l
k
Aijxixj.
(13)
i,j=l
The mean flux 5 is defined by Aijx,(
T) xj( T) dr.
If the replicator equation is permanent
[ll, 32, 33]- or at least if the trajectory remains in the interior of
75
W. Schnabl et al. /Full characterization of a strange attractor
the simplex-we
have
Since the last line is the equation fulfilled by the coordinates Xj = lim f/‘xj(r) t-m
0
of the interior fixed point we have
dr =Zj,
where f = (a,, E,, . . . , i,)
(14)
is the interior fixed point if it exists. Its coordinates
fulfil the equations
Vi=l,...,n
iAijaj=s,
and hence we have
h
Aijij
=
ns.
i,j=l
Accordingly we find lim :itdivY(_r(r))dr=
t+m
-2s.
(15)
Finally we are left with the following relation derived from eq. (12) for the sum
kAi=--$, $ Aij. i=l
( 124
l,J=l
Since we applied the burycentric transformation we obtain
s=;$
to our model, .i?j= l/n
holds for all coordinates
and
Aij.
t,J-1
Using the numerical values of (9) the sum of the four Lyapunov exponents becomes Cf= rhi = - i. The eigenvector pointing in the direction (1, 1, 1,l) is denoted by
1 e=2:
1 1 1 . 1
il
76
W. Schnabl et al. /Full
The corresponding Lyapunov the linearized model starting
we obtain
exponent (A,) can be calculated analytically. Let us consider from x(0) = e, which to first order in time is given by
the orbit of
@(T*).
X(T) =e+TdFe+ Then
characterization of a strange attractor
for the Lyapunov
exponent
The last equation is derived straightway by explicit evaluation of the integrand as the sum over all elements of the Jacobian. One obtains the same result -Q(T) for every time T and hence A, = -5. Thus we have in case of our model:
;hi= -a
and
A,=
-i.
(16)
i=l In the chaotic regime it is sufficient exponents are obtained from A,=O,
A,=
--i-A,
and
In fig. 5 we show the largest parameter CL. From the Lyapunov
.,=2-2=2(
therefore
A,=
ll++h’n:‘)
the largest
Lyapunov
exponent
A,. All other
-i.
Lyapunov
exponents
to compute
exponent
we can calculate
(17) (A,) of the strange an estimate
attractor
for the fractal
as a function
dimension
of the
of our attractor:
ifA,>O,
=l
if A, =O,
=o
if A, < 0.
(18)
From fig. 5 we obtain A, = 8 x lo-” for the highly developed chaos at p = Y = 0 and compute a dimension of D, = 2.06 with an estimated numerical uncertainty of kO.02 at this point. In order to study the fractal nature of the attractor quantitatively we computed also generalized or RCnyi dimensions, D,, which correspond to the scaling exponents of the qth moments of the attractor’s measure [34-371. The attractor is covered by M(1) cells with edges of length 1. By p, we denote the are probability to find a point of the attractor in cell number i (i = 1,. . . , M(l)). Then Rtnyi dimensions obtained from
log cT”=‘:ap=
D, = lim 140 (q-lllogl
lim
log P(I)
1+0 (4 - 1) log 1’
(19)
W. Schnabl et al. /Full characterization of a strange attractor
-.2
Fig. 5. The largest through the chaotic
-.15
-.05
-.l
Lyapunov exponent A, of the vector regime we chose v = 0.
0
.os
flow defined
.l
by eq. (9). In order
77
.15
.2
to show a representative
CI
cross-section
where we introduced the correlation integral defined by Cq(l) G Cfl’?pp. In order to circumvent an indeterminant expression in the denominator the RCnyi dimension D, is obtained by taking the limit lim q ~ 10, = D,. For the actual calculations we rewrite the sum over the pi’s in terms of the probability ji of the trajectory to cross a box of size I around the iterate x(tj). We make use of the equation
(20) and find that the latter probability is given by [38, 391:
pj = $+rfm
k:
o( I - llXi - Xjll))
(21)
i=l
where 0 is used for the Heaviside function: @(x-x,)=0
for
x
and
0(x -xt,)
= 1 for x2x,.
Thus we are left with the easier to compute expression
(22)
W. Schnabl et al. /Full
78
for the correlation
characterization of a strange attractor
integral [40]. In the limit q -+ 1 we find the information
dimension
(23) The RCnyi dimensions for q = 0 are equivalent to the Hausdorff dimensions, for q = 2 they correspond to the correlation dimensions. The evaluation of RCnyi dimensions turned out to be extremely time consuming and therefore we restricted computations to one point right in the centre of the chaotic regime (CL= v = 0). In order to calculate the correlation integrals P(1) the trajectory was decomposed into a time series X=(&=x&)
Vi=1,2,...,N)
(24)
with t, < t, c . . . -c t,. Some ten thousand points are required which have to be distributed independently along a trajectory. This condition can be met by collecting subsequent points which are separated by at least one tenth of a full revolution on the attractor. The data set used in our computation was obtained by sampling of 120 different trajectories with randomly chosen initial conditions. A transient phase of t, = 3000 was discarded and then 1768 points were recorded in equal time intervals of At = 20 for each trajectory. The mean revolution time on the attractor is of the order 70 < t, < 80. Correlation integrals were computed for individual trajectories and for the whole data set of about 200000 points. Superposition of about 5 trajectories was already sufficient to achieve deviations from the final result in the percent range. Convergence for q 2 - 2 was generally good. Inserting the fractal set obtained in the described way into eqs. (22) and (23) we found D, = 2.02,
D, = 1.91
and
D, = 1.77.
The estimated numerical uncertainty of our values is certainly less than &0.05. The data apparently are not in perfect agreement with the Kaplan-Yorke conjecture according to which the information dimension of the attractor, D,, should be equal to the Lyapunov dimension D,. The latter is about 15% larger than the former in our case. We cannot claim with ultimate certainty that this is not a result of some hidden numerical errors, but the discrepancy exceeds our estimated error bounds. The cause of such deviations may lie in the fact that higher-order terms in the box length 1 are non-negligible even at the smallest numerically accessible I-values. In addition, the Lyapunov dimension is unlikely to be too large since it has to fulfil D, > 2 by the nature of our attractor. For a recent discussion of systematic errors in computations of the correlation integral see several contributions in ref. [41]. Fractal sets can be characterized by the scaling properties of normalized distributions or measures lying upon them [42]. We briefly repeat the basic expressions and relations. Scaling properties are expressed in terms of singularity strengths (Y and a function f(o) which reflects the density of their distribution. The spectrum of singularities is given by the possible range of (Yvalues and the density f(a). In order to define these quantities we cover the attractor with boxes of size 1 and denote the number of points inside the ith box by Ni. Then the probability that this box contains a point of the attractor’s time series (24) is given by Pi(l) = limN_JNi/N). The 1ocal singularity strength (Y~is now given as the scaling exponent Pi(Z) a PI. The singularity strength cy, thus varies with the region of the attractor on which it is taken. Its probability density - the probability to encounter a value of (Ybetween (Yand 6 + d& - is given
79
W. Schnabl et al. /Full characterization of a strange attractor
Prob{ & < (Y< & + d&l = da p( (Y) I-fC’), where f(a) is a continuous function. The set of singularity strengths (Y~itself is a fractal and f(a) can be understood as its - fractal - dimension. In other words, fractal measures are modeled by interwoven sets of singularities of strength (Y-each one having its own fractal dimension f(a). The singularity spectrum f(a) is related to the RCnyi dimensions D, since P(Z)
= C(q,Z)
=/dGp(&)
Iqail-f(@.
The dimensions are defined in the limit I + 0 and accordingly we may approximate the integral for nonzero densities p(G) by the term with the smallest exponent q& -f(h). This value a(q) is defined by
and the existence of a minimum requires
which finally yields f’(a(q)) = q and f”(a(q)) < 0. We have two relations:
D, =
&
[w(q)
and 4s) =
-f(hl))l
&[(4- l)Dq]
(25)
which allow to compute the singularity spectrum from known RCnyi dimensions or vice versa. It can be shown straightway that the Hausdorff dimension of the attractor fulfills D, =f(ama,), where (Y,, is the position of the maximum of f(a). Feigenbaum [43] recently pointed out that this formalism is related to equilibrium statistical mechanics. A new function T(q) A (q - l>D, is introduced. The singularity spectrum f(a) is the Legendre transform of T. In Feigenbaum’s formalism T(q) and q are the conjugate thermodynamic variables to f(a) and D4. The singularity spectrum of the strange attractor with I_L= Y = 0 was calculated from a series of RCnyi dimensions. For a given set of values q = qk the singularity strengths were obtained by partial differentiation: a(q) = i3~/d(u. The partial derivatives were computed by means of a five-point differentiation formula. Then we derived f(a) for discrete points (Ye= dq,) from f(%)
&fk =
akqk
-
(qk
-
l)D,k.
(26)
Limiting dimensions can be derived straightway from the singularity spectrum: D,, = cdq = +m) and D_, = a(q = -m). Both limits are related to the singularity spectrum by f(a(q = km)) = 0. In the case of our chaotic attractor we find 1.3 I (YI 3. This result can be interpreted as follows: In case the trajectory visits a volume element in the three-dimensional concentration simplex S,, it occupies a fractal set whose measure has a Hausdorff dimension of cx= 1.3 or larger. For negative values of q the estimates of the RCnyi dimensions are rather poor, because the limit in eq. (22) converges very slowly. In our case convergence was poor for q = - 1.5 and q = -2, computed results with still smaller q-values
80
U! Schnabl et al. /FUN characterization of a strange attractor
Fig. 6. The singularity spectrum of the strange attractor in the centre of the chaotic regime (p = 0, v = 0). Note that the maximum of f(a) represents the Hausdorff dimension CD,,) of the attractor. Computation of the points in the branch for positive q-values (left-hand branch) converged well. Numerical uncertainties due to incomplete convergence lie within the size of the full circles. Points at the branch for negative q-values (right-hand branch) are inaccurate because of poor convergence. The curve extrapolates towards the physically meaningful limit a = 3.
are completely unreliable. Thus the right-hand branches of the computed singularity spectra are inaccurate and we expect rather large errors for the extrapolated values of D_,. Indeed straight extrapolation of the points obtained would yield an upper limit larger than (Y= 3 (fig. 61. The correct upper limit, however, has to fulfil (Y2 3 since the trajectory is ultimately restricted to the simplex S,.
5. The routes to chaos Routes into chaotic dynamics were studied by means of Poincare maps in the plane defined by x4 = $. Distributions of points at which a long-time trajectory crosses this plane are shown in fig. 7. A Cartesian coordinate system is defined with the origin in the central fixed point c = i(l, 1, 1,ll (fig. 7) and the branch starting near c and extending towards the projection of the [l, 21-edge onto the plane x4 = f is chosen as Poincare cross-section. Its shape is very close to a straight line. Time series % = (&,, _E2,. . . , E’) are defined as successive points at which the trajectory passes the Poincart cross-section. They were recorded as a function of the parameter /L (- 0.21 < p < 0.21) with the second parameter held constant at the value of the ACT model (Y = 0). The corresponding plot is shown in fig. 8. The first part ( - 0.21 < p < - 0.01) reminds of the analogous plot of the discrete logistic map on the unit interval: X,, +, = rX,(l -X,1. Indeed there is a monotonic - almost linear - relation between the values of the constant r and the parameter p at which the characteristic period doublings and periodic windows in the chaotic regime occur (fig. 9). The sequence of events into chaos and the structure of the chaotic regime along this route show perfect correspondence to the well-known Feigenbaum sequence. A more precise characterization of the individual periodic regimes is contained in table 1.
W. Schnabl et al. /Full characterization of a strange attractor
Fig. 7. Points at which a long-time trajectory crosses (A) and three examples with v = 0 and /L = -0.1001
of the attractor the plane x 4 = $. The orlentatlon (B), 0.0001 (C) and 0.1 (D) are shown.
relative
to the cross-section
W. Schnabl et al. /Full characterization of ‘a strange attractor
82
-4
I
I! i \
0.
F’N’t
I ! ‘I’(““’
-.2
-.15
! -.l
I”‘!““!
-.05
““‘I’ 0
.05
1’ .l
!
I”“’
.15
.2
Fig. 8. A PoincarC map of the attractors with Y = 0 as a function of the parameter g. A projection x-axis defined in fig. 7 is shown. The origin (x = 0) coincides with the hxed point in the centre (c).
0.
P
of the Poincare
map onto the
-
-.02
-
-.04
-
-.06
-
-.08
-
-.lO
-
-.12
-
-.14/-‘,,, 3.4
I ,,,,,,,,,,,,,,‘,,,,,,,11, 3.5
3.6
~. 3.7
3.8
3.9
r 4.0
Fig. 9. Relation between the critical values of the parameter r in Feigenbaum’s sequence of bifurcations and periodic windows in the chaotic regime of the logistic map on the unit interval and the values of the parameter p at which the analogous events occur in the dynamics of the replicator equation.
W. Schnabl et al. /Full
Table 1 Periodic windows Periodb’
in the chaotic Range
regime (-0.105
83
of a strange attractor
< /.L < 0.1684).@ Comments
of p values end
begin 12 6 5 10 3 ff. 5 7 6 4 ff. 6 7
characterization
begin of Feigenbaum
-0.1034 - 0.0956 - 0.0770
sequence
- 0.0948 - 0.760 - 0.0632
- 0.0587
- 0.0527 - 0.0398 - 0.0359 - 0.0326
- 0.0255
- 0.0250 -0.0182 -0.0116 end of Feigenbaum
9
sequence
0.0082
3 ff. 4 6
0.0118
2 ff. 12 7
0.0544
12 8 ff. 4 6 12 6
0.0865 0.0883 0.0922 0.1057
2 ff. 6 10 4
0.1207
0.0127 0.0227 0.0448 interior
crisis
0.0772 0.0784 0.0844 0.0868 0.0892 0.0980 0.1060 0.1090 0.1102 0.1291 0.1312 0.1489 0.1543
a)In the case of very narrow windows only one parameter value is given. “Windows in which at least one period doubling event was identified are indicated
by “ff’.
As in the logistic map, the range on the PoincarC section which is visited by the chaotic trajectory grows with increasing values of the parameter p. In contrast to the former, however, a certain region around the central fixed point c is always left blank. Beyond the end of the Feigenbaum sequence we observe a broad continuation of the chaotic regime up to a sharp interior crisis at p = 0.045. Characteristic accumulation zones of points on the cross-section which continue as lines through periodic windows appear in the range above /..L= -0.01 as well. They look as if they were reflected at the lower boundary of the chaotic regime near the point c. In order to search for multiple attractors in the interior of S, several different initial values- all chosen in the plane x4 = f -were applied in the parameter range 0.05 < p < 0.17. No indication of a second attractor in the interior of S, was found. There are, however, regions in this plane which have the fixed points P,, or P3 as w-limits (see fig. 2B). At p values above the interior crisis we observe a second period doubling cascade ending in a small chaotic regime which is again interrupted by periodic windows. Chaos expands with increasing p-values and eventually ends abruptly at a limit crisis around k = 0.17. Beyond the second crisis we observe in
W. Schnabl et al. /Full characterization of a strange attractor
84
Table 2 Chaotic windows in the essentially periodic in the range 0.1684 < /1 < 0.2085).“’ Range
of p values
begin 0.1687
Range
end
begin
0.1690
0.1981
regime
end 0.1993 0.2029 0.2041 0.2062
0.1807 0.1828 0.1891 0.1954
“‘In the case of very value is given.
2
of p values
0.1714 0.1801 0.1822 0.1873 0.1945
(period
0.2081
0.2074
narrow
windows
only one
parameter
essence a double periodic regime interrupted by small windows with significantly positive Lyapunov exponents. The chaotic windows at parameter values p < 0.16 are summarized in table 2. Approaching the chaotic regime in opposite direction, i.e. from high to low values of p, we observe a periodic attractor appearing first at p = 0.21. This attractor changes - at some points discontinuously - and small chaotic windows occur in the periodic regime. Eventually chaos is reached in one step at p = 0.17. For few selected values of the parameter p we constructed return maps on the unit interval which are shown in fig. 10. As in fig. 8 the downward Poincare cross-sections were mapped on a straight line. The 1.0
.a
.6
.4
.2
!I 0
.2
.4
.6
.8
1.0
0
.2
.4
.6
.a
1.0
0
.2
.4
.6
.a
1.0
0
.2
.4
.6
.6
1.0
0
.2
.4
.6
.B
1.0
0
.2
.4
.6
.B
1.0
1.0
.6
.4
.2
on the unit interval corresponding to the Poincare cross-section for different parameter values Fig. 10. Return maps f(,,, (,+,) - 0.029, -0.002, 0.016, 0.088 and 0.142. The points 5 = 0 and 5 = 1 correspond to the innermost and v = 0 and p = -0.101, outermost points of the (downward) PoincarC sets defined in fig. 7, respectively.
W. Schnabl et al. /Full characterization of a strange attractor
85
innermost points of the Poincare sets were defined as 5 = 0, the outermost points were chosen to represent the values 5 = 1. As expected the maps resemble closely the logistic map in the parameter range where we observe the Feigenbaum sequence. At higher parameter values (EL> 0) the map becomes more sophisticated, and resembles a cubic map. At still higher values of p it looks like a quartic map.
6. Mutation-caused
changes in the chaotic regime
Mutation is an unavoidable by-product of replication. Quantitative answers to the questions if and how complex dynamics is changing under the influence of mutation are of obvious importance for biological applications. Therefore we performed a study of the influence of mutation on the chaotic attractor in the replicator equation. Eq. (2) with rate constants Aij according to (9) was studied by a version of perturbation theory [44] under the simplifying assumption of equal mutation rates
+Aqx,e)
(27)
with L.& = xk
c Akixj
i
-
c Aijxixj
i,i
and
-kk =
C (Q/c,(e) AijXiXj - Qik(c> AkjXkXj). i,i
The first part-the unperturbed equation-is identical with eq. (4). Mutation A is considered as perturbation with E being the perturbation parameter. The replicator part is insensitive to the addition of constants to the reaction rate parameters: Aij -+Aij + A [45]. In order to make the perturbation approach applicable to eq. (9) we had to apply a sufficiently large value of this additive constant which was chosen to be A = 8. The influence of mutation on the position of the kth fixed point of the unperturbed replicator equation (4) (~~~‘1is expressed by X/JE) =x$+Edk+
@(E’).
(28)
As shown in ref. [44] the shift d, is obtained from
J( x~~‘)d, = -wk’(xi”), where J(x) represents the Jacobian of the unperturbed system: Jij = aSi/axj. In table 3 we summarize the shifts of all important fixed points of our model replicator equation determined by (9). Fig. 11 summarizes the changes in the chaotic attractor due to mutation. The main result of our study is that-at least this particular type of -chaotic dynamics is very sensitive to mutation: mutation
86
W. Schnabl et al. /Full characterization of a strange attractor Table 3 Shifts if individual fixed points caused by mutation assumed to be equal (Q,, = E Vi fj) and the results ka’
Fixed points
in the replicator model (9). All mutation rate were were obtained by means of perturbation theory [44]. d,
Shift vector
x~“’
16
1 ll(lOK
+ 17)
-2b + 1) 3j.L + 2
2
-(p+l) 1 0 0 1 00
3
’
3(18/.~ + 13) 89 176(3/~ + 4)
12
210~ + 181 - 88(3~
+ 4)
132
I
( - (9Oj.l + 73) 95 42(5/.~ + 6)
13
14
1 2(5/.~ + 9)
15(5/.~ + 6)
I
15/L + 88
\
-105
810~ + 1457 2(5/~ + 3)(5~
+ 9)(5/.~ + ll)(lO/.~
(5~ + ll)(
I
\
100/.~* + 265~ + 167)
(S/L + 9)(5/L + 11)(10/L
X
+ 17)
(5/k + 3)(5P
+9)(10/L
+ 17) + 17)
\- 1000/~~ - 4875~~ - 7730~ - 3979,
43
123
124
112(3/~ + 5)
651~ + 779
1 8(5~ + 6)
8(~ + 1)(3~ ’ X
+ 4)(3~
+ 5)(5fi
+ 3)(5~
+ 6)
(~+1)(5P+7)(30~z+25~-31)
’
150~’ + 395~~ - 55~’ - 767~ - 459 4(~ + 1)(3~
+ 4)(51.~ + 3)(5*
+ 6)
\ - 4( 150~~ + 530~~ + 544~~ + 65~ - 97)
W. Schnabl et al. /Full characterization of a strange attractor
87
Table 3. Continued Shift vector d,
Fixed points XI”’
ka’
134
2(5:+ 9)
160~ + 287 2(3~ + 5)(5j~ + 6)(5p + 9)(5/~ + 11)
( X
(5~ + ll)( loo/~’ + 265~ + 167) (5~ + 6)(5~ + 9)(5~ + 11) 1875~’ + 13125~~ + 33825~’ + 37815~ + 15392
\ - 1875cL4- 13750~’ - 36900$ - 42660~ - 17823 “The individual fixed points are denoted by indicating the non-vanishing coordinates as subscript. The tixed point I$‘) is non-hyperbolic - its Jacobian has one zero eigenvalue - and hence cannot be subjected to the analysis presented here. All other fixed points not included in the table lie outside the simplex S,.
105,
\
I
90-
765-
Fig. 11. Extension of the chaotic regime into the range of small mutation rates. The dynamics of the ACT model (V = 0) in the range -0.25. < r.~< 0.5 and 0 5 c < 1O-4 is shown. Typical shapes of the attractors are sketched. In the vertically hatched range trajectories converge to a stable point on the [1,4] edge. Periodic attractors are characterized by period multiplicities (1 = simple limit cycle). Chaotic dynamics is observed in the zone with slanted hatching. Note, that a Feigenbaum-sequence-like path with successive period doublings (1 + 2 + 4 + 8 + . .) is observed in the approach towards chaos from high to low mutation rates E.
88
W. Schnabl et al. /Full
characterization
of a strange attractor
6
106.E 015
020
iJ.025
Fig. 12. An enlargement of the bifurcation diagram near the limit crises of the unperturbed replicator equation (E = 0, F = 0.20). As in fig. 11 the numbers denote period multiplicities of the attractors vertical hatching indicates convergence to the boundary and slanted hatching refers to the chaotic regime.
probabilities as small as 2 X lo-” are sufficient to destroy the strange attractor. At higher mutation rates - E larger than some lop3 -the interior fixed point is stable. A decrease in the mutation rate leads to a sequence of period doublings which eventually ends up in chaos. A highly complex bifurcation structure arises on the introduction of mutation near the limit crisis of the unperturbed attractor. Fig. 12 shows an enlarged portion of fig. 11 with a resolution of approximately 0.005 in p and 5 x lo-’ in E. At the present stage of resolution we cannot say whether the bifurcation set forms a fractal or not. It is not surprising that even small mutation rates destroy the chaotic attractor since the mutation term causes the saddle focus P,,, to move into the unphysical range outside the simplex S, [44]. Therefore the trajectories can no longer come sufficiently close to the saddle and Silnikov’s mechanism breaks down. Furthermore the saddle point P,,, moves into the interior of S, and thus the limit crisis occurs earlier.
7. Conclusions The strange attractor in Lotka-Volterra and replicator systems of low dimension was fully characterized in a two-dimensional subspace of the twelve-dimensional parameter space. We were able to show that both chaotic attractors reported for these dynamical systems so far, lie within the same chaotic regime. There is no doubt that this regime extends also into other directions in parameter space. Because of technical limitations we have to leave the systematic multi-dimensional search as a project for the future. Still unanswered is also the question whether or not other chaotic regimes exist in areas which are not connected to the investigated one. One cross-section through the chaotic regime was studied in great detail. Along this path chaos developes in an almost perfect Feigenbaum scenario and then it breaks down in several interior crises - a major one and a few minor events - and ends in the final limit crisis. On further continuation along the path one encounters a periodic regime which is intermitted by a few, very small chaotic windows. The chaos is never fully developed: a small hole around the central - unstable -fixed point remains always
W. Schnabl et al. /Full characterization of a strange attractor
89
left out by the long-time trajectories. Right in the centre of the chaotic regime the fractal set produced by the attractor is of multi-fractal type. The distribution of dimensions covers a range from 1.3 to 3. The replicator equation allows straight inclusion of mutation, which inevitably occurs in nature. The most relevant result we obtained in the present study is a substantial reduction in the dynamical complexity with increasing mutation rates. Strange attractors are replaced by stable closed orbits, limit cycles are converted into stable attractors. It is interesting to note that we found a Feigenbaum sequence in opposite direction: reduction of the error-rate yields a cascade of bifurcations which eventually ends in the chaotic attractor of the unperturbed system. We cannot claim - nor do we intend to -that increased mutation will always leads to simpler dynamics but general intuition favours this view: mutation can be understood as a force pushing trajectories towards the centre of the concentration simplex. Wide excursions of trajectories which are characteristic for sophisticated closed orbits and strange attractors are less likely if the attractor lives near the centre.
Acknowledgements
Financial support was provided by the Austrian Fonds zur Forderung der wissenschaftlichen Forschung (Project No. 6864) and by the Stiftung Volkswagenwerk (BRD). Stimulating discussions with Professor Harald Posch and Professor Paul Phillipson are gratefully acknowledged. This work required extensive use of computers and we are particularly grateful for generous supply with CPU time at the IBM 3090-400 VF mainframe computer of the EDV-Zentrum der Universitat Wien within the EASI project of IBM and at the NAS 9160 mainframe computer of the IEZ Wien.
References 111 P.E. Kloeden and A.I. Mees, Bull. Math. Biol. 47 (1985) 697. [2] A. Chabra and R.V. Jensen, Phys. Rev. Lett. 62 (1989) 1327. [3] A.M. Fraser, Physica D 34 (1989) 391. 141 E. Gutierrez and H. Almirall, Bull. Math. Biol. 51 (1989) 785. 151 H.G. Schuster, Deterministic Chaos -An Introduction, 2nd Ed. (VCH, Weinheim, 1988). [61 R.W. Leven, B.-P. Koch and B. Pompe, Chaos in Dissipativen Systemen (Vieweg, Braunschweig, 1989). [7] M. Eigen, Naturwissenschaften 58 (1971) 465. [8] M. Eigen and P. Schuster, The Hypercycle - A Principle of Natural Self-organization (Springer, Berlin, 1979). [9] P. Schuster and K. Sigmund, J. Theor. Biol. 100 (1983) 533. [lo] P. Schuster, Physica D 22 (1986) 100. Ill] J. Hofbauer and K. Sigmund, The Theory of Evolution and Dynamical Systems- Mathematical Aspects of Selection (Cambridge Univ. Press, Cambridge, 1988). 1121 P. Schuster, Physica Scripta 35 (1987) 402. [13] J. Hofbauer, Nonlin. Anal. Theor. Meth. Appl. 5 (1981) 1003. [14] S. Smale, J. Math. Biol. 3 (1976) 5. [15] A. Arneodo, P. Coullet and C. Tresser, Phys. Lett. A 79 (1980) 259. [16] A. Arneodo, P. Coullet, J. Peyraud and C. Tresser, J. Math. Biol. 14 (1982) 153. [17] R.R. Vance, Am. Natur. 112 (1978) 797. [18] M.E. Gilpin, Am. Natur. 113 (1979) 306. 1191 J.F.C. Kingman, Mathematics of Genetic Diversity, CBMS NSF Regional Conf. Ser. Appl. Math., Vol. 34 (SIAM, Philadelphia, 1980). [20] J.F.C. Kingman, J. Appl. Prob. 15 (1978) 1. [21] P. Schuster, K. Sigmund and R. Wolff, J. Math. Anal. Appl. 78 (1980) 88.
90
W. Schnabl et al. /Full
characterization of a strange attractor
[22] P. Schuster, P. Stadler, W. Schnabl and C. Forst, Dynamics of Nonlinear Autocatalytic Reaction Networks (Akademie-Verlag, Berlin, 1990). [23] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields (Springer, Berlin, 1983) pp. 318-325. 1241 N. Samardzija and L.D. Greller, Bull. Math. Biol. 50 (1988) 465. [25] C. Grebogi, E. Ott and J.A. Yorke, Physica D 7 (1983) 181. [26] I.M. Bomze, Biol. Cybern. 48 (1983) 201. [27] P.F. Stadler and P. Schuster, Bull. Math. Biol. 52 (1990) 485. [28] S. Kim, S. Ostlund and G. Yu, Physica D 31 (1988) 117. [29] J.D. Farmer, E. Ott and J.A. Yorke, Physica D 7 (1983) 153. [30] G. Bennetin, L. Galgani, A. Giorgilli and J.-M. Strelcyn, Meccanica 15 (1980) 9. [31] A. Wolf, J.B. Swift, H.L. Swinney and J.A. Vastano, Physica D 16 (1985) 285. [32] P. Schuster, K. Sigmund and R. Wolff, J. Diff. Eqs. 32 (1979) 357. 1331 J. Hofbauer, P. Schuster and K. Sigmund, J. Math. Biol. 11 (1981) 155. [34] J.L. McCauley, Phys. Rep. 189 (1990) 225. 1351 P. Grassberger, Phys. I.&t. A 107 (1985) 101. [36] P. Grassberger, Phys. Lett. A 97 (1983) 227. [37] H.G.E. Hentschel and I. Procaccia, Physica D 8 (1983) 453. [38] P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50 (1983) 346. [39] P. Grassberger and I. Procaccia, Physica D 9 (1983) 189. [40] H. Atmanspacher, H. Scheingraber and W. Voges, Phys. Rev. A 37 (1988) 1314. [41] N.B.Abraham, A.M. Albano, A. Passamante and P.E. Rapp, eds., Measures of Complexity and Chaos (Plenum Press, New York, 1989). [42] T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia and B.I. Shraiman, Phys. Rev. A 33 (1986) 1141. [43] M.J. Feigenbaum, J. Stat. Phys. 46 (1987) 919. 1441 P.F. Stadler and P. Schuster, J. Math. Biol. submitted for publication. [45] J. Hofbauer, P. Schuster, K. Sigmund and R. Wolff, SIAM J. Appl. Math. 38 (1980) 282.
FULL CHARACTERIZATION
OF A STRANGE A’ITRACTOR
CHAOTIC DYNAMICS IN LOW-DIMENSIONAL REPLICATOR SYSTEMS
Wolfgang SCHNABL, Peter F. STADLER, Christian FORST and Peter SCHUSTER Institut fiir Theoretische Chemie der Uniuersitiit Wien, WiihringerstralJe 17, A-1090 Wien, Austria Received 17 January 1990 Revised manuscript received 8 September 1990 Accepted 8 September 1990 Communicated by A.T. Winfree
Two chaotic attractors observed in Lotka-Volterra equations of dimension n = 3 are shown to represent two different cross-sections of one and the same chaotic regime. The strange attractor is studied in the equivalent four-dimensional catalytic replicator network. Analytical expressions are derived for the Lyapunov exponents of the flow. In the centre of the chaotic regime the strange attractor was characterized by the Lyapunov dimension (D, = 2.06 f 0.02) and the Rtnyi fractal dimensions. The set corresponding to the attractor represents a multifractul. Its singularity spectrum is computed. One route in parameter space leading into the chaotic regime and crossing it was studied in detail. It leads to a Feigenbaum-type cascade of bifurcations. Before the chaos is fully developed the dynamic system passes an internal crisis through a period-two tangent bifurcation. Then a second sequence of period-doubling bifurcations leading to another chaotic regime is observed, which eventually ends in a crisis and then after a rather complicated periodic regime the attractor finally disappears. Trajectories in the intermediate periodic regime show transient chaotic bursts of intermittency type. A series of one-dimensional maps is derived from a properly chosen Poincare cross-section which illustrates structural changes in the attractor. Mutations are included in the catalytic replicator network and the changes in the dynamics observed are compared with the predictions of an approach based on perturbation theory. The most striking result is the gradual disappearance of complex dynamics with increasing mutation rates. The parameter space spanned by the mutation rate and selected parameters boundaries
of the replication network is split in regions of complex periodic behaviour near the chaotic regime. of these regions are likely to occur but the data available do not allow definitive conclusions yet.
Fractal
1. Introduction The search for chaotic dynamics became an issue in many scientific disciplines from physics to biology. The relevance of strange attractors in modeling real systems and the problem how to trace chaos in the experimental data were important questions already five years ago [l] and remained in the centre of interest since then [2-41. Numerical evidence for strange attractors was found in many model equations but unlike simple attractors they cannot be characterized by a few coordinates or coordinates and frequencies. In the pioneering works of several groups efficient techniques to analyse chaotic dynamics and fractal sets were developed. (For a recent review of the available repertoire of methods see refs. [5, 61.) In order to provide a sufficiently detailed basis for the comparison of model calculations and experimental data full characterization of strange attractors is required. We present such a study on a chaotic attractor which appears in the solutions of a differential equation which is of interest in biophysics and in theoretical ecology as well as in other disciplines. Autocatalytic reaction networks which involve two classes of catalytic processes - autocatalytic instruction for reproduction as well as material- and process-specific catalysis- were postulated as models for 0167-2789/91/$03.50
0 1991-
Elsevier
Science
Publishers
B.V. (North-Holland)
66
W. Schnabl et al. /Full characterization of a strange attractor
studies on prebiotic and early biological evolution scenarios [7, 81. Now they serve also as simple models showing the characteristic features of nonlinear dynamics in very different fields like molecular biology, population genetics, theoretical ecology or dynamical game theory [9]. Qualitative analysis and systematic numerical investigations revealed a very rich dynamics of the corresponding kinetic differential equations [lo]. The name replicator equations was coined for these systems because they reveal several unique - in the standardized form to be discussed in section 2 - fairly easy to prove features [ll]. The four-species replicator system shows chaotic dynamics for certain parameter values [121. According to Hofbauer [13] this autocatalytic network with four species is - apart from a transformation of the time scale - equivalent to a three-species Lotka-Volterra equation. Possible existence of very complex dynamical behaviour in Lotka-Volterra models was predicted by Smale [14]. Indeed, two different strange attractors were reported for the Lotka-Volterra system [15-181, which describes communities of one predator species living on two prey species. In this contribution we shall show that the two strange attractors lie on different cross-sections through one and the same chaotic regime. The number of degrees of freedom is minima1 for the occurrence of strange attractors and hence the chaotic regime in question houses the simplest strange attractors possible in Lotka-Volterra and replicator equations - replicator systems have normalized variables and hence the number of independent degrees of freedom is the number of species but one.
2. The replicator model Replicator equations are based on the simplest possible molecular mechanism of catalysed replication and mutation. The basic reactions are of the form (A) +Ij+Ii%
I,+I,+Ij,
i,j,k=l,...,
n.
(1)
Herein Ij is the template which is replicated and Ii the catalyst, (A) represents an overall notation for substrates needed for replication. It is assumed that the substrates are present in excess- their concentrations are buffered and hence do not change during the course of reactions. There are n3 possible reactions producing the n different species I, (k = 1,. . . , n) of the autocatalytic reaction network. Clearly it is very hard-if not impossible -to handle such a large number of parameters. In order to be able to analyze the reaction network simplifying assumptions are inevitable. It is straightforward to treat catalysed replication and mutation as independent superpositions: mutation frequencies depend on template (Ii) and target (I,), but not on the catalyst (Ii). The remaining 2n2 rate parameters are properly understood then as elements of two n x n matrices: a replication rate matrix A k (AjJ and a mutation matrix Q A (Qkj). The rate constants of the individual reactions are products of a replication and a mutation factor. Qkj .Aji, expressing the following sequence of processes: Ii catalyses the replication of I j which yields I k as an error copy. Error-free replication of species Ii- described by eq. (1) with I, = Ij - is assumed to occur with frequency Q,,. A mutation from Ij to I, (k #j) occurs with frequency Qkj. Every copy has to be either correct or a mutant and thus Q is a stochastic matrix (Xi =, Qkj = 1 Vj = 1,. . . , n). Concentrations of individual species are denoted by ck = [Ik]. As variables we use relative-or normalized - concentrations, xk = c&j”= icj, which obviously fulfil X; = ixk = 1.
U! Schnabl et al. /Full
characterization
Mass action kinetics applied to the autocatalytic reaction-mutation form dxk -
dt
61
of a strange attractor
network (1) yields an ODE of the
n =f,=
QkiAijxixj-x,@,
c
n,
k = l,...,
i,j=l
where an automatically adjusted dilution flux
k Aijxixj
CD=
i,j=l
is introduced in order to keep the sum of relative concentrations meaningful domain of eq. (2) is the simplex
s, A
(x,,...,
normalized
to unity. The physically
txi=l .
xn)l xi20,
i
I
i=l
If AijrO, i,j= l,..., n then S, is a compact, foreward invariant set for the ODE (2). An intensity of mutation can be measured by the mean mutation rate E: 1 E = n(n - 1)
5 i~i
Qij*
(3)
i,j=l
Throughout
this paper we use a particularly simple type of the mutation matrix
Qij= 1 - (n - l)e,
if i=j
and
Qij=e,
if
i#j.
It represents a special case of the house of cards model of mutation [19, 201, frequently used in theoretical population genetics. The important special case of error-free replication is represented by a unit mutation matrix, Q A IQ,, = S,,) with 6,, being the Kronecker symbol. Then eq. (1) is turned into the well known replicator equation :
~ AkiXi i=l
~
AijXiXj
(4)
i,j=l
The linear replicator equation (4) for n is flow-equivalent to the n - 1 species Lotka-Volterra
,
k=l,...,
n-l.
The variables are related by means of the transformation yk=xk/x,,
fork=l,...,n-1.
equation:
(5) introduced by Hofbauer [13]:
(6)
68
W. Schnabl et al. /Full
3. The chaotic regime in parameter
characterization
of a strange attractor
space
Two distinct chaotic attractors were found in Lotka-Volterra equations for three species: Vance [17] discovered a “quasi-cyclic” trajectory in a one-predator two-prey model and Arneodo, Coullet and Tresser [15] found a one-parameter family of strange attractors (ACT attractor). We used Hofhauer’s transform to construct the replicator equations which are equivalent to the models mentioned above. Then a barycentric transform [211 was applied to Vance’s model in order to shift the interior fixed point into the centre of the simplex S, - the point c = f
e
bkizi
-
i=l
k
bijziz/
i,j=l
has an interior fixed point i = (i,, i,, . . . , 2,). It can be shown [221 that there exists a di’omorphism which relates this replicator equation to another one in the notation of eq. (4) with ~~~~~
and
f=i(l,l,...,
1).
I
Both the ACT family and the Vance attractor (V) correspond to nonrobust phase portraits because of zero or almost zero (see e.g. A,, or A,, in A,) off-diagonal elements in the reaction matrices A. In order to make the numerical computations easy to reproduce all parameter values were rounded to three digits. Then the Vance model has a replication matrix 0
A,=
0.537 -0.535 0.536
0.063 0 0.38 - 0.032
0 - 0.036 0 - 0.004
0.437 - 0.001 0.655 0
I
(7)
and the ACT attractors are found at
J 0 A ACT=
1.1 -0.5 \ 1.7+/L
0.5 0 1 -1-p
-0.1 - 0.6 0 -0.2
0.1 '
0 0. 0,
(8)
Note that we use a parameter p which corresponds to p - 1.5 in ref. [15]. In order to find out whether or not these two attractors are related in the twelve-dimensional parameter space we performed a systematic numerical survey in the two-dimensional subspace (p, VI which connects the two models:
!
0
A(P, v) =
1 .l - 0.563~ 0.035v - 0.5 1.7 +/_L- 1.164~
0.5 - 0.437v 01 - 0.62~ - 1 - /L + 0.968~
-0.1 + O.lu 0.6 0 + 0.564~ - 0.2 + 0.196~
(9)
W. Schnabl et al. /Full
69
characterization of a strange attractor
1.0
.a .6 .4
.2
0.
-.2
-.4
-.6
-.a
t -1.0
1
! -1.5
’
’
’
’
! -1.0
’
’
’
’
! -.5
’
’
’
’
!
’
’
’
0
’ .5
v
Fig. 1. A two-dimensional cross-section through the chaotic regime of eq. (4). The cross-section is defined by the two parameters w and Y in the replication matrix A (9). The density of hatching refers to the dynamics of the attractor. In the high-density zone chaotic dynamics is observed. The lowest density indicates simple limit cycles, which undergo sequences of period doublings in the intermediate region. Numbers 1 to 14 indicate regions of qualitatively different dynamics on the boundary of the simplex S,. The corresponding phase portraits are shown in fig. 2.
We note that all three matrices of replication constants fulfil E;= ,Aij = 0.5 - a property which will be used together with the existence of a fixed point in the centre of the simplex later on. In the (IL, v) plane the ACT model lies on the line - 0.2 I p 5 0.05 and v = 0. The Vance attractor is close to the point (0,l). Both models are embedded in a Z-shaped two-dimensional manifold of strange attractors (fig. 1). In scanning through the manifold we could not find any sudden changes in the shape of the attractor. The mean revolution time also varies in apparently continuous manner all over the chaotic regime. A common feature of all phase portraits exhibiting chaotic behaviour (attractive or only transient) is a saddle focus Pi,, on the [l, 2,3] plane, which represents the stable manifold. The unstable manifold thus points into the interior of the simplex. We found that Silnikov’s condition [23] is always satisfied: the real positive eigenvalue is larger in magnitude than the real parts of the two complex conjugate eigenvalues. Although there is no homoclinic orbit possible for Plz3, a close relation to !%lnikov’s condition is evident, as was pointed out by Arneodo et al. [161: there is a heteroclinic orbit involving three saddles instead of
70
W. Schnabl et al. /Full
characterization
of a strange attractor
one and this orbit replaces the homoclinic one. For generalized Lotka-Volterra models with additional cubic terms Samardzija and Greller [24] found a different type of chaos arising from a fractal torus. For small values of p the fixed point in the centre of S, (c) is stable. Chaos is reached via a Hopf bifurcation and a subsequent cascade of period-doubling bifurcations. At first the limit cycle, and then the chaotic attractor, both grow steadily with increasing p values. At a parameter value around p = 0.21 the chaotic attractor collides with the stable manifold of a saddle on the [l, 2,41 plane and disappears suddenly. This behaviour is known as limit crisis [251: at the bifurcation point we find a heteroclinic orbit connecting the [l, 41 saddle, a saddle on the [l, 21 edge, and the saddle focus in the [l, 2,3] plane (see
A
a
Fig. 2. A sketch of fourteen flows on the boundaries of the simplex S, which can coexist with the chaotic attractor The numbers 1 to 14 are used as markers in fig. 1 to indicate where the phase portraits occur in parameter space.
in the interior.
71
W. Schnabl et al. /Full characterization of a strange attractor
fig. 7 below). The evidence for the orbit is derived from numerical integration. We cannot provide an analytical proof for the saddle connection between the saddle focus and Pi,. The saddle connections on the boundary follow from the invariance and the subsimplices and the classification of the flow of three-component replicator equations [26, 271. Beyond the limit crisis we find a field of transient chaos: the trajectory seems to draw an image of the chaotic attractor which in parameter space lies on the other side of the limit crisis before it ends up in a stable fixed point. For small values of v the frequency of the oscillations increases with decreasing parameter values-revolution times become shorter and shorter. Eventually the attractor collides with the stable manifold of the saddle on the [1,2] edge. The invariant hyperplane Ci,ixk = 1 is a repellor in this area and therefore integration is extremely sensitive to numerical errors. The chaotic attractor in the interior of the simplex S, coexists with a rather large variety of robust flows on the surface. Fig. 2 shows fourteen boundary flows consistent with chaos on the (II, v) plane. Fig. 3 sketches the phase portraits on the ACT and the Vance model together with a typical chaotic trajectory for each system. Direct evidence from the inspection of computed trajectories was supplemented by Fourier transformation [28]. Some typical examples of Fourier spectra in periodic and chaotic regimes of our attractor are shown in fig. 4. They are used as a sensitive indicator of the occurrence of chaos. In the four spectra chosen we observe a complete transition from a periodic window into the chaotic regime caused by a small change in the parameters CAP < 0.02).
4. Measures of the chaotic attractors The most common measure for a chaotic attractor is the spectrum of Lyapunov exponents [291. Let J(t) = dF(x(x(O), t)) be the Jacobian of the vector field 9 along a trajectory with starting point x(O).
x3
x2
A Fig, 3. Phase portraits and chaotic ACT model (p = 0, v = 0).
trajectories in the Vance
x3
x2
0 (7) and the ACT model (8). (A) Vance
model (CL= 0, v = 1) and (B)
W. Schnabl et al. /Full
characterization
of a strange attractor
I
u=-0.0536
8
6
I
4
2
0
-2 ~ 01
.02
.03
.04
.05
.06
.a7
.oa
.09
.I0
Frequency
A log[Al u=-0.0529
.Ol
.02
.03
.04
.05 6
.06
.07
.oa
.09
Frequency
Fig. 4. Fourier spectra of selected trajectories of the replicator model (9) with the parameter values v = 0 and (A) p = -0.0536, period-3 window, (B) -0.0529 period-3 window after period doubling, (C) -0.05 at the borderline of chaos and (D) -0.04, well developed chaos.
73
W. Schnabl et al. /Full characterization of a strange attractor
logCAl
logCAl a
.01
.02
.03
.04
.05
D Fig. 4.
Continued
.06
.07
.OB
.09
Frequency
74
W. Schnabl et al. /Full
characterization
of a strange attractor
Now we define an operator -Y(t)
A’? exp/oj(r)dT).
The time ordering operator ? has to be introduced since Jacobian matrices evaluated at different times - J(T) and J(T’) - usually do not commute. Let pi(f) be the eigenvalues of 1(t), then the Lyapunov exponents are given by *_= 1
lnh-%(t)ll
lim
f--to=
t
(11)
*
In most cases it is not possible to calculate Lyapunov exponents for the flows of differential equations analytically. An appropriate numerical algorithm has been given by Bennetin et al. 1301, a well tested FORTRAN program is available in the literature [31]. In the four-species replicator system there are three physical dimensions imbedded in a fourdimensional state space. The hyperplane
is invariant. Accordingly d(t) has three eigenvectors gi lying in this plane and their components fulfil C,$) = 0. The fourth eigenvector is pointing out of the hyperplane. Thus we expect to find three physically meaningful Lyapunov exponents. If all three of them are negative the attractor is a sink. For a periodic, quasiperiodic or chaotic attractor at least one exponent has to be zero: this exponent corresponds to an eigenvector which points in the direction of the trajectory. In particular, we have for periodic orbits A, I A, < 0, A, = 0, for quasiperiodic orbits A, < 0, A, = A, = 0 and finally for a chaotic attractor one Lyapunov exponent is negative, one is zero and one is positive, A, < 0, A, = 0, A, > 0. The sum of all Lyapunov exponents fulfils the equation
(12) For error-free divF=
replication - if the mutation matrix Q is the unit matrix - the divergence div F reduces to e
Akjxj - (n + 2)
k,j=l
k
Aijxixj.
(13)
i,j=l
The mean flux 5 is defined by Aijx,(
T) xj( T) dr.
If the replicator equation is permanent
[ll, 32, 33]- or at least if the trajectory remains in the interior of
75
W. Schnabl et al. /Full characterization of a strange attractor
the simplex-we
have
Since the last line is the equation fulfilled by the coordinates Xj = lim f/‘xj(r) t-m
0
of the interior fixed point we have
dr =Zj,
where f = (a,, E,, . . . , i,)
(14)
is the interior fixed point if it exists. Its coordinates
fulfil the equations
Vi=l,...,n
iAijaj=s,
and hence we have
h
Aijij
=
ns.
i,j=l
Accordingly we find lim :itdivY(_r(r))dr=
t+m
-2s.
(15)
Finally we are left with the following relation derived from eq. (12) for the sum
kAi=--$, $ Aij. i=l
( 124
l,J=l
Since we applied the burycentric transformation we obtain
s=;$
to our model, .i?j= l/n
holds for all coordinates
and
Aij.
t,J-1
Using the numerical values of (9) the sum of the four Lyapunov exponents becomes Cf= rhi = - i. The eigenvector pointing in the direction (1, 1, 1,l) is denoted by
1 e=2:
1 1 1 . 1
il
76
W. Schnabl et al. /Full
The corresponding Lyapunov the linearized model starting
we obtain
exponent (A,) can be calculated analytically. Let us consider from x(0) = e, which to first order in time is given by
the orbit of
@(T*).
X(T) =e+TdFe+ Then
characterization of a strange attractor
for the Lyapunov
exponent
The last equation is derived straightway by explicit evaluation of the integrand as the sum over all elements of the Jacobian. One obtains the same result -Q(T) for every time T and hence A, = -5. Thus we have in case of our model:
;hi= -a
and
A,=
-i.
(16)
i=l In the chaotic regime it is sufficient exponents are obtained from A,=O,
A,=
--i-A,
and
In fig. 5 we show the largest parameter CL. From the Lyapunov
.,=2-2=2(
therefore
A,=
ll++h’n:‘)
the largest
Lyapunov
exponent
A,. All other
-i.
Lyapunov
exponents
to compute
exponent
we can calculate
(17) (A,) of the strange an estimate
attractor
for the fractal
as a function
dimension
of the
of our attractor:
ifA,>O,
=l
if A, =O,
=o
if A, < 0.
(18)
From fig. 5 we obtain A, = 8 x lo-” for the highly developed chaos at p = Y = 0 and compute a dimension of D, = 2.06 with an estimated numerical uncertainty of kO.02 at this point. In order to study the fractal nature of the attractor quantitatively we computed also generalized or RCnyi dimensions, D,, which correspond to the scaling exponents of the qth moments of the attractor’s measure [34-371. The attractor is covered by M(1) cells with edges of length 1. By p, we denote the are probability to find a point of the attractor in cell number i (i = 1,. . . , M(l)). Then Rtnyi dimensions obtained from
log cT”=‘:ap=
D, = lim 140 (q-lllogl
lim
log P(I)
1+0 (4 - 1) log 1’
(19)
W. Schnabl et al. /Full characterization of a strange attractor
-.2
Fig. 5. The largest through the chaotic
-.15
-.05
-.l
Lyapunov exponent A, of the vector regime we chose v = 0.
0
.os
flow defined
.l
by eq. (9). In order
77
.15
.2
to show a representative
CI
cross-section
where we introduced the correlation integral defined by Cq(l) G Cfl’?pp. In order to circumvent an indeterminant expression in the denominator the RCnyi dimension D, is obtained by taking the limit lim q ~ 10, = D,. For the actual calculations we rewrite the sum over the pi’s in terms of the probability ji of the trajectory to cross a box of size I around the iterate x(tj). We make use of the equation
(20) and find that the latter probability is given by [38, 391:
pj = $+rfm
k:
o( I - llXi - Xjll))
(21)
i=l
where 0 is used for the Heaviside function: @(x-x,)=0
for
x
and
0(x -xt,)
= 1 for x2x,.
Thus we are left with the easier to compute expression
(22)
W. Schnabl et al. /Full
78
for the correlation
characterization of a strange attractor
integral [40]. In the limit q -+ 1 we find the information
dimension
(23) The RCnyi dimensions for q = 0 are equivalent to the Hausdorff dimensions, for q = 2 they correspond to the correlation dimensions. The evaluation of RCnyi dimensions turned out to be extremely time consuming and therefore we restricted computations to one point right in the centre of the chaotic regime (CL= v = 0). In order to calculate the correlation integrals P(1) the trajectory was decomposed into a time series X=(&=x&)
Vi=1,2,...,N)
(24)
with t, < t, c . . . -c t,. Some ten thousand points are required which have to be distributed independently along a trajectory. This condition can be met by collecting subsequent points which are separated by at least one tenth of a full revolution on the attractor. The data set used in our computation was obtained by sampling of 120 different trajectories with randomly chosen initial conditions. A transient phase of t, = 3000 was discarded and then 1768 points were recorded in equal time intervals of At = 20 for each trajectory. The mean revolution time on the attractor is of the order 70 < t, < 80. Correlation integrals were computed for individual trajectories and for the whole data set of about 200000 points. Superposition of about 5 trajectories was already sufficient to achieve deviations from the final result in the percent range. Convergence for q 2 - 2 was generally good. Inserting the fractal set obtained in the described way into eqs. (22) and (23) we found D, = 2.02,
D, = 1.91
and
D, = 1.77.
The estimated numerical uncertainty of our values is certainly less than &0.05. The data apparently are not in perfect agreement with the Kaplan-Yorke conjecture according to which the information dimension of the attractor, D,, should be equal to the Lyapunov dimension D,. The latter is about 15% larger than the former in our case. We cannot claim with ultimate certainty that this is not a result of some hidden numerical errors, but the discrepancy exceeds our estimated error bounds. The cause of such deviations may lie in the fact that higher-order terms in the box length 1 are non-negligible even at the smallest numerically accessible I-values. In addition, the Lyapunov dimension is unlikely to be too large since it has to fulfil D, > 2 by the nature of our attractor. For a recent discussion of systematic errors in computations of the correlation integral see several contributions in ref. [41]. Fractal sets can be characterized by the scaling properties of normalized distributions or measures lying upon them [42]. We briefly repeat the basic expressions and relations. Scaling properties are expressed in terms of singularity strengths (Y and a function f(o) which reflects the density of their distribution. The spectrum of singularities is given by the possible range of (Yvalues and the density f(a). In order to define these quantities we cover the attractor with boxes of size 1 and denote the number of points inside the ith box by Ni. Then the probability that this box contains a point of the attractor’s time series (24) is given by Pi(l) = limN_JNi/N). The 1ocal singularity strength (Y~is now given as the scaling exponent Pi(Z) a PI. The singularity strength cy, thus varies with the region of the attractor on which it is taken. Its probability density - the probability to encounter a value of (Ybetween (Yand 6 + d& - is given
79
W. Schnabl et al. /Full characterization of a strange attractor
Prob{ & < (Y< & + d&l = da p( (Y) I-fC’), where f(a) is a continuous function. The set of singularity strengths (Y~itself is a fractal and f(a) can be understood as its - fractal - dimension. In other words, fractal measures are modeled by interwoven sets of singularities of strength (Y-each one having its own fractal dimension f(a). The singularity spectrum f(a) is related to the RCnyi dimensions D, since P(Z)
= C(q,Z)
=/dGp(&)
Iqail-f(@.
The dimensions are defined in the limit I + 0 and accordingly we may approximate the integral for nonzero densities p(G) by the term with the smallest exponent q& -f(h). This value a(q) is defined by
and the existence of a minimum requires
which finally yields f’(a(q)) = q and f”(a(q)) < 0. We have two relations:
D, =
&
[w(q)
and 4s) =
-f(hl))l
&[(4- l)Dq]
(25)
which allow to compute the singularity spectrum from known RCnyi dimensions or vice versa. It can be shown straightway that the Hausdorff dimension of the attractor fulfills D, =f(ama,), where (Y,, is the position of the maximum of f(a). Feigenbaum [43] recently pointed out that this formalism is related to equilibrium statistical mechanics. A new function T(q) A (q - l>D, is introduced. The singularity spectrum f(a) is the Legendre transform of T. In Feigenbaum’s formalism T(q) and q are the conjugate thermodynamic variables to f(a) and D4. The singularity spectrum of the strange attractor with I_L= Y = 0 was calculated from a series of RCnyi dimensions. For a given set of values q = qk the singularity strengths were obtained by partial differentiation: a(q) = i3~/d(u. The partial derivatives were computed by means of a five-point differentiation formula. Then we derived f(a) for discrete points (Ye= dq,) from f(%)
&fk =
akqk
-
(qk
-
l)D,k.
(26)
Limiting dimensions can be derived straightway from the singularity spectrum: D,, = cdq = +m) and D_, = a(q = -m). Both limits are related to the singularity spectrum by f(a(q = km)) = 0. In the case of our chaotic attractor we find 1.3 I (YI 3. This result can be interpreted as follows: In case the trajectory visits a volume element in the three-dimensional concentration simplex S,, it occupies a fractal set whose measure has a Hausdorff dimension of cx= 1.3 or larger. For negative values of q the estimates of the RCnyi dimensions are rather poor, because the limit in eq. (22) converges very slowly. In our case convergence was poor for q = - 1.5 and q = -2, computed results with still smaller q-values
80
U! Schnabl et al. /FUN characterization of a strange attractor
Fig. 6. The singularity spectrum of the strange attractor in the centre of the chaotic regime (p = 0, v = 0). Note that the maximum of f(a) represents the Hausdorff dimension CD,,) of the attractor. Computation of the points in the branch for positive q-values (left-hand branch) converged well. Numerical uncertainties due to incomplete convergence lie within the size of the full circles. Points at the branch for negative q-values (right-hand branch) are inaccurate because of poor convergence. The curve extrapolates towards the physically meaningful limit a = 3.
are completely unreliable. Thus the right-hand branches of the computed singularity spectra are inaccurate and we expect rather large errors for the extrapolated values of D_,. Indeed straight extrapolation of the points obtained would yield an upper limit larger than (Y= 3 (fig. 61. The correct upper limit, however, has to fulfil (Y2 3 since the trajectory is ultimately restricted to the simplex S,.
5. The routes to chaos Routes into chaotic dynamics were studied by means of Poincare maps in the plane defined by x4 = $. Distributions of points at which a long-time trajectory crosses this plane are shown in fig. 7. A Cartesian coordinate system is defined with the origin in the central fixed point c = i(l, 1, 1,ll (fig. 7) and the branch starting near c and extending towards the projection of the [l, 21-edge onto the plane x4 = f is chosen as Poincare cross-section. Its shape is very close to a straight line. Time series % = (&,, _E2,. . . , E’) are defined as successive points at which the trajectory passes the Poincart cross-section. They were recorded as a function of the parameter /L (- 0.21 < p < 0.21) with the second parameter held constant at the value of the ACT model (Y = 0). The corresponding plot is shown in fig. 8. The first part ( - 0.21 < p < - 0.01) reminds of the analogous plot of the discrete logistic map on the unit interval: X,, +, = rX,(l -X,1. Indeed there is a monotonic - almost linear - relation between the values of the constant r and the parameter p at which the characteristic period doublings and periodic windows in the chaotic regime occur (fig. 9). The sequence of events into chaos and the structure of the chaotic regime along this route show perfect correspondence to the well-known Feigenbaum sequence. A more precise characterization of the individual periodic regimes is contained in table 1.
W. Schnabl et al. /Full characterization of a strange attractor
Fig. 7. Points at which a long-time trajectory crosses (A) and three examples with v = 0 and /L = -0.1001
of the attractor the plane x 4 = $. The orlentatlon (B), 0.0001 (C) and 0.1 (D) are shown.
relative
to the cross-section
W. Schnabl et al. /Full characterization of ‘a strange attractor
82
-4
I
I! i \
0.
F’N’t
I ! ‘I’(““’
-.2
-.15
! -.l
I”‘!““!
-.05
““‘I’ 0
.05
1’ .l
!
I”“’
.15
.2
Fig. 8. A PoincarC map of the attractors with Y = 0 as a function of the parameter g. A projection x-axis defined in fig. 7 is shown. The origin (x = 0) coincides with the hxed point in the centre (c).
0.
P
of the Poincare
map onto the
-
-.02
-
-.04
-
-.06
-
-.08
-
-.lO
-
-.12
-
-.14/-‘,,, 3.4
I ,,,,,,,,,,,,,,‘,,,,,,,11, 3.5
3.6
~. 3.7
3.8
3.9
r 4.0
Fig. 9. Relation between the critical values of the parameter r in Feigenbaum’s sequence of bifurcations and periodic windows in the chaotic regime of the logistic map on the unit interval and the values of the parameter p at which the analogous events occur in the dynamics of the replicator equation.
W. Schnabl et al. /Full
Table 1 Periodic windows Periodb’
in the chaotic Range
regime (-0.105
83
of a strange attractor
< /.L < 0.1684).@ Comments
of p values end
begin 12 6 5 10 3 ff. 5 7 6 4 ff. 6 7
characterization
begin of Feigenbaum
-0.1034 - 0.0956 - 0.0770
sequence
- 0.0948 - 0.760 - 0.0632
- 0.0587
- 0.0527 - 0.0398 - 0.0359 - 0.0326
- 0.0255
- 0.0250 -0.0182 -0.0116 end of Feigenbaum
9
sequence
0.0082
3 ff. 4 6
0.0118
2 ff. 12 7
0.0544
12 8 ff. 4 6 12 6
0.0865 0.0883 0.0922 0.1057
2 ff. 6 10 4
0.1207
0.0127 0.0227 0.0448 interior
crisis
0.0772 0.0784 0.0844 0.0868 0.0892 0.0980 0.1060 0.1090 0.1102 0.1291 0.1312 0.1489 0.1543
a)In the case of very narrow windows only one parameter value is given. “Windows in which at least one period doubling event was identified are indicated
by “ff’.
As in the logistic map, the range on the PoincarC section which is visited by the chaotic trajectory grows with increasing values of the parameter p. In contrast to the former, however, a certain region around the central fixed point c is always left blank. Beyond the end of the Feigenbaum sequence we observe a broad continuation of the chaotic regime up to a sharp interior crisis at p = 0.045. Characteristic accumulation zones of points on the cross-section which continue as lines through periodic windows appear in the range above /..L= -0.01 as well. They look as if they were reflected at the lower boundary of the chaotic regime near the point c. In order to search for multiple attractors in the interior of S, several different initial values- all chosen in the plane x4 = f -were applied in the parameter range 0.05 < p < 0.17. No indication of a second attractor in the interior of S, was found. There are, however, regions in this plane which have the fixed points P,, or P3 as w-limits (see fig. 2B). At p values above the interior crisis we observe a second period doubling cascade ending in a small chaotic regime which is again interrupted by periodic windows. Chaos expands with increasing p-values and eventually ends abruptly at a limit crisis around k = 0.17. Beyond the second crisis we observe in
W. Schnabl et al. /Full characterization of a strange attractor
84
Table 2 Chaotic windows in the essentially periodic in the range 0.1684 < /1 < 0.2085).“’ Range
of p values
begin 0.1687
Range
end
begin
0.1690
0.1981
regime
end 0.1993 0.2029 0.2041 0.2062
0.1807 0.1828 0.1891 0.1954
“‘In the case of very value is given.
2
of p values
0.1714 0.1801 0.1822 0.1873 0.1945
(period
0.2081
0.2074
narrow
windows
only one
parameter
essence a double periodic regime interrupted by small windows with significantly positive Lyapunov exponents. The chaotic windows at parameter values p < 0.16 are summarized in table 2. Approaching the chaotic regime in opposite direction, i.e. from high to low values of p, we observe a periodic attractor appearing first at p = 0.21. This attractor changes - at some points discontinuously - and small chaotic windows occur in the periodic regime. Eventually chaos is reached in one step at p = 0.17. For few selected values of the parameter p we constructed return maps on the unit interval which are shown in fig. 10. As in fig. 8 the downward Poincare cross-sections were mapped on a straight line. The 1.0
.a
.6
.4
.2
!I 0
.2
.4
.6
.8
1.0
0
.2
.4
.6
.a
1.0
0
.2
.4
.6
.a
1.0
0
.2
.4
.6
.6
1.0
0
.2
.4
.6
.B
1.0
0
.2
.4
.6
.B
1.0
1.0
.6
.4
.2
on the unit interval corresponding to the Poincare cross-section for different parameter values Fig. 10. Return maps f(,,, (,+,) - 0.029, -0.002, 0.016, 0.088 and 0.142. The points 5 = 0 and 5 = 1 correspond to the innermost and v = 0 and p = -0.101, outermost points of the (downward) PoincarC sets defined in fig. 7, respectively.
W. Schnabl et al. /Full characterization of a strange attractor
85
innermost points of the Poincare sets were defined as 5 = 0, the outermost points were chosen to represent the values 5 = 1. As expected the maps resemble closely the logistic map in the parameter range where we observe the Feigenbaum sequence. At higher parameter values (EL> 0) the map becomes more sophisticated, and resembles a cubic map. At still higher values of p it looks like a quartic map.
6. Mutation-caused
changes in the chaotic regime
Mutation is an unavoidable by-product of replication. Quantitative answers to the questions if and how complex dynamics is changing under the influence of mutation are of obvious importance for biological applications. Therefore we performed a study of the influence of mutation on the chaotic attractor in the replicator equation. Eq. (2) with rate constants Aij according to (9) was studied by a version of perturbation theory [44] under the simplifying assumption of equal mutation rates
+Aqx,e)
(27)
with L.& = xk
c Akixj
i
-
c Aijxixj
i,i
and
-kk =
C (Q/c,(e) AijXiXj - Qik(c> AkjXkXj). i,i
The first part-the unperturbed equation-is identical with eq. (4). Mutation A is considered as perturbation with E being the perturbation parameter. The replicator part is insensitive to the addition of constants to the reaction rate parameters: Aij -+Aij + A [45]. In order to make the perturbation approach applicable to eq. (9) we had to apply a sufficiently large value of this additive constant which was chosen to be A = 8. The influence of mutation on the position of the kth fixed point of the unperturbed replicator equation (4) (~~~‘1is expressed by X/JE) =x$+Edk+
@(E’).
(28)
As shown in ref. [44] the shift d, is obtained from
J( x~~‘)d, = -wk’(xi”), where J(x) represents the Jacobian of the unperturbed system: Jij = aSi/axj. In table 3 we summarize the shifts of all important fixed points of our model replicator equation determined by (9). Fig. 11 summarizes the changes in the chaotic attractor due to mutation. The main result of our study is that-at least this particular type of -chaotic dynamics is very sensitive to mutation: mutation
86
W. Schnabl et al. /Full characterization of a strange attractor Table 3 Shifts if individual fixed points caused by mutation assumed to be equal (Q,, = E Vi fj) and the results ka’
Fixed points
in the replicator model (9). All mutation rate were were obtained by means of perturbation theory [44]. d,
Shift vector
x~“’
16
1 ll(lOK
+ 17)
-2b + 1) 3j.L + 2
2
-(p+l) 1 0 0 1 00
3
’
3(18/.~ + 13) 89 176(3/~ + 4)
12
210~ + 181 - 88(3~
+ 4)
132
I
( - (9Oj.l + 73) 95 42(5/.~ + 6)
13
14
1 2(5/.~ + 9)
15(5/.~ + 6)
I
15/L + 88
\
-105
810~ + 1457 2(5/~ + 3)(5~
+ 9)(5/.~ + ll)(lO/.~
(5~ + ll)(
I
\
100/.~* + 265~ + 167)
(S/L + 9)(5/L + 11)(10/L
X
+ 17)
(5/k + 3)(5P
+9)(10/L
+ 17) + 17)
\- 1000/~~ - 4875~~ - 7730~ - 3979,
43
123
124
112(3/~ + 5)
651~ + 779
1 8(5~ + 6)
8(~ + 1)(3~ ’ X
+ 4)(3~
+ 5)(5fi
+ 3)(5~
+ 6)
(~+1)(5P+7)(30~z+25~-31)
’
150~’ + 395~~ - 55~’ - 767~ - 459 4(~ + 1)(3~
+ 4)(51.~ + 3)(5*
+ 6)
\ - 4( 150~~ + 530~~ + 544~~ + 65~ - 97)
W. Schnabl et al. /Full characterization of a strange attractor
87
Table 3. Continued Shift vector d,
Fixed points XI”’
ka’
134
2(5:+ 9)
160~ + 287 2(3~ + 5)(5j~ + 6)(5p + 9)(5/~ + 11)
( X
(5~ + ll)( loo/~’ + 265~ + 167) (5~ + 6)(5~ + 9)(5~ + 11) 1875~’ + 13125~~ + 33825~’ + 37815~ + 15392
\ - 1875cL4- 13750~’ - 36900$ - 42660~ - 17823 “The individual fixed points are denoted by indicating the non-vanishing coordinates as subscript. The tixed point I$‘) is non-hyperbolic - its Jacobian has one zero eigenvalue - and hence cannot be subjected to the analysis presented here. All other fixed points not included in the table lie outside the simplex S,.
105,
\
I
90-
765-
Fig. 11. Extension of the chaotic regime into the range of small mutation rates. The dynamics of the ACT model (V = 0) in the range -0.25. < r.~< 0.5 and 0 5 c < 1O-4 is shown. Typical shapes of the attractors are sketched. In the vertically hatched range trajectories converge to a stable point on the [1,4] edge. Periodic attractors are characterized by period multiplicities (1 = simple limit cycle). Chaotic dynamics is observed in the zone with slanted hatching. Note, that a Feigenbaum-sequence-like path with successive period doublings (1 + 2 + 4 + 8 + . .) is observed in the approach towards chaos from high to low mutation rates E.
88
W. Schnabl et al. /Full
characterization
of a strange attractor
6
106.E 015
020
iJ.025
Fig. 12. An enlargement of the bifurcation diagram near the limit crises of the unperturbed replicator equation (E = 0, F = 0.20). As in fig. 11 the numbers denote period multiplicities of the attractors vertical hatching indicates convergence to the boundary and slanted hatching refers to the chaotic regime.
probabilities as small as 2 X lo-” are sufficient to destroy the strange attractor. At higher mutation rates - E larger than some lop3 -the interior fixed point is stable. A decrease in the mutation rate leads to a sequence of period doublings which eventually ends up in chaos. A highly complex bifurcation structure arises on the introduction of mutation near the limit crisis of the unperturbed attractor. Fig. 12 shows an enlarged portion of fig. 11 with a resolution of approximately 0.005 in p and 5 x lo-’ in E. At the present stage of resolution we cannot say whether the bifurcation set forms a fractal or not. It is not surprising that even small mutation rates destroy the chaotic attractor since the mutation term causes the saddle focus P,,, to move into the unphysical range outside the simplex S, [44]. Therefore the trajectories can no longer come sufficiently close to the saddle and Silnikov’s mechanism breaks down. Furthermore the saddle point P,,, moves into the interior of S, and thus the limit crisis occurs earlier.
7. Conclusions The strange attractor in Lotka-Volterra and replicator systems of low dimension was fully characterized in a two-dimensional subspace of the twelve-dimensional parameter space. We were able to show that both chaotic attractors reported for these dynamical systems so far, lie within the same chaotic regime. There is no doubt that this regime extends also into other directions in parameter space. Because of technical limitations we have to leave the systematic multi-dimensional search as a project for the future. Still unanswered is also the question whether or not other chaotic regimes exist in areas which are not connected to the investigated one. One cross-section through the chaotic regime was studied in great detail. Along this path chaos developes in an almost perfect Feigenbaum scenario and then it breaks down in several interior crises - a major one and a few minor events - and ends in the final limit crisis. On further continuation along the path one encounters a periodic regime which is intermitted by a few, very small chaotic windows. The chaos is never fully developed: a small hole around the central - unstable -fixed point remains always
W. Schnabl et al. /Full characterization of a strange attractor
89
left out by the long-time trajectories. Right in the centre of the chaotic regime the fractal set produced by the attractor is of multi-fractal type. The distribution of dimensions covers a range from 1.3 to 3. The replicator equation allows straight inclusion of mutation, which inevitably occurs in nature. The most relevant result we obtained in the present study is a substantial reduction in the dynamical complexity with increasing mutation rates. Strange attractors are replaced by stable closed orbits, limit cycles are converted into stable attractors. It is interesting to note that we found a Feigenbaum sequence in opposite direction: reduction of the error-rate yields a cascade of bifurcations which eventually ends in the chaotic attractor of the unperturbed system. We cannot claim - nor do we intend to -that increased mutation will always leads to simpler dynamics but general intuition favours this view: mutation can be understood as a force pushing trajectories towards the centre of the concentration simplex. Wide excursions of trajectories which are characteristic for sophisticated closed orbits and strange attractors are less likely if the attractor lives near the centre.
Acknowledgements
Financial support was provided by the Austrian Fonds zur Forderung der wissenschaftlichen Forschung (Project No. 6864) and by the Stiftung Volkswagenwerk (BRD). Stimulating discussions with Professor Harald Posch and Professor Paul Phillipson are gratefully acknowledged. This work required extensive use of computers and we are particularly grateful for generous supply with CPU time at the IBM 3090-400 VF mainframe computer of the EDV-Zentrum der Universitat Wien within the EASI project of IBM and at the NAS 9160 mainframe computer of the IEZ Wien.
References 111 P.E. Kloeden and A.I. Mees, Bull. Math. Biol. 47 (1985) 697. [2] A. Chabra and R.V. Jensen, Phys. Rev. Lett. 62 (1989) 1327. [3] A.M. Fraser, Physica D 34 (1989) 391. 141 E. Gutierrez and H. Almirall, Bull. Math. Biol. 51 (1989) 785. 151 H.G. Schuster, Deterministic Chaos -An Introduction, 2nd Ed. (VCH, Weinheim, 1988). [61 R.W. Leven, B.-P. Koch and B. Pompe, Chaos in Dissipativen Systemen (Vieweg, Braunschweig, 1989). [7] M. Eigen, Naturwissenschaften 58 (1971) 465. [8] M. Eigen and P. Schuster, The Hypercycle - A Principle of Natural Self-organization (Springer, Berlin, 1979). [9] P. Schuster and K. Sigmund, J. Theor. Biol. 100 (1983) 533. [lo] P. Schuster, Physica D 22 (1986) 100. Ill] J. Hofbauer and K. Sigmund, The Theory of Evolution and Dynamical Systems- Mathematical Aspects of Selection (Cambridge Univ. Press, Cambridge, 1988). 1121 P. Schuster, Physica Scripta 35 (1987) 402. [13] J. Hofbauer, Nonlin. Anal. Theor. Meth. Appl. 5 (1981) 1003. [14] S. Smale, J. Math. Biol. 3 (1976) 5. [15] A. Arneodo, P. Coullet and C. Tresser, Phys. Lett. A 79 (1980) 259. [16] A. Arneodo, P. Coullet, J. Peyraud and C. Tresser, J. Math. Biol. 14 (1982) 153. [17] R.R. Vance, Am. Natur. 112 (1978) 797. [18] M.E. Gilpin, Am. Natur. 113 (1979) 306. 1191 J.F.C. Kingman, Mathematics of Genetic Diversity, CBMS NSF Regional Conf. Ser. Appl. Math., Vol. 34 (SIAM, Philadelphia, 1980). [20] J.F.C. Kingman, J. Appl. Prob. 15 (1978) 1. [21] P. Schuster, K. Sigmund and R. Wolff, J. Math. Anal. Appl. 78 (1980) 88.
90
W. Schnabl et al. /Full
characterization of a strange attractor
[22] P. Schuster, P. Stadler, W. Schnabl and C. Forst, Dynamics of Nonlinear Autocatalytic Reaction Networks (Akademie-Verlag, Berlin, 1990). [23] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields (Springer, Berlin, 1983) pp. 318-325. 1241 N. Samardzija and L.D. Greller, Bull. Math. Biol. 50 (1988) 465. [25] C. Grebogi, E. Ott and J.A. Yorke, Physica D 7 (1983) 181. [26] I.M. Bomze, Biol. Cybern. 48 (1983) 201. [27] P.F. Stadler and P. Schuster, Bull. Math. Biol. 52 (1990) 485. [28] S. Kim, S. Ostlund and G. Yu, Physica D 31 (1988) 117. [29] J.D. Farmer, E. Ott and J.A. Yorke, Physica D 7 (1983) 153. [30] G. Bennetin, L. Galgani, A. Giorgilli and J.-M. Strelcyn, Meccanica 15 (1980) 9. [31] A. Wolf, J.B. Swift, H.L. Swinney and J.A. Vastano, Physica D 16 (1985) 285. [32] P. Schuster, K. Sigmund and R. Wolff, J. Diff. Eqs. 32 (1979) 357. 1331 J. Hofbauer, P. Schuster and K. Sigmund, J. Math. Biol. 11 (1981) 155. [34] J.L. McCauley, Phys. Rep. 189 (1990) 225. 1351 P. Grassberger, Phys. I.&t. A 107 (1985) 101. [36] P. Grassberger, Phys. Lett. A 97 (1983) 227. [37] H.G.E. Hentschel and I. Procaccia, Physica D 8 (1983) 453. [38] P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50 (1983) 346. [39] P. Grassberger and I. Procaccia, Physica D 9 (1983) 189. [40] H. Atmanspacher, H. Scheingraber and W. Voges, Phys. Rev. A 37 (1988) 1314. [41] N.B.Abraham, A.M. Albano, A. Passamante and P.E. Rapp, eds., Measures of Complexity and Chaos (Plenum Press, New York, 1989). [42] T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia and B.I. Shraiman, Phys. Rev. A 33 (1986) 1141. [43] M.J. Feigenbaum, J. Stat. Phys. 46 (1987) 919. 1441 P.F. Stadler and P. Schuster, J. Math. Biol. submitted for publication. [45] J. Hofbauer, P. Schuster, K. Sigmund and R. Wolff, SIAM J. Appl. Math. 38 (1980) 282.