Two-dimensional spectrophotometry of spiral wave propagation in the Belousov-Zhabotinskii reaction

Physica 24D (1987) 87-96 North-Holland, Amsterdam

TWO-DIMENSIONAL SPECTROPHOTOMETRY OF SPIRAL WAVE PROPAGATION IN THE BELOUSOV-ZHABOTINSKII REACTION II. GEOMETRIC AND KINEMATIC PARAMETERS

Stefan C. MULLER, Theo PLESSER and Benno HESS Max-Planck-lnstitut fiir Erniihrungsphysiologie, Rheinlanddamm 201, 4600 Dortmund, Fed. Rep. Germany Received 9 June 1986

Spectrophotometric data of the two-dimensional concentration distribution of the reaction catalyst ferroin/ferriin are analyzed for the spiral wave of chemical activity propagating in the Belousov-Zhabotinskii reaction. The Archimedian spiral and the involute of a circle fit equally well the iso-intensity levels corresponding to the highest and lowest values of ferriin concentration. Initially, the spiral pitch is 1.19 mm and the revolution time is 17.0 s, From these data a wave propagation velocity of 70 / t m / s is calculated. During 5 min the pitch and the period increase slightly, whereas the velocity remains approximately constant. The time course of these parameters is described by second order polynomials. The time evolution of the parameters agrees with results that were previously obtained by local evaluation methods. Ferriin concentration profiles are fitted by splines. From their derivatives concentration gradients are estimated to reach 5 mM/mm. The transition from a quasi-stationary state at the spiral center to the full wave amplitude outside the core is quantified.

1. Introduction

The Belousov-Zhabotinskii reaction displays a large variety of structural and dynamic features. Among these the propagation of spiral-shaped waves of chemical activity in thin excitable reaction layers has received special attention, ever since it was discovered about 15 years ago [1, 2]. On the basis of the experimental data given in the literature various theoretical approaches for the explanation of this phenomenon have been put forward by general analysis of the underlying reaction-diffusion partial differential equations [3-10] as well as by elaboration of models specific for the Belousov-Zhabotinskii reaction [11-14]. For reviews see [15-17]. Recently, the experimental knowledge of the properties of such spiral waves has been increased by combining optical methods with computerized video techniques [18], as reported by the authors [19]. In a preceding publication [20] a detailed two-dimensional description of concentrations of the reaction catalyst and indicator ferroin (or its

oxidized form ferriin) and of the methods for representing the two-dimensional digital data were given. In this article an image sequence covering a 5 min time interval of spiral wave propagation is evaluated by fitting spiral functions to the measured data. These functions do not imply any assumptions about mechanistic models. The analysis remains descriptive in the sense that the evaluated parameters of the applied functions are to be considered as a global representation of relevant geometric and kinematic properties of the measured data set. After briefly introducing the methods for obtaining global parameters from the recorded time sequence of digitized images, it is shown in the result section that, within experimental error, a fit of an Archimedian spiral equals that of an involute of a circle. The time dependence of its parameters is determined and detailed information about the concentration gradients occurring in the spiral structure is derived from splines fitted to the extracted concentration profiles.

0167-2789/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

88

S.C. MiJller et aL / Geometric and kinematic parameters of spiral waves

2. Evaluation of image data 2.1. The data set and the fitting procedures The experimental data used in the following evaluations are provided by the sequence of 100 digital images that were taken at 3s intervals during the propagation of a spiral wave described in detail in [20]. An example is shown in fig. 1. Each image represents a 4.5 x 4.5 mm 2 section of the reactive solution layer of 1 mm thickness and consists of 450 x 450 pixels with a grey level resolution of 256 digital units each. The techniques of representing the information contained in these i m a g e s - averaging methods, three-dimensional perspectives, pseudo-color transformations-are given elsewhere [20, 21]. For the current fitting purposes, the image data are taken in their original non-averaged form but normalized such that an average grey level of 128

is obtained [20]. One of the main objectives is to extract the geometric form from those sets of pixels having one specific grey level, preferentially the maximum or minimum level. The Archimedian spiral and the involute of a circle will be compared with these data. The time evolution of the parameters of both spiral functions is derived for the complete sequence of images taken during the 5 min interval. Furthermore, the distribution of measured concentration values along a specific direction in the observation rectangle is fitted by splines [22] from which the concentration gradients are calculated.

2.2. Spiral functions In a first approach to a geometric description of the experimental data an Archimedian spiral [23] is considered. It is given by the equation r = a(qg-- q%),

