BioSystems 7 (1975) 9 2 - - 1 0 0 © N o r t h - H o l l a n d Publishing C o m p a n y , A m s t e r d a m - - Printed in The N e t h e r l a n d s

ANALYTIC PROCEDURES FOR L A R G E DIMENSION N O N L I N E A R BIOCHEMICAL O S C I LLATOR S PAUL RAPP

Department of Applied Mathematics and Theoretical Physics, Cambridge University, Cam bridge, England

A n analytic m e t h o d is p r e s e n t e d which can be used to d e t e r m i n e if the following s y s t e m o f n o n l i n e a r differential e q u a t i o n s has p e r i o d i c s o l u t i o n s Xl = h ( x n ) - - b l X l

xi = gi--lXi--1 - - bjxi

J = 2, ... n

A s y s t e m a t i c dual i n p u t describing f u n c t i o n p r o c e d u r e is given for c o n s t r u c t i n g a f u n c t i o n o f t h e reaction c o n s t a n t s R, w h e r e if R > 1 a p e r i o d i c s o l u t i o n exists and if R < 1 t h e r e is no periodic solution. The f o r m o f R c o n s t r u c t e d generalizes i m m e d i a t e l y to an arbitrarily large d i m e n s i o n . The m e t h o d generalizes to cover s y s t e m s displaying h y s t e r e s i s kinetics, s y s t e m s s u b j e c t t o chemical noise, and s y s t e m s c o n t a i n i n g delay c o m p o n e n t s . The m e t h o d has been applied to a well k n o w n b i o c h e m i c a l p r o b l e m w h e r e h ( x n ) = k/(1 + 0/xnP ). For p = 1, for all n, t h e r e are no stable limit cycles such t h a t xi(t ) > 0, t >/ 0. For p = 2, n/> 8 it is possible to c o n s t r u c t a p a r a m e t e r set such t h a t stable oscillations appear.

1. I n t r o d u c t i o n The biologically i m p o r t a n t properties of a biochemical system are the dynamic properties such as oscillations and transitions between steady states, n o t the steady states themselves. Theoretical studies of dynamic behavior have been impeded by a lack of mathematical tools which can be successfully applied. Accordingly most researchers have resorted to simulations f o r small sets of parameter values. Since the systems are large and highly nonlinear the simulation approach can fail to detect import a n t features o f the system, or may give unwarranted emphasis to properties which are limited to a very small domain in the parameter space. An area o f theoretical biology where this has been especially true has been in the study of biochemical oscillators. In this paper an analytic m e t h o d for locating periodic solutions to a system o f nonlinear differential equations describing a biochemical system is developed. Specific examples of its use are given in the analysis o f a nonlinearity which has been investigated b y Goodwin, Walter, Griffith and others.

An approxi m at e m et hod, the dual input describing function, is given for constructing a funct i on of the reaction constants, R, where if R > 1, t he system will oscillate and if R < 1, the system will n o t oscillate. The error in R is less than 5% for the problem considered here. However, the accuracy improves as the num ber of chemical reactants in the system increases. Thus the m e t h o d complements existing small dimension techniques such as Poincar~-Bendixson t h e o r y which is restricted to two dimensions. Other advantages of the describing function should be noted. First, and most important, the funct i on R is derived in a form that is valid for all dimensions. It is unnecessary to reconstruct R each time an additional step is introduced. Most techniques must be reapplied for each dimension. Second, it is possible to state a systematic procedure for the use of the m e t h o d unlike Liapounov techniques where it is impossible to specify a systematic program which will lead to a successful Liapounov function. The m e t h o d can be e x t e n d e d to cover more complicated systems. These extensions are considered in the final section.

93 2. The dual input describing function

(P+bl)(P+b2)(p+ba)...(p+bn)xn

