Simulation Practice and Theory 7 (2000) 657±689

www.elsevier.nl/locate/simpra

Construction, mathematical description and coding of reactive lattice-gas cellular automaton Jon-Paul Voroney a, Anna T. Lawniczak b,c,* a

Faculty of Medicine, University of Toronto, 1 King's College Circle, Toronto, Ont., Canada M5S 1A8 b The Fields Institute for Research in Mathematical Sciences, 222 College Street, Toronto, Ont., Canada M5T 3J1 c Department of Mathematics and Statistics, Guelph-Waterloo Program for Graduate Work in Physics, University of Guelph, Guelph, Ont., Canada N1G 2W1 Received 18 January 1999; received in revised form 7 October 1999

Abstract We construct a reactive lattice-gas cellular automaton (LGCA) for reaction±di€usion systems and provide extensive discussion of its software coding aspects. The software coding aspects provide rationale for some choices in the construction of LGCA which has been inspired by molecular dynamics. Portability of the C language source code, of the data structures, and of the data formats is discussed and explained. We illustrate the ideas behind the development of LGCA and its code by considering a particular reacting system, the Sel'kov model with immobile complexing species. We demonstrate usefulness of LGCA modelling of reactive systems by presenting various simulation results. We compare these results with the standard numerical simulations of reaction±di€usion equations. We conclude the paper by discussing how LGCA methodology can be applied and extended to other contexts. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Reactive lattice-gas; Gas dynamics; Cellular automata; Turing patterns

1. Introduction The vast literature on lattice-gas cellular automata (LGCA) has been devoted to their construction and application to study physical phenomena, in particular in ¯uid dynamics. [1,4,9,16,24] LGCA was developed with the purpose of creating models

*

Corresponding author. E-mail addresses: [email protected] (J.-P. Voroney), alawnicz@®elds.utoronto.ca; nuptek@ sympatico.ca (A.T. Lawniczak). 0928-4869/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 8 - 4 8 6 9 ( 9 9 ) 0 0 0 2 9 - 4

658

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

that are conceptually simple, well-suited for calculation on computers, and inspired by micro-dynamics of real systems. While many features of LGCA are chosen for reasons of computational feasibility, there has been little presentation in the literature of LGCA codes of how the software coding aspects a€ect LGCA construction. The LGCA codes use bit manipulations and hence avoid errors associated with ¯oating point arithmetic, while still retaining code portability. This makes LCGA di€erent from traditional numerical methods. In this paper, we emphasise software coding aspects, and illustrate their importance in LGCA development. LGCA was introduced to study turbulence without solving Navier±Stokes equations [10,11]. In simple ¯uid dynamics LGCA, particles move among the nodes of a hexagonal lattice through the six edges leading from each node to its nearest neighbours. The edge through which a particle leaves is determined by its velocity. At each node independently of the other nodes, the update rule probabilistically chooses between possible changes in the particles' velocities that model elastic collisions. These collisions obey three laws of physics: conservation of mass, energy, and linear momentum. In LGCA for simple ¯uid ¯ows, particles change direction only through simple two- and three-body collisions. This is sucient to reproduce essential aspects of ¯uid dynamics [10,11]. The success of LGCA for simple ¯uid ¯ow has stimulated rapid development of LGCA for ¯uid phase separation, hydrodynamic interfaces, multiple phase ¯ows, ¯ows through porous media, magnetohydrodynamics, and chemically reacting systems [1,4,9,16,24]. In this paper, our focus is on the development and coding of LGCA for reactive systems. Before presenting the LGCA code, we will ®rst illustrate ideas behind LGCA construction for reactive systems and suggest methods for analyzing the dynamics of the automaton. How to include in LGCA construction only essential features of microscopic dynamics will be explained by considering a particular reacting system, the Sel'kov model with complexing species. For this model we will describe the reactive LGCA mathematical algorithm and outline a derivation of the LGCA macroscopic reaction±di€usion equations. The LGCA constructed in this paper is an extension of a previous work [12,14,15]. It includes the addition of immobile species for studying dynamical e€ects of heterogeneity in the reacting medium. This extension requires a new calculation of the reactive transition probabilities. By comparing the coecients of equal powers of LGCA reaction±di€usion equations with those of the phenomenological description of the Sel'kov model we establish the relationships between LGCA reactive transition probabilities and the rate coecients of the Sel'kov model. From these relationships we calculate the reactive transition probabilitites. We demonstrate the usefulness of LGCA modelling by presenting various simulations. We compare these with simulation results obtained by traditional numerical methods for reaction±di€usion equations and rate laws. Extensive discussion of software coding aspects of LGCA forms the core of the paper. It explains rationales for some choices in LGCA modelling. Additionally, we discuss how the reactive LGCA can be adapted to other reacting or interacting models in biological or social contexts. In this way, the paper not only ®lls a gap in the literature but also provides a blue print for the LGCA modeller.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

659

In view of the importance of portability, the C language code has been implemented to be totally portable across compilers compatible with the ANSI C standard. The only C language libraries used are MATH.H and STDIO.H. A random number generator has been coded in order to avoid the C library standard random number generator and to increase the repetition period of the pseudorandom sequences. The code has been successfully compiled and run on an SGI Indigo II workstation, and an HP A9000/735 workstation, and also on a Toshiba Satellite T2100CS using the Borland compiler (vers. 3.1). The only CPU dependency is the use of a register as a consequence of the use of the declaration REGISTER INT I, J; it is expected that all UNIX workstations and all modern MS Windows 95/98/NT platforms will actually allocate registers to hold those integer variables. LGCA is often referred to as a ®ctitious micro-dynamical world of particles, the ¯ow of which is similar to the ¯ow of sand or water or gas or ball bearings. The notion of a regular lattice is dual to the notion of a tiling. Thus, we have dual interpretations of (1) particles interacting inside cells ± which correspond to tiles ± and moving to other cells through channels associated with the edges of the cells or (2) particles moving on edges of a lattice and interacting at the nodes. Both are possible interpretations of the micro-dynamics of the world in which our particles reside. Which of them is applied depends on the type of physical system that LGCA models. In ¯uid dynamics LGCA momentum, particle numbers and energy are conserved in collisions, while in reactive LGCA neither momentum, nor energy, nor particle numbers are conserved. Thus, the interpretation provided in this paper of particles reacting in cells and exiting through channels is more appropriate for the context of chemical systems than the traditional ¯uid dynamics interpretation of LGCA as a system of particles moving on the edges and colliding at the nodes of the lattice. This interpretation permits the extension of LGCA in other contexts, for example reaction±di€usion systems arising in ecology or epidemiology. In computational implementation of LGCA, the state of a cell is speci®ed by an ordered string of bits, typically an integer; the value of a bit represents the presence (1) or absence (0) of a particle, and the order (or position within the integer)

Fig. 1. The hexagonal cell and the six velocity directions. The arrows indicate the presence of two particles and their velocity directions. The value of the six bit state is 000110 for this con®guration of particles.

660

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

of the bit corresponds to the particle's velocity; see Fig. 1. This could have been implemented also by means of bit ®elds. However, it was decided not the use bit ®elds because several compilers implement them in ways that are processor-speci®c and not ANSI-compliant. While the states of the cells represent the number of particles and their velocities, the local operation step represents collisions between the particles. The operation modelling collisions can be deterministic but is more often probabilistic, and corresponds to local changes in the number and/or velocities of the particles. In LGCA for ¯uid ¯ows, only the velocities of the particles are changed; whereas in reactive LGCA the local operation probabilistically changes both the velocity directions of the particles and the number of particles according to a reaction scheme. LGCA for ¯uid ¯ows conserves the number of particles while reactive LGCA does not [8]; thus, the number of possible changes that can occur in a cell is much greater in LGCA for reactive systems than for ¯uid ¯ows. In an LGCA, the update rule usually consists of two parts: · a local operation which changes the state of a cell according only to its present state · a propagation step, where portions of all cells' states are updated synchronously according to matching portions of their neighbours' states. Hence, the value of the state represents the con®guration of the particles, and the matching portions of the states of the cells represent particles with identical velocity. The local operation corresponds to particle collisions, either elastic or reactive. Propagation corresponds to movement of particles according to their velocities and models their free streaming between collisions. The LGCA dynamics arises from the repetition of the local operation applied synchronously at each node and independently of the other nodes, followed by a propagation step. Particles in LGCA could represent individual molecules, groups of molecules, animals of speci®c species, quanta of energy, or other physical or abstract entities, depending on the application. Similarly, cells could represent sites in a rigid molecular lattice, regions of two- or three-dimensional space, or could have other interpretations which involve position-dependent information. Systems in the biological, physical, and social sciences often have similar abstract structure. Consider systems composed of a large number of components which move about and interact locally in simple, well-understood ways. Examples of such systems are ¯uid ¯ows, chemically reacting systems, gas dynamics, immunological systems, trac ¯ows, spread of epidemics in various populations, information propagation in social and political contexts, and cellular biochemical pathways. Hence, LGCA methodology developed here can be applied to study spatio-temporal dynamics arising from simple interactions in such systems. The interpretation of LGCA as micro-dynamical worlds makes these models distinct from models designed to solve PDEs. A PDE is already an abstract description of a system seen in the physical world, and may be highly phenomenological with little correspondence between the PDE variables and measurable quantities. Additionally, discretisation of space and time is necessary to numerically solve PDEs. Through discretisation an arti®cial microscopic description is introduced that

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

661