withqD>qo 0.

(1)

The polar coordinates r and ~0 have their origin at the center of the spiral with the cartesian coordinates x o and Yo, as sketched in fig. 2A. The angle ~, the independent variable, is incremented by 2~r for each whorl of the spiral. The parameter a determines the distance between successive whorls or the pitch A = 2~ra.

Fig. 1. Digital image of the spiral wave in the BelousovZhabotinskii reaction as already shown in fig. 2A of [20]. In this nonaveraged but normalized version of the original measurements, the grey levels of transmitted light intensity vary between 92 and 168. The maximumvalue is indicated by black dots along the crest of the spiral, the minimum value by white dots along the trough between the spiral crests.

(2)

The angle ¢Po describes the angular position of the spiral. This parameter varies with time ( % = % ( t ) ) because the spiral rotates. In this instance qOo(t) tends to more and more negative values, since for data collection the clockwise rotating spiral was selected from a pair of spirals. A second approach, supported by some plausible reasoning [2], is given by the involute of a circle [23] with radius a according to - (a/r)2

= a {• _

+ arccos (a/r)

- - > }, (3)

withr>a,

q~>~o +½~r.

S. C. Miiller et al. / Geometric and kinematic parameters of spiral waves

yT

A

/,"," /"

well as 2B are drawn with the same set of spiral parameters. The quantitative estimation of the spiral parameters from the experimental data was performed by the nonlinear least squares fitting procedure VA05A from the Harwell Subroutine Library [24]. Since the pixel locations extracted from an image are given by the cartesian coordinates x and y (fig. 1), the transformation into polar coordinates has to be included into the optimization process for eqs. (1) and (3),

r= ~(X- Xo)2 + (y-yo) 2, (Y-Yo]

Ixo

89

(4)

cp = 2~rn + arctg _--:-~o x 1.

yl

The integer n = 0,1, 2 . . . counts the whorl close te which the pixel with the coordinates (x, y) is located.

[3

3. Results 3.1. /"

txo X t-

Fig. 2. Sketch of an Archimedian spiral (A) and of an involute of a circle (B) according to eqs. (1) and (3), respectively.For explanation of the parameters see text.

This equation of a spiral for the region r > a has the same set of parameters as the Archimedian spiral. By addition of the term - ½~r in the bracket, eq. (3) is formulated such that for large r it transforms into eq. (1). For r = a one finds from eq. (3) cp = Cpo+ ½~r. Consequently, there is a difference of ½~r between the angle % and the angle of the radius vector for r = a being a tangent to the spiral, as sketched in fig. 2B. Both figs. 2A as

Fitting of spiral functions to iso-intensity levels

The digital image shown in fig. 1 was taken during the first revolution of the recorded spiral motion. The image is identical to that shown in fig. 2A in the preceding publication except that the spatial average of its grey level distribution is normalized to a value of 128 [20]. This operation does not affect the geometric properties of the structure. In addition, the pixels with maximum and minimum intensity are marked as black dots on a bright background and as white dots on a dark background, respectively. An operational definition of the minimum level is given by the grey level up to which the integral of the intensity distribution accumulates 1% of the 4 5 0 x 4 5 0 pixels of the full image. Correspondingly, the maximum intensity is defined as the grey level up to which 99% of the total number is accumulated. These levels, as marked in fig. 1, are 92 and 168, respectively.

90

S.C. Miiller et a l l Geometric and kinematic parameters of spiral waves

B

A

3

3 E E

v

~2

0

1

2

3

2

0

4

x

x (ram)

C

3

4

3

4

mm )

D

3

3 E E

E E

v

2

I

2

~2

0

1

2

3

4

x (ram)

0

1

2 x (ram)

Fig. 3. Fits of spiral functions (solid lines) to the pixels with maximum or minimum grey levels in the spiral pattern of fig. 1 (dots): (A) Archimedian spiral for maximum level; (B) Archimedian spiral for minimum level; (C) involute of a circle for maximum level; (D) involute of a circle for minimum level.

Since no spatial smoothing procedures were applied to the image, the pixel noise of the measurem e n t s is fully reflected b y the scattering of the pixel locations. Although the scattering is significantly smaller at the steep slopes of the wave fronts, the m a x i m u m and m i n i m u m levels have b e e n selected for the following evaluations, be-

