Turing systems as models of complex pattern formation

Turing systems as models of complex pattern formation

Teemu Leppänen^I; Mikko Karttunen^I; R.A. Barrio^I,II; Kimmo Kaski^I

^ILaboratory of Computational Engineering, Helsinki University of Technology, P.O. Box 9203, FIN-02015 HUT, Finland

^IIInstituto de Fisica, UNAM, Apartado Postal 20-364 01000 México, D.F., México

ABSTRACT

Half a century ago a reaction-diffusion system of two chemicals was introduced by Alan Turing to account for morphogenesis, i.e., the development of patterns, shapes and structures found in nature. Here we will discuss the formation of patterns and structures obtained through numerical simulation of the Turing mechanism in two and three dimensions. The forming patterns are found to depend strongly on the initial and boundary conditions as well as system parameters, showing a rich variety of patterns, e.g. stripes and spots (2D), and lamellae and spherical droplets (3D) arranged in structures of high symmetry, with or without defects or distortions.

1 Introduction

In the quest for understanding biological growth, it was Alan Turing (one of the key scientists of 20th century) who first demonstrated how a simple model system of coupled reaction-diffusion equations could give rise to spatial patterns in chemical concentrations through a process of chemical instability [1]. He showed that this kind of system may have a homogeneous stationary state which is unstable against perturbations, such that any random deviation from the stationary state leads through diffusion to a symmetry break. This process is called diffusion-driven instability. Since complex spatial patterns are commonly found in nature [2], for example in animal skins and also in some polymer systems [3], it is quite natural to think that such pattern formations could be caused by some general physico-chemical process. Inspired by this notion, Turing systems have been proposed to account for pattern formation in various biological systems [4], e.g. patterns on fish [5,6], butterflies [7] and lady beetles [8]. Furthermore, the first experimental evidence of a Turing structure was observed quite recently in a single-phase open chemical reactor [9].

The forms and variability of patterns generated by Turing systems have been studied by assuming inhomogeneous diffusion coefficients [10], or by introducing domain curvature [11] or growth [12]. In addition, the Turing systems have been investigated [10,13] from the point of view of symmetry breaking due to its potential relevance in understanding the mechanisms responsible for symmetries found in nature. Sofar most of the computational modelling work on Turing systems has been done in two-dimensions (2D), while three-dimensional (3D) systems have remained little studied, as pointed out by De Wit et al. [14]. They investigated the so called Brusselator reaction-diffusion model in 3D and demonstrated the possibility of additional high symmetry structures, i.e. body centered cubic (BCC) and hexagonally packed cylinders (HPC), and also a distorted lamellar structure which appears as a twist grain boundary. This latter structure embeds a minimal Scherk surface. Thus it is evident that 3D systems show a much richer set of morphologies than 2D systems.

Recently we have also studied the effect of dimensionality by simulating a three-dimensional Turing system, which clearly shows an increased complexity in pattern formation [15]. While in 2D we obtain structures with spots organized hexagonally or with stripes, or some type of labyrinthine patterns, in 3D complex lamellar, spherical droplet patterns or their combinations are seen. We have also observed that initial conditions play a crucial role as to which stationary or long living pattern the system evolves. For example with certain initial introduction of one of the morphogens the system seems to evolve to a hexagonally arranged spherical droplet pattern or to a perfect lamellar structure, depending on the parameters of the model system. Also we have studied situations in which the choice of system parameters indicate competition between spherical droplet and lamellar structures, and they give rise to complex hybrid patterns. In addition we have studied connectivity and clustering in the patterns [16] and the effect of noise on Turing structures [17].

These findings are in line with the observations by De Wit et al. [14] that 3D reaction-diffusion systems can show a wide variety of possible perfect and distorted patterns. However, much remains to be explored and studied in more detail. Thus the purpose of this study is to explore the morphologies of patterns within the framework of generic Turing system, in which it is possible to control pattern formation with the parameters of the system. Also the effects of initial and boundary conditions are going to be investigated.

In the next section we briefly present the generic Turing model and discuss its characteristics from the point of view of mode instability. Then we present the results of numerical simulations in both 2D and 3D systems. Finally we draw conclusions.

