Mathematical Biosciences 156 (1999) 255±269

Probabilistic lattice models of collective motion and aggregation: from individual to collective dynamics Andreas Deutsch a

a,* ,

Anna T. Lawniczak

b,1

Theoretical Biology, Botanical Institute, University of Bonn, Kirschallee 1, D-53115 Bonn, Germany b Department of Mathematics and Statistics, Guelph-Waterloo Program for Graduate Work in Physics, University of Guelph, Guelph, Ont., Canada N1G 2W1 Received 17 March 1998; accepted 4 May 1998

Abstract We construct lattice gas automaton models to study collective dynamics and aggregation processes in biological systems. Corresponding transport equations are derived in appropriate space and time scaling under mean-®eld assumption and complemented with automation simulations. Ó 1999 Elsevier Science Inc. All rights reserved. Keywords: Collective motion; Aggregation; Cellular automaton model; Transport equation; Mean-®eld assumption

1. Introduction Biological motion typically comprises social aspects since individuals (cells or organisms) can change direction, speed etc. according to the motion characteristics of their spatial neighbours. Here, we study e€ects of local interactions which might lead to swarming or aggregation, i.e. stationary pattern formation. We provide mathematical descriptions of two simple biological interaction * Corresponding author. Tel.: +49-228 732 081; fax: +49-228 735 513; e-mail: [email protected] 1 Tel.: +1-519 824 4120; fax: +1-519 837 0221; e-mail: alawnicz@®elds.utoronto.ca

0025-5564/99/$ ± see front matter Ó 1999 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 9 8 ) 1 0 0 6 9 - X

256

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

automata (BIA), namely concentration-dependent and orientation-dependent automata. These automata are applied to study emergence of collective motion in biological systems. Other attempts to model aggregation can be found in Refs. [1±3]. Street formation in myxobacteria [4,5] is an example of orientationdependent motion (alignment) while cell migration driven by (di€erential) adhesion provides indication for concentration-dependent phenomena [6]. In the automata, motion is modelled as a discrete velocity jump process. This means that the movement of individuals of a modelled system is approximated by a piecewise linear motion which consists of runs at constant speed lasting for constant time intervals. In the automaton all individuals choose a velocity from a discrete space of allowable velocities, and continue with the corresponding direction for the next time interval. At the end of this interval all individuals simultaneously choose new directions and repeat the process. Several mechanisms of motion can then be analyzed by introducing various dependencies of the direction choice, in particular of neighbour density and orientations. Note that the automata can be viewed as discrete versions of continuous velocity jump processes as introduced by Othmer et al. [7,8]. The construction of the automata is inspired by the reactive lattice gas automata introduced in Refs. [9±12] as well as the interaction automata in Ref. [13]. Contrary to the automata considered in Refs. [9±12] the biological interaction automata are characterized by interactions that are non-local, they are neighbourhood-dependent. Automata with neighbourhood-dependent interactions were previously analyzed in Refs. [14±16]. They were used to study alignment in biological systems [16] and the e€ects of anti-di€usion [15]. Here, we will continue the investigation of self-organization in collective biological motion [17]. First, the automata are introduced in terms of microdynamical variables. Next, we derive the corresponding Boltzmann equations and from them the macroscopic partial di€erential equations. Finally, we compare simulations of the individual (automaton) and the population (discrete macroscopic Boltzmann equation) model. 2. Biological interaction automata: construction In the BIA automata biological organisms are represented by particles. For simplicity reasons, we consider biological systems with only one type of moving organisms, hence, one type of particles in BIA. However, the introduced BIA can be easily generalized to model biological systems with various types of moving and/or immobile organisms hence, various types of particles. We assume that in the BIA particles move on a lattice L with coordination number m ˆ 2 and periodic boundary conditions (extensions to other boundary conditions are straightforward). Furthermore, the BIA can be naturally extended to lattices with coordination number greater than 2.

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

257

