Dalibor Štys [ stys@jcu.cz    Petr Jizba [    Anna Zhyrova [    Renata Rychtáriková [    Kryštof M. Štys [    Tomáš Náhlík [

Mesoscopic dynamics of self-organized structures is the most important aspect in the description of complex living systems. The Belousov–Zhabotinsky (B–Z) reaction is in this respect a convenient testbed because it represents a prototype of chemical self-organization with a rich variety of emergent wave-spiral patterns. Using a multi-state stochastic hotchpotch model, we show here that the mesoscopic behaviour of the non-stirred B–Z reaction is both qualitatively and quantitatively susceptible to the description in terms of stochastic multilevel cellular automata. This further implies that the mesoscopic dynamics of the non-stirred B–Z reaction results from a delicate interplay between a) a maximal number of available states within the elementary time lag (i.e. a minimal time interval needed for demise of a final state) and b) an imprecision or uncertainty in the definition of state. If either the number of time lags is largely different from 7 or the maximal number of available states is smaller than 20, the physicochemical conditions are inappropriate for a formation of the wave-spiral patterns. Furthermore, a white noise seems to be key for the formation of circular structures (target patterns) which could not be as yet systematically explained in existing models.

Belousov–Zhabotinsky reaction, hotchpotch machine, physical model, thermodynamics

ICPF]University of South Bohemia in České Budějovice, Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Institute of Complex Systems, Zámek 136, 373 33 Nové Hrady, Czech Republic \phone+420 38777 384 FJFI]Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University, Břehová 7, Prague 1, Czech Republic ICPF]University of South Bohemia in České Budějovice, Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Institute of Complex Systems, Zámek 136, 373 33 Nové Hrady, Czech Republic ICPF]University of South Bohemia in České Budějovice, Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Institute of Complex Systems, Zámek 136, 373 33 Nové Hrady, Czech Republic ICPF]University of South Bohemia in České Budějovice, Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Institute of Complex Systems, Zámek 136, 373 33 Nové Hrady, Czech Republic ICPF]University of South Bohemia in České Budějovice, Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Institute of Complex Systems, Zámek 136, 373 33 Nové Hrady, Czech Republic Hotchpotch model of the Belousov–Zhabotinsky reaction] Multi-state stochastic hotchpotch model gives rise to the observed mesoscopic behaviour in the non-stirred Belousov–Zhabotinsky reaction

1 Introduction

During the last two decades, a number of new paradigms for understanding complex living systems have emerged. These include, e.g. theory of dynamical systems, theory of complexity, nonlinear dynamics, evolutionary physics, and critical phenomena1, 2. Among these, chaotic attractors, (multi-)fractals, self-assemblies, dissipative structures and self-organization represent some of the most promising recent concepts. A particularly important testbed for a conceptualization of pattern formation in self-organizing systems is the B–Z reaction3, 4, 5, 6, 7, 8, 9 which belongs among the most extensively studied examples of chemical self-organization. However, despite decades of intensive research, there are still ongoing controversies over the actual chemical kinetics (i.e. details of rates of chemical reactions involved) and the mesoscopic dynamics (i.e. exact nature and mechanism of patterns formation at the mesoscopic scale) of the B–Z reaction10, 11, 12.

The B–Z reaction is considered as a textbook example of the so-called excitable medium13, 16, 14, 15. Majority of available chemical models which aim to explain the chemical self-organization are based on the standard reaction-diffusion analysis and the law of mass action applied to a few selected reactions17, 18. On the other hand, self-organization is a hallmark of a far-from-equilibrium dynamics which appears difficult to reconcile with a common-sense chemical reaction scheme based on the law of mass action.

Turing patterns that appear in some reaction-diffusion models are often considered as a theoretical embodiment of the B–Z patterns19. This is wrong for at least three reasons: Turing patterns (a) can explain only wave B–Z patterns but not spirals10, (b) appear only at specific parameter values in the reaction-diffusion equations, therefore they are unstable under fluctuations of parameters (in contrast to experiments where patterns are observably robust), (c) are stable solutions of Turing’s reaction-diffusion equations while, in the experiment, we observe a dynamic system on a trajectory through the state space towards a limit cycle with alternating spirals and wave fragments.

The mesoscopic description of the B–Z reaction is typically modeled with a cyclic cellular automaton which often generates patterns similar to the B–Z wave patterns found near to the final stage of the reaction. The morphological characterization of patterns is pivotal in these approaches while chemical aspects are often of a secondary interest.