cause for these extrema one can extract a data set for only one spiral. Each intermediate level produces a d a t a set which is c o m p o s e d of a spiral at the excited wave front, of a spiral at the relaxing back, and d a t a of the region connecting both. T h e results of a nonlinear least squares fit of an A r c h i m e d i a n spiral (eq. (1)) to the m a x i m u m and

S.C. Miiller et al./ Geometric and kinematic parameters of spiral waves

minimum level of fig. 1 are shown in figs. 3A and B. The spiral parameters for the maximum level having 671 data points are: a = 0.190 mm, Cpo-1.187r, x o = 2.23 mm, Yo = 1.98 nun. The standard deviation of the fit is + 0.038 mm. Those for the minimum level containing 831 data points are: a = 0.190 mm, Cpo= - 1.76~r, x o = 2.24 mm, Yo --" 1.98 mm. The standard deviation of the fit is _ 0.051 ram. The difference between the standard deviations of the two levels becomes clear from inspection of fig. 1. The data at the minimum are more scattered since this level is located at the b o t t o m of a broad trough. Obviously, the calculated curves fit the measured data points of constant intensity well. Systematic errors are caused by distortions in the camera target already mentioned in [20]. These are responsible for small deviations (up to 2%) between the data and the fitted spirals with respect to the x- and y-direction. (In fact, the measured spiral points appear to be somewhat "squeezed" in the x-direction, but this can be considered to be of negligible influence on the fitted parameter values.) The fits of an involute (eq. (3)) to the same data are shown in figs. 3C and D. The fitted curves are almost indistinguishable from the corresponding curves in plots A and B of the same figure. The involute parameters are a = 0.188 mm, Cpo= -1.231r, x 0 = 2.23 mm, Y0 = 1.98 mm (fig. 3C) and a = 0.189 mm, cp0 = - 1 . 8 0 ~ r , x 0 = 2.24 mm, Yo = 1.98 mm (fig. 3D). The standard deviations of the fits for the maximum (+0.038 mm) and for the minimum level (+0.051 mm) agree with those obtained for the corresponding Archimedian spirals. The quality of the fits of Archimedian spirals and involutes of a circle shows that with the currently available experimental techniques no decision can be made which of the two alternatives is the more appropriate one for the investigated system. At this point it should be emphasized that all of the parameters of the two spiral functions are fitted independently, neglecting the constraint that the functions should coincide for large spiral radii

91

5

v

2

1

2

3

4

x(mm) Fig. 4. Comparison of an Archimedian spiral and an involute of a circle with the same coordinates x o and Yo of the center and the same parameters a and ¢Po.

r. In fig. 4 it is illustrated how the involute differs from the Archimedian spiral if the numerical values of the parameters are the same in both functions (x o = 2.25 mm, Y0 -- 1.98 mm, a = 0.190 mm, Cpo= -1.25~r). The systematic deviations between the two curves within a region close to the center cannot be resolved by the experimental equipment. 3.2.

Time dependence of the spiral parameters

The data evaluation procedures presented in the preceding section was applied to the entire series of one hundred images. This yields the temporal evolution of the four spiral parameters a(t), Cpo(t), xo(t ), and yo(t). The four parameter sets describe the time dependence of the Archimedian spiral and the involute of a circle for both the minimum (1%) and the maximum (99%) intensity levels. The results are presented in figs. 5 to 7. Fig. 5 indicates that the pitch of the spiral, A(t)= 2~ra(t), is a slightly increasing function of time. Ignoring th~ plateau between 0 and about 60 seconds a crude linear approximation of the values for the maxi. m u m (symbol O ) is given by the expression

A(t)=l.188+O.994xlO-4t

(mm).