2 Generic Turing Model

As mentioned above a Turing system describes temporal behaviour of the concentrations of two reacting and diffusing chemicals or morphogens. This can be represented in general by the following coupled reaction-diffusion equations

where and are the morphogen concentrations, and D_U and D_V the corresponding diffusion coefficients setting the time scales for diffusion. The reaction kinetics is described by the two non-linear functions f and g.

In this paper we use the generic Turing model introduced by Barrio et al. [6], in which the reaction kinetics was developed by Taylor expanding the non-linear functions around a stationary solution (U_c,V_c), defined by f(U_c,V_c) = g(U_c,V_c) = 0. If terms of the fourth order and higher are neglected, the reaction-diffusion equations can be written as

where u = U – U_c and v = V – V_c. The parameters r₁ and r₂ set the amplitudes of the non-linear cubic and quadratic terms, respectively, serving as control for favouring a certain pattern formation. In fact by gradually changing the non-linear parameters we observe a transition from spotty (2D) or spherical droplet (3D) patterns to striped (2D) or lamellar (3D) patterns. In the equation the quantity D is the ratio of the diffusion coefficients of the two chemicals, and d acts as a scaling factor fixing the size of the system. In the analysis of the model system we can without loss of generality fix a = –g in order to have only one stationary state in the system. For details about the modes of instability and the linear stability analysis of the model we refer the reader to Barrio et al. [6,15]. Although the conditions for diffusion-driven instability to occur are widely known [4], it should be mentioned that setting D ¹ 1 is a necessary but not sufficient condition for the diffusion-driven instability in 2D and 3D systems [18], and that by restricting the parameter selection such that a Î (0,1) and b Î (–1,0) one is left with only two additional conditions, namely.

Now from the linear stability analysis [6,15] we obtain the following dispersion relation

Based on this equation one can find a region in k-space with positive growth rate, i.e., the eigenmodes u = u₀e^{l t} and v = v₀e^{l t} with eigenvalues l(k) > 0. In addition, we can analytically derive the modulus of the critical wave vector

For the numerical integration of Eq. (2) the (finite) model system needs to be discretized. Then in a three-dimensional cubic system, the wave number is of the form

where L is the system size and n_x, n_y, n_z are the wave number indices (Note that in a two-dimensional system n_z = 0). By adjusting the parameters and allowing only a few unstable modes, one can end up with several different parameter sets. As in our earlier work [15] we chose the parameters D = 0.516, a = –g = 0.899, b = –0.91 and d = 2 corresponding to a critical wave vector k_c = 0.45, and D = 0.122, a = 0.398, b = –0.4 and d = 2 corresponding to k_c = 0.84.

3 Simulation Results

For numerical simulations the spatial dimensions of the system were first discretized into a square or cubic lattice mesh with lattice constants dx = dy = dz = 1.0. The equations of motion, Eq. (2) were iterated in time using the Euler scheme with time step dt = 0.05. Both zero-flux and periodic boundary conditions were used.

In Fig. 1 we present the results for the 2D Turing system. In the left column we show the simulation results of Turing patterns for four sets of non-linear parameters (r₁ and r₂ in Eq. (2)). By slowly changing the strength of the quadratic nonlinear interaction term (r₂) one can observe transition from stripes to spots. Our preliminary results suggest that in the limit of infinite simulation time the transition between stripes and spots could become discontinuous. However, with a finite simulation time one always observes intermediate transient states (Figs. 3B and C) for fixed critical amplitudes of the nonlinearities. In the proximity of these critical values one observes critical slowing down of the dynamics: the closer one is to the critical values of the nonlinear terms the longer it takes for the system to reach either regular striped or spotty morphology.

One can get further insight into the transition by taking the discrete Fourier transform of the concentration data into the wave vector space, , i.e.,

This approach has been used earlier, e.g. for reaction-diffusion systems [19], Turing patterns [20], and to characterize the evolution of patterns [21]. The quantity corresponds to a diffraction pattern. The results are shown in the right column of Fig. 1, which shows a sequence of diffraction patterns corresponding to the concentration space Turing patterns on the left.

