Sufficient conditions for wave instability in three-component reaction-diffusion systems
aa r X i v : . [ n li n . PS ] F e b Sufficient conditions for wave instability in three-component reaction-diffusion systems
Shigefumi Hata, Hiroya Nakao, and Alexander S. Mikhailov Department of Physical Chemistry, Fritz Haber Institute of the Max Planck Society, Faradayweg 4-6, 14195 Berlin, Germany. Department of Mechanical and Environmental Informatics,Tokyo Institute of Technology, Ookayama 2-12-1, 152-8552 Tokyo, Japan
Sufficient conditions for the wave instability in general three-component reaction-diffusion systemsare derived. These conditions are expressed in terms of the Jacobian matrix of the uniform steadystate of the system, and enable us to determine whether the wave instability can be observed as themobility of one of the species is gradually increased. It is found that the instability can also occurif one of the three species does not diffuse. Our results provide a useful criterion for searching waveinstabilities in reaction-diffusion systems of various origins.
The wave instability provides an important mechanism for pattern formation in nonequilibrium chemical systems.When it takes place, a critical mode corresponding to a traveling wave with a certain wavenumber and oscillationfrequency begins to grow, destabilizing the uniform steady state. Although being less known, the wave instabilityhas already been considered in 1952 by A. Turing in his pioneering publication [1], where the classical (i.e. static)Turing instability, leading to the establishment of a periodic stationary pattern has also been introduced. Therefore,it may be also appropriate to describe it as the oscillatory Turing bifurcation. Moreover, it was already noticed byA. Turing [1] that at least three interacting species are needed for this instability to occur.Because of the spatial reflection symmetry, waves traveling in the left and right directions have the same growthrates and both of them begin to spontaneously develop above the instability threshold. Nonlinear interactions betweensuch modes determine whether one of the modes gets suppressed, so that a wave traveling in a certain direction isestablished, or standing waves, representing superpositions of left and right traveling waves, are instead formed [2, 3].The wave patterns resulting from such instability can also exhibit further instabilities, and wave turbulence may seton.In contrast to the classical Turing bifurcation, which has been extensively discussed for both biological and chemicalsystems [1, 4–8], the wave bifurcation has so far attracted less attention. It has been considered for special chemicalmodels [9] and its existence was suggested in the experiments with Belousov-Zhabotinsky microemulsions [10]. Thereare also publications in which this instability was discussed for special ecological models [11].Because at least three species are needed for the wave instability to occur, the linear stability analysis is morecomplex in this case, as compared with the classical Turing bifurcation in two-component activator-inhibitor systems.The complexity of the stability analysis, which has been previously performed separately for individual chemicalsystems, has probably also been responsible for the fact that the wave instability has not been broadly investigatedfor reaction-diffusion media.In this article, we present a general derivation of the sufficient conditions for the wave (or oscillatory Turing)bifurcation in arbitrary three-component reaction-diffusion models. The final conditions are formulated in terms ofthe elements of the Jacobian matrix of the uniform steady state. They tell whether the wave bifurcation is possiblewhen the mobility of any chosen species is gradually increased, while diffusion coefficients of other species are keptconstant. As we show, the instability is possible even if one of the three species is immobile. Using our generalexpressions, the analysis of the wave bifurcation in specific reaction models will be largely simplified.We consider reaction-diffusion systems with three chemical reactants
U, V and W . Local densities of the reactantsare denoted as u = [ U ] , v = [ V ] and w = [ W ]. All reactants diffuse over the space and undergo local chemicalreactions. Generally, such systems are described by equations dudt = f ( u, v, w ) + D u ∇ u,dvdt = g ( u, v, w ) + D v ∇ v,dwdt = h ( u, v, w ) + D w ∇ w, (1)where functions f, g and h represent the local reactions. Diffusion coefficients of the reactants are D u , D v and D w .We assume that a uniform steady state ( u, v, w ) = (¯ u, ¯ v, ¯ w ) determined by f (¯ u, ¯ v, ¯ w ) = g (¯ u, ¯ v, ¯ w ) = h (¯ u, ¯ v, ¯ w ) = 0exists and that this state is stable in absence of diffusion.We introduce small perturbations to the steady state as ( u, v, w ) = (¯ u, ¯ v, ¯ w ) + ( δu, δv, δw ). Substituting this intoEqs. (1), the following linearized differential equations for the perturbations are obtained: ddt δu = f u δu + f v δv + f w δw + D u ∇ δu,ddt δv = g u δu + g v δv + g w δw + D v ∇ δv,ddt δw = h u δu + h v δv + h w δw + D w ∇ δw, (2)where f u = ∂f /∂u | (¯ u, ¯ v, ¯ w ) , f v = ∂f /∂v | (¯ u, ¯ v, ¯ w ) , f w = ∂f /∂w | (¯ u, ¯ v, ¯ w ) , ... are partial derivatives at the steady state. Thefollowing rescaled variables are introduced for convenience: δ ˜ u = δu, δ ˜ v = s(cid:12)(cid:12)(cid:12)(cid:12) f v g u (cid:12)(cid:12)(cid:12)(cid:12) δv, δ ˜ w = (cid:12)(cid:12)(cid:12)(cid:12) f v h v (cid:12)(cid:12)(cid:12)(cid:12) δw, ˜ t = p | f v g u | t. (3)We substitute these variables into Eqs. (2) to obtain a set of equations dd ˜ t δ ˜ u = m δ ˜ u + αδ ˜ v + nδ ˜ w + D u µ ∇ δ ˜ u,dd ˜ t δ ˜ v = βδ ˜ u + p δ ˜ v + qδ ˜ w + D v µ ∇ δ ˜ v,dd ˜ t δ ˜ w = rδ ˜ u + γδ ˜ v + s δ ˜ w + D w µ ∇ δ ˜ w, (4)where m = f u p | f v g u | , n = f w | h v | q(cid:12)(cid:12) f v g u (cid:12)(cid:12) ,p = g v p | f v g u | , q = g w (cid:12)(cid:12)(cid:12)(cid:12) h v f v g u (cid:12)(cid:12)(cid:12)(cid:12) ,r = h u | h v | s(cid:12)(cid:12)(cid:12)(cid:12) f v g u (cid:12)(cid:12)(cid:12)(cid:12) , s = h w p | f v g u | ,µ = 1 p | f v g u | . (5)The coefficients α, β, γ are determined by signs of f v , g u , h v , α = sign( f v ) , β = sign( g u ) , γ = sign( h v ). The perturba-tions ( δ ˜ u, δ ˜ v, δ ˜ w ) are expanded over plane waves as δ ˜ u = Z d~k ˜ U ( ~k ) exp h λ ( ~k )˜ t − i~k · ~x i ,δ ˜ v = Z d~k ˜ V ( ~k ) exp h λ ( ~k )˜ t − i~k · ~x i , (6) δ ˜ w = Z d~k ˜ W ( ~k ) exp h λ ( ~k )˜ t − i~k · ~x i , where ~k = ( k , k , · · · ) is the wave vector, and λ ( ~k ) is a growth rate of the plane wave with wave vector ~k . Thus, weobtain the following equations for each wave vector ~k : λ ( k ) ˜ U ( k ) ˜ V ( k ) ˜ W ( k ) = m − D u µk α nβ p − D v µk qr γ s − D w µk ˜ U ( k ) ˜ V ( k ) ˜ W ( k ) , (7)where k = (cid:12)(cid:12)(cid:12) ~k (cid:12)(cid:12)(cid:12) is the wave number, the magnitude of the wave vector ~k . Because only the magnitude of the wavevector is important, here and below we drop the vector symbols.The condition det m − D u µk − λ ( k ) α nβ p − D v µk − λ ( k ) qr γ s − D w µk − λ ( k ) = 0 (8)should be satisfied for Eq. (7) to have non-trivial solutions. Thus, the linear growth rate λ ( k ) is determined by thecharacteristic equation λ − ( m + p + s ) λ + ( mp + ps + sm − nr − γq − αβ ) λ − ( mps − npr + αqr − γmq + βγn − αβs ) = 0 , (9)where m = m ( k ) = m − k µD u ,p = p ( k ) = p − k µD v , (10) s = s ( k ) = s − k µD w . The growth of each plane-wave mode is determined by the real part of λ ( k ) . The uniform steady state is stable ifRe( λ ( k ) ) is negative for all k . The instability occurs if Re( λ ( k ) ) becomes positive for at least one wave number k = k c .Then the uniform steady state is destabilized, leading to spontaneous development of wave patterns with critical wavenumber k c . If the imaginary part Im( λ ( k c ) ) of the unstable mode is zero, the first critical mode represents a stationaryplane wave and the Turing instability occurs. On the other hand, if Im( λ ( k c ) ) = 0, the critical mode is oscillatory intime and periodic in space, so that the critical mode represents a traveling wave and the wave instability takes place.The complex conjugate root theorem tells that a characteristic equation with real coefficients λ + aλ + bλ + c = 0 (11)has either three real roots or one real root and a pair of complex conjugate roots. In the latter case, the three rootscan be written as λ , = ψ ± iω, λ = φ, (12)so that the coefficients are represented as a = − (2 ψ + φ ) , b = ψ + ω + 2 ψφ, c = − ( ψ + ω ) φ. (13)Combining these three equations, we obtain c − ab = − ( ψ + ω ) φ + (2 ψ + φ )( ψ + ω + 2 ψφ )=2 ψ (cid:2) ( ψ + φ ) + ω (cid:3) . (14)At the threshold of the wave instability, we have ψ = 0, so that c − ab = 0. At the threshold of the Turing instability,we would have φ = 0, and therefore c = 0. Note that the Turing instability is also possible when the characteristicequation has three real roots and one of them becomes positive. It can be easily checked that, also in this case, theinstability threshold corresponds to c = 0.Thus, the wave instability first takes place when, for one wavenumber k = k c , the equation I wav =( mps − npr + αqr − γmq + βγn − αβs ) + ( m + p + s )( mp + ps + sm − nr − γq − αβ ) = 0 (15)becomes satisfied, where m, p and s are given by equations (10) with k = k c .The Turing instability first takes place when, for one wavenumber k = k c , the equation I st = mps − npr + αqr − γmq + βγn − αβs = 0 (16)becomes satisfied, where the coefficients m, p and s are again given by equations (10) with wavenumber k c .It is convenient to introduce the three-dimensional m - p - s space in order to represent these conditions graphically.As illustrated in Figures 1 and 2, each of the conditions (15) and (16) defines a boundary surface and Eqs. (10)determine a straight line Γ which is parameterized by the wave number k . If the line Γ touches the boundary surface I wav ( m, p, s ) = 0 or I st ( m, p, s ) = 0, then the conditions (15) or (16) are satisfied and the corresponding instabilitytakes place. (a) (b) P Γ PP P Γ FIG. 1: Boundary surface for the wave instability I wav ( m, p, s ) = 0. The line Γ touches the surface at a point P . Panel (b) isa closeup of (a). Parameters are fixed at n = q = r = α = β = γ = 1, m = 0 . , p = 0 . s = − .
8. Diffusion constantsare D u = 1 , D v = 1 and D w = 14 . I st ( m, p, s ) = 0. Parameters are fixed at n = q = r = α = β = γ = 1. The uniform steady state (¯ u, ¯ v, ¯ w ) should be stable if the diffusion is absent. In terms of of the Jacobian matrix atthe steady state J = f u f v f w g u g v g w h u h v h w , (17)this implies that det( J ) < J ) < . (18)Thus, the initial point P with coordinates ( m , p , s ) on the line (10) should lie inside the stable region in m - p - s space.The wave instability takes place if the first critical mode is oscillatory. At the threshold of the wave instability, theline Γ defined by Eq. (10) should touch the boundary surface I wav = 0 without having intersections with the surface I st = 0.Suppose that we want to check whether the wave instability can occur when the diffusion constant D w is varied.Let us denote A = D v /D u and B = p − m D v /D u and consider a plane p = Am + B parallel to the s -axis. Theline Γ is always lying on this plane irrespective of D w , because the conditions p = Am + B and D v /D u = A aresatisfied. A coordinate x is introduced on the plane in such a way that we have x = 0 when m + p = 0 holds. Theslope of the line Γ in the x - s space is D u D w / p D u + D v and is nonnegative. The initial point P of the line Γ isat ( x , s ) on the plane, where x = ( m + p ) / ( A + 1). Because we assume that the necessary conditions (18) aresatisfied, the point ( x , s ) is in the stable region of the x - s space. Two surfaces I wav = 0 and I st = 0 intersect theplane along the boundary curves ˆ I wav = 0 and ˆ I st = 0,ˆ I wav ( x, s ; A, B ) = I osc (cid:18) x − BA + 1 , Ax + BA + 1 , s (cid:19) = 0 , (19)ˆ I st ( x, s ; A, B ) = I st (cid:18) x − BA + 1 , Ax + BA + 1 , s (cid:19) = 0 . (20)Below we examine the dependences of the boundary curves ˆ I wav = 0 and ˆ I st = 0 on the parameters m , p , s , n, q, r, α, β and γ . This allows us to obtain the parameter conditions under which the wave instabilitycan occur.Let us examine the shape of the boundary curve ˆ I wav = 0 at large values of | s | . If | s | ≫
1, we can neglect O ( s )terms and obtain ˆ I wav ( x, s ; A, B ) ≈ s (cid:2) (1 + A ) x + s (1 + A ) x − nr − γq (cid:3) . (21)Then, the boundary curve ˆ I wav = 0 is given by the equation(1 + A ) x + s (1 + A ) x − nr − γq = 0 (22)and we have x = 12 − sA + 1 ± s(cid:18) sA + 1 (cid:19) + 4 nr + γq ( A + 1) ≃ s A + 1) (cid:20) − ± (cid:18) nr + γqs (cid:19)(cid:21) . (23)As follows from Eq. (21), the boundary curve approaches x = 0 in the limit of large | s | as( A + 1) xs = 0 . (24)Then, the plus sign should be chosen in Eq. (23). Thus, the asymptotic boundary curve in the limit s ≫ xs = nr + γqA + 1 . (25)If the coefficients of this hyperbola is negative, i.e. nr + γq <
0, the boundary curve lies in the region of x > s <
0. In such cases, given that x >
0, the line Γ always touches the boundary curve under increasing the diffusionconstant D w , as shown in Figs. 3 and 4. Therefore, we require nr + γq < x > f w h u + g w h v < , (26) f u + g v > . (27) (a) ( x, s ; A, B ) = I wav s x x ( x, s ; A, B ) = I ST increasing w D (c) s ( x, s ; A, B ) = I wav x x ( x, s ; A, B ) = I ST increasing w D x − x + P Γ (b) s ( x, s ; A, B ) = I wav x x ( x, s ; A, B ) = I ST increasing w D x − x + P Γ P Γ FIG. 3: Boundary curves ˆ I wav = 0 (red) and ˆ I st = 0 (blue) when (a) the condition (30) and (b,c) the condition (34) aresatisfied. All species are mobile. Positions of the line Γ defined by (10) are shown for two different values of D w . The boundarycurve ˆ I st = 0 can behave as shown in (b) or (c), depending on the model parameters. The boundary curve for the Turing instability ˆ I st ( x, s ; A, B ) = 0 is given by s ( x ) = − ( A + 1) ( Anr + γq ) x − ( A + 1) ( βγn + αqr ) + ( A + 1) B ( nr − γq ) A ( A + 1) x − ( A − Bx − αβ ( A + 1) − B , (28)which is a single-valued function of x .Let us assume that all three reactants diffuse over the space, D u = 0 , D v = 0 and D w = 0, which leads to A = 0. Ifthe denominator on the right hand side of Eq. (28) is not zero for all values of x , then s ( x ) is a continuous function.Figure 3 (a) illustrates the qualitative shape of the boundary curves ˆ I wav = 0 and ˆ I st = 0 in such a situation. Inthis case, the line Γ touches the boundary curve ˆ I wav = 0 without intersecting ˆ I st = 0 as D w is increased. As aconsequence, the wave instability takes place. The conditions are given by( A − B + 4 A ( A + 1) (cid:2) αβ ( A + 1) + B (cid:3) < f u g v − f v g u > ( f u D v + g v D u ) D u D v . (30)If the denominator on the right hand side of Eq. (28) vanishes at x − and x + , the boundary curve ˆ I st = 0 divergesin two ways depending on the values of n, q, r, α, β and γ (Fig. 3 (b) and (c)). In both cases, if x ≤ x − , the line Γtouches ˆ I wav = 0 first. Thus, if the condition x ≤ x − = ( A − B − ( A + 1) p B + 4 αβA A ( A + 1) (31)which is equivalent to f u g v − f v g u > f u D v + g v D u ≤ D w is increased.Next, we assume that the reactant V does not diffuse, that is D v = 0, leading to A = 0. In this case, the boundarycurve ˆ I st = 0 is given by s ( x ) | A =0 = − qx + κ ( n − qr ) + B ( nr − q ) Bx + 1 − B , (34) ( x, s ; A, B ) = I wav s (b) x x ( x, s ; A, B ) = I ST p p s ( x, s ; A, B ) = I wav (a) x P x ( x, s ; A, B ) = I ST p p increasing w D increasing w D Γ P Γ FIG. 4: Boundary curves ˆ I wav = 0 (red) and ˆ I st = 0 (blue). The reactant V does not diffuse ( D v = 0). The line Γ definedby (10) is shown as black lines for two different values of D w . The boundary curve ˆ I st = 0 can behave as shown in (a) or (b)depending on the parameters n, q, r, α, β and γ . which diverges at x = B − /B = p − /p [Figs. 4 (a) and (b)]. In this case, if x ≤ p − /p , the line Γ touchesˆ I wav = 0 without intersecting ˆ I st = 0 when D w is increased. Thus, the instability condition for a system with twodiffusible reactants is given by x ≤ p − p , (35)which implies g v ( f u g v − f v g u ) ≤ . (36)On the other hand, if the reactant U does not diffuse, D u = 0, one can choose a plane m = A ′ p + B ′ where A ′ = D u /D v and B ′ = m − p D u /D v and derive the condition f u ( f u g v − f v g u ) ≤ . (37)Thus, sufficient conditions for the wave instability under increasing diffusion constant D w of the reactant W hasbeen derived for the four cases depending on the diffusion constants. The first two conditions (26) and (27) arecommon to all cases. The last condition depends of the diffusion constants, i.e. we have f u g v − f v g u > ( f u D v + g v D u ) D u D v ( D u = 0 , D v = 0) , (38) f u g v − f v g u > f u D v + g v D u ≤ D u = 0 , D v = 0) , (39) f u ( f u g v − f v g u ) < D u = 0) , (40) g v ( f u g v − f v g u ) < D v = 0) . (41)These different requirements can however be expressed in a single equation;det (cid:18) f u + D u Λ f v g u g v + D v Λ (cid:19) = 0 for any Λ < . (42)with at least two non-vanishing diffusion constants.The wave instability may take place also under the variation of other two diffusion constants D u or D v . Thecorresponding sufficient conditions for each case can be obtained by permutating three variables u, v and w .In summary, the wave instability occurs under increasing diffusion constants D u , D v or D w , if the conditions g v + h w > g u f v + h u f w < (cid:18) g v + D v Λ g w h v h w + D w Λ (cid:19) = 0 for any Λ < , (43) h w + f u > h v g w + f v g u < (cid:18) h w + D w Λ h u f w f u + D u Λ (cid:19) = 0 for any Λ < , (44)or f u + g v > f w h u + g w h v < (cid:18) f u + D u Λ f v g u g v + D v Λ (cid:19) = 0 for any Λ < . (45)are satisfied.We have derived sufficient conditions for the wave instability in general three-component reaction-diffusion systems.The conditions are formulated in terms of the Jacobian matrix elements at a steady state and of the diffusion constantsThey do not depend on model details. Once these conditions are satisfied, the wave instability occurs as we increasethe diffusion mobility of one of the reacting species.Our general results are applicable for systems of various origins, including biological, chemical, physical and ecolog-ical systems. Our analysis has revealed that the wave instability may occur even if one of three reactants is immobile.This result can be important in a variety of applications involving both diffusible and non-diffusible reactants.Authors acknowledge the financial support through the DFG SFB 910 program “Control of Self-Organizing Non-linear Systems” in Germany, through the Fellowship for Research Abroad, KAKENHI and the FIRST Aihara Project(JSPS), and the CREST Kokubu Project (JST) in Japan. [1] Turing, A. M. The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B , 37-72 (1952).[2] Knobloch, E. & De Luca, J. Amplitude equations for travelling wave convection.
Nonlinearity , 975 (1990).[3] Walgraef, D. Spatio-Temporal Pattern Formation, with Examples in Physics, Chemistry and Materials Science. (Springer,New York, 1997).[4] Sick, S., Reiniker, S., Timmer, J. & Schlake, T. WNT and DKK determine hair follicle spacing through a reaction-diffusionmechanism.
Science , 1447-1450 (2006).[5] Kondo, S. & Asai, A reaction-diffusion wave on the skin of the marine angelfish Pomacanthus.
Nature , 765-768 (1995).[6] Nakamasu, A., Takahashi, G., Kanbe, A. & Kondo, S. Interactions between zerafish pigment cells responsible for thegeneration of Turing patterns.
Proc. Natl. Acad. Sci. USA , 8429-8434 (2009).[7] Castets, V., Dulos, E., Boissonade, J. & De Kepper, P. Experimental evidence for a sustained standing Turing-typenonequilibrium chemical pattern.
Phys. Rev. Lett. , 2953-2956 (1990).[8] Ouyang, Q. & Swinney, H. L. Transition from a uniform state to hexagonal and striped Turing patterns. Nature ,610-612 (1991).[9] Yang, L., Dolnik, M., Zhabotinsky, A. M. & Epstein, I. R. Pattern formation arising from interactions between Turing andwave instabilities.
Journal of Chemical Physics , 7259 (2002).[10] Vanag, V. K., & Epstein, I. R. Pattern formation in a tunable medium: The Belousov-Zhabotinsky reaction in an aerosolOT microemulsion.
Phys. Rev. Lett. , 228301 (2001).[11] Wang, W., Liu, Q. X. & Jin, Z. Spatiotemporal complexity of a ratio-dependent predator-prey system. Phys. Rev. E75