has no relationship with the physical microscopic world. On this micro-world, computations are performed. Because computer simulations only approximate solutions of continuous equations, to implement a PDE description on a computer is to abstract it even further from the physical world. LGCA take the unusual step of trying to simulate the physical world directly in the CA format designed for computer implementation. Additionally, since the dynamics of LGCA automatically include e€ects of ¯uctuations arising from probabilistic rules modelling microscopic local interactions, it captures aspects of the spatio-temporal dynamics that cannot be studied with continuous partial di€erential equations. Hence, the presented LGCA provides an alternate methodology to study spatio-temporal dynamics of complex systems. The paper is organised as follows. In Section 3, we construct an LGCA for the Sel'kov model with complexing species. Next, we discuss derivation of LGCA macroscopic reaction±di€usion PDEs and relate them to those of the Sel'kov model in Section 4. The core of the paper is a presentation of the LGCA code in Section 5. In Section 6 we present LGCA simulation results and compare them with the results obtained by solving PDEs and ODEs using traditional numerical methods. Finally, in Section 7 we show how the presented ideas and the constructed LGCA can be easily modi®ed for other applications. 2. Sel'kov model with complexing species In this section we describe an extension of the Sel'kov model [27] which includes complexing species. Our emphasis is on the essential e€ects of the local interactions between molecules of the reactive system which must be incorporated into the LGCA in order to reproduce macroscopic dynamics of the modelled system. Consider a system containing reacting species A; B; X ; Y ; C; P dispersed in some background medium, like a solvent or thin gel. Molecules of species A and B are uniformly and continuously fed in from the boundaries of the medium to maintain their concentrations nearly constant throughout the gel. In reactive collisions with the medium the species A reversibly produces the X species, and B reversibly produces the Y species. The species X and Y are mobile and they di€use through the gel. Their diffusion processes arise as a net e€ect of free streaming and collision processes. Additionally, in reactive collisions the species X and Y reversibly produce the Y species in a cubic autocatalytic reaction. Bound in the gel is a further species C, which is immobile. It can reversibly form an immobile complex product with the Y species which is called P. The described reaction scheme expands the Sel'kov model [26] k1

A X;

…1†

kÿ1

k2

X ‡ 2Y 3Y ; kÿ2

…2†

662

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689 k3

Y B;

…3†

kÿ3

by including a complexing reaction k4

Y ‡ C P:

…4†

kÿ4

The extended Sel'kov model can be described by phenomenological reaction±di€usion equations as follows. Since the concentrations of the species A and B are almost spatially homogeneous and time independent, we can assume that they are equal to constants qA and qB , respectively. Additionally, since the species C and P are immobile and one is converted to the other, the sum of their concentrations qS …x† does not change during the time evolution of the system. Hence, qS …x† ˆ qC …x; t† ‡ qP …x; t†, where qC ˆ qC …x; t† and qP ˆ qP …x; t† are concentrations of C and P species, respectively, and x denotes a vector of space variables and t a time variable. If qX ˆ qX …x; t† and qY ˆ qY …x; t† denote concentrations of the species X and Y, then assuming mass-action kinetics the macroscopic PDEs of the extended Sel'kov model have the form ot qX ˆ k1 qA ÿ kÿ1 qX ÿ k2 qX q2Y ‡ kÿ2 q3Y ‡ DX DqX ; k2 qX q2Y

…5†

kÿ2 q3Y

ÿ ot qY ˆ kÿ3 qB ÿ k3 qY ‡ ÿ k4 qC qY ‡ kÿ4 …qS ÿ qC † ‡ DY DqY ;

…6†

ot qC ˆ ÿ k4 qC qY ‡ kÿ4 …qS ÿ qC †:

…7†

In the above, DX and DY are the di€usion coecients of the species X and Y. The Sel'kov model with a complexing reaction has a complicated bifurcation structure with rich dynamical behaviours, including Turing patterns, steady states, oscillations, bistabilities and excitability [26,30]. Additionally, for heterogeneous distributions of complexing species interesting dynamical e€ects arise from interactions of non-linear local chemical kinetics with heterogeneous transport [30]. Traditionally these dynamical behaviours have been investigated by analysing the phenomenological macroscopic reaction±di€usion Eqs. (5)±(7). These equations are sucient for studying dynamical phenomena on large space and long time scales. However, if the aim is to understand emergence of various spatio-temporal structures from complex molecular interactions and e€ects of molecular ¯uctuations on the dynamics, or evolution of these structures on cellular space and time scales, then the description of the reacting system in terms of macroscopic PDEs is inadequate and it is necessary to look for alternatives. One alternative is a mesoscopic model: reactive LGCA. 3. Construction and mathematical description of reactive LGCA In this section, we construct an LGCA for the Sel'kov model with complexing reaction. The construction is inspired by microscopic dynamics yet is constrained by the realities of large scale simulations. We point out which of the features in the construction of the automaton result from computational considerations and which are

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

663

the essential features arising from an abstraction of the microscopic dynamics of a reacting system. From the presented physical picture in the ®rst paragraph of Section 2 it follows that the essential features of the Sel'kov reaction scheme which must be modelled by LGCA are changes in directions of molecular trajectories resulting from collisions, free streaming between collisions, and changes in number of molecules according to the given reaction scheme. This leads to the following simpli®cation of the microscopic dynamics. 3.1. Space and time discretization We discretise the space and time variables. For simulation and mathematical description purposes we use the term particles for our representation of molecules of the reacting species X ; Y ; C and P. We will perform all computations on bits, and therefore avoid round-o€ errors. A bit represents the presence of a given particle, hence a molecule. Other molecules, like the A and B species or solvents which are present in the real system are not explicitly modelled. Their e€ects on the dynamics of the X ; Y ; C and P species are included through the automaton rules. In LGCA, time moves forward in discrete uniform steps k 2 1; 2; 3; . . . of order of d. We tile the space by regular cells w…r† with centres at a discrete space variable r. Here, we consider two cases: the cells w…r† can be regular hexagons h ˆ h…r† or squares s ˆ s…r†. 1 Thus, depending on the context w…r† can stand for a hexagon h…r† or a square s…r†. We assume that the distance between centers of adjacent cells is 1. The cells w…r† represent some small area or thin volume element of the physical space. Consider the case when the space is tiled by hexagons h…r†. Let ci ˆ hcos…i ÿ 1†p=3; sin…i ÿ 1†p=3i

…8†

be a unit vector, for each i ˆ 1; . . . ; 6: If we connect the center of every hexagon h…r† with the centers of the neighbouring hexagons r ‡ ci , where i ˆ 1; . . . ; 6; then we obtain a hexagonal lattice structure, L ˆ Lh . Fig. 2 shows that the centers of the hexagons h…r† become the nodes of the hexagonal lattice Lh : We denote the space of all the hexagons h…r† by H, i.e., H ˆ fh…r† : r 2 Lh g. When the space is tiled by squares s ˆ s…r†, we denote by ci the Cartesian unit vectors, ci ˆ hcos…i ÿ 1†p=2; sin…i ÿ 1†p=2i

…9†

for i ˆ 1; . . . ; 4. Similarly, if we connect the center of every square s…r† with the centers of the neighbouring squares r ‡ ci , where i ˆ 1; . . . ; 4; then we will obtain a square lattice structure, L ˆ Ls : In this case the centers of the squares s ˆ s…r† become the nodes of the square lattice Ls . We denote the space of all the squares s…r† 1 By tiling the space with hexagons instead of squares we avoid spurious invariants in the dynamical behaviour of the automaton and in simulations.

664

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

Fig. 2. The tiling of hexagonal cells with a dashed lattice superimposed upon it. The lattice depicts the channels through which the particles are transported from cell to cell.

by S, i.e., S ˆ fs…r† : r 2 Ls g, and the space of all the cells w…r† by W; i.e., W ˆ fw…r† : r 2 Lg, where L ˆ Ls or Lh . Unless noted otherwise, we assume that space has periodic boundary conditions, i.e., that the boundary conditions of a modelled physical medium can be approximated through such an assumption. This implies that the lattice L has periodic boundary conditions. However, other boundary conditions like Dirichlet or Neumann boundary conditions can be considered and easily implemented in the constructed LGCA. The particles of each of the reacting species X ; Y ; C, and P are constrained to lie inside the cells w…r†; for r 2 L and move among the cells only along the channels corresponding to edges of the lattice L. Except where confusion might arise, we will identify a cell w…r† with its center r, hence with the node r of the lattice L. For a square lattice Ls , the lattice coordination number m ˆ 4, while for a hexagonal lattice Lh , m ˆ 6. In the next section, m will thus stand for 4 or 6, depending on whether the lattice is square or hexagonal. LGCA can be generalised to other types of tilings of two-dimensional space, or to three-dimensional space, for example two-dimensional triangular tillings or cubes in three-dimensional space, and thus the value of m can be other than 4 or 6. 3.2. Distribution of particles At the initial time, k ˆ 0, particles of the species X ; Y ; C, and P are distributed randomly and independently, according to probabilities given by initial concentrations, among the cells w…r†; r 2 L. This initial distribution is such that there are at most four or six particles of each species per cell w…r†, when the cells w…r† are,

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

665

