For a model glycolytic diffusion-reaction system, an amplitude equation has been derived in the framework of a weakly nonlinear theory. The linear stability analysis of this amplitude equation interprets the structural transitions and stability of various forms of Turing structures. This amplitude equation also conforms to the expectation that time-invariant amplitudes in Turing structures are independent of complexing reaction with the activator species, whereas complexing reaction strongly influences Hopf-wave bifurcation.
I. INTRODUCTION
Nonlinear diffusion-reaction (D-R) systems exhibit self-organized concentration patterns when the steady states remain stable to homogeneous perturbations, but become unstable to inhomogeneous perturbations. These patterns are called Turing structures.1 The inhibitor species must diffuse sufficiently faster than the activator species, which facilitates Turing structure formation. The structural transitions and the stability of various forms of Turing structures can be well-interpreted from the linear stability analysis of an amplitude equation (AE), which can be derived in the framework of a weakly nonlinear theory.
The AE derivation is, itself, a lengthy process. Only those D-R systems have been chosen in the past for AE derivation, for which the kinetic steady state is analytically solvable as well as simple-looking in algebraic form. As for example, AE2,3 for Turing pattern selection4–6 has been derived in the past for a few D-R systems – examples include the Brusselator7,8 and other models.9–13 AE derivation for the Swift-Hohenberg model14 has been widely used for a qualitative description of the convective structures originated from Benard-Marangoni instability or non-Boussinesq Benard convection.15,16
In this paper, we have derived the AE for the reversible Sel'kov model17,18 a kinetic model of glycolytic19,20 D-R system, in which the steady state from the reaction kinetics is required to be obtained by the analytical solution of a cubic equation using Cardon's method,21 which is rather complicated. For such a D-R model,17,18 we have attempted to derive the AE, which interprets the stability of various forms of Turing patterns as well as the structural transitions between them. In a recent paper,22 we have made a preliminary report of this derivation, but could not present in detail, particularly the nonlinear analysis and the amplitude equation derivation steps after incorporating the two solvability conditions – this has been accomplished in this manuscript. Also, a linear stability analysis of the amplitude equation together with some numerical results on Hopf and Turing bifurcations has been reported in detail, which we could not present before.22 The other related sections have been included in this manuscript in condensed form.
Section II has discussed the reversible Sel'kov model in brief. Section III has presented the results of the linear stability analysis of the D-R equations (2.4) and (2.5). Section IV has undertaken the decomposition of the D-R equations into linear and nonlinear parts as shown in Eqs. (4.1). Section V has discussed the nonlinear analysis – this is based on the expansion of the evolution of the inhomogeneous perturbations in terms of the amplitude of the branching solution for the bifurcation parameter in the neighborhood of its critical value at the Turing bifurcation point, and the two solvability conditions from second and third order terms have been obtained. In section VI, the two solvability conditions have been incorporated into the expanded form of the partial time derivative of amplitude (see Eq. (5.5)), which has finally derived the AE. In section VII, a linear stability analysis of the AE has been undertaken and has discussed how Turing structure selection problem can be interpreted from the linear stability analysis of the amplitude equation. Section VIII has reported numerical results on Hopf and Turing bifurcations. The new results have been summarized in section IX.
II. THE MODEL
A simple kinetic model of enzyme catalysis with product activation of the enzyme has been proposed by Sel'kov,19 which exhibits limit cycle oscillations. In this model, a substrate S (ATP) is supplied by a certain source at a constant rate (v1), which is converted irreversibly into a product P (ADP). The product P (ADP) is removed by an irreversible sink at another constant rate (v2.). The free enzyme E (PFK1), itself, is inactive but becomes active in the complex form EPγ after combination with γ product molecules. Richter et al.17 have proposed the following three reversible kinetic steps (2.1 to 2.3), based on Sel'kov model, known as the reversible Sel'kov model,
where, A and B are source and sink concentrations respectively, which could be controlled from outside — k±i (i = 1, 2, 3) are, respectively, the rate constants of the three steps (+ and – subscripts are, respectively, for forward and reverse reactions). It is required to note that the reversible Sel'kov model steps account for only the oscillatory part of the mechanism - the complete glycolytic mechanism, itself, is rather complicated. The autocatalyst P is the activator and S, the inhibitor, when these are described in terms of activator/inhibitor model. The D-R equations18,20 in dimensionless form are given by
where,
where, |$\hat p$| and |$\hat s$| are scaled concentrations, ρ, the scaled geometrical coordinate and τ, the scaled time as given below.
where, P and S represent the concentrations of the respective species, D's are diffusion coefficients, x is the geometrical coordinate and K′ is a measure of the degree of complexing reaction with the activator species P (ADP) such that
where, K is the equilibrium constant for the complex formation reaction.23,24
where, the activator P is captured partially producing the complex PC by reaction with the chemical species C. Also, the complexing agent is assumed to be in large excess, such that the initial concentration (co) of the complexing agent is almost equal to its concentration (c) at chemical equilibrium represented by equation (2.11b). The possible nature of the complexing species is an unexplored problem which is left open for future investigation
III. LINEAR STABILITY ANALYSIS
From the dimensionless D-R equations (2.4) and (2.5), one obtains the elements of the Jacobian matrix given as
where, so, po are the dimensionless steady states. The characteristic equation for the dimensionless partial differential equations (2.4) and (2.5) may be written as:
where ωk′s are eigenvalues, k is the wavenumber and D's are diffusion constants.
From the characteristic equation (3.2) for nonzero k mode, the condition for Hopf-wave bifurcation is given by
For Turing instability to set in,25,26 the constant term of the characteristic equation (3.2) must vanish.
Turing instability condition in terms of the dimensionless steady state concentrations of P and S is given by,
Equation (3.4), on algebraic manipulation for κ = κcbecomes
which gives the value of κc to be satisfied at the Turing bifurcation point.
The critical wavenumber (kc) at the Turing bifurcation point, which maximizes the constant term of Eq. (3.2) may be given by the expression
IV. SEPARATION INTO LINEAR AND NONLINEAR PARTS
One obtains the D-R equations (2.4) and (2.5) divided into linear and non-linear parts as given below.
where
V. NONLINEAR ANALYSIS
Let us consider the simplest possible case that the Turing structures constitute the superposition of three pairs of modes (|$\vec k_i, - \vec k_i)_{i = 1,2,3}$| making angles of 2π/3 between each pair, such that ||$\vec k_i$|| = kc and |$\sum\nolimits_{i = 1}^3 {\vec k_i }$| = 0: the inhomogeneous perturbations, (p, s)T unfold in terms of three amplitudes (|$A_i ^{\prime} s$|) according to the relation given by
where c.c. indicates the complex conjugate terms. In nonlinear analysis the following expansion is made in terms of the amplitude of the branching solutions for κ (say) in the neighborhood of κc (the critical value for Turing bifurcation).
Such that,
where A is the amplitude. Substituting Eqs. (5.2) to (5.4) into Eq. (4.1), one obtains first, second and third order balance equations as given below in Eqs. (5.6)–(5.8) respectively.
Order 1:
From Eq. (5.6), (p1, s1)T for the critical mode (at κ = κc, ||$\vec k_i$|| = kc) is given below (for calculations see Appendix C)
Order 2:
First, the right eigenvectors of LT with zero eigenvalue is to be obtained from Eq. (5.10) given below (see Appendix D for its solution)
This right eigenvecor, on transposing, gives the left eigenvector with zero eigenvalue as given below in Eq. (5.11).
We incorporate Eq. (5.9) into Eq. (5.7) and calculate (Fp1, Fs1)T, the coefficient of exp (i|$\vec k_1 \cdot \rho$|) in the form as given in equation (5.13).
Therefore,
Eqs. (5.11) and (5.13) give the first solvability condition as given in Eq. (5.14) [see Appendix E for its derivation).
Now, we express (p2, s2)T in Eq. (5.7) as the sum of a series of Fourier terms as follows,
The Fourier coefficient values have been calculated in Appendix F.
Order 3:
The coefficient of exp (i|$\vec k_1 \cdot \rho$|) i.e., (Gp(1), Gs(1))T from Eq. (5.8) is given as
Substituting the values of the Fourier coefficients from Appendix F into Eq. (5.16) and making long algebraic manipulations, one obtains
where, C1 and C23 are the expressions included in Appendix G. The second solvability condition from equation (5.8) will be such that the inhomogeneous term (Gp(1), Gs(1))T be orthogonal to the left eigenvectors of L with zero eigenvalue. Therefore, the second solvability condition should satisfy the relation:
Substituting the values of Gp(1) and Gs(1) from Eqs. (5.17) into Eq. (5.18), one obtains the second solvability condition as given in Eq. (5.19).
VI. AMPLITUDE EQUATION
The amplitude A1 from Eq. (5.5) is given by
The term (∂A/∂|$\tau ^{\prime} _0$|) from Eq. (5.5) is neglected since the slow variable at zero time is time-invariant. Substituting the values of (p1, s1)T and (p2,s2)T from Eqs. (5.9) and (5.15) respectively into Eq. (5.2), we equate the coefficients of exp(i|$\vec k_i \cdot \rho$|) from equations (5.1) and (5.2). One obtains using Eq. (F11) from Appendix F,
Substituting the value of A1 from equation (6.3) into equation (6.1), one obtains (neglecting terms of order higher than ε3),
We incorporate the solvability conditions, Eqs. (5.14) and (5.19), into equation (6.4) to obtain the amplitude equation as given below in Eq. (6.5) with the help of equations (5.3) and (6.3) – neglecting all terms involving ε including those of order ε3 and higher (Derivation of Eq. (6.5) including the coefficients τo, μ,|$\bar v$|, g, and h are given in Appendix H).
We can obtain similar equations for A2 and A3 by circular permutation of the indices from Eq. (6.5).
VII. LINEAR STABILITY ANALYSIS OF THE AMPLITUDE EQUATION
We make a linear stability analysis of the amplitude equation:12,27,28 as follows. Amplitudes in Eq. (6.5) may be expressed as
where,
is the phase. Therefore, one obtains from Eqs. (7.1) and (7.2),
Considering the real phase only, we have
Considering the real phase ϕ only, one obtains
We have considered both the positions29 ϕ = 0 and ϕ = π for which the value of |$\bar v$|becomes positive (>0) and negative (<0) respectively. Therefore, Eq. (7.6) takes the form
Equation (7.7) may possess three types of solutions:
A. Stability of Stripe Structures
The lower limit of the stability of stripe structures is obtained by a linear stability analysis of Eq. (7.7) around Rj = Rs + δRj (j = 1, 2, 3 respectively). Substituting R1 = Rs + δR1, R2 = δR2, and R3 = δR3 in Eq. (7.7) for the steady state value Rs, one obtains (neglecting terms of order δ2):
Substituting the steady state condition from Eq. (7.7), one obtains
Therefore, the characteristic equation governing the stability of stripe structures takes the form,
where, λ's are the corresponding eigenvalues. Substituting Rs = (μ/g)1/2 for stripe structures from Eqs. (7.8), one obtains
Therefore,
Stable stripe structures will grow, only if all the eigenvalues are negative. This implies
Therefore, the stripe structures are stable to all infinitesimal perturbations for the bifurcation parameter μ > μS as given in Eq. (7.14).
B. Stability of Hexagonal Structures
Substituting R1 = R2 = R3 = RH in Eq. (7.7), one obtains
For real values of RH±, the discriminant must be positive, which gives the lower limit of stable hexagonal structures in the form
The upper limit of the stability of hexagonal patterns is calculated by a linear stability analysis of Eq. (7.7) around Rj = RH + δRj (j = 1, 2, 3 respectively). Substituting R1 = RH + δR1, R2 = RH + δR2, and R3 = RH + δR3 into Eq. (7.7), one obtains (neglecting terms of order δ2 or higher):
where,
Substituting the steady state condition from Eq. (7.7) into Eq. (7.17), we have
The characteristic equation from Eqs. (7.19a), (7.19b), and (7.19c) for the amplitude variations δR1, δR2 and δR3 is given by
where, λ′s are the eigenvalues. The solution of this cubic equation (see Appendix I) is given by
where, a1 and a2 values are those given in Eq. (7.18). Let us now investigate the signs of the eigenvalues for the two cases separately namely RH = RH+ and RH = RH− respectively (see Eqs. (7.15))
Case RH = RH+:
Substituting the value of RH+ from Eqs. (7.15), one obtains
Now, we find out the condition for which the other two equal eigenvalues (λ2, λ3) become negative to ensure the stability of hexagonal structures. Imposing that λ2 = λ3<0, we have
which, on simplification gives
Therefore, the hexagons for RH = RH+ are stable if μ < μH2.
Case RH = RH−:
The expression given by Eq. (7.29) becomes positive by comparison with Eq. (7.24). Now, all the three eigenvalues for the case RH = RH− become positive ensuring that the hexagonal structures are unstable in this case.
The bifurcation scenario may be summarized schematically in Fig. 1 in the form of amplitude as a function of the bifurcation parameter μ (the normalized distance from onset) for stripe state(S) and hexagon state (H) respectively. For Turing pattern amplitude equation (6.5), the stable hexagonal branch is Ho (or Hπ) if |$\bar v$| > 0 (or |$\bar v$| < 0). A subcritical hexagonal branch starts at μH1 < 0, and remains stable until μH2>0. The supercritical stripe branch is unstable near the critical point, but it becomes stable when μ > μS. Both branches (stripes and hexagons) are stable in the range μS<μ <μH2(bistable region).
VIII. RESULTS AND DISCUSSION
Equation (3.3) indicates that Hopf bifurcations of the homogeneous (k = 0) and inhomogeneous (k ≠ 0) model systems are strongly affected by the complex formation with the activator P – the Hopf boundary would be shifted to a lower value of κ if K’ has values greater than zero. Hopf wavenumber (kH) as a function of a for four different values of d have been reported in Figure 2 of reference 24 which suggests that, due to the complexing reaction of the activator (K’ = 0.5), there is an approximately two-fold decrease (increase) of Hopf-wavenumber (wavelength) values. Also, the areas of phase diagrams in the parameter space (see Figs. 2 and 3 in the following paragraphs) and the corresponding maximum values of the possible wavenumbers (shortest wavelengths) are reduced drastically at high values of d and K’. When translated to real lengths, our model calculations have predicted24 inhomogeneous Hopf-waves in the order of a few millimetres in length (between 3 and 0.6 mm for Ds ≈ 10−5cm2s−1 and |$k_{^3} \approx 10\, s^{-1}$|). However, reliable experimental data are not currently available for direct comparison with the computed values.
It is interesting to note in figures 2(a)–2(c) that, for a very low value of d = 2.5, the Hopf-domain in the parameter space decreases with the increase of K’ (degree of complex formation with the activator P), but the tiny Turing region remains practically unchanged in size. Since a very low value of d, for the parameters chosen, does not help Turing domain expand for K’ > 0, a gap between these two domains is noticeable in figures 2(b) and 2(c). This gap of separation between Hopf and Turing domains at low values of d seems to be a new phenomenon which has not been reported before – but, logically not unexpected since the inhibitor species must diffuse much faster than the activator species to facilitate Turing structure formation. For large values of d ≥ 10, however, the Turing domain expands in the parameter space at the cost of the Hopf domain for K’ > 0, and apparently there exists no gap of separation between these two domains as presented in figures 3(a)–3(c) – these figures demonstrate that, while the parameter κis being decreased gradually, Turing region precedes the Hopf bifurcation line in the (a - κ) phase plane both in the absence and presence of complexing reaction with the activator.
But the complex formation reaction with the activator species P neither changes the wavenumber (wavelength) of Turing patterns (see equation (3.6)) nor does it affect the basic relation to be maintained for Turing instability to occur (see equation (3.4)). Turing wavenumber kc as a function of a is presented in figures 4(a) and 4(b) for different values of d (see equation (3.6)) in the absence and presence of the complexing reaction of the activator, respectively, to verify that kc, indeed does not depend on K’. It is interesting to note from these two plots that K’ = 1 does help to obtain Turing wavenumbers (kc) at a very low value of d = 5 which, otherwise could not produce any Turing structures in the absence of any complexing reaction with the activator. This result together with those presented in figures 2(a)–2(c) is in full agreement with the observations of De Kepper et al.6 b on the chlorite-iodide-malonic acid (CIMA) reaction with varying concentrations of starch. Also, by increasing K’ in Eq. (3.3) to an appropriate value, it is possible in principle to arrest the arrival of Hopf bifurcation in this model system and the Hopf domain, which will disappear by this technique, may be utilized to generate Turing patterns (see figures 2 and 3) by inducing inhomogeneous perturbations of nonzero k mode.
Different higher values of d can be realized in experiments by adopting a regulatory mechanism based on immobilization of the activator P in a real glycolytic system. The computed Turing wavenumbers in figures 4(a) and 4(b), when translated to real unit of length, give the Turing wavelength values in the range 0.21– 0.15 mm, which is left for comparison with the real experiments of glycolytic Turing patterns in the future. The linear stability analysis27–31 of the amplitude equation (6.5), as presented schematically in Fig. 1 in the last section interprets the structural transitions and stability of stripes and hexagonal patterns.
IX. CONCLUSION
The present manuscript has reported an amplitude equation (6.5) for a model D-R system of glycolytic oscillations in the framework of a weakly nonlinear theory. This amplitude equation derivation is also interesting from mathematical point of view, since the analytical steady state solution of this model system is rather complicated21 in algebraic form, whereas only the D-R systems with kinetic steady states simple-looking in algebraic form have been chosen in the past2,3,7–13 for amplitude equation derivation.
Figure 2 of reference 24 and Eq. (3.3) indicate that complexing reaction with the activator species P strongly influences the wavenumber (wavelength) of Hopf-wave bifurcations, whereas Fig. 4 and Eq. (3.6) show that Turing pattern wavenumber (wavelength) is independent of the complexing reaction with the activator species P. We can obtain the above-cited computed results immediately just by having a look at the amplitude equation (6.5) itself. Time dependent amplitudes in Eq. (6.5), such as in Hopf-wave bifurcations strongly depend on K’ (degree of complexing reaction with the activator species P), whereas time independent amplitudes such as in Turing bifurcations are independent of the degree of complexing reaction with the activator species P : K’ is included only in the expression for τo.
Complexing reaction arrests the arrival of Hopf bifurcation and the Hopf domain which will disappears due to complexing reaction may be utilized for Turing structure generation by inducing inhomogeneous perturbations (of nonzero k mode) – this is illustrated in Fig. 3. We obtain a big gap of separation between Hopf and Turing domains if the value of the ratio d (equal to Ds/Dp) is not large enough – this is well illustrated in Fig. 2 for d = 2.5, which could generate only a tiny Turing domain even in the presence of complexing reaction.The collective message from Figs. 2 and 3 is – unless the ratio d is sufficiently large, complexing reaction can not transform the lost Hopf domain to Turing domain resulting in a gap of separation between these two domains and the size of the Turing domain is left practically unchanged. Figure 1 interprets the structural transitions and the stability of stripes and hexagonal patterns qualitatively, which is obtained from the linear stability analysis of the amplitude equation (6.5) as described schematically in Fig. 1.
ACKNOWLEDGMENTS
I am thankful to Dr. J. Boissonade of Bordeaux for helpful suggestions. This research was supported in part by EPSRC (U.K).
APPENDIX A
Let the inhomogeneous perturbation (p, s) around the stable steady state (po, so) be represented by the following relations
Since (po, so) is the dimensionless steady state solution of the D-R equations (2.4) and (2.5), we have
Substituting the values of (|$\hat p$|, |$\hat s$|) from equations (A1) into the D-R equations (2.4) and (2.5), one obtains the D-R equations in the form as given below.
Therefore,
Incorporating the steady state conditions, Eqs. (A2) into Eqs. (A5) and (A6), one obtains the D-R equations divided into linear and non-linear parts as given below in Eq. (A8), which derives Eq. (4.1).
APPENDIX B
Equating the coefficients of the order of ε, ε2, and ε3 on both sides, one obtains,
which derives Eq. (5.6).
Therefore,
which derives Eq. (5.7)
which derives Eq. (5.8)
APPENDIX C
We have
Let. |$(\bar p_{1,}$||$\bar s_1$|)T be the solution for non-degeneracy such that
Since ||$\vec k_i$|| = kc, we have
Substituting the value of |$k_c ^2$| from Eq. (3.6), we have
which gives the solution for non-degeneracy as
For degeneracy as given in Eqs. (C3), we obtain the solution (p1,s1)T as shown in Eq. (5.9). We could also get the same solution starting from the second equation (C2) and substituting the values of κc and |$k_c ^2$| from Eqs. (3.5) and (3.6) respectively.
APPENDIX D
We have
Substituting the values of L and κc from Eqs. (4.2) and (3.5) respectively into Eq. (5.10) and considering that |$(\bar p_{2,}$||$\bar s_2$|)T be the solution for non-degeneracy such that
one obtains Eq. (5.10) divided in component forms as given below.
Using ||$\vec k_i$|| = kc, Eqs. (D3) and (D4) give
which gives the solution
APPENDIX E
Since the left eigenvector with zero eigenvalue given by Eq. (5.11) is orthogonal to (Fp(1), Fs(1))T as given by Eq. (5.13), we have
Therefore,
Substituting the values given in Eqs. (5.11) and (5.13), one obtains
Or,
Or,
which derives the first solvability condition given by Eq. (5.14).
APPENDIX F
Calculation of (Po, So)Texp (0):
Substituting Eq. (F1) into Eq. (5.7) and equating the coefficients of exp(0) (i.e. the constant term) on both sides, one obtains
Solving Eqs. (F2) and (F3), one obtains
where,
Calculation of (P1, S1)T:
Substituting Eq. (F8) into Eq. (5.7) and equating the coefficient of exp (i|$\vec k_1 \cdot \rho$|) from both sides, one obtains
Solving Eqs. (F9) and (F10), one obtains
Calculation of(P11, S11)T:
Substituting Eq. (F12) into Eq. (5.7) and equating the coefficient of exp (2i|$\vec k_1 \cdot \rho$|) from both sides, one obtains
Solving Eqs. (F13) and (F14), one obtains
Calculation of(P12, S12)T:
Substituting Eq. (F17) into Eq. (5.7) and equating the coefficient of exp (i|$(\vec k_1 - \vec k_2) \cdot \rho$|) from both sides, one obtains
Solving Eqs. (F18) and (F19), one obtains
APPENDIX G
and
APPENDIX H
Substituting the value of (∂W1/∂|$\tau ^{\prime} _1$|) from the first solvability condition, Eq. (5.14) into Eq. (H1) and using Eq. (6.3), one obtains
Dividing both sides byκc, one obtains
Neglecting all the terms involving ε including those of order ε3and higher, one obtains
where,
which, is a normalized distance from onset,
APPENDIX I
The equation (7.20) is of the form
where A = 1, B = −a1, C = a12 − a22, D = 3a1a22 − a13 − 2a23
The cubic equation (I1), when reduced to the standard form by the transformation Z = Aλ + B, takes the form
where H = AC − B2, G = A2D − 3ABC + 2B3. As a solution of Eq. (I2), we suppose that
Comparing Eqs. (I2) and (I3), one gets
Thus m3and n3are the roots of the equation
such that
Since there are three cube roots of any number, let them in this case be m, mω, mω2 and n, nω, nω2 respectively. Any value of m may be associated with any value of n giving altogether nine values of Z = m + n. But, these must be so chosen that their product is real as suggested by mn = −H. Therefore, we get three pairs of admissible values of m and n − (m, n), (mω, nω2), (mω2, nω), since the product of any other pair is imaginary. Thus, we get three roots of Z of Eq. (I2) given as
From the relation Z = Aλ + B, the roots of Eq. (I1) are given by
which, on substitution of the values of A, B, m and m takes the form (a1 + 2a2), (a1 − a2), (a1 − a2) respectively.