Limit Cycles in the Plane An Equivalence Class of Homogeneous Systems BY JOHN TEXTER Department of Chemistry and Center for Surface and Coatings Research Lehigh University Bethlehem Pennsylvania 1801 5 U.S.A. Receiued 30th July 1974 A set theoretic presentation of the necessary and sufficient conditions for the Occurrence of limit cycle behaviour in open homogeneous systems is given for the specific case of two time-dependent components. The results are obtained by using several theorems due to Zubov in the framework of Lyapunov's stability theory. The starting point in the analysis is the definition of a closed invariant set which depicts the hypothesized or experimentally observed limit cycle. Classes of kinetics which evolve to the defined limit cycle may then be generated.When kinetics of a particular polynomic form are presumed necessary conditions on the rate constants and external constraints are generated. Time reversal generates an equivalence class of kinetics which has the previously defined limit cycle set as an asymptotic stability domain boundary. This concept may be of use in the consideration of the turning on and off of biological clocks. Introductory results regarding the admissible kinetics giving rise to elliptical limit cycles are presented and the use of the theory in treating experimentally observed limit cycles is discussed. 1. INTRODUCTION Renewed and increasing interest has over the past decade been shown in the study of physicochemical systems which give rise to sustained oscillatory behaviour.The number of such systems examined experimentally seems to be steadily increas- ing.1-4 In formulating a firm physical basis for such oscillatory phenomena it is required that the observed behaviour be reproduced quantitatively in terms of descriptive mathematical models. Attention herein will be focussed on open homo- geneous systems with time-independent boundary conditions (external constraints) and two time-dependent components. The dynamics of such systems may be suitably modelled in terms of two coupled nonlinear ordinary differential equations. Since the sustained oscillations which have been reported in the literature are usually preceded by an induction period attention herein is further restricted to oscillations of the limit cycle type.It should be noted that many hypothetical chemical mechanisms have been formu- lated which (as proven by numerical integration procedures) admit limit cycle behavio~r.~-' Although numerical integration of the rate equations for a proposed mechanism may certainly be sufficient to prove that a given mechanism admits limit cycle behaviour such models of real systems are often stiff.14-15 In these circum- stances numerical integration becomes prohibitive. The problem at hand then is to develop necessary and sufficient criteria (exclusive of numerical integration) by which it may be judged whether or not a given model will admit limit cycle oscillations. The ultimate test of any proposed model remains one of numerical simulation however. Several necessary conditions for chemical limit cycle oscillations have been proposed Physically these include the requirement of autocatalytic or crosscatalytic interactions among the chemical species.These requirements follow naturally from 254 JOHN TEXTER the necessary condition that such systems be described by nonlinear ordinary differential equations. Another necessary condition usually considered is the requirement that the linear part of a system’s rate equations (as obtained by an expansion of the original rate expressions about the critical point or steady state whose stability character is being investigated) have eigenvalues with positive real parts. This guarantees that the critical point (any limit cycle surrounds at least one critical point) will be repulsive to small perturbations.This convention will be adopted herein in that the describing differential equations (in perturbed coordinates) associated with a given limit cycle will be presumed to have a linear part which is positive definite. The important point here is that the critical point surrounded by the limit cycle (for the sake of discussion cases involving multiple critical points surrounded by a limit cycle will be excluded herein) must be repulsive with respect to infinitesimal perturbations. To the best of this writer’s knowledge no firm physical basis has so far been established to exclude the possibility of limit cycles occurring for systems whose linear part is negative or positive semidefinite. The theory of indices provides another (qualitative) necessary condition for the occurrence of limit cycles and Bendixson’s negative criterion gives sufficient conditions for the non- existence of limit cycles.The Poincare-Bendixson theorem gives both necessary and sufficient conditions for the occurrence of limit cycles but these criteria are not easily implemented. In this paper it will be attempted to implement a more quantitative theory of limit cycle oscillations than has heretofore been utilized in considering physico- chemical oscillations. Two general questions (not unrelated to one another) are investigated. First does a given pair of coupled rate expressions admit a limit cycle as an asymptotic solution? Secondly what is the general class of kinetics which will give rise to an a priori defined limit cycle solution? The basic theory utilized herein is that of Zubov 7* ’ in the framework of Lyapunov’s theory of stability of motion.Both necessary and sufficient conditions for limit cycles are included. A set theoretic presentation of Zubov’s theory is given in sections 2 and 3. In section 2 the general class 8,of kinetics is constructed which has an a priori defined asymptotic stability domain. Usage of Zubov’s theory in this context has been noted previously 2o and has successfully been applied to the characterization of bistability in two-mode lasers 21 and the characterization of thermal explosion limits.22 In section 3 the construction of the general class of kinetics having a given a priori defined limit cycle solution is discussed.Various applications of the theory are illustrated and discussed in section 4. The results reported herein are only to be considered introductory and do not provide means for complete system identification ;however sufficient criteria for constructing physicochemical limit cycle oscillations are obtained. 2. CONSTRUCTION OF gS In this section the construction of the general class 8,of two-dimensional open homogeneous reaction kinetics X = Fl(X Y) P = F2(X Y) is considered where the system (2.1) admits a finite domain of asymptotic stability about a critical point (Xc Yc)contained in this domain. For the sake of what follows the righthand sides of (2.1) are assumed to be holomorphic (i.e. expressible as a convergent Taylor series about any point) in some domain embedding the domain of asymptotic stability.To facilitate discussion the following definitions are introduced. 256 LIMIT CYCLES IN THE PLANE Definition 1. The coordinates (X Y) are restricted to some open (simply connected) set GI of the real euclidean plane according to their type. For example if X and Yare both chemical concentrations or activity variables CI is the positive quadrant of the plane. DeJinition 2. The critical point (X, Y,) comprises an invariant set MS,ASc C! and JlS= {(X, Y,) I Ff(Xc,Y,) = 0 i = 1,2). Asis assumed to be unique in its domain of attraction d,. DeJinition 3. The domain of attraction of A%’, is an open bounded and simply connected set d where dScR and &isc ds.Asis further assumed to be an interior point of ds. Definition 4. The boundary of the asymptotic stability domain d,in a closed invariant set Al c CI which is a limit set of dS.Unless noted otherwise constitutes a closed curve in Q. DeJinition 5. The set d1is a closed set where dl= dsuAl. Note that every neighbourhood of every point in Atl contains points of ds. Definition 6. The set of a priori admissibility constraints is denoted by G = {Gi i = 1,2 . . .>where the Giare related to the parameters which appear in the righthand sides of (2.1). For example if X and Y denote chemical concentrations certain of the Gi would ensure the positivity of X and Y. Others of the Gfmay reflect what is already known about or presumed for the system (2.1). For example if the right- hand sides of (2.1) are presumed to have a given polynomic form certain of the Gi may denote values of certain of the rate constants or pseudo rate constants and some of the G1may denote the sign appearing in front of certain of the terms accordingly as to whether or not restrictions are imposed on the molecularities of the interactions between Xand Y.The transformation is now made to perturbed coordinates (x,y) = (X- X, Y-Y,) so that (2.1) becomes 2 =fdX7 Y> (2.2) P =fZ(T Y) wheref,(O 0) = 0 i = 1,2. The linear part of (2.2) will be assumed to be negative definite. As this exposition is fundamentally based on Zubov’s results his principal theorem for the cases considered herein is stated as follows Theorem. The system (2.2) has dsas its domain of asymptotic stability.A Lyapunov function v(x,y) exists for (2.2) such that for all (x,y) E ds,0 < ZJ < 1 and for all (x,y) EA’l,ZJ = 1. A continuous positive definite function 4(x,y) exists for all (x y) ES2 such that v is defined by the condition v(0,O) = 0 and the equation au at, fi -+ fi -= -q5(1-v)(1+ f;+f;)+. ax ay Proofs may be found in Zubov,17* l8 Hahn,23or in Bhatia and Szego 24 where d,is an open invariant set. The construction now proceeds by defining All through a choice of u such that v(x,y) = 1 for all (x,y) E4,and such that MIdescribes a closed curve. 8,is thus obtained as the solution (fl,fz) to the inhomogeneous system where q5 (satisfies the hypothesis of the Theorem) a b and c (all scalar functions of x and y or constants) are selected in conformity with the desired spectral properties required of the linear part of (2.2) and such that the determinant of the coefficient matrix of the lefthand side of (2.4) may not vanish in d,except possibly at the origin JOHN TEXTER ('in which case the righthand side of (2.4) vanishes identically at the origin also).The a b and c may be restricted in such a manner that the solutions have a particular polynomic or nonpolynomic form (i.e. they must be chosen in conformity with the admissibility constraints G). As will be subsequently illustrated the solution of (2.4) is simplified somewhat if the substitution 4 = 4(1 +ff+f;)-* is made. It should also be noted that the set A may or may not (depending on the choices of a b c and G) contain additional critical points of the system (2.2) obtained from (2.4).To illustrate the flexibility of the method the class b of kinetics giving rise to global asymptotic stability is constructed. The construction proceeds analogously as that for Q except that d,= ICZ and the solutions are obtained from the inhomogeneous system rather than from (2.4). In addition u must be chosen positive definite such that v(x,y)-+ GO as x2+y2 -+ GO. The extension of the construction to higher dimen- sional systems is straightforward. Both polynomial and nonpolynomial kinetics are generated depending on the selection of the functions 21 4,a band c. Nazarea 25-27 has illustrated a different generating function method far obtaining these results.3. CONSTRUCTION OF 81 In this section the construction of the general class Q of admissible kinetics having an invariant set A1as an asymptotic stability boundary is outlined. Proceeding as before in section 2 the class gSis constructed under the additional constraint that A contain no critical points of the solutionsf andf of (2.4) i.e. .fi andf do not vanish identically for all (x y)E Ar. Under the previously presumed condition that the righthand side of (2.1) (and hence of (2.2)) is holomorphic in some domain embedding dl,the continuation of integral curves from points in d,to the negative semi-trajectory tE (0 -a),is assured as is their boundedness. As a practical consideration it may additionally be required that there exist some sufficiently large open set $4' embedding d1such that there exist no critical points of (2.2) in the set 9 where 9 = %-&',.This can be an especially important consideration in systems admitting multiple steady states. The class is then obtained by the substitution of t + -1 in the solutions (2.2) obtained from (2.4). The added restrictions of this section are where necessary implemented in the selection of the functions 4,a band c and the admissibility constraint G. Under such considerations as above gSand g1are seen to be equivalent classes under time reversal. 4. APPLICATIONS AND EXAMPLES The procedures outlined in sections 2 and 3 have several foreseeable applications in the identification and analysis of temporal ordering in nonlinear open systems.The generation of the classes cTs and 8,is a very general procedure and therefore highly flexible. The flexibility is inherent in the selection of the functions v 4 a b and c (cf. (2.4)) which are chosen in accordance with various constraints. Since Zubov's theory provides necessary and sufficient stability criteria it may be used as a double-edged sword in practical applications. As the ensuing discussions will indicate these procedures provide us with a powerful analytical tool. Implementa-tion of the theory is not devoid of difficulties and restrictions however. LIMIT CYCLES IN THE PLANE ANALYSIS OF STIFF SYSTEMS In cases where a proposed model of a limit cycle is sufficiently stiff to make numerical dynamical simulation prohibitive the following procedure may be used to establish the existence of a limit cycle.If the limit cycle exists its phase portrait is also obtained. The substitution t + -t is first made in the differential equations suspected of admitting a limit cycle. The equations are then inserted into (2.3)where 4(1+f +fi)3 is replaced by a quadratic form in conformity with the linear parts of fi andf,. The solution u(x,y) is then obtained necessarily by the method of un-determined coefficients as a power series of homogeneous forms until convergence is achieved. Numerical details of this method have been given elsewhere.28*29 After convergence has been obtained the indicated boundary set Aflis checked to see if it contains any critical points of the systemf andf,. If there are no critical points in the existence of the cycle is proven.Unfortunately the methods of section 3 are restricted to two-dimensional systems so that higher-dimensional limit cycles cannot be constructed in the manner for two- dimensional cases. However three-dimensional cases may be treated in the following manner but only necessary conditions are generated. The three-dimensional limit cycle model is first transformed as before (t -,-t). One of the coordinates is treated as a parameter and two-dimensional cross-sections of the asymptotic stability domain boundary are generated in a step-wise manner proceeding along both the positive and negative semiaxes of the parameterized coordinate. The three-dimensional stability surface is thus generated sequentially.If the original three-dimensional system considered admits a limit cycle it will necessarily describe a closed curve on the three dimensional asymptotic stability surface. The locus of points observed experimen- tally may be compared with the numerically generated stability surface. TIME REVERSAL AND BIOLOGICAL CLOCKS The mapping of time reversal has been illustrated to generate an equivalence class of kinetics. This concept may be of some use in developing a fuller appreciation of biological clocks and rhythms which exhibit intermittent periodic behaviour. For the sake of discussion consider the equivalence class of kinetics which admits the invariant set A’ illustrated in fig. 1 respectively as a domain of asymptotic stability and as a limit cycle.The subspaces k and k’ denote the external constraints for the respective situations. The desired time reversal may thus be realized as a invertible transformation T operating on the external constraints where $k = k’ and T-’k’ = k. The case illustrated in fig. 1 is highly idealized. From a practical standpoint the respective invariant sets Asneed not coincide with each other and similarly the invariant sets JA%‘~ may also be disjoint. A physical example of a system which admits this kind of switching is the Teorell membrane oscillator 30-33 which exhibits damped oscillatory behaviour below a given current threshold and sustained oscillations above the threshold. SYSTEM IDENTIFICATION The present theory should prove to be useful in furthering our understanding of the dynamics of physicochemical systems exhibiting limit cycle behaviour.It will complement the more physical methods used in studying oscillations. Its ultimate JOHN TEXTER 259 utility will only be realized when considered in conjunctjon with exhaustive experi- mental data. The reason for this is that the theory focusses on invariant sets. A single application of the theory to an isolated data set describing an observed limit cycle would probably be of marginal utility as far as system identification is concerned. This is because an infinity of kinetics may exist which admits a given invariant set A1as a limit cycle. In other words induction period behaviour is not considered explicitly in the theory although certain implicit constraints may be included in the construction as indicated in section 2.I-Iy Y FIG.1.-Schematic representationof the concept of time reversal in the turning on and offof biological clocks. The small open circles depict the critical point associated with the illustrated invariant sets A]. In the lefthand figure A1is an asymptotic stability domain boundary and in the right figure .MI is a limit cycle. The coordinates k and k’ denote the external constraints associated with the respec- tive situations. Best use of the theory is made when an appreciable amount of information is already known (or suspected) about the physical system under investigation. Thus if kinetics of a given polynomic or nonpolynomic form are suspected necessary conditions on the rate constants proportionality coefficients and external constraints may be generated.These considerations then provide a severe test of admissibility that must be satisfied by any finally accepted model for a given system. EXAMPLE CONSTRUCTIONS To illustrate the methods of section 3 a class of kinetics which admits the unit circle as a limit cycle is constructed. The invariant set A = ((x y) 1 u(x y) = x2+y2 = 11 is illustrated in fig. 2 in arbitrary nonperturbed coordinates. Poly-nomial kinetics of the form 2 = alx+a2y+a3x2 +a4xy+a,y2 +a6x3 +a7x2y+a8xy2 +agy3 (4.2) j = blx + b,y+ b3x2+ b.Q~y+ bsy2+ b6x3+ b7x2y+ b8Xy2 + b9y3 are presumed where the variables x and y are perturbed coordinates and the coefficients at bi are to be determined.A particular solution is obtained by solving (2.4) with +(l +f +f,”>* = dlx2+d2y2(dl,d2 > 0) a = 0,b = 1 and c = 3. The result is obtained as three inhomogeneous equations in the unknowns a, b, a2 and b2,four homogeneous equations in the unknowns a3 b3 a, b, a and b5,and five inhomo-geneous equations in the unknowns as b6 a7 b7 as b8 ag and b9. As all three sets of equations are undetermined an infinity of solutions exists. Unfortunately these solutions will not necessarily admit limit cycle behaviour following time reversal as LlMIT CYCLES IN THE PLANE 2 Y I 0 0 I 2 3 X FIG.2.-The unit circle as a limit cycle. Coordinates X and Y are arbitrary. The small open circle represents the origin of the perturbed coordinates (x y) of (4.3).The solid lines approaching the cycle illustrate the attraction of the cycle for the indicated initial conditions and were obtained by numerical integration of (4.3)with a1 = b1 = 0.5 and a2 = -b2 = 0.1. The dashed curve illustrates the type of behaviour that may be obtained when critical points are not explicitly excluded from the proposed limit cycle set. The dashed curve was obtained by numerical integration of the system (4.2) with a1 = -a6 = 1.0 b2 = -b9 = 0.5 a8 = 1.5 and the remaining coefficients set equal to zero and is seen to come to rest at the critical point “X” at the top of the circle. I 0 -0.2 -0.1 0.0 Ql a2 V FIG.3.-Qualitative illustration of the treatment of experimental data reported for the Teorell membrane oscillator.E is the potential in V and V is the net volume rate in ml min-’. Curves a and b were defined by making a least squares fit of the function ae2 +yv2 = 1to the experimental data (open circles) where e = E-E and v = Vfor two different choices of & (a & = 1.85 V ; b Ec = 2.25 V). The curves a and b were generated by numerical integration of the system (4.4) with dl = d2 = 1.0 I = 1.0 (ar,~) = (0.545 35.3) for a and (a,y) = (0.639 27.7) for b. JOHN TEXTER critical points have not necessarily been excluded from the boundary. In fact ten arbitrarily initiated distinct solutions were found to have critical points on the boundary. To obviate this difficulty the following less general kinetics were considered i= alx(l-x2-y2)+a2y 3 = bly(1-x2-y2) +b2~.(4.3) The general solution (d = d2 = 1) is easily demonstrated to be given as al = bl = 0.5 and a = -b # 0. The form of (4.3) was selected so as to exclude critical points from the unit circle. In view of the fact that the determinant 0,5(1 -x2-y2) I b2 0.5(1 -x2-y2) I -b2 is positive for all (x,y) and finite b2 it is assured that the origin is the only critical point admitted by (4.3) so that the limit cycle is globally attractive. It should be noted that (4.3) is not admissible as a purely chemical model if interactions among the species are limited to trimolecular and lower order interactions (i.e. certain of the admissibility constraints G1for chemical reaction kinetics are never satisfied by (4.3)). The results for (4.3) may immediately be extended to the case where ~(x, y) = orx2+yy2 = 1 describes the general ellipse in standard position and $(l +f +Sf)) is chosen as the positive definite quadratic form dlx2+d2y2.The general solution is obtained as (4.4) d2 3 = -y( 1-ax2-yy2)-zx 2Y where I is finite. For illustrative purposes the results (4.4) were generated (d = dz = 1) for two ellipses obtained by a least squares fit to the limit cycle (in the potential-net volume rate plane) reported by Teorell 34 for his membrane oscillator. The results are depicted in fig. 3 where the open circles are the data to which the ellipses were fitted. The limit cycles a and b were defined by arbitrarily presuming the location of the critical point to be given by (E, V,) = (1.85 0) and (2.25 0) respectively and then carrying out the least squares fit in perturbed coordinates (e v) where e = E-E and v = V.The cycles a and b are respectively defined by the relations u(e v) = 0.545e2 +35.3v2 = 1 and v (e v) = 0.639e2+27.7v2 = 1 ; the illustrated curves were generated by numerical integration of (4.4). The purpose here is not to confirm or dispute the empirical model proposed by Teorell for this system as the experimental data do not depict an ellipse but rather to stress the importance of having at least limited a priuri knowledge of the model being investigated. For example what is the actual location of the critical point encircled by the experimental cycle? In this instance a good initial guess might be obtained by extrapolation of the observed sequence of critical points observed in the relaxation of the system at various driving currents below threshold.5. CONCLUSIONS Necessary and sufficient conditions for the admissibility of limit cycle oscillations in physicochemical systems have been presented. Hypothetical reaction kinetics giving rise to an a priori limit cycle solution are constructed for the first time in a LIMlT CYCLES IN THE PLANE deterministic non-trial-and-error fashion and use of the theory in treating experi- mental results has been outlined. Extensions of the current results are currently being attempted with the hope of determining whether or not purely chemical kinetics (with the accompanying constraints G) will admit a quadratic form (as well as higher order forms) as a limit cycle.The author would like to thank Prof. G. Stenglefor several criticisms and suggestions concerning this paper. B. Has and A. Boiteux Ann. Rev. Biochent. 1971 40,237. G. Nicolis and J. Portnow Chem. Rev. 1973 73 365. B. Chance E. K. Pye A. K. Ghosh and B. Hess (eds.) Biological and Biochemical Oscillators (Academic Press New York 1973). D. E. F. Harrison and H. H. Topiwala Adv. Biochem. Eng. 1974 3 167. I. Prigogine and R. Lefever J. Chem. Phys. 1968,48,1695. R. Lefever J. Chem. Phys. 1968 49,4977. ’R.Lefever and G. Nicolis J. Theor. Biol. 1971 30 267. G. Nicolis Adv. Chem. Phys. 1971 19 209. B. Lavenda G. Nicolis and M. Herschkowitz-Kaufman,J. Theor. Biol. 1971,32,283. lo M. Herschkowitz-Kaufman and G.Nicolis J. Chem. Phys. 1972 56 1890. l1 J. J. Tyson J. Chem. Phys. 1973,58,3919. l2 J. J. Tyson and J. C. Light J. Chem. Phys. 1973,59,4164. l3 M. L. Smoes and J. Dreitlin J. Chem. Phys. 1973,59 6277. l4 R. J. Field and R. M. Noyes J. Chem. Phys. 1974 60,1877. l5 E. M. Chance ref. (4) p. 177. l6 V. V. Nemytskii and V. V. Stepanov Qualitative Theory of Diferential Equations (Moscow 1949 English translation Princeton University Princeton 1960). V. I. Zubov Prikl. Mat Mekh. 1955 19 179. V. I. Zubov Methods of A. M. Lyapunov and Their Application (Noordhoff Groningen Netherlands 1964). l9 A. M. Lyapunov General Problem of Stability of Motion Comm. SOC.Math. Kharkov 1892 (Probkme Generale de la Stabilitk de Mouvement Ann. Math. Studies 17 Princeton 1949).2o J. Texter J. Chem. Phys. 1973 58,4025. 21 J. Texter and E. E. Bergmann Phys. Rev. A 1974 9 2649. 22 D. Dwyer and J. Texter unpublished. 23 W. Hahn Stability of Motion (Springer-Verlag Berlin 1967). 24 N. P. Bhatia and G. P. Szego Stability Theory of Dynamical Systems (Springer-Verlag Berlin 1970). 2s A. D. Nazarea Biophys. 1971 7 85. 26 A. D. Nazarea Biophys. 1972 8 96. 27 A. D. Nazarea Biophys. 1973 9 93. 28 S. G. Margolis and W. G. Vogt IEEE Trans. Automat. Contr. 1963 AC-8 104. 2g J. J. Rodden in 5th Joint Automatic Control Conf. (Stanford Electronic Laboratories Stanford California 1969) p. 261. 30 T. Teorell Disc. Faraday SOC.,1956 21 9. 31 T. Teorell J. Gen. Physiol. 1959 42 831. 32 T. Teorell J. Gen. Physiol. 1959 42 847. 33 T. Teorell Biophys. J. 1962 2 27. 35 See ref. (31) fig. 3.