So far, mesoscopic studies of the B–Z reaction have been limited to a low-level cellular automata20 which can reasonably well account only for some of the observed B–Z wave patterns, while it is as yet unclear how to generalize these approaches to obtain a full-fledged evolution of wave patterns together with dynamics of spiral patterns. To the best of our knowledge, the influence of the number of levels has been systematically studied only in one case20 while most of the systematic studies in the literature have been confined to maximally 8 levels21.

In order to numerically implement a hotchpotch model, we have adopted an approach based on the version of Wilensky NetLogo model15. In our case the model is limited to 200 achievable state levels and simulated on a square 50 50 grid. After a random setup of the space distribution of initial centers of as


where is the maximally achievable number of levels of the cell state. The model at each time step may proceed in four possible ways:

  1. When a cell is at the , so-called quiescent, it may be “infected” by surrounding cells according to the equation


    where and is a number of cells at the and , respectively, and are characteristic constants of the process and is a maximum allowed level of state/excitation.

  2. When a cell is at the , its new state is calculated as


    where is a state of the -th cell in the Moore neighbourhood, which directly surrounds the examined cell, and is another arbitrary constant. For simplicity’s sake, we denote in the following text and as 2a and 2b, respectively.

  3. When a cell is at the , then

  4. When a cell achieves the , then


Here we posit that our improved stochastic version of the multilevel hotchpotch model not only faithfully represents the dynamics of wave-spiral patterns of the B–Z reaction but also provides insight into underlying chemical reactions in terms of the Eyring activated complex theory22. The emergence of correct spatial structures in our model depends on the ratio between the number of available states within the elementary time lag and on the rate of the internal increase in excitation process. At appropriate levels of noise, we observe a dynamical evolution of the system through circular waves up to the final stage — spiral-wave interchange. This structure-stabilizing mechanism via the noise is a typical trademark of a mesoscopic description and is not simply attainable via microscopic deterministic rules. By the presence of the noise, we can explain many controversies in self-organization in the Nature.

2 Results

2.1 Modelling the Belousov–Zhabotinsky reaction in excitable media and constructive role of noise

The B–Z reaction is not easily comprehensible in terms of the standard law of mass action (which represents the “canonical method” of interpretation of the chemical reactivity). This is due to the fact that the reaction space is separated into regularly evolving/travelling structures and that one has to consider a large number of interlocked chemical processes involved. In this work, we report on a new cellular-automaton based stochastic model of the B–Z reaction. The model retains some of the key features of the multi-level hotchpotch machine which, however, outperforms both in its ability to faithfully mimic the onset stage of the B–Z reaction and its potential to correctly describe the morphology of the interacting wave-spiral patters during their evolution.

Figure 1a compares a late stage of the B–Z reaction (full data are accessible in Supplementary Material 1) in our least possible spatial constraining (a 200-mm Petri dish) and roily conditions (gentle mixing at 1400 rpm using an orbital mixer) with simulation of the modified Wilensky model of excitable medium15. The structures are, for some parameter ranges, astonishingly similar. However, the most regular spirals and waves, best comparable to the model, are expected to arise in a very gently mixed, homogenous solution of thin layer in a vessel of the unlimited size which does not spatially constrain the evolution of waves.

In order to achieve the manifest morphological similarity between the B–Z experiments and our simulation, we implemented the following changes into the Wilensky model:

  1. the enlargement of the cellular grid to 1000 1000,

  2. start from a very few points which enabled to analyze the behaviour of individual centers of emanation,

  3. a sequence of switching the values of cell states from natural to decimal numbers which extended the span of each cellular state,

  4. the addition of a uniform white noise to each automaton step which compensated for our limited knowledge of precise underlying mechanism, and

  5. the extension of the number of achievable states and rate of the internal cell excitation up to 2000 and 280, respectively, to smooth the model waves.

The first modification — usage of the finer grid — suppressed the influence of the nonidealities of the periodic boundaries on the evolution of the model system. The second intervention into the Wilensky model was performed using random-exponential function for generation of the starting (ignition) points. Multiplication of each cell state by the of the exponential distribution, i.e.


ensured that the simulation started with a small number of the ignition points.

At this initial condition, we first tested the influence of the of natural numbers due to the rounding in processes 1–2 (see Eqs. (2)-(3)), thus new states in time were calculated as