We assume that each node of the lattice is labelled by a discrete vector r and is connected to its nearest neighbours r ‡ c; r ÿ c, by a unit vector c. With each node r we associate two channels (r, c) and …r; ÿc†. We denote them by r; ci , where i ˆ 1, 2 and c1 ˆ c and c2 ˆ ÿc. The particles reside at these lattice channels. They can leave their channels as a result of a probabilistic local interaction step or a deterministic propagation step. These steps are described below. The exclusion principle prevents two particles to be at the same channel. Hence a channel can have only two different con®gurations, and a node can have only four di€erent con®gurations. Due to the exclusion principle the presence or absence of a particle in a given channel …r; ci † at a given discrete time k can be represented by a Boolean random variable gi …r; k† de®ned as follows: gi …r; k† ˆ 1, if at time k there is a particle in the channel (r; ci ) and gi …r; k† ˆ 0, if at time k the channel (r; ci ) is empty. A con®guration of particles at a node r at time k can be described by a vector g…r; k† with components which are the Boolean random variables gi …r; k†, where i ˆ 1, 2. The random vector g…r; k† has values in a state space S ˆ fs ˆ hs1 ; s2 i where si ˆ 0 or 1; for i ˆ 1; 2g: The distribution of the particles on the lattice at a time k is described by a Boolean random ®eld fg…r; k†: r 2 Lg: The BIA dynamics arises from repetitive applications of superpositions of local interaction steps followed by a propagation step. In a propagation step at time k a particle residing in a channel …r; ci † moves to a neighbouring channel …r ‡ ci ; ci † for i ˆ 1,2. This step can be described easily in terms of microdynamical variables. Namely, if gPi …r ‡ ci ; k† denotes the contents of a channel …r ‡ ci ; ci † after the propagation step took place at time k, then gPi …r ‡ ci ; k ‡ 1† ˆ gi …r; k†;

…1†

for i ˆ 1, 2, respectively. The local interaction step is performed at each lattice node r independently of the other nodes. The rules governing this step for the concentration-dependent automaton are di€erent from those of the orientation-dependent automaton. They depend on the con®guration of particles at a given node and in its neighbourhood. The local interaction can take place at a given node r only if at this node there is one particle, if there are two particles or none, the local interaction does not take place. The con®gurations of particles in a neighbourhood of a node r determining local interactions at a node r are di€erent for the concentration and the orientation-dependent automaton. In order to characterize them, we introduce the following notation. We reserve a superscript c for the concentration-dependent automaton and a superscript o for the orientation-dependent automaton. We de®ne Dci …r; k† ˆ g1 …r ‡ ci ; k† ‡ g2 …r ‡ ci ; k†;

…2†

258

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

and Doi …r; k† ‡ gi …r ‡ c1 ; k† ‡ gi …r ‡ c2 ; k†:

…3†

Dci …r; k†

The variable gives the total number of particles at node r ‡ ci at time k, for i ˆ 1,2. The variable Doi …r; k† counts the total number of particles in the channels …r ‡ c1 ; ci † and …r ‡ c2 ; ci † at time k, for i ˆ 1,2. The local interactions in the concentration and orientation-dependent automaton (i.e., where s ˆ c or o) are characterized as follows. If at a node r there is only one particle, i.e. if g1 …r; k†…1 ÿ g2 …r; k†† ˆ 1 or …1ÿ g1 …r; k††g2 …r; k† ˆ 1, then the particle residing in the channel …r; c1 † can jump randomly to the channel …r; c2 † in the ®rst case, or from the channel …r; c2 † to the channel …r; c1 † in the second case. The probabilities of these transitions depend on the neighbourhood con®gurations. Let gi …r; k†…1 ÿ gj …r; k†† ˆ 1

…4†

Dsj …r; k† ÿ Dsi …r; k† P 1;

…5†

and where s ˆ c or o, and if i ˆ 1 then j ˆ 2, and vice versa. In this case a particle can jump from the channel …r; ci † to the other (empty) channel …r; cj † at a node r with large probability pijs , where s ˆ c,o, and as before, if i ˆ 1 then j ˆ 2, and vice versa. The probabilities pijs , for ij ˆ 12 and 21, are the transition probabilities, such that if Eqs. (4) and (5) are satis®ed for i ˆ 1 and j ˆ 2, then s s ˆ p12 …10; 01† p12

…6†

is a transition probability from the state h1; 0i to the state h0; 1i, and if Eqs. (4) and (5) are satis®ed for i ˆ 2 and j ˆ 1, then s s ˆ p21 …01; 10† p21

…7†

is a transition probability from the state h0; 1i to the state h1; 0i. If Dsj …r; k† ÿ Dsi …r; k† 6 0;

…8†

then the particle can jump from the channel …r; ci † to the other (free) channel …r; cj † at a node r with a very small probability qsij , where s ˆ c, o, and if i ˆ 1 then j ˆ 2, and vice versa. If Eqs. (4) and (8) are satis®ed for ij ˆ 12 and 21, respectively, then the probabilities qsij , for ij ˆ 12, or 21, are the transition probabilities, such that qs12 ˆ qs12 …10; 01†;

…9†

and …10† qs21 ˆ qs21 …01; 10†: They describe transitions from the state h1; 0i to the state h0; 1i, and from the state h0; 1i to the state h1; 0i, respectively, for s ˆ c, o. The local transitions are

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

259

governed at time k at each node r, independently of the other nodes and a past evolution of the automaton by independent Bernoulli type random variables nsij …r; k† and csij …r; k†, such that P …nsij …r; k† ˆ 1† ˆ pijs ;

…11†

P …csij …r; k† ˆ 1† ˆ qsij ;

…12†