(51

S.C. Mfiller et al./ Geometric and kinematic parameters of spiral waves

92 1.24

1.22

1.2 -

<

involute:

ooo-

I • "

0

00~

0 00

°'o

. . . . . "°2 °~* °* * ~ " .~%"o° °%,°° % ~°.~ .~ ,,, ~°~. ,÷ ,

0

. 60

.

. 120

. 180

. 240

300

Time ( s )

Fig. 5. Temporal evolution of the pitch A = 2~ra, determined by fitting the Archimedian spiral to the maximum (O) and minimum intensity levels (+) in the sequence of 100 spiral images. Previously reported data obtained by local methods (fig. 10 of [20]) are included (11).

A has the same trend as previously found by local methods [20] (symbol II). The time dependence of the spiral orientation given by the parameter ~0 in eqs. (1) and (3), respectively, is shown in fig. 6. For each full spiral rotation an increment of -2~r is added to the value of ¢Po found by nonlinear optimization. For convenience the absolute value of ¢p0/~r is plotted. Its variation with time can be well described by a parabola. The following relationships were found by the fitting routine VA05A of [24]. For maximum intensi, ty level: Archimedian: Cpo/~r = - 1.75 - 0.1176t + 1.063 X lO-St 2

Here again it turns out that the Archimedian anc the involute give almost identical results. Frorr eqs. (6) and (7) follows a rotation period of T--17.0 s at the beginning of the experiment. Afte~ 300 seconds the period reaches 17.9 s. The quadratic term in the equations accounts for second order effects, presumably caused by the chemical reaction driving the solution from the reducec to the finally fully oxidized state. The accuracy ir determining the small second order term reflect,, the systematic precision of the experiment. The velocity of the moving wavefront is given for an Archimedian spiral by the time derivative of eq. (1) resulting in

o(t)

d--7

With eqs. (5) and (6a), this expression gives for the

40 32

_

involute:

dr =

(8)

(6a)

standard deviation: + 1.88 × 10-2;

~Oo/~r = - 1.795 - 0.1176t + 1.06 × 10-5t 2

24

--16

(6b)

standard deviation: __ 1.86 x 10 -2.

8 0 60

F o r minimum intensity level: Archimedian: ~Oo/~r = - 0 . 3 6 - 0.1176t + 1.0 × 10-5t 2 standard deviation: _+2.29 x 10-2;

(7b'

standard deviation: +_2.7 × 10 -2.

o"

1.18

1.16

~o/~r = - 0 . 3 9 7 - 0.1177t + 1.0 × 10-st 2

120

180

240

300

Time ( s )

(7a)

Fig. 6. Temporal evolution of the parameter q~o, determined by fitting the Archimedian spiral to the maximum ( O ) and minimum intensity levels ( + ) in the sequence of 100 images. The solid lines are fitted parabola (see eqs. (6a) and (7a) in the text).

S. C. Miiller et al. / Geometric and kinematic parameters of spiral waves

maximum intensity level

v(t)

secutive whorls of the spiral dA

= 69.9 × 10 -3 + 1.58 × 10-sip - 9.39 ×10-7t-1.58×10-9/2

02~

(mms-1).

-- V2~r(n-1) =

dt

(9)

= 0.994 × 10 - 4 Eqs. (8) and (9) show that the velocity is a function of time and of the direction of motion given by the angle ~. The direction dependence is a consequence of the time dependence of the pitch A. The sequence of wavefronts of the spiral along a fixed direction is counted by multiples of 2~r. With the expression % = ~0 + 2~rn one finds from eq. (8) for the difference in velocity of two con-

120

% @**~+********

A

I00

**** **

8O

•

60

• +.*

I,+ *

~ooo

ooo

oo •

aoo •

oooooo •

•

o

40 '°

~'~..,.*'/÷ **

20

• •

u*

+

°°o°

o o ~o

.+++ooo#o o o

0 -20 • -40 • 0

60

120 Time

240

180

300

( s )

120 100

B

8O 60 4.0

• +~ . ÷o**

~o 20 o~

+Lo

93

o *

** + o%00 + ÷ ~oo°o,.~°.°* .

•

(mm

s-l).

(10)

This shows that the leading wavefront moves slightly faster than its follower. The velocity of any wavefront remains approximately constant. Its value is lowered by about 0.6% within the 300 s time interval of the experiment as calculated by eq. (9). This result differs from observations on target patterns in the excitable [25] and oscillatory regime [26]. The time dependence of the spiral pitch A and the second order time effect of % indicate that the applied spiral functions, the Archimedian spiral as well as the involute of a circle, are reasonable approximations of the experimental data only as far as the time scales of overall chemical changes in the solution are much larger than the rotation period of the spiral. The location of the spiral center shows also some variation with time (fig. 7). The drift during a 300 s time interval is less than one tenth of the spiral pitch. The forces responsible for this movement are unknown but variations in surface tension due to temperature or chemical composition, leading to convective motion is one of the prospective candidates [27, 28]. Note that there is a distinct gap between x 0 - 20 values found for the maximum intensity ((3) and the minimum intensity level ( + ) . This difference may be due to the geometric distortion of the camera target.

0 -20 -40 0

3.3. 60

120 Time

180

240

300

( s )

Fig. 7. M o t i o n o f the center o f the spiral. Its coordinates x o (A) a n d Yo (B) are derived from fits of any of the two spiral f u n c t i o n s to the m a x i m u m ( O ) and m i n i m u m intensity levels ( + ) . T h e y are plotted relative to a reference position at )70 = 2.25 m m a n d )70 = 1.98 m m previously obtained for the first spiral r e v o l u t i o n b y local m e t h o d s [20]. The results for the s u b s e q u e n t revolutions, as s h o w n in fig. 10 of [20] are included

(-,).

Intensity profiles

In parts I and 2 of the result section the locations of those pixels were analysed which have a preselected intensity level or concentration of ferdin, respectively. For the analysis of intensity profiles the locations of the pixels are preselected by a given curve in the 2D plane. In a first instance the size of the core region can be estimated by a plot of the intensity or concentration

94

S.C. Miiller et al./ Geometric and kinematic parameters of spiral waves

~L

/

_L =EC_

Et

o~

T

T i

0.

.15

.45

.3

.6

.75

r(mm)

6

Fig. 8. Average variation of ferriin concentration along the route of fitted Archimedian spirals from the center to the crest as well as to the trough. The averageis calculated from the first fifteen images of the recorded sequence (see also fig. 9 in [20]).

4

~2 ×

0

~-2

values along an Archimedian spiral on its way f r o m the center to the outer region. The curves in fig. 8 represent the ascent along the Archimedian f r o m the ferriin concentration at the very center to m a x i m u m ferriin concentration and the descent f r o m the center to minimum ferriin concentration in the trough between two fronts. The independent variable of the plot is the length of the radius vector r describing the spiral (eq. (1)). Both curves in fig. 8 are obtained by a spline fit of the data extracted f r o m the first fifteen images of the recorded sequence. For comparison see fig. 9 in [20]. A radius of 0.35 m m is a reasonable estimate of the core size, a number which was also found by other means [19, 20]. The crossover of the two curves in the region close to r equal zero is due to systematic instrumental errors. The shape of the ascending curve supports theoretical results [29] but a quantitative characterization, for example in terms of relaxation constants, is not yet meaningful. An example of a concentration profile along a straight line is shown in fig. 9A (see also fig. 3 in [20]). The data are obtained from fig. 1 after a 3 × 3 moving average and are fitted by cubic splines with the routine VC03A of [24]. Fig. 9B displays the concentration gradients calculated f r o m the splines. The peak values of the steepness of the wavefronts amount to about 5 m M / m m .

-4 1 -6 O.

/5

~

~15 2

2Z5 3

315 ;

415

x ( rnrn }

Fig. 9. (A) Profile of ferriin concentration extracted from fig. 1 and passing in the x-direction through the center of the spiral at Y0= 1.98 mm as already shown in fig. 3 of [20]. The solid line is a spline fit of the measured data. (B) Derivative of the curve in (A) as calculated from the splines.

They exceed the values estimated by graphical methods [20] by a factor of 3. Disregarding the small wiggles there is a dearly detectable asymmetry in the derivative indicating an at least biphasic leading wave front. The wave backs show no distinct structural details as reported for circular waves [30]. The differences may be explained by the fact that the wave amplitude of the circular waves was b y a factor of three larger than that of the spiral waves analyzed in this paper:

4. Conclusions An analysis of the geometry and the kinematics of the spiral shaped chemical wave in the Belousov-Zhabotinskii reaction is presented. It uses digitized experimental data of the spatio-tem-

S.C. Miiller et al. / Geometric and kinematic parameters of spiral waves

poral distribution of the reaction catalyst and indicator ferroin for finding optimal values of the parameters of two spiral functions: the Archimedian spiral and the involute of a circle. Furthermore, profiles of ferriin concentration obtained along a line of given direction through the spiral center are represented by fitting splines to the measured data. The resolution of the image data allows for a significant estimation of the relevant spiral parameters. It turns out that no decision can be made yet whether an Archimedian spiral is more appropriate for the description of the spiral geometry than the involute of a circle. This calls for obtaining a still better spatial resolution of the core region which is one of the major objectives of current work by the authors. Both functions were used without further implications for specific model assumptions. As a second order effect we find that the spiral parameters are time dependent. Further analysis of this observation may lead to deeper insights into the relationship between the progress of the reaction and the excitability properties of the medium. While there exist various mathematical models on the geometric form of chemical spirals, much less has been proposed yet for the shape of the profiles of concentration waves. Our measurements provide experimental material on the detailed features of these profiles, and the fitting of their shape by splines is a first step for further theoretical analysis. It is noteworthy that information could be gained about local concentration gradients. They amount to 5 m M / m m at the steep wave fronts. With regard to the core, the applied evaluation procedures corroborate previous results for its spatial extent (maximum diameter = 700 /~m) and the small drift in space of its center location. Generally, the improvement of the experimental knowledge regarding the geometry and the kinematics of chemical spiral patterns by one or even more components of the reacting species is expected to motivate further theoretical efforts for

95

quantitative wave solutions of the appropriate sets of reaction-diffusion equations.

Acknowledgements We thank Dr. W. Kramarczyk for programming support and Ms. B. Plettenberg for typing the manuscript.

References [1] A.M. Zhabotinsky, Investigation of homogeneous chemical auto-oscillating systems, Thesis (in Russian), 1970, Institute of Biological Physics of the Academy of Science, USSR, Pushchino. [2] A.T. Winfree, Science 175 (1972) 634. [3] N. Wiener and A. Rosenblueth, Arch. Inst. Cardiol. Mex. 16 (1946) 205. [4] I.S. Balakhovskii, Biofizika 10 (1965) 1063. [5] J. Guckenheimer, Lect. Notes in Mathematics, vol. 525 (Springer, Berlin, 1976). [6] D.S. Cohen, J.C. Neu and R.R. Rosales, SIAM J. Appl. Math. 35 (1978) 536. [7] J.M. Greenberg, SIAM J. Appl. Math. 39 (1980) 301. [8] N. Kopell and L.N. Howard, Adv. Appl. Math. 2 (1981) 417. [9] J.P. Keener, SIAM J. Appl. Math. 39 (1980) 528. [10] P.S. Hagan, SIAM J. Appl. Math. 42 (1982) 762. [11] J.A. DeSimone, D.L. Beil and L.E. Scriven, Science 180 (1973) 946. [12] A.S. Mikhailov and V.I. Krinskii, Physica 9D (1983) 346. [13] A.B. Rovinsky, J. Phys. Chem. 90 (1986) 217. [14] J.P. Keener and JJ. Tyson, Physica 21D (1986) 307. [15] Springer Series in Synergetics, H. Haken, ed., Vols.: 5 (1979); 6 (1980); 12 (1981); 19 (1984); 27 (1984); 28 (1984); 29 (1985). [16] J.J. Tyson, The Belousov-Zhabotinskii Reaction, Lect. Notes Biomath., vol. 10, (Springer, Berlin, 1976). [17] Oscillations and Traveling Waves in Chemical Systems, R.J. Field and M. Burger, eds. (Wiley, New York, 1985). [18] S.C. Miiller, Th. Plesser and B. Hess, Anal. Biocliem. 146 (1985) 125. [19] S.C. Mtiller, Th. Plesser and B. Hess, Science 230 (1985) 661. [20] S.C. Miiller, Th. Plesser and B. Hess, Physica 24D (1987) 71. [21] S.C. Miiller, Th. Plesser and B. Hess, Naturwissenschaften 73 (1986) 165. [22] C. de Boor, A Practical Guide to Splines, Applied Mathematical Sciences, vol. 27 (Springer, Berlin, 1978). [23] I.N. Bronshtein and K.A. Semendyayev, Handbook of Mathematics (Harri Deutsch, Frankfurt, 1985).

96

S.C. Mfiller et al./ Geometric and kinematic parameters of spiral waves

[24] Harwell Subroutine Library (AERE Harwell, Oxfordshire, England). [25] K. Showalter, J. Phys. Chem. 85 (1981) 440. [26] M.L. Smoes, in: Dynamics of Synergetic Systems, H. Haken, ed. (Springer, Berlin, 1980), p. 80. [27] J.K. Platten and J.C. Legros, Convection in Liquids (Springer, Berlin, 1984).

[28] S.C. Mtiller, Th. Plesser and B. Hess, Ber. Bunsenges. Phys. Chem. 89 (1985) 654. [29] D. Walgraef, G. Dewel and P. Borckmans, J. Chem. Phys. 78 (1983) 3043; T. Yamada and Y. Kuramoto, Progr. Theor. Phys. 55 (1976) 2035. [30] P. Wood and J. Ross, J. Chem. Phys. 82 (1985) 1924.