This modification, which was originally implemented to start the process from these few centers (ignition points) quite surprisingly increased the structural-patterns similarity between the B–Z experiment and the simulation. The results are shown in Figure 1b–c and Supplementary Materials 2–4 and were discussed in some detail in23. In Figure 1b, we present early simulation steps 2, 4, 14, and 16 after the ignition in process 1. For and , at least two non-zero points in a proper configuration , were required for the evolution of the waves in the simulation, since at least one addend in process 1 has to be equal to 1. In this case, the early evolution gave octagons (Figure 1b, i) and final state was formed by spirals (Figure 1c, i). In contrast, if and , then, e.g. and the non-zero cell was surrounded by evolving wave of 8 cells in . This early evolution resulted in squares with central circular objects (Figure 1b, ii) which further led to the filamentous structures (Figure 1c, ii).

The next step softened the definition of the state by allowing 1 decimal place in Eqs. (2)-(3), i.e.




and kept and (Supplementary Material 5). If , the first layer of 8 cells in evolves in the Moore neighbourhood. This modification provided the trajectory as observed for and with natural number states. More decimal places did not change the state space trajectory any more. This proves that the trajectory is determined by the ratio of the constant to the number of levels . The ratio close to 7 (28/200) showed the highest structural similarity to the experiment. There seems to be also a lower limit for the proper development of the system trajectory. In case of and , we observed evolution of square spirals and the simulation grid was never fully covered by spirals and waves. At the 2/14 ratio, the grid was already completely covered. At the 3/20 ratio, the circular structures similar to those observed in the experiment were firstly observed. At keeping the ratio, the next increase of the value smoothed only the edges of the structures. Obviously, 7 state levels does not allow to cover the Moore space with 8 neighbouring cells and 14 state levels are not yet sufficient enough to create an appropriate curvature.

The course of the simulation and the type of the limit set (late state) depends on the type of the Garden of Eden as follows: a) the emergence of the ignition point needs one or two neighbouring non-zero points and determines the overall type of the trajectory and b) the number and distribution of ignition centers determines the duration of each trajectory phase. There are as many Gardens of Eden as possible geometrical set-ups of the ignition points, however, the thickness of square and circular waves, shapes of the central objects and the structure of the limit set remain the same. The multilevel cellular automata exhibit less versatile trajectory than the low-level ones21.

In order to compensate our lack of knowledge of the details of the process, we introduced noise into each sub-process in Eqs. (2)-(3) by multiplication of the terms , , and , respectively, by the factor . The noise was generated by a computer random number generator and its level was upper limited to be the uniform noise within certain range. At , and precision of 10 decimal places, levels of noise were 6%, 12%, and 30% for , , and , respectively. When ignition led to circular structures with filled interior, the noise restored spiral formation (Figure 1b–c, iii and Supplementary Material 5) and the simulation was the most similar to the experiment (Figure 2 and Supplementary Material 1). In other experiments at and (data not shown), we observed that the respective noise levels gave the same trajectory. A formation of circular structures which arose as watersheding of regular structures was also the same as observed for the natural-number states.

The proper combination of values and uniform noise range (Figure 2a) generates a sequence of simulated structures as follows:

  • The simulation grid is filled with systems of square dense waves. This has not been observed in the experiment and we interpret it as a lag phase, which precedes the observed formation of circular waves.

  • Circular structures emanate from the centre of square waves.

  • At the certain state, the simulation grid is nearly covered by large circular structures. A few spirals occur and break into a new generation of spirals.

  • Final state is similar to that in the simulation where the states are natural numbers, and (Supplementary Material 4), however, the waves are about 2 grid elements thicker.

  • In the terminology of the multilevel cellular automata21, the uniform (white) noise shifts the system to the basin of attraction similar to the system where only natural numbers are allowed and open a novel system trajectory through the state space.

Let us mention further key similarities between our simulation and actual experiments (Figure 2b–c):

  • The chemical waves do not interfere like material waves but merge.

  • The chemical waves do not maintain the shape (as is the case, e.g., for solitons27).

  • The morphology of interacting patterns (merger patterns) in simulations has comparable traits as in real experiments.

  • Quantitative features of the limit sets, i.e., the last evolutionary stage of the wave-spiral patterns can be set as close as possible to actual experimental data by an appropriate choice of the parameter range.

It should be stressed that the possibility to replace partly the elementary anisotropy by the uniform noise (see Figure 4) — the constructive role of the noise — is a novel observation, which can be viewed as a superposition of starting conditions.

2.2 Decay and re-assembly of transition state complex and theory of mesoscopicity in chemical systems

When noise matters, an observed process is typically mesoscopic. It does follow neither the deterministic rules of the microscopic (or purely mechanical) system nor the statistical-physics tenet of Boltzmann that only the most frequent events are observed17. So, in case of the B–Z reaction, a chemical origin of mesoscopicity should be sought.