for s ˆ c,o and ij ˆ 12 and 21. Let 1Asji …r; k† be the characteristic function of the event Asji …r; k† ˆ fDsj …r; k† ÿ Dsi …r; k† P 1g and 1Bsji …r; k† be the characteristic function of the event Bsji …r; k† ˆ fDsj …r; k† ÿ Dsi …r; k† 6 0g; where s ˆ c,o, and if i ˆ 1 then j ˆ 2, and vice versa. It is easy to see that 1Asji …r; k† ‡ 1Bsji …r; k†  1:

…13†

Depending on the type of the automaton the functions 1Asji …r; k† have the following representations. For the concentration-dependent automaton, 1Acji …r; k† ˆ

2 2 Y Y gl …r ‡ cj ; k†…1 ÿ gl …r ‡ ci ; k†† ‡ gl …r ‡ cj ; k† lˆ1

lˆ1

 ‰g1 …r ‡ ci ; k†…1 ÿ g2 …r ‡ ci ; k†† ‡ …1 ÿ g1 …r ‡ ci ; k††g2 …r ‡ ci ; k†Š ‡ ‰g1 …r ‡ cj ; k†…1 ÿ g2 …r ‡ cj ; k†† ‡ …1 ÿ g1 …r ‡ cj ; k††g2 …r ‡ cj ; k†Š 

2 Y lˆ1

…1 ÿ gl …r ‡ ci ; k††;

…14† where if i ˆ 1 then j ˆ 2, and vice versa. For the orientation-dependent automaton, 1Aoji …r; k† ˆ

2 2 Y Y gj …r ‡ cl ; k†…1 ÿ gi …r ‡ cl ; k†† ‡ gj …r ‡ cl ; k† lˆ1

lˆ1

 ‰gi …r ‡ c1 ; k†…1 ÿ gi …r ‡ c2 ; k†† ‡ …1 ÿ gi …r ‡ c1 ; k††gi …r ‡ c2 ; k†Š ‡ ‰gj …r ‡ c1 ; k†…1 ÿ gj …r ‡ c2 ; k†† ‡ …1 ÿ gj …r ‡ c1 ; k††gj …r ‡ c2 ; k†Š 

2 Y lˆ1

…1 ÿ gi …r ‡ cl ; k††; …15†

260

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

where if i ˆ 1 then j ˆ 2, and vice versa. The results of local interactions at a node r at time k can be described in s terms of microdynamical variables as follows. If gIi …r; k† denotes a microdynamical variable describing a content of the channel …r; ci † after a local interaction took place, then s

gIi ˆ gi gj ‡ gi …1 ÿ gj †1Asji …1 ÿ nsij † ‡ gi …1 ÿ gj †1Bsji …1 ÿ csij † ‡ …1 ÿ gi †gj 1Asij nsji ‡ …1 ÿ gi †gj 1Bsij csji

…16†

for s ˆ c,o. After some simple algebra s

gIi ˆ gi ÿ gi …1 ÿ gj †‰1Asji nsij ‡ 1Bsji csij Š ‡ …1 ÿ gi †gj ‰1Asij nsji ‡ 1Bsij csji Š;

…17†

where s ˆ c,o and if i ˆ 1 then j ˆ 2, and vice versa. In the preceding two formulas we dropped the dependence on …r; k†. We will follow this convention if confusion is unlikely to occur. 3. Macroscopic equations Let us recall that the full automaton dynamics arises from repetitive applications of superpositions of local interaction steps, applied at each node r of the lattice independently of the other nodes, followed by a propagation step. By Eqs. (1) and (17) the superposition of a local interaction step with a propagation step can be expressed in terms of microdynamical variables as follows: s

…18† gi …r ‡ ci ; k ‡ 1† ˆ gIi …r; k† for s ˆ c,o, and i ˆ 1,2. Since all the random variables ns0 s and cs0 s governing local interactions are mutually independent, then by construction the entire Boolean ®eld fg…r; k†: r 2 L; k ˆ 0; 1; 2; . . .g is a Markov stochastic process, with some initial distribution ls , for s ˆ c,o. Let Els be the expected value with respect to the distribution of this Markov stochastic process, then …19† fis …r; k† ˆ Els gi …r; k† is the average occupation number, for s ˆ c,o, and i ˆ 1,2. By Eq. (17) taking the expected value of Eq. (18) we get fis …r ‡ ci ; k ‡ 1† ÿ fis …r; k† o n ˆ Els ÿ gi …1 ÿ gj †‰1Asji nsij ‡ 1Bsji csij Š ‡ …1 ÿ gi †gj ‰1Asij nsji ‡ 1Bsij csji Š ; where the right hand side is evaluated at …r; k†:

…20†

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

261

