We show that a model reaction-diffusion system with two species in a monostable regime and over a large region of parameter space, produces Turing patterns coexisting with a limit cycle which cannot be discerned from the linear analysis. As a consequence, Turing patterns oscillate in time, a phenomenon which is expected to occur only in a three morphogen system. When varying a single parameter, a series of bifurcations lead to period doubling, quasi-periodic and chaotic oscillations without modifying the underlying Turing pattern. A Ruelle-Takens-Newhouse route to chaos is identified. We also examined the Turing conditions for obtaining a diffusion driven instability and discovered that the patterns obtained are not necessarily stationary for certain values of the diffusion coefficients. All this results demonstrates the limitations of the linear analysis for reaction-diffusion systems

Aragón, J. L.

Baker, R. E.

Barrio, R. A.

Maini, P. K.

Woolley, T. E.

Physical Review E

English

Oxford University Research Archive (ORA)

Non-linear effects on Turing patterns: time oscillations and chaos.

We show that a model reaction-diffusion system with two species in a monostable regime and over a large region of parameter space produces Turing patterns coexisting with a limit cycle which cannot be discerned from the linear analysis. As a consequence, the patterns oscillate in time. When varying a single parameter, a series of bifurcations leads to period doubling, quasiperiodic, and chaotic oscillations without modifying the underlying Turing pattern. A Ruelle-Takens-Newhouse route to chaos is identified. We also examine the Turing conditions for obtaining a diffusion-driven instability and show that the patterns obtained are not necessarily stationary for certain values of the diffusion coefficients. These results demonstrate the limitations of the linear analysis for reaction-diffusion systems

Aragón, JL

Barrio, RA

Woolley, TE

Baker, RE

Maini, PK

Oxford University Research Archive

Nonlinear effects on Turing patterns: time oscillations and chaos.

We show that a model reaction-diffusion system with two species in a monostable regime and over a large region of parameter space produces Turing patterns coexisting with a limit cycle which cannot be discerned from the linear analysis. As a consequence, the patterns oscillate in time. When varying a single parameter, a series of bifurcations leads to period doubling, quasiperiodic, and chaotic oscillations without modifying the underlying Turing pattern. A Ruelle-Takens-Newhouse route to chaos is identified. We also examine the Turing conditions for obtaining a diffusion-driven instability and show that the patterns obtained are not necessarily stationary for certain values of the diffusion coefficients. These results demonstrate the limitations of the linear analysis for reaction-diffusion systems. © 2012 American Physical Society.This work was supported by CONACyT and DGAPA-UNAM, Mexico, under Grants No. 79641 and No. IN100310-3, respectively, and was based on work supported in part by Award No. KUK-C1-013-04, made by King Abdullah University of Science and Technology (KAUST)

KAUST Repository (King Abdullah University of Science and Technology)

Nonlinear effects on Turing patterns: Time oscillations and chaos

Nonlinear effects on Turing patterns: time oscillations and chaos

ORA - Oxford University Research Archive