respectively, squares or hexagons. The LGCA discrete dynamics is constructed in such away that the upper bound of four or six particles per each cell, imposed in the initial distribution, is maintained during the time evolution of the automaton. This condition on the upper bound of the number of particles is imposed for computational purposes. It allows calculations to be performed on bits which results in the elimination of ¯oating point round-o€ errors. Depending on whether the tiling is composed of square or hexagonal cells, we use four or six bits for each mobile species to specify the presence (1) or absence (0) of a particle about to move into one of the adjacent cells. In other words, con®gurations of particles of the mobile species X and Y at each cell w…r† at time k can be described by m-dimensional random vectors gm …r; k† ˆ hgmi …r; k†iiˆ1;...;m , where m 2 fX ; Y g, of Boolean random variables gmi …r; k†. These variables are de®ned as follows: 8 < 1 if at time k there is a particle of species m that exits cell w…r† through the channel determined by ci ; …10† gmi …r; k† ˆ : 0 if not: The Boolean random variables specify not only the number of particles of the mobile species X and Y at each cell w…r† but also the channels through which these particles will leave a cell w…r† at time k. The number of particles of the species X and Y at a cell w…r† at time k are given by random variables gX …r; k† and gY …r; k†, respectively, where gm …r; k† ˆ

m X gmi …r; k†

…11†

iˆ1

for m 2 fX ; Y g. The presence or absence of particles of the immobile species C and P at each cell w…r† at time k is described by discrete random variables gj …r; k†, where j 2 fC; P g: At each cell w…r† the sum of the variables gj …r; k† is time independent, i.e., gC …r; k† ‡ gP …r; k† ˆ gS …r†: The function gS …r† is determined at k ˆ 0, i.e., by the initial distributions gj …r; 0† of particles of the species C and P among cells w…r†. In the simulations presented later for the hexagonal lattice, gS …r† was chosen to be 0 or 3 to minimise storage requirements. 3.3. Reactive transformations of particle numbers At each time k and at each cell w…r†, we model changes in molecular numbers resulting from reactive collisions taking place in a small volume element represented by w…r† during a short time d. These changes are modelled by probabilistic rules which create and/or annihilate particles in the cells w…r†. Hence, in the reactive transformation step, the number of particles of the species X ; Y ; C; and P change at each cell w…r† according to the reaction scheme aX X ‡ aY Y ‡ aC C ‡ aP P ! bX X ‡ bY Y ‡ bC C ‡ bP P ;

…12†

666

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

from an input state a ˆ haX ; aY ; aC ; aP i of numbers as of particles to an output state b ˆ hbX ; bY ; bC ; bP i of numbers bs of particles, where s 2 fX ; Y ; C; P g. Each reactive transition at node r, from a particles at time k to b particles at time k ‡ 1; takes place with probability P …a; b; r; k; k ‡ 1†. We assume that the process of creation and annihilation of particles in the chemical transformation step is completely homogeneous, i.e., is cell and time independent, hence for every r and k, P …a; b; r; k; k ‡ 1† ˆ P …a; b†. These probabilities are some of the parameters of LGCA and their values a€ect the behaviour of the automaton. The probability of no reactive transition taking place is X P …a; b†: …13† P …a; a† ˆ 1 ÿ b6ˆa

The process of creation or annihilation of particles is carried out independently in each cell at each time step and is consistent with the reaction scheme (1)±(4). The creation and annihilation transition probabilities are calculated from the relationships obtained from the comparison of the coecients of equal powers of the phenomenological rate law with the rate law obtained from LGCA. The transition probabilities are expressed in terms of the rate coecients of the Sel'kov reaction scheme in Section 4. 3.4. Advection and channel selections At each time step after completion of the reactive transformation step, particles of the species X and Y residing in a cell w…r† move simultaneously to other neighbouring cells, not necessarily nearest neighbours. They move through the channels determined by the vectors ci for i ˆ 1; . . . ; m, a distance lm given by a positive integer for m 2 fX ; Y g. Hence, if gmi …r; k† ˆ 1, then at time k a particle moves along a channel of the lattice L from a node r to a node r ‡ lm ci , where i 2 f1; . . . ; mg, and m 2 fX ; Y g. The advection of particles is deterministic and such that at most one particle of a given mobile species can leave a cell through each channel. At each time k and at each cell w…r† independently of the other cells two con®gurations of exit channels, one for particles of the species X and the other one for particles of the species Y, are selected randomly and independently. The probabilities with which these channels are selected depend only on the number of particles that require exit channels but not on prior or current channel con®gurations. Let bm be the number of particles of species m in a cell w…r†. Then the probability p…bm † of selecting an exit channel con®guration for the species m, where m 2 fX ; Y g is the inverse of the number of di€erent ways of distributing the bm particles in the m channels, i.e.,  ÿ1 m : …14† p…bm † ˆ bm The described probabilistic selection of con®gurations of channels for each species X and Y at each cell w…r† re¯ects that the physical system is at local equilibrium in a small volume element represented by a cell w…r†:

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

667

3.5. Reactive collisions and di€usive transport The probabilistic selection of con®gurations of exit channels for each species X and Y at each cell w…r† at time k re¯ects changes in trajectories due to the e€ects of a mixing process arising from elastic and reactive collisions in a small volume element of the physical space during a short time d. Hence, the process of creation and annihilation of particles followed by the random selection of channels mimics main features of reactive collisions: changes in the number of particles and their trajectories. When particle numbers are not changed in a cell then the random selection of channels mimics only elastic collisions. The mixing process followed by the advection of the mobile particles models di€usive transport in the physical system. 3.6. LGCA dynamics The LGCA dynamics arise from the repetition of the following three steps: (i) probabilistic creation or annihilation of particles performed synchronously and independently at each cell, (ii) independent probabilistic selection of exit channels through which particles of each species X and Y leave their cells, and (iii) simultaneous deterministic movement of particles to other cells. Since the rules governing changes in particles' numbers and their redistribution among exit channels are probabilistic, LGCA incorporates automatically a€ects of ¯uctuations into its dynamics. This is an important feature of LGCA. Other important aspects of LGCA are the possibility of controlling ¯uctuations, the separation of the reactive time scale from the di€usive one, and the modelling of the di€usion processes of each species X and Y in various ways and independently of each other. For example, the separation of time scales and di€erent di€usion coecients can be achieved by applying the transport step a di€erent number of times for each species per each reactive step. Additionally, this automatically a€ects the ¯uctuations in particle numbers arising from reactive events. The di€erence in di€usion coecients for the species X and Y can also be obtained by allowing their particles to move a distance lX and lY , respectively, in one time step. This results in the ratio of their di€usion coecients proportional to lX : lY . 4. Calculation of transition probabilities for the Sel'kov kinetics with a complexing reaction In order to apply LGCA to study dynamics of a given reaction±di€usion system it is necessary to establish relationships between LGCA parameters and those of the modelled system. We establish these relationships by comparing the coecients of equal powers of the LGCA macroscopic reaction±di€usion equations with those of the phenomenological reaction±di€usion equations of the physical system. Full derivation of the LGCA reaction±di€usion equations can be found in [18]. Here, we only brie¯y outline how LGCA equations can be derived. First, we describe

668

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

the automaton dynamics in terms of LGCA micro-dynamical equations. Since LGCA is a Markov stochastic process, by taking expectations of the micro-dynamical equations we obtain the LGCA discrete Boltzmann equations. These equations provide a discrete deterministic description of the automaton dynamics in contrast with the discrete probabilistic description given in Section 3. In appropriate space and time limits, called di€usive limit, the automaton Boltzmann equations are reduced to the LGCA macroscopic reaction±di€usion equations for the concentrations qs ˆ qs …x; t† of the species s for s 2 fX ; Y ; C; P g, where x is a space variable and t time variable. Hence, the transition is made from the discrete to the continuous description of the automaton. Table 1 summarises the properties of the automaton equations and shows how they emerge from LGCA dynamics of reacting and di€using particles. The LGCA macroscopic reaction±di€usion equations of the Sel'kov kinetics with a complexing reaction are of the form XX oqm sl2m 2 …bm ÿ am †P …a; b†P …a† ‡ r qm ˆ ot 4 a b

…15†

for v 2 fX ; Y g and oqj XX ˆ …bj ÿ aj †P …a; b†P …a† ot a b

…16†

for j 2 fC; P g [18]. In Eq. (15) the parameter s arises from the di€usive limit of the space and time variables and P …a† gives a probability of having aX ; aY ; aC ; aP particles of the species s 2 fX ; Y ; C; P g at a spatial point x at time t, when all the correlations are neglected, hence mÿas Y  m  q …x; t† as  qs …x; t† s 1ÿ : …17† P …a† ˆ as m m s2fX ;Y ;C;P g In the derivation of the above formula we have assumed for simplicity that the maximum number of particles of each of the immobile species at each node is also given by m, however, this assumption is not necessary in the derivation of LGCA Table 1 LGCA micro-dynamical equations variable type: discrete, probabilistic variable values: binary + by taking the expected value of the micro-dynamical equations LGCA Boltzmann equations variable type: discrete, deterministic + by taking the di€usive limit LGCA reaction±di€usion equations variable type: continuous, deterministic

variable values: reals in ‰0; 1Š

variable values: non-negative reals

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

669

equations or calculation of LGCA transition probabilities of the Sel'kov kinetics with a complexing reaction. Let X …bs ÿ as †P …a; b† …18† E…gs ÿ as ja† ˆ b