Assuming molecular chaos, i.e. mean-®eld approximation in which the expected value of a product of random variables is the product of their expectations, we obtain from Eqs. (14), (15) and (20), that fis …r ‡ ci ; k ‡ 1† ÿ fis …r; k† s

s

ˆ ÿfis …r; k†…1 ÿ fjs …r; k††‰A~ji …r; k†pijs ‡ B~ji …r; k†qsij Š ‡ …1 ÿ

fis …r; k††fjs …r; k†

…21†

s s  ‰A~ij …r; k†pjis ‡ B~ij …r; k†qsji Š

for s ˆ c,o and i ˆ 1,2, where as usual, if i ˆ 1 then j ˆ 2, and vice versa. In s s s s Eq. (21) A~ij …r; k†; A~ji …r; k†; B~ij …r; k†; B~ji …r; k† are expected values of the random variables 1Asij …r;k† ; 1Asji …r;k† ; 1Bsij …r;k† ; 1Bsji …r;k† ; respectively, which are calculated under molecular chaos assumption. They are obtained from Eqs. (13)±(15) by replacing in these formulas gi s; gj s and gl s by fis s; fjs s; and fls s, respectively, for s ˆ c,o, and if i ˆ 1, then j ˆ 2, and vice versa. Note that in Eqs. (14) and (15) gi ; gj ; and gl are evaluated at the nodes r ‡ c1 and r ‡ c2 . Hence, the terms s s s s fis s; fjs s; and fls s in the expressions A~ij ; A~ji ; B~ij ; B~ji ; are also evaluated at the nodes r ‡ c1 and r ‡ c2 . In Eq. (21) the ®rst term expresses the probability in mean ®eld approximation of loosing a particle from the channel …r; ci †, due to its jump to the empty channel …r; cj †; the second term expresses the probability in mean ®eld approximation of gaining a particle at the empty channel …r; ci †, due to its jump from the occupied channel …r; cj †. The probability of loosing or gaining a particle from or in the channel …r; ci † is neighbourhood-dependent. This dependence is expressed through the transition probabilities pijs ; pjis ; qsij ; qsji s s s s and the terms A~ij ; A~ji ; B~ij ; B~ji ; which give probabilities in mean ®eld approximation of various neighbourhood con®gurations at the nodes fr ‡ c2 ; r ‡ c1 g. The form of these terms depends on the type of the automaton. Eq. (21) can be also expressed in terms of turning frequencies ksij …r; k† and s kji …r; k†, de®ned by s

s

ksij …r; k† ˆ …1 ÿ fjs …r; k††‰A~ji …r; k†pijs ‡ B~ji …r; k†qsij Š;

…22†

s s ksji …r; k† ˆ …1 ÿ fis …r; k††‰A~ij …r; k†pjis ‡ B~ij …r; k†qsji Š;

…23†

which are equivalents of the turning frequencies of Ref. [4]. They depend, of course, on the average concentration pro®le in the perception set fr ‡ c2 ; r; r ‡ c1 g, and using them Eq. (21) becomes of the same form as in Ref. [4], namely fis …r ‡ ci ; k ‡ 1† ÿ fis …r; k† ˆ ÿksij …r; k†fis …r; k† ‡ ksji …r; k†fjs …r; k†:

…24†

s s A~ij …r; k† ‡ B~ij …r; k†  1;

…25†

Since

262

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

if i ˆ 1 and j ˆ 2, and vice versa, observe that, if pijs ˆ qsij for ij ˆ 12 and 21, then the rules of the BIA are neigbourhood-independent and Eq. (21) becomes fis …r ‡ ci ; k ‡ 1† ÿ fis …r; k† ˆ ÿfis …r; k†…1 ÿ fjs …r; k††pijs ‡ …1 ÿ fis …r; k†fjs …r; k†pjis :

…26†

Eqs. (21), (24) and (26) are discrete space and time Boltzmann equations of the concentration- and orientation-dependent BIA. These equations show that, if pijs 6ˆ qsij for ij ˆ 12 and 21, then the space and time evolution of the automaton concentrations fis and fjs is neighbourhood-dependent. The discrete Boltzmann equations (21), (24) and (26) in appropriate space and time scaling, can be reduced to the macroscopic partial di€erential equations of the concentration- and orientation-dependent BIA. These equations, derived under molecular chaos assumption, describe the `limiting behaviour' of the BIA when the number of sites in the lattice L goes to in®nity. The `limiting behaviour' emerges on macroscopic scales from the BIA microscopic dynamics when the distance between neighbouring sites, time between biological interactions and ¯ipping rates go to zero. This `limiting behaviour' depends on the relationships between considered time, space and rate scalings. It is re¯ected through the types of emerging macroscopic PDEs which can be di€erent for di€erent scalings. Here, we consider Euler scaling whereby time and space scale at the same rates as opposed to di€usive scaling at which time scales as a square root of the spatial length (discussed in Ref. [18]). In order to derive these macroscopic equations we assume the following. We assume that the BIA dynamics takes place on a one-dimensional lattice L ˆ L  Z, centered at the origin, with number of sites equal to two times the integer part of ÿ1 and periodic boundary conditions. This implies that various BIA parameters become -dependent. Furthermore, we assume that the initial distribution ls ˆ ls is such that Els …gi …r; 0†† ˆ gis …r†; gis …x†,