The key process at the end of the reaction scheme is the breakage of a crucial chemical bond in the Eyring transition state complex, possibly the decomposition of brominated malonic acid to formic acid and carbon dioxide24, 25. A sketch of molecules involved in the last reaction step is outlined in Figure 3. This process is hidden in the experiment in which we observe only the reduction of the Fe-phenanthroline complex. The rest of, known and unknown, reactions in all their reaction intermediates contributes to the restoration of this activated complex. The activated complex consequently reduces Fe which enables our measurement. The scheme described here still includes breakage of a few individual bonds, possibly in certain orders. Moreover, it is not certain that this redox process is the bottleneck process. It should be also noted that there is a deep controversy about the actual details of the mechanism including the role of molecules involved in the carbon dioxide evolution. Thus the proposed chemical mechanism must be understood more as illustration of the complexity and of the origin of time extent of key reaction steps. The chemical reaction is certainly not occurring at timeless instant as it is assumed in the reaction–diffusion model.

We propose that the state changes of the activated complex are described by processes 2–3 (see Eqs. (3)-(4)). Process 4 (see Eq. (5)) represents the breakage of the complex, which vacates the space for the restoration of a new activated complex in process 1 (see Eq. (2)). The mesoscopicity may be explained from the multitude of processes involved in these elementary steps.

To examine properties of the elementary spatial unit (pixel in the simulation), we analyzed the properties of process 2. We assume that process 2 (see Eq. (3)) reflects all reorganization processes and reactions of all parts of the complex, down to the level of a bond which breaks upon the disruption of the complex. We split the process into sub-processes 2a and 2b. Process 2a is the sum of all processes which occur between the elementary spatial unit and the immediate surroundings. Larger surroundings were ignored, since their influence was assumed to be mediated by neighbouring elementary units. Process 2b represents all processes within this unit. The relatively high noise level of 0.14 and 0.30 for process 2a and 2b, respectively, includes our ignorance of the details and randomness of the process.

Figure 4a and Supplementary Material 6 show the results of the noise-free simulation at , , and different . At the , the process exhibited the same trajectory as at the . At the , the process was in its initial part similar to the course at the , and , later, it tended to give round structures. At the , the course of the simulation was very different from any previous simulation. The final state of this simulation was characterized by spiral doublets — ram’s horns26. At the , the process was already very diffusive. Obviously, the ratio between the constants of processes 2a and 2b determines the geometry of the process. For the course of the simulation similar to the experiment, the suitable ratio is equal to 7. The values of constants and themselves do not significantly influence the course of the simulation.

The time needed for the decay of the active state complex and its restoration in process 4 determines the characteristic timescale of the whole simulation. The ratio determines properties of waves and the size of the spatial unit. This ratio represents an aliquot of energy/entropy needed for the decay of the transition state complex in process 2b. In order to observe waves, process 2b needs to be significantly faster than the transition from the neighbourhood in process 2a and, also, the neighbourhood has to be asymmetric. In other words, noise represents geometric asymmetry (2a noise) and internal kinetic irregularity (2b noise) of the neighbourhood.

2.3 Size of the elementary spatial unit

Figure 5 shows the analysis of wave profiles in the B–Z experiment. The striking similarity of the simulation to real experiment intensity profiles of dense waves (Figures 1, 2 and 5) motivated us to guess the number of molecules per an elementary spatial unit (i.e. the pixel of experimental wave). The number of elementary units per the width of the wave was in the range of 10–20. Since the average width of the wave was 1.5 mm, the elementary unit had 0.07–0.15 mm. The solution above the elementary unit had thickness and volume of 0.5 mm and 10 mm, respectively. Then, the solution contained ca. 3 10 and 10 molecules of water and reactants per elementary unit, respectively. This number lies within the thermodynamic limit. The source of the mesoscopicity has to be sought in the physico–chemical dynamics. It means that only a few energetic/re-organizational events occur within a given time. Since an elementary spatial unit contains roughly 10 molecules of reactants, it is probably unsatisfactory to explain the energetic transition only due to one Eyring-type reaction complex. It is likely that we are dealing with a phase separation which gives rise to structures of an analogous type as, e.g., in liquid crystals27.

3 Discussion

