key: cord-0039853-lqrssqqg authors: Szakály, Tamás; Lagzi, István; Izsák, Ferenc; Roszol, László; Volford, András title: Stochastic cellular automata modeling of excitable systems date: 2007 journal: nan DOI: 10.2478/s11534-007-0032-7 sha: 5ee2297e23415f9650285ecdbd7085c88f44ae8d doc_id: 39853 cord_uid: lqrssqqg A stochastic cellular automaton is developed for modeling waves in excitable media. A scale of key features of excitation waves can be reproduced in the presented framework such as the shape, the propagation velocity, the curvature effect and spontaneous appearance of target patterns. Some well-understood phenomena such as waves originating from a point source, double spiral waves and waves around some obstacles of various geometries are simulated. We point out that unlike the deterministic approaches, the present model captures the curvature effect and the presence of target patterns without permanent excitation. Spontaneous appearance of patterns, which have been observed in a new experimental system and a chemical lens effect, which has been reported recently can also be easily reproduced. In all cases, the presented model results in a fast computer simulation. Experimental and numerical study of excitable systems has a long tradition. In life sciences, spatiotemporal spread of epidemic diseases or migration of some species can also be successfully simulated by considering the population or the territory as an excitable system. In the study of reaction-diffusion systems, propagation of chemical waves can be put into this framework. For a broad overview on excitable systems and for further literature we refer to [1, 2] . The underlying reaction-diffusion mechanism is usually modelled using a system of partial differential equations. In general, they perform well [1] , but programming a proper numerical solver can be rather complicated and the computational simulations are time consuming (especially, in case of 3 dimensions). The basic framework is the cellular automata (CA) models which have significantly improved in the last decades. Besides the relatively low computational cost of these models, they are popular also due to their generality and simple setup. The development of the CA models can be directed to the qualitative simulation of some phenomena such as diffusion [4] , some reactions [5] and pattern formation [4, 6, 7] , neural signal transport [5] , highway traffic [8, 9] , "pedestrian dynamics" [10] , propagation of forest fires [11] or epidemic diseases [12, 13] . A stochastic CA model has been initiated in [14] , and developed in many aspects. For a broad overview of stochastic CA we refer to [5] . Using another model, the conservation of mass was ensured in [15] and an adaptive modeling was developed in [16] . Some computer programs are available in [17, 18] , which can be used to simulate excitation waves in Belousov-Zhabotinsky type (BZ) reactions or the spread of infectious diseases. In this paper, we present an improved cellular automata model and begin with the simulation of two phenomena: propagation of excitation waves from a single point and formation of double spirals. Moreover, we point out the curvature effect. These present our model as an appropriate tool for the simulation of special phenomena. The aim of our investigation is to simulate the propagation of excitation waves around some obstacles, the formation of target patterns without external excitation and a lens effect for excitation waves. A simple and useful algorithm to simulate excitation waves in BZ reactions is called the "hodgepodge machine" [19] , which can be used with a random initial data resulting in wave phenomena [20, 21] . Note that the transition rules in these models remain deterministic. We combine here the basic ideas of [14] and some details in the systematic study in [6, 7] : we insert a fully probabilistic approach into a detailed model of excitable systems which results in a successful simulation of the above phenomena within a single framework. In the presentation of CA, we follow the approach in [4] and make use of the detailed study in [23] by extending it to include stochastic elements. A cellular automaton consists of a set L of cells, which are usually represented with unit squares or even with their midpoints in a coordinate system. Each cell l ∈ L should have a fixed state at a given time; the set of the states is denoted by Q. In our case, L ⊂ {(c 1 , c 2 ) : c 1 and c 2 are positive integers}. The time evolution of the system is driven by the interaction between the cells and their neighbors l 1 , l 2 , . . . , l n . Basically, in the classical cases, the time step is characterized by a function f : Q × Q n → Q, which gives the state of a given cell in the subsequent time step depending on its present state (first component) and that of its neighbors (last n components). Note that the number n can be different for some cells: it can be reduced near the boundaries or walls depending on the geometry of the investigated media. The cells can represent individual species [24] , small segment of a forest [25] , a micro-volume in a reaction space [7] . In our model, the possible states are given as: In real-life situations, active cells describe small segments where reaction occurs or an infected individual is present. Resting state is the simple interpretation of those segments where the above phenomena are not present but could turn into an active state by the neighbors. In the literature, instead of "active" one frequently uses the term "excited", and accordingly, refractory cells are often called non-excitable ones. In this context, the boundary of the region formed by the active cell is called "front". For a cell l = (c 1 , c 2 ) we define the neighbors as follows: l 2 , l 4 , l 6 and l 8 denote the neighbors which share an edge with l and l 1 , l 3 , l 5 and l 7 denote those which has only one common vertex with l. The time step f : Q d × Q n → Q depends now also on the "past" of a given cell (over d time steps). Corresponding to real life cases, we prescribe in this way that each cell remains in the refractory state over d 1 time steps. Similarly, the number of the time steps while a cell resides in the active state is denoted by d 2 and with these, d := max{d 1 , d 2 }. We suppose that any cell can turn into the active state (loosely, say, it can be infected) by its neighbors in the square lattice with a certain probability: l can infect l 2 , l 4 , l 6 and l 8 (any of its edge neighbors) with probability p, -l can infect l 1 , l 3 , l 5 and l 7 (any of its vertex neighbors) with probability αp, with an appropriate parameter α such that for the front of the active cells, the uniform propagation velocity, or, in mathematical terms, uniform expected propagation length is ensured in all directions. A straightforward computation gives the possible pairs (p, αp): Table 1 Parameters for the propagation probabilities over the edges (p) and vertices (αp) of an active cell. We define the step function depending on the actual state of cells. Let a i denote the state of the cell l at the preceding i-th step (i = 1, 2, . . . , d) and let b j denote the state of its neighbor l j at the actual time step (j = 1, 2, . . . , 8) . • An active cell will turn into · refractory state if it was active during the preceding d 2 time steps. · active state if it was not turning into active state d 2 time steps ago. • A refractory cell will turn into a · resting state if it was refractory during the preceding d 1 time steps. · refractory state if it was not turning into active state d 1 time steps ago. • A resting cell will turn into a · resting state if none of its active neighbors infects it, which occurs with a certain probability. Formally: If a 1 = q 0 then where k 1 is the number of its active neighbors which has a common edge with l: k 1 = #{j ∈ {2, 4, 6, 8} : b j = q 1 } and k 2 is the number of its active neighbors which has a common vertex with l: k 2 = #{j ∈ {1, 3, 5, 7} : b j = q 1 }. · active state if one of its neighbors infects it, which occurs with a certain probability. Formally: If a 1 = q 0 then with k 1 , k 2 defined above. In this section, we discuss the results from a set of test cases motivated by some interesting phenomena studying excitation waves in BZ reactions. First we point out that using the above framework, some simple geometries of excitation waves can be simulated such as a circle shaped single wave and formation of double spirals. Moreover, the variation of propagation speed (curvature effect) is investigated. Afterwards, we present simulations of more unique phenomena which have been also observed in real experiments. For a comparison, we refer to these in each test case. The simulations have been executed with dimensionless quantities. In a concrete situation, one can adjust the parameters using the following considerations: • The choice of the parameters p and αp should lead to an isotropic propagation. For some suitable pairs, see Table 1 . • The quotient of d 1 /d 2 gives the ratio of the active and the refractory zones in the simulations. In the experiments, it is a measurable quantity (depending mainly on the initial concentration setup), and accordingly, one can find an appropriate pair (d 1 , d 2 ) for the simulations. According to the formal transition rules defined in the previous section, the states of the cells were determined in consecutive time steps using a Monte Carlo simulation. Test case 1. First we simulated an excitation wave induced by a point source and the evolution of the active state was investigated, which corresponds to an excitation wave (Figure 1 (a), (b) ). The probabilities p and αp in the transition rules were p = 0.35 and αp = 0.06, respectively, with the same length d 1 = d 2 = 8 of the refractory and the active period, respectively. Using the stochastic approach, the shape of the reaction front corresponds to the observations. Apart from some noise, the reaction front will be radially symmetric after a transient period, such that the propagation velocity is the same in all directions. For the corresponding simulation see Figure 2 . The average of the positions of a planar front depending on the time step are depicted in Figure 3 ( symbols). The case of constant propagation velocity corresponds to the slope of the dotted line which has been fitted to the observed data. In a similar way, we analyzed the evolution of a radially symmetric front of active cells according to Figure 1 (a), (b) . When we computed the average distance from the midpoint, we took into consideration those active cells which have at least two resting neighbors. The curvature effect is demonstrated in the way that a linear curve was fitted to the first 6, 12 and 24 observations, respectively. The slope of these lines approaches the propagation velocity of the planar front. These results are shown in Figure 3 . Note that for these qualitative results, it is essential that we used a stochastic model with the parameter set in Table 1 . Otherwise, the geometry of the cell network will highly influence the shape of the active zone. For a visible comparison, we also simulated the above excitation waves using a deterministic approach, which can be interpreted as a special probabilistic approach with p = αp = 1 (see Figure 1 (c), (d) ). These results correspond to the ones in [23] . We have also simulated the formation of double spiral waves (Figure 4 (a), (b) , (c)). "Spiral waves" have been observed in many experimental setups, and simulated using deterministic [27] and stochastic [28] models. For a solid theoretical framework see [29] . The probabilities p and αp in the time stepping function were chosen to be 0.25 and 0.056, respectively, with the active and refractory period d 1 = 20 and d 2 = 100, respectively. The simulation was initialized by a single active cell. After some steps, one half of the formed active and refractory region was removed and substituted with resting cells. The shape of the active zone is again in good accordance with the results of the real experiments [26] . In this test case, it is essential when half of the formed pattern is removed (see Figure 4 (a) ). If the diameter of the initial half-circle is smaller than d 1 + d 2 , no double spiral evolves. In this case, the active front returning to the midpoint cannot make the cells active here, as they are still in refractory state. Test case 2. In this case, we modeled the propagation of the active zone around obstacles of various shapes ( Figure 5 Initially, a strip perpendicular to the obstacle was filled with active cells which initiate fronts propagating both in forward and backward directions. After some time, the back section was removed. According to the theoretical [29, 30] and experimental [31] investigations, the zone of the active cells are involutes of the obstacle (arcs in case of polygonal obstacles). Some of the theoretical results refer to constant velocity field. Indeed, if the curvature is small enough then its effect can be neglected [32] (see Figure 2 therein). For the corresponding experimental studies, we refer to [33, 34, [36] [37] [38] and for numerical simulations see [39] . Our simulations are in good accordance with these experiments, see Figures 5 (a), (b) and (c) . We also provided a statistical analysis to point out that a reaction front according to Figure 5 (c) consists of arcs centered at the vertices of a slightly bigger square. We considered the cells being in anactive state and we analyzed their distance measured from the vertices. For the data see Figure 6 (a) and (b). Here cells which belong to the flat region of Figure 6 (a) and (b) were considered as the cells in the smaller and the bigger "arcs", respectively. The number of these cells is 9780 and 26000, respectively. We computed the corresponding sample variancesŝ 1 andŝ 2 of the distances for both "arcs". We obtainedŝ 1 = 9.81 for the smaller arc andŝ 2 = 12.91 for the bigger one. Since the length of the active state is 40 time steps in the simulation, this confirms the visual impression that the reaction front consists of arcs. Note that the distance data should not be normally distributed since the outward periphery of the active zone consists of more cells than the inward one. Test case 3. In the framework of the above model, we could reproduce a unique phenomenon which can hardly be simulated using any deterministic approach. In some real experiments, "target patterns" can be detected, in more precise terms: the whole active region consists of concentric circles which propagate outward. These are usually considered as degenerated double spirals. Under some circumstances, such a pattern can appear without any externally forced initial excitation in its center [26] . Using a deterministic approach, only one front of active cells appears, followed by the zone of refractory cells (cf. Figure 1 (d) ). Therefore, in the deterministic simulations, usually a permanent artificial excitation is applied to simulate target patterns [40] . In our approach, however, the excitation can arise in a natural way without any external forcing. Due to the probabilistic transition rules, some cells become active with a time delay. If the length of the refractory period is short compared to the active period, then a refractory cell can become resting at the back of the active zone. This resting cell has a neighboring active one which can infect it again. Usually, it occurs as soon as the first generation of the refractory cells becomes resting again such that the same initial cells will be active. One can thus observe a periodically active centrum without any external forcing. Note that using an external forcing may result in a rather complex behavior: presence of twisted spirals and front reversal has been predicted, which can also be controlled [35] . In course of the simulation, we used again p = 0.35 and αp = 0.06, with the parameters d 1 = 10, d 2 = 3. The results are shown in Figure 7 (a), (b) . Test case 4. We have also simulated the spontaneous appearance of spirals and target patterns [36] using the parameter set p = 0.35, αp = 0.06, d 1 = 10 and d 2 = 3, respectively. In several situations (e.g. in supersaturated or excitable media), physical and chemical processes are initiated by a small random perturbation of the system. Accordingly, in our simulation, a randomly generated new excited cell was placed in the domain in each time step with uniform spatial distribution. In this way, after some time, target patterns, spirals and double spirals can appear. A simulation result is depicted in Figure 8 . Pattern formation phenomena in a precipitation system may produce similar structures. Our new experimental investigations [41] have shown the coexistence of precipitation process and excitability without any external forcing. Spontaneous appearance of travelling waves and spiral formation inside of the precipitation front was thus described for the first time. The dynamics and spatial structure of the observed traveling waves suggest the similar origin of BZ waves and the phenomena in [41] . The result in a real experiment is depicted in Figure 9 . Test case 5. Recently, it has been reported that precipitation fronts and excitation waves may produce a refraction-like behavior [42] . It has been shown that precipitation patterns generated by diffusion front travel through spatial discontinuities similar to their optical counterparts, and obey a Snell-like law. A chemical lens based on a refraction phenomenon, which corresponds to the geometrical wave theory was realized experimentally in a BZ system. For details, see [43] . The propagation velocities are different inside and outside the lens: the wave starting from one given point of the outer "faster" medium is refracted and forms a circular front traveling inwards, collapsing at a given point of the inner "slower" medium. Refraction of chemical waves are explained by simple geometrical arguments [43] , and this effect can be successfully reproduced in numerical simulation using our CA model. Accordingly, we have chosen smaller parameters for the propagation velocities inside of the lens, by multiplying both probabilities with 0.45 such that the ratio of the propagation velocities inside and outside of the lens corresponds to the experimental observation in [43] . This reduces the thickness of the active and refractory zones, whenever the lengths of the refractory and active states were unchanged. The results of the experiment and simulation are shown in Figure 10 (a) and Figure 10 (b), respectively. We may Summarize the work reported in this paper as follows. Combining a standard model of excitable systems with a fully probabilistic approach, we improved the existing models in some aspects: • In our model we can prescribe the length of the active or refractory states for the cells. • We can avoid the rather unrealistic property of many simulations namely, that the shape of the active zone inherits the geometry of the cell network. After a transient time, the shapes of these zones will be smooth apart from some fluctuations. • We combined the simple geometry of the cell network with the isotropic wave propagation property using probabilistic parameters. • We could reproduce the curvature effect on the propagation speed and rotating excitation waves around some obstacles within the cellular automata framework, which allows fast simulation. • We could simulate the formation of target patterns without a permanent external excitation and a lense effect for excitation waves. Varieties of spiral wave behavior: an experimentalist's approach to the theory of excitable media Foundations of Synergetics I. Distributed Active Systems A model for fast computer-simulation of waves in excitable meadia Cellular Automata Modeling of Physical Systems Theory and Applications of Cellular Automata A cellular automaton model of excitable media. 2. Curvature, dispersion, rotating waves and meandering waves A cellular automaton model of excitable media. 3. Fitting the Belousov-Zhabotinskii reaction Statistical physics of vehicular traffic and some related systems A stochastic cellular automaton model for traffic flow with multiple metastable states Simulation of evacuation processes using a bionics-inspired cellular automaton model for pedestrian dynamics A forest fire model and some thoughts on turbulence On the spread of drug-resistant diseases Clustering model for transmission of the SARS virus: application to epidemic control and risk assessment Equivalence of cellular automata to Ising-models and directed percolation Probabilistic cellular automata with conserved quantities Adaptive stochastic cellular automata -Theory A cellular automaton describing the formation of spatially ordered structures in chemical systems Computer recreations: The hodgepodge machine makes waves Spontaneous evolution of spatiotemporal patterns in materials Class of cellular automata for reaction-diffusion systems Phenomenology of excitation in 2-D cellular automata and swarm systems A simple cellular automaton model for influenza A viral infections Formation of space-times structure in a forest-fire model Concentration wave propagation in 2-dimensional liquid-phase self-oscillating system Langevin molecular dynamics of interfaces: Nucleation versus spiral growth On a forest fire model with supposed self-organized criticality A geometrical theory for spiral waves in excitable media Geometric theory of trigger waves -A dynamical system approach Involutes -The geometry of chemical waves rotating in annular membranes Signal transmission in chemical systems -Propagation of chemical waves through capillary tubes Refraction of chemical waves propagating in modified membranes Waves of excitation on nonuniform membrane rings, caustics, and reverse involutes Front Reversals, Wave Traps, and Twisted Spirals in Periodically Forced Oscillatory Media Oscillations, Waves and Chaos in Chemical Kinetics Asymmetric wave propagation between zones loaded with different catalyst concentrations Rotating chemical waves: theory and experiments A fast method to simulate travelling waves in nonhomogeneous chemical or biological media Target patterns in a realistic model of Belousov-Zhabotinsky reaction Pattern Formation and Self-Organization in a Simple Precipitation System Wave Optics of Liesegang Rings Chemical lens This research has been supported by the Hungarian Academy of Sciences (OTKA K60867 and D048673) and theÖveges Research Fellowship of the National Office for Research and Technology. We also ackowledge the support of the national program BSIK of the Dutch government, Netherlands: ICT project BRICKS, Theme MSV1.