…27†

with x 2 ‰ÿ1; 1Š, is a non-negative smooth function independent of where , bounded by one and with periodic boundary conditions, i.e. gis …ÿ1† ˆ gis …1†, for s ˆ c,o and i ˆ 1,2. For r 2 L a microscopic point, x ˆ r is the corresponding macroscopic point. In other words, each -discretization of the interval [ÿ1, 1] is mapped onto its corresponding lattice L . In this way, through the functions gis , we impose an initial concentration pro®le on each lattice L . The functions gis vary on a distance of order of  on the interval [ÿ1, 1]. Hence, after rescaling each -discretization of the interval [ÿ1, 1] by ÿ1 , on the corresponding lattice L the initial concentration pro®le varies on average on a distance of order of 1, i.e. from the node to node. Since biological interactions in BIA take place only at the lattice nodes this implies that the e€ective mean free path of the BIA, i.e. the distance free of biological interactions, is of order 1 on a lattice; hence, of order  on macroscopic scales, i.e. on the interval

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

263

[ÿ1, 1]. Since the microscopic space scale was rescaled by  and there are ®nite velocities in the system we must also rescale time. Namely, the microscope time k corresponds to the macroscopic time t ˆ k in the considered kinetic regime characterized by a Knudsen number (the ratio of the mean free path and the typical distance on which densities vary) of order one [19]. Additionally, we assume that transition probabilities, representing ¯ipping rates, scale as pijs ! pijs and qsij ! qsij , for ij ˆ 12 and ij ˆ 21. Under the above scaling Eq. (21) becomes fis …x ‡ ci ; t ‡ † ÿ fis …x; t† s s ˆ ÿfis …x; t†…1 ÿ fjs …x; t††‰ A^ji …x; t†pijs ‡ B^ji …x; t†qsij Š s s ‡ …1 ÿ fis …x; t††fjs …x; t†‰ A^ij …x; t†pjis ‡ B^ij …x; t†qsji Š;

…28†

for s ˆ c, o, i ˆ 1, 2, and as usual, if i ˆ 1 then j ˆ 2 and vice versa. The terms s s s s s s s s A^ij ; A^ji ; B^ij and B^ji are obtained from the terms A~ij ; A~ji ; B~ij and B~ji by replacing in them the arguments …r ‡ ci ; k† and …r ‡ cj ; k† of the average occupation numbers fis ; fjs ; fls by their macroscopic corresponding points …x ‡ ci ; t† and …x ‡ cj ; t†, respectively, for i, j, l ˆ 1, 2. By taking the truncated Taylor expansions of fis0 s; fjs0 s; and fls0 s up to order 1 in Eq. (28), for s ˆ c, o, and i, j, l ˆ 1, 2, we reduce Eq. (28), in a limit  ! 0, to two coupled partial differential equations de®ned on the interval [ÿ1, 1], namely ot f1s …x; t† ‡ ox f1s …x; t† s ˆ ÿf1s …x; t†…1 ÿ f2s …x; t††‰As21 …x; t†p12 ‡ Bs21 …x; t†qs12 Š s ‡ …1 ÿ f1s …x; t††f2s …x; t†‰As12 …x; t†p21 ‡ Bs12 …x; t†qs21 Š

…29†

and ot f2s …x; t† ÿ ox f2s …x; t† s ‡ Bs12 …x; t†qs21 Š ˆ ÿf2s …x; t†…1 ÿ f1s …x; t††‰As12 …x; t†p21 s ‡ …1 ÿ f2s …x; t††f1s …x; t†‰As21 …x; t†p12 ‡ Bs21 …x; t†qs12 Š;

…30†

with the initial conditions f1s …x; 0† ˆ g1s …x† and f2s …x; 0† ˆ g2s …x†:

…31†

This initial boundary value problem is describing the macroscopic dynamics of the BIA, for s ˆ c, o. The terms As12 ; As21 ; Bs12 and Bs21 are obtained from the s s s s expressions A^12 ; A^21 ; B^12 and B^21 of Eq. (28) by replacing in them the arguments …x ‡ ci ; t† and …x ‡ cj ; t†, for i, j ˆ 1, 2, of the average occupation numbers fis s; fjs s; and fls s by (x, t), for i, j, l ˆ 1, 2. This is the result of taking the limit  ! 0, of the truncated Taylor approximations up to order 1 of fis s; fjs s; and fls s: Hence, the terms As12 ; As21 , for s ˆ c, o become