The kinetic structure of the system analysed here has been pre,;ented in its most general form by Walter (1969a, 1969b, 1970, 1971). It consists of a slow reaction, formation of compound x l , followed by a sequence of fast reactions. The product of a fast reaction, Xn, inhibits the slow reaction. This scheme is probably the most c o m m o n feedback loop found in biological systems. As it is particularly simple it was probably one of the first feedback reaction schemes to appear in a prebiotic system. It was assumed that no important features of the system were lost by linearizing the equations describing fast reactions. Since it is usually the case that the relaxation times of the fast and slow reactions differ by several orders of magnitude, this linearization is justifiable. If it is true that the source material used for the construction of xl is slowly varying in time, then the rate of synthesis depends only on the value of Xn, the concentration of the feedback metabolite. The differential equations describing the scheme are given in equation 1, where bj's and gj's are positive reaction constants and x i's are positive concentration functions.

gn-- lh(Xn)

-- f ( X n )

= g l ... (2)

This system can be represented by a closed feedback loop of the form shown in the transfer function block diagram given in Fig. 2, where f(Xn) is the nonlinear element in equation 2 and G(p) is the linear filter n

G(p) = 1-I 1/(b~ + i¢o) j=l

The dual input describing function technique will now be used to analyse this feedback system. The method is a standard one in the control theory literature (Gelb and van der Velde, 1968; Hsu and Meyer, 1968). As it is applied here, the method differs from that usually encountered in that it is used in its analytical form rather than as a numerical tool. The limit cycle and o u t p u t of the nonlinear element are approximated up to first order harmonics. The limit cycle is approximated by the expression y + x cos ~0t where one requires y > x > 0 in order to insure positive concentrations. The o u t p u t of the nonlinear element is approximated by the first terms of its Fourier series (Note that a2 = 0 if f is single valued).

dx 1

dt = h(xn ) -- bl xl

f(xn ) = a0/2 + a 1 cosc0t + a2 sincot

dx2

dt

2 -

gl x1 - - b 2 x 2

dx n dt - g " - 1 Xn-

1

(1)

--

1

ao = ~- f o 1 2n

bn x . al

Let p denote the differentiation operator. Using the linearity of p gives the following equivalent expression governing the variation in Xn ;

f(y + xcosO)dO

=;f

f(y + xcosO)cosOdO

o

S'T x1~*x2---" . . . . . . . --~ 'xn - - T Fig. 1. T h e x I f o r m a t i o n r e a c t i o n is i n h i b i t e d b y f e e d b a c k m e t a b o l i t e xa.

Fig. 2. Transfer f u n c t i o n block d i a g r a m c o r r e s p o n d ing to c h e m i c a l r e a c t i o n s y s t e m in Figure 1. G ( p ) is a linear filter, f is a n o n l i n e a r f u n c t i o n .

94 1 a2 = ~- f 0

Q(x,y)G(o) = 1

2~

Z e r o t h Balance E q u a t i o n (3)

f(y + x c o s 0 ) s i n 0 d 0