be the conditional expected change in the number of particles of species s when an input state has a numbers of particles. Assuming that the automaton dynamics is spatially homogeneous, Eqs. (15) and (16) become the LGCA rate laws which by (17) and (18) are expressed in terms of integer powers of qs of at most m mÿac Y  m  qc …t† ac  qc …t† dqs X ˆ E…gs ÿ as ja† 1ÿ …19† ac m m dt a c2fX ;Y ;C;P g for s 2 fX ; Y ; C; P g. For some phenomenological mass action rate laws, for example those with polynomials with powers of concentration higher than m or those of nonpolynomial form, the form (19) of LGCA rate laws is not compatible and transition probabilities P …a; b† cannot be determined. Since the phenomenological rate laws of the Sel'kov kinetics with the complexing reaction deduced from (5)±(7) are polynomials of order at most three, it is possible to determine LGCA transition probabilities P …a; b† so that (19) are the rate laws for the Sel'kov reaction scheme with a complexing species. By comparing the coecients of equal powers of (19) with those of the rate laws deduced from (5)±(7) the following relationships are obtained:   E…gX ÿ aX ja† ˆ r1‡ …a† ÿ r1ÿ …a† ÿ r2‡ …a† ‡ r2ÿ …a† ;   E…gY ÿ aY ja† ˆ r2‡ …a† ÿ r2ÿ …a† ‡ r3‡ …a† ÿ r3ÿ …a† ‡ r4‡ …a† ÿ r4ÿ …a† ;   E…gC ÿ aC ja† ˆ r4‡ …a† ÿ r4ÿ …a† ; E…gP ÿ aP ja† ˆ ÿ E…gC ÿ aC ja†; where r1‡ …a† ˆ k1 qA ;

r1ÿ …a† ˆ kÿ1 aX ; m ; r2‡ …a† ˆ k2 aX aY …aY ÿ 1† mÿ1 r2ÿ …a† ˆ kÿ2 aY …aY ÿ 1†…aY ÿ 2†

m2 ; …m ÿ 1†…m ÿ 2†

…20† …21† …22†

r3‡ …a† ˆ kÿ3 qB ;

r3ÿ …a† ˆ k3 aY ;

…23†

r4‡ …a† ˆ kÿ4 aP ;

r4ÿ …a† ˆ k4 aY aC :

…24†

This comparison speci®es only the conditional expected change E…gs ÿ as ja† in the number of particles of species s when the input number of particles is a. Thus, there are various choices for the transition probabilities P …a; b†. The selection of the transition probabilities P …a; b† is related to an interpretation of the automaton space

670

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

and time scaling [12]. Assuming that (1) the LGCA cells w…r† represent small volume elements of a physical medium, (2) the reactive species X and Y are dilutely dispersed in the reacting medium, and (3) the automaton time step represents a short time, then it is likely that at most one reactive collision happens in the small volume represented by the cell w…r† during the time d corresponding to the automaton update. Under these assumptions the guiding principle in the selection of transition probabilities P …a; b† is to assign non-zero probabilities only to transitions which can be interpreted as single reactions in the reaction scheme (1)±(4), when possible. Applying this rule we can derive the following set of transition probabilities 2 P …a; b†:   P …a; aX ÿ 1; aY ; aC ; aP † ˆ kÿ1 aX ÿ k1 qA ‡ m3 …kÿ2 ÿ k2 †daY ;m daX ;m P …a; aX ‡ 1; aY ; aC ; aP † ˆ k1 qA …1 ÿ daX ;m †  P …a; aX ; aY ÿ 1; aC ; aP † ˆ k3 aY ÿ kÿ3 qB ‡ m3 …k2 ÿ kÿ2 †daX ;m  ‡ kÿ4 aP daC ;0 daY ;m P …a; aX ; aY ‡ 1; aC ; aP † ˆ kÿ3 qB …1 ÿ daY ;m †  P …a; aX ÿ 1; aY ‡ 1; aC ; aP † ˆ k2 aX aY …aY ÿ 1†

m mÿ1  …aY ÿ 1†…aY ÿ 2† 2 m daX ;m …1 ÿ daY ;m † ÿ kÿ2 aY …m ÿ 1†…m ÿ 2†  …aY ÿ 1†…aY ÿ 2† 2 m P …a; aX ‡ 1; aY ÿ 1; aC ; aP † ˆ kÿ2 aY …m ÿ 1†…m ÿ 2†  ÿ k2 m2 aX daY ;m …1 ÿ daX ;m † P …a; aX ; aY ‡ 1; aC ‡ 1; aP ÿ 1† ˆ kÿ4 aP …1 ÿ daY ;m † P …a; aX ; aY ÿ 1; aC ÿ 1; aP ‡ 1† ˆ k4 aY aC ÿ kÿ4 aP daY ;m …1 ÿ daC ;0 † P …a; aX ; aY ; aC ‡ 1; aP ÿ 1† ˆ kÿ4 aP daY ;m daC ;0 with the probabilities of no transition P …a; a† assigned the values X P …a; b†; P …a; a† ˆ 1 ÿ

…25†

b6ˆa

and all other transitions assigned probability zero. In the above transition probabilities, the ®rst eight model single reactive collisions of the chemical micro-dynamics and the ninth models two reactive collisions in order to obtain the correct mass action rate law. This reaction, C ! C ‡ 1, P ! P ÿ 1, is needed because the exclusion principle limits the maximum number of molecules at a node. Reactions that would push the total number of particles of any species past the limit of m must be assigned a probability zero, and reaction nine must be added to compensate for this. With the above transition probabilities, ¯uctuations in particle numbers may be adequately described if the density of mobile particles in the system is not too close to m. 2

The symbol daY ;m represents the delta function, de®ned to be one if aY ˆ m and zero otherwise.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

671

Note that with the given transition probabilities, aC ‡ aP ˆ bC ‡ bP ; hence P …gC …r; t† ‡ gP …r; t† ˆ qS …r†jgC …r; 0† ‡ gP …r; 0† ˆ qS …r†† ˆ 1;

…26†

i.e., the total number of C and P particles is speci®ed at the ®rst time step. This implies that strong correlations exist between gC and gP , which we neglected in formulating P …a† in (17). In order to ensure that the probabilities P …a; b† are not negative, the rate coecients ki must satisfy 1. kÿ1 m ÿ k1 qA P 0, 2. kÿ1 m ÿ k1 qA ÿ m3 …kÿ2 ÿ k2 † P 0, 3. k3 m ÿ kÿ3 qB ÿ kÿ4 m P 0, 4. k3 m ÿ kÿ3 qB ÿ kÿ4 m ÿ m3 …k2 ÿ kÿ2 † P 0, 3 5. k2 m ÿ kÿ2 mÿ2 P 0, 6. kÿ2 m ÿ k2 …m ÿ 1† P 0, 7. k4 m ÿ kÿ4 …m ÿ 1† P 0. In addition, the units of time and space in the macroscopic rate law must be such that all P …a; b† are not greater than one and that P …a; a† are non-negative.

5. Reactive LGCA algorithms and codes The mathematical formulation of LGCA in Section 3 described the stochastic evolution of the ®elds of occupation numbers gm …r; k† ˆ hgmi …r; k†iiˆ1;...;m , where m 2 fX ; Y g. The form of these ®elds was chosen to be Boolean in order to eliminate round-o€ errors, increase speed, and achieve bit democracy in large-scale simulations. Bit democracy means that all bits in a data structure have equal importance; it permits optimal memory utilisation. The Boolean ®elds are stored in a data structure which is described in Section 5.1. This section will provide the link between the computer implementation in terms of bit arithmetic operations on a data structure and the model of particles reacting and di€using on a lattice presented in Section 3. As such, Section 5.1 complements Section 4 where LGCA was related to continuous PDEs and the natural phenomena of molecules reacting and di€using in space. To model chemical reactions, in Section 5.2, we will show how to determine the number of particles of each species in a hexagonal cell 3and how to change this number according to the transition probabilities. We then describe in Section 5.3 an algorithm that distributes the particles into channels in preparation for movement to other hexagonal cells. The movement of the particles among hexagonal cells is implemented eciently through bits shifts within the data structure, Section 5.4. In Section 5.5, we will describe various manipulations of the data used to obtain useful information about the simulation. In the ®nal section, we will provide an over-

3 In this section, we will concentrate on hexagonal tilings of space. From the LGCA code description for a hexagonal tiling, it is easy to deduce how LGCA for square lattices could be implemented.

672

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

view of the structure of the LGCA code. The complete code is available on the WWW at http://www.fields.utoronto.ca/ alawnicz. 5.1. The state data structure Ecient storage of the particles' numbers and velocities is essential to LGCA. In computational implementation, the information contained in the Boolean ®elds about the number and direction of movement of the mobile molecules and the number of immobile molecules is stored in a n-bit word. This section describes a code for which the maximum of the sum of C and P particles at a node is three, and there are at most six particles of the mobile X and Y species at each node; we will see that for these maximums, the word is 16-bits long. The two highest order bits of the word are reserved to indicate the number of C molecules. These two bits are binary representations of the numbers zero through three; thus there is a maximum of three C particles in each hexagonal cell. 4 Similarly, the next two bits represent between zero and three P particles. The next six bits represent the con®guration of X particles in the hexagonal cell. The six lowest order bits are reserved for the particles of the Y species. We label the neighbours of a hexagon from one through six, beginning with the neighbour to the east, as in Fig. 3. The six bits for each of X and Y, from highest order to lowest order, correspond to the presence (1) or absence (0) of a particle about to move into the neighbouring hexagon from highest number to lowest number. The meaning of the bits in the word can be summarised as CC PP X 6X 5X 4X 3X 2X 1 Y 6Y 5Y 4Y 3Y 2Y 1. Thus sixteen bits are needed to determine the state of the hexagonal cell; that is the numbers of C and P particles and the distributions of X and Y particles. Fig. 3 is a graphical representation of a possible state of a hexagonal cell, and gives the value of the 16-bit word for that state. We will need to perform various operations on the 16-bit word to simulate reaction and di€usion of particles; in many computer languages, such as C which is used in the codes presented in this section, the natural data type on which bitwise operations are performed is the integer. The total number of 16-bits in our word is convenient, since an unsigned short integer is 16-bits long. De®ne the number of bits for immobile species as NBI ˆ 2 and the number of bits for mobile species as NBM ˆ 6, since for other simulations a user may wish to change the number of bits needed for immobile species and mobile species storage. For the code in this paper, parameters are denoted with upper case letters. Four macros implement bit masks, each of which has 1's in the positions corresponding to one of the four species, and 0's elsewhere and each of which represents the bits reserved for one of the four species: #define ALLC … ……1  NBI† ÿ 1†  CSHIFT ) /* 11 00 000000 000000 */ #define ALLP … ……1  NBI† ÿ 1†  PSHIFT ) /* 00 11 000000 000000 */ #define ALLX … ……1  NBM† ÿ 1†  XSHIFT ) /* 00 00 111111 000000 */ #define ALLY ……1  NBM† ÿ 1† /* 00 00 000000 111111 */

