- Email: [email protected]
Pattern formation in enzyme inhibition and cooperativity with parallel cellular automata
Parallel Computing 30 (2004) 741–751 www.elsevier.com/locate/parco
Pattern formation in enzyme inhibition and cooperativity with parallel cellular automata Xin-She Yang School of Engineering, University of Wales Swansea, Swansea SA2 8PP, UK Received 10 November 2003; accepted 15 December 2003 Available online 13 May 2004
Abstract Enzyme reactions with inhibition and cooperativity are modelled in terms of a pair of coupled nonlinear reaction–diffusion equations. The governing equations are solved using stochastic cellular automata with local rules derived from the corresponding nonlinear partial differential equations. The parallel cellular automaton is implemented using domain decomposition according to the nature of the locality of its update rules. Numerical simulations show stable 2-D and 3-D pattern formation, and complex patterns have the interesting feature of self-organized criticality. The numerical results of cellular automata are also compared with results obtained from finite difference and finite element methods. 2004 Published by Elsevier B.V. Keywords: Cellular automata; Enzyme reaction; Pattern formation; Reaction–diffusion; Self-organized criticality
1. Introduction Biological reactions involve enzymes which are catalysts to convert proteins and other molecules called substrates into products, while enzymes themselves are not changed by the reaction. Enzymes are highly specific in the sense that one enzyme usually can only catalyze one particular substrate or closely related substrates [7,11]. Due to the complexity of regulation of reactions and feedback in enzyme activities, enzymes can interact, inhibit and cooperate to show complicated
E-mail address: [email protected] (X.-S. Yang). 0167-8191/$ - see front matter 2004 Published by Elsevier B.V. doi:10.1016/j.parco.2003.12.013
742
X.-S. Yang / Parallel Computing 30 (2004) 741–751
behaviors. There have been intensive interesting and important studies such as experimental investigation and observation of pattern formation in catalytic reactions (e.g. [1,12]) in the last decades. However, most of the theoretical models of enzyme kinetics mainly concern the time behavior without inclusion of space coordinates [7,19]. Furthermore, since Turing’s pioneer work on morphogenesis [21], pattern formation has received much attention in theoretical and mathematical biology [15,16,20]. Turing discovered that the primary condition for pattern formation is the interaction of two substance with different diffusion rates [21]. Meinhardt and his colleagues have shown a similar but more profound generation mechanism of pattern formation via local self-enhancement and long range inhibition [15,17]. Many phenomena such as animal skin coating, shell structure and cell morphogenesis have been found to associate with reaction–diffusion pattern formation. Since the generic nonlinearity in the governing equations of reaction–diffusion pattern formation, mathematical analysis is usually very difficult, while their numerical simulations are also difficult because the stiffness of the related matrix and no general numerical scheme is universally robust for the problem. In addition, self-organized criticality has been found in many systems in nature pioneered by Bak and his colleagues [2,3]. Thus, most numerical simulations have mainly concerned the 2-D pattern formation using explicit schemes (e.g., [15,22]). It can be expected that 3-D pattern formation will show more complex and interesting features than its 2-D patterns though it is also much more difficult in both 3-D simulations and the selection of parameters. Enzyme reactions usually involves enzyme inhibition and cooperativity in order to control accurately the direction, extent and speed of the reactions concerned. An enzyme inhibitor is a substrate that inhibits the catalytic action of the enzyme. There are many different types of inhibitors such as irreversible or catalytic poisons, competitive inhibitors and allosteric inhibitors [7,11]. The cooperativity is the effect of one enzyme can bind more than one substrate molecule but the binding of one substrate molecule affects the rate of subsequent binding of other molecules. This is essentially what occurs in hemoglobin and oxygen binding with cooperative behavior [11,19]. The transport process of oxygen (O2 ) and carbon dioxide (CO2 ) in blood show the effect of low level of CO2 concentration on the increase of O2 transport via the allosteric effect of Hþ on oxygen binging to hemoglobin and the enhancement of CO2 transport by low levels of oxygen. Monod–Wyman–Changeux model [19] is the earliest model for cooperative effect in terms of enzyme’s conformation, however, they did not consider the competitive inhibition and there was no spatial diffusion involved [19]. In fact, there are many other reaction and transport processes have the inhibition and cooperativity such as the glycolysis and glycolytic oscillations, cellular homoestasis and intercellular calcium transport (e.g., [9,11]). This paper aims to formulate a new reaction–diffusion model for competitive enzyme inhibition and cooperativity. Then, we non-dimensionalize the model equations and present the typical values of parameters to be used for the computer simulations in the following sections. The stability analysis provides a guideline to simulate different pattern formations in the selection of different parameter combinations. Numerical methods such as stochastic cellular automata and finite element
X.-S. Yang / Parallel Computing 30 (2004) 741–751
743
method will be given in detail. Numerical simulations for 2-D and 3-D pattern formation are then given and the self-organized criticality is tested from the simulated results. We also present some results concerning the spiral waves. The comparison of results using cellular automata and finite element method will be also discussed. 2. Inhibitive-cooperative enzyme reactions Most enzyme reactions involve two steps: first, enzyme E combines with a substrate (S) to form an intermediate complex (C1 ), then the complex breaks down into product (P) releasing E in the process. In addition, another substrate (I) combines with E forms a competitive intermediate complex (C2 ) as an inhibitor, then this complex C2 breaks down into the product (P) and another complex C1 in a cooperative manner as in hemoglobin and oxygen binding [11,19]. This competitive-cooperative enzyme reaction mechanism can be schematically written as k1
k2
S þ E C1 !P þ E; k1
k3
k4
E þ I C2 !C1 þ P; k3
ð1Þ
where k1 , k2 , k3 , k4 are forward reactions, and k1 , k3 are reverse reaction rates. Let S, I, E, C1 , C2 denote the concentrations of [S], [I], [E], [C1 ], [C2 ], respectively. The Michaelis–Menton [18] quasi-steady state theory of enzyme kinetics implies that the intermediate species are in quasi-steady state (i.e. dC1 =dt 0; dC2 =dt 0) and E þ C1 þ C2 ¼ E0 (initial enzyme concentration). Thus, we only need to consider the concentration changes for the two substrates (S and I). After some algebraic simplifications, we have h i E0 ðka k1 k1 ÞS kk1b kk14 I S_ ¼ DS r2 S þ AðS0 SÞ; ð2Þ ka þ S þ ðka þkkb4 =k1 Þ I E0 ðka k3 k3 ka =kb ÞI I_ ¼ DI r2 I ; ka þ S þ ðka þkkb4 =k1 Þ I
ð3Þ
where _ ¼ oto , ka ¼ ðk2 þ k1 Þ=k1 ; kb ¼ ðkb þ k3 Þ=k3 . The term þAðS0 SÞ is the feed rate for adding more substrate S starting at an initial concentration S0 to maintain the process. 2.1. Governing equations in dimensionless form By choosing the typical concentration E0 ka for S and E0 kb =ð1 þ k4 =ka k1 Þ for I so that u ¼ S=E0 ka , v ¼ Ið1 þ kak4 =k1 Þ=E0 kb and using the length scale L and time scale ðk2 þ k1 Þ=k1 k2 , the governing equations in the dimensionless form become ou bu cv ¼ Du r2 u þ að1 uÞ ; ot 1þuþv
ð4Þ
744
X.-S. Yang / Parallel Computing 30 (2004) 741–751
ov dv ¼ Dv r2 v ; ot 1þuþv
ð5Þ
where Aka a¼ ; k2
k1 ka k1 b¼ ; k2
ka k1 k4 ; c¼ k2 ðk1 ka þ k4 Þ
ka k3 d¼ k3 ; k2 kb
ð6Þ
and Du ¼
DS ðk2 þ k1 Þ ; k1 k2 L2
Dv ¼
DI ðk2 þ k1 Þ : k1 k2 L2
ð7Þ
From typical values given in [7], we can estimate that Du ¼ 0:01–0.2, Dv ¼ 0:05– 0.25, and a; c ¼ 0:01–0.5 and b; d ¼ 0:1 to 0.4. 2.2. Stability and pattern generation The reaction diffusion system Eqs. (4) and (5) have thep property offfi diffusion-driffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ven instability. The steady state solutions are u0 ¼ bð 1 þ 4a2 =b2 1Þ=2a and T v0 ¼ 0. Let w ¼ ðu u0 ; v v0 Þ be the small perturbation, then w satisfies ow ¼ Dr2 w þ Mw; ot
ð8Þ
where D ¼ diagðDu ; Dv Þ, and M ¼ ½ð2au0 þ bÞ=ð1 þ u0 Þ;P ðað1 u0 Þ þ cÞ=ð1 þ u0 Þ; 0; d=ð1 þ u0 Þ is a 2 · 2 matrix. Writing w in the form of expðktÞWk where summation is over all wavenumber k, we have jM kI Dk 2 j ¼ 0;
ð9Þ
where I is a 2 · 2 unity matrix. Since RðkÞ > 0 implies that instability, this requires that Dv =Du < d=ð2au0 þ bÞ. The range of unstable wavenumbers between the two roots of k 2 at bifurcation point is given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dDu Dv ð2au0 þ bÞ 1 1 þ 4Du Dv d 2 k ¼ ; ð10Þ Du Dv ð1 þ u0 Þ 2 2
where d ¼ ð2au0 þ bÞ=½dDu Dv ð2au0 þ bÞ . When the unstable criteria is satisfied, any small random perturbation can generate complex patterns.
3. Numerical methods The model Eqs. (4) and (5) for enzyme reactions can be generally written in a system of reaction–diffusion equations in the form /t ¼ Dr2 / þ f ð/Þ;
/ ¼ ½u
T
v ;
ð11Þ
where D ¼ diagðDu ; Dv Þ. The rate f ð/Þ is a general nonlinear function coupling the different components u and v. This nonlinear reaction–diffusion can be solved using
X.-S. Yang / Parallel Computing 30 (2004) 741–751
745
conventional numerical methods such as finite difference (FD) method and finite element (FE) method or the state-of-the-art method such as stochastic cellular automata. As the governing equations are generally nonlinear, some approximations and iterations are needed to obtain the convergent numerical results in the case of FD and FE methods. Alternatively, the nonlinear systems can be simulated by using cellular automata [10,22–24]. The cellular automata for reaction–diffusion systems can be formulated to correspond to the solutions of the related partial differential equations in a qualitative manner or quantitative manner [10]. The macroscopic approach via cellular automata as demonstrated by [22] is usually more efficient than explicit numerical method, and can be more efficient than better numerical technique in many cases. 3.1. Stochastic cellular automaton Using the central difference scheme for space and explicit time-stepping, we can discretize Eq. (11) as " n # n /nþ1 /iþ1;j 2/ni;j þ /ni1;j /ni;jþ1 2/ni;j þ /ni;j1 i;j /i;j ¼D þ þ f ð/ni;j Þ: ð12Þ 2 2 Dt ðDxÞ ðDyÞ By taking Dt ¼ Dx ¼ Dy ¼ 1, this becomes h i n n n n n n n /nþ1 i;j ¼ D ð/iþ1;j þ /i1;j þ /i;jþ1 þ /i;j1 Þ 4/i;j þ /i;j þ f ð/i;j Þ;
ð13Þ
which can be written as a generic form r X /nþ1 ck;l /niþk;jþl þ f ð/ni;j Þ; i;j ¼
ð14Þ
k;l¼r
where summation is over the 2r þ 1 neighborhood. The coefficients ck are determined from the discretization of the governing equations, and for this special case, c1;0 ¼ cþ1;0 ¼ c0;1 ¼ c0;þ1 ¼ D, c0;0 ¼ 1 4D. In fact, this is equivalent to a finitestate cellular automaton [10,22] with a transition rule G ¼ ½gij ði; j ¼ 1; 2; . . . ; N Þ from one state Un ¼ ½/nij ði; j ¼ 1; 2; . . . ; N Þ at time level n to a new state Unþ1 ¼ ½/nþ1 ij ði; j ¼ 1; 2; . . . ; N Þ at time level n þ 1. That is, G : Un 7!Unþ1 ;
gij : /nij 7!/ijnþ1 :
ð15Þ
The basic assumption for the rule inference for stochastic automata is that the function gij ð/nij Þ ¼ /ijnþ1 ; t 6 Cð/nij ; f ; r; N Þ, and C is a function with a range of ½0; 1 , and t is a random variable. At every time step, a random number t is generated for each automaton ði; jÞ. The new state will only be updated if the generated random number is greater than C, otherwise, it remains unchanged. Following a similar derivation for ordinary differential equation [10], the probability function C can be written as Cð/nij ; f ; r; N Þ ¼ Cðef Þ ¼ a þ bef . The coefficients a and b depend on the number of finite states and other factors such as the accuracy of the approximation to the partial differential equations. The parameters of the continuum reaction–diffusion model
746
X.-S. Yang / Parallel Computing 30 (2004) 741–751
shall be calibrated to fit the results given by the stochastic cellular automaton model using a least-squared procedure so as to get the related accurate transition rules. 3.2. Finite difference method The cellular automata method is an explicit method as we can see its equivalence in the discretization. As the governing equations are nonlinear and timestepping is time-consuming, the implicit method usually works well. The nonlinear reaction–diffusion Eq. (11) can be solved numerically using the two-stage scheme [14] nþ1=2 2ð/i;j /ni;j Þ nþ1=2 ¼ D½d2x þ d2y /i;j þ f ð/ni;j Þ; Dt
ð16Þ
where dx ¼ ð/iþ1;j /i1;j Þ=2Dx, dy ¼ ð/i;jþ1 /i;j1 Þ=2Dy and d2x / ¼ ð/iþ1;j 2/i;j þ 2
/i1;j Þ=ðDxÞ . The next stage gives /nþ1=2 from /inþ1 /ni nþ1 nþ1=2 ¼ D½d2x þ d2y ð/i;j þ /ni;j Þ=2 þ f ð/i;j Þ: Dt The convergence is second-order in space and OðDtÞ
ð17Þ 2
in time for any < 1=2.
3.3. Finite element method Finite element method is a very successful classic method with rigorous mathematical PM foundation. For a nonlinear parabolic equation, after substituting / ¼ j¼1 uj ðtÞNj ðx; yÞ into Eq. (11) and some elaborate calculations, the finite element formulation can be written as Mu_ ¼ Ku þ gðuÞ;
ð18Þ R
T Rwhere u ¼ ½ u1 u2 u3 uM . M ¼ X Ni Nj dX is the mass matrix, K ¼ rNi DrNj dX is the stiffness matrix, and gðuÞ is a nonlinear function due to the X nonlinear source term f ð/Þ and natural boundary conditions. Using the implicit time-stepping and Newton–Raphson iteration procedure, the above equation can be solved efficiently.
3.4. Parallel implementation As the cellular automaton for enzyme inhibition and cooperativity involves a large grid and time evolution, thus the parallel implementation is essential for the overall speed of the simulations. The parallel implementation used in this paper is similar to the method of domain decomposition with ghost cells at interacting domain boundaries. Cappuccio et al. gave an excellent implementation and discussion of the ghost cell technique for the cellular automata-based model for coffee percolation [6]. From the implementation point of view, the advantage of cellular automata over the conventional numerical methods is the locality of the updating rules, and its
X.-S. Yang / Parallel Computing 30 (2004) 741–751
747
domain can be easily decomposed into many subdomains for parallelization. With the ghost cell technique and period boundary conditions, the domain of simulation is subdivided among processors with balanced computation load.
4. Simulations and results Using the three numerical methods described in the above section, we can simulate the enzyme inhibition and cooperativity process. We first compare the results using different methods of finite difference, finite element and cellular automata to see the accuracy of the cellular automata. We also present the parallel performance using multiprocessors. Then, we show some simulation results for pattern formations and spiral waves in enzyme activities. 4.1. Comparison of numerical methods To test the accuracy of cellular automata (CA) formulated in the above section, we first compare the results of cellular automata with the conventional finite difference (FD) and finite element (FE) methods for the same system with the same initial configurations and boundary conditions. The left graph in Fig. 1 shows the comparison of the three different methods for a chosen line in the middle of the simulation domain. The value of concentration is plotted versus x coordinate. The dots, solid and dashed correspond to the CA, FD and FE results, respectively. We can see clearly that these three methods give essentially the same results. 4.2. Parallel performance The efficiency of the parallel implementation using domain decomposition and ghost cell technique is estimated by the speed-up of the algorithm over a single processor. The right graph in Fig. 1 shows the parallel performance of speed-up versus
0.8
12
CA FD FE
0.7
10
Speedup
0.6
0.5
8 6
0.4
4
0.3
2
0.2
0 20
30
40
50
60
70
80
2
4
6
8
10
12
14
16
Number of Processors
Fig. 1. Comparison of solutions using cellular automat (CA), finite difference (FD) and finite element (FE) methods (left). Speed-up of the parallel implementation of multi-processors (right).
748
X.-S. Yang / Parallel Computing 30 (2004) 741–751
the number of processors used. We can see that the speed-up of parallel implementation is quite satisfactory even for a small number of processors, and this suggests the domain decomposition for cellular automata is an efficient way for parallelization. 4.3. Pattern formation From the initial random configuration, nonlinear reaction–diffusions of enzyme inhibition and cooperativity can lead to complex patterns. The left picture in Fig. 2 shows an example of the substrate concentrations evolving to the interesting patterns in 2-D configuration with 200 · 200 cells at t ¼ 500. The right picture in Fig. 2 shows the 3-D pattern formation starting from an initial random distribution with 50 · 50 · 50 at t ¼ 250. The patterns and the general characteristics of patterns shown in Fig. 2 only change slowly with time. This implies that patterns and structures formed by local transition rules are relatively stable. Our present results are consistent with the other well-known studies in pattern formation such as [16] for plants and sea shells by using nonlinear reaction–diffusion equations. In addition, our present simulations suggest that patterns arise naturally from the local interactions either through rule-based/agent-based system in terms of relationships among the neighborhood of individuals in stochastic cellular automata or partial differential equations in terms of system variables such as substrate concentrations. 4.4. Spiral waves Spiral and scroll waves may occur in reaction–diffusion systems with simple reaction kinetics in excitable media [4,5,8,13]. Since our present system is also a reaction– diffusion system, we can expect that complex spiral wave patterns may arise under certain boundary and initial conditions. Using the similar technique proposed by Barkley and his colleagues [4,13], we investigate the spiral wave formation in enzyme inhibition and cooperativity. Fig. 3 shows the 2-D simple spiral waves (left) and complex 3-D pattern (right) formed under the conditions of Du ¼ 0:1, Dv ¼ 0:005,
Fig. 2. Pattern formation of substrate concentration from initial random configuration. 2-D patterns of spots start to merge for Du ¼ 0:05, Dv ¼ 0:1 and d ¼ a ¼ 0:02, b ¼ c ¼ 0:05 (the left figure), and 3-D patterns form for Du ¼ 0:15, Dv ¼ 0:05, a ¼ b ¼ 0:1, c ¼ d ¼ 0:05 (the right figure).
X.-S. Yang / Parallel Computing 30 (2004) 741–751
749
Fig. 3. Spiral waves in 2-D (left) and 3-D (right) configurations. The values are Du ¼ 0:1, Dv ¼ 0:005 for 2-D pattern, Du ¼ 0:1, Dv ¼ 0:001 for 3-D case.
d ¼ a ¼ 1, b ¼ c ¼ 0:5 for 2-D pattern, and Du ¼ 0:1, Dv ¼ 0:001, a ¼ b ¼ 1, c ¼ d ¼ 2 for 3-D patterns. The snapshot is taken at t ¼ 500. This suggests that some of the important enzyme functions may be related to wave characteristics [7]. 4.5. Self-organized criticality For a lattice N ¼ 50 50 50, the complex pattern of substrate concentration in enzyme inhibition and cooperativity can be measured by grouping or counting the number of cells for a given value of concentration. The results display a power law in the distribution of the number of cells ðnÞ versus discrete concentration ðuÞ. A least-square fitting of n / uc [2,3], leads to the exponent of c ¼ 1:25 0:04. This implies the self-organized criticality in the complex patterns as generally found in many areas of physical sciences [3]. These results may have important implications to the enzyme reactions in biological systems. Although the detail of enzyme inhibition and cooperativity in some real multi-cellular organisms or biological systems may be not clear, it is likely that the enzyme activity and interactions mainly occur locally so that the same cell or substrate may take different functions at different locations. The blood clotting resulting from a trauma is one good example of localization of enzyme chain reactions. 3-D pattern formation may result in the positional information. Thus reaction–diffusion dominates the process without much contribution from the convection mechanism.
5. Discussion The enzyme reactions with inhibition and cooperativity have been modelled using a nonlinear reaction–diffusion system. A finite state stochastic cellular automaton in 2-D and 3-D configurations has been formulated to simulate the nonlinear system of enzyme pattern formation using the transition rules derived from the partial
750
X.-S. Yang / Parallel Computing 30 (2004) 741–751
differential equations. By using the proper stochastic and transition rules between different states, the finite state automaton can simulate the complexity of enzyme activities. Numerical experiments show that complex patterns can arise from the initial random configuration due to the local transition rules between entities of certain nearest neighborhood while the detail of initial configuration is not important. The stochastic cellular automaton can provide a good alternative way for pattern formation if the local rules are derived from the related reaction–diffusion equations and yet give the same results as numerical solutions of parabolic PDEs. In addition, the power-law relationship between the number of cells and substrate concentrations implies self-organized criticality in the complex patterns. Furthermore, the parallelization of cellular automata using domain decomposition and ghost cell technique is a good way for the speedup of the simulations because of the nature of the locality of the updating rules of the cellular automata.
References [1] K. Asakura, J. Lauterbach, H.H. Rotermund, G. Ertl, Modification of spatiotemporal pattern formation in an excitable medium by continuous variation of its intrinsic parameters co-oxidation on Pt (1 1 0), Phys. Rev. B 50 (1994) 8043–8046. [2] P. Bak, C. Tang, K. Wiesenfeld, Self-organized criticality: an explanation of 1/f noise, Phys. Rev. Lett. 59 (1987) 381–384. [3] P. Bak, How nature works: the science of self-organized criticality, Springer, 1996. [4] D. Barkley, M. Kness, L.S. Tuckerman, Spiral-wave dynamics in a simple model of excitable media: the transition from simple to compound rotation, Phys. Rev. A 42 (1990) 2489–2492. [5] D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D 49 (1991) 61–70. [6] R. Cappuccio, G. Cattaneo, G. Erbacci, U. Jocher, A parallel implementation of a cellular automata based model for coffee percolation, Parallel Computing 27 (2001) 685–717. [7] M. Dixon, E.C. Web, Enzyme, Academic, New York, 1979. [8] M. Gerhardt, H. Schuster, J.J. Tyson, A cellular automaton model of excitable media including curvature and dispersion, Science 247 (1990) 1563–1566. [9] A. Goldbeter, Biochemical Oscillations and Cellular Rhythms: the Molecular Bases of Periodic and Chaotic Behaviour, Cambridge University, Cambridge, 1996. [10] V. Guinot, Modelling using stochastic, finite state cellular automata: rule inference from continuum model, Appl. Math. Model. 26 (2002) 701–714. [11] J. Keener, J. Sneyd, Mathematical Physiology, Springer, 1998. [12] R. Imbihl, Catalysis on microstructured bimetallic surfaces, Chaos 12 (2002) 182–189. [13] D. Margerit, D. Barkley, Selection of twisted scroll waves in three-dimensional excitable media, Phys. Rev. Lett. 86 (2001) 175–178. [14] P.C. Meek, J. Norbury, Twos-stage, two-level finite-difference scheme for nonlinear parabolic equations, IMA J. Number. Anal. 2 (1982) 256–335. [15] H. Meinhardt, Models of Biological Pattern Formation, Academic, London, 1982. [16] H. Meinhardt, The algorithmic beauty of sea shells, Springer, New York, 1995. [17] H. Meinhardt, Biological pattern formation as a complex dynamic phenomenon, Int. J. Bifurcatoin and Chaos 7 (1997) 1–26. [18] L. Michaelis, M.I. Menton, Die Kinetik der Invertinwirkung, Biochem. Z 49 (1913) 333–369. [19] J. Monod, J. Wyman, J.P. Changeux, On the nature of allosteric transition: a plausible model, J. Mol. Biol. 12 (1965) 88–118. [20] J.D. Murry, Mathematical Biology, Springer, New York, 1989.
X.-S. Yang / Parallel Computing 30 (2004) 741–751
751
[21] A. Turing, The chemical basis of morphogenesis, Philos. Trans. B 237 (1952) 37–72. [22] J.R. Weimar, Cellular automata for reaction–diffusion systems, Parallel computing 23 (1997) 1699– 1715. [23] S. Wolfram, Cellular automata as models of complexity, Nature 311 (1984) 419–424. [24] S. Wolfram, Cellular automata and complexity, Addison-Wesley, Reading, Mass, 1994.
Pattern formation in enzyme inhibition and cooperativity with parallel cellular automata Xin-She Yang School of Engineering, University of Wales Swansea, Swansea SA2 8PP, UK Received 10 November 2003; accepted 15 December 2003 Available online 13 May 2004
Abstract Enzyme reactions with inhibition and cooperativity are modelled in terms of a pair of coupled nonlinear reaction–diffusion equations. The governing equations are solved using stochastic cellular automata with local rules derived from the corresponding nonlinear partial differential equations. The parallel cellular automaton is implemented using domain decomposition according to the nature of the locality of its update rules. Numerical simulations show stable 2-D and 3-D pattern formation, and complex patterns have the interesting feature of self-organized criticality. The numerical results of cellular automata are also compared with results obtained from finite difference and finite element methods. 2004 Published by Elsevier B.V. Keywords: Cellular automata; Enzyme reaction; Pattern formation; Reaction–diffusion; Self-organized criticality
1. Introduction Biological reactions involve enzymes which are catalysts to convert proteins and other molecules called substrates into products, while enzymes themselves are not changed by the reaction. Enzymes are highly specific in the sense that one enzyme usually can only catalyze one particular substrate or closely related substrates [7,11]. Due to the complexity of regulation of reactions and feedback in enzyme activities, enzymes can interact, inhibit and cooperate to show complicated
E-mail address: [email protected] (X.-S. Yang). 0167-8191/$ - see front matter 2004 Published by Elsevier B.V. doi:10.1016/j.parco.2003.12.013
742
X.-S. Yang / Parallel Computing 30 (2004) 741–751
behaviors. There have been intensive interesting and important studies such as experimental investigation and observation of pattern formation in catalytic reactions (e.g. [1,12]) in the last decades. However, most of the theoretical models of enzyme kinetics mainly concern the time behavior without inclusion of space coordinates [7,19]. Furthermore, since Turing’s pioneer work on morphogenesis [21], pattern formation has received much attention in theoretical and mathematical biology [15,16,20]. Turing discovered that the primary condition for pattern formation is the interaction of two substance with different diffusion rates [21]. Meinhardt and his colleagues have shown a similar but more profound generation mechanism of pattern formation via local self-enhancement and long range inhibition [15,17]. Many phenomena such as animal skin coating, shell structure and cell morphogenesis have been found to associate with reaction–diffusion pattern formation. Since the generic nonlinearity in the governing equations of reaction–diffusion pattern formation, mathematical analysis is usually very difficult, while their numerical simulations are also difficult because the stiffness of the related matrix and no general numerical scheme is universally robust for the problem. In addition, self-organized criticality has been found in many systems in nature pioneered by Bak and his colleagues [2,3]. Thus, most numerical simulations have mainly concerned the 2-D pattern formation using explicit schemes (e.g., [15,22]). It can be expected that 3-D pattern formation will show more complex and interesting features than its 2-D patterns though it is also much more difficult in both 3-D simulations and the selection of parameters. Enzyme reactions usually involves enzyme inhibition and cooperativity in order to control accurately the direction, extent and speed of the reactions concerned. An enzyme inhibitor is a substrate that inhibits the catalytic action of the enzyme. There are many different types of inhibitors such as irreversible or catalytic poisons, competitive inhibitors and allosteric inhibitors [7,11]. The cooperativity is the effect of one enzyme can bind more than one substrate molecule but the binding of one substrate molecule affects the rate of subsequent binding of other molecules. This is essentially what occurs in hemoglobin and oxygen binding with cooperative behavior [11,19]. The transport process of oxygen (O2 ) and carbon dioxide (CO2 ) in blood show the effect of low level of CO2 concentration on the increase of O2 transport via the allosteric effect of Hþ on oxygen binging to hemoglobin and the enhancement of CO2 transport by low levels of oxygen. Monod–Wyman–Changeux model [19] is the earliest model for cooperative effect in terms of enzyme’s conformation, however, they did not consider the competitive inhibition and there was no spatial diffusion involved [19]. In fact, there are many other reaction and transport processes have the inhibition and cooperativity such as the glycolysis and glycolytic oscillations, cellular homoestasis and intercellular calcium transport (e.g., [9,11]). This paper aims to formulate a new reaction–diffusion model for competitive enzyme inhibition and cooperativity. Then, we non-dimensionalize the model equations and present the typical values of parameters to be used for the computer simulations in the following sections. The stability analysis provides a guideline to simulate different pattern formations in the selection of different parameter combinations. Numerical methods such as stochastic cellular automata and finite element
X.-S. Yang / Parallel Computing 30 (2004) 741–751
743
method will be given in detail. Numerical simulations for 2-D and 3-D pattern formation are then given and the self-organized criticality is tested from the simulated results. We also present some results concerning the spiral waves. The comparison of results using cellular automata and finite element method will be also discussed. 2. Inhibitive-cooperative enzyme reactions Most enzyme reactions involve two steps: first, enzyme E combines with a substrate (S) to form an intermediate complex (C1 ), then the complex breaks down into product (P) releasing E in the process. In addition, another substrate (I) combines with E forms a competitive intermediate complex (C2 ) as an inhibitor, then this complex C2 breaks down into the product (P) and another complex C1 in a cooperative manner as in hemoglobin and oxygen binding [11,19]. This competitive-cooperative enzyme reaction mechanism can be schematically written as k1
k2
S þ E C1 !P þ E; k1
k3
k4
E þ I C2 !C1 þ P; k3
ð1Þ
where k1 , k2 , k3 , k4 are forward reactions, and k1 , k3 are reverse reaction rates. Let S, I, E, C1 , C2 denote the concentrations of [S], [I], [E], [C1 ], [C2 ], respectively. The Michaelis–Menton [18] quasi-steady state theory of enzyme kinetics implies that the intermediate species are in quasi-steady state (i.e. dC1 =dt 0; dC2 =dt 0) and E þ C1 þ C2 ¼ E0 (initial enzyme concentration). Thus, we only need to consider the concentration changes for the two substrates (S and I). After some algebraic simplifications, we have h i E0 ðka k1 k1 ÞS kk1b kk14 I S_ ¼ DS r2 S þ AðS0 SÞ; ð2Þ ka þ S þ ðka þkkb4 =k1 Þ I E0 ðka k3 k3 ka =kb ÞI I_ ¼ DI r2 I ; ka þ S þ ðka þkkb4 =k1 Þ I
ð3Þ
where _ ¼ oto , ka ¼ ðk2 þ k1 Þ=k1 ; kb ¼ ðkb þ k3 Þ=k3 . The term þAðS0 SÞ is the feed rate for adding more substrate S starting at an initial concentration S0 to maintain the process. 2.1. Governing equations in dimensionless form By choosing the typical concentration E0 ka for S and E0 kb =ð1 þ k4 =ka k1 Þ for I so that u ¼ S=E0 ka , v ¼ Ið1 þ kak4 =k1 Þ=E0 kb and using the length scale L and time scale ðk2 þ k1 Þ=k1 k2 , the governing equations in the dimensionless form become ou bu cv ¼ Du r2 u þ að1 uÞ ; ot 1þuþv
ð4Þ
744
X.-S. Yang / Parallel Computing 30 (2004) 741–751
ov dv ¼ Dv r2 v ; ot 1þuþv
ð5Þ
where Aka a¼ ; k2
k1 ka k1 b¼ ; k2
ka k1 k4 ; c¼ k2 ðk1 ka þ k4 Þ
ka k3 d¼ k3 ; k2 kb
ð6Þ
and Du ¼
DS ðk2 þ k1 Þ ; k1 k2 L2
Dv ¼
DI ðk2 þ k1 Þ : k1 k2 L2
ð7Þ
From typical values given in [7], we can estimate that Du ¼ 0:01–0.2, Dv ¼ 0:05– 0.25, and a; c ¼ 0:01–0.5 and b; d ¼ 0:1 to 0.4. 2.2. Stability and pattern generation The reaction diffusion system Eqs. (4) and (5) have thep property offfi diffusion-driffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ven instability. The steady state solutions are u0 ¼ bð 1 þ 4a2 =b2 1Þ=2a and T v0 ¼ 0. Let w ¼ ðu u0 ; v v0 Þ be the small perturbation, then w satisfies ow ¼ Dr2 w þ Mw; ot
ð8Þ
where D ¼ diagðDu ; Dv Þ, and M ¼ ½ð2au0 þ bÞ=ð1 þ u0 Þ;P ðað1 u0 Þ þ cÞ=ð1 þ u0 Þ; 0; d=ð1 þ u0 Þ is a 2 · 2 matrix. Writing w in the form of expðktÞWk where summation is over all wavenumber k, we have jM kI Dk 2 j ¼ 0;
ð9Þ
where I is a 2 · 2 unity matrix. Since RðkÞ > 0 implies that instability, this requires that Dv =Du < d=ð2au0 þ bÞ. The range of unstable wavenumbers between the two roots of k 2 at bifurcation point is given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dDu Dv ð2au0 þ bÞ 1 1 þ 4Du Dv d 2 k ¼ ; ð10Þ Du Dv ð1 þ u0 Þ 2 2
where d ¼ ð2au0 þ bÞ=½dDu Dv ð2au0 þ bÞ . When the unstable criteria is satisfied, any small random perturbation can generate complex patterns.
3. Numerical methods The model Eqs. (4) and (5) for enzyme reactions can be generally written in a system of reaction–diffusion equations in the form /t ¼ Dr2 / þ f ð/Þ;
/ ¼ ½u
T
v ;
ð11Þ
where D ¼ diagðDu ; Dv Þ. The rate f ð/Þ is a general nonlinear function coupling the different components u and v. This nonlinear reaction–diffusion can be solved using
X.-S. Yang / Parallel Computing 30 (2004) 741–751
745
conventional numerical methods such as finite difference (FD) method and finite element (FE) method or the state-of-the-art method such as stochastic cellular automata. As the governing equations are generally nonlinear, some approximations and iterations are needed to obtain the convergent numerical results in the case of FD and FE methods. Alternatively, the nonlinear systems can be simulated by using cellular automata [10,22–24]. The cellular automata for reaction–diffusion systems can be formulated to correspond to the solutions of the related partial differential equations in a qualitative manner or quantitative manner [10]. The macroscopic approach via cellular automata as demonstrated by [22] is usually more efficient than explicit numerical method, and can be more efficient than better numerical technique in many cases. 3.1. Stochastic cellular automaton Using the central difference scheme for space and explicit time-stepping, we can discretize Eq. (11) as " n # n /nþ1 /iþ1;j 2/ni;j þ /ni1;j /ni;jþ1 2/ni;j þ /ni;j1 i;j /i;j ¼D þ þ f ð/ni;j Þ: ð12Þ 2 2 Dt ðDxÞ ðDyÞ By taking Dt ¼ Dx ¼ Dy ¼ 1, this becomes h i n n n n n n n /nþ1 i;j ¼ D ð/iþ1;j þ /i1;j þ /i;jþ1 þ /i;j1 Þ 4/i;j þ /i;j þ f ð/i;j Þ;
ð13Þ
which can be written as a generic form r X /nþ1 ck;l /niþk;jþl þ f ð/ni;j Þ; i;j ¼
ð14Þ
k;l¼r
where summation is over the 2r þ 1 neighborhood. The coefficients ck are determined from the discretization of the governing equations, and for this special case, c1;0 ¼ cþ1;0 ¼ c0;1 ¼ c0;þ1 ¼ D, c0;0 ¼ 1 4D. In fact, this is equivalent to a finitestate cellular automaton [10,22] with a transition rule G ¼ ½gij ði; j ¼ 1; 2; . . . ; N Þ from one state Un ¼ ½/nij ði; j ¼ 1; 2; . . . ; N Þ at time level n to a new state Unþ1 ¼ ½/nþ1 ij ði; j ¼ 1; 2; . . . ; N Þ at time level n þ 1. That is, G : Un 7!Unþ1 ;
gij : /nij 7!/ijnþ1 :
ð15Þ
The basic assumption for the rule inference for stochastic automata is that the function gij ð/nij Þ ¼ /ijnþ1 ; t 6 Cð/nij ; f ; r; N Þ, and C is a function with a range of ½0; 1 , and t is a random variable. At every time step, a random number t is generated for each automaton ði; jÞ. The new state will only be updated if the generated random number is greater than C, otherwise, it remains unchanged. Following a similar derivation for ordinary differential equation [10], the probability function C can be written as Cð/nij ; f ; r; N Þ ¼ Cðef Þ ¼ a þ bef . The coefficients a and b depend on the number of finite states and other factors such as the accuracy of the approximation to the partial differential equations. The parameters of the continuum reaction–diffusion model
746
X.-S. Yang / Parallel Computing 30 (2004) 741–751
shall be calibrated to fit the results given by the stochastic cellular automaton model using a least-squared procedure so as to get the related accurate transition rules. 3.2. Finite difference method The cellular automata method is an explicit method as we can see its equivalence in the discretization. As the governing equations are nonlinear and timestepping is time-consuming, the implicit method usually works well. The nonlinear reaction–diffusion Eq. (11) can be solved numerically using the two-stage scheme [14] nþ1=2 2ð/i;j /ni;j Þ nþ1=2 ¼ D½d2x þ d2y /i;j þ f ð/ni;j Þ; Dt
ð16Þ
where dx ¼ ð/iþ1;j /i1;j Þ=2Dx, dy ¼ ð/i;jþ1 /i;j1 Þ=2Dy and d2x / ¼ ð/iþ1;j 2/i;j þ 2
/i1;j Þ=ðDxÞ . The next stage gives /nþ1=2 from /inþ1 /ni nþ1 nþ1=2 ¼ D½d2x þ d2y ð/i;j þ /ni;j Þ=2 þ f ð/i;j Þ: Dt The convergence is second-order in space and OðDtÞ
ð17Þ 2
in time for any < 1=2.
3.3. Finite element method Finite element method is a very successful classic method with rigorous mathematical PM foundation. For a nonlinear parabolic equation, after substituting / ¼ j¼1 uj ðtÞNj ðx; yÞ into Eq. (11) and some elaborate calculations, the finite element formulation can be written as Mu_ ¼ Ku þ gðuÞ;
ð18Þ R
T Rwhere u ¼ ½ u1 u2 u3 uM . M ¼ X Ni Nj dX is the mass matrix, K ¼ rNi DrNj dX is the stiffness matrix, and gðuÞ is a nonlinear function due to the X nonlinear source term f ð/Þ and natural boundary conditions. Using the implicit time-stepping and Newton–Raphson iteration procedure, the above equation can be solved efficiently.
3.4. Parallel implementation As the cellular automaton for enzyme inhibition and cooperativity involves a large grid and time evolution, thus the parallel implementation is essential for the overall speed of the simulations. The parallel implementation used in this paper is similar to the method of domain decomposition with ghost cells at interacting domain boundaries. Cappuccio et al. gave an excellent implementation and discussion of the ghost cell technique for the cellular automata-based model for coffee percolation [6]. From the implementation point of view, the advantage of cellular automata over the conventional numerical methods is the locality of the updating rules, and its
X.-S. Yang / Parallel Computing 30 (2004) 741–751
747
domain can be easily decomposed into many subdomains for parallelization. With the ghost cell technique and period boundary conditions, the domain of simulation is subdivided among processors with balanced computation load.
4. Simulations and results Using the three numerical methods described in the above section, we can simulate the enzyme inhibition and cooperativity process. We first compare the results using different methods of finite difference, finite element and cellular automata to see the accuracy of the cellular automata. We also present the parallel performance using multiprocessors. Then, we show some simulation results for pattern formations and spiral waves in enzyme activities. 4.1. Comparison of numerical methods To test the accuracy of cellular automata (CA) formulated in the above section, we first compare the results of cellular automata with the conventional finite difference (FD) and finite element (FE) methods for the same system with the same initial configurations and boundary conditions. The left graph in Fig. 1 shows the comparison of the three different methods for a chosen line in the middle of the simulation domain. The value of concentration is plotted versus x coordinate. The dots, solid and dashed correspond to the CA, FD and FE results, respectively. We can see clearly that these three methods give essentially the same results. 4.2. Parallel performance The efficiency of the parallel implementation using domain decomposition and ghost cell technique is estimated by the speed-up of the algorithm over a single processor. The right graph in Fig. 1 shows the parallel performance of speed-up versus
0.8
12
CA FD FE
0.7
10
Speedup
0.6
0.5
8 6
0.4
4
0.3
2
0.2
0 20
30
40
50
60
70
80
2
4
6
8
10
12
14
16
Number of Processors
Fig. 1. Comparison of solutions using cellular automat (CA), finite difference (FD) and finite element (FE) methods (left). Speed-up of the parallel implementation of multi-processors (right).
748
X.-S. Yang / Parallel Computing 30 (2004) 741–751
the number of processors used. We can see that the speed-up of parallel implementation is quite satisfactory even for a small number of processors, and this suggests the domain decomposition for cellular automata is an efficient way for parallelization. 4.3. Pattern formation From the initial random configuration, nonlinear reaction–diffusions of enzyme inhibition and cooperativity can lead to complex patterns. The left picture in Fig. 2 shows an example of the substrate concentrations evolving to the interesting patterns in 2-D configuration with 200 · 200 cells at t ¼ 500. The right picture in Fig. 2 shows the 3-D pattern formation starting from an initial random distribution with 50 · 50 · 50 at t ¼ 250. The patterns and the general characteristics of patterns shown in Fig. 2 only change slowly with time. This implies that patterns and structures formed by local transition rules are relatively stable. Our present results are consistent with the other well-known studies in pattern formation such as [16] for plants and sea shells by using nonlinear reaction–diffusion equations. In addition, our present simulations suggest that patterns arise naturally from the local interactions either through rule-based/agent-based system in terms of relationships among the neighborhood of individuals in stochastic cellular automata or partial differential equations in terms of system variables such as substrate concentrations. 4.4. Spiral waves Spiral and scroll waves may occur in reaction–diffusion systems with simple reaction kinetics in excitable media [4,5,8,13]. Since our present system is also a reaction– diffusion system, we can expect that complex spiral wave patterns may arise under certain boundary and initial conditions. Using the similar technique proposed by Barkley and his colleagues [4,13], we investigate the spiral wave formation in enzyme inhibition and cooperativity. Fig. 3 shows the 2-D simple spiral waves (left) and complex 3-D pattern (right) formed under the conditions of Du ¼ 0:1, Dv ¼ 0:005,
Fig. 2. Pattern formation of substrate concentration from initial random configuration. 2-D patterns of spots start to merge for Du ¼ 0:05, Dv ¼ 0:1 and d ¼ a ¼ 0:02, b ¼ c ¼ 0:05 (the left figure), and 3-D patterns form for Du ¼ 0:15, Dv ¼ 0:05, a ¼ b ¼ 0:1, c ¼ d ¼ 0:05 (the right figure).
X.-S. Yang / Parallel Computing 30 (2004) 741–751
749
Fig. 3. Spiral waves in 2-D (left) and 3-D (right) configurations. The values are Du ¼ 0:1, Dv ¼ 0:005 for 2-D pattern, Du ¼ 0:1, Dv ¼ 0:001 for 3-D case.
d ¼ a ¼ 1, b ¼ c ¼ 0:5 for 2-D pattern, and Du ¼ 0:1, Dv ¼ 0:001, a ¼ b ¼ 1, c ¼ d ¼ 2 for 3-D patterns. The snapshot is taken at t ¼ 500. This suggests that some of the important enzyme functions may be related to wave characteristics [7]. 4.5. Self-organized criticality For a lattice N ¼ 50 50 50, the complex pattern of substrate concentration in enzyme inhibition and cooperativity can be measured by grouping or counting the number of cells for a given value of concentration. The results display a power law in the distribution of the number of cells ðnÞ versus discrete concentration ðuÞ. A least-square fitting of n / uc [2,3], leads to the exponent of c ¼ 1:25 0:04. This implies the self-organized criticality in the complex patterns as generally found in many areas of physical sciences [3]. These results may have important implications to the enzyme reactions in biological systems. Although the detail of enzyme inhibition and cooperativity in some real multi-cellular organisms or biological systems may be not clear, it is likely that the enzyme activity and interactions mainly occur locally so that the same cell or substrate may take different functions at different locations. The blood clotting resulting from a trauma is one good example of localization of enzyme chain reactions. 3-D pattern formation may result in the positional information. Thus reaction–diffusion dominates the process without much contribution from the convection mechanism.
5. Discussion The enzyme reactions with inhibition and cooperativity have been modelled using a nonlinear reaction–diffusion system. A finite state stochastic cellular automaton in 2-D and 3-D configurations has been formulated to simulate the nonlinear system of enzyme pattern formation using the transition rules derived from the partial
750
X.-S. Yang / Parallel Computing 30 (2004) 741–751
differential equations. By using the proper stochastic and transition rules between different states, the finite state automaton can simulate the complexity of enzyme activities. Numerical experiments show that complex patterns can arise from the initial random configuration due to the local transition rules between entities of certain nearest neighborhood while the detail of initial configuration is not important. The stochastic cellular automaton can provide a good alternative way for pattern formation if the local rules are derived from the related reaction–diffusion equations and yet give the same results as numerical solutions of parabolic PDEs. In addition, the power-law relationship between the number of cells and substrate concentrations implies self-organized criticality in the complex patterns. Furthermore, the parallelization of cellular automata using domain decomposition and ghost cell technique is a good way for the speedup of the simulations because of the nature of the locality of the updating rules of the cellular automata.
References [1] K. Asakura, J. Lauterbach, H.H. Rotermund, G. Ertl, Modification of spatiotemporal pattern formation in an excitable medium by continuous variation of its intrinsic parameters co-oxidation on Pt (1 1 0), Phys. Rev. B 50 (1994) 8043–8046. [2] P. Bak, C. Tang, K. Wiesenfeld, Self-organized criticality: an explanation of 1/f noise, Phys. Rev. Lett. 59 (1987) 381–384. [3] P. Bak, How nature works: the science of self-organized criticality, Springer, 1996. [4] D. Barkley, M. Kness, L.S. Tuckerman, Spiral-wave dynamics in a simple model of excitable media: the transition from simple to compound rotation, Phys. Rev. A 42 (1990) 2489–2492. [5] D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D 49 (1991) 61–70. [6] R. Cappuccio, G. Cattaneo, G. Erbacci, U. Jocher, A parallel implementation of a cellular automata based model for coffee percolation, Parallel Computing 27 (2001) 685–717. [7] M. Dixon, E.C. Web, Enzyme, Academic, New York, 1979. [8] M. Gerhardt, H. Schuster, J.J. Tyson, A cellular automaton model of excitable media including curvature and dispersion, Science 247 (1990) 1563–1566. [9] A. Goldbeter, Biochemical Oscillations and Cellular Rhythms: the Molecular Bases of Periodic and Chaotic Behaviour, Cambridge University, Cambridge, 1996. [10] V. Guinot, Modelling using stochastic, finite state cellular automata: rule inference from continuum model, Appl. Math. Model. 26 (2002) 701–714. [11] J. Keener, J. Sneyd, Mathematical Physiology, Springer, 1998. [12] R. Imbihl, Catalysis on microstructured bimetallic surfaces, Chaos 12 (2002) 182–189. [13] D. Margerit, D. Barkley, Selection of twisted scroll waves in three-dimensional excitable media, Phys. Rev. Lett. 86 (2001) 175–178. [14] P.C. Meek, J. Norbury, Twos-stage, two-level finite-difference scheme for nonlinear parabolic equations, IMA J. Number. Anal. 2 (1982) 256–335. [15] H. Meinhardt, Models of Biological Pattern Formation, Academic, London, 1982. [16] H. Meinhardt, The algorithmic beauty of sea shells, Springer, New York, 1995. [17] H. Meinhardt, Biological pattern formation as a complex dynamic phenomenon, Int. J. Bifurcatoin and Chaos 7 (1997) 1–26. [18] L. Michaelis, M.I. Menton, Die Kinetik der Invertinwirkung, Biochem. Z 49 (1913) 333–369. [19] J. Monod, J. Wyman, J.P. Changeux, On the nature of allosteric transition: a plausible model, J. Mol. Biol. 12 (1965) 88–118. [20] J.D. Murry, Mathematical Biology, Springer, New York, 1989.
X.-S. Yang / Parallel Computing 30 (2004) 741–751
751
[21] A. Turing, The chemical basis of morphogenesis, Philos. Trans. B 237 (1952) 37–72. [22] J.R. Weimar, Cellular automata for reaction–diffusion systems, Parallel computing 23 (1997) 1699– 1715. [23] S. Wolfram, Cellular automata as models of complexity, Nature 311 (1984) 419–424. [24] S. Wolfram, Cellular automata and complexity, Addison-Wesley, Reading, Mass, 1994.