P ELSEVIER
CA )
Physica D 99 (1996) 303-317
Turing pattern formation in heterogeneous media J.-E Voroney a., A.T. Lawniczak a, R. Kapral b a Department of Mathematics and Statistics and Guelph-Waterloo Programfor Graduate Work in Physics, Universily of Guelph, Guelph, Ont., Canada NIG 2W! b Chemical Physics Theory.Group, Department of Chemistry, Universi~ of Toronto, Toronto, Ont., Canada M5S IAi
Received 24 October 1995; revised 20 March 1996; accepted 9 April 1996 Communicated by J.D. Meiss
Abstract The effects of spatial inhomogeneities in the initial distribution of immobile complexing species on the spatiotemporal dynamics of activator-inhibitor systems are investigated. If the initial concentration of the complexing species is below a critical value and homogeneously distributed, the stable attractor is a limit cycle; for higher concentrations the stable attractor is a Turing pattern. When the complexing agent is inhomogeneously distributed interactions between the oscillatory dynamics and Turing pattern formation are possible. Such phenomena are investigated in a model system, the Sel'kov model with a complexing reaction. The effects of the symmetries and lengths that characterize the initial distribution of complexing agent on the spatiotemporal attracting states are studied. The study provides insight into the types of behavior that are possible in
heterogeneous media undergoing Turing-like bifurcations. PACS': 82.20.Mj; 82.20.wt: 82.65.Jv
Keywords: Turing systems; Spatial pattern formations; Heterogeneous media; Chemical waves
1. Introduction The nature of the medium in which reaction and diffusion processes occur can influence the character of the spatiotemporal dynamics that is observed. For example, the rates at which reactions occur can depend on temperature or light intensity which can vary throughout a medium; furthermore diffusion can be affected by the make-up of the medium and its geometry. Spatial dependence of diffusion and rate coefficients can have interesting consequences in far-from-equilibrium systems where a variety of chemical patterns are observed as a result of
symmetry-breaking bifurcations. In this context, wave propagation in heterogeneous excitable and oscillatory media has been investigated (see, for instance [ 1]), reaction-diffusion fronts in periodically layered media have been studied [2], and diffusion-driven instabilities ha"e been explored when the diffusion coefficients are spatially dependent [3]. In this article we examine how an inhomogeneous distribution of an immobile reactive species can affect Turing pattern [4] formation. This is an especially interesting case to consider since the symmetries and length scales associated with the heterogeneous nature of the system can interact with the characteristic symmetries and lengths of the Turing pattern to produce new structures.
* Corresponding author. 0167-2789/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved Pll SO ! 67-2789(96)00132-7
304
J..P. Voroney et aL/Physica D 99 (1996) 303-317
Recent experiments on the chlorite-iodide-malonic acid: (CIMA) reaction have been carded out in gel media with an immobilized complexing species [5,6]. A natural source of system inhomogeneity arises from the m~hanism which involves the binding of the mobile activator species to an immobile complexing agent. This mechanism is presumed to modify the reaction and diffusion processes [7,8] that lead to Turing bifurcations in systems with diffusion coefficients that would not otherwise give rise to such bifurcations, i In addition, the immobilized activator is no longer free to react as it would if it were not immobilized. Thus both transport and chemical kinetics of the activator species are affected by the presence of the ;.mmobile species bound to the medium. If the complexing species is inhomogeneously distributed with various symmetries and characteristic lengths, these inhomogeneities can affect the subsequent Turing and "luring-like pattern formation processes. This is the problem we address here. "The paper is organized as follows: an outline of the linear stability analysis for an inhibitor-activator system, the Sei'kov model, supplemented with a complexing reaction, is presented in Section 2. A more general version of this analysis is presented in Appendix A. This section, while devoted to a standard analysis of a system with no underlying heterogeneity, serves to provide background for the discussion of heterogeneous media that is presented in the following section. Section 3 describes the effects of spatial inhomogeneities in the initial distribution of the immobile complexing species on the one- and two-dimensional dynamics of the Sel'kov model with a complexing reaction. Various spatially inhomogeneous initial distributions of complexing species are considered and their effects on Turing and Tunng-like patterns are investigated. The results show that a wide variety of complex attracting states I Actually, in the general case, one should carry out a bifurcation analysis of the reaction-diffusion system with the complexing reaction included, rather than simply adiabatically eliminating this step from the reaction mechanism to obtain renormalized reaction and diffusion terms. See, for example, [9] and Appendix A of this paper.
are possible as a result of the interactions between the symmetries and characteristic lengths of the initial distribution of complexing agent and those of the underlying reaction-diffusion system. The conclusions of the study as well as comments on the prospects for experimental observations of the phenomena are presented in Section 4.
2. Sel'kov model with complexing reaction We consider reaction-diffusion systems with two mobile chemical species X and Y and an immobile complexing species C whose only reaction is a reversible binding with the Y species to form an immobile product P. Let the concentrations of the reacting chemical species X, Y, C, and P be Px, Pr, Pc, and PP, respectively. The types of complexing reactions that we allow have the property Opc/at = -~pp/Ot so that the number of concentration variables can be reduced from four to three through the introduction of a (possibly) spatially dependent parameter r/: r/(r) = pc(r, t) + pp(r, t). It is straightforward to carry out a linear stability analysis for a general model with these features and such an analysis is sketched in Appendix A. For the present purposes we outline the results for a modification of the Sel'kov model [10] in order to investigate effects on the spatiotemporal dynamics due to heterogeneity in the distribution of the complexing species. The reaction scheme considered here is
A k~1_~__.X,
X + 2Y k--~2~3y,
k-t
"v+c k-3
k-2 k4
(1)
e.
k-4
Under the assumptions that the concentrations of the A and B .~I~i~ : ~ ; ~i,,ed by external constraints, the variables are scaled as in [ l l ], where time is in units of k3, concentrations are in units of ~/rk2-/k3. With the parameters defined as a = (kl/k3)PA, b = (k-3/k3)pB, y = (k-i/k3), ( = ( k 4 / ~ ) , p = (k-4/k3), and k2 = k-2, the reaction-diffusion equations for the Sel'kov model (1) are of the form
J.-P. Voroney et aL/Physica D 99 (1996) 303-317
Op.__.x.x = f ( P x , Or) + DV2px, Ot Opr = g ( P x , PY) + h ( p y , Pc, rl) at
(2)
+ DTZpy,
Opc = h ( p r , Pc, rl), Ot
where f , g and h are functions describing the local chemical kinetics and are given by f ( P x , Py) = a - FOx - p x p 2 + p3, g(Px, Pv) = b - pv + p x p 2 -- 9 3. h(py, Pc, rl) = - ( P Y P c + vpp
(3)
= vrI - ( p Y p C - v p C .
The diffusion coefficients Dx and Dr of the X and Y species, respectively, have been taken equal, Dx = Dy=D.
The spatially homogeneous, steady-state concentration of the species Y, p~,, can be found from the solution of P~,3( 1 + F) -- P~,2(a + b) + YP~' - bF = 0,
(4)
while p~(, p~ and p~, can be obtained from I'P~ = a + b - p~, vq
pb = v + ¢ p ; ' p; =.
(5)
-
Descartes rule of signs [12] implies that (4) has one o~"three positive solutions and no negative solutions, so the system has one or three fixed points. We investigate changes in behavior caused by heterogeneity of the reaction medium, not bifurcations due to varying rate parameters; therefore, we vary the initial distribution of the immobile species, and in all the simulations described in the next sections we use the same non-dimensional parameters: a=0.51,
b=0.06,
¢=0.15,
v=0.1,
305
to zero, so r/(r) = pc(r, t = 0). Thus spatial inhomogeneity arises in the system (2) through .'h~ conservation of the total amount of immobile species, and any inhomogeneity present due to the initial distribution of complexing species will remain for all time. The unique fixed point of the system of equations (4) and (5) with parameters (6) is given by ,o]:=3.7,
p~, = 0.2,
p ~ __ i0
(7)
The linear stability analysis about this fixed point is straightforward if 1/does not depend on r. One may verify that the steady state is an unstable node when r / = 0. As r/increases it becomes stable when a complex conjugate pair of eigenvalues crosses the imaginary axis. The steady state is unstable to spatially uniform perturbations when 1/< 1.797, and stable to spatially uniform perturbations when r/ > 1.797, but unstable to perturbations with wave number 0.02 _< k 2 = ( m / L ) 2 D < 0.2 where L is the system size and L / m is a positive integer. The stable solution is a limit cycle when rl < 1.797 and a Turing pattern when r/> 1.797. Graphs of the real and imaginary parts of the eigenvalues of the matrix of the system iinearized about its fixed point as a function of r/are given in Fig. I. Using the Poincar6-Andronov-Hopf bifurcation theorem (see for e.g. [13]), it can be shown that a periodic orbit exists only for r / < 1.797. With this background we are now in a posmon to study heterogeneous media. Most of our results will rely on numerical simulations of the reactiondiffusion equations since analytical studies are difficult for arbitrary r/(r). All the simulation results presented below were obtained for systems with periodic boundary conditions using Pvt explicit Euler integration scheme with time step At and space step Ax given for each simulation in the figure captions. The time and space variables are reported in these units.
y=0.1, (6)
except for the parameter 1/which depends on the distribution of the immobile species. The spatial average of r/was 5.0 in all simulations. The initial concentration of the immobile product species P was always set
3. Spatiotemporal dynamics in heterogeneous media The nature of the observed spatiotemporal state of the system depends on the character of the
J.-P. Vomney et al.IPhysica D 99 (1996) 303-317
306
0.2
i
I
I
I
II
.
0.1 0
Re(.~)
.
-o.l -0.2
-0.3 -0.4 -0.5 0
I
I
l
!
I
1
2
3
4
5
!
!
6
77 0.15 0.1 0.05
Im(,~)
!
i,i
!
!
0
-0.05 -0.1 -0.15 0
l
I
I
a
|
]
2
3
4
5
r/
Fig. !. The top graph shows the real parts of the eigenvalues as a function of O. The bottom graph shows the imaginary parts of the complex conjugate eigenvalues as a function of O. Note the third eigenvalue -0.13 at 0 - 0. The one-dimensional stable manifold associated with this eigenvalue does not exist in the Sel'kov model without the complexing reaction. However. since 0 -- 0 implies that p¢(r. t) = 0 and pp(r. t) - 0 the disappearance of the stable manifold in the two-species Sel'kov system is not a physically unrealistic discontinuity in the dynamical behavior because there exists only one physical initial condition for ~1= 0 on the stable manifold.
inhomogeneous distribution of complexing agent. The initial spatial distributions may be regular or random and for each of these categories there are various ways in which the concentrations of complexing agent may vary. Naturally, we shall examine only few examples to illustrate the types of phenomena and focus on particular classes of regular and random distributions. Before beginning this discussion we remark that in most experiments carded out in gel reactors the distribution of complexing agent (e.g. starch in the CIMA reaction) is presumably random on a scale that is small compared to the characteristic lengths that enter the reaction-diffusion description. Diffusion will be able to homogenize the system on these small microscopic or mesoscopic length scales and the distribution of complexing agent will appear to be effectively homogeneous on macroscopic scales. In this instance one
expects to observe standard Turing pattern formation in the appropriate r/regime. This is demonstrated in Fig. 2 which shows Turing pattern formation for an initial, fine-scale, random distribution of species C. Consequently, the remainder of this article will focus on distributions whose characteristic length scales are comparable to those which are intrinsic to the reaction--diffusion system. System symmetries play an important role in determining the nature of the attracting state or states. 2 If 17 is spatially dependent then there exists a steady state that is spatially homogeneo~Js in X and Y, but spatially inhomogeneous in C and P. The functional form of this solution is given by (4) and (5), with 17 replaced by ~(r). Consider a two-dimensional system with orthogonal directions denoted by rl and r2. Because the kinetics and dynamics are local and the rate coefficients and diffusion coefficient are not spatially dependent then, for the complexing species initially distributed homogeneously, the symmetry group of the reaction--diffusion equations is E2, the Euclidean group of all translations, rotations, and flips in twodimensional space. The symmetry of the boundary value problem is thus a subgroup of E2. Due to the conservation law' r/(r) = pc(r, t) + pp(r, t), initial conditions of the immobile species influence bifurcations in the Serkov system. Thus the symmetry of the problem is determined by the function ~7(r) and the boundary conditions. Since the boundary conditions are periodic in both the rl and r2 directions, the symmetry of the reaction-diffusion equations is a subgroup of 0(2) x 0(2) C E2, the group of translations and flips in the rl and.r2 directions. Results of Golubitsky et al. [ 14], indicate that there exist two solutions that bifurcate from the spatially homogeneous solution, one of which has symmetry group of at least Zz(rl), the other of which has symmetry group of at least Zz(K2), where Xl and ~2 are flips in the rl and r2 directions. When the reaction-diffusion boundary value problem has no symmetry, steady state solutions II
2 Even if the spatial distribution of C occurs on scales small compared to those intrinsic to the system the symmetry of the distribution can influence the transient evolution of the system to its attracting state or reduce the number of attracting states.
J.-P. Voroney et a L / P h y s i c a D 99 (1996) 303-317
X
C -.~'.' .~,-~-~:: ,,~s~.r,.'..,~':~.".
•. ~ . ~ .
~.
:~-~.-...-." .:.~,.>~.%
,-:,-.
•. > . . . . . , . . ? ~ . . . ~ :
;.>.-~.~.-~
~.~:.
• ~..
..~..~.:..,¢.~.~,+..,., .~,,.
. .
Y
~'.%:
• ...
..~.~..,~,
307
.~,,.
,.~ ,.-~.,~.~,
~,~:.:.~:i~~:;~.,~ .~:~¢~!~~:~.~:.,.,~zS~':,~, Fig. 2. Snapshots of the concentrations of C, X and Y after 200 000 time steps. The initial distribution of C was taken to be binomial on the 128 x 128 spatial discretization grid with values ranging between zero and six and average equal to five. The time step and the scaled diffusion coefficient are At = 0.01 and D A t / A x 2 = 0.0156, respectively.
will have no global symmetry; however, solutions can in some approximation have a symmetry, or can have a local symmetry with defects that are visible on larger scales. In the remainder of this section we consider the complexing agent to be confined to n-dimensional disks or stripes and suppose that the initial concentration of C is constant throughout each disk or stripe. System parameters may be chosen such that Turing pattern formation would be possible in a large system for the t/(r) values within the disks or stripes, but in regions without the complexing species, i.e. when ~/(r) = O, the dynamics would be oscillatory in a large system. Given this situation we may then examine how the spatiotemporal dynamics depend on the various characteristic lengths in the system: ~-T, the Turing wavelength; s, the disk diameter or width of the stripe; I, the characteristic length of a gap between disks or stripes; and L, the system size. The distributions of disks may be regular or random, and the characteristics and symmetry of the distribution play a role in the pattern formation process as discussed above. The following sections present a variety of the spatiotemporal attracting states for different inhomogeneous distributions of C.
in space with period ~.d = S -q- I. The nature of the steady state will depend on s and l as well as ~.d in relation to the characteristic Turing length, kT. Note that the Turing length kT is independent of 17provided that 17 is everywhere above the critical value, 1.797, where a bifurcation from a limit cycle to a Turing pattern occurs; therefore, ~.T is taken to be that of a Turing pattern which forms in a large system with 17 = const. > 1.797.3 Consider the evolution when ~.d = ~.T = 30. We observe that at low values of s, s < 13, when there are large gaps between regions with complexing species, most realizations of the evolution starting from random initial conditions for X and Y concentrations contain traveling waves. As s increases, standing waves also begin to appear, and if s increases further timeindependent Turing patterns emerge. After s reaches a high enough value, s = 17, no oscillatory behavior is observed, and all realizations evolve to a Turingpattern state. As an illustration of the complicated nature of the dynamics for values of s intermediate between these extremes we present in Fig. 3 several examples of the system dynamics for s = 15 and I = 15 for different initial values of the X and Y concentrations. Sometimes a Turing pattern forms as in Fig. 3(a). Figs. 3(b) and (c) show that the same initial
3.1. Regular distributions: one-dimensional systems We begin with a discussion of one-dimensional systems where the one-dimensional disks (line segments) containing complexing agent are regularly distributed
3 It is possible that a regular Turing pattern forms when the underlying distribution of C is inhomogeneous. This Turing pattern can in general have a different Turing length than that for a homogeneous C distribution if 17 is not greater than the bifurcation value everywhere in the domain.
Z-P. Vorrmey et al./Physica D 99 (1996) 303-317
308 P¥ 1.0 0.8 0.6 0.4 0.2 0.0
r1
600 0
"
0.0
1.01Pv
(b)
0.8 i
0.6 0.4 0,2 0.0' 00.0 11,0
r
(c) 0.8 0.6 0.4 11.2
0.0 , ~
800 0"
3.0
Fig. 3. Concentration of Y for six different times, separated by 258 time units, beginning at time 5 000000 as a function of r for s "= 15 and I = 15 with L = 900. (a) Turing pattern: the concentration is time independent and all six curves are identical. The "luring wavelength is ~-T = 30. For reference, the period of the limit cycle is 1290 time units when the C species is absent. (b) Wave of excitation in the concentration of the Y species. (c) Standing wave: The spatial modulations in the heights of the waves that occur in different regions that lack complexing species decrease in time after many cycles of the standing wave. The time step and the scaled diffusion coefficient are At = 0.04 and D A t / A x 2 = 0.0833, respectively•
J.-P. Voroney et aL/Physica D 99 (1996) 303-317
309
I :"
.
:"
'
.:i
0.8 0.6 0.4 0.2
Cvv 0 -0.2 -0.4 -0.6 -0.8 "1
Fig. 4. Correlation functions Cyy(r ! ) versus rl for three sets of simulations with ~ - d - " 30. Twenty realizations of the dynamics were used to compute the correlation functions. Dotted line: s = 17; dashed line: s = 15; solid line: s = 13.
distribution of the species C that produced the steady patterns in (a) can produce traveling (b) or standing (c) chemical waves of excitation for different initial distributions of the X and Y species. Since the system exhibits a high degree of multistability for a given inhomogeneous initial distribution, the spatial patterns are best analyzed in statistical terms. For this reason we compute the spatial correlation functions by averaging over the initial distributions (of X and Y) that characterize the system state. (If the initial distributions of the C species are random then averages over the random distribution must also be considered.) The spatial correlation function is defined as T+to
Crr(r) m(T-i f
dtL -n
/0
xf
dr'Spr(r' + r, t)dp~,(r', t)) ,
(8)
I"
where (.. ")r denotes an average over the initial distribution and 8pr(r, t) = pr(r, t) - / S t is the deviation of Pr (r, t) from its space-time and ensemble average. Here to is a transient time and T is a time long compared to typical decay or oscillation times in the system. For the one-dimensional system, Fig. 4 plots the
normalized correlation function Crr for the (s, I) pairs (13,17), (15,15), and (17,13). For s = 17 all realizations show Turing pattern formation and the correlation function possesses long-range order with period ~.T. The correlation functions for s = 13 and s = 15 show damped oscillatory decay indicating disruption of the Turing pattern by traveling wave states and consequent destruction of the long-range order. One sees for s = 13 that positions separated by a multiple of ~.T are positively correlated at small separations and negatively correlated at large separations. Notice that a phase shift in the correlation length at large separations has taken place due to the traveling waves. Since the steady state is linearly unstable to perturbations with wavelengths between 0.64~.T and 2.03~.T, one might expect that Turing patterns with wavelengths other than the intrinsic one may form if the complexing species is distributed so as to favor the growth of perturbations of some other wavelength in this range. When the complexing species was distributed in an array with ~ d "- 2~-T and $ / ~ d - - 0.8, simulations show that a Turing pattern forms which has a wavelength that varies from 0.96~.T to 1.16~.T. We now consider a number of representative cases where s and ! vary. These results are shown in Figs. 5(a)-(c). Fig. 5(a) illustrates oscillatory dynamics when ~.d = 35 and ! = 20. Note the two pairs
J.-P Vomney et a t / P h y s k ' a D 99 (1996) 303-317
310
t' ....
(a)
0, 0.8 0.7 '
i
0~6
p,y
o.5 O.4 0.3 0.2 0.1
of steady peaks; other realizations have one or three pairs of steady peaks. Oscillations and Turing patterns for s = 35 and I = 5 are shown in (b). The number of oscillating regions varies from one realization to another. In the same realization, the positions of the oscillations do not change in time, but their amplitudes do. The transition to a Turing pattern is very sharp when ~.d = 40. At s = 36 no simulations exhibit oscillations, whereas at s = 35, all simulations contain oscillations. Fig. 5(c) illustrates a regular distribution of oscillatory and steady regions for a simulation with s = 15 and ! = 5.
0.0 0.0
(b)
200.0
400.0
600.0
800.0
o..~
3. 2. Regular distributions: two-dimensional systems 3.2.1. Striped distributions
0.4
0,3
Pv
oo 0.0
200.0
1it[!ti !j, lIL!Ij!,1 400,0
800,0
600,0
(c) 0.4
0.3
Pv
When the complexing species is distributed on m stripes with the same width and distance apart, oriented along the r2 direction, the symmetry group of the reaction-diffusion boundary value problem is Din~2x 0(2) where Dm/2 is the dihedral group of order m generated by an element that corresponds to a shift in the rl direction of 2L/m and the flip xl in the rl direction. The symmetry of solutions cannot be greater than Din~2x O(2), and a hexagonal pattern of dots with the Turing wavelength does not have a symmetry that is a subgroup of Dm/2 x 0(2) for most system sizes; however, the patterns can still have a local hexagonal symmetry. If the stripes are approximately one quarter of the intrinsic Turing wavelength, as in Fig. 6, then a striped pattern forms with twice the wavelength of ihe C distribution. Initially, a Turing l~attern begins to form on the stripes containing C, but as the local regions of activation grow larger than the width of the stripe, they each excite a chemical wave. This wave splits and
0.2
0.1
:HI'r'! y ~f rfrili:'I!~ [~ [] ii ~"i nr rIl,!,"-!
Id,1,~!!itllttl! !iul!~I1!i If!H~~riIi[l~!liII~N~rJI~!!,ii I!~ ~I,!I o.olHl,Jt~lllIt!,ttltl,~lj!ltlttl!lj!mj,!~!tl,~il 'i
0,0
li u '
200.0
: ,,
ili li 400.0
q
i
t
'~ Ilil i ~ ,
600.0
800.0
Fig. 5. (a) Oscillations for s = 15, ! = 20, and L = 840. (b) Oscillations and Turing patterns for s = 35, l = 5, and L = 840. (c) Regular distribution of oscillatory and steady regions for a simulation with s = 15, I = 5, and L = 840. The bottom line in the figures shows the location of the complexing species and the top lines show the Y concentration for six different times separated by 258 time units beginning at time 5 000000. The time step and the scaled diffusion coefficient are At = 0.04 and D A t / Ax 2 = 0.0833, respectively.
J.-P. Voroney et al./Physica D 99 (1996) 303-317
C
311
y
Fig. 6. The concentration of ¥ as gray shades is shown for five times separated by 55 time units beginning at time 1000000. Time increases from left to right and top to bottom. The first panel (upper left) shows the initial distribution of C on 16 equally spaced stripes (m = 32) with widths s = I = 4 ~ ~ kT. The simulation shows periodic dynamics with a small drift of the horizontal line of "luring dots. The time step and the scaled diffusion coefficient are At = 0.04 and D A t / A x 2 = 0.0226, respectively, and L = 128.
travels in two directions along the stripe without C, and locally activates the adjacent stripes which contain C. Eventually each of the two waves meets another wave or is absorbed by the stripes containing C. After several oscillations, the waves on the stripes become almost in phase, and occur only on every second stripe that lacks C. The Turing structures disappear except along one horizontal line that drifts slowly upwards in time. This interesting type of behavior occurs because the thickness of a stripe containing C is not large enough to support the formation of an entire Turing "dot", nor is the separation between stripes small enough to suppress oscillations. This phenomenon is a two-dimensional version of that presented in Fig. 5(c) for a one-dimensional system. Fig. 7(a) shows the results of a simulation where m = 18 and s = l "~ ½~.T. A Turing pattern has formed with its symmetries and wavelengths adjusted so that areas of activation fit in locations where there is complexing agent. The chemical pattern forms along each stripe containing C. Since I ~ ½~.T, the patterns on each stripe are in "communication" with each other,
and a global symmetry develops. This symmetry is not /)6 because the pattern is not symmetric with respect to flips of l~r or ~a'. After a transient period, where the system evolves to a state which has some defects in the global hexagonal symmetry seen in Fig. 7{a), oscillations suddenly begin. The oscillations arise from Turing dots which are near defects in the pattern and. spread into the regions with no C. The oscillations have the same period as the limit cycle which is stable in the absence of the complexing species. Fig. 7(b) shows a simulation where there are four stripes and s = i is large compared to kT. Now Turing patterns form within the large stripes containing C and chemical waves form in regions without C. The chemical waves are initiated from the edges of the Turing pattern where distortions in the regular Turing structure exist. 3.2.2. Hexagonal and square arrays Since a Turing pattern with hexagonal symmetry is selected for the above parameters when I/ = const., it is interesting to consider the modifications in the
312
Z-P Voroney et aL/Physica D 99 (1996) 303-317
Fig. 7. (a) The concentration or c at time 0 (left panel)and ¥, at times 150000 and 1000000 for a simulation with m = 18 (9 stripes with C), s -- ! - 7, and L - 126. (b) The concentration of C at time 0 (left panel) and ¥ at times 250000 and I 000000 for a simulation with m -- 4 (two stripes with C), s - I - 32, and L - 128.
Fig. 8. Formation of a spir ted by ½~.T. In front of the wave a Turing pattern with half the ~ r i n g wavelength f o ~ s on the disks with complexing species~ and in the wake of the wave a pattern with the Turing wavelength appears in the regions without complexing species. The panels show the Y concentration as gray shades for several times.
Z-P. Voroney et al./Physica D 99 (1996) 303-317
pattern when the underlying distribution of C is comprised of a hexagonal array of two-dimensional disks with separ,~tion different from that of ~.T. Fig. 8 presents the results of a simulation where the disk spacing in the hexagonal array is ½kT. A spiral wave has formed in the regions with no C. In front of the wave, a Turing pattern with half the intrinsic wavelength emerges on the disks containing C. In the wake of the wave, a pattern with the intrinsic Turing wavelength emerges in the regions without C. In much the same way that Turing patterns with symmetries other than hexagonal can form in systems whose dimensions are not much larger than the Turing wavelength, see [12] for example, when the length scales associated with the geometry of the complexing species distribution are near or not mach larger than the Turing wavelength, the pattern which forms in systems where the complexing species is distributed inhomogeneously can locally have non-hexagonal symmetries. In Fig. 9 the size of the disks with complexing species is not large enough to support an entire hexagonal structure. The Turing pattern is forced locally to have pentagonal symmetry; this symmetry permits the greatest number of dots of activation of the Y species to be present while maintaining the separation required by the long-range inhibition. Chemical waves pass through the regions without complexing species as in previous simulations. If the disks with complexing species are arranged in a square array with the separation between centers of the disks equal to ~.T, then because the points of activation of the Y species are forced to occur on disk containing C a Turing pattern forms which has square symmetry, as in Fig. 10. A wave of excitation with a well defined direction of propagation moves in the region with no C. 3.3. RasMom distributions
Finally, we briefly consider random distributions of disks of diameter s. To generate the initial distribution of C we randomly seed the system with disks with a given probability p. The disks may overlap and, depending on p, a certain fraction of the medium is randomly covered with complexing agent. We may
313
~!ii!!!~ii!i~iii~!)~i/!il ! i i~ i¸~ ~ ~i)!~i~!~i!!~i!~i~i!~!i~ii iil) ii!ili ~iii!~ ii¸¸ii~i~i~i~ i!i ¸:¸¸¸¸¸~¸¸
Fig. 9. Snapshot of the concentration of V at time 250 000 on a 252 square lattice for a simulation where the complexing species was initially distributed on disks in a square array.. The size and separation of the disks was large compared to ~.T, but each of the disks containing C could not support the formation of an entire hexagonal structure, so structures with pentagonal symmetry develop. The time step and the ,scaled diffusion coeffioent are At = 0.04 and DAt/Ax" = 0.00565, respectively.
then investigate the nature of the pattern formation process as a result of the fractional coverage of the system. Clearly, if the coverage is low one expects the spa- • tiotemporal dynamics to be dominated by the traveling waves in the oscillatory medium with no C, while if the fractional coverage is large enough that most of the system is covered and the gaps between disks are small then normal Turing pattern formation should result. In view of the nature of the seeding process the initial distribution has no special symmetry properties. Figs. I l(a) and (b) show two such simulations, one for a relatively low area fraction, (a) 0.41, and the other for a considerably larger area fraction, (b) 0.71. In (a) the pattern is dominated by complex traveling waves with small local regions showing incipient Turing pattern formation. In contrast, in (b) a recognizable Turing pattern with defects has formed. Wave propagation occurs in the regions where the species C is absent.
J.-P. Voroney et aL/Physica D 99 (1996) 303-317
314
Fig. I0. Snapshots s was initially distributed on disks in a square array. The separation ~ t w ~ n centers of the disks was approximately equal to ~-T and s = 18, I = 7.
X
(b)
Fig. !1. Snapshots of the concentrations of C'. X and Y for two different seeding probabilities of randomly placed disks cCntaining C. The fraction of the area covered by disks is (a) 0.41 and (b) 0.71.
4. Conclusions The simulations presented in this paper demon~trated that the distribution of the immobile species can influence the Turing and 'luring-like pattern formation processes. The (possibly)externally controlled lengths and symmetries can interact with those intrinsic to the reaction-diffusion system to give rise to spatiotemporal dynamics that are different from that seen in homogeneous media. Our results suggest that some of the interactions between oscillatory behavior and Turing patterns presented in this paper could be
observed experimentally in the CIMA and other reactions. Different gels, used in the expeliments to damp out convection, have different complexing properties. The transition to oscillatory behavior has already been detected [15] in agarose gels with low concentrations of starch, used as an indicator of iodide and a complexing reagent. In addition it has been observed that waves form and spread symmetrically from defects in the gel or in some cases from Turing spots, believed to be the results of a Turing-Hopf bifurcation. Thus, with the correct spatial distributions of agarose, polyacrylamide, and starch, results similar to
315
J.-P. Voroneyet aL/Physica D 99 t!996) 303-317
those presented here may be experimentally observable. In principle, one should be able to design gels with distributions of complexing agent to probe the phenomena reported in this paper. While our simulation results were confined to a modification of the Serkov model, the l:kenomena described above should occur in a wide class of activator-inhibitor systems in view of the general features required for their appearance, namely, a complexing reaction to an immobile complexing agent. Most naturally occuring media are heterogeneous. The present work shows if the scales of the inhomogeneities are macroscopic, albeit small, and comparable to those intrinsic to the system then Turing pattern formation different from that in homogeneous media can result. If the scales of the inhomogeneities are small compared to these macroscopic lengths and occur on microscopic or mesoscopic scales then diffusion can compensate for. these spatial variations and normal Turing pattern formation occurs. This investigation demonstrates that the distribution of complexing agent can serve as useful probe and source of complex spatiotemporal dynamics in chemical pattern formation processes. Inhomogeneous distributions can be designed in the laboratory to experimentally examine these phenomena. They are also surely relevant for some large-scale biological pattern formation processes where the systems are intrinsically heterogeneous, assuming that Turing bifurcations play a role.
Acknowledgements This work was supported in part by grants from the Natural Sciences and Engineering Research Council of Canada (RK and AL), by a Killam Fellowship (RK) and a Natural Sciences and Engineering Research Council of Canada postgraduate scholarship (JV).
Appendix A Assume that r/(r) = const, and (2) has a unique, spatially independent, fixed point, p*, a solution of f(P*x, P~') = g(P~, P~,) = 0 and h(p~,, p~., rl) = O.
The addition of a complexing reaction does not change the equilibrium values of the concentrations of X and Y species, but it can alter their stability. General discussions of the type of ~,tability analysis performed here can be found in [16,17]. We shall show that if the fixed point of (2) is stable or a s~.~dle, then its stability :emains unchanged by the addition of complexing species, but that if the fixed point of (2) is an unstable node or focus, then under certain conditions a Hopf bifurcation occurs with a Turitig pattern expected for i/above a critical value and a limit cycle expected for 17below this critical value. Let L be the matrix obtained fi'om linearizing (2) about the fixed point p*,
raf
of
Opx L=
Og
Opr ~
Og
0
\o
Oh
Opr
[}PC
Oh
Oh
Opy
OPc p=p,
a12 a22+bl bt
Oh
+
~PX Opr
all t --,a21
o
0 ) b2 ,
(A.I)
b2
then the characteristic equation of L is ~.3 _ t~2 -1- Sk - 8 = 0,
(A.2)
where r is the trace, 8 is the determinant, and S is the sum of the minors of the diagonal elements of the matrix L: r = ra + bl + b2, 8 = b28a, and S = 8a + alibi + v:ab2, where ra = all + a22 and 8a = a ! Ia22 - ai2a21 • For complexing reactions of the type considered here bl and b2 are always negative, and only bl depends on r/. The variables ra and 8a were chosen because the eigenvalues of the two-species system in the absence of immobile species are solutions of ~.2 _ ra~. +80 = 0. The two species system is stable when ra < 0 and 8a > 0, and the three-species system is stable when, using the Routh-Hurwitz conditions, r < 0, Sr - 8 < 0, and 8 < 0, are satisfied. It is straightforward to prove tha' if the two species system is stable then adding complexing species, which corresponds to decreasing bl, will not change the dynamics
316
Z-P. Voroney et al./Physica D 99 (1996) 303-317
of the system. 4, Decreasing bl can cause r to become negative when ra is positive so the addition of a sufficient amount of complexing species can stabilize a node, but since 8a < 0 implies 8 > 0 a saddle cannot be stabilized by adding C. If an unstable node of (2) changes its stability and becomes stable to spatially homogeneous perturbations, then L either has a zero or complex conjugate pair of eigenvalues. Since the determinant of L is independent of the bifurcation parameter 11 and is generically non-zero, the constant term of (A.2) cannot be zero which implies that zero is not an eigenvalue. Hence, by increasing r/the fixed point of (2) becomes stable when a cemplex conjugate pair of eigenvalues -+-ioJ of the matrix L crosses the imaginary axis and the third eigenvalue r is negative• In this case, the characteristic equation of L is of the form ~.3 _ r~.2 + oj2~. _ rw2 = 0. This implies that t S = 8, r < 0, S > 0, and 8 < 0 at the point where the system becomes stable to spatially homogeneous perturbations. Since this occurs with a single complex conjugate pair of eigenvalues crossing the imaginary axis, a limit cycle is the expected steady state when the concentration of the complexing species is below a critical value determined by r S = ~. A Turing pattern can form in the region of parameter space where a fixed point p* that is stable to spatially homogeneous perturbations loses its stability to spatially varying perturbations. The matrix La = L - k 2 D must have a positive eigenvalue, where D is the diagonal matrix whose elements are Dx, Dr, and 0. (Here we allow for unequal diffusion coefficients.) In order for an eigenvalue of Ld to be positive one of the inequalities, rd < 0, 8,t < 0, or Sd r d -- 8d < 0,
(A.3)
must be violated, where ra = r - k 2 Dx - k 2Dr, 8d = 8 - k 2Dxa22b2 - k 2Dra,,b2 + k 4Dx Di, b2, and Sd = 4 We know ra < 0, 8a > 0, bl < 0, and b2 < 0 and must show r < 0, Sr - d < 0 and 8 < 0, Proof: 8 < 0 and r < 0 tbllow directly from their definitions. Sr - 8 < 0 is equivalent to (r~ +bl)Sa + (alibi + a l l b 2 +a2262)r < 0 which can only be violated if alibi + a l l b 2 +a2262 = ra(bl +b2) - a l l b 2 = alibi -4-tab2 > 0 which is not possible since ra < 0 implies at least one of all and a22 is positive.
S - k 2Dx (a22 -t-bl -I-b2) - k 2Dy(all A-b2) h-k 4 Dx Dr • When the fixed point p* is stable to homogeneous perturbations then r < 0, implying that rd < 0; thus, the first of the inequalities in (A.3) is always satisfied. If 8,/> 0, then ~ = b28d and b2 < 0 imply 8a - k2Drall - k2Dxa22 + k4DxDr < 0.
(A.4)
This condition is exactly the same as one of the conditions for the formation of a Turing pattern in a two-species activator-inhibitor model (see [12] for example). If ra < 0 and 8a > 0 a Turing pattern can form in the system with or without the complexing species, and no bifurcations occur with an increase of 11. If ra > 0 and 8a > 0 then the spatially homogeneous steady state of the subsystem without the complexing species is an unstable node that can become stable to spatially homogeneous perturbations upon the addition of sufficient complexing species but remains unstable to certain spatially inhomogeneous perturbations when diffusion is present. When ra > 0, Eq. (A.4) can be satisfied even when the diffusion coefficients Dx and Dr are equal. The last inequality in (A.3), Sdrd -- 8d < 0 may be violated if for some k
kODxDy(Dx + D r ) - k 4 l D 2 ( r - - a l l ) + 2DxDI, r + D2(all +b2)] +k2[Dx(S +a2262 + r 2 - a l l r ) + Dr(S + r" - a 2 2 r - b l r - allb2)] -St
+8=0.
(A.5)
The coefficient of k 6 and the constant term are always positive when the fixed point is stable. If either or both the coefficients of k 4 and k 2 are positive then there is exactly one negative root, - r , and either two positive roots or a complex conjugate pair of roots. If N, A, B, and C are the coefficients of the above polynomial, then there are three real roots if (2A 3 - 9 A B N + 27N2C) 2 _< 4(A 2 - 3NB) 3. This is equivalent to - 1 8 A N B C + 27N2C 2 + 4 N B 3 < AeB 2 - 4 A 3 C . If we think of Dr as a bifurcation parameter, in systems where ell + b2 > 0 two positive roots emerge when the imaginary part of complex conjugate pair goes to zero, leaving a real positive double root d. When
J.-P. Voroney et aL/Physica D 99 (1996) 303-317
this happens, the characteristic equation has the form ).3 _ (2d - r)). 2 + d ( d - 2r)). + r d 2 = O. In the above discussion we have considered systems where all is negative. One may observe a bifurcation from a steady state to a Turing pattern by choosing all > 0 and D r much larger than D x . This corresponds to long-range inhibition with the inhibitor complexed. If the rate of the complexing reaction is increased, a critical value of the rate may be found where the Turing pattern disappears and a stable steady state is obtained.
References
Ill A.M. Zhabotinsky, in: Chemical Waves and Patterns, eds. R. Kapral and K. Showalter (Kluwer, Dordrecht, 1995) p.401; J.A. Sepulchre and A. Babloyantz, ibid, p.191 and references therein; O. Steinbock, A. Trth and K. Showalter, Science 267 (1995) 868; M.D. Graham, I.G. Kevrekidis, K. Asakura, J. Lauterbach, K. Krischer, H.-H. Rotermund and G. Erti, Science 264 (1994) 80; M.D. Graham, M. Biir, I.G. Kevrekidis, K. Asakura, J. Lauterbach, H.-H. Rotermund and G. Ertl, Phys. Rev. E 51(6) (1995); Rotermund and G. Ertl, Science 264 (1994) 80; A.T. Winfree, When Time Breaks Down: The Three Dimensional Dynamics of Electrochemical Waves
317
and Cardiac Arrhythmias (Princeton University Press, Princeton, NJ, 1987). G. 121 Papanicolaou and X. Xin, J. Slat. Phys. 63 (1991) 915; J.X. Xin, J. Slat. Phys. 73 (1993) 893; J.X. Xin and J. Zhu, Physica D 81 (! 995) 94. 131 D.L. Benson, J.A. Sherratt and P.K. Maini, Bull Math. Biol. 55 (1993) 365. 141 A. Turing, Phil. Trans. R. Soc. London B 237 (1952) 37. I51 V. Casters, E. Dulos, J. Boissonade and P. DeKepper, Phys. Rev. Lett 64 (1990) 2953. 161 Q. Ouyang and H.L. Swinneyo Nature 352 (1991) 610. [71 I. Lengyel and I.R. Epstein, Science 251 (1991) 650. 181 A. Hunding and P.J. Sorenson, J'. Math. Biol. 26 (1988) 27. 191 J.E. Pearson, Physica A 188 (1992) 178. [10l E.E. Sei'kov, Eur. J. Biochem. 4 (1968) 79. !111 P.H. Richter, P. Rehmus and J. Ross, Prog. Theor. Phys. 66 (1981) 385. !121 J.E~. Murray, Mathematical Biology (Springer, New York, 1989). II31 S. Wiggins, Introduction to Applied Non-Linear Dynamics and Chaos (Springer, New York, 1990). I141 M. Golubitsky, I. Stewart and D.G. Schaeffer, Singularities and Groups in Bifurcation Theory, Vol. 2 (Springer, New York, 1988). [151 K. Agladze, E. Dulos, and P. De Kepper, J. Phys. Chem. 96 (1992) 2400; P. De Kepper, J.-J. Perraud, B. Rudovics and E. Dulos, Int. J. Bifurc. and Chaos 4 (1994) 1215. !16l P. Borckmans, G. Dewel, A. De Wit and D. Walgraef, in: Chemical Waves and Patterns, eds. R. Kapral and K. Showalter (Kluwer, Dordrecht, 1995) p.323. !171 J.E. Pearson and W.J. Bruno, Chaos 2 (1992) 513.