Both wavy and spiral patterns often occur in the Nature. However, despite being usually formed in continuous media, they can also be generated artificially in the threshold-range of cyclic multilevel cellular automata21. The generation of waves and spirals in cellular automata has been initially perceived by the scientific community with scepticism28. The reasons why the formation of spirals by cellular automata is not easily accepted by general chemical community is a) the enforcement of the natural number states (i.e., a coarsening of the configuration space) b) and the need of discrete dynamics on a regular grid which is driven by deterministic albeit not differential rules. In our simulation, the former limitation is completely suppressed by a (near–)continuous set of levels. The latter limitation seems to be only a prejudice. Whereas the formation of Bénard cells in Rayleigh–Bénard convection is driven by temperature gradient29 and the speed gradient generate cells in the Taylor–Couette flow, the Gibbs energy gradient drives the discrete cells formation in the B–Z reaction.

The basic principles of the formation of waves and spirals, i.e. the specific mesoscopic behaviour, are in our simulation following:

  • Proper ratio for a given spatial discretization23. A high and low ratio promotes structural asymmetry and diffussiveness, respectively.

  • Presence of sufficient amount of irregularities due to other processes of comparable rates and leads to a noise, which is specific to each (sub)-process. Only at the appropriate combination of noise in processes 2a and 2b, we observed concentric circular waves followed by spirals. When the specific noise, given by decimal places in state levels, by the pulsed white noise, we may conclude that an increase of the noise in process 2a – i.e. increase of the diffusiveness – leads to the faster formation of spirals without any circular waves. The asymmetry of the surroundings of the elementary unit, as observed in the simulation with natural number state levels, and , is probably a sub-set of the appropriate specific noises modulating the wavefront. We call the whole range of state-level distributions, which attenuate the specific mesoscopic behaviour, a mesoscopic noise.

Main differences between the experiment and the simulation are following:

  1. In the simulation, formation of a few dense waves precedes the emanation of circular waves (never observed in the experiment).

  2. In the experiment, four different wave frequencies are observed (not yet achieved in the simulation).

  3. In the simulation, dense waves/spirals are arising from remnants of broken initial dense waves. In the experiment, we observe a formation of new centers of dense waves, often at place of evolving micro-bubbles. But scenario of broken waves may also be found by careful inspection of the experiment.

  4. In the experiment, after the state of dense waves/spirals, the system still evolves. The waves broaden due to the exhaustion of chemical energy (see a re-started experiment in Supplementary Material 1).

Many differences may be attributed to our ignorance of the proper character of the underlying mesoscopic noise. The introduction of two noise levels into the simulation indicates that the mesoscopic noise originates from two different physico-chemical processes which lead to the geometrical and kinetic asymmetry.

Even if we use the square tessellation and neighbourhood which allows to cover the whole simulation grid — the Moore neighbourhood, the simulation exhibits the formation of circles. A symmetrical ignition rule gives rise to circles and filamentous structures, whereas the simplest asymmetrical ignition leads to interchanges of spirals and waves (with fractal-like structure). The noise provides a proper condition for switching between these two tendencies. The question remains whether the same rules of ignition and noise introduction give rise to the same geometrical structures in any tessellation or whether the neighbourhood of the active complex is, for some reason, due to the physico-chemical, or a more fundamental physical or geometric, principle, always a Moore-like neighbourhood. We can imagine some organization of molecules into a near-regular square lattice which is maintained by chemical dynamics (i.e., by movement of the molecules), similar to that by which the Bénard cell arises in thin layers of viscous liquids.

We assume that the time is discretized by a step needed for the bottleneck process, i.e. for the restoration of the original concentrations of molecules in spatial element after their depletion. It is described by processes 3–4 (see Eqs. (4)-(5)) and represented by 2 time elements. All other processes are related to this time measure. Thus, for the formation of the proper structures, a near regular structure (lattice) where the processes of depletion and repletion take 7 and 2 time elements, respectively, is needed. In the noise-free hotchpotch machine, mixtures of spirals and waves are formed in a relatively broad range of the ratio23 – between and . All these arguments, including the trajectory of any described multilevel cellular automaton, explain the great similarity to the reaction.

Let us conclude that a computer-based reaction-diffusion model of the B–Z reaction is calculated on a regular grid and all but the fastest process are represented by a stepwise increase/decrease. From this point of view, any reaction-diffusion model is a special kind of the cellular automaton26, 30, 21. Apart from a few special cases when an analytic solution of Turing patterns may be used19, the threshold-range cellular automaton offers promising and very realistic tool for description and modelling of the chemical self-organization.

