摘要:

Unified Theory of Temporal and Spatial Instabilities Institute of Physics AND H. DEGN Institute of Biochemistry Odense University Denmark Receiced 29th July 1974 Chemical reaction-diffusion systems are discussed with the view of clarifying the relation between spatial and temporal concentration variations. The predicted absence of oscillating or propagating structures in closed two component systems which are not oscillatory when stirred points to the necessity of studying systems with more than two components. 1. INTRODUCTION The occurrence of concentration oscillations in homogeneous chemical systems was a matter of controversy for several decades until the sixties when there came a boom in experimental and theoretical work on the subject. The existence of homogeneous chemical oscillations then became firmly established.This had hardly happened before the logical extension of temporal periodicity namely spatial periodicity took over the role of the flying saucers of chemistry. The possibility of the emergence of a spatial concentration pattern in an initially homogeneous chemical system seems to have first occurred to Turing.' He was led to the idea when searching for an explanation of the stage in morphogenesis when a small spherical array of identical cells assigns different parts of its surface to develop into different organs. Turing's work published in 1952 only became famous after 1967 when the Brussels school began the publication of a series of papers on spatial instability in homogeneous chemical systems.The report by Zhabotinskii on spatial phenomena in the Belousov reaction at the Prague Symposium on Oscillating Reactions in 1968 inspired experimental work on spatial patterns in many labora- tories. Despite the parallel courses of the theoretical and experimental approaches it is not always evident that the theoreticians and the experimentalists talk about the same thing. The pretentious title of this paper reflects the fact that temporal and spatial instab- ilities are closely related phenomena as they both require reaction schemes with feed- back. The aim of the paper is to clarify this relationship and to give operational interpretations of theoretical results regarding spatial instabilities and spatial period- icity. The concept of concentration oscillations in a homogeneous phase is so straight-forward that it needs no special introduction.Spatial periodicity on the other hand is an equivccal expression covering different phenomena which may be easily confused. Before we proceed to theoretical analyses we shall therefore classify these phenomena in operational terms. 233 TEMPORAL AND SPATIAL INSTABILITIES Imagine a solution in an ordinary vessel made of non-reactive material. We stir the solution and it does not exhibit time periodicity as long as it is stirred. When we stop the stirring conceivably we may observe that concentration gradients form spontaneously in the solution. These gradients may tend to form a time-invariant structure or they may be periodic in time either as propagating or as oscillating standing concentration waves.A reaction solution which forms spatial structures after it was homogenized by stirring and did not exhibit concentration oscillations whilc stirred is said to be spatially unstable. Current theory predicts spatial instability in reaction systems which comprise diffusion of chemical components and certain types of chemical reactions with feed-back loops. In the discussion below systems which produce time-invariant structures are said to have aperiodic spatial instability and systems giving rise to coilcentration waves are said to have oscillatory spatial instability. Let us now return to an imaginary experiment where a solution in an ordinary vessel with non-reacting walls exhibits concentration oscillations when stirred.When we stop the stirring conceivably we may observe the emergence of concentration structures which are oscillating with the same period as the previous homogeneous oscillation. In this case spatial instability is not involved. Time-dependent structur:can be explained by small local heterogeneities or externally induced gradients causing different parts of the oscillating solution to come out of phase. Such patterns in oscillating systems will tend to be destroyed by diffusion of the reaction components of the chemical system. This paper is concerned with structures and waves which result from spatial instability and thus grow and are maintained with the assistance of diffusion. 2. TWO-COMPONENT SYSTEMS We consider a system in which two substances A and B are present with concen- trations A(r t) and B(r t) as function of space P and time t.The relevant reactions are assumed to be fully described by two rate equations 8AS -= fA(A B)+DAV2A at dB -= fn(A B)+DBV2B at wheref andf describe the chemical reaction rates DA and DBare diffusion cozfTi-cients and V is the gradient operator. We shall first neglect the spatial dependence of A and B by assuming stirring or large diffusion coefficients. As shown in several works a linear theory provides important information on the dynamics near singular points (As,B,) in the A-B phase plane where For sniall deviations from the singular point a(t)= A(t)-A b(r) = B(t)-R one obtains da -= mAAa+mABb dt I.BALSLEV AND H. DEGN 235 171 are partial derivatives where the coefficients a f (3) in the singular point. Solations to eqn (2) may be written a = a. exp((,u-iio)t) b = bo exp((p-ico)t) (4) where p is the growth rate and co is the angular frequency. The complex growth rate 3. = ji-ico is then determined as the roots of the secular equation (?nAA -l)(mBB -?L) -mABn1BA = 0. (5) Let us denote these roots by A1 = pl+io and l2= p2-io. For LO + 0 we have pl = p2 and we shall in this case denote the singular point as oscillatory unstable if ,ul > 0 and oscillatory stabZe if p1 < 0. For CD = 0 the singular point is aperiodically stabZe if both ,ul and p2 are negative. Two positive values correspond to aperiodical instability and different signs lead to a saddle type instability.In the later discussion of spatial effects this notation is slightly more suitable than the conventional terms focus node and saddle. Asit turns out that the case mABmBA 2 0 is uninteresting in connection with spatial effects included later we shall assume mABmBA < 0. In a qualitative discussion it is convenient to scale the time unit and the unit for the concentrations a and b. In such rationalized units the complex growth rate is obtained from solutions to (.-A)(p-n)+l = 0 (6) where The solutions of eqn (6) are conveniently represented in the ap parameter space with domains for different types of dynamics according to the values of and &. The domains are shown in the ap diagram in fig.1. This graphical representation is useful in the analysis of specific chemical reaction schemes. For fixed rate constants the scheme defines a single point (a,p). For example in the reaction scheme of Lotka (1910) the points obtained for different values of rate constants fall on the negative p-axis and the Volterra-Lotka scheme has always (E p) = (0,O). In the reaction scheme of Prigogine and Lefever it is possible to obtain points in the region between the hyperbola branches and the axes. This latter scheme and reaction systems involving forward inhibition or an enzymatic step were analyzed previously by transforming the ap diagram into a parameter space directly related to the chemical rate constants. Including the spatial dependence of A and B of eqn (1) and assuming a spatially periodic solutioii we find in the linear theory a(r t) = a exp(,ut -iot + ik.r) b(r,t) = boexp(pt -iot + ikx) where k is thc wave vector defining a propagation direction and a wave length of 2x/(kl.Introducing a scaled unit for length we obtain the secular equation (u-DAk2-A(k))(p-D,k2-i(k))+1 = 0 (9) where k is the magnitude of k and I@) = p(k)-ico(k) is the scaled complex growth rate for wave vector k. 236 TEMPORAL AND SPATIAL INSTABILITIES FIG. l.-Domains in the ap-plane for different types of dynamics near singular points in two. component systems. 0s and AS correspond to oscillatory and aperiodic stability respectively- The character of instabilities are abbreviated to 01 A1 and SI corresponding to oscillatory aperiodic and saddle type instability respectively.Comparing eqn (6) and (9) we note that the dynamics of periodic patterns are obtained by subtracting DAk2and DBk2from a and /I,respectively corresponding to a displacement in fig. 1 towards negative parameter values along a line with slope D,/DA. These considerations lead to a classification of spatio-temporal behaviour near singu- lar points shown as domains 1-111 in fig. 2. On the left hand side of these domains there is global stability and on the right hand side the long wave limit is aperiodically unstable. In the latter case the non-linear terms are generally unable to prevent a FIG.Z.-Instability domains used in the discussion of spatio-temporal development beyond the linear approximation.The slope of the straight line is DB/DA. Consequently no lines with this slope through domain I intersects the domain for saddle type instability to the left of the hyperbola. In the figure DslD~m-9. I. BALSLEV AND H. DEGN global drift away from the singular point. In domains I and I11 the long wave limit* is oscillatory unstable while domains I1 and 111 have aperiodically growing solutions in a finite wave vector interval. Domain I1 and 111 vanish unless the self- activating component (here the A-component) is diffusing slower than the other component. The necessity of having D,/DA > 1 for obtaining aperiodical growth for finite wave vectors has been clearly stated by Edelstein and by Segel and Jackson.' The discussion of the further development of unstable systems beyond the linear approximation is limited to systems which do not drift away to other singular points and to closed systems having zero gradient of all concentrations perpendicular to the to the surface of the system.Such boundary conditions allow all space-independent solutions to the rate equations. Under these conditions the nature of the final state is closely related to the position in the parameter space of the point (a,/3) (see fig. 2). In Domain I-systems the homogeneous solution (k = 0) has the largest growth rate and is allowed by the boundary conditions. The final state of such systems is a homogeneous limit cycle. In Domain I1 the long wave limit is stable and there is saddle type instability in a limited wave vector range k < k < k,.If the boundary conditions allow the pres- ence of Fourier components in this range the system develops into states with large deviatims of the steady state. So far the ultimate destiny of such systems has not been analysed in general. Numerical integrations 6* on specific reaction schemes seem to indicate that the above situation generally leads to stable structures. How-ever at present it is not established if a system with specific spatial and chemical parameters allows a multiplicity of different stable patterns and to what extent these are accessible from a a slightly perturbed homogeneous steady state. From the linear theory it is obvious that below a critical size L,xn/k of the largest linear dimension there is spatial stability k2 being the largest wave vector with positive growth rate.In the case of niuch larger systems it would be interesting to search for general trends concerning the favoured spatial wave length of the final pattern. Oscillatory instability in the long wave limit and aperiodic growth of patterns is simultaneously present in Domain I11 (see fig. 2). Numerical integration in one dimension on such systems indicates that homogeneous temporal oscillations are strongly favoured if the concentration gradients are zero at the boundaries." In fact the existence of any form of stable patterns in closed Domain 111 systems is questionable. We can summarize that in closed two-component systems spontaneous oscillations tend to be homogeneous and spatial instabilities will develop into time-independent structures.For the experimentalist the most pertinent result of the above analysis is first that the formation of stable spatial structures may occur from an initially homogeneous solution. Secondly the interpretation of observed propagating structures in systems which are not oscillating when stirred cannot be based on two-component models. Thirdly in order to obtain spatial instability the self-activating component must diffuse considerably slower than the other one. In the search for systems capable of forming time-independent patterns one should look for reaction systems where an autocatalytic component is considerably more bulky than its precursor. This condition may not be fulfilled in the Belousov reaction where the autocatalytic species almost certainly is a fragment such as BrO of the bromate ion." Since bromate is present in excess it should not be considered as a variable component.The functional precursor of the autocatalytic species in the * The long wave limit corresponds to homogeneous solutions. TEMPORAL AND SPATIAL INSTABILITIES Belousov reaction is cerous ion. We have no data on the diffusion coefficients of these substances. However cerous ion is probably the bulkier one because of its tendency to form complexes. In the Bray reaction the autocatalytic substance is probably a fragment of the iodate ion and its functional precursor is molecular iodine. The large mass of molecular iodine and possibly its tendency to form complexes suggest a large diffusion coefficient compared to that of an 10 species.Thus neither of the systems mentioned above seem to have ratio of diffusion Coefficients required for spatial instability in a two-component model. 3. THREE- AND FOUR-COMPONENT SYSTEMS The absence from two-component systems of oscillatory spatial instabilities is believed to be a consequence of their simplicity rather than a general trend in reaction- diffusion systems. This suggests that models with more than two components should be investigated. We have studied a special class of three- and four-component systems nameiy those containing a chain of first order reactions with a feed-back on the first reaction from the last intermediate in the chain. Such systems can be symbolized as folIows s-I -+A-+B-+C-+ -1 I -+A+B-+C-+D-+ where the feed-back is activation or inhibition.The structure of the matrix of partial derivatives (cf. eqn (3)) for a chain of three components is given by where AfP = IF],+ Dpk2 t 12) and niA,mBand rn are positive (no self-activation). This system can be shown to be globally stable ifm c 0 (backward inhibition). For nz0 > 0 the system may have oscillatory or aperiodic instabilities. However as shown in the Appendix the oscillatory solutions have increasing growth rate for increasing wave lengths in the vicinity of marginal stability (dp(k2)/d(k2) < 0 for A(k) = +iu) whereas aperiodic spatial instability may occur. Thus this three-component system does not offer any spatial phenomena unknown in two-component systems.Preliminary studies lo of four-component systems by the methods shown in the Appendix indicate that backward inhibition may lead to oscillatory spatial instability. This requires a stable long wave limit (all four roots of the secular equation have negative real parts for k = 0) and the existence of wave vectors for which two rwts with p(k) = 0 and ~(k) + 0 have dp(k)/d(k2) > 0. From eqn (12) we have Thus oscillatory spatial instability requires that at least one of the derivatives dpld :VIP is positive for a marginally stable oscillatory solution. This is the case for large diffusion coefficients of the A-component and a large rate constant of the reaction which is subject to inhibition.I. BALSLEV AND €1. DEGX The above preliminary results for models with more than two components indicate that such models should be explored in more detail. The results should be related to chemical reaction schemes involving backward inhibition such as the model for the Belousov reaction. If models predicting oscillatory spatial instabilities turn out to be chemically realistic then a quite new field is opened not only in connection with instabilities. It is known in general linear wave mechanics that globally stable steady states in systems which are close to oscillatory spatial instability allow wave propagation with little spatial attenuation. APPENDIX Considering the matrix in eqn (11) of a three component system with backward inhibition we find that the complex growth rate A = p-iw is the solution to (MA+A)(MB+?J(Mc+A)-m,m,MA i-inOmAmB= 0.(14) With reference to eqn (13) we concentrate on the signs of ap/aMp. Investigating first the behaviour near marginally stable aperiodic solutions (E. 2~ 0) we require for a finite wave vector ko MAMBMC -(MA-mA)mOmB= 0. In the vicinity of ko Evaluating dp/dMp and inserting in eqn (12) and (15) we find and ap/akfA = (M&f,m~3/(DAk2N) ap/dh!f = -MAhfJN (17) where ap/aM = -M~M~/N N = MAMB+MAMc- MBhfC172A/DAk1e (18) The different signs of the derivatives ap/aMpallow the existence of aperiodic spatial instab- ility. For oscillatory solutions we find in the case of marginal stability In the evaluation of @/dMp we consider the left hand side of eqn (14) as a function MA MB,Mc A) and use the fact that Inserting in this expression the above value of A = +iowe find that all derivatives &/dMp are negative.Consequently oscillatory spatial instabiIity is absent in this three-component system. ' A. M. Turing Phil. Trans. B 1952 237 37. A. M. Zhabotinskii in B. Chance E. K. Pye A. K. Ghosh and R. Hess (eds.) Biologied and Riochemical Oscillators (Academic Press New York 1973) p. 89. P. Ortoleva and J. Ross J. Chem. Phys. 1973 58 5673. A. Lotka J. Plzys. Chem. 1910 14 271. TEMPORAL AND SPATIAL INSTABILITIES V. Volterra Lepw sur la Theorie Mathematique de In Lutte pour la Vie (Gautier-Villars Paris 193 1). I. Prigogine and R. Lefever J. Chem. Phys. 1968,48 1695. I. Balslev and H. Degn J. Theor. Biol.,1974. B. B. Edelstein J. Theor. Bid 1970 26 227. L. A. Segel and J. L. Jackson J. Theor. Biol. 1972 37 545. lo I. Balslev unpublished work. R. J. Field E. Koros and R. M. Noyes J. Amer. Chem. Soc. 1972 25 8649. l2 H. Degn J. Chem. Ed. 1972,49,302.

