Reaction-Diffusion a la Turing in MGS


Turing was the first to suggest (in 1952) that the mechanisms of morphogenesis can be explained by reaction-diffusion equations. The idea is that each cell in a tissu (e.g. the skin) is characterized by two numbers, representing concentrations of substances (morphogens) a and b. The substances diffuse and react with each other following a partial differential equations of the form

where the vector-valued at each time t describes the changing spatial concentrations of a and b. The function specifies the reaction and the laplacian term describes the diffusion. Numerous variants have been studied and widely simulated. If one of the morphogen is colored, these systems have been shown to produce regular patterns in the tissu, reminiscent of zebra stripes and leopard spots. See here various example of reaction-diffusion patterns that can be produced. Turing's analysis stimulated considerable theoretical research on mathematical models of pattern formation, but Turing-type patterns were not observed in controlled laboratory experiments until 1990. This kind of model have also been used not only to explain the marking on animals, but reaction-riffusion processes have been used in areas such as:

Whether or not reaction-diffusion mechanisms are an actual phenomena in such systems, simulation does tend to generate regular patterns that appear as standing-wave solutions to the equations. This offers an explanation of the breakdown of symmetry and homogeneity or the emergence of a pattern in an originally homogenous medium.

We propose two examples coded in MGS in 1D and 2D.

Diffusion-reaction in a ring of cells (1D)

reaction-diffusion in a ring of cell

click for a larger view (3Mb animated gif)

This bubling torus is a visual representation of a diffusion-reaction process "à la Turing" on a ring of cells. A very simple MGS program is used to compute the density of two chemicals that diffuse and react on a ring of cells. The data are then written in a file together with a simple Mathematica program used to generate this plot. The diameter of the torus is proportional to the density of one of the chemicals.

The corresponding MGS program

specifying a reaction-diffusion in 1D

gbf RING = < d; 40 d = 0 >;;

NbCell := 40;;
Definition of a GBF type. GBF means Group Based Fields. A GBF specifies a uniform topology where the elements are indexed by a mathematical group. Here, all elements in an instance of the GBF  type RING has two neighbors: one following the direction d and the other following the inverse direction. If we go 40 steps following the direction d, we reach again our starting point, hence the equation  40 d = 0.

This is juste the way in MGS to define a circural buffer (a ring of cells) and a displacement nammed d.
rsp := 1.0/16.0;;
diff1 := 0.25;;
diff2 := 0.0625;;

fun da(a, b, la, ra) =
      rsp * (16.0 - a * b) + diff1*(la + ra - 2.0*a);;
fun db(a, b, beta, lb, rb) =
      rsp*(a*b - b - beta) + diff2*(lb + rb - 2.0*b);;
Definition of some constant parameters.