A balance e q u a t i o n c o r r e s p o n d i n g t o the l o o p e q u a t i o n (fG)xn = Xn is w r i t t e n separately f o r each h a r m o n i c c o m p o n e n t o f the f(xn) o u t p u t . If it is possible to satisfy these balance e q u a t i o n s s i m u l t a n e o u s l y t h e n a periodic solut i o n exists. The f o r m o f the e q u a t i o n corres p o n d i n g to t h e z e r o t h h a r m o n i c is established in the following manner. T h e n o n l i n e a r elem e n t is replaced b y Q, its equivalent f u n c t i o n in t h e zero h a r m o n i c . Given i n p u t equal to the c o n s t a n t y, t h e e l e m e n t Q has o u t p u t equal to t h e c o n s t a n t ao/2. Thus Q is itself a gain c o n s t a n t d e f i n e d b y the relation Qy = a 0/2, (Q = a 0 / 2 y ) . T h e c o n s t a n t o u t p u t o f Q is o p e r a t e d o n b y G(p). Since a 0 / 2 is a c o n s t a n t , its derivative is zero; thus G ( p ) ( a o / 2 ) = G ( O ) (a0/2). T h e c o n d i t i o n which must be satisfied to have o u t p u t y result f r o m i n p u t y at e l e m e n t Q is t h e z e r o t h balance e q u a t i o n ; QG(O) = 1. The n e x t balance e q u a t i o n is established b y a similar a r g u m e n t . T h e o u t p u t o f the l o o p in the first h a r m o n i c is xcoscot. Nonlinear e l e m e n t f is t o be replaced b y e l e m e n t P where i n p u t xcoscot t o c o m p o n e n t P gives o u t p u t al coscot + a2sincot. Substantial simplification can be achieved if t h e p r o b l e m is r e s t a t e d in terms o f c o m p l e x quantities. xcoscot = [:{(xeiw t) al coscot + bl sincot = R [(al -- ia2 )e iw t ] In c o m p l e x f o r m P is d e f i n e d b y the relation, i n p u t x e i¢° t gives o u t p u t (al -- ia2 )e i~ t. This fixes P as (al -- ia2 )/x. P will be r e f e r r e d to as t h e describing f u n c t i o n . The objects o p e r a t e d u p o n b y G(p) all have the f o r m Ce i~ t. C is a c o n s t a n t . T h e derivative is icoCeiC°t; thus app l y i n g t h e d i f f e r e n t i a t i o n o p e r a t o r is t h e same as m u l t i p l y i n g b y ico. G(p) can be replaced b y G(ico). T h e balance c o n d i t i o n requires t h a t i n p u t x e i ~ t t o e l e m e n t P results in l o o p o u t p u t x e ic°t. T h e p r o d u c t o f P and G(ico) m u s t be one. T h e resulting e q u a t i o n s are

P(x,y)G(ico) = 1

First Balance E q u a t i o n (4)

Q(x,y) = ao/2y

(5)

P(x,y) = (al -- ia2 )/x

(6}

T h e z e r o t h e q u a t i o n is used t o fix y as a f u n c t i o n o f x. Eliminating y d e p e n d e n c e , the n e x t e q u a t i o n b e c o m e s a f u n c t i o n o f x, the a m p l i t u d e o f the oscillation, and co, the freq u e n c y . An a t t e m p t t o find a solution to this e q u a t i o n is m a d e b y considering the behavior o f t h e G(ico) and 1/P c o n t o u r s in the c o m p l e x plane. If an i n t e r s e c t i o n occurs t h e n the balance e q u a t i o n s are satisfied. It is i m p o r t a n t to observe t h a t G(ico) is the N y q u i s t c o n t o u r o f the linear e l e m e n t . Thus b y a p p l y i n g the Nyquist stability criterion to analysis o f the l / P , G(ico) intersections, one has i m m e d i a t e l y a p r e d i c t i o n o f the stability o f t h e limit cycles. Stability is d e t e r m i n e d b y considering perturbations o f the p o i n t o f i n t e r s e c t i o n in each d i r e c t i o n along t h e 1/P line. T h e s y s t e m will r e t u r n to, or m o v e away f r o m , the limit cycle a c c o r d i n g to w h e t h e r or n o t the p e r t u r b e d p o i n t is encircled b y G(ico) and b y the d i r e c t i o n o f increasing x on the 1/P line. T h e third diagram gives t h r e e examples o f c o m m o n l y e n c o u n t e r e d behavior for systems possessing real describing f u n c t i o n s . The arrows on the 1/P lines indicate t h e d i r e c t i o n o f increasing x. T h e arrows on the G(icn) c o n t o u r s indicate the d i r e c t i o n o f increasing co. In t h e first case the describing f u n c t i o n is m o n o t o n i c a l l y increas-

Fig. 3. The first example gives an intersection corresponding to an unstable limit cycle. The second example corresponds to a stable limit cycle. The two intersections in the third example each correspond to a limit cycle. One is stable, the other is unstable.

95 ing. The intersection corresponds to an unstable limit cycle. In the second example the describing function is monotonically decreasing and the intersection corresponds to a stable limit cycle. In the third case each of the two intersections represents a different limit cycle. The first is unstable and the second is stable. Results concerning the accuracy of the conventional single input describing function carry over intact to the dual input form presented here (Gelb and van der Velde, 1968; Hsu and Meyer, 1968). Justification of the expansion of the expressions for the limit cycle and the o u t p u t of the nonlinear element to only first harmonics follow from the properties of G(p) which acts as a low pass filter (suppresses higher harmonics). Several properties of the G(p) function encountered in this analysis indicate that the procedure will be especially accurate for this problem. ]:n general the G of most systems is expressible as a rational function in p. The filtering properties of a well behaved G(p) increase geometrically with A, where A is the difference in the: order of p in the denominator of G and the order of p in the numerator of G. For the G(p) encountered in the analysis of equation 1, the numerator is zeroth order in p. Thus A = n where n >~ 3. Most biological systems have a large number of reactants. This means that A will generally be quite large and that G(p) will have marked low pass filtering properties. It has been pointed out that in almost all cases where the describing function fails G(ico) has high resonant peaks (Hsu and Meyer, 1968). Such peaks are not present in the well behaved G(p) studied in the paper. Also in almost all counter examples A = 1, where here it is at least t h r e e If the m e t h o d fails it is usually because the system possesses an oscillation that has a large c o m p o n e n t in the second or third harmonic. Thus borderline cases can be resolved by deriving the balance equation for these harmonics and solving them simultaneously with those; for the zeroth and first harmonics. Accuracy of the describing function is generally ~aken as being 10%. Here extensive numerical tests have found the error to be less than 5%.

3. Applications The m e t h o d as presented will now be applied to the system with h(x n ) as shown for the cases p = 1 and p = 2. As the details have appeared elsewhere (Rapp, 1974) only an outline is given here. This system of equations has already been examined theoretically by several workers, notably Goodwin (1965), Walter (1969a, 1969b, 1970, 1971), Griffith (1968a, 1968b), Morales and McKay {1967), and Tiwari and Fraser (1973, 1974). All of these papers have been concerned with periodic solutions. Griffith (1968a, 1968b) has shown that for n = 1 and n = 2 no limit cycles exist. Walter has approached the question via Liapounov techniques (1969a, 1969b) and extensive simulation work {1970). These studies have yielded information about the singular points but did not succeed in establishing a set of conditions which would ensure the presence of limit cycles. Since elementary techniques indicate that the system can not oscillate for n < 2, we consider here only systems where n >~ 3. k h(xn) - 1 + ux~p f(Xn) = glg2 "'" g . - l k

_

1 + ax.p

S = {k, a, g, ... g n -

1, bl

d=

(7)

1 + d2xnP ,..

bn }

S is said to be permissible o n l y if each element is positive. The positive constants dl and d2 are defined as shown in equation 7. Also the constant Co, e0 = bl b2 ... bn = 1/G(O) is positive. Some important properties of the integrals ao,al, and a2 can be stated for general p. First since f is single valued a2 = 0. This means that the describing function is real. It can also be shown that al is always negative and thus 1/P is negative. The zeroth balance equation takes the form 1

1 {1 2~ 61 dO 2CoY f 1 + d2(y + xcos0) °

}

(8)

0

A solution y(x) to this equation is permissible if

96 y ( x ) > x > 0. It can be s h o w n t h a t e q u a t i o n 8 has a p e r m i s s i b l e r o o t if a n d o n l y if x is in t h e interval ( 0 , x * ) w h e r e x * is d e f i n e d b y t h e equation

1{1/

0=1

2X'Co

0

1 + d 2 x * ( 1 + c o s 0 ) p d0

}

T h e analysis f o r p = I n o w follows. Using t h e f o r m o f t h e z e r o t h b a l a n c e e q u a t i o n gives e q u a t i o n 9 w h i c h fixes y as a f u n c t i o n o f x . T h e e q u a t i o n 1/P is f o u n d b y evaluating a l .

(c0y)2(1 + d2Y)2 --(Coy)2(d2x)2 - - d l2 = 0

(9) 1 P(x)

Cj = t a n - 1 (¢o/b~)

(12)

Since all bj's are positive, as ~ increases f r o m zero, b o t h the m a g n i t u d e o f G(ico) a n d its a r g u m e n t are m o n o t o n i c a l l y decreasing. G(ico) spirals in t o t h e origin in a c l o c k w i s e d i r e c t i o n f r o m G ( O ) = 1/Co. This m e a n s t h a t if an i n t e r s e c t i o n exists it is o f t h e f o r m o f t h e first e x a m p l e in Fig. 3 a n d t h u s t h e c o r r e s p o n d i n g limit cycle is an u n s t a b l e one. Define COc as t h e f r e q u e n c y value at w h i c h G ( i w ) first intersects t h e negative real axis. An i n t e r s e c t i o n will exist o n l y if i n e q u a l i t y 13 is satisfied. F o r p = 1 a n d p = 2 t h e m i n i m u m value o f 1 / P is necessarily less t h a n G ( i w c ) f o r all p e r m i s s i b l e p a r a m e t e r sets. T h u s satisfying i n e q u a l i t y 13 a n d its p = 2 a n a l o g ensures t h e e x i s t e n c e o f an intersection. 1 _ --IG(ic~c)i = G(icoc) < ~¢¢LimP(x)

d2 x2 [(1 + d2y) 2 -- (d2x) 2] i/2 2dl[(l + d2y) 2 - (d2x) 2] 1/2 __ 2di(I + d2Y)

(13)

(lO) 1/P(x} is a m o n o t o n i c a l l y increasing f u n c t i o n o f x a n d it has t h e limits s h o w n below. T h e limit o f 1 / P ( x ) as x a p p r o a c h e s i n f i n i t y is a valid limit a n d c o n s t i t u t e s t h e m a x i m u m value o f 1 / P ( x ) f o r p = 1 even t h o u g h t h e value o f y is less t h a n x f o r x is g r e a t e r t h a n x*.

' LG(0)I

D e f i n i n g R as s h o w n we see t h a t an intersect i o n can n e v e r o c c u r if R < 1. T h e s e c o n d expression for R follows from the definition o f Cj. n

R = 21G(iwc)/G(0}l = 21-[ cos Cj

(14)

j=l

1 Lim x-~0 P(x)

--(I + d2Yo)2 dl de

--c o + [c 2 + 4 c o d l d 2 ] 1/2 Yo =

2cod 2

1 1 = - - ! IG(0)I 2 Lim p(x~-)= --2C--o T o c o m p l e t e t h e analysis f o r p = 1, o n e m u s t consider the behaviour of G(i~). n

G(ico) = 1-] 1 _ j= 1 (bj + ico)

lel 1 ) e-i(ol+'''+~n) j=l (bj2 + w2)112

I t is n o w n e c e s s a r y to find t h e p a r a m e t e r set which gives the optimum value of IG(i~c)/G(O)I s u b j e c t t o the c o n s t r a i n t co = coc (argG = --~r). This can be d o n e either b y t h e m e t h o d o f Lagrange m u l t i p l i e r s or b y c o n s i d e r i n g t h e e q u i v a l e n t p r o b l e m o f optimising log I G ( i ~ ) / G ( O ) I and using t h e concave p r o p e r t i e s o f log(cos ~bj). T h e o p t i m u m is r e a c h e d w h e n bj = b* > 0 f o r all j. T h e o p t i m u m value o f R is (n = d i m e n s i o n o f system) Rop t = 2{cos(Tr/n)} n

(11)

N o t e t h a t t h e o n l y r e s t r i c t i o n o n b* is t h a t it be positive. Ro p t is d e t e r m i n e d i n d e p e n d e n t l y

97 o f t h e individual values o f e l e m e n t s o f S. It d e p e n d s o n l y on t h e d i m e n s i o n o f t h e system. F o r n >~ 3 Ro p t m o n o t o n i c a l l y increases with increasing n. For n ~< 7, Ropt < 1 and for n ~> 8 R o p t > 1. Thus for these higher dimensions it is possible to c o n s t r u c t a p a r a m e t e r set which yields a s y s t e m possessing an unstable limit cycle. T h e analysis for the case where p = 2 follows in a similar manner. G(ico) is the same as f o r t h e p = 1 case as are t h e definitions o f d~, d2 and Co. T h e zel:oth balance e q u a t i o n and 1 / P ( x ) are as follows

It can be s h o w n t h a t Yo = ( x . ) e , t h e equilibr i u m value o f Xn (t). T h e oscillation c o n d i t i o n is analogous to t h a t c o n s i d e r e d in t h e p = 1 case. T h e f u n c t i o n R m u s t be greater t h a n one if an oscillation is to occur.

X~U P(x)

--dI 2d2co2y3 R=

dl

0=1

2

2 d l ( s i n ~ ' / 2 - - x / ~ 2 Y cos 7 / 2 )

d2x2 --d2Y2 + 1) 1/2 K

K = [(d2 x2 - - d 2 Y 2 + 1) 2 + 4d2YZ] l/2 1 + COS "7/2 = +

sin 7 / 2 = +

cos ,,/)1/2 2

1 -- cos. 7)1/2 2

T h e f u n c t i o n 1 / P ( x ) is m o n o t o n i c a l l y decreasing with increasing positive x. T h u s if t h e r e is an i n t e r s e c t i o n it is o f t h e f o r m which w o u l d indicate the e x i s t e n c e o f a stable limit cycle. T h e limit o f 1/P as x a p p r o a c h e s zero is: 1

Lim x-~0 P(x)

+

2d2co

(xn)

) G(ic°c)[ G--i

=

(15)

Again it is necessary to find the p a r a m e t e r set t h a t yields t h e o p t i m u m value o f R. T h e fact o r U and t h e f a c t o r IG(icOc)/G(O)l b o t h dep e n d o n t h e values o f bj(co = bl ... bn) and t h u s in general one c o u l d n o t o p t i m i z e the t w o factors i n d e p e n d e n t l y . H o w e v e r it was f o u n d t h a t the o p t i m u m o f IG(icOc)/G(O)[ was r e a c h e d w h e n all t h e bi's were equal. T h e value o f Co was c o m p l e t e l y u n s p e c i f i e d at the e x t r e m u m . Thus o n e can consider the maxim u m values o f U and o f IG(icOc/(G(O)L separately. Defining t h e variable A as A = d 1 V ~ 2 / x / 2 c 0 o n e can s h o w t h a t the f o l l o w i n g holds:

U(dl, d 2 , 2-

~+

Co) =

U(A)=

~+27As

!

2d2co2yoa

+ 1)1,2}1,

\4d2c~

dl 2d2 Co

12d2co

-- (~2 -- (~4 "{'i7~)1/2} 1/3

--dl

dl +Id YO =

--dI 2(Xn)ea 2d2Co

t U ( d l , d2, Co ,) I G(icdc) / ~

x/d2 x2 K1/2

COS7 =

-

(1 "t-cos "{} 1/2

Co yK1/2

1 P(x)

1

-IG(ic~c)l = G(icOc) < l i r a

27d~ +

\4d.2 %2

27d~/

T h u s t h e m a x i m u m value o f U is 2. This means t h a t Rop t is t h e same as f o u n d for p=l. R o p t = 2 { c o s (Tr/n)} n T h e r e are no oscillations f o r n ~< 7, b u t it is possible t o c o n s t r u c t a permissible p a r a m e t e r

98 set giving oscillations for p = 2 if n/> 8. These p r e d i c t i o n s have been c o n f i r m e d b y n u m e r i c a l studies. It should be n o t e d t h a t Walter's extensive set o f simulations (1970, 1 9 7 1 ) failed to d e t e c t a n y oscillations for p = 2. Given the results p r e s e n t e d here it is possible t o conclude t h a t had a higher d i m e n s i o n been t e s t e d or a slightly larger value o f k used, t h e n t h e oscillations w o u l d have been d e t e c t e d . This c o n s t i t u t e s an e x c e l l e n t e x a m p l e o f the p o w e r o f an analytic t e c h n i q u e , even an a p p r o x i m a t e one, over a trial and e r r o r search t h r o u g h a large p a r a m e t e r space b y n u m e r i c a l simulations. We n o w s u m m a r i z e the results t h a t have been o b t a i n e d f o r p = 1 and p = 2. (1) F o r all p, for all n and for all permissible p a r a m e t e r sets, if t h e r e is an oscillation in positive space t h e n its a m p l i t u d e is less t h a n x*, where x * is established u n i q u e l y b y the following equation:

0 = 1

2X'Co

1 + d 2 x * ( 1 + cos 0) p d0

(2) F o r p = 1, for n and f o r all permissible p a r a m e t e r sets there are no stable limit cycles such t h a t xj(t) is always positive. (3) F o r p = 1, f o r all n less t h a n or equal t o seven and f o r all permissible p a r a m e t e r sets, t h e r e are n o limit cycles in positive n-space e i t h e r stable or unstable. (4) F o r p = 1, and for all n greater t h a n or equal to eight, it is possible to cons t r u c t a permissible p a r a m e t e r set which yields a s y s t e m possessing an unstable limit cycle. The limit cycle will exist if R > 1, where R = 2 IG(i~c )/G(0){ (5) For p -- 2, for all n less t h a n or equal t o seven and f o r all permissible p a r a m e t e r sets, t h e r e are n o limit cycles in positive n-space e i t h e r stable or unstable. (6) F o r p = 2, for all n greater t h a n or equal to eight it is possible t o c o n s t r u c t a permissible p a r a m e t e r set which yields a s y s t e m possessing a stable limit cycle where xj(t) > 0 for all t t> 0. A stable limit cycle will be p r e s e n t if R is greater t h a n one.

A = d I x/~2/x/~c 0 G(icoc) R =

X 2--

~-~ +

~-

27A6 ]

j

8 'i "27"3/ j

J

Defining (Xn)e as the equilibrium value o f xn,R becomes: R - 2 d 2 c ° (x~) 3 IG(icoc)/G(0)l dl If bi = b > 0 for all j t h e n R b e c o m e s : R -

2d2 Co dl (xn)3 {cos (zr/n) }n

4. Conclusions and Discussion It has been seen t h a t the distinction b e t w e e n an oscillating and a non-oscillating system rests principally with the s t o i c h i o m e t r i c c o e f f i c i e n t p. This means t h a t the b e h a v i o r is d e t e r m i n e d by t h e kinetic s t r u c t u r e o f the chemical netw o r k r a t h e r t h a n exclusively b y t h e n u m e r i c a l values o f the r e a c t i o n constants. G o o d w i n ( 1 9 6 5 ) originally p r e s e n t e d these e q u a t i o n s as a m o d e l f o r the c o n t r o l o f p r o t e i n synthesis. In t h a t m o d e l p r e p r e s e n t e d t h e n u m b e r of molecules o f i n h i b i t o r t h a t w o u l d have to c o m b i n e with t h e o p e r a t o r to p r e v e n t transcription. If t h a t m o d e l is valid t h e n the p r e s e n t results indicate t h a t t h e specification o f the oscillatory p r o p e r t i e s o f a p r o t e i n synthesis c o n t r o l l o o p is c o n t a i n e d in t h e genetic material o f the operator. E x p e r i m e n t a l l y one observes t h a t p = 1 is the m o s t c o m m o n value o f p. This w o u l d m e a n t h a t oscillatory b e h a v i o r w o u l d be an exceptional m o d e o f behavior. This observation is c o n s i s t e n t with a c o n c e p t u a l m o d e l o f metabolic regulation in which a small n u m b e r o f aut o n o m o u s oscillators c o o r d i n a t e a large n u m -

99 ber of reactions by coupling through c o m m o n reactants. It is possible to determine the higher harmonic c o n t e n t in the resulting oscillations. This c o n t e n t is found to be very low especially in a biologically plausible system which contains a large number of variables. This means that the oscillations can be approximated by a simple trigonometric function for most purposes. Accordingly, when considering the effect of coupling this oscillator to another pathway, it is possible to make use of the large body of mathematical results concerning sinusoidally forced nonlinear systems. An example of such a coupled system is one in which the present oscillator serves as an amplitude dependent state-transition trigger. The qualitative behavior of the envelope of the oscillation during 1;he run up to the stable oscillation varies with initial conditions. In some cases the osciUation increases slowly and monotonically to the final limit cycle, while for a different set of initial values the system may very quickly spike to a large amplitude and then decay to the same stable oscillation pattern. Clearly the triggering properties of the two cases are very different. Since most models of developing systeras are built on the hypothesis that differences in position are expressed as different concentration values, this special property of nonlinear oscillators may be of importance. Well established mathematical tools (Gelb and van der Velde, 1968; Voronov, 1957) for examining the run up to a stable oscillation exist and can be applied to this problem. The describing function can give us information about the minimum number of reactants that can establish an oscillating reaction network. This question has been asked in connection with prebiotic ,;ystems. As we have already pointed out, since it is incorrect to attribute sophisticated kinetic mechanisms to prebiotic systems, the simple feedback loop examined here may have been one of the first cyclic regulatory structures to appear. For the most general form of the equations as given in equation 1, it is seen that the smallest system

which can oscillate is the three dimensional one, since that is the smallest dimension which can possess a nonzero value for ~c- If time lags between the reactions are considered however, smaller systems can have a nontrivial ¢Oc. The describing function can easily be extended to include systems which contain delays between reactions since the Laplace transform of the time delay operator is particularly simple. Such delay components are essential to the accurate modeling of several important biological processes; notably those of macromolecule synthesis. In general the presence of these time delays tends to increase the value of R, and thus makes oscillations more likely. The nonlinear function f need not be continuous for the describing function m e t h o d to be applicable. It only needs to be integrable. It is also possible to extend the m e t h o d to cover cases where f is not a function but rather a closed hysteresis loop. For an f of this form the integral a2 is not zero, but instead it is the area enclosed by the loop. The describing function becomes complex, but the m e t h o d carries over with no further variation. Several examples of hysteresis type kinetic systems have been considered in the chemical literature (Selkov, 1974; Bunov, 1974; Tanner, 1 9 7 4 ) r e c e n t l y . The describing function m e t h o d is the only approach to theoretically studying these systems in large dimensions except for simulation techniques. Recently it has been suggested, on the basis of experimental evidence (Hess, 1974) that biochemical control circuits when subjected to a randomly varying input may begin to oscillate. This effect is c o m m o n l y observed in control systems. It is possible to derive describing functions for nonlinear systems with random inputs. Several tables of these functions have been compiled (Gelb and van der Velde, 1968). The m e t h o d provides a means of analytically predicting a control system's response to chemical noise. It is possible to extend the m e t h o d to cover multiple loop systems. Since a control circuit of this type will have a frequency dependent describing function, the graphical analysis pro-

100

cedure used here is no longer applicable. Alternative procedures have been derived. The method is being applied by the author to a system in which both transcription and translation are inhibited by a feedback metabolite. In summary we wish to suggest that several problems encountered in the theoretical study of metabolic regulation, that cannot be resolved by linear control theory or simulation techniques, can be solved by the methods of nonlinear control theory.

References Bunow, B., 1974, Paper given at the 9th FEBS Meeting, Budapest,. Gelb, A., W. van der Velde, 1968, Multiple Input Describing Functions and Nonlinear System Design (McGraw-Hill Book Company, New York). Goodwin, B., 1965, in Advances in Enzyme Regulation Weber, Ed., (Academic Press, New York).

Griffith, J., 1968a, J. Theoret. Biol., 20,202. Griffith, J., 1968b, J. Theoret. Biol., 20, 209. Hess, B., 1974, Paper given at the 9th FEBS Meeting, Budapest. Hsu, J. and A. Meyer, 1968, Modern Control Principles and Applications, (McGraw-Hill Book Company, New York). Morales, M. and D. McKay, 1967, Biophysical J., 7, 621. Rapp, P., 1974, J. Mathematical BioSciences, 23, 289. Selkov, E., 1974, Paper given at the 9th FEBS Meeting, Budapest. Tanner, R., 1974, Chem. Eng. Sci., 29, 169. Tiwari, J. and A. Fraser, 1973, J. Theoret. Biol., 39, 679. Tiwari, J., A. Fraser and B. Beckman, 1974, J. Theoret. Biol., 45, 311. Voronov, A., 1957, Automatika i Telemekhanika, 18, 253. Walter, C., 1969a, J. Theoret. Biol., 23, 23. Walter, C., 1969b, J. Theoret. Biol., 23, 39. Walter, C., 1970, J. Theoret. Biol., 27,259. Walter, C., 1971, in: Biochemical Regulatory Mechanisms in Eukaryotic Cells, Kun and Grisolia, Eds., (Wiley Interscience, New York).