B dletin ol Mathematical Biolob'v VoI. 53, No. 3, pp. 469 485, 1991. Printed in Great Britain
0092 8240/9153.00 + 0 Pergamon Press pie © 1991 Society for Mathematical Biology
SMALL AUTOCATALYTIC REACTION NETWORKSIII. MONOTONE GROWTH FUNCTIONS BARBEL M . R . STADLER AND PETER F . STADLER I n s t i t u t ffir theoretische Chemie der Universitfit Wien, WfihringerstraBe 17, A 1090 W i e n , A u s t r i a (E.mai~bitnet
:
[email protected] W1UNI11)
The classification of the dynamical behaviour of first order replicator equations is extended to models with monotonical growth rates. It is shown that for two species there is a general classification independent of the particular form of the growth function. For three species a common dynamical behaviour for all power laws can be found and the existence of limit cycles is disproved. F o r more general growth functions, however, limit cycles may occur.
1. Introduction. In the last decade it turned out that a single type of nonlinear ODE, in the beginning designed to describe the replication of macromolecules, is appropriate to model systems in very different fields of science, as in molecular biology, population ecology or evolution theory. Some examples are template induced replication of RNA (Cech, 1986; Doudna and Szostak, 1989), idiotypic networks in theoretical immunology (Perelson, 1988), game dynamics of MaynardSmith models of altruism in animal behaviour (MaynardSmith, 1982) and the hypercycle model ofprebiotic evolution (Eigen and Schuster, 1979). The kinetics of autocatalytic, i.e. RNA template instructed RNA synthesis has been shown to be an extremely complicated reaction network consisting of thousands of elementary steps (Biebricher et al., 1983, 1984, 1985). Nevertheless, under simple experimental conditionsexcess of monomer educts and polymerizing enzymethe overall rate of RNA synthesis is proportional to the concentration of the template (I~). Biebricher et al. (1984) showed that the experiments on in vitro replication of singlestranded RNA species agree well with the following MichaelisMententype model if the concentration [E] of the enzyme, Qflreplicase, is sufficiently small : d kl [E][So] dt [I] = [I+1 k2 +k3[So] ' 469
470
B . M . R . STADLER and P. F. STADLER
where [I+] and [I] denote the concentrations of complementary RNA strands and [So] is the effective substrate concentration, the harmonic mean of the concentrations of the mononucleotidetriphosphates. Clearly, the reaction rate is a nonlinear but monotone function of [So]. The replication of RNA phages in vitro consists of similar but much more complicated reactions (Weissmann, 1974; Biebricher and Eigen, 1988). Modelling RNA replication and catalysis inevitably leads to autocatalytic reaction networks. In general a reaction of the form (A) + I T + I c + 2 I T + I c + (B), where the same molecule may act as a catalyst, Ic, and as a template, Iv, gives raise to reaction rates of the form d[¢] d t  fk([I1], [ I 2 ] , . . . ,
which depend on the concentrations of all interacting species. The most simple attempt is based on mass action kinetics (Eigen, 1971 ; Schuster et al., 1980) ; the overall rate of RNA synthesis is therefore proportional to the concentrations of both the catalyst and the template : d[I~] d~ oc [IT] [ic]A cleavage reaction of an RNA molecule by a specific RNA catalyst follows the same rate law with a negative rate constant (Cech, 1986). Recent attempts to model the dynamics of the immune system (Perelson, 1988) involve differential equations of the form d[B ] dt  f([&]'
 fl,[&]
where [Bk] denotes the concentration of B cells of type k, [B] is the total concentration of B cells. The constants ak and flk involve the fraction of stimulated ktype B cells and their degradation rates. The function f is assumed to be very rapidly decreasing in [B] and monotonously increasing (self stimulation) or decreasing (self inhibition) in [Bk]. Some 20 years ago models were designed to explain the origin of altruistic behaviour in conflicts between animals of the same species (MaynardSmith, 1982). Similar equations have been used to investigate the relation between evolutionary stable strategies (ESSs) and the approach towards stable rest points of more complicated attractors of the underlying dynamical systems. (Taylor and Yonker, 1978; Hofbauer et al., 1979; Zeeman, 1981). These
SMALL AUTOCATALYTIC REACTION NETWORKSIII
471
game dynamical models yield an ODE of the form xk = Xk([AX]k  x A x ) .
For a review of this topic see, for example, Schuster and Sigmund (1985). There is also a close relationship between dynamical systems describing autocatalytic processes and Lotka Volterra type systems of differential equations used in predator prey models (Hofbauer, 1981). Because of their widespread applicability as model systems, nonlinear autocatalytic reaction networks have been studied previously and attempts on a systematic classification of their long time behaviour have been made (Hofbauer et al., 1980; Bomze, 1983). A complete classification has been done recently for small networks based on mass action kinetics in a flow reactor (Stadler and Schuster, 1990): £Ck = Xk
ak./Xj ~ .
1
i,j=
aiyxixj
.
1n
In this contribution we extend these results to systems in which reaction rates are monotone functions of the concentration of the catalyst : ~
= xk
[email protected](Xj)j
I
aiixi~'(XJ)
'
e x ~ > O,
i,j~
2. The Model. Let us consider reactions of the form (A)+IT+IK
+ 2 I T + I x + ( B )
where IT denotes the template, Ix stands for the catalyst and (A), (B) denotes small monomers, the concentrations of which are assumed to be held constant and thus do not occur in the kinetic equations. Let us assume furthermore that the growth rate depends on the concentration of the catalyst only. For convenience we will use ck instead of [lk] for the concentration of a species Ik. The kinetic equation then becomes
We take these "elementary reactions" to be independent and in order to keep the total concentration finite, we introduce an unspecific dilution flux (I)(x). A well known experimental setup in chemical kinetics, the evolution reactor (Schuster, 1981) is designed to keep the total concentration of macromolecules constant. This special choice of the dilution flux is referred to as constant organisation. Independence of the reactions allows us to simply add up the contributions of all possible catalysts.
472
B . M . R . STADLER and P. F. S T A D L E R
The kinetic equations under the constraint of constant organization n o w read
= ok E a.q i(oj)ok*.
(1)
i=1
The dilution flux is determined by c,(t) = Co ~ Co(t). t1
If we are interested in the long time behaviour, we may eliminate Co f r o m our set of equations and rewrite it in terms of relative concentrations xi = cffco: 2k  ck  x k
(2)
~ akj¢,(Cj)Xk¢(X).
("0
1=1
F r o m ~0 = 0 we get
I,/=
1
Using this expression for the dilution flux and rescaling the time by a factor 1~co equation (2) becomes
fk

dcot

Co
Y k
1
a/c/:j(C/)i,j= 1a(/Xi~)j(Cj) 
•
We m a y introduce functions q~k defined on the unit interval by 

C0
The above definition of Ck has been chosen such that the identity on [0, Co] corresponds to the identity on [0, 1]. Our final result is 2k = xk
ak:4)/(x:)1
a,:x,49:(x:) i,j=
.
(3)
1
Thus the flux function in terms of relative concentrations reads • (x) =
i
i,j~
a~/x,¢:(X/). 1
E q u a t i o n (3) is defined on the nsimplex
i~ 1
(4)
SMALL AUTOCATALYTIC REACTION NETWORKS III
473
We restrict our investigations to strictly m o n o t o n e growth rates q~j(cj) on [0, Co]. Obviously qbk is strictly m o n o t o n e on [0, 1] iff q~k is on [0, Co]. Formally we assume that the following three conditions are fulfilled : (i) ~bj: [0, 1] * [0, 1] is continuously differentiable ; (ii) @(0) = O, and q~j(1) = I; (iii) 4~.}> 0 for all x e (0, 1). Clearly the special case 4~k(Xk) = X~ fulfils the above conditions. It is known as the first order replicator, equation. For this type of autocatalytic network we know, that there is some redundancy in the matrix o f rate constants A :
We can add an arbitrary constant to each column of A without changing the r.h.s, of equation (3) Proof. LEMMA 1.
(a~jb~)~j(xj) ~ (aijbj)xigPj(X,) j=l
i,j~l
j= l
]~1
i,j1
i,j=l
= ~akjdpj(xj) ~ aijxi49j(xj)(1~x~)'~bjq~j(xj) j=l
=
i,j=l
i=1
j=l
a~j4j(xj) ~ aijx/pj(xj). j= I
i,j= 1
As for the firstorder replicator systems we can choose a normal form for the matrix A such that ajj = 0 for all j. We will generalize Hofbauer's transform, which has proved to be very useful in the analysis of firstorder replicator networks, to our class of growth functions. Let us perform the following change of coordinates : Yi=
Xj

for" i 
1,...,t~
I.
Xn
Then we get
.gk = Yk ~ (akjan/)(~j(Xj(y))
(5)
j= 1
without any further change in velocity. This set of n  1 ODEs may be considered as a generalized version of the famous Lotka Volterra equation. Therefore we will refer to equation (3) as the replicator form and to equation (5) as the Lotka Volterra (LV) form of our model. For further calculations we have to take into account that the n variables XJ are functions of all n  1
474
B . M . R . STADLER and P. F. STADLER
variables Yr We have x~
Yk ,,_~
for
k = I,...,n1
and
x,,
1 + Z YJ
1 n1 1 + Z YJ
j=J
j~l
We remark t h a t   w i t h y,  1   f o r j = 1 , . . . , n and k = 1 , . . . , n  1 we find
~3Yk
OYk
Yi+ ~Sk,/ 1

y~
= qS}(xj(y)) •
(6)
This relation will be useful for the investigation of the case n = 3.
3. Existence and Uniqueness of an Interior Equilibrium. Let us denote the determinant of A as A := det A. F o r an interior fixed point we need (A'~b(2))k = qb(2)
for all
k = 1,...,n.
(7)
If A is regular, i.e. A ¢ 0, the linear equation A  ( ~ , . . . . , ~ , ) ~ = @  ( 1 , . . . , 1) T has a unique solution. Using Cramer's rule we get (,Ok
=
(8)
where co~ is the determinant of the matrix obtained by substituting the k t h column by the vector ( 1 , . . . , 1) r. LEMMA 2. I f there is an interior f i x e d point, then all numbers o9k are nonzero and have the same sign. Proof. Suppose 2 is an equilibrium. 2 e int S, iffgk > 0 for all k and @/A is independent of k. Therefore we have to satisfy ~ok" @/A > 0 for all k. • It remains to determine the value qb(2). LEMMA 3. For all k let COk be nonzero and let all these numbers have the same sign. Then there is a number 2 ~ ~ such that =
j=l
1.
SMALL AUTOCATALYTIC REACTION NETWORKSIll
2 Proof. L e t p =
2
for for
475
coj>0 c o j < 0 and let co0 = maxj{m/} . Furthermore we
introduce j=l
This function is defined for 0 ~
j=l
and h(O) = O. Since q~j is m o n o t o n e so is its inverse ~bf ~. Thus h is a monotonically decreasing function of #. By the theorem o f Rolle there is exactly one #* such that h(/~*) = 1 with 0
qSk(Xk) = (I)0?) " ( A  l " ( 1 , . . . , 1)r)k = 0(2)" ~ a~,
(9)
where A* denotes the inverse of A. We can substitute this into equation (4) and if ~ 0 ?) is nonzero we may divide by this term, otherwise equation (4) is satisfied trivially. In the nontrivial case we have n i,j= I
l= 1
i,l
j= l
i,I
Thus the definition of the flux function q~(x) does not give rise to an additional condition for the interior fixed point 2. Summarizing, we have proven thus far'
n
THEOREM 1. L e t A be invertible. For :~k = xk(akjY~j=l qbj(x)(I)(x)) there is an equilibrium in intSn if/" all numbers (ok have the same sign and are nonzero. Then it is unique. Remark. If A = 0 then q)(2) = 0. If d i m k e r A = 1 then (I)(2)  0 is just the missing condition for a unique determination of the rest point 2 provided the
476
B . M . R . STADLER and P. F. STADLER
determinants COk are nonzero. U n d e r this condition Theorem 1 is still valid because q)(2i/A remains finite. If, on the other hand, dim ker A ~> 2 there are at least three linearly dependent columns in A so that after substituting ( 1 , . . . , 1) 7`for one of them, the resulting matrix contains at least two linearly dependent columns, thus ogj = 0. In this case there m a y be lines, surfaces or in general hyper surfaces which consist entirely of rest points. 4. Two Species. In the special case of two interacting species it is not convenient to use the general formulae of section 3, because we can reduce our problem to a rather simple O D E on the unit interval. Let us introduce the notation
a ;), Our system reduces to
0o)
2 = x(1 x)(bq~2(1  x )  ~qS, (x)).
LEMMA 4. Let qS" [0, 1] ~ ~+ be strictly monotonically increasing, 4 [ 0 , 1] ~+ be strictly monotonically decreasing and let 4)(0) = 4(1) = 0. Then h(x) := ~(x)  ~ ( x )
(11)
has exactly one zero in the interior o f (0, 1) (ff~ > 0. It has exactly one zero on the boundary o f [0, 1] and none in the interior of this interval iff c~ = O, and has no zero in [0, 1] tff~ < 0. Proof. Let ~ < 0. Then  ~ > 0 and h(x) is the sum of two positive functions and thus positive everywhere on [0, 1]. F o r ~ = 0, h(x) reduces to O(x) and thus has a unique zero for x = 1. In the third case c~ > 0 we find h(0) = 4(0) > 0 and h(1) = ~4~(1) < 0. Since h(x) is continuous it has at least one zero in the interior of the unit interval. Since both ~(x) and ~q~(x) are monotonically decreasing this zero is unique. • Noticing that x = 0 and x = 1 are always zeros of equation (10) we conclude directly: THEOREM 2. The two species model admits exactly two f i x e d points iff sgn a ¢ s#n b, has exactly three f i x e d points iff sgn a = sgn b ¢ O. Every point in [0, 1] is rest point iff a = b = O. The stability of the equilibria is given by the derivative of the r.h.s, of equation (10) : 02 J(x)  Ox  ( 1  2 x ) [ b O ( x )  aqS(x)] + x ( 1  x ) [ b ~ ' ( x ) 
a~b'(x)]
(12)
where we used qS(x) = 01(x) and O(x) = q52(1  x ) . F o r the end points x = 0
SMALL AUTOCATALYTIC REACTION NETWORKsIII
477
a n d x = 1 we find J(0) = b,
J ( 1 ) = a.
(13)
A t the interior fixed p o i n t x*, if it exists, we find
J(x*) =  x*(1  x*)[a~', (x*) + b~b'2(1  x*)].
(14)
I f x* exists, a a n d b h a v e the same sign a n d the derivatives o f ~bl and q~2 are positive e v e r y w h e r e in the interior o f the unit interval. T h u s sgn J(x*) =  sgn a =  sgn b. W e can state LEMMA 5.
If the two species model admits three fixed points, the interior is
stable iff a, b > 0 and it is unstable iff a, b <~O. The equilibria on the endpoints of the interval have the opposite stability. I f there are exactly two rest points they have opposite stability. F i g u r e 1 shows a c o n s t r a i n t  c o n s t r a i n t p l o t (a p h a s e d i a g r a m ) for the two species model. ft. T h r e e Species. W e a d o p t the same n o t a t i o n as in Stadler and Schuster (1990) for the r e a c t i o n m a t r i x A •
A =
0
0
.
(15)
•
)
a
~C~
Figure 1. Phase diagram for the two species model. Sinks are denoted by a bullet, circles denote sources.
478
B.M.R. STADLER and P. F. STADLER
For det A and the determinants ok we find
coj
=
fl~+&/?~
02 = e y + e ~   ? ~ o03 = ~ f l + ¢ a  ~ a .
(16)
Applying Hofbauer's t r a n s f o r m   e q u a t i o n ( 5 )   w e obtain the following system on ~2+. = x(~
l (Xl) + ( a  g)~2(x2) + yq~,(x,)) ,=
xq'(x, y)
? = y((o~~)dp,(x,)fl(02(Xz)+e43(x3)),=y®(x,y).
(17)
At the interior fixed point (2j/23, 22/23) we find ~F = ~ = 0. The Jacobian then reads
(2Udx
2~,~.
(18)
A rather lengthy calculation yields an explicit expression for the determinant of the Jacobian at the position of the interior fixed point. After transforming back to the simplex coordinates we find det J(2,fO = 2~1x223 • [~1 (21)(~2 (22)(03 "~ ~ l (21)q~; (23)(D2 "3i~ 2 ( 2 2 ) ~ t 3 (23)(D 1] • Note that we reduced our system by H o f b a u e r ' s transformation to a flow equivalent system in the plane and then transformed back to the simplex. We have done this to restrict our analysis to the two physically meaningful directions of our simplex, i.e. to exclude the external eigenvalue. The same procedure will be done for the trace of the Jacobian. F r o m Section 3 we know, that if there is a unique fixed point in int $3 then the numbers cOjhave equal signs and are nonzero ; by definition the derivatives ~b;~(2k) are positive. If det J(2, f~) is negative then the two relevant eigenvalues have opposite signs and therefore the interior fixed point is a saddle. If on the other hand det J(2, .~) is positive the eigenvalues must have equal signs and the interior rest point is a source, a sink, or a centre depending on the sign of the trace o f J(2, 39). We have : THEOREM 3. The interior fixed point & a unique saddle point iff all numbers (o~
are negative. THEOREM 4. I f C01,~02, O03 < 0 then there are no periodic orbits in $3. Proof. If there are periodic orbits then they must surround an equilibrium which is not a saddle. For $3 and all subsimplices on its boundary, i.e. all
SMALL AUTOCATALYTIC REACTION NETWORKSIll
479
edges and corners are invariant under the flux of our model, this fixed point must lie in the interior of $3. We know that the interior equilibrium is unique and that it is a saddle. • The calculation of the trace of J(2, f~) yields after rewriting in simplex coordinates Tr J(2, )9) =  [2101 (21) ((x22 ~~23) }220'2 (22) (621 }/~2~3)
There is another representation of this expression which might be useful for practical calculations : Tr Y(2, 9)) = (
[email protected](~))  @(2).
(19)
Let us now have a look at the edges of the state space. Without losing generality we can choose the ( 1 , 3 ) edge, i.e. x2 = 02(x2) = 0. In one dimensional formulation the restriction to this subsimplex reads ~ = X (   ~ 0 1 ( X ) ' q  ~ 0 2 ( 1   X ) ) .
F r o m section 4 we know that there is an interior equilibrium iff7 and ( have the same sign. The eigenvalue along the edge has the sign of 7. Let us now calculate the trace of the Jacobian. After rewriting in simplex coordinates we find 1 TrY(x*,O)=xTx~({O'1(x~)+70"3(x*))+ co201(x*). (20) 7 It is easy to calculate the eigenvalue along the edge since the eigenvector, i.e. the direction of the edge, is known. We find that
21 =  x*x~( ~O"~(x'~) q 743 (x*) )
(21)
belongs to the edge. Thus the transversal eigenvalue is given by 1 ~2 =
1 (D201(Xl~) =
('0303(X~)"
(23)
7
Table 1 summarizes the behaviour of fixed points in the interior of edges. It Table 1. Signs of the eigenvalues for fixed points on the edges of the $3 Edge
(1,2)
(1,3)
(2, 3)
Condition
sgne=sgn6
sgnT=sgnff
sgnfl=sgne,
2~ 22
sgnc~ sgn (c~co3)
sgn7 sgn (?(D2)
sgnfl sgn (rico,)
480
B. M R. S T A D L E R a n d P. F. S T A D L E R
Table 2. Eigenvalues for corner fixed points of the S~ Eigenvalue
e~
e2
e3
along (1,2) along (1,3) along (2, 3)
~ ~ 
6 fl
7 e,
remains to analyse the corners. W i t h o u t losing generality we select e3. The Jacobian is diagonal : J(e3) = (~
0)
(24)
The stability of the corners is tabulated in Table 2. We emphasize that the flux on the b o u n d a r y matches exactly with the behaviour of the three species first order replicator equation. This statement also holds for int $3 if there is no interior fixed point or if the interior fixed point is a saddle. In s u m m a r y we have proved so far THEOREM 5. I f not all three numbers o2 ~, co2, o93 are positive then there is a homeomorphism which maps the phase portraits o f equation (3) onto the phase portraits o f the first order repiicator equation with the same matrix Of rule coejficients A. 6. Power Laws. Let us n o w restrict our consideration to @(xj) = xJ'
forall
p e r +.
(25)
It is easy to get a closed expression for the coordinates of the interior fixed point. We will restrict ourselves to the case o91, o92, co3 > 0 for we already know the behaviour in the other cases. We find (26)
=
P
P
P '
A n explicit formula for the trace o f the Jacobian can also be found : Tr J(2) =  p
A (P,f _ _ +px/co +p~)p
p.
(27)
Thus we have sgn Tr J =  sgn k. A H o p f bifurcation occurs for A = 0. We will see that it is always degenerate, i.e. it does not give rise to stable or unstable limit cycles.
SMALL AUTOCATALYTIC REACTION NETWORKSIll
481
THEOREM 6. The three species selection model equation (3) with ~bj(xj) = x~ does not admit isolated periodic orbits f o r any p ~ R + . There is a continuum oJ periodic orbits iff A = 0 and co t, co2, o~3 > O. Proof. The LVform of this model reads =
x
(1 + x + y)P " (  ~xp + (6  fl)yP + ~)
Y )) = (1 + x + y)P " ((c~ ( ) x p  flyP + e).
(28)
We can get rid of the factor (1 + x + y ) p by a change of velocity. With a second change of coordinates, Y = x p and fi = yP, we find x = p x p  1 ~ and = pyP ~ . If we change velocity once m o r e   b y a factor 1/pwe are left with x = Y (  (2 + (3  fl) )~+ 7) (29) ¢ = y.
=
This is the well k n o w n Lotka Volterra equation. The first part of theorem 6 is a classical result on this equation; see Coppel (1966) or Andronov et al. (1973) for a proof. The proof of the existence of a continuum of periodic orbits in equation (29) iff the trace of the Jacobian at an interior nonsaddle equilibrium vanishes, can be found in H o f b a u e r and Sigmund (1988). A trivial consequence of theorem 6 is the l\~llowing. For n <~ 3 there is a homeomorphism mappit~q the pha,w portraits of the replicator equation (3) with the growth [unctions ~/~ d~ftined iH equation (25) onto the P P s o f the first order replicator equation with the same matrix o f reaction constants A. CORO1,LARY.
7. Numerical ExperimentsLimit Cycles. A generalisation of the results of the previous section seems to be quite difficult, for there is no general recipe to find Dulac functions for models with arbitrary growth rates ~k(xk). Trajectories have been integrated using the I M S L library subroutine D I V P A G with Gear's algorithm for stiff ODEs. All calculations have been performed on an IBM 3090 mainframe. Since there are no limit cycles in the case of power laws, provided there is a single power p for all three functions 4~k, we expected that there are no limit cycles if we choose 4 , ( x ) = q 2(x) = 43(x)  4(x).
(30)
We have tested this conjecture for a few functions, for example q~(x) = a sin x or 4~(x) = b In (1 + x ) . We could not find limit cycles in these cases.
482
B . M . R . STADLER and P. F. STADLER
If we deal on the other hand with three arbitrary m o n o t o n e functions ~bk, then the cubic coefficient a" determining the stability of the interior fixed point at a H o p f bifurcation contains so m a n y independent parameters, i.e. second and third derivatives of the three functions, that it is hard to imagine that 1
a" fxxx+J~,y+gpyy+gy~x+ 09 [fxy(f~+ fyy)gxy(gxx+gyy).f~xgx~+ fyyg~.~] =
vanishes at all H o p f bifurcations point. In the above formula f and g denote the nonlinear part of the H o f b a u e r  t r a n s f o r m e d system after a second (linear) transformation of coordinates has been performed in order to bring equation (17) to the form
\g(x,y)/ f~x etc. denote partial derivatives with respect to these coordinates x, y at the position of the interior fixed point. Clearly the first three derivatives o f qS~, q52 and q~3 at the position of the interior fixed point enter as parameters into a'. A nonzero value o f a" would prove the existence of a limit cycle on one side o f the H o p f bifurcation (see e.g. Gluckenheimer and Holmes, 1986). In fact we obtained a limit cycle for power laws with different exponents : (~1 (X1) :
.X.1,
(~2(X2) = X4,
(~3(_~.3) = .¥7
and symmetric reaction constants e = fl = 7 =  l, 6 = e = ~ = 0.88 Figure 2 shows a phase portrait for this case as well as two examples of centres occurring at degenerate H o p f bifurcations in the case of power laws. We close this investigation with a remark on where one can expect to find H o p f bifurcations in parameter space, i.e. in the space of the six independent reaction constants. O f course there are no H o p f bifurcations if the interior fixed point is a saddle or if there is no interior fixed point. F r o m the remaining part of the parameter space there are two cases which we can rule out: Tr J(2, ))) cannot change sign in a region where all reaction constants have a c o m m o n sign. For functions which do not differ too m u c h from linear functions we expect the hypersurfaces {qs ~6[TrJ(20; q) = 0} to be close to the hypersurfaces {q e N6[A(q) = 0}. (The additional argument q e N6 in J(20; q) has been introduced in order to emphasize that we deal with hypersurfaces in parameter space.) F o r a map of this parameter space for firstorder replicator equat i o n s  a n d thus also for the power laws described in section 6   w e refer to Stadler and Schuster (1990). 8. Conclusions. The studies reported here show that we can expect the firstorder replicator equation to be a good model to describe the evolution of a
t%l
483
x
:2 7
c~ : o= ~
~ ;.= N II
,..= <. N ~ ~: I1 ~ 
J,
v
L~
'
o
II = ~ m c~
X
o ~
~
"~
C22
,.
<
0
I"'
~
il~
t",I il ~
II "~
))) '"
~,O~
,,.i
v
ca
<, ¢' ~
X
©
II
484
B, M. R. STADLER and P. F. STADLER
small autocatalytic network even if the growth rates are not exactly linear but only monotone. A large class of rate functionsincluding MichaelisMenten (or Holling) type interactionyields the same asymptotic behaviour as networks with up to three species based on linear growth functions. This implies that the firstorder ansatz is sufficient, e.g. to decide whether such a system o f interacting species can be cooperative or not and to find out number and type of equilibria. On the other hand the firstorder ansatz is not appropriate to model systems in which the response o f a species to the concentrations o f the different catalysts does not have the same functional dependence for all catalysts at least in a substantial part of the parameter space. In this case limit cycles may occur also for n = 3, contrary to the model with linear growth functions. Stimulating discussions with Prof. Dr P. Schuster, Doz. Dr J. Hofbauer and Dr J. Swetina, and the generous supply with computer time by the "EDVZentrum der Universitfit Wien" are gratefully acknowledged. The authors wish to thank Mr J. K6nig for drawing the figures.
LITERATURE Andronov, W. W., E. Leontovich, I. Gordon and A. Maicr. 1973. Qualitatit;e 17woJ7 ojSecot~d Order Dynamic Systems. Halsted Press, New York. Biebricher, C. K., M. Eigen and W. C. Gardiner. 1983. Kmetic~ of RNA replication. Biochemistry 22, 2544~..2559. Biebricher, C. K_, M. Eigen and W_ C. Gardiner. 1984. Kinetics of R N A replication : plus minus asymmetry and double strand formation. Biochemislry 23, 31863 i 94. Biebricher, C, K., M. Eigen and W. C. Gardiner. 1985. Kinetics of R N A replication: Competition and selection among selfreplicating R N A species. Biochemistry 24, 6550 6560. Biebricher, C. K. and M. Eigen. 1988. Kinetics of R N A Replication by Qflreplicase. In: RNA Genetics (Vol_ 1), E. Domingo, J. J. Holland and P. Ahlquist (eds), pp. 122. CRC Press, Boca Raton, Florida. Bomze, I. M. 1983. Lotka Volterra equation and replicator dynamics. A twodimensional classification. Biol. Cybern. 48, 201211. Cech, T. R. 1986. R N A as an enzyme. Sci. Am. 255(5), 7684. Cech, T. R. 1988. Conserved sequences and structures of group I introns : building an active site for R N A catalysis. Gene 73, 259271. Coppel, W. 1966. A survey of quadratic systems. J. Diff. Eqns. 2, 293304. Doudna, J, A. and J. W. Szostak. 1989. R N A catalysed synthesis of complementarystrand RNA. Nature 339, 519522. Eigen, M. 1971. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften 58, 465523. Eigen, M. and P. Schuster. 1978. The Hypercycle. A principle of natural selforganisation. Part B : The abstract hypercycle. Naturwissenschq/?en 65, 7 41. Guckenheimer, J. and P. Holmes, 1986. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Applied Mathematical Sciences Vol. 42. Berlin : Springer. Hofbauer, J., P. Schuster and K_ Sigmund. 1979. A note on evolutionary stable strategies and game dynamics. J. theor. Biol. 81,609612.
SMALL AUTOCATALYT1C REACTION NETWORKSIll
485
Hofbauer, J_, P. Schuster, K_ Sigmund and R. Wolff. 1980. Dynamical systems under constant organization II : homogeneous growth functions of degree p = 2. S I A M J. appl. Math. 38, 282304. Hofbauer, J. 1981. On the occurrence of limit cycles in the LotkaVolterra equation_ Nonlinear Analysis 5, 1003 1007_ Hofbauer, J_ and K. Sigmund. 1988. The Theory of Evolution and Dynamical Systems. Cambridge University Press. MaynardSmith, J. 1982. Evolution and the Theory of Games. Cambridge University Press. Perelson, A. S. (ed.). 1988. Theoretieal Immunology. Vols. I and II_ Redwood City, C A : AddisonWesley. Schuster, P., K. Sigmund and R. Wolff. 1978. Dynamical systems under constant organisationI. Topological analysis of a family of nonlinear differential equationsa model for catalytic hypercycles. Bull. math. Biol. 40, 743769. Schuster, P., K. Sigmund and R. Wolff. 1981. Mass action kinetics of selfreplication in flow reactors. J. Math. Anal. Appl. 78, 88112. Schuster, P. 1981. Selection and Evolution in Molecular Systems. In Nonlinear Phenomena in Physics and Biology, R. H. Enns, B. L. Jones, R. M. Mimura and S. S_ Rangnekar (eds), pp. 485543. New York : Plenum Press. Schuster, P. and K. Sigmund. 1983. Replicator Dynamics. J. theor. Biol. 100, 533 538. Schuster, P. and K. Sigmund. 1985. Towards a dynamics of social behaviour : strategic and genetic models for the evolution of animal conflicts. J. social_ biol_ Struct. 8, 255277. Stadler, P. F. and P. Schuster_ 1990. Dynamics of small autocatalytic reaction networks 1. Bifurcations, permanence and exclusion. Bull Math. Biol. 52, 485508. Taylor, P. D_ and L. B. Jonker. 1978. Evolutionary stable strategies and game dynamics. Math. Biosci. 40, 145156. Weissmann Ch. 1974. The making of a phage. FEBS Letters 40, S I0 SI 8. Zeeman, E. C. 1981. Dynamics of the evolution of animal conflicts. J. theor Biol. 89, 249270.