In Fig. 1A, the diffraction intensity is to k_x-direction because the stripes are aligned in y-direction. The distance of these diffraction peaks from the origin gives the length of the wave vector of the unstable mode, while its width in the perpendicular direction describes deviations of stripes from the principal directions due to stripe being tilted and bent. In Fig. 1B the diffraction intensity becomes more spread around the diagonal in (k_x k_y)-space, as a result of nucleating defects and spots. After that in Fig. 1C the diffraction intensity starts splitting into separate peaks due to more spots forming, then developing into six separate equidistant (from the origin and from each other) intensity peaks as evident in Fig. 1D. This hexagonal symmetry in the reciprocal space is because the system evolves towards regular spotty pattern with predominantly triangular symmetry in the real concentration space. However, these diffraction peaks are somewhat spread, which indicates that the triangular symmetry is not perfect over the whole system.

In three dimensions the pattern selection is not as straightforward as in two dimensions. Instead of perfectly organized spherical structures and planar lamellae, one may obtain structures that establish order in more general ways. In the case of spherical droplet structures, the droplets are not always organized on a perfect lattice. The lamellae are rarely planar, but more complex minimal surfaces or twisted grain boundaries arise, c.f. [14]. In addition, the intermediate structures are not separate domains of spherical droplet structures and lamellae as spots and stripes in 2D, but a variety of stable structures is obtained. For typical 3D Turing structures we refer the reader to [15].

Some of the 3D structures generated by the generic Turing model resemble those that have been predicted and simulated in diblock copolymer melts [22,23] and also those obtained by De Wit et al. [14] using the Brusselator model. In Fig. 2 we show specifically three Turing structures that are also observed in polymer melt systems: spherical droplets organized into BCC lattice, planar lamellar structure, and HPC structure. These structures were obtained by numerical simulation of the generic Turing model of Eq. (2) after 400000 iterations in a 50 × 50 × 50 system. The initial conditions and the system parameters had to be adjusted carefully to obtain these structures. However, the structures proved to be stable against noise once established.

For the BCC structure in Fig. 2 the initial concentration field had the BCC symmetry, which made the system find the BCC structure of larger spherical droplets. This structure we could not, however, get by starting from a random initial configuration of morphogens. For the planar lamellar structure, chemical U was introduced only in the midplane of the system and V was distributed uniformly over the whole system. For the same set of parameters but starting from random introduction of both chemicals we find a lamellar minimal surface, which is characterised to have both positive and negative local principal curvature. As for the HPC structure, it was obtained starting the simulation from random introduction of both chemicals. The used parameter values can be found from the caption of Fig. 2.

There is an optimal minimal surface for connecting perpendicular planes such that the local principal curvatures (c₁ and c₂) have opposite signs yielding zero local mean curvature (H = (c₁ + c₂)/2 = 0). This surface is called Scherk surface, which is a well-known minimal surface solution for twisted lamellar surfaces. The stability of a Scherk surface has previously been observed in the case of the Brusselator model [14] by initializing the minimal surface into the system. In Fig. 3 we show the final structure obtained as a result of a simulation run of a system initializing one of the chemicals into two planar domains set perpendicularly. Many twisted connections of lamellar surfaces are found. Based on the close-up of one of the connections we can observe that also spontaneous formation of the Scherk surface is possible in the generic Turing system. The parameters in this case were r₁ = 3.5 and r₂ = 0 for mode k_c = 0.84.

4 Summary

In summary we have studied the pattern formation in a generic Turing model in two and three dimensions and under different initial and boundary condition. The generic Turing model is computationally very interesting and rich because one can change the nonlinear parameters to control the pattern formation and explore the space of pattern the system produces. In 2D we have observed striped and spotty patterns and a transition between them when the nonlinear parameters are changed. This transition turns out to be very sharp, or possibly even discontinuous at the limit of infinite time. In 3D the number of possible stationary patterns turned out to be much larger than in 2D. We have observed spherical droplet pattern (organized either somewhat irregularly or in hexagonal arrangement), regular BCC pattern of spherical droplets, regular HPC structure, perfect lamellar structure, distorted lamellar structures which carry characteristics of minimal surface, and Scherk surface structure of minimal surface solution for twisted lamellar surfaces.