264

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

Ac12 …x; t† ˆ Ac21 …x; t† ˆ f1c f2c …1 ÿ f1c †…1 ÿ f2c † ‡ ‰f1c …1 ÿ f2c † ‡ …1 ÿ f1c †f2c Š‰f1c f2c ‡ …1 ÿ f1c †…1 ÿ f2c †Š; 2

2

…32† 2

2

Ao12 …x; t† ˆ …f1o † …1 ÿ f2o † ‡ 2…f1o † f2o …1 ÿ f2o † ‡ 2f1o …1 ÿ f1o †…1 ÿ f2o † ; …33† 2

2

2

2

Ao21 …x; t† ˆ …f2o † …1 ÿ f1o † ‡ 2…f2o † f1o …1 ÿ f1o † ‡ 2f2o …1 ÿ f2o †…1 ÿ f1o † ; fis

where the average occupation numbers From Eqs. (25) and (32) it follows that

ˆ

fis …x; t†,

…34† for s ˆ c, o and i, j ˆ 1, 2.

Bc12 …x; t† ˆ 1 ÿ Ac12 …x; t† ˆ 1 ÿ Ac21 …x; t† ˆ Bc21 …x; t†;

…35†

Bo12 …x; t†

…36†

ˆ1ÿ

Ao12 …x; t†;

and Bo21 …x; t† ˆ 1 ÿ Ao21 …x; t†: …37† In the macroscopic partial di€erential equations (29) and (30) the neighbourhood-dependence of the BIA rules is manifested through the above As and Bs terms. The fact that the average occupation numbers f1s ; f2s , appearing in As and Bs terms lost their neighbourhood-dependence in their argument and are evaluated at the point (x, t) is the result of taking the limit,  ! 0, of their Taylor approximations. Observe that, if pijs ˆ qsij for ij ˆ 12 and 21, Eqs. (29) and (30) become s s ot f1s …x; t† ‡ ox f1s …x; t† ˆ ÿf1s …x; t†…1 ÿ f2s …x; t††p12 ‡ …1 ÿ f1s …x; t††f2s …x; t†p21

…38† and s s ot f2s …x; t† ÿ ox f2s …x; t† ˆ ÿf2s …x; t†…1 ÿ f1s …x; t††p21 ‡ …1 ÿ f2s …x; t††f1s …x; t†p12 :

…39† Eqs. (38) and (39) could be obtained directly from Eq. (26). If in addition pijs ˆ pjis ˆ ps , i.e. if the ¯ipping rates are equal, then Eqs. (38) and (39) are reduced to ot f1s …x; t† ‡ ox f1s …x; t† ˆ …ÿf1s …x; t† ‡ f2s …x; t††ps

…40†

and …41† ot f2s …x; t† ÿ ox f2s …x; t† ˆ …ÿf2s …x; t† ‡ f1s …x; t††ps : The transport equations (38)±(41) depend only on the probabilities of a particle of being or not at a given channel at a macroscopic point x at time t and on the ¯ipping rates. This is in agreement with that fact, observed earlier in Eq. (26), that in these cases the BIA rules become neighbourhood-independent, i.e. the

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

265

jumping frequencies are the same regardless of the neighbourhood con®guration. Since the neighbourhood-dependence has been lost, the form of the macroscopic partial di€erential equations is the same for the concentrationand orientation-dependent automaton. Alternatively, if in the derivation of Eqs. (29) and (30) we scale the ¯ipping rates as 2 instead as , i.e. pijs ! 2 pijs and qsij ! 2 qsij , then we obtain equations ot f1s …x; t† ‡ ox f1s …x; t† ˆ 0;

…42†

ot f2s …x; t†

…43†

ÿ

ox f2s …x; t†

ˆ 0:

In this case, the biological interactions happen on a slower scale and their e€ects can not be detected by the macroscopic equations of the BIA. Further analysis of various space and time scales under which macroscopic PDEs of BIA can be derived is under investigation [18]. 4. Simulations We present simulation results for concentration- and orientation-dependent BIA. Since neighbourhood-dependence of the dynamics is only manifested if pijs 6ˆ qsij for ij ˆ 12 and 21, both in the rules of the automata and in the corresponding macroscopic Eqs. (29) and (30) (through the terms As and Bs ), we direct our attention to this situation. First, we consider the orientation- dependent case …s ˆ o†. Starting out from a spatially homogeneous initial state automaton simulations clearly exhibit formation of oriented patches collectively moving over their lattice habitat (Fig. 1, top). Coexistence of di€erently oriented patches persists for long times with alternating predominance of one orientation or the other (left or right). Contrary, simulations with the same parameter set of the corresponding discrete macroscopic Boltzmann equations (28) in the long run always yield growth of a single direction at the cost of the other (Fig. 2, left). Presumably, correlations and ¯uctuations prevent convergence to a single predominant orientation (mode) in the automaton. In the case of the concentration-dependent automaton …s ˆ c† retaining the same parameters and a similar initial distribution of particles as for the orientation-dependent model, simulations show formation of spatial patterns visible as vertical bands (Fig. 1, bottom). Note that these bands do not have a distinguished spatial frequency. In contrast, in a typical simulation of the corresponding macroscopic equations (28), illustrated in Fig. 2 (right), we always observe periodic pattern formation. In the long run f1 …x; t† ˆ f2 …x; t† is 2-periodic where the period length corresponds to the neighbourhood range. Obviously, since correlations have been neglected in the derivation of the macroscopic equations (28), we can not expect that these equations will capture correctly all qualitative behaviours manifested by the underlying