4 Using more than two bits for immobile species would allow this maximum to increase, for example to six as in the theoretical development of LGCA in the previous sections of this chapter.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

673

Fig. 3. The C particles are shown as large, gray balls. The Y particles are smaller white balls, and the X particles are yet smaller black balls. The mobile X and Y particles have arrows attached to them, pointing to the neighbouring hexagon into which they are to be moved. There is one immobile Y particle: it is attached to a C particle, forming a P complex. The value of the computer word describing this hexagonal cell is 10 01 011111 010010, and there are 2 C, 1 P, 5 X, and 2 Y particles.

CSHIFT is de®ned as (2*NBM+NBI) (or 14), PSHIFT is de®ned as (2*NBM) (or 12), and XSHIFT is de®ned as NBM (or 6). They correspond to the lowest order bits' positions of the C, P, and X species. No shift is needed for the Y species, since the Y bits are already in the lowest order position. The compact structure of the word described above is necessary for large scale simulation; all information about the particles and their directions of motion for a hexagonal cell is stored in a minute amount of data, just four bytes. This is one sixteenth the amount of space needed for a typical ODE solver, where the the amount of space needed would be a 64 bit double for each of the four species. The word structure is bit democratic: each bit has a speci®c purpose and no bit is signi®cantly more important than another. This contrasts with most numerical techniques for solving PDEs, where ¯oating point operations are used, and thus lower-order bits lose any meaning (due to round-o€ error) towards the end of a simulation. For each hexagonal cell in a tiling of NX by NY hexagons, we need one word to describe its modelled particles. This structure is two-dimensional, but is implemented as a linear array for computational eciency; the internal implementation is linear regardless, but we make it explicit because the compiler may not do this transformation as eciently as our explicit instructions. Let asState denote this linear array. In the code we use Hungarian notation for variables: the ®rst one or two letters are lower case and indicate the data type (as for an array of short integers); this is followed by a capitalised variable name (in this case State). The element asState[i*NX+j] where i ˆ 0::NX ÿ 1; j ˆ 0::NY ÿ 1 holds the information about C, P, X, and Y species at the hexagon numbered i*NX+j, corresponding to row i and column j. See Fig. 4 for the relationship between the two-dimensional

674

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

Fig. 4. The array of hexagonal cells is shown on the right, with a dashed lattice superimposed upon it. The lattice depicts the channels through which the particles are transported from cell to cell. The twodimensional graph on the left shows the numbering of the nodes of the linear array of words. If periodic boundary conditions are used, then there are additional channels, for example from node 7 to nodes 0, 4, and 8.

array of hexagonal cells, the linear array of computer words, and the index i*NX+j. If NX is chosen to be a power p of two, then a multiplication by NX in binary becomes a simple bit shift left by p, resulting in a faster code; this is analogous to the ease of multiplying by ten in decimal notation. Many compliers automatically convert an integer multiplication by a power of two to a bit shift, but it is best to check this; the conversion from *NX to  p is a simple search and replace. On the right side of Fig. 4, nodes on the graph of computer words which are connected with a line correspond to hexagons which are adjacent. Rows have alternate connections: either connections E, NE, SE, and W or else connections E, W, NW, and SW. Hexagons have identical adjacent neighbours directly to the east and west, but NOT directly north or south. A particle moving north passes through a sequence of hexagons the centers of which are not on a line pointing directly north. Half of the centers are east of this line, the other half are west. A particle moving north must move alternately east and west as well. Thus a particle moving north in Fig. 4 from hexagon 9 to hexagon 1 can take two paths: 9, 5, 1 or 9, 4, 1. We have seen the manner in which the Boolean ®elds are stored for computational purposes, and provided a graphical image of this data structure. Now we are ready to see how bits in each word can be added or removed, corresponding to chemical reaction, and how bits can be moved around the data structure, corresponding to particles moving among the cells of a tiling of space. 5.2. Reaction Simulating a reactive collision involves counting the number of particles of each species in a hexagonal cell, and changing this number according to a probabilistic reaction scheme. An ecient way to ®nd the number of ones in the Boolean representation of an integer is to realise that passing an integer through a bitwise AND gate with the number that is one less than the integer produces an output that is

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

675

identical to the ®rst integer, except that its lowest order non-zero bit (i.e., the rightmost one) is replaced by a zero bit. Thus, if we pass an integer through an AND gate with the number that is one less than the integer, and repeat this process using the output of the AND gate as the next input, then the number of times this needs to be done to make the output zero is the number of ones in the integer. The following C function, 5 popcnt(i) returns the number of ones in a word i by eciently implementing this algorithm: popcnt(i) register int i; f register int j ˆ 0; while (i) f j++; i & ˆ i ÿ 1; g =  this is equivalent to i ˆ i&…i ÿ 1†;  = return(j); g The ®rst step in changing the number of particles involves counting the number of particles of mobile species X and Y, which are called the populations of X and Y or sPopx and sPopy, respectively, and shifting the number of immobile species C and P into the correct space for interpretation as integer numbers between 0 and 3 sPopy ˆ popcnt(asState[i*NX+j]& ALLY); sPopx ˆ popcnt(asState[i*NX+j]& ALLX); sPopc ˆ (asState[i*NX+j]&ALLC)CSHIFT; sPopp ˆ (asState[i*NX+j]&ALLP)PSHIFT; Before discussing probabilistic changes associated with the step simulating chemical reaction, a short note on random number generation is in order because, with the exception of a small number of computers with specialised architectures, most computers do not have probabilistic mechanisms for generating a random real number of a speci®ed number of decimal places. Instead, to produce the same sequence of random numbers on a variety of computational formats and for reasons of speed and eciency, many uniform random number generators use binary bit operations on a sequence of 32-bit integers. Speed is essential in this subroutine, since LGCA requires the generation of many random numbers. An example of a generator using a fast binary bit operation is the lagged Fibonacci generator [22]. This generator produces a sequence fxn g of 32-bit integers by the recursion xn ˆ xnÿr  xnÿs , where  is a binary operation, typical addition modulo 232 , and r and s are lags. For appropriate choices of r and s the period of the sequence is 232 …2r ÿ 1†. As suggested in [22], we use a subtract-with-borrow generator, xn ˆ …xnÿ22 ÿ xnÿ43 ÿ c† mod …232 ÿ 5†;

…27†

where c is 1 (or 0, respectively) if …xnÿ23 ÿ xnÿ44 ÿ c† was positive (or negative, respectively). With 43 initial seeds and an initial value of c equal to 0 or 1, this

5 This function was created by Shiyi Chen, at Los Alamos National Laboratory. Note the shortness of the code, compared to the length of the text describing the algorithm.

676

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

generator produces a sequence of integers in the inclusive range from 0 to 232 ÿ 6. We de®ne the constant MAXINT to be the number 232 ÿ 5 corresponding to the number one larger than the highest 32-bit integer the generator produces. In the code presented here, the random number lRanNum is generated by a call to the subroutine ranlong_( ). Random numbers are needed for two purposes in the code: to choose between a ®nite number of options with equal probabilities (for example to choose between possible exit channel con®gurations with sPopx number of particles), and to make a change based on some probability (for example in selecting which chemical reaction will occur). To choose between k equally likely options, take the option numbered the (32-bit random integer) mod k, which is uniformly distributed for values of k which are small compared to MAXINT. To make a change based on some probability, ®rst multiply the probabilities by MAXINT, and then compare them to the random integer. This eliminates ¯oating point round-o€ errors, through forcing the probabilities initially to be written as long integers, with the interpretation they are numerators of fractions with the large denominator MAXINT. This has the advantage of distributing the random numbers evenly on ‰0; 1† in intervals of about 10ÿ10 . This even distribution contrasts with random ¯oats: ¯oats near zero must come up less frequently than those near one because the density of ¯oating point numbers decreases away from zero, typically from intervals of about 10ÿ77 near zero to intervals of about 10ÿ7 near one. 6 Returning to the problem of simulating reactive transformations, to change the number of particles according to the reaction probabilities we ®rst determine the number of possible ways in which the particle numbers can change. For the Sel'kov kinetics with a complexing reaction, in Section 4 nine transition probabilities for nine possible changes corresponded to nine reactions; thus the parameter NREA ˆ 9. If we de®ne a population state as the number of particles of each of the species (this corresponds to a ˆ haX ; aY ; aC ; aP i), then there are (NCOMP+1)2 (NBM+1)2 possible population states. NCOMP is the maximum number of immobile species (i.e., three when two bits are reserved for each immobile species). De®ne the population state as sPops: sPops ˆ …sPopc  …NBM ‡ 1†  …NBM ‡ 1†  …NCOMP ‡ 1† ‡ sPopp  …NBM ‡ 1†  …NBM ‡ 1† ‡ sPopx  …NBM ‡ 1† ‡ sPopy†  NREAC; Initially a reaction table is created such that the table entry in the sPops+i position is the probability (times MAXINT) that reaction i or any reaction number lower than i occurs for population state sPops. The entry with number i ˆ 0 is interpreted as no reaction occurring. The ninth reaction occurs with conditional probability 1, if a reaction has occurred and none of the previous eight reactions