PHYSICAL REVIEW E 86, 026201 (2012)Nonlinear effects on Turing patterns: Time oscillations and chaosJ. L. Aragón,1 R. A. Barrio,2 T. E. Woolley,3,4 R. E. Baker,4 and P. K. Maini4,51Departamento de Nanotecnologı́a, Centro de Fı́sica Aplicada y Tecnologı́a Avanzada, Universidad Nacional Autónomade México, Apartado Postal 1-1010, Querétaro 76000, México2Departamento de Fı́sica Quı́mica, Instituto de Fı́sica, Universidad Nacional Autónoma de México, Apartado Postal 01000,76230 Distrito Federal, México3Oxford Centre for Collaborative Applied Mathematics, Mathematical Institute, University of Oxford, 24-29 St. Giles’,Oxford OX1 3LB, United Kingdom4Centre for Mathematical Biology, Mathematical Institute, University of Oxford, 24-29 St. Giles’,Oxford OX1 3LB, United Kingdom,5Oxford Centre for Integrative Systems Biology, Department of Biochemistry, University of Oxford,South Parks Road OX1 3QU, United Kingdom(Received 12 April 2012; published 8 August 2012)We show that a model reaction-diffusion system with two species in a monostable regime and over a largeregion of parameter space produces Turing patterns coexisting with a limit cycle which cannot be discerned fromthe linear analysis. As a consequence, the patterns oscillate in time. When varying a single parameter, a series ofbifurcations leads to period doubling, quasiperiodic, and chaotic oscillations without modifying the underlyingTuring pattern. A Ruelle-Takens-Newhouse route to chaos is identified. We also examine the Turing conditionsfor obtaining a diffusion-driven instability and show that the patterns obtained are not necessarily stationary forcertain values of the diffusion coefficients. These results demonstrate the limitations of the linear analysis forreaction-diffusion systems.DOI: 10.1103/PhysRevE.86.026201 PACS number(s): 05.45.Pq, 82.40.Bj, 87.18.HfI. INTRODUCTIONReaction-diffusion systems can produce spatially periodicstationary patterns, usually known as Turing patterns. In hisseminal work [1], Turing suggested the possibility of obtainingperiodic patterns, oscillating both in time and space, and thatthey only occur with three or more morphogens. More recently,it has been shown that Turing patterns oscillating in timecan be generated in a two-species system when Turing andHopf instabilities interact [2], by coupling stationary Turingstructures with oscillating layers [3], or in some bistablesystems [4] in which neither Hopf nor wave instabilities arerequired.Diffusion-induced chaos has been investigated due to itsconnection with chemical chaos [5] and to the existenceof some biological phenomena displaying this behavior [6].Spatiotemporal chaos in reaction-diffusion systems can arisein different ways, for instance from the interaction betweenTuring and Hopf bifurcations [7], by homoclinic explosions[8], by interaction with propagating fronts [9], or behindtransition fronts simulated by moving boundaries [10].In this paper we study a simple reaction-diffusion systemwith two species in a monostable regime and with a fixeddomain size, which produces temporal oscillations and chaosin a region where linear analysis predicts only Turing patternswith no temporal oscillations. As a single parameter varies, theoscillating patterns undergo successive bifurcations, produc-ing period doubling, tori, and chaos. Furthermore, when theratio of diffusion coefficients fulfills the conditions for a Turinginstability, we find the unexpected result that in a large regionof parameter space, only Turing patterns oscillating in timeemerge from the homogeneous steady state. The importance ofthese results is that they not only show the stringent limitationsof linear analysis but also open up new perspectives on theapplication of reaction-diffusion systems to real biological andchemical systems.II. THE MODELThe reaction-diffusion model that we study in this workis the so-called Barrio-Varea-Aragon-Maini (BVAM) model[11]:∂u∂t= D∇2u + η(u + av − Cuv − uv2),(1)∂v∂t= ∇2v + η(bv + Hu + Cuv + uv2),with zero flux boundary conditions. This model presents arichness of behavior, making it suitable as a good “laboratory”to gain insight into the mechanisms of pattern formation [12].We will keep parameters D, C, and H free while the remainingparameter values will be η = 1, a = −1, and b = −3/2.According to linear analysis, this choice ensures that theequilibrium point (0,0), in the absence of diffusion, is a stablespiral for the values of H considered in this paper and alsothat the trace of the Jacobian is negative as we shall see below.A. Linear analysisThe diffusionless system has the following equilibriumpoints:(u0,v0) = (0,0) ,(u∓,v∓) =(−5C ∓ √5√4H + 4 ,−5C ∓ √5√10),where  = 5C2 − 8H + 12. We observe that if  > 0, thenthree real equilibrium points are present, but if  < 0, then026201-11539-3755/2012/86(2)/026201(4) ©2012 American Physical SocietyARAGÓN, BARRIO, WOOLLEY, BAKER, AND MAINI PHYSICAL REVIEW E 86, 026201 (2012)2.5 2.3 2.1 1.9 1.7 1.5 1.3 1.1 0.9 0.7 0.5 0.3 0.11.01.52.02.53.03.5CHFIG. 1. (Color online) The bold line is the saddle-node bifurcationcurve (2) of system (1) without diffusion. (Only the region C < 0 isshown.) In the shaded region, (0,0) is the only equilibrium point,which is linearly stable. The dashed curve is the locus of a Hopf bi-furcation calculated numerically from Eq. (1) as described in Sec. III.(0,0) is the only real equilibrium point, hereafter denoted as O.The particular case  = 0 defines a saddle-node bifurcationcurve given by the equation (see Fig. 1)H SN = 18 (5C2 + 12). (2)Along this curve, two equilibrium points are present, and tothe left of H SN, the three equilibrium points are recovered. Inthe following we consider parameter values C and H insidethe shaded region in Fig. 1. The Jacobian evaluated at O isJO =(1 −1H − 32).From this we obtain the eigenvalues λ∓ = 14 (−1 ∓√25 − 16H ) from which it is concluded that O is linearlystable provided H > 3/2. The value H = 3/2 corresponds tothe minimum of H SN. Thus, the shaded region bounded byH SN in Fig. 1 corresponds to a monostable regime.Since tr (J 0) = −1/2 is always negative, a Hopf bifurca-tion, which requires tr (J 0) = 0, is not feasible.B. The Turing instabilityBy considering diffusion in the linearized system, followingthe standard analysis for a Turing instability [13], we constructthe matrix  = J 0 − k2D, whereD = diag (D,1), from whichwe find tr() = −1/2 − (1 + D) k2. Since D  0, it followsthat tr() < 0; thus, neither oscillations with k = 0 (wavebifurcation) are allowed, and consequently, linear analysispredicts no oscillations at all. Since O is linearly stable,there remains the possibility of having a Turing (diffusion-driven) instability. Using D as the bifurcation parameter, theconditions for a Turing instability [13] yield the critical valuesk2c = H −32+√H(H − 32), (3)Dc = 89[H − 34−√H(H − 32)]. (4)Notice that both critical quantities are independent of C.III. RESULTSThe model in Eq. (1) was solved numerically in onedimension with zero flux boundary conditions and randominitial conditions around the equilibrium point (0,0). Thesecond partial derivative terms in Eq. (1) are discretized usingfinite differences, including zero flux boundary conditions inthe discretization. Using a mesh with N nodes yields a set of2N ordinary differential equations which were solved usingthe software library CVODE [14] with a numerically calculatedJacobian matrix. Scalar absolute and relative tolerances of10−9 were used, and t = 0.01. We set N = 50 and x =0.2; thus, the system size is L = 9.8.Since the critical parameters depend only on H , we firstobtain a Dc valid for the shaded region in Fig. 1, thatis, the range 1.5 < H < 3.8. In this interval Dc decreasesmonotonically; thus, we need to consider only the values ofDc when H is evaluated at the boundaries of the inequality.From Eq. (4) we see that for H = 1.5, Dc ≈ 0.66 and thatfor H = 3.8, Dc ≈ 0.083. Thus, if we set D = 0.08, Turingpatterns should be produced in the range 1.5 < H < 3.8. Wenow consider C as a control parameter to explore its effect onthe generated patterns. In what follows and unless otherwiseindicated, we set D = 0.08 and H = 3. For this value of H ,Eq. (3) gives kc ≈ 1.9029, and from the dispersion relation, itis found that the admissible (taking into account boundaryconditions), most rapidly growing mode is k = 7π/9.8 ≈2.244. Therefore, we expect a wavelength of 2π/2.244 ≈ 2.8or 9.8/2.8 ≈ 3.5 stripes in the simulations.In Fig. 2 numerical simulations for three different values ofC are displayed; in the left column, the space-time plots foru up to 15 000 time iterations only (t = 150) are shown forclarity. In the right column, the corresponding time series ofthe central point u (25) are shown. As expected, for C = −0.6,as in Fig. 2(a), a stationary Turing pattern is obtained. WhenC = −1 [Fig. 2(b)], a stationary pattern is generated initially,but eventually regular oscillations emerge. If C is decreasedfurther, approaching the line H SN, the pattern oscillates in anapparently chaotic fashion as shown in Fig. 2(c) for C = −1.5.Since oscillations are not predicted by the linear analysisand the conditions for a Turing instability are valid for theseparameter values (that is, there are admissible modes growing),it seems plausible to assume that these oscillations are aconsequence of the nonlinear terms and, in particular, of therelative strength of the quadratic and cubic nonlinearities,given by the magnitude of C. Since an analytical treatmentof nonlinearities is often prohibitive, we used the bifurcationsoftware AUTO [15] to detect possible oscillatory instabilities.In Fig. 3(a), we show the bifurcation diagram. As Cdecreases from −0.1, the stable stationary (Turing) stateloses stability through a supercritical Hopf bifurcation atC = −0.7546515 (successive bifurcations will be discussedafterwards). Bifurcation diagrams can be generated for differ-ent values of H so that a curve of Hopf bifurcations in a C vsH plane can be calculated numerically. By using steps of 0.01in H , the resulting curve is plotted in Fig. 1 as a dashed line.This curve separates the monostability region where conditionsfor a Turing instability are satisfied into two regions: on theright-hand side of the Hopf curve, Turing patterns are realized,and on the left, patterns oscillate in time.026201-2NONLINEAR EFFECTS ON TURING PATTERNS: TIME . . . PHYSICAL REVIEW E 86, 026201 (2012)−1.0−0.50.00.51.01.50 50 100 150−1.0−0.50.00.5Timeu(25)−0.50.00.51.0−2−101230 50 100 150−2.0−1.5−1.0−0.50.00.51.0Time0 50 100 150−0.6−0.4−0.20.00.2Timeu(25)(a)(b)(c)0150t0150t0150tx0.0 9.8FIG. 2. (Color online) Numerical simulations of the system (1) inone dimension for H = 3 and D = 0.08. The left column shows thespace-time plot of u for three values of C. (The plots for v are similar.)On the right, the time series registered at the central point u (25)of the corresponding space-time plots are shown. (a) C = −0.6,(b) C = −1, and (c) C = −1.5.It is worth mentioning that in the interval 2.1 < H < 2.5,no Hopf bifurcation occurs, and in this region only oscillatoryTuring patterns emerge from the homogeneous steady stateinstead of the expected stationary patterns. We observe fromEq. (4) that Dc increases as H diminishes such that when,for instance, H = 1.8, the critical value Dc = 0.2801  0.08.According to linear analysis, the uniform steady state becomesunstable to many modes, and one expects that the most rapidly-0.50.00.51.01.52.02.53.03.5u(25)-2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2CHBPD TR (a) (b)-2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2C-0.10-0.050.000.050.100.150.200.25λFIG. 3. (a) Numerically calculated bifurcation diagram of system(1) using AUTO. Supercritical Hopf (HB), period doubling (PD), andtorus (TR) bifurcations are detected at C = −0.7546515, −1.200656,and −1.281831, respectively. (b) Numerically calculated maximalLiapunov exponents.−1.0 −0.5 0.0 0.5−1.5−0.50.00.51.0u(25)v(25)−1.0 −0.5 0.0 0.5u(25)v(25)u(25)v(25)1.01.0u(25)v(25)−1.5−0.50.00.51.0−1.5−0.50.00.51.0−1.0 −0.5 0.0 0.5−1.50.50.0−1.0−0.51.5−1.0 −0.5 0.0 0.5−1.5(a)(c)(b)(d)FIG. 4. Phase portraits of the central points for different valuesof C after 50 000 time iterations (t = 500). (a) C = −1.15, (b) C =−1.21, (c) C = −1.285, and (d) C = −1.5.growing mode will eventually dominate. However, numericalsimulations show that this is not so. If D is just below 0.2801,stationary Turing patterns are obtained as expected, but whenD  0.2801, only Turing patterns oscillating in time emergefrom the homogeneous steady state.A. ChaosOnce the Hopf bifurcation line is crossed, if C decreasesfurther, the oscillations become irregular, and near the H SNline they appear to become chaotic. In Fig. 3(b), maximalLyapunov exponents (λ) are plotted versus C. According toour numerical results, λ becomes zero at C = −1.285 andpositive after this value except for a small region of stability(−1.33  C  −1.35). Consequently, chaos is expected when−1.285 < C < −1.33 and when C < −1.35.In Fig. 4, phase portraits are displayed for different valuesof C. After crossing the Hopf bifurcation line, a limit cycleis formed, as shown in Fig. 4(a), for C = −1.15. WhenC = −1.21, the cycle doubles its period as shown in Fig. 4(b).At smaller values of C a torus is generated at C = −1.285as in Fig. 4(c), and chaotic orbits, shown in Fig. 4(d),appear for C = −1.5. All these results are consistent with thesuccessive bifurcations displayed in the bifurcation diagramof Fig. 3(a): after the Hopf bifurcation (at C = −0.7546515),the limit cycle loses stability through a period doubling forC = −1.200656, and this cycle itself loses stability through atorus at C = −1.281831.The quasiperiodicity of the torus in Fig. 4(c) can be verifiedby calculating its power spectrum. In Fig. 5(a), the powerspectrum for the torus (C = −1.285) is shown. The limitcycle after the Hopf bifurcation has frequency ω1 = 0.3382,the period doubling bifurcation introduces the frequencyω1/2, and the torus bifurcation introduces the new frequencyω2 = 0.2732. The ratio of frequencies ω1/ω2 = 1.2379 and,as it is impossible to discern numerically in this way whetheror not a number is irrational, we cannot distinguish betweenquasiperiodic orbits or orbits with very long periods; thus, weverified that the Poincare plots are typical of quasiperiodicorbits (ringlike). The power spectrum for C = −1.455, shown026201-3ARAGÓN, BARRIO, WOOLLEY, BAKER, AND MAINI PHYSICAL REVIEW E 86, 026201 (2012)0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5Frequency0.000.050.100.150.200.25Amplitude (arb. units)ω1ω1/22ω1ω23ω1/23ω15ω1/20.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4Frequency0.000.050.100.150.20Amplitude (arb. units)FIG. 5. (Color online) Power spectrum for C = −1.285 after500 000 iterations (t = 5000). The frequency ω1 = 0.3382 and itsinteger multiples appear after the Hopf bifurcation. The perioddoubling bifurcation introduces the frequency ω1/2 and the torusbifurcation ω2 = 0.2732. Inset: Power spectrum for C = −1.455showing a continuous band of frequencies, typical of a chaoticoscillation.in the inset of Fig. 5, exhibits a continuous band of frequencies,confirming the chaotic nature of the oscillation.From all of these numerical results, we observe the follow-ing sequence of bifurcations: stationary → Hopf → two-torus.Numerically, when a third frequency is about to appear(three-torus), small perturbations transform the quasiperiodicorbit into a chaotic orbit, or a strange attractor. This resultis consistent with the well-known Ruelle-Takens-Newhouseroute to chaos [16].B. Two-dimensional oscillating patternsIn two dimensions, the parameter C controls the typeof Turing pattern, spots or stripes [11], so by varying thisparameter, we can generate breathing spots or stripes in aTABLE I. Parameter values used in supplementary movies.Movie Name D H CD008H3C-1.mp4 0.08 3.0 1.0D008H3C-15.mp4 0.08 3.0 1.5D008H16C-02.mp4 0.08 1.6 0.2D008H35C-12.mp4 0.08 3.5 1.2regular, quasiperiodic, or chaotic fashion. Some illustrativemovies of the numerical simulations in two dimensionsare included as supplementary material [17]. Equation (1)was implemented in the commercial finite element packageCOMSOL MULTIPHYSICS with zero flux boundary conditionsand random initial conditions. The parameter values for thesupplementary movies are shown in Table I.IV. CONCLUSIONSIn summary, we have studied a simple reaction-diffusionsystem that exhibits oscillations and chaos in a novel way.Turing patterns lose stability through a Hopf bifurcation whichcannot be discerned from present linear analysis techniques.The bifurcation parameter C is absent in the linear analysis,and as it varies, quasiperiodic and chaotic oscillations aregenerated in a large region of parameter space. A Ruelle-Takens-Newhouse route to chaos is identified. The mainconclusion of this work is that linear analysis may be a stringentlimitation for studying reaction-diffusion systems and that thenonlinearities play an important role, not only on stabilizing apattern, but also in producing unsuspected bifurcations.ACKNOWLEDGMENTSThis work was supported by CONACyT and DGAPA-UNAM, México, under Grants No. 79641 and No. IN100310-3, respectively, and was based on work supported in partby Award No. KUK-C1-013-04, made by King AbdullahUniversity of Science and Technology (KAUST).[1] A. M. Turing, Philos. Trans. R. Soc., B 237, 37 (1952).[2] L. Yang, M. Dolnik, A. M. Zhabotinsky, and I. R. Epstein, Phys.Rev. Lett. 88, 208303 (2002).[3] L. Yang and I. R. Epstein, Phys. Rev. Lett. 90, 178303 (2003).[4] V. K. Vanag and I. R. Epstein, Phys. Rev. E 71, 066212(2005).[5] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence(Springer, Tokyo, 1984).[6] L. F. Olsen and H. Degn, Q. Rev. Biophys. 18, 165 (1985).[7] A. De Wit, G. Dewel, and P. Borckmans, Phys. Rev. E 48, R4191(1993); V. Petrov, S. Metens, P. Borckmans, G. Dewel, andK. Showalter, Phys. Rev. Lett. 75, 2895 (1995).[8] J. Elezgaray and A. Arneodo, Phys. Rev. Lett. 68, 714 (1992).[9] J. H. Merkin, V. Petrov, S. K. Scott, and K. Showalter, Phys.Rev. Lett. 76, 546 (1996).[10] J. A. Sherratt, Dynam. Stabil. Syst. 11, 303 (1996).[11] R. A. Barrio, C. Varea, J. L. Aragón, and P. K. Maini, Bull.Math. Biol. 61, 483 (1999); T. Leppänen, Ph.D. thesis, HelsinkiUniversity of Technology, 2004 (unpublished).[12] C. Varea, J. L. Aragón, and R. A. Barrio, Phys. Rev. E 60, 4588(1999); D. A. Striegel and M. K. Hurdal, PLoS Comput. Biol.5, e1000524 (2009); R. A. Barrio, R. E. Baker, B. Vaughan,K. Tribuzy, M. R. de Carvalho, R. Bassanezi, and P. K. Maini,Phys. Rev. E 79, 031908 (2009); T. E. Woolley, R. E. Baker,P. K. Maini, J. L. Aragón, and R. A. Barrio, ibid. 82, 051929(2010).[13] J. D. Murray, Mathematical Biology (Springer, Berlin, 2002),Vol. 2.[14] S. D. Cohen and A. C. Hindmarsh, Comput. Phys. 10, 138(1996).[15] B. Ermentrout, Simulating, Analyzing, and Animating Dynami-cal Systems: A Guide to XPPAUT for Researchers and Students(Society for Industrial Mathematics, Philadelphia, 2002).[16] S. Newhouse, D. Ruelle, and F. Takens, Commun. Math. Phys.64, 35 (1978); J. P. Eckmann, Rev. Mod. Phys. 53, 643 (1981).[17] See movies in the Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevE.86.026201 for numerical sim-ulations of two-dimensional oscillating patterns.026201-4

Mathematical Institute Eprints Archive

Non-linear effects on Turing patterns: time oscillations and chaos.

Abstract

Similar works

Full text

Available Versions

Oxford University Research Archive (ORA)

Oxford University Research Archive

KAUST Repository (King Abdullah University of Science and Technology)

Oxford University Research Archive (ORA)

ORA - Oxford University Research Archive

Oxford University Research Archive

Mathematical Institute Eprints Archive

Oxford University Research Archive

ORA - Oxford University Research Archive