266

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

Fig. 1. Temporal evolution of biological interaction automata starting from homogeneous initial  ˆ 0.5). Motion tracks condition (left vertical line t ˆ 0, averaged density of particles per node q distinguish right and left moving particles. Only initial …0 6 t 6 300† and ®nal evolution …9700 6 t 6 tmax ˆ 10 000† are shown. Top: orientation model …s ˆ o†, street formation in both directions (collective motion); bottom: concentration-dependent model …s ˆ c†, aggregation in the s s ˆ p21 ˆ 0:9; qs12 ˆ qs21 ˆ 0:01, lattice size ˆ 100, form of stripe pattern formation. Parameters: p12 periodic boundary conditions.

Fig. 2. Numerical simulation of Eq. (28) starting from randomly perturbed homogeneous initial condition (f1 …x; 0†  f2 …x; 0†  0:5†. Left: orientation model …s ˆ o†; right: concentration model s s ˆ p21 ˆ 0:9; qs12 ˆ qs21 ˆ 0:01;  ˆ 1=100, periodic boundary …s ˆ c†, Parameters (as in Fig. 1): p12 conditions.

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

267

concentration-dependent BIA. Further analysis of these patterns, including linear stability analysis is presented in Ref. [18]. 5. Discussion We have introduced biological interaction automata and derived (in mean®eld approximation) corresponding macroscopic equations under certain assumptions with regard to space, time and ¯ipping rates scaling. Thereby, it could be demonstrated that di€erences between the individual and the population level model description may be signi®cant, particularly in the concentration-dependent model. Here, the macroscopic equations always yield periodic banding patterns with a well-de®ned frequency depending on the transport/interaction range of the model. On the contrary, at the automaton level, we observe formation of stripes of various widths and various life spans. These unregularities are due to strong correlations induced by the low dimensionality of the state space. In particular, this suggests that the mean-®eld description of the one-dimensional concentration-dependent automaton does not properly describe its spatio-temporal dynamics. It seems that the molecular chaos assumption is too crude, due to the strong correlations in the one-dimensional model. However, it may be appropriate in higher dimensions, i.e. systems with less in¯uence of correlations. Such systems can be easily generated by considering higher spatial dimensions or, alternatively, by allowing more particles per node. For example, it is known from two-dimensional versions of the concentration-dependent (anti-di€usion) automaton (Alexander model) that stationary patterns are often even periodic [15]. This result is also represented by linear stability analysis under mean-®eld assumption, i.e. in two dimensions the corresponding instabilities are strong enough even under the assumption of uncorrelated particle distributions. On the other hand, in the orientation-dependent automaton we observe pattern formation both at the individual and the population description of the model. In Ref. [20] three di€erent scenarios characterizing orientation-dependent pattern formation are distinguished. Starting from an initially disordered state instabilities might lead to: (A) formation of angular order in a spatially homogeneous distribution; (B) aggregation and formation of regular clusters without common direction and (C) development of orientation waves in patches of aligned entities. Both (A) and (C) are possible at the individual level while at the population level always scenario (A) is present. If we include a concentration-dependence in the automaton dynamics, most probably, also scenario (B) would become possible. Scenario (B) has been already observed both in simulations and in the corresponding stability analysis of a model with rest particles (i.e. zero velocity) [21].

268

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

A systematic analysis and comparison of stability results is currently undertaken for both the concentration and the orientation model [18]. So far, we have only shown simulations of the discrete macroscopic Boltzmann equations (28), systematic examination of the continuous Boltzmann equations (29) and (30) remains to be done. Comparison of the resulting hyperbolic transport equation with other reaction transport equations (e.g. those considered in Ref. [22]) is also a topic of future investigation. In this paper, only concentration and orientation-dependent non-local rules were considered. Nevertheless, it is straightforward to extend the analysis to models in higher dimensions, with more particle types and with more complex interactions. An example with three cell types and hybrid concentration/orientation interactions has already proven useful in simulating pigment pattern formation of certain salamander larvae which can be explained by di€erential adhesion and local alignment dynamics of certain moving and resting cells [23]. This example nicely demonstrates how elementary model modules of the type in the focus of this article may be combined in order to explain key steps in the embryology even of multi-cellular organisms.

Acknowledgements The authors greatfully acknowledge valuable discussions with Michael Wurzel and Wolfgang Alt, Bonn, as well as ®nancial support from SFB 256 under the heading Nonlinear partial di€erential equations. Additionally, A.T.L. acknowledges partial support from the Natural Sciences and Engineering Research Council of Canada and The Fields Institute for Research in Mathematical Sciences.

References [1] D. Gr unbaum, A. Okubo, Modelling social animal aggregations, manuscript. [2] K.P. Hadeler, Random walk systems modeling spread and interaction, CIME lectures at Martina Franca, June, 1997. [3] W. Alt, A. Deutsch, G. Dunn (Eds.), Dynamics of Cell and Tissue Motion, Birkh auser, Basel, 1997. [4] B. P®stner, A one-dimensional model for the swarming behaviour of myxobacteria, in: W. Alt, G. Ho€mann (Eds.), Biological Motion, Springer, Heidelberg, 1990, p. 556. [5] B. P®stner, Simulation of the dynamics of myxobacteria swarms based on a one-dimensional interaction model, J. Biol. Syst. 3 (2) (1995) 579. [6] M.S. Steinberg, Does di€erential adhesiveness govern self-assembly processes in histogenesis? Equilibrium con®gurations and emergence of a hierarchy among populations of embryonic cells, J. Exp. Zool. 173 (1970) 395.

A. Deutsch, A.T. Lawniczak / Mathematical Biosciences 156 (1999) 255±269

269

[7] S.R. Dunbar, H.G. Othmer, On a nonlinear hyperbolic equation describing transmission lines, cell movement, and branching random walks, in: H.G. Othmer (Ed.), Nonlinear Oscillations in Biology and Chemistry, Springer, Berlin, 1986, p. 274. [8] H.G. Othmer, S.R. Dunbar, W. Alt, Models of dispersal in biological systems, J. Math. Biol. 26 (1988) 263. [9] A. Lawniczak, D. Dab, R. Kapral, J.P. Boon, Reactive lattice gas automata, Physica D 47 (1991) 132. [10] R. Kapral, A. Lawniczak, P. Masiar, Oscillations and waves in a reactive lattice-gas automaton, Phys. Rev. Lett. 66 (1991) 2539. [11] R. Kapral, A. Lawniczak, P. Masiar, Reactive dynamics in a multi-species lattice-gas automaton, J. Chem. Phys. 96 (1992) 2762. [12] J.P. Boon, D. Dab, R. Kapral, A. Lawniczak, Lattice gas automata for reactive systems, Phys. Rep. 273 (2) (1996) 55. [13] A. Deutsch, Towards analyzing complex swarming patterns in biological systems with the help of lattice-gas cellular automata, J. Biol. Syst. 3 (1995) 947. [14] A. Deutsch, Orientation-induced pattern formation: swarm dynamics in a lattice-gas automaton model, Int. J. Bifruc. Chaos 6 (1996) 1735. [15] H.J. Bussemaker, Analysis of a pattern-forming lattice-gas automaton: mean ®eld theory and beyond, Phys. Rev. E 53 (1996) 1644. [16] H. Bussemaker, A. Deutsch, E. Geigant, Mean-®eld analysis of a dynamical phase transition in a cellular automaton model for collective motion, Phys. Rev. Lett. 78 (1997) 5018. [17] L. Edelstein-Keshet, G. Ermentrout, Models for contact-mediated pattern formation: cells that form parallel arrays, J. Math. Biol. 29 (1990) 33. [18] A. Deutsch, A. Lawniczak, Characterization of instabilities in lattice models of collective motion, manuscript in preparation. [19] C. Carcignani, The Boltzmann Equation and its Applications, Springer, Berlin, 1988. [20] A. Mogilner, A. Deutsch, J. Cook, Models for spatio-angular self-organization in cell biology, in: A. Deutsch, W. Alt, G.Dunn (Eds.), Dynamics of Cell and Tissue Motion, Ch. III.1, Birkhauser, Basel, 1997, p. 173. [21] A. Deutsch, A new mechanism of aggregation in a lattice-gas cellular automaton model, Math. Comp. Mod. to appear. [22] T. Hillen, A Turing model with correlated random walk, Report No. 9 at Sonderforschungsbereich 382, Universities of T ubingen and Stuttgart, 1995. [23] A. Deutsch, M. Spielberg, L. Olsson, Cell-based modeling of salamander pigment pattern formation, manuscript in preparation.