6 This is relevant in the simulations for the Sel'kov model because probabilities for some reactive transitions approach 10ÿ7 , and thus when ¯oating points are used for probabilities these reactions can occur with the wrong probability or, worse, not occur at all.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

677

occurred. The number of particles in each hexagonal cell is changed by means of the random number generator ranlong() and the look up reaction table alRtable with the algorithm described by the following C-code, used at each site of the array asState: lRanNum ˆ ranlong_(); if (lRanNum< ˆ alRtable[sPops])f; g else if (lRanNum< ˆ alRtable[sPops+1]) fsPopx ÿ ÿ; g else if (lRanNum< ˆ alRtable[sPops+2]) fsPopx ‡ ‡; g else if (lRanNum< ˆ alRtable[sPops+3]) fsPopy ÿ ÿ; g else if (lRanNum< ˆ alRtable[sPops+4]) fsPopy ‡ ‡; g else if (lRanNum< ˆ alRtable[sPops+5]) fsPopx ÿ ÿ; sPopy ‡ ‡; g else if (lRanNum< ˆ alRtable[sPops+6]) fsPopx ‡ ‡; sPopy ÿ ÿ; g else if (lRanNum< ˆ alRtable[sPops+7]) fsPopy ‡ ‡; sPopc ‡ ‡; sPopp ÿ ÿ; g else if (lRanNum< ˆ alRtable[sPops+8]) fsPopy ÿ ÿ; sPopc ÿ ÿ; sPopp ‡ ‡; g else fsPopc ‡ ‡; sPopp ÿ ÿ; g 5.3. Selection of exiting channels To complete the reactive collision, at each node the particles need to be distributed among the channels in preparation for movement to neighbouring hexagonal cells. This step involves probabilistically choosing a word from all possible word values with the number of particles given at the end of the reaction step. First, we place the number of C and P particles in asState at the correct bit positions: asState‰i  NX ‡ jŠ ˆ …sPopc  CSHIFT† ‡ …sPopp  PSHIFT†; Next we probabilistically choose one from the possible six bit values which have sPopx number of X particles. Since if there are zero (000000) or six (111111) particles, no random number generation is required, we rule out these possibilities ®rst. There are six choose sPopx number of population states, or 6, 15, 20, 15, or 6 depending on whether sPopx is 1, 2, 3, 4, or 5. Since the least common multiple of these numbers is LCM ˆ 60, we create a mixing table asMtable with 5 rows and 60 entries per row. The number of the row corresponds to the number of particles, and values of the entries in each row are evenly distributed between the six bit binary numbers with the row's number of ones. Thus, in the following code we probabilistically choose a short integer sPick to pick between the LCM ˆ 60 entries of the row with sPopx number of particles: if (sPopx>0) f if (sPopx
678

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

asState[i*NX+j] | ˆ asMtable[sPick]  XSHIFT; g else asState[i*NX+j] | ˆ ALLX; g Similarly, a short integer with the number of non-zero bits equal to the number of particles of the Y species must be probabilistically chosen to distribute the particles of the Y species in preparation for movement to neighbouring cells. 5.4. Movement In this section, we will consider movement of the mobile species on the graph representing the hexagonal cells. This can be done eciently by sweeping through the lattice four times: sweeping west and moving particles east; sweeping north and moving particles south; sweeping east and moving particles west; sweeping south and moving particles north. Particles moving NE, NW, SE, and SW will be moved by two of the sweeps. For example, particles moving NW will be moved by the south and east sweeps. The sweep is done in the reverse direction of the movement of the particles so that information is not copied onto information that is still needed. The exception is that one row or column of asState may need to be stored: this is the ®rst row or column that is copied over, and corresponds to particles exiting the lattice (or in the case of periodic boundary conditions, moving to the opposite side of the lattice). The code presented below 7 is for moving the particles north and assumes periodic boundary conditions, so that particles moving north from row asState[0*NX+j], j ˆ 0..NX-1 are placed in row asState[(NY-1)*NX+j], j ˆ 0..NX-1. for(j ˆ 0;j< ˆ NX-1;j++) f sTempSite ˆ (asState[j]) & MASKN; for(i ˆ 0; i< ˆ (NY-2); i++) asState[i*NX+j] ˆ (asState[i*NX+j] & NOTN) | (asState[(i+1)*NX+j] & MASKN); asState[(NY-1)*NX+j] ˆ (asState[(NY-1)*NX+j] & NOTN) | sTempSite; g; For a model with two mobile species, MASKN is de®ned to be …6  26 ‡ 6†10 ˆ 000110 0001102 . See Fig. 5. The parameter NOTMASKN is de®ned as …57  26 ‡ 57†10 ˆ 111001 1110012 . This algorithm does not check if it is moving a 1 or a 0, since such a check would take as much time as actually moving the 1 or 0. It is no more computationally expensive to move two species than it is to move one species. Movement of particles in a southward direction is done in a similar fashion, but 7 The functions in this section and the previous section were partially based on codes provided by Karen Diemer.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

679

Fig. 5. The hexagonal cell with its six channels to nearest neighbours and the direction in which the channels point for nodes on odd and even rows. The mask 000110 000110 is used to pick out the mobile X and Y particles which have a northward velocity component.

sweeping through the lattice in the opposite direction and using di€erent bit masks to pick out the south-moving particles. Other boundary conditions can be implemented by changing the two lines containing sTempSite. For example, for zero Dirichlet boundary conditions, remove line three and remove | sTempSite from the end of line seven. Moving particles east is slightly more complicated, since either one or three directions point eastwardly depending on whether the particle is on an odd or even row. Thus, on odd rows, MASKE3ˆ …35  26 ‡ 35†10 ˆ 100011 1000112 represents particles moving east, while on even rows MASKE1ˆ …1  26 ‡ 1†10 ˆ 000001 0000012 represents particles moving east. The following code describes the movement of particles in the east direction: for(i ˆ 0;i<(NY-1);i++) f sTempSite ˆ (asState[ i*NX + NX-1 ]) & MASKE1; for(j ˆ (NX-1);j>0;j- -) asState[i*NX+j] ˆ (asState[i*NX+j] & NOTE1) | (asState[i*NX+j-1] & MASKE1); asState[i*NX] ˆ (asState[i*NX] & NOTE1 ) | sTempSite; i++; sTempSite ˆ asState[i*NX+NX-1] & MASKE3; for(j ˆ (NX-1);j>0;j- -) asState[i*NX+j] ˆ (asState[i*NX+j] & NOTE3) | (asState[i*NX+j-1] & MASKE3); asState[i*NX] ˆ (asState[i*NX] & NOTE3) | sTempSite; g Moving particles west is done in a manner similar to moving particles east; for details of the code, consult the WWW at

680

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

http : ==www:fields:utoronto:ca= alawnicz 5.5. Analysis of simulations To complete the LGCA code requires that data be stored, though not generally at every time step, and that analysis be performed on the data. To produce the pictures shown in this paper, the asState data structure must be manipulated and stored as a graphical representation of the spatially dependent concentrations, for example as a
gif (graphical image ®le). Care must be taken in interpreting the concentrations: concentrations are particles per unit cell. If the area of the hexagonal cell is one, then q p the distance between centers of adjacent hexagons is 2= 3. Thus, an NX by NY array of hexagons covers an approximately rectangular region of space with size q q p p  2= 3NX by 3=2NY. Also, many standards for images store data from left to right and bottom to top; care must be taken so that e€ects due to possible re¯ections of the images are correctly interpreted. 5.6. Code overview The main module of the LGCA code has the following architecture: create or load the mixing table create or load the reaction table initialise the lattice or load an initial asState initialise the random number generator iterate (over time): iterate (over space): perform reaction and distribution perform transport for some times, analyse data, store results, or produce a picture In addition, a program is needed to create the reaction and scattering tables. Complete versions of the LGCA code used to create the simulations displayed in this paper are available on the WWW at http : ==www:fields:utoronto:ca= alawnicz 6. Comparison with other methods LGCA for reactive systems have several advantages over PDE models that are solved using traditional numerical methods; however, care must be taken in interpreting LGCA results. LGCA do not approximate solutions to PDEs in the traditional sense: standard schemes, such as ®nite di€erence methods derived using Euler's approximation for time and space derivatives, discretise space, time, and the operation of taking a derivative, whereas a LGCA simulation is a realisation of discrete, probabilistic micro-dynamics which can be approximated by a solution

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

681