Figure 1: Illustration of the key aspects of the B–Z model. (a) Comparison of the B–Z experiment (i) with the simulation at optimal levels of noise for each process, i.e. 9%, 14% and 30% for process 1, 2a, and 2b, respectively, and and (ii). Images were expanded so as to have comparable widths of traveling waves. (b) Starting points of the simulations (steps 2, 4, 14, 16). The noise-free simulation with natural number states, and in step 2,000 (i), the noise-free simulation with natural number states, and in step 2,596 (ii) and the process described under a in step 18,400 (iii). (c) Final states (limit sets) of processes defined in b. For all processes, and . In the simulation, the black and white corresponds to 0 and , respectively. Original datasets supplied in Supplementary Material 3.
The intensity is colour-coded in the blue colour, the darkest blue represents value 200, the white colour the intensity 0. The unquestionably inspection of the data has to be done using the original data matrices as demonstrated in Fig. 4.
Figure 2: Similarities between the trajectories of the simulation at optimal noise and those of the Belousov–Zhabotinsky experiment. (a) Selected states of the simulation (i) and corresponding images from the course of the experiment (ii). The early stage of the experiment corresponds to the lag phase of the experiment when no waves evolve. For the later stages of the simulation, corresponding structures were found in the experiment. (b) Sections of images which show wave merging. Similar behaviour was not found for material waves and another wavelike structures and indicates that threshold-range cellular automata (i) are proper models for observed phenomena in the Belousov–Zhabotinsky experiment (ii). (c) States of spiral formation. In the simulation (i), the distortion of the dense waves leads to their merging which is the source of formation of spirals. In the experiment (ii), the source of the distortion is often a bubble of carbon dioxide. Otherwise, the formation of spirals is similar to the experiment. For all processes, , , and . In the simulation, the black and white corresponds to 0 and , respectively.
Figure 3: Eyring transition state complex in the Belousov–Zhabotinsky reaction. Selection of crucial molecules involved in oxidation of a final bond is shown. During the reaction, a few bonds are broken. Vibration mode of the bond which is broken first is responsible for the conversion of the activated complex to the product22. In the experiment, the reduction of Fe-phenanthroline complex to the Fe-phenanthroline complex is monitored. The absorbance of the Fe complex is negligible in comparison to the Fe form31.
Figure 4: Ratio of processes 2a and 2b determines the size of elementary unit and the type of state trajectory. (a) Images of later states of simulation at different ratios, , and . At (i), spirals evolve into forms of ram’s horns. To the opposite, (ii) does not form spirals. At (iv), the process is fully diffusive. At (iii), the trajectory is almost identical to the experimental trajectory. (b) The intensity profiles of waves at different ratios. Decrease of the ratio leads to the broadening of waves. The intensity profile of the circular structure is very noisy. In the simulation, the black and white corresponds to 0 and , respectively.
Figure 5: Analysis of traveling waves in the Belousov–Zhabotinsky experiment. (a) Figure with identified wave profiles. (b) Intensity profile of the early circular wave (1) and later dense wave (2). Three colours represent camera channels.

4 Experimental

4.1 Performance of the chemical reaction

For experiments, the oscillating bromate-ferroin-bromomalonic acid reaction (the B–Z reaction recipe32 provided by Dr. Jack Cohen) was chosen. The reaction mixture included 0.34-M sodium bromate, 0.2-M sulphuric acid, 0.057-M sodium bromide (all from Penta), 0.11-M malonic acid (Sigma-Aldrich) as substrates and a redox indicator and 0.12-M 1,10-phenanthroline ferrous complex (Penta) as a catalyst. All reagents were mixed by hand directly in a dish – a circular Petri dish with diameter of 35, 90, 120 and 200 mm and square dish with side of 75 and 30 mm, respectively — in the above-mentioned sequence for 1 min. A special thermostat, which was constructed from a Plexiglas aquarium and a low-temperature circulating water bath-chiller, was keeping a reaction temperature at C.

The chemical waves were recorded by a Nikon D90 camera in the regime of Time lapse (10 s/snapshot) with exposure compensation EV, ISO 320, aperture and shutter speed s. The original 12-bit NEF raw image format was losslessly transformed to the 12-bit PNG format. See Supplementary Materials. All results and protocols are freely available (Supplementary Material 3).

4.2 Model of the Belousov–Zhabotinsky reaction

The model was derived from that provided in the official release15 of Netlogo 5.1.0 and run there. The full model is provided and described in Supplementary Material 7.


This work was financially supported by the Ministry of Education, Youth and Sports of the Czech Republic – projects CENAKVA (No. CZ.1.05/2.1.00/01.0024), CENAKVA II (No. LO1205 under the NPU I program) and The CENAKVA Centre Development (No. CZ.1.05/2.1.00/19.0380). P.J. was supported by the GAČR Grant No. GA14-07983S.