ISSN:0301-5696    

DOI:10.1039/FS9740900233

出版商:RSC    

年代:1974

数据来源: RSC

摘要:

Nucleation in Systems with Multiple Stationary States BYA. NITZAN,P. ORTOLEVA AND J. Ross Department of Chemistry Massachusetts Institute of Technology Cambridge Massachusetts 021 39 Received 26th July 1974 We consider a reaction diffusion system far from equilibrium which has multiple stationary states (phases) for given ranges of external constraints. If two stable phases are put in contact then in general one phase annihilates the other and in that process there occurs a single front propagation (soliton). We investigate the macroscopic dynamics of the front structure and velocity for two model systems analytically and numerically and for general reaction-diffusion systems by a suitable perturb- ation method. The vanishing of the soliton velocity establishes the analogue of the Maxwell con- struction used in equilibrium thermodynamics.The problem of nucleation of one phase imbedded in another is studied by a stochastic theory. We show that if the reaction dynamics is derived from a generalized potential function then the macroscopic steady states are extrema of the probability distribution. We use this result to obtain an expression for the critical radius of a nucleating phase and confirm the prediction of the stochastic theory by numerical solution of the deterministic macro- scopic kinetics for a model system. 1. INTRODUCTION Chemical reaction mechanisms with macroscopic rate laws of sufficient non- linearity in systems maintained far from equilibrium may have multiple stationary stable states for given external constraints.We refer to each such state as a phase and transitions between phases are possible. The analogy of the theory of phase transitions and critical phenomena to transitions between stable states and critical (marginal stability) points has been discussed in some detail both at the macro- scopic and statistical mechanical level of anal~sis.~-~ In this paper we investigate the nucleation of one phase within another phase as well as the conditions of coexistence of two phases (the analogue of the Maxwell construction). These problems have been considered by Kobatake who showed for a particular case the similarity of the behaviour of the generalized entropy production in transitions between stable branches of steady states and the Gibbs free energy in equilibrium phase transitions ; by Schlogl,* who took into account reaction and diffusion and treated coexistence of phases and the analogue of the Maxwell construction for one variable systems with the help of a mechanical analogy; and by Nicolis Malek-Mansour van Nypzlseer and Kitahara,lO who analyzed the onset of instability as a nucleation process derived a ncn-linear master equation for that purpose and applied that result to some ex-amples.We approach these problems from two diff'erent points of view. First in section 2 we investigate the macroscopic dynamics of two stable stationary states of semi-infinite extent at given external constraints placed in contact with each other. Except for one value of the constraints one phase annihilates the other and in that process there occurs a single-front propagating wave a soliton.We investigate the behaviour of the concentration profile and velocity of such a wave for two model systems analytically and numerically and show that the values of the external con- straints for zero soliton velocity establish the Maxwell construction. 241 242 NUCLEATION Further we study these quantities for general reaction-diffusion systems by a suitable perturbation method and find that the soliton velocity vanishes linearly in the deviation of the external constraint from its value at the Maxwell construction. Next in section 3 we discuss some aspects of a stochastic theory of instability phenomena on the basis of a Fokker-Planck equation assumed for a reaction-diffusion system.That equation can be solved for the steady state probability distri- bution for which the reaction dynamics is obtainable from a generalized potential function. In that case we show that the macroscopic steady states are extrema of the probability distribution maxima for stable steady states. We apply this theory to analyze the behaviour of a system consisting initially of a nucleus of one phase embed- ded in an infinite bulk of the other phase ; we do so in a manner similar to the conven- tional treatment of nucleation in first order phase transitions.". l2 We thus obtain an expression for the critical radius for such a nucleus above which the nucleus grows and below which it disappears in time. This prediction agrees with numerical solutions of the deterministic diffusion-reaction equation for a simple model system.2. MACROSCOPIC DYNAMICAL THEORY We consider a diffusion-reaction system for which the niacroscopic deterministic equation of motion is The symbols denote $ a column vector of concentrations and possibly other state variables such as temperature; Y spatial coordinates; t time; B a matrix of diffusion coefficients ; F,the variations due to chemical reaction ; and X the set of external constraints (boundary concentrations of species light intensity in an illurn-inated ~ystern,~ etc.). We assune that 1. the diffusion matrix is constant and sym- metric and the diffusion process is stable in that in the absence of reaction eqn (2.1) with F = 0) the steady state solution is stable; 2.P is analytic in 9 and 1. We consider systems for which there exist two stable stationary states (and one unstable stationary state) in a given range of k. Thus if in a one-dimensional system (co- ordinate x) we prepare two systems under identical external constraints one system in one stable state (one phase) and the other system in the other stable state (the other phase) then the boundary at x = 0 between the two systems will move in one direction or another depending on X. There occurs a single front propagation,13 a soliton in which one phase annihilates the other. We shall find the condition for coexistence of two such systems (phases) to be zero soliton velocity or a standing front at a given X. We therefore seek wave solutions of (2.1) which are most readily studied in a reference frame moving with the front.W7edefine the phase 4 = x-ct (2.2) where the soliton velocity o is a function of k. With (2.2) the partial differential eqn (2.1) is converted into an ordinary differential equation which is a nonlinear eigenvalue equation for u at a given X. Coexistence occurs when u = 0 at h = AM. The value of hMis obtained from eqn (2.4) which may be inrer-preted as an eigenvalue equation for AM (2.4) A. NITZAN P. ORTOLEVA AND J. ROSS We first consider a soluble model system and then present an analysis for more general {D F}systems. The model system consists of a single species with third-order kinetics which obey the equation a2t,b all/-= D 7-q(t,b-a)(t,b-b)(ll/-c).at dx-Homogeneous steady states occur at $* = a,b c where these parameters are arbitrary functions of 2.. We take b 2 c > a q > 0,and thus only states a and b are stable to small homogeneous perturbations. If we define a new variable u such that *-a UE- 6-a' then we obtain dU a2u -= D-2-q(b-a)2u(u-l) at ax This is identical in form to a model considered by Montro11.12 In terms of our notation his solutions for soliton fronts are = -(b-a). (2y The front velocity is given by v = (+qD)*(a+b-2~). (2.9) Coexistence occurs at u = 0 i.e. c = +(a+b) at 2 = AM (2.10) which is obvious from the symmetry of the trilinear kinetics for which the two stable states a and b are located symmetrically about the unstable state c.It is interesting to take a particular case and analyze it further. Let F = -q@($2-2$ + 1 /A). (2.1 1) Then v = (+qD)3[3(1-1 /A)+-I] and = (q/2D)+[1+(1-1 /A)+] from(2.8,9). The uniform steady states and the soliton velocity are shown in fig. I(a) and (b). Co-existence is found at AM = 918. With this we obtain the structure of the coexistence region to be $coex .( x = $( +exp(kPMMX))-(2.12) where & = 3(8q/D)+;the structure is shown in fig. 2. We note from fig. l(b) that the velocity passes through the origin at AM with a finite but non-zero slope l(du/dL)AMl# co. This general property provides a useful aid in determining AM experimentally. Next we consider the general system (2.1) by analyzing the solutions of the front eqn (2.3)in the vicinity of AM.To do so we introduce a parameter A to be defined at a later stage in the development which measures the deviation of the system from coexistence conditions. In general A is a function of A 2 = &A" (2.13) n=O 244 NUCLEATION 0 5 10 A FIG.1 .-Steady states and soliton velocity versus the parameter A for the model system eqn (2.1 1). and we expect that A approaches zero proportional to the factor (A-AM) some positive power. Furthermore 1,6 and v depend on A raised to 03 (2.14) rm v = C vnAn. (2.15) n=O Such expansions have been employed elsewhere for the study of chemical waves.lS On substitution of the expansions (2.13-15) into (2.3) we recover to zero'th order in A FIG.2.4tructure of the coexistence region versus distance eqn (2.12) for the model system eqn (2.11).A. NITZAN P. ORTOLEVA AND J. ROSS 24s eqn (2.4) with the boundary conditions of the system being one stable stationary state at x = -a,and the other stationary stable state at x = +a. We also have by definition in this order uo = u(R& = 0. To first order in the parameter A we find (2.16) with the operator 2’ defined by d2 9 = D-+R (2.17a) d242 (aF/a$) (2.17b) and the subscript 0 following differentiation implies evaluation in lowest order. We note that differentiation of (2.4) with respect to 4 yields (2.18) We assume that d$/d@ is the only eigenfunction of 9 with zero eigenvalue. P-l(d$/d#) does not exist and hence the solution of (2.8) for t+h1cannot be found for arbitrary (vl R1).The problem is overcome by a proper choice of these coefficients. We then define a vector 0 such that 9+0= 0 (2.19) where P+is the adjoint operator of 9.If we define an inner product (LIP)= c d4 44 (2.20) 1 -53 where the sum on i runs over all species then from (2.17,20) L?+ is given by d2 9+= D+-+a+. (2.21) d+2 (By assumption we have D+ = Dt Onsager’s relations). We can be assured of the solubility of the first order eqn (2.16) if we can choose v1 and A1 such that the r.h.s. of (2.16) is orthogonal to the zero eigenvalue function d$,/d4. Thus we obtain by use of the adjoint eigenfunction to d$o/d4 that is 0,see (2.19) that u,(@I(d$o/d4)) +M@l(aFlwo)= 0. (2.22) The coefficient of v may be chosen to be unity by proper normalization of 0.Hence if the coefficient of A1 is finite then this relation can be satisfied for finite vl R1. Finally (aF/d;l) approaches constant values as 1&-+ co. A necessary condition for al solubility is 1 d$lO(4)l < 00 in order to find R1 vl # 0. For systems where -0 R+ = i-2 9is self adjoint and therefore 0= d$o/d&. Since d$o/d$ is localized to the region of the soliton front it satisfies the necessary condition of solubility. However we have not investigated this condition for systems with R+# s1. It is sufficient to all orders to choose Rn32 = 0 ; thus with the convenient choice R1 = 1 we have A R-RM. (2.23) To first order the perturbation solution is b Y-Y 246 NUCLEATION + ' VI4M -~(ol(aF/a~)M)/(ol(d~O/d#))](~-~~) ' (2.25) For the one variable system (2.5-12) we have verified these relations using the equation 0 = dt,bo/d4.I I --4.0 5.0 6.0 So FIG.3.-Steady state concentration and soliton velocity versus Sofor the model system described in Appendix A. For curve (a) the diffusion coefficients Ds Dp are both 1 and for curve (b)DS = 1 and Dp = 0.05. 0.0 8.0 16.0 X FIG.4.-Soliton structure against distance at different times for the model system Appendix A. The diffusion coefficients are Ds= 1 Dp = 0.05. The external concentration S" = 5. A. NITZAN P. ORTOLEVA AND J. ROSS Finally we note that the structure and velocity of the soliton depends in general not only on the reaction term P but also on the diffusion matrix D.If D is equal and diagonal then it may always be eliminated from the equation for the soliton front by a suitable scaling of the length 4-+(JD)4'. However for more complex diffusion matrices this simple scaling no longer holds. In Appendix A we consider a two species system with rate mechanism & S"$S+P+P" (2.26) where So,Poare maintained constant and the activity of the enzyme 8depends on P. In fig. 3 we show the steady state concentration P* as a function of s"and the results of the numerical integration of the partial differential equations yielding the variation of the soliton velocity against the concentration of s" for various choices of diffusion coefficients of S and P. For this system the value of Soat coexistence (AM of the general theory) does not appear to depend on the choice of diffusion coefficients to within the accuracy of the numerical computations; the results of the stochastic theory section 3 indicate no dependence on that choice for systems for which the kinetics is derivable from a generalized potential function.In fig. 4 we show the soliton stucture for the same system (Appendix A) for the choice of diffusion coefficients Ds = 1 Dp = 0.05 and So = 5. 3. STOCHASTIC THEORY Our treatment is based on a model in which to the set of reaction diffusion rate equations we add a stochastic source term to account for fluctuations in the system. The set of concentrations at each spatial point becomes then our set of stochastic variables. It is convenient to make this set discrete by dividing the system into cells of homogeneous composition.We thus consider the system of volume V to be partitioned into N cells each of volume ij = V/N. The macroscopic reaction kinetics is again given by with a similar equation for $n in the nth cell. The rate of change of numbers of particles of any species x, = fiq,,,due to diffusion is where the sum extends over the 2d nearest neighbour cells (d is the dimensionality of the system). The factor of Ui-2/d comes from multiplying by the surface area ijl-l'd of the cell and dividing by the characteristic length of a cell such that ($n+i-$,,)/fi*/d is an approximation to the concentration gradient. Now we write a Langevin equation for this system as (3.3a) or (3.3b) For simplicity we take the stochastic term to be Gaussian and thus to satisfy the usual relations (A,(?))= 0and (&,,(~')h,Jt)) = 6,p6,,,$(t-t') (a,/3 denote different chemical NUCLEATION components).Eqn (3.3b)is equivalent to the following Fokker-Planck equation for the probability density of finding the concentration distribution $,, P[Y], in the system where Mn(Y](with components MaJ3!]) is given by 2d and where Y denotes the “ supervector ” *n (n = 1 . . . N) (= (It = 1 . . . N a goes over the different chemical species)). In (3.4) K was taken for simplicity as a scalar matrix K = KI where I is the unit matrix. The integration of this many- variable Fokker-Planck equation for steady state conditions (Plat) = 0 is possible provided the reaction kinetics are such that the curl of’M is zero,I6 aMan -aMpn (3.6) 89s wa This condition should be modified if D is not a symmetric matrix or if K is not a scalar matrix.If (3.6) holds then we obtain P[W] = Cexp(-U[Y]) (3.7) where C is a normalization constant. For simplicity we assume now that we have a one-component system for which (3.6)always holds. The results for a more general case with a symmetric diffusion matrix are given in Appendix C. We write U[+]as which in the continuum limit becomes The function G[$(v)]is defined by (3.10) which in case of many chemical species should be understood as a line integral in $ space. Condition (3.6) ensures that this line integral is well defined and does not depend on the integration path.It is interesting to note that (3.9) has the same form as the free energy in the Landau-Ginzburg theory (there G($) is usually taken as a quartic form). Also critical phenomena have been discussed with Hamiltonians of that form.17 Having arrived at an expression for the steady state probability distribution P[$] we show next that extrema in that distribution are given by the steady states of the deterministic equation. Differentiation of (3.7) yields (3.11) A. NITZAN P. ORTOLEVA AND J. ROSS 249 or hp= dr{F[$]G$-DV$ VS$) (3.12) which on partial integration becomes K dr(F[1,b]+DV~$}6$. (3.13) For P to be an extremum for the arbitrary variation S$ we must have the integrand vanish ; that is the macroscopic condition for steady states.It can be further shown that stable steady states correspond to maxima of the probability distributions (see Appendix B). The system may be in a single stable stationary state (although there exist other such states at the same external constraints) and in that case the system is homogeneous i.e. F[$] = 0. However two (or more) stable stationary states may coexist and the distribution of concentrations is the solution of DV'++F[$] = 0 for given boundary conditions. Note that this is exactly the equation for the zero velocity soliton structure (eqn (2.3) with u = 0). The stochastic theory for the probability distribution of a system in a steady state may be applied to the problem of nucleation. Consider a spherical nucleus of radius R of one stable stationary state labelled phase B immersed in an infinite bulk of another stable stationary state labelled phase A both under the same external constraints.Thus we have F[@J= F[i,bB]= 0 for the individual homogeneous phase. For phase B immersed in A however we may inquire about the expected stability of the spherical nucleus that is the nucleus is expected to evolve in the direc- tion of increasing the probability of the overall structure (nucleus +bulk). The Iogarithm of the ratio of the probability distribution for B phase immersed in A phase, PBA,to that of pure A phase PAis where I is the width of the interface region and where we have approximated the gradient by ($B-$A)/Z. In eqn (3.14) we have used the following expression for the surface area S and the volume V of a d-dimensional sphere (ris the r-function) The derivative (with respect to R) of the r.h.s.of (3.14),set equal to zero gives the radius of the sphere of phase B for which the 1.h.s. of (3.14) is a minimum that is (for a given I) (3.15) If G[#,] > G[@*] then there exists a distribution (characterized by R,) of maximum "free energy " (minimum probability) such that for R > R the nucleus will grow (as this is the direction of increasing the probability) and for R < R it will disappear in time. If G[I,~~] < G[$A] then no positive R exists. Finally if G[$J = G[t,bA]then R = a,and phases A and B coexist with a planar (zero curvature) surface of contact. This condition is clearly the analogue of the Maxwell construction in equilibrium 250 NUCLEATION phase transitions.Furthermore we expect that this condition is equivalent to the one of zero soliton velocity and have confirmed the expectation for the example given by (2.5). In a many component curl-free system with a symmetric non-scalar diffusion matrix the result (see eqn (C.l) Appendix C) shows that the Maxwell construction is again given by the condition G[31g]= G[$J. We stress that the arguments presented here are based on probabilistic consider- ations especially on the postulate that a system will evolve in the direction of increasing the probability. The result (3.15) is similar in many respects to that for the radius of critical nucleus in equilibrium first order phase transitions.' * l2 The contribution of the diffusion process to the " free energy " U eqn (3.9) takes here the place of the surface tension in the equilibrium nucleation treatment.The appearance of the parameter I io eqn (3.15) takes into account the non-zero thickness of the interphase layer. It should be remembered that I is not an independent parameter but is deter- mined by the external parameters. Also the simplifying assumption of linear varia- tion of the concentration t,b across the interphase layer limits the applicability of the result (3.15). We compare the results of the stochastic theory to the solution of the deterministic dynamical equations. To this end we numerically integrated the one variable partial differential equation !?!! = DV2t+b-J13+bt,b+~ (3.16) at in which the cubic polynomial F($) = -G3+b$+c plays the role of a chemical source term.For b > 0 and Jcl < (2b/3)(,/b/3) this equation has two stable and one I .€ 4 O.( -1.t 1 t=O I J 0.0 4.5 91 1.0 4.5 9.0 r r FIG.5.-Evolution of a nucleus as a function of radial coordinate for various times as obtained from numerical solution of eqn (3.16) with b = 3 c = 1.S and D = 1. Figures at left and right are for initial radius smaller and larger than & respectively. A. NITZAN P. ORTOLEVA AND J. ROSS 25 1 unstable homogeneous steady states. Denoting the two stable homogeneous steady states by (phase A) and t,kB (phase B) we put our initial condition to be a sphere of radius R centred around the origin and consisting of phase B imbedded in a spherical region of phase A which in the computation process has to be taken finite but with radius R1much larger than R.Since the initial condition is spherically symmetric and the coefficients of (3.16) are constants it is sufficient to consider the radial coordinate only. The boundary conditions are taken as aIrl//ar = 0 at r = 0 and at the larger bounding sphere. Fig. 5 shows the time and space dependent solution to eqn (3.16) in three dimen- sions for the choice of coefficients D = 1 b = 3 c = 1.5 and for two different choices of the initial radius R R < R and R > R,. The nucleus shrinks and disappears for R < R,and grows in time for R > R,. It is interesting to note that the growth process (for R > R,)accelerates for small R but the velocity approaches a constant for larger radii as the front becomes more planar.We calculate the velocity of a planar front (from eqn (2.9) using the steady state solutions t,ka = -1.38 t,hb = 1.94 (stable) and Irl/c = -0.558 (unstable)) to be 1.2 while from fig. 5b the velocity for the largest radius obtained is 1.0. TABLECRITICAL RADIUS FOR NUCLEATION FOR THE SYSTEM GIVEN BY EQN (3.16) Rccalculated from eqn (3.15) Rc from numerical computations with 1 = 2 -C 2 dimensions 3 dimensions 2 dimensions 3 dimensions 0.50 2.6 < R < 3.2 3.4 0.75 3.5 < R < 5.5 4.6 1.oo 1.5 < R < 1.75 2.7 < R < 3.7 1.7 3.4 1.25 2.1 c R < 2.3 2.6 1SO 0.85 < R < 1.0 1.65 < R < 1.95 1.05 2.1 1.75 0.8 < R < 1.0 1.35 -c R < 1.5 0.9 1.75 1.875 1.2 < R < 1.4 1.6 2.0 R < 0.5 R < 0.8 0.7 1.35 In table 1 we present the results of numerical computation of R,.as well as results based on eqn (3.19 for b = 3 and for different values of c in two and three dimensions.The calculations based on eqn (3.15) were done with I = 2 which is an estimate based on the numerical results such as given by fig. 5b. The agreement between the values of R,as obtained from (3.15) and from numerical integration of (3.16) is fairly good for all the values of the parameter c that were tried excluding the bifurcation value c = 2. (Our failure to provide a more accurate numerical answer for R near the bifurcation point is caused by the slow time evolution of the system near that point which makes the numerical integration very expensive.) This agreement is remark- able as the two ways of obtaining R,are different; one is based on solutions of the time-dependent deterministic equation and the other on the steady state probability distribution.To conclude this section we should note that even though the dynamics of nucleus growth as given by eqn (3.16) or fig. 5 give important information on the rate of nucleation of one phase within another the rate of this process depends also on the rate of formation of nuclei by spontaneous fluctuations which is not discussed in the present work. We thank John M. Deutch for helpful discussions. This work was supported in part by the National Science Foundation and Project SQUID Office of Naval Research.252 NUCLEATION APPENDIX A For the reaction mechanism (2.26) the reaction-diffusion equations are as -= 82s -2D +(SO-S)-G(P)S at ax ap a2p --D -+(Po-P)+d(P)S at ax2 with all rate coefficients taken to be unity. The assumed form for the dependence of 8 on P is P2 (A.3) &(P) = 1+P+P2 for which chemical examples are available.'* The steady state solutions of (A.l,2) are P* = 0; t(So-1)+$[(So-1)2-8J* 64.4) which are shown in fig. 3a. APPENDIX B Here we show for the probability distribution given by eqn (3.7) and (3.9) that stationary points which correspond to stable macroscopic steady states are local maxima of the distribution or minima of the potential U[$] given by (see eqn (C.l)) U[$l = -1lK 1dr{G[$(r)l-4D :(V$m2}.(B.1) Let the steady state variable $*(r) satisfy the equation F[$*(r)]+DV2$*(r) = 0 (B.2) and let $(r) = $*(r)+a$(r). We assume that \II*(r)corresponds to a stable steady state of the system. The first order (in S$) term of U then vanishes and the second order term is given by a2u = -3 JdrslC/(n+ Dv2)S$ where i2 = (iW/a$)+*. We have used the identity JdrVWeV V= -J drWV2V+1 WneVVdo (B.4) and the fact that $*(r) and \II(r)must satisfy the same boundary conditions and there- fore = 0. In eqn (B.4) E is the boundary of the system do is an element of this boundary and n is a unit vector perpendicular to the boundary. Since we assume that $*(r) is a stable steady state the operator +DV2 must be negative definite and therefore d2 U in eqn (B.3) is positive.This establishes the fact that U[$] is a mini-mum (and P[$]is a maximum) for $ = $*. APPENDIX C Here we write the modified relations of section 3 for a many variable system and for a general symmetric diffusion matrix. The potential U[$] (eqn (3.9)) is given by wwl = -1 lq dww91-3D (V$(rN2} (C.1) A. NITZAN P. ORTOLEVA AND J. ROSS where The expression for the critical radius R is now replaced by N. Rashevsky Mathematical Biophysics (Dover New York 1960) Vol. I. B. B. Edelstein,J. Theoret. Bid 1970 29 57. P. GIansdorfT and I. Prigogine Thermodvnamic Theory of Structure Stability and Fluctuations (Wiley-Interscience New York,1971). G. Nicolis Adv. Chem. Phys.1971 19 209. A. Nitzan and J. Ross J. Chem. Phys. 1973 59 241. H. Hahn P. Ortoleva and J. Ross J. Theoret. Biol. 1973 41 503. ’Y. Kobatake Physica 1970,48 301. F. Schlogl 2.Phys. 1972 253 147. A. Nitzan P. Ortoleva J. Deutch and J. Ross J. Chem. Phys. 1974,61 1056 see also other references therein. lo G. Nicolis M. Malek-Mansour A. van Nypolseer and K. Kitahara private communication (preprint). L. D. Landau and E. M. Lifshitz Statistical Physics (Addison-Wesley Reading Massachusetts 1958.) l2 see for example B. Widom in Phase Transitions and Critical Phenomena Vol. 2 ed. C.Domb and M. S.Green (Academic Press New York 1972). l3 R.Davidson Methods in Nonlinear Plasma Theory (Academic Press New York 1972). l4 E. W. Montroll in Statistical Mechanics ed.S. A. Rice K. F. Freed and J. C. Light (Univ. of Chicago Press 1972) p. 69. P. Ortoleva and J. Ross J. Chem. Phys. 1974 60,5090. l6 R. L. Stratonovich Topics in the Theory of Random Noise (Gordon and Breach New York 1963). see for example S. Ma Rev. Mod. Phys. 1973 45 589. l8 H. Hahn A. Nitzan P. Ortoleva and J. Ross Pruc. Nat. Acad. Sci. (submitted).

ISSN:0301-5696    

DOI:10.1039/FS9740900241

出版商:RSC    

年代:1974

数据来源: RSC

摘要:

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.

ISSN:0301-5696    

DOI:10.1039/FS9740900254

出版商:RSC    

年代:1974

数据来源: RSC

摘要:

GENERAL DISCUSSION Dr. P. Ortoleva and Prof. J. Ross (MZT)said We have investigated the problem of the establishment of gradients of concentrations in a chemical system far from equilibrium which is originally homogeneous.' For this purpose we chose a volume surrounded by a membrane and that composite system is immersed in a bath of fixed concentrations. Note that the boundary condition across the membrane for each M -Memb raw-> M distance FIG.1. concentration is a flux proportional to the concentration difference across the mem- brane. For cells this is likely to be a physically better condition than one of constant concentration or constant flux. We formulated a general theory for symmetric and asymmetric perturbations and then derived for a model system the equations for the symmetric and asymmetric stable states.Thus we showed that from an initially stable symmetric state in a system immersed in a bath of constant concentrations (dotted line in the figure) an asymmetric stable state (solid line) evolves in certain ranges of kinetic and transport coefficients and bath conditions. A different problem of spatial order is that of Liesegang rings referred to earlier in this Discussion. A theory and some experiments on this topic have been presented by Flicker and Ross2 Dr. A. Babloyantz and Prof. G Nicolis (Brussels) (partly communicated) :To our knowledge the paper by Ortoleva and Ross refers to a discontinuous system of two compartments subject to a flow of matter from the outside. One cannot really refer in terms of " boundary condition " to the flux across the external " membrane " as this flux is incorporated into the kinetic equations themselves.Asymmetric states of concentration under similar conditions were first obtained by Prigogine and Lefe~er.~ In spite of its interest and its impact as the first model for dissipative structures a discontinuous two-box model suffers from the inconvenience of imposing artificially a wavelength on the system. This point is discussed in detail in Prigogine's paper at the second Versailles conference From Theoretical Physics to Biology (Versailles 1969). When diffusion is taken into consideration either in a P. Ortoleva and J. Ross Biophys. Chem. 1973 1 87. M. Flicker and J. Ross J. Chem. Phys. 1974,60,3458.I. Prigogine and R. Lefever J. Chem. Phys. 1968,48 1695. 263 GENERAL DISCUSSION continuous formalism or by dividing the space into at least 3 cells then the system can evolve to a structure whose wavelength is no longer imposed but is instead determined by the system's parameters the boundary conditions and the size. Under these conditions it has recently been shown that polarity can appear spontaneously in an initially homogeneous morphogenetic field. Dr. P. Ortoleva and Prof. J. Ross (MIT) (communicated):Prigogine and Lefever' showed with a model reaction in a two-box enclosure that gradients of chemical species may be established under the conditions that these species (X Y)do not inter- change with the surroundings. We showed that in a similar system (eqn A.7 8) a gradient of species can be established under more restrictive conditions such that the species for which the gradient is formed may exchange with homogeneous surroun- dings.2 Thus despite the driving force by the constant concentrations in the sur- roundings (bath) to erase any gradient such a gradient arises. This is of importance in those cases of morphogenesis taking place in homogeneous surroundings for which the morphogens may move across the boundary of the system. The two-box model is indeed a simplification but nonetheless is related to the continuum partial differential reaction diffusion equations with boundary conditions as may be shown by a mode-mode coupling theory. The lowest order approximation in that approach yields a cosine distribution in position with an amplitude obeying a nonlinear equation which is essentially that for the two-box model.Dr. J. S. Turner (Azutilz Texas) said In Section 2 of the paper by Goodwin a model for membrane excitability is proposed with the qualification that full experi- mental justification for a phase transition interpretation of the phenomenon is not yet available. Insofar as the underlying chemical kinetics of an appropriate model for membrane excitability is found to exhibit multistationary states then interpretation as a phase transition may well be accurate. The question of phase transitions in systems for which macroscopic chemical kinetics predicts several stable states has been studied recently by Lefever Prigogine and myself.Analyzing spontaneous localized fluctuations in a simple chemical model we have verified the existence of metastable states in the coexistence region and have obtained the minimum coherence length L for spontaneous fluctuations necessary to induce a transition from such a state. (For details see my discussion remark following the paper of Nicolis and Prigogine in this Symposium.) In this view for example the phenomenon of latent excitability referred to at the end of Section 2 would correspond to naturally occurring fluctuations having a characteristic spatial extent L < L, a consequence possibly of physical constraints on membrane-related processes. Dr. A. Babloyantz (Brussels) said In his paper Goodwin considers the sudden appearance of spontaneous action potential during the regeneration of the unicellular marine alga Acetabularia.But the differential equations constructed to describe his model are linear ones and cannot give a spontaneous generation of inhomogeneity in a given physico-chemical parameter starting from a homogeneous situation. The inhomogeneous result found by the author is due to his special choice of the so-called pick-up function. In particular for the special case of Acetabularia any theory has to explain how in anucleated fragments one can have a generation of polarity giving rise either to one cap formation at one end of the fragment or two caps at either end I. Prigogine and R. Lefever J. Chem. Phys. 1968 48 1695. P. Ortoleva and J. Ross Biophys. Chem.1973. 1 87. GENERAL DISCUSSION of the fragment. Can the author account for these experimental facts? And if so how is the number of caps determined within his linear theory? Dr. B. C. Goodwin (Sussex) said Babloyantz and Auchmuty questioned whether the model of wave propagation and gradient formation on a membrane was spatially uniform initially and whether it is dynamically stable. It is evident from fig. 1 that all the molecular components of the model are distributed uniformly in the membrane so that there is no initial asymmetry. Once a wave begins to propagate from some point however there is first dynamical asymmetry (i.e. at any point in space one process is followed by another in a defined manner as determined by the pattern of interactions described by the kinetics of the reaction scheme) and subsequently spatial asymmetry in the bound metabolite arising from the dynamical asymmetry together with wave propagation.This leaves open the question how the wave origin is determined which was done deliberately in order to cover the various cases to which the model could be applied. These include (1) an external environmental event determining where the symmetry-breaking wave is initiated as in the case of light acting upon a fertilized Fucus egg; (2) an initial gradient of some kind which establishes the wave origin as in an amphibian embryo or a decapitated Hydra; (3) a system with neither initial polar order nor an environmental asymmetry which can break symmetry spon- taneously as a result of spontaneous fluctuations which exceed the level for initiation of the wave.In this instance more than one wave origin can arise so that com- peting organising centres may result as occurs in bipolar Acetabularia segments. The model is applicable to all these instances but it does not of course explain them. Regarding stability the basic equations of the model as described for example by Goodwin are necessarily stable when P(x,t)is a propagating delta function (the case analysed) because a finite amount of substance is produced over any finite spatial region of the membrane. Thus there cannot in this instance be any question of instability and the problem considered is simply the spatial distribution of the morphogen resulting from the passage of a wave along a membrane.Dr. J. F. G. Auchmuty and Prof. G. Nicolis (Brussels) (partly communicated) In their paper Balslev and Degn appear to overlook the effect of boundary conditions on the solutions of reaction-diffusion equations. Their conclusions that (D,/D,) > 1 is a necessary condition for obtaining an aperiodic instability and that in domain I of fig. 2 the system will necessarily evolve to a uniform limit cycle only hold under special boundary conditions. An analysis of the one-dimensional trimolecular model in the case where the concentrations are held fixed and equal to their uniform steady state values at the boundaries (see the paper of Nicolis and Prigogine at this Symposium) shows that the first instability may give rise to a space-dependent and time-periodic structure.More-over such structures arise before aperiodic solutions provided where we are using the notation of Nicolis and Prigogine. Here D,,D are the dif- fusion coefficients of X (the ‘‘ self-activating ” component) and of Y. In particular D,/Dx need not be greater than 1 to lead to an aperiodic stability when A is suffi- ciently small. B. C. Goodwin Theoretical Physiology in Selected Topics in Physics Astrophysics dnd Bio-physics. ed Laredo and Juristic (D. Reidel Publ. Co. 1973) pp. 381-420. GENERAL DISCUSSION The reason for this is the fact that linear stability analysis for one dimensional systems of length Z with fixed concentrations at the boundaries allows modes of the form for n = 1,2 . . .. The state with n = 0 is now excluded because of the boundary conditions and so one gets these different results.Fixed concentration boundary conditions may arise if the reaction takes place inside semipermeable membranes and there is a reservoir of reactants outside. Prof. D. Thomas at Compihgne has been conducting experiments with these boundary conditions. Dr. A. Winfree (Indiana) said Balslev and Degn point out that in the linear approximation near a homogeneous stationary state propagating waves cannot occur in a 2-component reaction/diffusion system unless it oscillates while well stirred (i.e. with V2 term omitted eigenvalues are now real) and that initiation of such spatiaZ instability requires that the diffusion coefficients be unequal. As usual with linear systems one must enquire how much non-linearity or how great an excursion from the stationary state is required to admit the forbidden phenomena.I submit by way of answer only an isolated example but it is an in- structive one I think aA -= V2A-A-B+H(A-0.05) at aB -= V2B+A/5 at where H(x) = 1 if x > 0 or 0 otherwise. It will be noted that the diffusion coeffi- cients are equal and the equation is piece-wise linear the two pieces meeting along A = 0.05. In the well-stirred case (V2 terms omitted) the stationary state (0,O)is a unique global attractor with negative real eigenvalues -0.72 -0.27. Yet this equation’s stable solutions (depending on initial conditions and boundary conditions) include not only the homogeneous stationary state but also propagating plane waves and concentric circular waves and spiral waves polarized to rotate clockwise or counterclockwise (unpublished numerical computations).Prof. I. Balslev and Dr. H. Degn (Denmark) said In answer to Auchmuty and Nicolis we remark that in our paper we have emphasized that the boundary conditions we wish to consider are the ones a chemist will always assume if nothing else is ex-plicitly stated namely the ones of “an ordinary vessel with non-reacting walls ”. We have abstained from discussing the consequences of other boundary conditions such as fixed concentrations at the walls because we would feel obliged to discuss also the operational meaning of such boundary conditions. We have not been able to think of any operations which will assure fixed concentrations at the walls.In answer to Winfree we affirm that the hypothesis of the existence in certain reaction systems of spatial instability towards infinitesimal perturbations is a fun- damental one. The aim of our paper is to contribute to the clarification of the con- ditions required for such spatial instability. The interesting observations of propagating patterns presented by Winfree in this symposium seem to be dependent on Jinite perturbations. Likewise Winfree’s example given in his comment has a globally stable steady state. Therefore a departure from the corresponding homo-genous state can be brought about only by ajinite perturbation. GENERAL DISCUSSION Dr. R. G. Gilbert and Mr. S. McPhail (University of Sydney) said In theoretical discussions of reactions which produce spatial oscillations an obvious point which is frequently overlooked is that a stable final stationary inhomogeneous state I/~(x) is required in addition to an unstable initial state.In a one-dimensional system described by a$pt = F(+) +Da*+/aXz (1) (where t) is the usual vector of concentrations etc.) then t,V is defined by putting the right-hand side of (1) equal to zero. An indication of the stability of +I is that the real parts of the eigenvalues of the matrix aF/a$I be negative for all x. We have written a general program for numerically determining $I requiring the solution of coupled differential equations with double-ended boundary conditions (which as the equations are very badly conditioned poses considerable numerical difficulty).esting this program over a small range of x on the Brusselator mechanism we found not only the solution derived by Lefever by complete numerical solution of (I) (which was of course stable) but also several unstable solutions one of which was asymmetric and hence relevant to Babloyantz’ previous remark. We feel that the search for stable stationary inhomogeneous states will prove most useful in under-standing morphogenesis. Dr. A. Babloyantz and Dr. J. Hiernaux (Brussels)(partry communicated) :We wish to present a mechanism explaining the spontaneous onset of polarity in an initially homogeneous morphogenetic field. Ever since Turing’s * important work it has been known that a reaction diffusion system of two or several morphogens may account for many phenomena of morphogenesis.In particular the spontaneous formation of inhomogeneous patterns was viewed as the result of the instability of the homo- geneous state triggered by random disturbances. However due to the particular boundary conditions used by Turing the patterns he obtained did not present any polarity. This was also the case for all subsequent work on dissipative struct~res.~-~ Consider a reaction-diffusion system of two morphogens nz and m2 described by the following kinetic equations fand g describe chemical reactions which must obey kinetics similar to those proposed by Turing. We impose zero flux at the boundaries of the system. Then it can be shown 6* that for a proper choice of parameters there is a critical length Lo below which the system remains homogeneous.R. Lefever J. Chem. Phys. 1968,49,4977. A. M.Turing Phil. Trans. R. SOC.B 1958,237 37. P. Glansdorff and I.-Prigogine Thermodynamics of Structure Stability and Fluctuations (Wiley-Interscience New York 1971). I. Prigogine and R. Lefever J. Chem. Phys. 1968,48 1965. J. F. G. Auchmuty and G. Nicolis Buff. Math. Biof. 1975. A. Babloyantz and J. Hiernaux Science 1975 to appear. A. Babloyantz and J. Hiernaux Bull. Math. Biol. 1975. GENERAL DISCUSSION When the length of the system is greater than Lo a one-wavekngth structure appears in the system and is of the form seen in fig. 1 which presents a polarity. For higher values of L there appear structures of several wavelengths.72 m cell position FIG.1. Dr. B. L. Clarke (AZberta) said Balslev and Degn have demonstrated the desir- ability of being able to study reaction-diffusion instabilities which cannot occur in systems with too few reactants. I would like to indicate how the steady state stability of large reaction-diffusion systems can be determined. The stability analysis of large realistic chemical systems is greatly simplified by the fact that the rate constants usually range over many orders of magnitude. The matrix elements of eqn (2) and (11) (of the paper by Balslev and Degn) are then also of different orders of magnitude and the equations in the Appendix can be simplified to a few dominant terms. For large reaction systems this simplification makes an otherwise unwieldy algebraic problem tractable.Since the stability domains are determined by the signs of various polynomials the curves separating the stability domains have approximate equations which state that two possible dominant terms of opposite sign are equal. Such equations are always hyperplanes when plotted logarithmically in the parameters. From this it follows that the stability domains of a reaction-diffusion system are asymptotically convex polyhedral cones in the log- arithmic parameter space. This approach is only valid when the parameters cannot change sign. If each of the quadrants in fig. 1 and 2 is mapped into the space with coordinates (logjcxl loglp() the hyperbolas and lines through the origin will be mapped into straight lines and all the stability domains (except region I1 of fig.2) will appear as convex wedges with distortions near their apexes. The cone-like structure of the stability domains breaks down near loglctl = loglpl = 0 because the parameters in eqn (2) are then all of similar magnitude. For chemical networks it can be shown that a complete breakdown of this approximation only occurs when all reactions have the same steady state currents and all evolving reactants have the same steady state concentration. The conditions for even a partial breakdown of this approx- imate approach are not expected to be encountered often. If the stability analysis is carried out for a particular reaction network rather than for the general matrix as in eqn (2) the parameters in the matrix will all be of the form nj/c where c is a steady state concentration and j is a current and n is a small integer which comes from kinetics and stoichiometry.If diffusion is included the parameters Dk2 also occur. Since these parameters are all positive the stability B. L. Clarke J. Chem. Phys. 1974,60,1481. GENERAL DISCUSSION domains in the analogue of fig. 1 are only of interest in the positive orthant and only one mapping into the logarithmic space is necessary. Furthermore each convex cone of instability can be related to the network features which produce it. The dominant term in the stability polynomials which is responsible for destabilizing a network is a product of matrix elements and can be factored into a set of feedback loops which are important in the instability.Thus each cone of instability is caused by recognizable combinations of feedback loops and the stable regions adjacent to each cone are stabilized by other recognizable network features. Dr. M. Herschkowitz-Kaufman (Brussels)said In the paper by Balsev and Degn it is stated that the existence of multiple stable patteriis is not yet established for a reaction-diffusion system. I should like to emphasize that the multiplicity of stable steady state structures has now been firmly established both by computer simulation and analytically. Consider the trimolecular reaction scheme. In a bounded medium subjected to zero flux or fixed boundary conditions computer simulations have shown the existence beyond a certain distance from the instability point of afinitr space (arbitrary units) space (arbitrary units) (4 (b) L J 0 space (arbitrary units) (4 FIG.1.-Steady state profiles of X obtained for the same values of the parameters but different initial conditions.The curves have been established for zero flux boundary conditions. Numerical values used A = 2 L = 1 B = 4.6 D = 1.6~ D = 8.0~ multiplicityof spatial steady state dissipative structures each one having its own domain of attraction as a function of the initial conditions. This is illustrated on fig. 1. Each one of the three different stable structures observed corresponds to a different G. Nicolis Proc. AMS-SIAM Symposium on Appl. Math. vol. 8 (American Math.SOC. Providence 1974). GENERAL DISCUSSION wavelength and amplitude. They emerge from the honiogeneous steady state for the same values of the parameters of the system but different initial conditions i.e. depending on the sign location and number of initial perturbations. This phenome- non can be inderstood on the basis of a bifurcation analysis which can provide an analytical approximation for the new steady state solutions bifurcating from the n FIG.2.-Real part of the eigenvalues wIIcorresponding to the unstable modes emerging from the thermodynamic branch for different values of B (zero flux or fixed boundary conditions). The other parameters have the same values as in fig. 1. homogeneous steady state not only in the neighbourhood of the first instability point but in the neighbourhood of any marginal stability point.2 However at the present time the stability of the successive bifurcating solutions cannot be guaranteed a priori.Computer simulations indicate that those solutions are stable which correspond to the fastest growing unstable modes according to the linear stability analysis (i.e. to the space(arbitrary units) space(arbitrary units) (4 (b) FIG.3,-Degenerate steady state dissipative structures obtained for the same values of the parameters by a single initial perturbation of the homogeneous steady state of the same strength but opposite sign. The arrows show for some points the sign of the perturbation leading to the corresponding spatial distribution. The numerical values of the parameters are the same as in fig.1. The homogeneous steady state concentrations X = A Y = B/A are maintained at the boundaries. first few successive bifurcations). This is seen explicitly in fig. 2 where n = 7 8 and 9 appear indeed as the leading modes for the numerical values corresponding to fig. 1. An additional mechanism of multiplicity operating even near the first bifurcation point is degeneracy as appears from the analytical expression (3.4) presented in the G. Nicolis and J. F. C. Auchmuty Proc. Nut. Acud. Sci. 1974 71 2748. GENERAL DISCUSSION 27 I paper by Nicolis and Prigogine. Fig. 3 gives a typical example of such degenerate states obtained by computer simulations. Dr. P.Hanusse (CNRS Talence) said I wish to mention to Balslev that we have established several theorems on the occurrence of instability especially in two corn- ponent systems.* The main result is that in a system with two components when only mono and bi-molecular steps are involved whatever the form of reaction scheme 1 no unstable node or focus may occur so no limit cycle type oscillation 2 when the system is stable with respect to homogeneous perturbations that is stable when stirred it can not be spatially unstable. Consequently models with two components can lead to oscillations only if they contain tri-molecular steps as does the Brussellator. Dr. J. F. G.Auchmuty (Brussels) said I would wish just to mention to Ross that for single component reaction-diffusion systems one often gets the single front propagation described in his paper in which one phase annihilates the other.In multicomponent systems however other phenomena besides coexistence and pro- pagating fronts can occur. These phenomena include the existence of fronts whose propagation velocity is not constant and of unstable propagating fronts as well as the nonexistence of such fronts. In such examples the boundaries between the states can become very com- plicated (as for example often occurs in the Stefan problem) and for certain model systems one can even get explosions when two phases are mixed in the manner described in this paper. Even the very concept of boundary between two phases is very difficult to define when the steady state solutions are space-dependent. Prof.J. Ross (MZT)said In the concluding paragraphs of our paper we comment on the agreement of two methods for analyzing aspects of nucleation one based on solutions of the time-dependent macroscopic (deterministic) equations the other on the steady state probability distribution obtained from a master equation. Since submission of this article we (K. Kithara H. Metiu and J. Ross) have found reasons for such agreement. In brief we have derived by path-integral methods solutions of the master equation for " curl-free " reaction systems and a certain form of the transition probability which show that the deterministic equations yield the trajectory for the most probable dynamic evolution. Hence the evolution of such non-linear system in the direction of increasing probability is on the average given by the deterministic equations.Dr. P. Fife (Tucson Arizona) said These remarks extend the results given by Nitzan Ortoleva and Ross in their Section 2. 1. Consider the case of a single dependent variable tl/. Let $*(A) and t.,bl(A) be the two stable stationary states and P,Hanusse Compt. Rend. l973,277C,(63 1247.l972,247C ,; GENERAL DISCUSSION If there is onIy one stationary state (unstable) between $o and 11/1 then the critical values R,(2.4) can be found from the simple relation l* J(RM) = 0. If there is more than one stationary state between $o and $ 1 then the extra condition F(IC/ nM)d$ > 0 for $O(AM) < k < #l(RM) 'hO(AM) has to be imposed. In the case of more than one dependent variable apparently there is no such simple criterion available for determining lbM.2. Again consider one dependent variable but an inhomogeneous medium with the constraints A held fixed. Thus dependence on x replaces dependence on 2. Define t,bo(x),$l(x) and J(x) as above. The coexistence problem is to find a solution of D9+F($ x) = 0, dx2 ~ with exhibiting a phase transition from valuesnear t,bo(x)to values near t,bl(x) as x crosses some transition region. If D is small a singular perturbation analysis '9 may be applied to determine the conditions for such a solution to exist. Suppose there is a value x* such that J(x*) = 0 and J(x) changes sign as x passes through x* (also assume an extra condition like that above if there is more than one stationary state between $o and Then a phase transition may occur in an interval centred at x* with width of order D*.Furthermore a heuristic argument shows this coexistence solution to be stable provided that J(x) > 0 for values of x "in phase lcfo "; J(x) < 0 for values of x "in phase ". This contrasts with the homogeneous case treated in Ross's paper in which a small perturbation of 1 will start the front moving with consequent annihiIation of one of the phases. If J is negative for some finite interval a < x < b only then can exist in that same interval with transitions to t,bo at both ends. 3. Consider now the case of two dependent variables u and v and a homogenous medium.2 Assume the diffusion coefficient D of u to be much smaller than that of u which we set equal to 1.The steady-state equations are d2u D -= f(zr 0) dx2 d2v -= g(u u). dx2 Suppose the curve f = 0 is S-shaped with t,bo and $1 the two stable branches P. Fife J. DiY. Qn. 1974 15 77. See esp. pp. 102 and 103. P. Fife J. Math. Anal. Appl. 1975 to appear. 'P. Fife and W. M. Greenlee Uspeki Mafem. Nuuk SSSR 1974 49 103 ; also to appear in Russ. Math. Surveys. A. B. Vasil'eva and V. A. Butusov Asymptotic Expansions of Solutions of Singularly Perturbed Equations (in Russian) (Nauka Moscow 1973). GENERAL DISCUSSION 1 V Again J(v) is defined and there will be a value v* at which J(u*) = 0 and J changes sign. Let If the reacting mixture is confined to a finite region say 0 < x < 1 with boun- dary conditions imposed on u and v at both ends we proceed by solving with the given boundary conditions for v imposed upon V.This problem is easily analyzed by phase-plane arguments despite the discontinuity on the right. If there exists a solution V(x) crossing the value v* at some point x* then one can expect a solution (u,v) of the original problem when D is small for which u(x) is close to V(x) and for which u exhibits a transition in a small interval of width of the order D* at x = x*. On one side of the interval u is approximated by $,,( V(x)),and on the other side by t,h1(V(x)). The proof needs some extra minor assumptions. Again the solution appears to be stable if J > 0 in phase $o and J < 0 in The following example is chosen for clarity rather than reference to any physical system.The function f is not S-shaped but the above results hold in any event. The fact that negative concentrations are assumed can easily be remedied. Primes denote differentiation in x. Example Du” = (u-u)(u2-1) v’’ = -au u’ = v’ = Oatx = Oand 1. Here $,,(v) = -1 ql(u) = 1 u* = 0 and V(x)is piecewise parabolic. On the other hand if the medium is infinite (or large compared to the square root of v’s diffusion coefficient) we niay seek a solution passing from one stable GENERAL DISCUSSION stationary state (uo,uo) to another (ul q). It is reasonable to assume that u0 = t,b0(q,) u1 = @,(u,). Then the problem reduces to finding a solution V(x) of (I) exhibiting a phase transition from vo to ul.This is like the one dimensional problem in the paper presented. It may not have a solution unless there is a dependence on some parameters A,and d is chosen appropriately. 4. The procedure in (2) and (3) above can be carried out formally when there is more than one space dimension and the region is bounded. One sometimes obtains a cell of one phase surrounded by the other phase (nucleation). The mathematical details and proof for case (2) are given in ref. (3). 5. The nucleation problem discussed in the presented paper is such that when D is small the radius of the cell Rc is of the order D3 as is the transition layer width 1. The cells constructed as in (4) however have size independent of D for small D. Dr. J. S. Turner (Austirz Texas) said I would like to point out that nucleation by spontaneous localized fluctuations in systems exhibiting multiple steady states has been studied by Lefever Prigogine and myself.We have discussed the occurrence of metastable states in a simple model and have obtained a critical coherence length for spontaneous fluctuations necessary to form an unstable nucleus in an initially homogeneous metastable phase. For details see my discussion remark following the paper of Nicolis and Prigogine in this Symposium.

ISSN:0301-5696    

DOI:10.1039/FS9740900263

出版商:RSC    

年代:1974

数据来源: RSC

摘要:

AUTHOR INDEX* Aarons L. J. 129. Auchmuty J. F. G. 98 265 271. Balslev I. 233 266. Babloyantz A. 162 263 264 267. Blank M. 218. Boiteux A.. 202,221 Bond J. R. 97 156. Burger M. 28. Busse H. G. 94. Clarke B. L. 79 80 83 161 162 215,268. Collins M. 77. Degn H. 233,266. DePuy P. E. 47. Ellis R. 153. Feinberg M. 163. Field R. J. 21. Fife P. 271. Franck U. F. 137. Friedrich U. 28. Gilbert R. G. 77 153 267. GoIdbeter A. 222. Goodwin B. C. 226 265. Gray B. F. 129 150 159 161. Gray P. 103 151 152 153 158. Griffiths J. F. 103. Hanusse P. 158 271. Hess B. 94 202. Hiernaux J. 267. Jones G. T. 192. Karvaly B. 182. Kaufman-Herschkowititz M. 70 269. Keleti T. 219. Kirsch L. J. 150 153. KGriis E. 28 92 94 96 99.Ladinyi L. 28. Lefever R. 160,215. Lewis T. J. 192. Liebhafsky H. A. 55. Mason D. M. 47 101. IMatsUzaki I. 55 101. McPhail S. 267. Meares P. 166,218. Moule R. J. 103. Nagy Zs. 28. Nakajima T. 55. Nicolis G. 7 68 69 70 74 98 263 265. Nitzan A. 241. Noyes R. M. 21 78 79 81 82 84 86 87 89 93 97 101 102 154 158. Orbrin M. 28. Ortoleva P. 69 83 97 241 263 264. Page K. R. 166. Perche A. 157. Plesser Th. 225. Portnow J. 72. Prigogine I. 7 68 69 70. Rapp P. 215. Ross J. 66 68,69 82 83,97,241,263,264,271. Rossler 0. E. 87 91 97. Shashoua V. E. 174. Smoes M-L. 73 74 85 86 98. Ssrensen P. G. 88. Texter J. 254. Tributsch H. 164 217 221. Turner J. S. 75 264 274. Winfree A, 38,83 100,224,266. Yang C. H. 114 153 154 155 157. *The references in heavy type indicate papers submitted for discussion. 275

ISSN:0301-5696    

DOI:10.1039/FS9740900275

出版商:RSC    

年代:1974

数据来源: RSC