to a continuous, deterministic PDE. As such, LGCA simulates a system's dynamics without solving a PDE. Nevertheless, in many situations we expect a LGCA simulation and a numerical approximation to the solution of a PDE to be very similar. In this section, we present simulations which show that LGCA realisations and numerical approximations to solutions of PDEs or LGCA realisations averaged over the whole lattice and numerical approximations to solutions of ODEs can be very close. However, we will see that LGCA dynamics show e€ects not exhibited by PDE models, due to the incorporation of ¯uctuations arising from the probabilistic description. Hence, LGCA can capture aspects of spatio-temporal dynamics of physical systems that cannot be captured by PDEs or ODEs models. We will end this section by summarising the di€erences and similarities between the LGCA method for simulating a reactive system and methods based on numerically approximating a solution to PDEs or ODEs. 6.1. LGCA simulations of excitable dynamics of the Sel'kov model Many spatially extended systems either quickly become spatially homogeneous despite an initially heterogeneous state, or remain in a spatially homogeneous state if they begin that way. Examples are the mixing of miscible ¯uids, bulk oscillations in well-stirred chemical systems, or systems of yeast cells which synchronise the oscillations of their glycolytic cycles. It is interesting to examine similarities and di€erences in the dynamics of LGCA and ODEs models of such systems, after transient spatial e€ects have damped out. As an example, we have chosen to simulate dynamics in an excitable regime of the Sel'kov model, without the complexing species. Fig. 6 shows

Fig. 6. Trajectories of the Sel'kov model in the excitable regime. The solid lines show solutions to a difference equation obtained using the Runge±Kutta approximation for four di€erent initial conditions. The dashed lines show the spatially averaged concentrations of a LGCA simulation on a hexagonal lattice.

682

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

the results of simulations with several initial conditions which compare LGCA realisations and simulations of the ODE rate equations using the Runge±Kutta approximation for the time derivative [2]. In Fig. 6, the parameters are k1 a ˆ 0:00001677, kÿ1 ˆ 0:00001, k2 ˆ kÿ2 ˆ 0:002, k3 ˆ 0:0001, and kÿ3 b ˆ 0:0000004652. There is no complexing species. For the LGCA simulations, the time step is chosen to be unity, the spatial averages of the concentrations of the X and Y species are plotted every 50 time steps, and the lattice size is 128 by 128. For the Runge±Kutta ODE simulation, the time step is 0.1 and points are plotted every 500 time steps. The trajectories for the LGCA simulations are very similar to those of the Runge±Kutta simulations. If the number of reactive steps is increased, say to three reactive steps per di€usive step, then spontaneous nucleation can arise on a square lattice, as in Fig. 7. On square lattice simulations, particles initially on ``red'' squares move to ``black'' squares, and thus never interact with particles initially on ``black'' squares. Thus, due to the checkercoard parity, only one of the two subsimulations is shown in Fig. 7. Note that when the ratio of reactive steps to di€usive steps is greater than one, the assumption that di€usion smoothes out local ¯uctuations after each reactive step is certainly very poor. To derive the appropriate PDE for these micro-dynamics, transition probabilities equivalent to performing three consecutive reactive steps should be calculated, and these used to determine the macroscopic rate law. However, spontaneous nucleation is a feature that is not found in solutions to deterministic reaction±di€usion equations regardless of the rate law used. 6.2. LGCA simulations of Turing pattern in the Sel'kov model Turing patterns belong to a special class of stable, spatially inhomogeneous solutions of reaction±di€usion equations [29]. They are self-organising phenomena resulting from spontaneous symmetry breaking associated with bifurcations from spatially homogeneous steady states. In an initially homogeneous medium, symmetry breaking instabilities arise from the interaction of local chemical kinetics with transport processes like di€usion [20] or ¯ows [25], and can lead to spatial patterns. A small peak in the concentration of one species does not always decay; hence, the system does not approach a spatially homogeneous state. Instead, a local, shortrange activation of one chemical species and a local, short-range inhibition of another chemical species is encouraged by the faster transport of the inhibitor species away from the peak of the activator species. The locations of the activated peaks (or inhibited valleys) undergo a spatial annealing process until a steady or periodic pattern develops. In chemically reacting systems the faster transport of the inhibitor species can result from the complexing reaction of the activator with an immobile species. Turing patterns in two-dimensions can be an hexagonal array of dots or a series of stripes. Many observed Turing structures are one wavelength thick [3,23]. Thus, their mathematical analysis can often be reduced to two spatial dimensions. Turing patterns are characterised by an intrinsic wavelength which depends on the reaction rates, input rates, and transport properties of the reactants. This wavelength does not depend on the geometry of the system. Parameter regimes where Turing

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

683

Fig. 7. The concentration of the Y species shows spontaneous nucleation on a square lattice. Three reactive steps per di€usive step are necessary to see this e€ect. The boundary conditions are periodic. The time between snapshots increases by 1000 time steps of the LGCA.

684

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

patterns form in the reaction±di€usion systems can be predicted by a linear stability analysis. For the Sel'kov model with complexing reaction this analysis has been performed in [30]. While Turing patterns have been investigated numerically and theoretically for many years, it was not until 1989 that these spatio-temporal structures were observed experimentally [3,23]. The ®rst experiment to show Turing pattern formation involved chlorite and iodide fed in at one end of the gel reactor and malonic acid fed in at the other. It has been suggested by Lengyel and Epstein [19] that in the experiments [3,23] the Turing patterns result from the e€ect of complexing species ± starch and the gel ± on the transport of the activator iodine. The original purpose of the gel was to damp out convection, and the purpose of the starch was to visualise optically triiodide. Yet starch and the gel have important transport e€ects; for some chemical concentrations the immobilization has created a sucient e€ective di€erence in di€usion coecients for Turing pattern to emerge. While such complexing reactions had not received much study in the inorganic chemical systems where experimental evidence of Turing patterns was long sought, membrane bound species and large, nearly immobile polymers are commonly seen to play various roles in cellular systems and they could play roles similar to those of the gel and starch of the above experiments. The LGCA simulations presented here show formation of Turing structures in the Sel'kov model with immobile complexing species. These patterns are stable to molecular ¯uctuations. The Sel'kov model was derived from analyzing the glycolytic pathway in order to study oscillations [26]. It describes qualitatively the basic dynamical properties of one step production of ATP which occurs naturally and universally in cell cytoplasm. This highly stylised reaction is typical of many enzyme catalysed reactions which occur in the highly heterogeneous medium of biological cells [12]. The presented LGCA simulations of the Turing pattern can be of importance to studies in molecular cell biology and biochemistry of the cell [13]. Speci®cally, since LGCA include ¯uctuations which have a physical interpretation: the ¯uctuations correspond to local increases and decreases in the number of molecules due to discrete chemical reactions. Since Turing's original work [29] attempted to describe the formation of morphological structures, it is pleasing that these new results on the importance of immobile complexing species may be particularly applicable to cellular systems. Fig. 8 shows LGCA simulation performed in the regime of parameters such that for the Sel'kov model without the complexing reaction the stable solution is a limit cycle but when the complexing reaction is included and the initial concentration of the species C is greater than a critical value of 1.797 the Turing pattern can form. These parameters were found by performing linear stability analysis [30] on the system of Eqs. (5)±(7) and they are the same as those used in [30] to simulate Turing pattern. In [30] the e€ects of heterogeneity in the initial distribution of the complexing species on Turing pattern formation were investigated and various simulations of PDEs (5)±(7) were performed using ®nite di€erence methods. In Fig. 8, complexing species is distributed on the lower two-thirds of the spatial domain. Boundary conditions are periodic, to approximate an in®nite domain.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

685

Fig. 8. Formation of waves and Turing patterns in a LGCA simulation. Two snapshots of the concentration sum, qY ‡ qP are shown for times separated by 2000 time steps. For reference the period of oscillations is approximately 16600 time steps. The rate coecients are k1 a ˆ 0:001785, kÿ1 ˆ 0:00035, k2 ˆ kÿ2 ˆ k3 ˆ 0:0035, kÿ3 b ˆ 0:0000306, k4 ˆ 0:000525 and kÿ4 ˆ 0:00035, and the lattice is composed of 1024 by 1184 hexagonal cells.

Turing patterns are expected to form on the lower two-thirds, and waves to propagate on the upper third [30]. One e€ect observed in the LGCA simulation, which has not been seen in the simulations in [30] is the apparently spontaneous formation of waves away from the Turing dots, i.e., in the center of the region with no complexing species [18,30]. In the region where there is no complexing species the stable attractor is a limit cycle and due to ¯uctuations in LGCA simulations, small regions in the automaton become phase shifted with respect to each other. This phase decoherence means that these small regions can be further ahead on the trajectory of the limit cycle than near by regions. In these regions the particles that make up the local peaks in the concentration of the Y species di€use outwards, and in doing so cause rapid chemical activation of the surrounding regions. The result is the spreading of the peaks. This explains the apparently spontaneous formation of waves [18,31]. In ®nite di€erence simulations such e€ects due to internal ¯uctuations are not observed [30]. In ®nite di€erence simulations [30], waves originated only from the Turing dots on the boundary between the region with complexing species and the region without complexing species were observed. These type of waves are also observed in LGCA simulations. 6.3. Summary Finite di€erence schemes are often faster than LGCA on conventional, serial machines because they converge on small lattices with large time steps. However, the architecture of serial machines was designed for algorithms such as ®nite di€erence schemes. On machines designed for CA simulation [21,28], LGCA codes can in many

686

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