The following files are available free of charge:

  • Supplementary Material 1: Video of the Belousov–Zhabotinsky reaction

    The movies videoBZ.avi and videoBZreshake.avi were made from sections of the image series BZexperiment and BZ reshake as described below. In the movies, every 5th frame of the series is shown with the frequency of 10 frames per second.

    The image series BZexperiment is a dataset from an experiment in a 200-mm Petri dish. The experiment was preceded by a 2-min preparation period, which included mixing of chemicals, shaking on an orbital shaker and lag period under which no travelling waves evolved. The image series BZreshake is a dataset from the re-shaken, i.e. re-started, experiment in a 200-mm Petri dish. Images were captured in 10s intervals.

  • Supplementary Material 2: Video of the noise-free simulation – original excitable medium simulation with spirals

    The movie video_rounding0_noise0.avi was made from the simulation when only natural numbers were allowed, no noise was included, = 28, = 200, and . In the movie, every 30th simulation step is shown with the frequency of 100 frames per second.

  • Supplementary Material 3: Original data

    All original data for experiments and simulations are available upon request to the authors. We apologise for not placing all datasets freely, because some simulations resulted in files of 250 GB size. Complete datasets of a few TB size will be sent upon request on respective hardware media.

  • Supplementary Material 4: Video of the noise-free simulation – original excitable medium simulation with filaments

    The movie video_rounding1_noise0_k3.avi was made from the simulation when one decimal place was allowed, no noise was included, = 28, = 200, and . In this case, similarly as in the case when and it is sufficient to have one non-zero point to start ignition. This further illustrates the conclusion depicted in Figure 4. In the movie, every 30th simulation step is shown with the frequency of 100 frames per second.

  • Supplementary Material 5: Video of the simulation at optimal noise

    The movie video_rounding1_noise_appropriate.avi was made from the simulation at and , 10 decimal places were allowed and optimal levels of noise for each process, i.e., 0.00, 0.12, and 0.30 for process 1, 2a, and 2b, respectively. In the movie, every 100th simulation step is shown with the frequency of 100 frames per second.

  • Supplementary Material 6: Video of the noise-free simulation at different ratios

    Files video_rounding0_noise0_g1_levels2000.avi, video_rounding0_noise0_g280_levels2000.avi, and video_rounding0_noise0_g1000_levels2000.avi show results of simulations at and , , , and as noted in the file name. Simulations exhibit distinct state space trajectories. In the movie, every 30th simulation step is shown with the frequency of 100 frames per second.

  • Supplementary Material 7: Netlogo 5.2 model of the B–Z reaction as an excitable medium


  • 1 Pikovsky, A.; Rosenblum, M.; Kurths, J. Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge University Press: Cambridge, UK, 2001.
  • 2 Tiezzi, E. Steps towards an Evolutionary Physics, WIT Press: Southampton, UK, 2006.
  • 3 Belousov, B. P. Collection of Short Papers on Radiation Medicine, Med. Publ.: Moscow, USSR, 145, 1959.
  • 4 Zhabotinsky, A. M. Periodical Process of Oxidation of Malonic Acid Solution (a Study of the Belousov Reaction Kinetics). Biofizika 1964, 9, 306–311.
  • 5 Biosa, G.; Masia, M.; Marchettini, N.; Rustici, M. A Ternary Nonequilibrium Phase Diagram for a Closed Unstirred Belousov-Zhabotinsky System. J. Chem. Phys. 2005, 308, 7–12.
  • 6 Belmonte, A. L.; Ouyang, Qi; Flesselles, J.-M. Experimental Survey of Spiral Dynamics in the Belousov-Zhabotinsky Reaction. J. Phys. II France, 1997, 7, 1425–1468.
  • 7 Kitahata, H. Spatio-temporal pattern formation in reaction-diffusion systems coupled with convection and geometrical effect, Open Access Thesis and Dissertation, 2006, available at http://hdl.handle.net/2433/64953.
  • 8 Liveri, M. L. T.; Lombardo, R.; Masia, M.; Calvaruso, G.; Rustici, M. Role of the reactor Geometry in the Onset of Transient Chaos in an Unstirred Belousov-Zhabotinsky System. J. Phys. Chem. A 2003, 107, 4834–4837.
  • 9 Rustici, M.; Branca, M.; Brunetti, A.; Caravati, C.; Marchettini, N. Inverse Ruelle-Takens-Newhouse Scenario in a Closed Unstirred Cerium Catalyzed Belousov-Zhabotinsky System. Chem. Phys. Lett. 1998, 293, 145–151
  • 10 Cross, M. C.; Hohenberg, P. S. Pattern Formation Outside of Equilibrium. Rev. Mod. Phys. 1993, 65, 851–1112.
  • 11 Budroni, M. A.; Rustici, M.; Tiezzi, E. On the Origin of Chaos in the Belousov-Zhabotinsky Reaction in Closed and Unstirred Reactors. Math. Model. Nat. Phenom. 2011, 6, 226–242.
  • 12 Holley, J.; Adamatzky, A.; Bull, L.; Costello, B.; Jahan, I. Computational Modalities of Belousov-Zhabotinsky Encapsulated Vesicles. Nano Commun. Netw. 2011, 2, 50–61.
  • 13 von Neumann, J. The General and Logical Theory of Automata; Jeffress,L. A., Ed.; In: Cerebral Mechanisms in Behavior – The Hixon Symposium, 1951, pp. 1–31.
  • 14 Dewdney, A.K. Computer Recreations: The Hodgepodge Machine Makes Waves. Sci. Am. 1988, 225, 104.
  • 15 Wilensky, U. NetLogo B–Z reaction model, available athttp://ccl.northwestern.edu/netlogo/models/B-ZReaction. Center for Connected Learning and Computer-Based Modeling, Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL, 2003.
  • 16 Kinoshita, S.; Iwamoto, M; Tateishi, K.; Suematsu, J. B.; Ueyama, D. Mechanism of spiral formation in heterogeneous discretized excitable media. Phys. Rev. E 2013, 87, 062815.
  • 17 Krapivsky, P. L.; Redner, S., Ben-Naim, E. A Kinetic View of Statistical Physics, Cambridge University Press: Cambridge, UK, 2010.
  • 18 Ortiz de Zárate, J. M.; Sengers, J. V. Hydrodynamics Fluctuations in Fluids and Fluid Mixtures, Elsevier: Oxford, UK, 2006.
  • 19 Turing, A. M. The Chemical Basis of Morphogenesis. Phil. Trans. R. Soc. London 1952, 237, 37–72.
  • 20 Hiltemann, S. Multi-Coloured Cellular Automata, Bachelor thesis, Erasmus Universiteit: Rotterdam, Netherlands, 2008.
  • 21 Wuensche, A. Exploring Discrete Dynamics, Luniver Press, 2011.
  • 22 Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107–115.
  • 23 Štys, D.; Náhlík, T.; Zhyrova, A.; Rychtáriková, R.; Papáček, Š., Císař, P. Model of the Belousov-Zhabotinsky Reaction, in press in LNCS.
  • 24 Field, R. J.; Körös, E.; Noyes, H.M. Oscillations in Chemical Systems. II. Thorough Analysis of Temporal Oscillation in the Bromate-Cerium-Malonic Acid System. J. Am. Chem. Soc. 1972, 94, 8649–8664.
  • 25 Rovinsky, A. B.; Zhabotinsky, A. M. Mechanism and Mathematical Model of the Oscillating Bromate-Ferroin-Bromomalonic Acid Reaction. J. Phys. Chem. 1988, 88, 6081–6084.
  • 26 Fisch, R.; Gravner, J.; Griffeath, D. Threshold-Range Scaling of Excitable Cellular Automata. Stat. Comput. 1991, 1, 23–39.
  • 27 Blasone, M.; Jizba, P.; Vitiello, G. Quantum Field Theory and Its Macroscopic Manifestations, Boson Condensation, Ordered Patterns and Topological Defects, Imperial College Press: Cambridge, UK, 2011.
  • 28 Lauzeral, J.; Halloy, J.; Goldbeter, A. Desynchronization of Cells on the Developmental Path Triggers the Formation of Spiral Waves of cAMP during Dictyostelium aggregation. Proc. Natl. Acad. Sci. USA 1997, 94, 9153–9158.
  • 29 Getling, A. V. Rayleigh-Bénard Convection: Structures and Dynamics. World Scientific, 1998.
  • 30 Kapral, R.; Lawniczak, A.; Masiar, P. Oscillations and Waves in a Reactive Lattice-Gas Automaton. Phys. Rev. Lett. 1991, 66, 2539–2542.
  • 31 Kinoshita, S. Pattern Formations and Oscillatory Phenomena, 1st edition, Elsevier: 2013.
  • 32 Cohen, J. at http://drjackcohen.com/BZ01.html.
  • 33 http://www.biowes.org.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description