The Mathematics of Patterning
One of the most important mathematical models in developmental biology has been that formulated by Alan Turing (1952), one of the founders of computer science (and the mathematician who cracked the German “Enigma” code during World War II). He proposed a model wherein two homogeneously distributed solutions would interact to produce stable patterns during morphogenesis. These patterns would represent regional differences in the concentrations of the two substances. Their interactions would produce an ordered structure out of random chaos.
Turing's reaction-diffusion model involves two substances. One of them, substance S, inhibits the production of the other, substance P. Substance P promotes the production of more substance P as well as more substance S. Turing's mathematics show that if S diffuses more readily than P, sharp waves of concentration differences will be generated for substance P (Figure 1). These waves have been observed in certain chemical reactions (Prigogine and Nicolis 1967; Winfree 1974).
Figure 1 Reaction-diffusion (Turing model) system of pattern generation. Generation of periodic spatial heterogeneity can come about spontaneously when two reactants, S and P, are mixed together under the conditions that S inhibits P, P catalyzes production of both S and P, and S diffuses faster than P. (A) The conditions of the reaction-diffusion system yielding a peak of P and a lower peak of S at the same place. (B) The distribution of the reactants is initially random, and their concentrations fluctuate over a given average. As P increases locally, it produces more S, which diffuses to inhibit more peaks of P from forming in the vicinity of its production. The result is a series of P peaks (“standing waves”) at regular intervals.
The reaction-diffusion model predicts alternating areas of high and low concentrations of some substance. When the concentration of such a substance is above a certain threshold level, a cell (or group of cells) may be instructed to differentiate in a certain way. An important feature of Turing's model is that particular chemical wavelengths will be amplified while all others will be suppressed. As local concentrations of P increase, the values of S form a peak centering on the P peak, but becoming broader and shallower because of S's more rapid diffusion. These S peaks inhibit other P peaks from forming. But which of the many P peaks will survive? That depends on the size and shape of the tissues in which the oscillating reaction is occurring. (This pattern is analogous to the harmonics of vibrating strings, as in a guitar. Only certain resonance vibrations are permitted, based on the boundaries of the string.)
The mathematics describing which particular wavelengths are selected consist of complex polynomial equations. Such functions have been used to model the spiral patterning of slime molds, the polar organization of the limb, and the pigment patterns of mammals, fish, and snails (Figures 2 and 3; Kondo and Asai 1995; Meinhardt 1998). A computer simulation based on a Turing reaction-diffusion system can successfully predict such patterns, given the starting shapes and sizes of the elements involved.
Figure 2 A photograph of the snail Oliva porphyria (left), and a computer model of the same snail (right) in which the growth parameters of the shell and its pigmentation pattern were both mathematically generated. (From Meinhardt 1998; computer image courtesy of D. Fowler, P. Prusinkiewicz, and H. Meinhardt.)
Figure 3 Pigment patterns of zebrafish homozygous for the wild-type allele (A) and for three different mutant alleles (B–D) of the leopard gene. Computer simulations of the pigment patterns are shown in the bottom row. Changing a single parameter of the reaction-diffusion equation results in changes in the pattern. As the values of this parameter become larger, the stripes break into spots and the spots get smaller and less dense. (From Asai et al. 1999; photographs courtesy of S. Kondo.)
One way to search for the chemicals predicted by Turing's model is to find genetic mutations in which the ordered structure of a pattern has been altered. The wild-type alleles of these genes may be responsible for generating the normal pattern. Such a candidate is the leopard gene of zebrafish (Asai et al. 1999). Zebrafish usually have five parallel stripes along their flanks. However, in the different mutations, the stripes are broken into spots of different sizes and densities. Figure 3 shows fish homozygous for four different alleles of the leopard gene. If the leopard gene encodes an enzyme that catalyzes one of the reactions of the reaction-diffusion system, the different mutations of this gene may change the kinetics of synthesis or degradation. Indeed, all the mutant patterns (and those of their heterozygotes) can be computer-generated by changing a single parameter in the reaction-diffusion equation. The cloning of this gene should enable further cooperation between theoretical biology and developmental anatomy.
The Development of Zebra Striping Patterns
Few patterns are more obvious than the alternating black-and-white stripes of the zebra (Figure 4). There are actually three extant species of zebra, and each has a different pattern of stripes. The imperial zebra (Equus grevyi) has some eighty stripes perpendicular to the long axis of its body. The common zebra (E. burchelli) has 26 wide caudal stripes, some of which extend towards the belly in the rear of the animal. The mountain zebra (E. zebra) has some 55 stripes, with three horizontal bands near the hindlegs. Each of these three species are members of the horse genus and can interbreed among themselves and other horses to produce infertile offspring.
Figure 4 Different species of zebra: (a) The imperial zebra (Equus grevyi), (b) the mountain zebra (Equus zebra), (c) the common zebra (Equus burchelli), and (d) the quagga (Equus quagga). From a drawing. Photograph from Bard (1977), with permission of the author.
How did the zebra get its stripes?
Ultimate (evolutionary) mechanism
It is generally believed that zebras are dark animals with white stripes where the pigmentation is inhibited. The pigment of the hair is found solely in the hair and not in the skin. The reasons for thinking that they were originally pigmented animals are that (1) white horses would not survive well in the African plains or forests; (2) there used to be a fourth species of zebra, the quagga (which was overeaten to extinction in the 1800s). The quagga had the zebra striping pattern in the front of the animal, but had a dark rump; (3) when the region between the pigmented bands becomes too wide, secondary stripes emerge, as if suppression was weakening.
Zebra stripes have often been thought to be an adaptation that prevents zebras from being seen by predators such as lions or hyenas. (This hypothesis goes back at least to Rudyard Kipling [1908]). The alternating stripes obscure the outline of the zebra. This may serve as camouflage, allowing the zebra to blend in with its backgound (Thayer 1909; Marler and Hamilton 1968) and/or it may serve to confuse a predator as to the distance of the fleeing animal (Cott 1957; Kruuk 1972). However, neither of these hypotheses can be easily confirmed. A different hypothesis (Waage 1981) contends that the stripes serve to obliterate a large single-colored region that is favored by biting insects such as the tsetse fly. These flies prefer large, dark, moving animals (Vale 1974).
How did the zebra get its stripes?
Proximate (developmental) mechanism
Jonathan Bard of Edinburgh has hypothesized a mechanism for the production of zebra stripes in the three species of extant zebras (1977, 1981). His model claims that while neural crest cells begin migration at week two of gestation (in the horse), the zebra striping patterns are generated between weeks three and five, depending upon the species. Moreover, Bard asserts that the three patterns of striping are precisely those predicted if the original pattern was the same in each zebra, but was established at different times within this three week period. In the case of the imperial zebra, all the stripes are perpendicular to the dorsal axis, but are thicker towards the neck. This would be expected if the striping pattern originated at week five (Figure 5a). At week five, most of the differential body growth has ceased, except for the neck region, which becomes extended, and the rump, which is slightly shortened. Thus, if the stripes were formed at week five, they should all be parallel, but slightly wider at the neck and narrower at the rump.
Figure 5 Bard's hypothesis for the generation of stripes in three species of zebras. The spacing and size of the stripes are the same. What differs is the time at which the stripes were generated. If generated during week 3, the stripes begin perpendicular to the anterior-posterior body axis, but become parallel to this axis in the rump, since the rear of the zebra is still growing. This generates the pattern of common zebra. If the striping pattern is generated during week 4, most of the rump has grown, and the hind stripes are more perpedicular to the body axis. This generates the pattern seen in the mountain zebra. If the striping pattern is generated during week 5, there is space for many more stripes, all of which are perpendicular to the body axis. This generates the striping pattern of the imperial zebra. (After Bard 1977.)
The stripes of the mountain zebra probably form towards the end of week four. If the stripes were originally parallel, those in the rear of the embryo would be pulled back towards the rump by the growth of the hindparts of the horse (Figure 5b). Similarly, if the stripes of the common zebra were generated during the third week of zebra gestation, the differential growth rate of the rump between weeks three and four would also pull the stripes posteriorly (Figure 5c).
Bard's hypothesis that all the stripes originally are the same width and are generated at different times in the three species also explains the numbers of stripes in each species. The common zebra has 26 stripes per side, and the 3-week Equus embryo is generally 11 mm long. This gives a spacing of about 0.42 mm per stripe. If the 43 stripes of the mountain zebra were generated in the 17 mm embryo of the 3.75 week zebra, the spacing is also 0.40 mm per stripe. At week 5, the embryo is 32 mm long, and the 80 stripes would yield the spacing of 0.40 mm per stripe. Therefore, the striping patterns of the common zebra, mountain zebra, and imperial zebra can be explained if the stripes are generated 0.4 mm apart in the 3-, 4-, and 5-week embryos, respectively.
It is not known how the pattern is initiated or what activators or inhibitors are being generated. It is difficult to imagine how such a pattern can be generated by preformed maternal instructions, responses to gradients, or regional inductions. It has been proposed that the Turing reaction-diffusion models could produce these alternative pigmented and non-pigmented bands. Murray (1981) has shown that the chevrons at the base of the zebra's limbs are the shape expected by the overlapping of two Turing-type reaction-diffusion systems.
Mathematical Background of Pattern Formation
The Mathematics of Reaction-Diffusion Mechanisms
Probably the best introduction to the mathematics of pattern formation in developmental biology is on Hans Meinhardt’s website, a remarkable and beautiful resource.
The classic paper Gierer, A. and Meinhardt, H. 1972. A theory of biological pattern formation, Kybernetic 12: 30-39 is a critically important paper in modeling developmental phenomena (PDF format).
I. The Basic Activator-Inhibitor System
(Based on Meinhardt, H. 1998. The Algorithmic Beauty of Sea Shells, Second ed., Springer-Verlag, New York.)
The mathematics of the reaction-diffusion system involves the interactions between an autocatalytic activator (P or a) and its antagonist, the inhibitor (S or b). The partial differential equations relate the concentration change per time of both these substances as a function of their concentrations:
Activator:
Inhibitor:
Where:
t is the time, x is the spatial coordinate, Da and Db are the diffusion coefficients for the activator and the inhibitor, respectively (and the diffusor inhibits faster than the activator).
The terms can then be understood as such:
sa2/b is the production rate of the activator. The activator has an autocatalytic rate such that the more a is present, the more a is made. The product formation is slowed down by the presence of b. The source density s describes the ability of the agents to perform autocatalysis.
-raa is the rate at which the activator is removed. As a general rule, the rate at which molecules disappear is proportional to the number of molecules present.
Da(Δ2a/Δx2) is the exchange by diffusion. It is proportional to the second derivative for the following reason: The net exchange of molecules by diffusion is 0 if all cells have the same concentration. The net exchange is also 0 if a constant concentration difference exists between neighboring cells (as in a linear concentration gradient). In that case, each cell loses and receives the same amount of substance. Thus, this term represents not the change of concentrations, but the change in the change of concentrations in space.
ba is the basic activator production. A small activator-independent production can initiate the system at low activator concentration levels. This is the term that was altered in the equations to produce Figure 20.22 in the textbook. It is required for sustained oscillations during growth.
bb is the basic inhibitor production. This function can change the stable state at which the system arrives. It is also important in producing "travelling waves."
II. More Detailed Mathematics of the Reaction-Diffusion System
One of the most interesting and widely studied topics in modern theoretical biology is the structure formation from a more or less homogeneous egg. Here some global field models of such processes are represented. A characteristic feature of these models is that they consider autonomous growth (autocatalysis) coupled to dissipative processes, such as diffusion.
In some cases one can consider these models as Turing's systems of the first kind, that is represented by differential equation:
More rarely Turing's systems of second kind are used:
As was noticed, before the stable pattern can be generated two conditions have to be fulfilled: a local deviation from an average concentration of a pattern forming substance should further increase (1), and this increase should not go to infinity (2). In supposition that increase in one part of the field is necessarily connected with a decrease in another part of it, (i.e., that the total amount of substances is roughly conserved), the process of the pattern formation should reach the stable steady state.
Thus the mechanism underlying pattern formation should be similar to local autocatalysis with strong positive feedback and lateral inhibition. We separately consider the model with diffusion (class 1 models) and without it (class 2 models). Here we will only discuss class 1 models.
Models with diffusion (Class 1 models)
Here we represent the application of principles of autocatalysis and lateral inhibition to biochemical reactions with diffusion. Let us consider a substance a, called an activator, which stimulates its own production (autocatalysis) and the production of its antagonist i, called an inhibitor. To carry out the necessary long-range inhibition, the inhibitor must diffuse more rapidly. In an extended field of cells, a homogeneous distribution of these substances is unstable, since any small local elevation of the activator concentration, resulting perhaps from random fluctuations, will be amplified by the activator autocatalysis. The inhibitor, which is produced in response to the increase in activator production, cannot halt the locally increased activator production, since it diffuses quickly into the surrounding tissue and suppresses activator production outside the activated center. Thus, the locally increased activator concentration will increase further, and with increasing concentration, the maximum becomes narrower and narrower until some limiting factor comes into play. For instance, the loss of activator from the narrow peak by diffusion become sequel to the net production.
A stable activator and inhibitor profile is ultimately obtained, although both the substances continue to be made, to diffuse, and to be broken down. Such a simple system of two interacting substances is, therefore, able to produce a stable, strongly patterned distribution from a nearly homogeneous initial distribution, as it occurs in biological pattern formation. The general representation of such a system is:
The first term in the right part represents the diffusion process, where Da and Di are diffusion constants for a and i correspondingly; the second term represents other processes such as production and decay, they may be written as:
here if c1=0, then c2=1 and vice versa. In these equations, production is represented by the first and third terms, the second term represents the decay. One of the most known and frequently used types of systems of such a kind can be derived if to set c1=1 (correspondingly c2=0), k=r1=r2=0, namely
Resulting distribution may be monotonous or periodical and changes as a function of parameters. The pattern formed in this way may be stable or oscillating with the time. The positioning of high concentrations is produced by small internal or external asymmetries or by local disturbance. This local high concentration can serve as a signaling system, for instance, to initiate head formation in hydra. A pattern formed in this way has strong self-regulatory properties. Let us also consider a special case of such a systems with a feedback loop.
To do this set c2=0, then
The small basic (activator-independent) activator production r1 can initiate the autocatalysis in areas of low activator concentration. As we will see, this term is important if new centers have to arise during growth or regeneration. In contrast, the basic inhibitor production can suppress the appearance of secondary maxima, a feature which is important if an ordered sequence of structures is to be specified by positional information. If the activator production saturates at a high concentration due to term 1/(1+k × a × a), the activated area is self-regulated.
Other molecular realization of autocatalysis and lateral inhibition principle may be the processes in which the inhibition effect is realized by depletion of a substrate consumed in autocatalysis. From a mathematical point of view this means:
Also one has to remember Sel'kov's model in developmental biology:
where n is known as Hill's number.
In this model, one assumes that the activator reproduction is compensated by self-regulated reproduction of the inhibitor. The stable pattern formation is also possible in this case.
For a deeper mathematical understanding, see Oscillating Reaction-Diffusion Spots.
III. Reaction-Diffusion in Dictyostelium
As documented on pages 676-679 of the textbook (and in Vade Mecum2), the slime mold Dictyostelium discoideum is a part-time multicellular organism. Most of the time, Dictyostelium is an amoeba, crawling around the forest floors of North America, eating bacteria. However, when the food supply is insufficient, they aggregate by the thousands. The molecule signaling the aggregation is cyclic AMP (cAMP). This molecule appears to be autocatalytic, in that it tells the cell receiving it to make more of it. In fact, when cAMP is sensed by the responding cell, it takes only 30 seconds for that cell to start making more cAMP (Gerisch 1968; Devreotes 1989; Goldbeter 1996). However, that same autocatalysis also signals inhibition as the receptor molecules become phosphorylated. Once phosphorylated, the receptor cannot bind cAMP. The refractor period lasts about five minutes and then the receptor can function again. If Dictyostelium amoebae are packed together densely enough, the activation spreads in a wave-like motion (see Vade Mecum2). Those cells that initiated the release of cAMP will become the centers of aggregation. Cells will move towards these centers since one side of their cytoplasm receives the signal before the other. In this way, the wave is propagated from the center to the periphery of the aggregate. The movement of the slug also seems to be under the regulation of the waves of cAMP.
Literature Cited
Devreotes, P. N. 1989. Dictyostelium discoideum: A model system for cell-cell interactions in development. Science 245: 1054–1058.
Gerisch, G. 1968. Cell aggregation and differentiation in Dictyostelium. Curr. Top. Devel. Biol. 3: 157–232.
Goldbeter, A. 1996. Biochemical Oscillations and Cellular Rhythm: The Molecular Bases of Periodic and Chaotic Behavior. Cambridge University Press, Cambridge.
Literature Cited
Asai, R., E. Taguchi, Y. Kume, Saito M. and S. Kondo. 1999. Zebrafish Leopard gene as a component of the putative reaction-diffusion system. Mech. Dev. 89: 87–92.
Bard, J. B. L. 1977. A unity underlying the different zebra striping patterns. J. Zool. (London) 183: 527–539.
Bard, J. B. L. 1981. A model for generating aspects of zebra and other mammalian coat patterns. J. Theoret. Biol. 19: 363–385.
Cott, H. B. 1957. Adaptive Colouration in Animals. John Dickens. Northampton.
Kipling, R. 1908. Just So Stories. Macmillan, London.
Kondo, S. and R. Asai. 1995. A reaction-diffusion wave on the skin of the marine angelfish Pomacanthus. Nature. 376: 765–768.
Kruuk, H. 1972. The Spotted Hyena. University of Chicago Press, Chicago.
Marler, P. and W. J. Hamilton. 1968. Mechanisms of Animal Behavior. Wiley, New York.
Meinhardt, M. 1998. The Algorhythmic Beauty of Sea Shells. Springer: Berlin.
Murray, J. D. 1981. A pre-pattern formation mechanism for animal coat markings. J. Theoret. Biol. 88: 161–199.
Prigogine I. and G. Nicolis. 1967. On symmetry-breaking instabilities in dissipative systems. J. Chem. Phys. 46: 3542–3550.
Thayer, A. H. 1909. Concealing Coloration in the Animal Kingdom. Macmillan, New York.
Turing, A. M. 1952. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. 237: 37–72.
Vale, G. A. 1974. The response of tsetse flies (Diptera, Glossinidae) to mobile and stationary baits. Bull. Entom. Res. 64: 545–588.
Waage, J. K. 1981. How the zebra got its stripes: biting flies as selective agents in the evolution of zebra coloration. J. Entom. Soc. South Afric. 44: 351–358.
Winfree, A. T. 1974. Rotating chemical reactions. Sci. Am. 230(6): 82–95.