Definition of the function da and db that computes the change in a and b. Arguments a and b are the concentration in morphogen a and b at somme cell xla is the concentration of morphogen a at the cell at left of x and ra the concentration at the cell at right of x. Identically for b, lb and rbbeta is a constant for each cell that introduces some fluctuations.
trans Turing =
  x => let xl = hd('neighbor(x, <d|))
       and xr = hd('neighbor(x, |d>)) in
    { a = x.a + da(x.a, x.b, xl.a, xr.a),
      b = max(0.0, x.b + db(x.a, x.b, x.beta, xl.b, xr.b)),
      beta = x.beta
The state of a cell is defined by a record x containing the concentration in a, in b and the parameter beta. The special function 'neighbor(x, dir) returns all the neighbors of the cell x following the direction dir. It returns a sequence. Here the sequence returned is a singleton (there is only one neighbor following a given direction in a GBF) and we take this element using hd (the head of the sequence).
So, xl and xr denotes the element at the left and at the right of x.
A new record is computed using the da and the db function to compute the new values of a and b.
trans init =
  x => { a = 4.0,
         b = 4.0,
         beta = 12.0 + random(0.1) - 0.05 };

ring0 := init(iota(NbCell, ():seq) following |d>::seq:());;
To compute an initial state, we apply the transformation init to an initial ring. We do not care the elements in this ring, we rather replace each element by a record containing some predefined value for the concentration in a and b and a beta value which fluctuates around 12 +/- 0.05.

The transformation init is applied to an initial ring obtained by coercing a sequence of 40 integer (computed by the function iota) into a ring. This coercion makes use of the function following. Following is an infix function that takes two arguments. The first argument is a sequences of elements and the second is a sequence of directions that specify the distribution of the elements into the ring.
This expression computes 150 time steps of the process.

Reaction-diffusion at the surface of a torus (2D)

diffusion reactiojn at the surface of a torus.

click for a larger view (1.5 Mb)

In this view, the radius of a point at the surface of the torus is proportional to one of the two morphogens that diffuse and react at the surface. The deformation of the torus is simply used to visualize the concentration. The diffusion does not take place in a deformable medium but at the fixed surface of an ordinary torus.
The reaction-diffusion equations are taken from a model proposed by G. Turk in 1991to generate textures on arbitrary surface (Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion,  Computer Graphics, Vol. 25, No. 4, pp. 289-298, July 1991 (SIGGRAPH '91)).

The corresponding MGS program

the 2D diffusion-reaction

gbf torus = < a, b; 180 a, 60 b>;;

NbCell_a := 180;;
NbCell_b := 60;;

Definition of a GBF type. GBF means Group Based Fields. A GBF specifies a uniform topology where the elements are indexed by a mathematical group. Here, all elements in an instance of the GBF  type torus has two neighbors: one following the direction a, the other following the direction b. (and also two others neighbors reached following the inverse directions). If we go 180 steps following the direction a, we reach again our starting point, hence the equation  180 a = 0 which is abreviated as 180 a. The same property holds for the direction b but after 60 steps.

This is juste the way in MGS to define a circural 2D grid (a torus) whose dimensions are named a and b.
beta_init := 12;;
beta_rand := 0.1;;

u_steady := 4;;
v_steady := 4;;

diffu := 0.25 / 2.0;;
diffv := 0.0625 / 2.0;;

p1 := 0.2;;
p2 := 0.0;;
p3 := 0.0;;

k := p1/16.0;;

dt := 0.4;;
Some consyant parameters are declared (see the article of G. Turk for an explanation of these parameters).

Each point if the torus is associated with a record with fields u, v and c. The first two fields record the concentration in morphogens.Think to this record as a vector. The field c is constant in timeand correspond to a perturbation term.
fun acc(r, suite) = { u = r.u + suite.u,  v = r.v + suite.v}

An auxilliary function is used to compute the laplacian. This function simply make the pointwise sum of the components u and v of two records.
trans Turing =
    x => let laplacian =
             neighborsfold(acc,{u= -4.0*x.u, v= -4.0*x.v},x)
         let du = k * (16.0 - x.u*x.v) + diffu*laplacian.u
         and dv = k * (x.u * x.v - x.v - x.c)
                  + diffv*laplacian.v
         in { u = x.u + dt*du,
              v = x.v + dt*dv,
              c = x.c }
For each point x of the torus, we compute the laplacian, which is simply a the sum of the neighbors minus the current state at x (times an adequate factor).The laplacian is then used to compute the increase in morphogens u and  v which are finally used to build the new state.

trans init =
    x => { u = u_steady,
           v = v_steady,
           c = beta_init + random(2*beta_rand) - beta_rand }

fun genere(n) =
      if n == 0
      then seq:()
      else (iota(NbCell_a, ():seq))::genere(n-1) fi

torus0 := init(genere(NbCell_b) following |b>, |a>);;

An auxilliary transformation init is used to define an initial state. Starting from a torus, each element is replaced by a record with fields u, v and c.

The function genere is used to build a nested list. We don't care the element in the nested lists (they are generated by the iota function that gives a sequence of integers)

The operator following is then used to convert this list into a torus and the transformation init is used to associate a record to each element of the toroidal grid.

output of the simulation data

itere := 1400;;
step  := max(1, itere/40);;
out   := "tmp.2turing-"

The simulation will evolve for 1400 time steps and a snapshot is taken every 40 iterations. The name of the output file include the size of the torus, the number of iteration done and the step.

fun Show(x) =
  if 0 == ('iteration -1) % step
  then show[t='iteration](x)
  else Turing(x)
The function Show is used to run the simulation. Parameter x is the current state. The variable 'iteration is a contextual variable, managed by the system and updated when the function Show is iterated (see below). If the current iteration is a multiple of step, then the current state is outputed to a file using the transformation show (Show and show are two different names because capital letters are meaningful in MGS).
Else, an evolution step is made by applying the transformation Turing.
trans show[t=0] =
    x => let pp = sequify(pos(x)) in
         let bx = take(pp, 0)
         and ax = take(pp, 1)
         in  out << "dataTuring[" << t
                 << "][" << ax
                 << "][" << bx
                 << "] = " << x.v
                 << ";\n";

The state of the torus is written to a file using the transformation show. This transformation takes an additional parameter t which represent the current time. Applied on the torus, each element is selected by the only rule of this transformation.
The position pos(x) of the element x is a word in the torus generator a and b. This word is transformed into a sequence of integers using the sequify function. The variables ax and bx represent the integer coordinates of this element in the circular grid along the axis a and b. And for each element, a line:
dataTuring[t][ax][bx] = v
(where v is the concentration of the morphogen of the element x at time t) is written in the output file.
This file is then processed by a Mathematica program to visualize the torus. This visualization program in Mathematica makes use of the Tuba Mathematica package written by Prof. M A Berger.
The function is iterated starting from the initial state torus0.

Back to top
MGS examples index
MGS home page

(currently under construction logo)
this site is under construction.
Pages started: May 2002. Last revision: 24 jully 2003.

Creative Commons License
Pictures, graphics and animations are licensed under a Creative Commons License.

English keywords for indexation: computer science, programming language, topological collections, transformation, declarative programming language, functional languages, simulation of biological processes, cell model, biological pathway, interaction network, gene regulation, signal transduction, morphogenesis, developmental biology, integrative simulation, biological organization, dynamical systems, dynamical structure, Gamma, CHAM, P system, L system, Paun, Lindenmayer, cellular automata, membrane computing, aqueous computing, artificial chemistry, GBF, Cayley graph, data fields, nested collections, rewriting, rule based programming, pattern-matching, intentional programming, compilation, interpretation, type, type inference, nested type, polytypism, catamorphism, static analysis, sequence, multiset, combinatorial algebraic topology, chain complex, chain group.