To conclude it is evident that Turing systems have very general ways of establishing order since we have observed spontaneous formation of a wide variety of stable minimal surfaces in addition to the Scherk surface. All these surfaces have the characteristic length fixed by the parameters and predicted by linear analysis, but the morphologies can be very different. This interesting problem will be addressed in more detail in future research.

5 Acknowledgements

One of us (R. A. B.) wishes to thank the Laboratory of Computational Engineering at Helsinki University of Technology for their hospitality. This work has been supported by the Finnish Academy of Science and Letters (T.L.) and the Academy of Finland through its Centre of Excellence Program (T. L. and K. K.).

Received on 9 August, 2003

• [1] A.M. Turing, Phil. Trans. R. Soc. Lond. B237, 37-72 (1952).
• [2] P. Ball, The self-made tapestry: pattern formation in nature, Oxford Univ. Press, Oxford (2001).
• [3] S.O. Kim et al., Nature 424, 411-414 (2003).
• [4] J.D. Murray, Mathematical Biology, 2nd. ed., (Springer Verlag, Berlin 1993).
• [5] S. Kondo, and R. Asai, Nature 376, 678 (1995).
• [6] R.A. Barrio, C. Varea, J.L. Aragón, and P.K. Maini, Bull. Math. Biol. 61, 483 (1999).
• [7] T. Sekimura, A. Madzvamuse, A.J. Wathen, and P.K. Maini, Proceedings Royal Society London B 267, 851 (2000).
• [8] S.S. Liaw, C.C. Yang, R.T. Liu, and J.T. Hong, Phys. Rev. E 64, 041909 (2002).
• [9] V. Castets, E. Dulos, J. Boissonade and P. de Kepper, Phys. Rev. Lett. 64, 2953 (1990).
• [10] R.A. Barrio, J.L. Aragón, M. Torres, and P.K. Maini, Physica D 168-169, 61 (2002).
• [11] C. Varea, J.L. Aragón, and R.A. Barrio, Phys. Rev. E 60, 4588 (1999).
• [12] C. Varea, J.L. Aragón, and R.A. Barrio, Phys. Rev. E 56, 1250 (1997).
• [13] J.L. Aragon, M. Torres, D. Gil, R.A. Barrio, and P.K. Maini, Phys. Rev. E, 65, 051913 (2002).
• [14] A. de Wit, P. Borckmans and G. Dewel, Proc. Natl. Acad. Sci. USA 94, 12765 (1997).
• [15] T. Leppänen, M. Karttunen, K. Kaski, R.A. Barrio, and L. Zhang, Physica D 168-169, 35 (2002).
• [16] T. Leppänen, M. Karttunen, R.A. Barrio, and K. Kaski, ''Morphological transitions and bistability in Turing systems", submitted 2003. cond-mat/0302101.
• [17] T. Leppänen, M. Karttunen, R.A. Barrio, and K. Kaski, Prog. Theor. Phys. 150 (2003).
• [18] J.A. Vastano, J.E. Pearson, W. Horsthemke, and H.L. Swinney, Phys. Lett. A 124, 6 (1987).
• [19] M. Karttunen, N. Provatas, T. Ala-Nissila, and M. Grant, J. Stat. Phys. 90, 1401 (1998).
• [20] L. Yang, M. Dolnik, A. M. Zhabotinsky, and I. Epstein, Phys. Rev. Lett. 88, 208303 (2002).
• [21] M. Karttunen, M. Haataja, K. R. Elder, and M. Grant, Phys. Rev. Lett. 83, 3518 (1999).
• [22] M.W. Matsen and F.S. Bates, Macromolecules 29, 1091 (1996).
• [23] R.D. Groot and T.J. Madden, J. Chem. Phys. 108, 8713 (1998).

Publication Dates

• Publication in this collection
  01 Sept 2004
• Date of issue
  June 2004

History

• Accepted
  09 Aug 2003
• Received
  09 Aug 2003