ways outperform standard simulation methods. Finite di€erence techniques also produce smoother solutions which are closer to analytical solutions of PDEs. The purpose of lattice-gas cellular automata for reactive systems is not to replace standard techniques for solving PDEs; instead LGCA naturally includes e€ects due to ¯uctuations arising from the probabilistic particle description. In PDEs, including ¯uctuations often involves appending noise terms which often lack an interpretation from the modelled system's micro-dynamics. Due to the compact storage of data, LGCA can be performed on larger lattices allowing for more detailed simulations. Due to the bit level description of lattice-gas cellular automata, LGCA can also model systems where traditional schemes fail due to convergence problems associated with round-o€ error. It is straightforward to relate LGCA to PDEs, and bifurcation theory performed on PDEs can be used to predict and classify much of the behaviour expected in LGCA simulations. The results of ®nite di€erence schemes are smooth, while LGCA pictures, even after averaging, show evidence of the noisy nature of molecular interactions. The detail of the pictures is similar; however, LGCA pictures can be magni®ed (averaged over fewer than 16 nodes) and more detail (as well as more noise) could be observed. No extra detail can be observed in ®nite di€erence pictures. In performing simulations, it should be noted that the rate coecients in the LGCA simulation have a maximum due to the requirement that transition probabilities be non-negative and not greater than one, and thus the rate law of a reaction±di€usion equation may have to be altered by rescaling time in order to satisfy these requirements. The di€usion coecient is determined by the spacing between lattice nodes and the update rules, for example the number of lattice sites the particles move between velocity randomisations. In the ®nite di€erence scheme, the Courant±Friedrichs condition requires that D 6 l20 =2t0 ; convergence also requires that the time step be small compared to the relaxation time of the system, a time which is roughly proportional to the inverse of the rate coecients. Thus, we see that both methods have restrictions on space and time steps. Alternatively, these restrictions can be seen as bounds on rate and di€usion coecients, if lattice spacing and time steps are the units of space and time. 7. Conclusions The code, algorithms, and derivations presented in the previous sections were for illustrative purposes, and in this presentation clarity took precedence over generality. In this section, we will see how LGCA can be generalised to other classes of problems. Increasing the maximum number of C or P particles past three simply requires more bits for each element of the array asState: a maximum of n immobile particles of one species requires log2 n, rounded up, bits. Having a maximum of two mobile particles, instead of one, moving through each of the channels to the six neighbours would require twelve bits per mobile species instead of six. In addition, bits can be reserved so that some of the X and Y species could remain in a hexagonal cell, with a certain probability. This naturally (and desirably) alters the

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

687

di€usion coecient. 8 In order to create a tuneable di€usion coecient in a model without rest particles, the random selection of channels in the distribution step can be replaced by probabilistic permutations of the original distribution of particles in channels as in [5,17]. For larger di€usion coecients, channels can lead to far away neighbours. To simulate a convective ¯ow in one direction, channels in that direction can be to hexagons many neighbours away [17]. Non-local, short-range interactions can also be included, and are useful for modelling non-linear e€ects especially in ecology, for example swarming of insects, schooling of ®sh, ¯ocking of birds, etc. [6,7]. Interaction of the particles with a ®eld can also be accomplished, for example by having position-dependent update rules that change in time. These can correspond to a force ®eld (for example in electromagnetics), or a concentration ®eld (in bacterial chemotaxis), or a changing environment (due to angiogenesis during cancerous tumour formation). Transport e€ects on boundaries, such as mechanical transport or facilitated di€usion, can be included with appropriate rules for nodes on boundaries. This is particularly relevant in biological cells, where many types of transport, for example in protein targetting and sorting, and selective passage of nutrients through a membrane, take place instead of simple di€usion. In addition, since PDE models for such systems have not in all cases been developed, the method presented here for relating LGCA micro-dynamics to macroscopic PDEs could aid in the development of continuous equations to describe such systems. We feel that the planned portability has been achieved. ANSI C compatability has been veri®ed by means of compiler settings (i.e., ANSI C violatation error enabled). Data format portability of the graphical images has been achieved by means of the many utilities that allow manipulation of GIF images. Readers intending to port this code into their applications, changing lattice size and number of species, should check that their application ®ts within the memory limits of their environment. We are interested in reading comments and suggestions from any who use this C code for their work. Please email to alawnicz@®elds.utoronto.ca and to [email protected].

Acknowledgements The authors gratefully acknowledge Karen Diemer for having taught Jon-Paul Voroney general LGCA coding in C, and Ray Kapral and Bruno Di Stefano for helpful discussions. Additionally, they acknowledge partial support from the Natural Sciences and Engineering Research Council of Canada and The Fields Institute for Research in Mathematical Sciences.

8 In LGCA with arrays of square cells instead of hexagons, a checkerboard parity results unless a stationary channel is included; this is due to particles initially in ``black'' square cells never interacting with particles initially in ``red'' square cells.

688

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

References [1] J.P. Boon, D. Dab, R. Kapral, A. Lawniczak, Lattice-gas automata for reactive systems, Phys. Rep. 273 (2) (1996) 55±147. [2] W. Boyce, R. DiPrima, Elementary Di€erential Equations and Boundary Value Problems, Wiley, Toronto, New York, 1986 (for the Runge±Kutta and Euler approximations). [3] V. Castets, E. Dulos, J. Boissonade, P. DeKepper, Experimental evidence of a sustained standing Turing-type non-equilibrium chemical pattern, Phys. Rev. Lett. 64 (1990) 2953±2957. [4] S. Chen, S.P. Dawson, G.D. Doolen, D.R. Janecky, A. Lawniczak, Lattice methods and their applications to reacting systems, Comput. Chem. Engrg. 19 (6/7) (1995) 617±646. [5] B. Chopard, M. Droz, Cellular automata model for the di€usion equation, J. Stat. Phys. 64 (1991) 859±892. [6] A. Deutsch, Orientation-induced pattern formation: swarm dynamics in a lattice-gas automaton model, Int. J. Bifurc. Chaos 6 (1996) 1735±1752. [7] A. Deutsch, A.T. Lawniczak, Probabilistic lattice models of collective motion and aggregation from individual to collective dynamics, Math. Biosci. 156 (1999) 255±269. [8] G. Doolen (Ed.), Lattice-Gas Methods for Partial Di€erential Equations, Addison-Wesley, New York, 1990. [9] G. Doolen (Ed.), Lattice-Gas Methods for PDE's Theory, Applications and Hardware, NorthHolland, Amsterdam, 1991. [10] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice-gas automata for the Navier±Stokes equation, Phys. Rev. Lett. 56 (1986) 1505±1509. [11] U. Frisch, D. d'Humieres, B. Hasslacher, Y. Pomeau, J.P. Rivet, Lattice-gas hydrodynamics in twoand three-dimensions, Complex Syst. 1 (1987) 649±707. [12] B. Hasslacher, R. Kapral, A. Lawniczak, Molecular Turing structures in the biochemistry of the cell, Chaos, Solitons & Fractals 3 (1993) 7±13. [13] B. Hess, A. Mikhailov, Self-organization in living cells, Science 264 (1994) 223±224. [14] R. Kapral, A. Lawniczak, P. Masiar, Oscillations and waves in a reactive lattice-gas automaton, Phys. Rev. Lett. 66 (1991) 2539±2542. [15] R. Kapral, A. Lawniczak, P. Masiar, Reactive dynamics in a multispecies lattice-gas automaton, J. Chem. Phys. 96 (4) (1992) 2762±2776. [16] A.T. Lawniczak, in: R. Kapral (Ed.), Pattern Formation and Lattice Gas Automata. The Fields Institute Communications, American Mathematical Society, 1996. [17] A.T. Lawniczak, Lattice-Gas Automata for Di€usive-Convective Dynamics, CNLS Newsletter, Los Alamos National Laboratory, vol. 136, 1997, pp. 1±15. [18] A.T. Lawniczak, J.P. Voroney, LGCA modeling of spatio-temporal dynamics of reactive systems, in preparation. [19] I. Lengyel, I.R. Epstein, Modeling of turing structures in the chlorite±iodide±malonic acid±starch reaction system, Science 251 (1991) 650±652. [20] J.D. Murray, Mathematical Biology, Springer, New York, 1989. [21] N. Margolous, T. To€oli, Cellular Automata Machines, A New Environment for Modeling, MIT Press, Cambridge, MA, 1987. [22] G. Marsaglia, B. Narasimhan, A. Zaman, A random number generator for PC's, Comput. Phys. Comm. 60 (1990) 345±349. [23] Q. Ouyang, H.L. Swinney, Transitions from a uniform state to hexagonal and striped Turing patterns, Nature 352 (1991) 610±612. [24] D.H. Rothman, S. Zaleski, Lattice-Gas Cellular Automata Simple Models of Complex Hydrodynamics, Cambridge University Press, Cambridge, UK, 1997. [25] A.B. Rovinsky, M. Mezinger, Chemical instability induced by a di€erential ¯ow, Phys. Rev. Lett. 69 (1992) 1193±1197. [26] P.H. Richter, P. Rehmus, J. Ross, Control and dissipation in oscillatory chemical engines, Prog. Theor. Phys. 66 (1981) 385±405. [27] E.E. Sel'kov, Self-oscillations in glycolysis. 1. A simple kinetic model, Eur. J. Biochem. 4 (1968) 79±86.

J.-P. Voroney, A.T. Lawniczak / Simulation Practice and Theory 7 (2000) 657±689

689

[28] T. To€oli, Cellular automata as an alternative to (rather than an approximation of) di€erential equations in modeling physics, Physica D 10 (1984) 117±127. [29] A. Turing, The chemical basis of morphogenesis, Philos. Trans. R. Soc. London B 237 (1952) 37±72. [30] J.P. Voroney, A. Lawniczak, R. Kapral, Turing pattern formation in heterogeneous media, Physica D 99 (1996) 303±317. [31] J.P. Voroney, Spatial and Temporal Patterns in Chemical Systems: Theoretical and Computational Approaches, PhD thesis, Department of Mathematics and Statistics, University of Guelph, 1998.