UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Mathematical models for biological pattern formation in two dimensions Lyons, Michael J. 1992

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1992_spring_lyons_michael.pdf [ 4.97MB ]
JSON: 831-1.0085628.json
JSON-LD: 831-1.0085628-ld.json
RDF/XML (Pretty): 831-1.0085628-rdf.xml
RDF/JSON: 831-1.0085628-rdf.json
Turtle: 831-1.0085628-turtle.txt
N-Triples: 831-1.0085628-rdf-ntriples.txt
Original Record: 831-1.0085628-source.json
Full Text

Full Text

MATHEMATICAL MODELS FOR BIOLOGICAL PATTERNFORMATION IN TWO DIMENSIONSByMichael J. LyonsB. Sc. (Physics) McGill UniversityM. Sc. (Physics) The University of British ColumbiaA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESPHYSICSWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1991© Michael J. Lyons, 1991Signature(s) removed to protect privacyIn presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)_______________________________Department of PscThe University of British ColumbiaVancouver, CanadaDate______________DE-6 (2/88)Signature(s) removed to protect privacyAbstractThe general problem of pattern formation in biological development is described. Onepossible type of general solution, reaction-diffusion theory, is outlined, together with thedifficulties involved in definite experimental confirmation or rejection of it. The methodsfor linear analysis of reaction-diffusion equations are presented. An approximate analytical technique for solving nonlinear reaction-diffusion equations close to the onset of thepattern-forming instability is discussed. The technique, known as the adiabatic approximation, is applied to the problem of stripe formation in two spatial dimensions. Thisanalytical result demonstrates the existence of a class of pattern-forming models whichpreferentially select striped patterns. This class of models is defined by the requirementthat the reaction terms contain only odd powers of the departures from the homogeneoussteady state.The nonlinear properties of reaction-diffusion models in two spatial dimensions havebeen extensively studied numerically using finite-difference methods. These computations support the results of the analysis and extend their validity to a wider range ofparameter values. Additional nonlinear effects, for example elimination of pattern defects, orientation of pattern using monotonic gradients, and the stripe/spot transition,have been investigated computationally.Two specific instances of biological stripe formation are considered: the expressionof segmentation genes in the early embryogenesis of Drosophila melanogaster, and thedevelopment of ocular dominance domains in the primary visual cortex of some highervertebrates. Mechanisms by which the segmentation patterns could arise by reactiondiffusion are discussed. Existing models for ocular dominance, which take the form ofiiintegro-differential equations, are analyzed using the adiabatic approximation. A fundamental connection between a symmetry in the visual system and stripe selection isestablished.111Table of ContentsAbstract iiList of Tables viiList of Figures viiiAcknowledgement xi1 Introduction 11.1 Biological Pattern Formation 21.2 Snmmary of Thesis 92 Reaction-Diffusion Theory 122.1 Reaction-Diffnsion Systems 122.2 Two Component Reaction-Diffnsion Mechanisms. 142.3 Experimental Evidence for Tnring Strnctllres . . . 202.4 Tnring Instability in the Ring of Discrete Cells 232.5 Conditions for Pattern Formation and the Chemical Wavelength 262.5.1 A4k. Space 322.6 Reaction-Diffnsion Mechanisms 362.6.1 Brnsselator 372.6.2 Hyperchirality Model 413 Stripe Formation in Two Dimensions by Reaction-Diffusion 45iv3.1 Reaction-Diffusion Theory in Two Dimensions 453.2 A Class of Reaction-Diffusion Mechanisms Which Preferentially SelectStripes 504 Numerical Solution of Reaction-Diffusion Equations 644.1 Finite-Difference Methods for the Numerical Solution of Reaction-DiffusionEquations 644.2 Forward Time Differencing - the Explicit Method 664.3 Implicit Finite-Difference Methods 714.4 Two Spatial Dimensions 724.5 Simulation of Two-Dimensional Pattern Formation 744.6 Description of Numerical Results 764.6.1 Animation of Pattern Formation Simulations 845 Biological Examples: Drosophila Segmentation and Ocular Dominance 1045.1 Background Material on Drosophila 1065.2 Reaction-Diffusion Models of Drosophila Segmentation 1095.3 Ocular Dominance Stripes 1135.4 Models for Ocular Dominance Pattern Formation 1185.5 Nonlinear Analysis 1196 Conclusion 126Appendices 129A Turing’s Two Cell Reactor 129A.1 Solution of Turing’s Example 129VB An Exact Result for Coupled Brusselators 135Bibliography 139viList of Tables2.1 Turing’s conditions in terms of dimensionless parameters 33viiList of Figures1.1 Citations of ‘The Chemical Basis of Morphogenesis’ per year since 1961. 31.2 Development of a Turing structure in an activator-inhibitor reaction-diffusionsystem 112.1 Qualitative stability diagram for two first order differential equations. . 182.2 Castets’ gel strip reactor 212.3 Microtubular dissipative structures 222.4 Lacalli’s k’k parameter space 342.5 Dispersion relation for reaction-diffusion systems 353.6 Dispersion relation and spectrum of modes on a narrow domain 483.7 Dispersion relation and spectrum of modes on a wider domain 483.8 Turing structure on a narrow domain 493.9 Turing structure on a wider domain 493.10 Arrangement of wavevectors in k-space 584.11 Finite difference lattice 674.12 Discretization of the Laplacian operator in two spatial dimensions . . . 734.13 The two basic types of singularities for ridge systems 804.14 Examples of computed reaction-diffusion patterns displaying the two basictypes of singularities 814.15 X and Y morphogen profiles in a brusselator computation 864.16 Figure 4.15 continued 87viii4.17 X and Y morphogen profiles in a computation with the hyperchirality model 884.18 Figure 4.17 continued 894.19 Comparison of brusselator and generic stripe making model 904.20 Effect of quadratic nonlinearity on stripe formation 914.21 Striped patterns observed with reaction-diffusion 924.22 Effect of a monotonic gradient on pattern selection 934.23 Calculation as in figure 4.22 but without gradient 944.24 Effect of a monotonic gradient on pattern selection far from the diffusiveinstability 954.25 Calculation as in figure 4.24 but without gradient 964.26 Effect of a monotonic gradient on pattern selection in the brusselator. . 974.27 Early stages of pattern formation by reaction-diffusion 984.28 Continuation of figure 4.27 994.29 Continuation of figure 4.27 1004.30 Later stages of pattern formation by reaction-diffusion 1014.31 Continuation of figure 4.30 1024.32 Continuation of figure 4.30 1035.33 The Drosophila body plan 107.5.34 Pattern of ftz expression 1115.35 Ocular dominance stripes in the visual cortex of the macaque monkey. 1155.36 Schematic of the visual pathway in the human brain 1165.37 Intracortical interaction kernels 1195.38 Reaction-diffusion calculation resembling ocular dominance stripes. . 1225.39 Ocular dominance stripes in the visual cortex of a monkey 1235.40 Ocular dominance patches in the visual cortex of a cat 125ixA.41 Phase plane representation of the Tnring instability in a two cell reactor. 134xAcknowledgementIt is a pleasure to thank Prof. Lionel Harrison for supervising this thesis, and for beinga source of inspiration through his devotion to science and sense of humour. Thanksalso to co-supervisor Prof. Geoff Hoffmann, for introducing me to theoretical biologyand for his stimulating approach to research. I am grateful to my thesis committee:Profs. Birger Bergersen, Myer Bloom, and Gren Patey, for encouragement and helpfulcomments. The following people also directly influenced this work: Prof. ThurstonLacalli (Saskatchewan), Dr. Ken Miller (Caltech), Prof. Zoltán Rácz (Budapest), andProf. Nick Swindale.•Thanks to David Holloway and Bernard Lakowski for commenting on the manuscriptand Florence Tam for helping with some of the figures. Thanks also to Robin MacQueenfor comments on a rehearsal of my talk.I learned much from my friends and colleagues during my years in graduate school. MarkLegros, Jacqueline Purtzki, and Rob Scharein deserve special mention.I am grateful to the members of my family for their emotional support.This thesis is dedicated to the memory of Enrico Kindl, whose graduate career was cutshort in summer 1987. Enrico’s enthusiasm continues to inspire all who knew him.xiChapter 1IntroductionFrom the point of view of the physical scientist, biology consists of a vast set of interestingand unexplained phenomena. One fundamental problem which is largely unsolved isdevelopment: how do the complex structures of an organism arise from a fertilized egg.Eggs are unicellular, contain a single nucleus, and are not structurally complex comparedto the adult form. From this initial state emerges an ordered society of differentiatedcells with various roles in the functions of the body. It is a deep mystery how structuresas complex as the human nervous system can self-construct.Several kinds of cellular activities may be distinguished in the process of biologicaldevelopment including: the regulation of gene expression, the movement, growth, differentiation, and death of cells, and intercellular signalling. A proper understanding ofhow the many cellular and molecular events are integrated to create a developmentalprogram is a still distant goal. It seems clear that such an understanding will requirenot only a description of the nature of the components of a developmental process, butalso a model for how these are organized. The central questions of developmental biology ask how order at one level emerges from the collective properties of the componentsat a lower level. Such problems require the theoretical techniques familiar to studentsof complex phenomena in the physical sciences: mathematical analysis and computersimulation. Attempts to analyze biological development must at this stage be contentwith semi-realistic models, and modest success. Part of the difficulty is in consequence ofthe immense complexity of living organisms. For the present, modelling efforts are most1Chapter 1. Introduction 2valnable in testing general hypotheses abont how developmental processes could work.1.1 Biological Pattern FormationBiological pattern formation may be defined to inclnde developmental phenomena leading to spatially ordered strnctnres. A few commonplace examples of pattern formation inbiological systems are: animal coat markings, the pigmentation of shells, insect segmentation, vertebrae, knnckles, and the digits of the vertebrate limb. Diagrams of the biologicalpatterns most relevant to this thesis may be found in chapter 5. How spatial symmetryis broken to produce regular patterns in tissue that is initially apparently homogeneous,is an interesting scientific question, independently of its importance in development. Thesingle greatest advance in the theory of biological pattern formation was published ina paper by the English mathematician Alan Turing 1 entitled ‘The Chemical Basis ofMorphogenesis’[136]. Turing’s fundamental contributions to several areas, most notablyto mathematical logic, the foundations of computing science, and artificial intelligence,in addition to his work on morphogenesis, established him as one of the great scientificgeniuses of the 20th century. Turing is also known for cracking the German ‘Enigma’code during the second world war. The influence of Turing’s work on morphogenesis mayperhaps be gauged by the record of citation of the paper (see figure 1.1). This graphshows clearly that interest in the paper is continuing to increase.Turing had a long term interest in biological growth. He collected fir cones exhibitingFibonacci series, studied spiral phyllotaxis, and was familiar with D’Arcy Thompson’sbook On Growth and Form. The morphogenetic theory is an attempt to understandthe mechanism by which biological pattern arises. The development of the theory was amajor preoccupation for the last years of Turing’s life, and there remains a considerable1An excellent account of Turing’s life and scientific work is given in Hodges’ biography [53]. Most ofmy statements about Turing are based on Hodges’ book.Chapter 1. Introduction 36040Cl,0:C-)2001960 1990Figure 1.1: Citations of The Chemical Basis of Morphogenesis’ per year since 1961. Datafrom Science Citation Index.1970 1980YearChapter 1. Introduction 4body of unpublished work in the Turing archives at King’s College, Cambridge.As mentioned in the last section, a large range of different types of processes areactive in even a single developmental event. Turing’s models concentrated on minimalsystems involving only chemical reactions and diffusion, often called reaction-diffusionsystems. Turing’s morphogenetic theory may be summarized as consisting of two mainaxioms (or hypotheses):1. Biological structure is preceded by a chemical pre-patterri of morphogens which controlgrowth or differentiation.2. Such chemical patterns of morphogens can self-organize from an initially homogeneousstate through reaction-diffusion.Turing demonstrated the validity of the second hypothesis mathematically, through thelinear analysis of reaction-diffusion systems. Recent experimental results [18], discussedin section 2.1.2, have confirmed the second hypothesis. With respect to the first hypothesis, spatial patterns of biochemicals have been observed in several developing systems(some of which I will discuss below), but in all cases the origin of the pattern is unknown. For example, in the fruitfiy Drosophila melanogaster a morphological pattern(the segments of the body plan) has been shown to be preceded by patterns of geneproducts which appear very early in embryogenesis; however the mechanisms controllingthe spatially inhomogeneous activation and inactivation of genes remain controversial(see chapter 5).Maynard Smith [93] has provided a qualitative explanantion for the initially counterintuitive result that reaction-diffusion may generate spatial patterns of chemical concentration (illustrated in Figure 1.2). Consider a pair of morphogens which are part of achemical reaction scheme such that one slow diffusing morphogen (labeled X) catalyzesits own departure from an equilibrium value X0, and the other morphogen (labeled Y)Chapter 1. Introduction 5inhibits its own departure from equilibrium value }‘, and diffuses rapidly. Further suppose that X catalyzes the formation of Y, but Y inhibits X formation. These catalyticinteractions may be represented abstractly with the following ‘community matrix’ 2[71]:HNow consider a one dimensional reaction vessel which contains an initially homogeneousdistribution of the morphogens, at a steady state of the dynamics. Thermal fluctuations,proportional to the square root of the morphogen concentrations, cause minute spatialinhomogeneities in the concentration profiles of the morphogens. Such inhomogeneties areamplified by the reaction-diffusion system described above. If the X has a small positiveincrease above the steady state concentration the autocatalysis will amplify this. Thegrowing X concentration stimulates Y production. The Y morphogen diffuses rapidlyoutwards from the site of the X peak and causes X and Y to drop below the steady statevalues in this region. The diffusion and antagonistic interaction of X and Y cause thespatial disturbance to propagate, forming alternating peaks and troughs of X and Y. Itis easy to see that the related system of catalytic interactions having community matrix:will also form patterns, but peaks of X will coincide with troughs of Y and vice-versa.Models having the first sign pattern are often termed activator-inhibitor models. Thesecond type are usually called depletion models as the kinetics may be thought of arisingfrom a depleted substrate and a self-enhancing substance which feeds on that substrate.One can guess that the patterns formed by the above instability have an intrinsic wavelength which will be determined by the rates of diffusion of the morphogen and the rates2A term from the theory of ecological networks.Chapter 1. Introduction 6at which they are formed. The simplest expression having dimensions of length one canform from a characteristic diffusivity, D and first order reaction rate constant, k, is:A (1.1)In section 2.3, it will be seen that this does in fact provide a useful estimate of the lengthscale of pattern formation by reaction-diffusion.Turing’s ideas on morphogenesis were slow to attract the attention of the biologicalcommunity, as evidenced by the infrequent citation of his paper in the nineteen sixties.While most developmental biologists are now aware of its existence, there remains awide variety of misconceptions about reaction-diffusion. This is partially due to themathematical nature of the theory, which is alien to many biologists. There have alsobeen a number of criticisms of the theory which I will mention briefly here; some of thesewill be discnssed in greater depth in chapter 5.One class of objections to reaction-diffusion theory is concerned with statements aboutwhat it cannot do. C. H. Waddington’s early (1956) and influential assessment of thecapabilities of Turing’s model for organizing biological patterns falls into this category[140]:It seems .. . improbable that fundamental rhythmic patterns, such as those ofthe somites of the vertebrate body, would be dependent on such an inherentlychancy mechanism as that investigated by Turing. Probably the processeswhich he .. . disdllssed play a part only in the quasi-periodic dapplings andmottlings which often fill np relatively unimportant spaces.Some of the computations presented in chapter 4 should dispel any such doubts thatthe reaction-diffusion mechanisms are incapable of generating very regular structures ina robust fashion. A related criticism of reaction-diffusion is that it is very sensitive toChapter 1. Introduction 7boundary conditions, or the exact geometry of the embryo, which can in some circumstances be highly variable. This is the problem of regulation in developmental biology:how do embryos of different sizes and shapes produce the correctly scaled structuresand organs? A paper by Bard and Lauder [6] came to the conclusion, based on theresults of some numerical simulations, that Turing’s model was not capable of regulation. Careful numerical and analytical work by Lacalli and Harrison [73] showed thatthis was not the case. Indeed, Lacalli and Harrison demonstrated that reaction-diffusionwas capable of explaining scale-invariance over a factor of about five in the slug of theslime mold Dictyostelium discoideurn. The paper by Bard and Lauder, a later paperby Bunow et. al.[12], and the studies of animal coat patterns by Bard [7] as well asMurray [101. 104] have all contributed to a misconception that reaction-diffusion modelsare incapable of producing striped pattern on two dimensional domains which are notnarrow. One consequence of this was the rejection of an early proposal for a mechanismto control segmentation in Drosophila [69]. A major goal of some of the work presentedin this thesis was to establish that reaction-diffusion is capable of stripe-formation, andto examine the conditions under which it occurs.Perhaps the most serious criticism facing Turing reaction-diffusion models is that theexistence of the morphogens remains hypothetical [44] Before addressing this objection,it should first be pointed out that the term ‘morphogen’ (which Turing invented), hasnow passed into the common biological vocabulary to describe the many form producingsubstances that have been discovered in the last few decades. This sense of the word isnot what was intended by Turing, who used ‘morphogen’ to describe the X and Y substances, which act not only as signallers of differentiation or growth but also as the keyparticipants in the self-organization of biological pattern. It is the interactions betweenthe morphogens, and their physical properties, that cause form to arise. The current‘morphogens’ are more accurately described as Wolpert positional signalers. With thisChapter 1. Introduction 8semantic point clarified, the question is: ‘Where are the Turing morphogens’? What isneeded to prove the existence of a Turing mechanism operating in a biological systemis to elucidate the dynamics of a chemical network for the synthesis and degradation ofan appropriate set of substances which are found in spatial patterns at some point indevelopment. Better still, one should know the time evolution of these patterns, andmeasure the diffusivities and reaction rates for the biochemical network under consideration. Such an experimental program faces many difficulties. The patterns are transitoryand must be recorded at the correct stage of development. Also the patterns may besensitive to boundary conditions, so an in situ method of detecting the morphogens isrequired, or spurious results are likely to be obtained. Further, if they are cell signallers,the morphogens are likely to be active at low concentrations, thus difficult to measure,especially in situ. Recent experimental work in Harrison’s laboratory on the marine algaAcetabularia [42, 43, 46] and in Cox’s on the cellular slime molds[13, 14], come the closestto satisfactory attempts to test reaction-diffusion models.This thesis is exclusively concerned with reaction-diffusion models of pattern formation and with a class of phenomenological models for ocular dominance bands in whichthe underlying mechanism is not explicitly specified. Discussion of the latter modelsis included because their dynamics closely resembles the behaviour of reaction-diffusionsystems. Many other models for biological pattern formation exist, but are not treatedhere. Most of the alternate models attempt to incorporate cellular activities beyond thekinetics of biochemical reactions. For example a natural extension of reaction-diffusiontheory is to use the morphogen concentrations to drive the growth of a domain, to modelthe morphological changes as tissue develops. Such a mechanism has been used to modeldichotomous branching in the cell wall growth of certain algae [45]. Other models consider the movement of cells in response to chemical or adhesive gradients, cell tractionforces and the response of the extracellular matrix and other mechanical effects. It isChapter 1. Introduction 9clear that such effects are important in some developmental processes. With many ofthese ‘mechanochemical’ models [103, 105, 112, 113, 114, 115] it is possible to obtaiupatterning behaviour through the coupling of mechanics and chemical reactions; in somecases the mathematics of these mechanisms is nearly isomorphic to the mathematics ofthe Turing instability. It may therefore be possible to apply some of the results of thework reported here to these mechanochemical models.1.2 Summary of ThesisChapter 2 introduces the Turing instability, presents techniques for the linear analysisof reaction-diffusion systems and discusses several hypothetical reaction-diffusion mechanisms. Chapter 3 starts with a discussion of the problems specific to reaction-diffusionsystems in two spatial dimensions and introduces the problem of biological stripe formation. A technique for the nonlinear analysis of pattern forming systems is appliedto stripe formation in two dimensions. Chapter 4 contains results of computations withreaction-diffusion equations. Finite-difference methods for the solution of the partial differential equations, and their implementation on digital computers are first discussed.Chapter 5 discusses two examples of biological pattern formation to which the resultsapply. Segmentation gene patterns in Drosophila are candidates for a biological exampleof reaction-diffusion, and were the initial motivation for my work on stripe formation.Ocular dominance patterns in the visual cortex are another example of self-organizingbiological stripes. T extend the analysis presented in chapter 3 to existing models for thedevelopment of ocular dominance patterns. Chapter 6 contains some concluding remarkson the thesis. In appendix A some work on models of two cells coupled by Fickian transport is presented. The problem was used by Turing to introduce the concept of diffusiveinstability. In appendix B I give the exact solution for the inhomogeneous steady statesChapter 1. Introduction 10in the case where each cell is modelled nsing the brnsselator.Chapter 1. Introduction 11abx—— YFigure 1.2: Development of a Turing structure in an activator-inhibitor reaction-diffusionsystem. After a sketch by Maynard Smith [931.CdChapter 2Reaction-Diffusion TheoryIn this chapter I discuss in detail the second axiom of Turing theory, that there exist physical systems capable of producing pattern in a self-organizing fashion. The mechanismwhich Turing presented for the self-organization of patterns of chemical concentrationis generally known as the Turing instability. A system exhibiting a Turing instabilityis stable to spatially homogeneous perturbations, but linearly unstable to heterogeneousperturbations having a wavelength within a certain range. In a reaction-diffusion systemthis type of instability is also called a diffusive or diffusion-induced instability becauseof the critical role played by diffusive transport of the system components. The mathematical origin of the Turing instability in reaction-diffusion mechanisms is derived anddiscussed with respect to a recent experimental confirmation of Turing’s predictions. Thestandard mathematical technique in the study of Turing models is linear stability analysis. This is equivalent to Turing’s exact solution of the linear reaction-diffusion equations.Linear stability analysis and related topics are presented in this chapter. Some of therelevant model chemical mechanisms used in this thesis, especially the brusselator andhyperchirality model, are also described.2.1 Reaction-Diffusion SystemsThe state of a reaction-diffusion system with N chemical species may be given at anytime by the N concentrations of the those chemicals at all positions in the spatial domainoccupied by the system. For the present purposes of discussion denote the state as the12Chapter 2. Reaction-Diffusion Theory 13set of N functions of the spatial variable i as:{X(i,t);i = 1,N} (2.1)Abbreviated henceforth as {X}. X is used for the amount of the component per unitvolume, as well as to stand for the actual chemical substance. We treat systems wherethe only processes affecting these N variables are chemical reactions occurring within thesystem, and physical diffusion. We assume that other effects such as convection, or theactive transport of molecules in a biological system are absent. For the time being, werestrict our attention to the interior of the system, returning later to the discussion ofthe boundaries where there may be an exchange of matter with the surroundings ( as inthe case of Dirichlet boundary conditions).The dynamical equations governing the time evolution of the state 2.1 may be foundfrom the conservation of matter at each position of the domain:rate of rate of production rate of influxincrease of = due to chemical + due to diffusionX reaction from surrounding volumeThen chemical reaction net production functions determine the first part of the RHS anddepend only on the state of the system as well as some control parameters. For reactionalone then:8Xt)= f({}). (2.2)The rate of influx due to diffusion at a given point which we write as Ij can be foundusing the divergence theorem. If Q is an infinitesimal ball at f’ and 6Q its boundingsurface, and Jj a flux, then:i=—f.n=jV.2, (2.3)Chapter 2. Reaction-Diffusion Theory 14where ñ is the outwards normal to SQ. So the rate of change of concentration at a pointdne to diffnsion is:8Xt)=—v.2. (2.4)We assume linear Fickian diffusive transport [33, 108]:= E (2.5)In the case of a dilute mixture of the N components (which should be the case formorphogens) cross diffusion may be neglected:= DS3 (2.6)and= —V. (DVX). (2.7)We neglect any concentration dependence of the diffusivities (again this is likely to be avery good approximation for dilute solutions) to obtain the usual form of the diffusionequation:= DAX. (2.8)A represents the Laplacian operator. The time-evolution equations for reaction-diffusionsystems are then:= f({X4) + DAX i = 1,N. (2.9)These coupled parabolic partial differential equations are known as reaction-diffusionequations. The reaction-diffusion mechanisms discussed in this thesis can all be writtenin this general form.2.2 Two Component Reaction-Diffusion MechanismsBy far the most common case of equations 2.9 in the literature is N = 2. Systems withN = 4 will also be encountered in this work. With N = 2 the two substances will beChapter 2. Reaction-Diffusion Theory 15denoted X and Y. Then the standard Turing eqnations are:X = f(X,Y)+Dx/XX,Y = g(X,Y)+DyzIY. (2.10)N = 2 is the minimum number of components of the system required for the Turing instability to be possible. The one component system does not exhibit a diffusive instability.For consider the following general one-dimçnsional single component equation:d2XX=f(X)+D_a_ (2J1)with a homogeneous steady state at X0:fix0) = 0 . (2.12)Linearizing the equation about this point:cPUU=kiU+D—- where U=X—X0. (2.13)dr2To test the stability of the homogeneous state consider perturbations of the form:&‘tsin(kx) , (2.14)where k is the wavevector of the sinusoidal perturbation. Then the dispersion relationfor 2.13 is:=— lc2D (2.15)Then if the homogeneous system is stable (k1 < 0), all inhomogeneous perturbations willalso be stable.Chapter 2. Reaction-Diffusion Theory 16Consider now the two morphogen system of equation 2.10. One is usually interestedin a case where there is a fixed point X0, Yo which is locally stable in the absence ofdiffusion:f(X0,1”) = g(X0,Y0) = 0 (2.16)The stability of this solution to local, spatially inhomogeneous perturbations is studied.First the equations are linearized by changing variables to departures from the steadystate, U = X — X0 and V = Y — Y, and discarding all terms nonlinear in U and V:U =k1U+k2V+DUV =3U+k4V+DyZXV (2.17)The constants k1 — k4, or marginal reaction rates, may be determined from:(&f kk1=2)= () , = () . (2.18)The derivatives are evaluated at the steady state. The set of coupled partial differentialequations 2.17 is the system Turing was able to solve exactly. Let the patterning domainbe a one-dimensional ring of continuous tissue of unit length. Then the appropriatebasis functions for looking at spatially patterned perturbations are plane waves withwavevector k = 2’irn where n E (—, cc):U, V o (2.19)This leads to the eigenvaltie problem:U k1—k2D U0 =. (2.20)V k3 k4—k2D VChapter 2. Reaction-Diffusion Theory 17The eigenvalues are:± 2 Tr+VTr2_4dei2where:Tr = kj—k2Dx+4---kydci = (k1—k2Dx)(k4— k2Dy) — (2.21)We now show that while the homogeneous system may be stable:a(k2 = 0) < 0 , (2.22)it may be unstable to perturbations with non-zero k2, that is spatial patterns.A phase diagram for the behaviour of the linear system 2.20 in terms of the traceand determinant of the linear stability matrix is plotted in figure 2.1. The homogeneoussystem is required to be stable, so:Tr =dci = k14—k23>O. (2.23)For certain values of the marginal reaction rates and diffusivities, the largest eigenvalue,can become positive for intermediate values of k2. Since Tr < 0, the appropriatebifurcation occurs when the determinant of eigensystem 2.20 changes sign from positiveto negative. For the critical value k0 of the wavevector for which this first occurs, asa control parameter is changed, the linear stability of the homogeneous state changesfrom a stable node to a saddle point. This corresponds to crossing from the third to thefourth quadrant in figure 2.1. Patterns of wavelength A = 2ir/k theu grow exponentiallyin amplitude. Farther from the intial onset of the instability, where there is a singleunstable wavelength, a range of modes is unstable.Chapter 2. Reaction - Diffusion T1 eorv 18Figure 2.1: Phase diagram showing regions of stability of system 2.20 as a function ofthe two invariants (trace and determinant) of the eigensystem. The parabola opening tothe right is defined by Tr2 = 4det.The two further required conditions for the occurrence of a diffusion induced instability are then:Tr = k1 + k4 — k2(D1 + D2) < 0det = (k1 —k2D1)(k4— k2D)— k23 < 0 for some range of k2 . (2.24)The first condition is satisfied if the system is stable to homogeneous perturbations. TheTr f4,det4,+tsecond condition may be simplified as follows: det(k2) is an parabola opening upwards,Chapter 2. Reaction-Diffusion Theory 19and will have negative values if the discriminant:(k1D2 +k4D1)2 — 4D1D2(kk— k23) > 0 (2.25)is positive. The wavelength at which the instability first occurs may be found by settingthe minimum value of the parabola to zero:1(k +k4mzn2D D2The dispersion relation, c+(k2), for the reaction-diffusion system, gives the rate of growthof a perturbation as a function of the magnitude of the wavevector of the perturbation.Typical forms of the dispersion relation for a Turing instability are given in figure 2.5. Perturbations of differing spatial wavelengths grow exponentially at different rates. Spatialpatterns having very short or very long wavelengths experience negative growth whilefor an intermediate range of wavelengths the patterns are amplified. The wavelengthhaving the largest value of cr will after some time dominate the pattern of morphogenconcentrations. The wavelength of the fastest growing pattern fitting the boundary conditions of the domain is close to the maximum in the dispersion curve. The wavelengthat the maximum of u+(k2) is determined by intrinsic properties of the reaction-diffusionmechanism, i. e. by the marginal reaction rates and diffusivities, and is consequentlyknown as the chemical wavelength.The existence of Turing structures in chemical systems is somewhat counter-intuitive.Diffusion normally smooths inhomogeneities in concentration profiles, and leads asymptotically to homogeneity. This is the ‘smoothing property ‘of the diffusion operator. Forconsider the diffusion equation:X = DAX. (2.26)Sinusoidal perturbations in the form of equation 2.19 are all damped, with the higherChapter 2. Reaction-Diffusion Theory 20frequency modes decaying the most rapidly:= —k2D (2.27)With the Turing instability, however, diffusion is found to have a destabilizing influenceon a chemical network which is otherwise stable. Spatial heterogeneities in the system,resulting from small initial heterogeneities or from thermal fluctuations in concentrationsmay become amplified to macroscopic pattern.2.3 Experimental Evidence for Turing StructuresThe year 1990 was a significant one for reaction-diffusion theory in that the first evidenceconfirming the experimental reality of the Turing instability was reported [18, 131]. Twogroups (both French) working with very different systems observed the sustained chemicalpatterns typical of Turing structures. Aside from the intrinsic interest these results holdfor modellers, it is widely anticipated that these results will lead to heightened interestin Turing’s theory amongst developmental biologists [111, 146].Castets et. at. [18] used a chlorite-iodide-malonic-acid (CIMA) reaction, alreadyestablished as a chemical system exhibiting complex kinetic phenomena. In Castet’sreactor (figure 2.2) reagents are fed into opposite edges of a thin polyacrylamide gel.This provides the flow-through conditions required for non-equilibrium pattern formation.The chemicals of the CIMA reaction diffuse from the two reservoirs to the middle of thegel, where the reaction takes place. A few hours after the beginning of an experimentstationary waves of the [I]/[I2] ratio (visualized using a starch indicator embedded inthe gel) are observed. The pattern is approximately a hexagonal array with an interpeakspacing of about 0.2mm, which was found to be truly intrinsic and not influenced bygeometric factors. More recent studies witb the CIMA-gel strip setup by Ouyang andSwinney [117], have established that striped patterns are also formed under conditionsChapter 2. Reaction-Diffusion Theory(a)21Figure 2.2: Experimental arrangement of Castets’ et. at., for the observation of standingwaves in the CIMA reaction. A and B indicate distinct reagents flowing along either edgeof a polyacrylamide gel strip. Figure courtesy of V. Castets.of high iodide or low malonic acid concentration 2• So far all of the observations arein accord with what is expected in Turing theory: stationarity, spontaneous symmetrybreaking, and intrinsic (chemical) wavelength are all features which characterize Turingstructures. In a discussion of the possible underlying mechanism for the experimentalobservations, Lengyel and Epstein [82] have suggested that the ion ClOg plays the roleof a fast-diffusing inhibitor, while 1, which binds to starch, has a reduced mobility inthe gel, and is the slower diffusing activator. Realistic modelling of the structures is stillat an early stage , and should be an area of productive research for the next few years.2These results have been extended to a wider variety of conditions, Harry Swinney, personalcommunication.3john Pearson, Los Alamos National Lab, personal communication.1et out1e3r A -inlet outle ‘\.(c)(b)Chapter 2. Reaction-Diffusion Theory 22-‘ T 1 •? •,,• 1/,,’ White0. 0. / • / •I •tL .,•-IIt t -\Dark% %%%%•%%%%% Ses’: n/,, ••, • ‘• White• •. • 0..__,1,p, f _.AFigure 2.3: Sketch of the microtubular dissipative structures observed by Tabony andJob [1311Tabony and Job [131] have reported stationary patterns in a quasi-biological system.Microtubules are polymers of a protein, tubulin, which are active in developmental processes in living cells. Tabony and Job polymerized tubulin in thin layers of solutionwith other reagents viz GTP (guanosine tn-phosphate) and acetyl phosphate present toprovide a thermodynamic driving force for polymerization. Using light scattering andsmall-angle neutron scattering, they showed that arrays of microtubules were formed asstripes of oriented nematic liquid crystals in which alternate stripes have orthogonal orientations of microtubules 45° to the stripe direction in one, 135° in the next, as sketchedin figure 2.3. The striped pattern had a wavelength of about 2mm and evolved over aperiod of a few hours displaying the types of defects seen with zebra stripes as well asocular dominance bands (as described in chapters 4 and 5). Using 31P NMR to followthe reactions involving phosphate, Tabony and Job demonstrated that formation andChapter 2. Reaction-Diffusion Theory 23maintenance of the structures depends on energy dissipation. The mechanism underlying pattern formation is at present unknown, but Tabony has suggested3 it may bereaction- diffusion.2.4 Turing Instability in the Ring of Discrete CellsBefore continuing with the linear analysis of the reaction-diffusion system introduced insection 2.1, I will discuss patterning domains consisting of discrete elements. The linearanalysis of the continuous and discrete systems is nearly identical and can be dealt withsimultaneously. I give an original discussion of the solution to the ring of cells using thediscrete Fourier transfomation, which permits a close analogy between continuous anddiscrete systems.The simplest discrete reaction-diffusion system is two well-stirred reactors (or cells)coupled by a linear exchange of morphogens. Appendix A contains a detailed discussionof this problem. A generalization of the two cell reactor is a ring of identical cellscoupled by linear exchange. This is the discrete one-dimensional system with periodicboundary conditions. Discrete systems are of interest in modelling problems where thecellular nature of tissue is important. In these circumstances the transport of morphogensmay occur by exchange through gap junctions, rather than diffusion in the extracellularmedium. Turing used the ring of discrete cells to discuss Hydra (a freshwater polyp) headmorphogenesis. Tn the context of continuous domains, the finite-difference approach tosolving the reaction-diffusion equations involves a discretization of space, resulting in aset of equations identical to the corresponding discrete-cell problem.The system is a ring of N cells. The reaction-diffusion equations are the 2N linear3J. Tabony, personal communication.Chapter 2. Reaction-Diffusion Theory 24ordinary differential equations:== k3X + k4Y + Dy(Yj1 + Y — 2c). (2.28)The morphogen profiles are represented by the 2 N-dimensional column vectors, X andY. The Fourier transforms of these are given by [130]:=FX and =FY, (2.29)where:1 1 1 ... 11 th ...1 2 ... 2O—1) where ii, = . (2.30)1 thThe inverse transform is computed with:1 1 1 ... 11 iv ... wn—11 in2 in4 w2&_1) where to = . (2.31)1 w1 w2(fl_1)In matrix form the reaction-diffusion equations 2.28 are:dx AB X— =, (2.32)dt Y CD YChapter 2. Reaction-Diffusion Theory 25where:kj—2D Dx U DxDx ki-2Dx Dy 0Dx U Dx ki—2Dxk4—2Dy Dy U DyDy k4—2Dy Dx U (2.33)U Dy k4—2Dyand:B—k], C=k31. (2.34)Factoring out the inverse transform:FO IL. AB FO (2.35)UFdtif CD OFwhence:d ( A*B (I = , (2.36)CD* i7in which:A*=FAF,B*=1PBF. (2.37)Use of the fact that B and C commute with the transform matrices, F, F has been made.A and B have circulant form, and are therefore diagonalized by the Fourier transform130, 54]. The diagonal entries are the Fourier components of the sequence:(k1— 2Dx,Dx,U, . ,Dx) , (2.38)Chapter 2. Reaction-Diffusion Theory 26given by:—4j2() k — 45jn2() •.. — 452((n —1))) (2.39)The equations for j and i’j then decouple:ci = — 4sin2() (2.40)cit Ic3 k2—4sinH)These may be solved separately for each i 0, 1, . N — 1.To solve the initial value problem, given20=(t=0), =(t=0) (2.41)we find the initial mode amplitudes:- 1--= —FX0,= 1P. (2.42)Then each mode amplitude equation 2.40 is solved. The soliltion in terms of morphogenconcentrations is given by the inverse transformation:(t) =i(t) = Fi7QL). (2.43)2.5 Conditions for Pattern Formation and the Chemical WavelengthThe discrete ring of cells and the continuous ring of tissue both gave rise to an eigenvalueproblem of the form:d c(t) = Ic1 — k2D Ic2 c(t) (2.44)d(t) Ic3 Ic4— k2D d(t)Chapter 2. Reaction-Diffusion Theory 27The same equation holds for any two morphogen reaction-diffusion problem, regardlessof the dimensionality of the patterning domain. For the continuous domains:= (2Km) , (2.45)whereas for discrete tissue:k2 = 4sin2Qirm/N). (2.46)The solution of equation 2.44 for the spectrum of growth rates and eigenvectors is identicalin either case. In the discrete case we bear in mind that only a finite number of modesare available: m taking values from 0 to N — 1. In the continuous case there is an infinityof modes with m ranging from negative to positive infinity.The general solution to 2.44 is straightforward:c(t) + -= c1ue +c2ve , (2.47)dQt)where o are roots of the secular equation:— k2D k2= 0 . (2.48)Ac3 Ac4 — k2DFor a fixed value of Ac2, u and v are eigenvectors:ki—k2Dx cf Ac2 U1= 0, (2.49)Ac3 k4_k2Dy_cr+with an analogous equation for v. c1 and c2 are to be determined from the initialconditions.The solution of linear two component reaction-diffusion systems is thus reduced tofinding the eigenvectors and eigenvalues of the Turing operator T:Ac1— Ac2Dx Ac2T = (2.50)I 12nChapter 2. Reaction-Diffusion Theory 28The eigenvalues, written in terms of two invariant quantities of T:TrT =defT = (k1—k2Dx)(k4— k2Dy) — k23 (2.51)are:= TT+ Discr = (TrT)2—4(detT). (2.52)The discriminant may be written in the simplified form [49, 75]:Discr = b2 + 41c2k3 , (2.53)with:b = (k4— — k2(Dy— Dx)) . (2.54)The larger of the two eigenvalues a±(k2) determines the linear pattern selection properties. The concept of ‘chemical wavelength’, as well as Turing’s conditions are all basedon analysis of the form of the function c+(k2), which is known as the dispersion relationfor the reaction diffusion system.2a(k) = (k1 + k4) — k2(Dx + Dy) + b2 + 4k2k3. (2.55)There is a single (non-zero) value of k2 for which(2.56)To see this, consider:2- = —(Dx -b Dy) + a(Dx — Dy), (2.57)where:a =_________. (2.58)b2 + 4k2k3Chapter 2. Reaction-Diffusion Theory 29Then the condition for the derivative to be zero is:D1a1p259Dxa+1at < 1 implies the unphysical result Dy/Dx < 0. It follows that:k23 < 0 , (2.60)for there to be a solution to equation 2.56 for non-zero k2. This solution is unique:/ k23bmax = +2a1/ 2 . (2.61)(1-a)The sign of b is already determined by the choice of parameters. By convention thenegative solution is used here; this is equivalent to the arbitrary choice of label Y for thefaster diffusing substance (Dy > Dx). It may be shown that this solution is a maximumin the dispersion curve by calculating the second derivative of the dispersion curve at theextremum. If a+(k2) becomes complex for some value of k2 then there is a discontinuityin the derivative at that point. The cusp is a local minimum of the growth rate (realpart of uj. There are only three possible local extrema of the dispersion curve (thethird at k2 = 0 which can be a minimum or a maximum). For two of these the derivative2.57 is undefined. Under certain conditions the solution to 2.56 is the global maximumof the dispersion curve, and selection of a particular wavelength of pattern occurs. This‘chemical wavelength ‘ may be calculated from the solution to equation 2.57:k2 —— + (Dx + 23x (2 2mary— xThe chemical wavelength is defined as:2K=. (2.63)‘m axLinear theory predicts that the modes closest to the chemical wavelength have the greatestgrowth rate and will dominate the pattern after a short time. If the chemical wavelengthChapter 2. Reaction-Diffusion Theory 30is much smaller than the dimensions of the patterning domain, then a linear system willselect a pattern with wavelength equal to AD.In reaction-diffusion systems, discrete or continuous, different choices of the marginalreaction rates, k1 — k4, and the morphogen diffusivities, D, D, lead to qualitativelydifferent dynamics of the pattern. Most of the asymptotic behaviour of the linear equations can be understood in terms of the dispersion relation derived above. Some of thefeatures of the nonlinear equations, especially wavelength selection, may be understoodin terms of the linearized equations. It is therefore useful to know how the functiona+(1c2) changes as the parameters are changed. Turing [136] listed a set of six conditionson the parameters by which one may divide the parameter space into regions with qualitatively similar behaviour. Turing provided little justification for the conditions, someof which have not received much attention in the literature of reaction-diffusion. Lacalliand Harrison [76] were the first to discuss all six inequalities in detail. They also provideda convenient graphical interpretation of Turing’s conditions, which is useful in computerwork with two morphogen models. In this section I present a derivation of the conditionswhich overlaps [75, 76] to an extent. Three conditions derived earlier, in the section onthe Turing instability, will not be rederived.The first condition is that the value of u(k2) at the critical point k2 given by equation 2.62 must be positive for wavelength selection to occur. This condition, the mostimportant of the six, was already derived in the chapter on the Turing instability. Forcompleteness, I list the formula again:(k1Dy +k4Dx)2 > 4DxDy(kik — k23) . (2.64)Condition two is that there must be a solution to equation 2.62; that is the maximumin the dispersion relation must occur for non-zero k2. From the expression for the chemicalChapter 2. Reaction-Diffusion Theory 31wavelength, it is sufficient that:k1 > k4 + (Dx + DY)1JD3 (2.65)The third condition divides parameter space into a region in which there are nooscillatory modes, and one in which damped or undamped oscillations may occur. Anecessary and sufficient condition for the eigenvalues of T to be real is:2+ 4k2k3 > U , (2.66)this simplifies to:k4 < k1 — 2—kk3. (2.67)If oscillations are possible, then it is important to distinguish between parameterchpices for which some oscillatory modçs may grow from those for which it must decayfor all values of k2. In the former case it may be difficult to sustain a stable pattern. Thefourth condition, as given in section 2.2, is merely a test on the sign of the real part ofthe eigenvalues of T (which are a complex conjugate pair):k1-j-k4<U. (2.68)As a further refinement on the previous two conditions, growing oscillations mayincrease in magnitude either faster or slower than the fastest growing stationary pattern(at kiax). The fastest growing oscillatory mode in the Turing equations has k2 = U anda growth rate of k1 + k4. So for pattern to grow faster than oscillations it is sufficientthat:b2 + 4k2k3 > k2(Dx + Dy), (2.69)evaluated at the turning point in the dispersion curve kax. Then using the same notationin equation 2.58:b/a> kax(Dx + Dy), (2.70)Chapter 2. Reaction-Diffusion Theory 32which simplifies to:kax = (k4 - k1). (2.71)Inserting and simplifying:k1 > k4 + (Dx+Dy) —k2k3. (2.72)The final condition that the homogenous system be stable, was found previously to be:k14 < k23 . (2.73)2.5.1 kk SpaceThe linear reaction-diffusion system:X k1X + k2Y + DzX,Y =3X+k4Y+DY, (2.74)can be written in terms of dimensionless parameters (see Lin and Segel [85] for a discussionof this technique in the context of the diffusion equation):= ki/I—k2s,=n = Dy/Dx. (2.75)The modified equations read:== (+)X’ + kY’ + nAY’. (2.76)Where:= /—k2k3,Chapter 2. Reaction-Diffusion Theory 33Table 2.1: Turing’s conditions in terms of dimensionless parameters.Label Condition Commenta k < nk — 2J max a positiveb k>14—(n+l)/i../7 A0realc k < 14 — 2 no oscillationsd 14 + 14 < U all oscillations dampede 14 < 14 — 4(n + 1)/4/ oscillations grow slower than patternf kA4 > —1 homogeneous system stable= r—k2ks/Dx (2.77)x’=(2.78)System 2.76 is completely equivalent to 2.75. One can study equations 2.76 with no lossof generality. This is beneficial as there are only three parameters to be varied in 2.76,where there were six in 2.75.With the above transformation, Turing’s conditions take a simpler form (see Table2.1). The conditions can be interpreted as planes curves in 1414 space. Five of theseconditions are straight lines, while the sixth is a square hyperbola. For a fixed valueof n (the ratio between diffusivities) the parameter space is determined. An exampleis given in figure 2.4. The conditions in the parameter space are labelled with theletters a-f. The dispersion diagrams corresponding to some of these regions are shown infigure 2.5. The regions bordered by these curves correspond to a given qualitative shapefor the dispersion plot. Specification of the three parameters determines the shape ofdispersion relation a(k2), although it leaves the units of the ordinate (time) and abscissa(inverse distance squared) unspecified. Lacalli and Harrison were the first to use suchdiagrams (known as Lacalli or 14k! plots) for exploring Turing models [76]. The regionsChapter 2. ReactionDiffusion TheoryI.1’i34Figure 2.4: Lacalli’s kk parameter space. Labels on Turing’s conditions are defined intable 2.1. Figure courtesy of L. G. Harrison.El10—1 b—24%4%4.—3Chapter 2. ReacL ion -Diffusion Th eorv 9—001—0•0100040—0004EFigure 2.5: Representative dispersion plots for four regions of kk space. Letters correspond to the labelled regions of figure 2.4. Here, k9 stands for the linear growth rate o.Figure courtesy of L. C. Harrison.S001SA CDChapter 2. Reaction-Diffusion Theory 36of primary interest for the formation of stationary Tnring structures are regions C andD. All of the computations in the present work have parameters from these two regions.Aside from a being convenient device for checking which form of dispersion curve a givenchoice of parameters leads to, the Lacalli parameter space has several other uses. Inchapter 4 the results of computations in which there is a gradient of reactants or reactionrates across the patterning domain are discussed. Such gradients may be used to limitthe spatial extent of pattern, mimicking the phenomenon of ‘localization’ seen in manybiological instances. For such studies one needs to know the characteristics of the modelat each point along the gradient. The switching on and off of the instability may becrudely interpreted as entering and leaving the regions C and D of Lacalli space. Suchinformation is viewed most easily by plotting the gradient in 14k space.2.6 Reaction-Diffusion MechanismsThere is little doubt that the ingredients required for pattern formation by a reaction-diffusion process are present in biological tissue. Biochemical pathways are sufficientlycomplex and nonlinear to generate the kind of kinetics required in Turing models. Theexperimental work on the in vitro microtubular dissipative structures, which is discussedin section 2.3, is an example of this principle: all of the components of the reactionare of biological origin. Another widely studied example of a quasi-biochemical schemeexhibiting spatio-temporal pattern-formation is the Belousov-Zhabotinski reaction whichhad its origin as a model of the citric acid cycle [145, 144]. A biological example ofchemical pattern formation by reaction-diffusion is the cyclic AMP signalling system inthe cellular slime mold Dictyostelium discoideum [139, 145, 144].Nonetheless, there is not yet a biological example of a Turing structure where thebiochemistry has been worked out sufficiently well for there to be a satisfactory agreementChapter 2. Reaction-Diffusion Theory 37between theory and experiment. In the past, Turing theory has been confined to the stndyof hypothetical chemical schemes. A recent exception is Lengyel’s realistic model of thegel strip reactor patterns [82] which is itself a highly simplified model of the CIMA reactor.Models of biological pattern formation have almost always made use of hypotheticalreaction mechanisms. Two such hypothetical biochemical schemes will be discussed inthis section. Many additional mechanisms have been reported in the literature, see forexample references [125, 30, 69, 95, 101]. The two presented have special relevance to thecurrent work. The brusselator, originally proposed by the Prigogine group in Brussels[108, 120, 121], is the most widely used reaction-diffusion model. Partly for this reason,it was used at an early stage of the project as a generic model. Also the brusselatorhas some non-generic properties that are useful in models of segmentation gene patternsin Drosophila (chapter 5). The second model discussed, the hyperchirality model, ispresented as a biochemically possible mechanism with a certain kind of symmetry presentin the dynamics which leads to special pattern formation properties. This model wasoriginally derived by Harrison and Lacalli [40] as a model for dichotomous branchingin algal cell wall morphogenesis, but may be more appropriate for instances of stripedpattern.2.6.1 BrusselatorPrigogine’s group considered the following set of chemical reactions [108, 120, 121]:aA-’XB+X -÷ Y+D2X+Y 2 3XX -* E (2.79)Chapter 2. Reaction-Diffusion Theory 38This is also known as the ‘trimolecular model’, after the stoichiometry of the third step.Lower case letters stand for forward reaction rates. Reactants A, B and products D,E areto be held at constant concentrations by providing a continuous source of the former, andremoving the latter at a sufficiently rapid rate that they do not build up to significantconcentrations. The reverse reactions can then be considered to occur at negligible ratecompared to the forward reactions. Such conditions are necessary to ensure that thesystem is held far from equilibrium which is a pre-requisite for pattern formation [108].The continuous flow of matter and energy through such a system is the origin of the termdissipative structure [108, 122, 123] to describe the material inhomogeneities developed.These conditions are duplicated in the various flow reactors that have been constructedfor experimental studies [18, 117, 131]. Such flow reactor conditions are easily envisagedfor biological systems: waste products can he removed by degradation, reactants suppliedby active transport or activation of precursors. The requirements that a system be opento the flow of matter and free energy are normal for living organisms.Rate equations can be derived from the chemical equation 2.79 above, using the massaction law:= aA—bBX+cX2Y--Xd,bY/Ut = bBX— cX2Y. (2.80)Justification of the use of mass action kinetics may be found in [108]. More complexkinetics may apply in biological reactions, but we ignore such considerations in view ofthe lack of detailed knowledge about the actual molecular events underlying biologicalpattern.The first step in the study of a mathematical model is to non-dimensionalize theequations so that the number of parameters to be varied is reduced. The rate constantsa and b may be absorbed into reactant concentrations A and B respectively. c and d mayChapter 2. Reaction-Diffusion Theory 39be set to unity with the following non-dimensionalizing transformation [27, 108]:1’ = dl=Y’ = (2.81)A’= vAB’=— DxLIX— dD = (2.82)Dropping the primes, the non-dimensionalized brusselator rate equations read:DX/8l = A-BX+X2Y- ,OY/Dt = BX — X2Y. (2.83)The brusselator has the homogeneous steady state solution:X0=A, Y0=. (2.84)The reaction kinetics may be written in terms of departures form this steady state x =X-A,y=Y-B/A:th B—i A2 x B 1= + (2Axy -b —x2 + x2y) . (2.85)—B —A2 y —1The nonlinear contributions to th and are equal and opposite. This is a special featureof the brusselator arising from the fact that autocatalytic creation of X takes place withthe destruction of Y. This property has the consequence that for a system with no-fluxboundary conditions, integration of x over the extent of the pattern shows that the totalareas of peaks (x > 0) is equal to the total area of troughs (x < 0) [48].Chapter 2. Reaction-Diffusion Theory 40Each pair of values A, B for the brusselator corresponds to a unique point in Lacalliparameter space:— B—i‘ AlE’— —A286Some points in the parameter space are not accessible to the brusselator equations.(2.87)so the area to the right of the hyperbola in figure 2.4 is excluded. Given values k, k, aunique pair A, B may be calculated, providing they exist:A_____—B= i’k’ (2.88)It is possible to draw a parameter space for the brusselator which is equivalent to Lacalli’sparameter space diagram”There are several variants of the brusselator model which will not be discussed indetail here. Notable is Tyson’s version [137] in which the roles of X and Y are reversedin the last two steps:A-÷XB+X -* Y+DX+2Y - 3YY —i F (2.89)This is mentioned as it is used by Thurston Lacalli in work on Drosophila segmentation[78]. This model is similar to the brusselator, with the important difference that the4M. J. Lyons, unpublished results.Chapter 2. Reaction-Diffusion Theory 41region of k k parameter space accessible to the model is more limited than for thebrusselator.Another model derived from the brusselator that will be mentioned later is Lacalli’sdual brusselator [79]:Precursors —* X1Precursors —* X2X1 -* X1+Y2X2—, X2+Y12X1+Y1 —* 3X1-j-Y12X2 + Y2 —, 3X2 -f- Y2X1 — ProductsX2 —* ProductsYl —* ProductsY2 —+ Products , (2.90)used as a model of the Drosophila segmentation gene network.2.6.2 Hyperchirality ModelHarrison’s hyperchirality mechanism [40, 41] was proposed in 1978 as a biochemicallypossible model for the kinetics of algal cell wall morphogenesis. The model possessesa dynamical symmetry that gives it a strong intrinsic capacity for generating stripedpatterns. The striping capacity of the hyperchirality equations was first demonstratedin attempts to model segmentation gene patterns in Drosophila [88], about a decadeafter it was devised. The hyperchirality model and Lacalli’s dual brusselator model 2.90remain the only stripe generating reaction-diffusion models based on explicit (thoughChapter 2. Reaction-Diffusion Theory 42hypothetical) biochemical schemes.The model is based on a mathematical model for spontaneous optical resolution, andinvolves two pairs of chiral enantiomers denoted XL, X, YL,, YD. These were called hyperchiral because the chirality was postulated as arriving from the two modes of attachmentto a membrane of proteins which, as usual, are not present as both enantiomers. TheX and Y substances are formed from two precursors A and B in reactions which arecatalyzed stereospecifically with bimolecular kinetics, the production terms of the rateequations being given by:7 42172 7 42172= ItXXn&D+ItXYn1L,XL = kxxAX+kxyAY,17 j n2 tr2ID = ‘X-°-“-D= kB2X. (2.91)is a rate constant for autocatalysis of the X morphogens; there is no such auto-catalysis for the Y system (kyy = 0). ky multiplies cross-catalysis of X by Y, witha cross-over in the chiralities of X and Y morphogens for X formation but not for Yformation. It is assumed that there is a fixed total amount of both X and Y substances:XL+XD =YL+YD = 14. (2.92)One can imagine such a situation arising on a saturated membrane or catalyst bed wherethere are a fixed number of binding sites for molecules of the X and Y specificities. Thisconstraint reduces the number of independent variables by two, and introduces into thekinetic equations 2.91 ‘displacement’ terms [40] of the form:— A2 [kXK(X + X) + kxy(Y + n)]. (2.93)Chapter 2. Reaction-Diffusion Theory 43The two remaining independent variables which are most convenient for the presentpnrposes are the optical asymmetries:X=XD—XL. Y}D—YL, (2.94)which represent departnres from the racemic state. In terms of these asymmetries thehyperchirality reaction-diffusion equations take the form:X = (kia[1 — (X2/P)] + k1)X + k2Y — k5XY2+ DxAX,Y =3X+k4Y—6+DY, (2.95)with the following parameter definitions:kia = kxxA2Fx, k1 kxyA2P/2Fx.—kA2F, k5 = kxyJI2/2Fx= kyxB2P , k4 —kyxBFj/2Py, Ic6 = kyxB2/2Py. (2.96)Ic1 = kia + kib, Ic2, Ic3, and Ic4 are already in the form of the linear Turing constantsas the homogeneous solution of equations 2.95 is X = 0, Y = 0. The special propertyof the hyperchirality equations mentioned in the first paragraph of this section, is thatit contains only cubic powers of the departures from the equilibrium state. Quadraticterms are found in the brusselator as well as in most other biochemically-derived reaction-diffusion mechanisms. Their absence in this model is a consequence of the symmetry ofthe the mechanism with respect to the interchange of the L and D morphogens.Any doubts about whether such seemingly contrived dynamics might exist in biological systems may be abated by considering some of the highly unlikely structures ofbiomolecules. One example I am familar with is the linear pentadecapeptide gramicidin5llaving L5 residues.Chapter 2. Reaction-Diffusion Theory 44A! found in the bacterium Bacillus brevis [87]. This peptide has a primary sequence consisiting of perfectly alternating L and D amino acids. The alternating sequence allows thepeptide to adopt a non-conventional secondary structure endowing it with ion channelproperties. Such highly ordered structures are common in biological molecules. It is notinconceivable that biological dynamics may also make use of high degrees of symmetry.Chapter 3Stripe Formation in Two Dimensions by Reaction-DiffusionBiological pattern is commonly formed in two-dimensional domains: single layers ofcells. sufaces of large single cells, single layers of nuclei in a multinucleate cell. In someinstances, such as mammalian coat markings [101, 105, 104], what is presumably a singlemechanism produces both striped and spotted patterns on various parts of the skin. Onemust then seek a mechanism which can be controlled by such influences as the shape ofthe region to behave in either of these ways.By contrast, when a pattern forms very early in development and hierarchically governs the essential body plan of the organism, one needs a mechanism with extreme reliability both as to the geometry of the pattern units (stripes or spots) and their precisetotal number. In this chapter I discuss models which are capable of producing stripedpatterns in a very robnst fashion, independent of the particnlar boundary conditions orsize of the system. The initial motivation for the investigation of the stripe forming models were the Drosophila pair-rule gene patterns treated in chapter 5. However, the resultsof this study apply to an entire class of models which possess a particular dynamicalsymmetry.3.1 Reaction-Diffusion Theory in Two DimensionsThe initial boundary-value problem of interest is the solution of:X = f(X,Y)+DxZXX,45Chapter 3. Stripe Fbrmation in Two Dimensions by Reaction-Diffusion 46Y = g(X,Y)+DyAY, (3.1)for the morphogen profiles X(i, 1) and Y(i, t) where iis the position in a two dimensionaldomain. I will consider rectangular domains of width Lr and length L8. The boundaryconditions may be cyclic, no-flux, or Dirichlet, or some combination of these.For a domain with cyclic boundary conditions the appropriate eigenfunctions of theLaplacian operator are plane waves:= kr + k3, (3.2)where:2irn 2Kmicr= —i-—, k = -i---, ri,m e (—oc,,o). (3.3)With no-flux conditions at all sides, the eigenfunctions are:cos(krr)cos(ks.s), (3.4)with:rnkrr 1c3=—, n,mE(O,). (3.5)The linear analysis of equation 3.1, was presented in section 2.2. For the two dimensionalproblems, the magnitude of the wavevector is found from:= k + k. (3.6)In models of one dimensional patterns the rate constants and diffusivities of thereaction-diffusion equations may be adjusted so that the chemical wavelength matchesthe experimentally observed periodicity. The linear behaviour of the model is usuallysufficient to determine pattern selection. A consequence of the larger number of modesavailable in a two dimensional domain is that linear pattern selection is no longer sufficientto describe the qualitative behaviour of models. Figure 3.6 shows the dispersion relationChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 47a(k2) for a set of parameters with the spectrum of modes available on a domain 66units long and 8 units wide with no-flux boundary conditions all around. In this casethe domain is narrow enough that. for large values of the chemical wavelength, it iseffectively one-dimensional. In Figure 3.7 the same plot of u(k2) is used, but with themode spectrum for a domain 66 by 16 units. In 3.6 the spacing of modes is sparseenough that, for small values of n, the linear behaviour may be adjusted so that onlya single mode is unstable. For the domain in 3.7 this is no longer practical: the modesare clustered even for small values of k2. Surface plots of Turing structures observed incomputations with the brusselator model on these domains are shown in figures 3.8, 3.9.In spite of the large number of closely spaced modes for non-narrow two-dimensionalregions, recent computational work [79 88, 106) indicates that there are reaction-diffusionmodels which have a strong intrinsic tendency to produce striped patterns, independentof the geometry of a region or its boundary conditions. This chapter employs an approximate analytical technique due to Haken [35, 36, 37, 38, 39] to establish that the stripeforming tendency of these models is a result of a dynamical symmetry of the reaction-diffusion equations. The models treated in this chapter have a spatially homogeneoussteady state labelled:Y(i)=Y0. (3.7)For this chapter a more compact notation for the state will be used. A transformationis made to variables which are departures from the spatially uniform fixed point:X(i,t)—X0=. (3.8)q2(r,t) Y(ñ,t)— YoThe reaction-diffusion equations may then be written:cctfa(q)+Dc4qci, cv=1,2, (3.9)Chapter 3. Stripe Formation in Two Dimensions0by Reaction-Diffusion 48Figure 3.6: Dispersion relation and spectrum of modes available on a domain 66 unitslong and 8 units wide with no-flux boundary conditions all around.a(1,1)Figure 3.7: Dispersion relation and spectrum of modes available on a domain 66 unitslong and 16 units wide with no-flux boundary conditions all around.(1,0) (1,1)(0.4)11(0,5)10(0,7) (0,9)Growth Rate of Linear Modes—58x66CeIIDoman (0,10)(1,3)(1,0)(1,6)1.2_______—-(1,7) (1,8) (1,9)—1—2—3—4(0,4)(0,5)2.5 10(0,7)Growth Rate of Linear Modes16 x 66 Cell Domain (0,1Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusionx7.3.X Morphogen Concentration49Figure 3.8: Turing structure computed using the brusselator model on the narrow domainfor which the dispersion relation is shown in figure 3.6.x9.X Morphogen ConcentrationFigure 3.9: Turing structure computed using the brusselator model on the domain forwhich the dispersion relation is shown in figure 3.7.16 x 66 Cell Domain4..—oce,,,,(‘#7%IChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 50where q = 0 is the spatially homogeneous steady state. The functions:fi(q) = f’(X,Y) f2(q) g’(),Y) , (3.10)may be calculated from the transformation 3.8. Efforts to model the self-organizationof Drosophila pair-rule gene expression products into striped patterns have led to theinvestigation of models anti-symmetric about q = 0:f(—q) = —f(q) (3.11)Then the operation:q —, —q (3.12)is a symmetry of the dynamics of equation 3.9. We are aware of three such reaction-diffusion models [106, 79, 40]. These three models have been found to have a strongtendency to produce stripes regardless of the boundary conditions [88, 79, 106] empiricallyby computation. The ability of a local property of the reaction kinetics to influence thepattern globally is at first sight surprising. By analysis in the weakly nonlinear regimeclose to the Turing instability, it will be shown that the only stable solutions of systems3.9 are stripes.3.2 A Class of Reaction-Diffusion Mechanisms Which Preferentially SelectStripesIn linear Turing theory (f(q) linear in q) the solution q(t) is a superposition of noninteracting modes whose amplitudes have exponential time dependence. The mode withthe largest positive growth rate dominates the pattern after a short time. With nonlinearf(q) the linear modes interact, and the pattern selection properties of a model are afunction of these interactions. In this section Haken’s slaving principle (which I will alsoChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 51refer to as the ‘adiabatic approximation‘ ) [35, 36, 37] is used to derive a set of equationsdescribing the time evolution of the amplitudes of the unstable modes when nonlinearitiesare present.The first step in the derivation of a set of Landau-Ginzburg equations correspondingto 3.9 is to separate the linear and nonlinear parts of the RHS:çj. = T,,pqp + Na(q) . (3.13)The Turing operator (my nomenclature) T contains the linear part of the time evolutionof the state:= Tq, (3.14)where:k2(3J5)The marginal reaction rates can be found from the net production functions as describedin the discussion of linear Turing theory. Close to the homogeneous steady state solutionthe noulinearities may be expanded in a Taylor series about q = 0:(1) (2)Na= gap7qq+ gasqaq7q +... . (3.16)fryThe coefficients gj’7,gj5 are constants to be determined from the net production functions of a particular model. For example the hyperchirality equations [40]:U =k1U+k2V—3-k6V+DzXU,V = U+k4V—7+Dv/XV, (3.17)lead to coefficients:= 0,.91111 = k5, 91122 = —k6, 92122 = —1c7, (3.18)Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 52all other gp are zero. For the brusselator:X = A-(B+1)X+X2Y+DxAX.Y = BX — X2Y + DRAY, (3.19)the coefficients are:gui = B/A 9112 = 2A 91112 = 19211 = —B/A 9212 = —2A 92112 = —1 . (3.20)For the above examples the polynomial expansion for the morphogen. net productionfunctions is exact and terminates at the cubic terms. For several reaction-diffusion modelsin the literature this is not the case. The Gierer-Meinhardt model [30, 95] for exampleyields an infinite series. This does not affect the conclusions of this chapter, as theanalysis is conducted close to q = 0 and only the leading order nonlinearities shouldmatter.Elsewhere it was shown that the linear Turing equations are exactly soluble. Thesolution may be compactly written in the form:q(t)= Eku3(k)e9, (3.21)which is a linear superposition of the eigenmodes of the Turing operator. Plane wavesare used to represent the spatial part of the modes as they are eigenfunctions of theLaplacian operator:=—k2e’. (3.22)Plane waves are utilized, as the case of periodic boundary conditions will be analyzed.The treatment of an alternate representation, for example the sines or cosines that arisewith Dirichlet or Neumann problems, is similar. The vectors u (Ic) are eigenvectors ofChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 53the Turing operator for a given value of k2:k1—DA Ic2 u(k) u(k)= a3(k) (3.23)Ic3 k4—D2z\ u(k) u(k)More compactly:= a(k)u(k)e. (3.24)j = +1— distinguishes between the two eigenvalnes a of T for a given . I have replacedIc with Ic = Ic when a quantity does not depend on the direction of the wavevector. Theeigenvectors will be normalized to unit magnitude:u(k) . u3(k) = 1. (3.25)T is not a self-adjoint operator (as Ic2 Ic3 in general) so its eigenvectors are not orthogonal. The left eigenvectors of T will be needed later:v(Ic)T = a(Ic)v(Ic), (3.26)with normalization chosen so:v’(Ic)u’(Ic) 6jt. (3.27)The solution to the linear equations gave mode amplitudes:(t) = eUJt. (3.28)The constants (O) can be found by Fourier analysis of the initial state. As the statevariable is real:= (t), (3.29)at t = 0 and at all later times.The eigenmodes of the linear Turing operator form a complete set of basis functionsfor the solution space of q(?, t). Since the time-evolution of patterns of a given spatialChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 54frequency (especially the steady state) is of interest, it is natural to attempt to write thesolution to the nonlinear equations in the same form as 3.28:q(t) = (3.30)j,kThe Turing operator is linear so:Tq > (3.31)j,kExpansion of N(q) into mode amplitudes causes a proliferation of indices reflecting thecomplexity of couplings between modes induced by the nonlinear terms. Separating thecubic and quadratic contributions:N(q) = N1(q) + N2(q) . (3.32)One finds:N’= Zo7 (z:(t)u(k’)e) (z t)u4”(k”)e) . (3.33)j’k’Similarly:= gaps (E t)uk1)eik’P) (z t)u(k”)e’) (z c(t)uc”(k”’)e”)76 j’k’(3.34)Taking the inner product on the left with v3(k)ek”9and integrating over the domain,using:J e1dA = (3.35)as well as the orthogonality between left and right eigenvectors 3.27 one finds:— if iNt Ji’i” t’t”R P’”” Ci’Ci”ti”. — -—‘ iS m,akk!kfls,s,,ukk1+k!f m.Ls kk’k”k’S,S,,’,,,, k,k’+k”+k”jnnj!fl ttjf!IChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 55where:ii’J”—.3 (IA i’(j,i\ j”(j,ll\ (1)akklk — ti1U 1u,. jg5afrv= v (k)u% (k’)u7 (k”)i4” (k”)g3j5. (3.37)876The reaction-diffusion equations 3.1 are thus transformed into an infinite set of ordinarydifferential equations for mode amplitudes. For the systems which terminate at thirdorder these ordinary differential equations are equivalent to the original partial differential equations. With higher nonlinearities present 3.36 represents the first stage of anapproximation which is only valid close to q = 0. In either case, this representationof the reaction-diffusion equations clearly displays the immense complexity of nonlinearsystems. The Kronecker deltas in the quadratic and cubic terms imply a special angularrelationship that mode wavevectors must fulfill in order to couple, greatly reducing thenumber of terms the series. Nevertheless, in this form, the nonlinear problem appears tobe intractable. It is difficult to know whether it will one day become possible to treatsuch problems analytically.An approximate analytical technique that allows one to reduce the complexity ofsuch coupled nonlinear equations is applicable very close to the onset of the Turinginstability. In this so-called ‘weakly nonlinear’ regime the modes separate readily into ahighly reduced set of unstable modes with positive growth rates, and stable modes withnegative eigenvalues. This situation is relevant to the Turing problem where a bandpasslimited range of spatial modes is unstable. We analyze the case where a single wavevectorlength k has a positive eigenvalue. We require further that the eigenvalues of the stablemodes be larger in absolute magnitude that the unstable modes. The exact conditionsfor the mode elimination procedure are:> OforChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 56a(k) < 0 otherwise,>> a(k (3.38)These conditions hold near the onset of the Turing instability, just inside regions C+D inLacalli’s Turing space. The mode elimination (or adiabatic approximation) takes place asfollows. From the outset it is expected that, as predicted by the linear analysis, the largestamplitudes will be those associated with planewaves with = k. The amplitudes of thestable modes should adjust quickly to the current amplitudes of the unstable modes, astheir linear relaxation times a’(k) are much shorter than the timescale for the unstablemode dynamics. We may solve for the stable mode amplitudes in terms of the unstablemode amplitudes by using the quasi-static approximation. That is for the stable modeODEs we use:0 for stable modes. (3.39)Then we may solve 3.36 for the stable mode amplitudes:—1 is sniff sunsc ci(k) E + E b%%ns , (3.40)JVfIVnkfffiii” jfjftjiffIt may be impossible to solve these coupled nonlinear algebraic equations exactly. However only the leading order terms are needed, since we are already making an approximation. Here, the analysis is considerably simplified if attention is restricted to the problemof immediate interest, i.e. systems for which the quadratic terms vanish:= 0 . (3.41)The leading order term in 3.40 is then third order in the unstable mode amplitudes. Sucha term does not contribute to the approximation for the unstable mode amplitudes up tofifth order, which is not needed here. The ordinary differential equations for the unstableChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 57mode amplitudes are, to third order:—+(j ‘ -- 3 42Sk c,iSJç + k,Ic’+k”+k”unstablemodesPresence of the Krouecker delta implies a geometrical relationship between modes contribnting to the snm (see figure 3.10). As there is a single value of = k0 the notationmay now be simplified:0- ==tb — 4is the angle that a given nnstable mode vector makes with a fixed axis. Note that asin equation 3.29(3.44)for the solution qQt) to be real-valued. Also it is clear that it is necessary to have:b<0, (3.45)for the global stability of the dyuamics.Equality 3.44 allows restriction of the variable to angles less than :0<’<K, (3.46)1”The coefficients 3 (self-coupling) and 6 (coupling with distinct modes) arise from combinatorial considerations when contribntions to the sum over modes are counted up.Equations 3.46 are the Landau-Ginzburg equations describing approximately the pattern formation close to the onset of the instability. The number of equations to be solveddepends on the boundary conditions of the problem. There is an equation for each ofChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 58V(a)—(b)Figure 3.10: Arrangement in k-space of wavevectors for modes having k k andsatisfying (a) = 1 in the quadratic term of (7), and (b) = 1 in thecubic term.Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 59the N degenerate modes fitting the boundaries. This represents an enormous reductionin the complexity of the original equations 3.36.Of primary interest are the steady-state and stability properties of the steady statesof equations 3.46 . Those steady states correspond to stationary Turing structures givenby a superposition of the unstable plane waves 3.22. The lack of any angle dependencein the coupling between modes and the absence of quadratic terms further simplify thesituation. It was a surprising result that it is in fact possible to write down all of thesteady states of 3.46 and describe completely the local stability of those steady states. Torender the analysis of 3.46 more concrete I will treat the simplest possible two-dimensionalcase: a square domain in which there are four plane wave modes unstable.(a) Simple example - four modes unstable. For concreteness we start with the simpleexample of the unit square with periodic boundary conditions. Suzpose that there areonly four unstable modes corresponding to striped patterns in mutually perpendiculardirections. The unstable modes have = 0, r/2, ir, 3K/2 and k = 2ir10,l E Z It isnecessary to restrict l to integers which cannot be written as the sum of two squares,or there will be other possible spatial modes. For example, the case l = 5 admits 4additional unstable modes. Then equations 3.46 are:= ao + 3be04o+ 6b0= a + 3br + 6bt0. (3.47)Equation 3.47 has fixed points:(I) = —a/3b,sj = 0 or = —a/3b,o = 0 (3.48)(II) o = = -a/9b. (3.49)Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 60The linear stability matrix=8/D is:A= ( + 6b(4 + 6b (3.50)6ba4 a + 6b(0 + Ievaluated at the fixed points. For solutions ITrA = —2a<0,detA = a2>0, (3.51)so both eigenvalues of the linear stability matrix will be negative. For solutions II:TrA =detA = _2 > o, (3.52)which means there will be a positive eigenvalue and a negative one. So the only stablesolution of this system is a set of stripes running in either direction x or y. The othersolutions, which are equal admixtures of four modes, are not stable and are saddle pointsof the dynamics.(b,) General Case - N unstable modes. It turns out that it is still possible to analize 3.46for the case of general boundary conditions. This is fortunate as it was the purpose of theanalysis to investigate in general terms the stripe selecting property of 3.9, independentof the given set of boundary conditions. We still treat cyclic boundaries, but it is alsopossible to carry out the analogous calculations for Dirichlet or Neumann cases.Suppose now there is an arbitrary number, N, of directions for which plane waveswith wavevector = k satisfy the boundary conditions. Then there are N coupledordinary differential equations to be solved. The mode amplitudes at a fixed point arelabelled so that:0 i=1M,= 0 i = (M+1).N. (3.53)Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 61The label i segregates the modes with non-zero amplitude from those with zero amplitude.Then 3.46 yields M coupled linear equations for the M unknownsa + + 6b > = 0 (3.54)j=1•MjiIt is possible to show that this set of equations is non-singular, and so has the uniquesolution:= 3b(2M— 1) j = . .. M. (3.55)Linear stability analysis is now used to examine the stability of these N solutions. Elements of the linear stability matrix A are given by:,1—A has structure:BOA = , (3.57)ocwhere B is MxM:ST••TTS•••T(3.58)TT...Swith:S = a+6(M+1)ba*a= (2M+1)’T = 6ba*a = (2M+ 1)• (3.59)C is (N — M)x(N— M) and diagonal:C=(a-)I. (3.60)Chapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 62Submatrix C has only negative eigenvalues. The eigenvalues of B may be found by firstsubtracting a multiple of the identity matrix I:B’=B—xI, (3.61)where x = T— S. The eigenvalues of B’ are then related to those of B in a simple fasion:(3.62)and it is straightforward to show that the eigenvalues of:1111 ...1B’ = T (3.63)1 •.. 1are TM and 0, from which it follows that the eigenvalues of B are:—2(M— 1)a a(M—l)T+S= 2M+1 (2M+1)’S-T2I- (3.64)The second eigenvalue is positive for M> 1 indicating that such solutions are not stable.Only the solution where there is a striped pattern M = 1 is stable. The mutual inhibitionof amplitude for each pair of modes is the cause of the instability of the mixed state.This is evident in the stability matrix A as:> 0 i 0 j (3.65)and< 0 . (3.66)This situation is described as ‘competition’ between modes. In the absence of cooperativeinteractions the only stable solution is one in which a single mode is populated andChapter 3. Stripe Formation in Two Dimensions by Reaction-Diffusion 63suppresses all other modes. In the above analysis periodic boundaries were assumed.It is straightforward to extend the result on stripe selection to rectaugular domainswith alternative boundary conditions. The striping propensity of this class of modelsis intrinsic to the mechanism: it depends on the form of the reactious alone. It isinteresting to note that the mode amplitude equations share similarities with the modelsfor optical resolutiou upon which the hyperchirality model is based. The competitionbetween optical enautiomers gives rise to the competition between patterning modes. Inthe models for optical resolution, the entities are two equally-matched molecular species(enantiomers) and one must ultimately destroy the other because of a particular kindof dynamics iu their formation and destruction. In the models for stripe selection, theentities are equally-matched patterns, and again only one can ultimately survive becauseof features of the dynamics.Chapter 4Numerical Solution of Reaction-Diffusion EquationsNonlinear reaction-diffusion eqnations are all bnt intractable with a purely analyticalapproach. Numerical experimentation with a model is essential to understanding it’spatterning behaviour beyond what is expected from the linear stability analysis. Computer simulation formed a significant portion of Turing’s work on morphogenesis in theearly 1950’s. Turing was a key figure in the development of digital computing machinesand programming techniques. Reaction-diffusion equations were therefore amongst tbefirst differential equations to be studied computationally. In spite of this, there does notappear to be an extensive study of numerical techniques for the simulation of reaction-diffusion models available in the literature. The methods commonly used are extensionsof algorithms for the solution of the diffusion equation, a more widely known parabolicpartial differential equation. Descriptions of these techniques may be found in standardbooks on the subject of numerical methods [55, 72, 81, 119], upon which this chapter ispartly based. This approach was continued in the present work.4.1 Finite-Difference Methods for the Numerical Solution of Reaction-DiffusionEquationsFor concreteness consider first one-dimensional reaction-diffusion equations for two morphogens with dirichlet boundary conditions. The initial boundary-value problem to besolved takes the form:64Chapter 4. Numerical Solution of Reaction-Diffusion Equations 65a2xX =82YY = (4.1)se(O,1), te(O,cx) (4.2)X(s,t = 0), Y(s,t = 0) given (4.3)The length of the domain has been scaled to 1 nnit, with no loss of generality. In mostcases the long-term behaviour of the solution is of interest. For a Turing instability thesolntion approaches a steady-state, and it is the concentration profiles of the morphogensat this state that are of main interest. To solve 4.3 numerically the spatial and temporalpartial derivatives are approximated using difference expressions. This converts the partial differential equations into a set of coupled difference equations which may be solvediteratively given the initial X and Y concentration profiles. To approximate the timederivative, first write X(.s, t) as a Taylor expansion from the current timestep:X(s, t + M) = X(s, t) + k(s, t)t + k(5,t)(t)2 . (4.4)The approximation:, X(.s,t + t) — X(s,t)45has truncation error:T.E.The error in 4.5 is first order in At, the timestep. Higher order expressions for X arealso possible, and will be mentioned below.Two distinct approaches to differencing the spatial derivative lead to two classes ofnumerical methods for solving the reaction-diffusion equations. In both cases the startingChapter 4. Numerical Solution of Reaction-Diffusion Equations 66point is again a Taylor expansion for X(s, t), this time in the spatial variable.X(s + As) = X(s) + X’(s)As + X”(s) (As)2+ X”(s) (As)3+ X’m(s) (As)4+“ , (4.6)then:X(s + As) + X(s — As) — 2X(s)— X”sX”(s)A N2 ... 4 7(As)2 ( sj . (.The approximate expression:82X X(s+As)+X(s—As)—2X(s)48)(As)2has a truncation error of:NT.E. =—(As)2. (4.9)so 4.8 is accurate to second order in As.With an explicit finite-difference method the approximation for the spatial derivativeis applied at the current timestep. With an implicit method, the approximation is appliedat the next timestep. These methods are discussed separately, below.4.2 Forward Time Differencing - the Explicit MethodWith the explicit method the spatial derivative is calculated at the current timestepusing the known solution values. Using the difference expressions listed above yields analgorithm that is first order in At and second order in As.The difference equation to be solved takes the form:X(s, t + At) — X(s, t) X(s + As, t) + X(s — As, t) — 2X(s, t)At = f(X(s,t),Y(s,t))+Dx (As 2(4.10)with a similar equation for Y. To render the notation more compact, the indices shownin the figure 4.11 are used to label the finite difference grid and timesteps. Then theChapter 4. Numerical Solution of Reaction-Diffusion EquationsSpace67Figure 4.11: Discretization of time andtions by finite-difference methods. Thein the explicit Euler technique.difference equations are:space for the solution of reaction-diffusion equabold line shows the computational molecule used(4.11)EH-____ n±1nn—ij— 1 j+1X’ = X +f(X,}7)t± [x1 +X71 - 2X],y+1 = +g(X,Yt + [1 +Y - 2].4.11 is a prescription for finding the state of the system at the (n + l)th timestep giventhe state at the timestep. Given X, Y one calculates YJ’. Given the initialstate (n = 0) the solution may be marched forward ‘in time.Chapter 4. Numerical Solution of Reaction-Diffusion Equations 68The three varieties of boundary conditions usually encountered are dealt with in astraightforward fashion:• Cyclic: define the indexj modulo N. Then j+i = Oifj = N—i and j—i = N—iif j = 0• Dirichlet: X, P7 do not change with time for j = 0, N — 1• No-Flux: define buffer nodes at j = —i,N withX1,Y = X07’,Yj’ and X,YJ =X_1,Yfl_1. This ensures that no matter flows beyond the boundaries.The most serious drawback of the explicit method is that it can exhibit numericalinstability for certain choices of L\t and z5. Numerical instability is a condition in whichthe approximation causes numerical errors to be amplified at each timestep. Such asituation is generally unacceptable for a computation. It is well known that the explicitmethod for solving the diffusion (or heat) equation:= , (4.12)becomes numerically unstable when:DAt 1(As)2 2 (4.13)A similar condition arises when the explicit method is applied to reaction-diffusion equations. The stability analysis presented here makes use of techniques originally devised byvon Neumann and is analogous to the derivation of 4.13. This analysis of the stability ofthe explicit method for reaction-diffusion systems has not been previously discussed.The eigenmodes of the difference equations 4.11 are of the form:=(4.14)Chapter 4. Numerical Solution of Reaction-Diffusion Equations 69k is a wavevector limited to a finite, discrete set of values determined by the boundaryconditions and the grid spacing As. The linearized reaction-diffusion difference equationsare:= k1U + k2Vf + ()2 [u;+1 + U1 - 2U$],At= k3U + k4 + (A)2 + V — 2] . (4.15)Substituting 4.14 into 4.15 yields a set of linear difference equations for , ,y:n+1= A , (4.16)n+1with:A=kiAt+1_4jjhisin2( ) k2At. (4.17)k3AtThis system will not suffer numerical instability provided:<1, (4.18)where the A are the eigenvalues of A. This condition leads to a complicated expressionfor the maximum allowable value of At/(As)2in terms of the parameters. This conditionmay most simply be expressed in a form originating with R. May described in [71]:2> 1 + Det(A) > Tr(A) . (4.19)This becomes too unwieldy to be of practical use here, although it could be convenientlyevaluated with a simple program. A more useful criterion may be derived by calculatingthe approximate value of the eigenvalue of A which has the largest magnitude. We assumethat the timestep is chosen sufficiently small that the reaction part of the equation isbeing accurately calculated i.e.:kAt s U where i = 1,4. (4.20)Chapter 4. Numerical Solution of Reaction-Diffusion Equations 70Then the eigenvalues of A will be close to the eigenvalues of the system:1 — 0(zXs) 2. (4.21)0 1 4D,t 2/ks()2 SZ7? \ 2The most stringent condition on the timestep follows from:Dv At1 4(AS)2 <1, (4.22)which implies:. (4.23)In practice, this formula is an adequate guide in the choice of timestep for the computation.Expression 4.11 is only first order accurate in time. Because the solution approachesa stable attractor this is sufficient accuracy to avoid incorrect conclusions on the qualitative behaviour of reaction-diffusion equations. In some circumstances a more accuraterepresentation of the time derivative may be required. One approach is first to applythe finite-difference method to discretize the spatial derivatives to generate a set of 2Ncoupled ordinary differential equations for the 2 morphogen concentrations at N gridpoints:= f(Xj,})+Dx (x3+i+x 2x)= g(Xj, ) + D (As)2- 21) (4.24)One may also think of this system as describing a discrete set of coupled cells. The manymethods for solving ordinary differential equations may then be used, such as the popular4th order Runge-Kutta which is considerably more accurate than the Euler or first ordertime-differencing method. With the Runge-Kutta method intermediate time-steps aretaken to obtain an estimate of the solution which is then used to calculate the solutionChapter 4. Numerical Solution of Reaction-Diffusion Equations 71at the next timestep. To use the Runge-Kutta method, closed form expressions for thederivatives of f(X, Y) and g(X, 1) with respect to X and Y are required. This is not aproblem for the models being studied. In general a timestep requires more computationwith this method, as several evaluations of these derivatives are needed. This methodcan be used when a more accurate representation of the steady state profile or time-dependence is required. It requires somewhat more computation (function evaluations)but may be worthwhile if the system is reaction-dominated (i.e. Dy/Xt/(As)2is small)and taking a larger timestep may improve efficiency without instability occuring in thespatial differencing.4.3 Implicit Finite-Difference MethodsThe numerical instability of the explicit scheme may be avoided by approximating thespatial derivative at timestep (n + 1) (where the solution is known up to timestep vi):—xrb r’ x1 2XTh+lAt = f(X7,}7) + Dx ( + ) . (4.25)This leads to a set of 2N coupled linear equations for the 2N unknowns (X’1,4n+1))which may in general be solved. To obtain the time-solution of the initial boundary valueproblem such a set of equations must be solved for each timestep of the computation.The equations have the form:1+2a —a 0 0 X1 X+f(X,YjjAt—a l+2a —a 0 X1 X+f(X,Y1jAt0 —a 1+2a 0 X1 X+f(X’,Y’)At0 0 0 1 + 2a - X1 + f(X_1,,_1)At(4.26)Chapter 4. Numerical Solution of Reaction-Diffusion Equations 72for the case of no-flux boundary conditions with:DAta = . (4.27)A similar equation for 1 must be solved at each timestep. With no-flux or Dirichletboundary conditions these equations are in tn-diagonal form, for which there is a veryefficient algorithm [119]. Stability analysis shows that the implicit method is unconditionally stable, so that there is no limit on the size of timestep that will give a stablesolution. However, the truncation error is of the same magnitude as the explicit method.In fact, the presence of reaction terms in this system limits the use of superior stabilityproperties of implicit methods. Estimates of the error indicate that in most circumstancesof concern here, the implicit methods are only of limited use. Specifically they are mostuseful close to the primary bifurcation where the computation must run for a long timefor the pattern to grow.4.4 Two Spatial DimensionsSolution of the reaction-diffusion equations on two-dimensional domains is possible througbgeneralizations of the one-dimensional techniques. The main additional requirement is arepresentation of the spatial derivatives as a difference expression on the finite differencegrid. The method presented here may be used on arbitrary finite difference grids. Whileonly rectangular grids with lattice spacing equal in the two dimensions were used in thiswork, I also give the discretization on a hexagonal network to show how to proceed moregenerally.First the Laplacian derivative is approximated using its integral representation:AXfVX.d. (4.28)This is a line integral around the boundary of a region surrounding a node point (seeChapter 4. Numerical Solution of Reaction-Diffusion Equations 73Figure 4.12: Discretization of the Laplacian operator in two spatial dimensions. A, squarelattice. B, hexagonal lattice.figure 4.12). i is a unit vector perpendicular to the boundary, and pointing outwards.In the case of a square lattice the resulting discretization is:1 7 V- ‘.. 17 17 171fj.A() 2-”O -“3’O 4-O— a± am a+ a= a2(X1+a2+ X3 + X — 4X0)a a(4.29)For an hexagonal network it is:2 a(Xi_XOX_ OX4XOXS_XO6)a, a a a a a(4.30)Or,+X23456-60). (4.31)In the case of a forward time difference generalization of the one-dimensional methodis straightforward. One can again simply calculate the state at timestep n ± 1 giventhe state at using the equations 4.11 but replacing the spatial part with expression1.. .A 3.6.a.2..5.B4 .3Chapter 4. Numerical Solution of Reaction-Diffusion Equations 744.29. The ease of implementation of the two dimensional method and its satisfactoryperformance for simulations of the Tnring instability have led to its widespread use. Inthe present work, this was the most heavily used technique.The equations for the backward time difference are no longer as straightforward aswas found for the one dimensional case. If one writes the state of the system as a vectorof dimension N x P (where N and P are the number of grid points in the two directions)then the matrix M of equation 4.26 is no longer tn-diagonal. It is, however, sparse andone may resort to a sparse solver to implement the implicit scheme. More commonly,the alternating-direction implicit algorithm is used. With the ADI algorithm, alternatetimesteps are used to calculate the matter fluxes in the two directions independently ofeach other. So at each timestep, one has N (or P) one-dimensional problems to solve,involving the solution of P (or N) tn-diagonal systems of equations.In most cases the reactions terms on the left side(4.32)require that we take a sufficiently small timestep to ensure accuracy and reduce the (1storder) truncation error due to this part of the discretization. This limits the benefitsgained over the explicit method, so that the increased stability of the implicit methodmay not be a significant advantage. While it is clear that in the case of heat conductionor diffusion problems the implicit (or semi-implicit as in the case of the Crank-Nicholsonmethod) are superior to the explicit method, this is not so clear for reaction-diffusionproblems.4.5 Simulation of Two-Dimensional Pattern FormationThe software used in this project was designed to be as simple and portable as possible.Fortran was the programming language used in nearly all of my numerical work, as itChapter 4. Numerical Solution of Reaction-Diffusion Equations 75is still the most widely supported language of the University Computing Services. Thecode relies on no external system routines, and can be compiled with little modificationon several different machines, from IBM AT compatibles through Unix workstations ofthe new ‘Unixg’ network at UBC, to the IBM 3090 150s running AIX. The latter machineproved to be the main workhorse of my two-dimensional calculations. At various stagesof my research the IBM 3090 was virtually unused, allowing me to carry out a largenumber of numerical experiments with a short turnaround time between jobs.Simplified one-dimensional versions of programs used to solve reaction- diffusion equations were developed for interactive use on IBM compatibles. The Waterloo Watcom Fortran77 compiler with GKS graphics library was used to compile and run the programs.The one-dimensional programs were used routinely to test new parameters or gradientsbefore running calculations in two dimensions. For a reasonable number of finite difference nodes (20-40) one-dimensional calculations are sufficiently rapid that the resultscan be viewed as an animated sequence as the calculation proceeds. I have found thatthese one-dimensional ‘cartoons’ of reaction-diffusion dynamics are of pedagogical valuein introducing the concept of self-organizing pattern formation to newcomers in the field.One-dimensional morphogen profiles can most easily be represented by line plots ofconcentration versus position in the patterning domain. With two dimensional dataone has a number of options for ways to represent the outcome of a calculation. Fourpossibilities that I considered and experimented with are: surface plots, contour plots,false colour mappings, and grayscaling. The most convenient options for most purposesare surface plots or grayscaled images. Contour plots are not as quickly interpretedas the other means of data representation, although they are an excellent way to codequantitative information on pattern amplitude. I used an hue-saturation-value colourmodel to code morphogen concentrations as hues in graphics programs to routinely viewthe results of my computations on an AT compatible. Such a model was also used in aChapter 4. Numerical Solution of Reaction-Diffusion Equations 76video of reaction- diffusion dynamics. Colour graphics, though often aesthetically moreattractive than grayscaled images, are not so readily or cheaply produced, and so werenot used for hardcopy records of the computer results.I found that the features of two dimensional patterns of greatest interest to me in thiswork (spots vs. stripes and other qualitative aspects of biological patterns) can be mosteasily discerned in grayscaled images. Surface plots were used in my earliest work on thisproject and found to be inferior for complex patterns. Elsewhere in this thesis surfaceplots are used to represent some simple two dimensional patterns. Grayscaled outputhas the further advantage over surface plots that it actually resembles the patterns beingmodelled. These include photographs of zebra coats, immunohistochemical visualizationof protein distfibution patterns in Drosophila and the autoradiograms of striate cortexshowing ocular dominance bands.-All of the grayscaled output at the end of the section was produced by first mappingthe double precision numbers in the range of concentrations at a given timestate to theintegers between 0 and 255. The resulting two dimensional arrays of integers were thenprinted by constructing simple Postscript programs to drive a laser printer.4.6 Description of Numerical ResultsThis section contains selected reaction-diffusion computations on two dimensional domains. Morphogen concentration profiles are displayed as grayscaled images at severalintermediate timesteps in the computation. The figures are collected together at the endof this section, preceded by a discussion summarizing our conclusions on the behaviourof the systems studied. The discussion is based on the figures presented here as well asa large number of similar numerical experiments with reaction-diffusion systems.Linear stability analysis predicts that modes with a wavelength close to the maximumChapter 4. Numerical Solution of Reaction-Diffusion Equations 77of the dispersion curve should dominate the the system after some time. Roughly similarbehaviour is observed in calculations with the nonlinear equations when the parametersare chosen to lie in regions (c) and (d) of Lacalli’s Turing space. That is, the steady statesare characterized by a periodicity as close to the chemical wavelength of the linearizedsystem as the boundary conditions permit. In two dimensions several degenerate modesmay share the same value of k2. Here the behaviour of nonlinear systems departs fromwhat is expected from linearization of the equations. Figures 4.15 and 4.17 illustrate thatthe qualitative nature of the patterns formed varies for different models. The parametersof the two models were adjusted so that they share identical marginal reaction rates anddiffusivities, and only the nonlinear terms differ. The simulation with the brusselatorleads to a steady state which is an ordered array of localized concentration maxima orspots. The simulation with the hyperchirality equations leads to a pattern of irregularstripes or elongated blobs (which I include in my definition of stripes). Further iterationsof such computations are observed to lead to a very regular pattern having periodicityin only one direction. Another comparison between two models is shown in figure 4.19.In this case there is only a single wavelength for which the growth rate is positive. Thebrusselator computation shows that the degenerate plane waves with this wavelengthsuperpose to form an array of spots (an hexagonal array for this calculation). Theother part of this figure shows stages in a calculation with stripe-selecting dynamics.One spatial mode predominates and suppresses the other modes. Which mode ‘wins’is arbitrary and depends on the stochastic elements of the simulation, i.e. the randominitial state, and the fluctuations during the time evolution of the system.In the previous chapter it was shown that models with odd net production functionsselect striped patterns. The above computations, and many other computations undera large variety of conditions, have confirmed this result and suggest it extends into theentire pattern forming region of parameter space (region (c) and (d) of Turing space).Chapter 4. Numerical Solution of Reaction-Diffusion Equations 78The numerical work has also suggested that the preseuce of quadratic terms induces atendency to form spotted patterns. The factors influencing the stability of spots appearto be more complicated than merely to involve the form of the net production functions.I have found that in many cases the brusselator can make very regular striped patterns.The empirical results indicate that the important factors in addition to the strength ofthe quadratic nonlinearities include the shape of the domain, the boundary conditions, aswell as the linear parameters (which affect the spectrum of eigenvalues of the patterningmodes). Figure 4.20 shows the effect of the strength of the quadratic term in the followingrate equations:X = kiX+k2YkgXYkrX3,= k3X + k4Y — k6XY2 , (4.33)for increasing values of kg. The system has an increasing tendency to produce spottedpatterns as kg is gradually increased. The sign of kq used in these experiments is unimportant, as one may replace kg with kg and the equations 4.33 are unchanged exceptthat X —* —X and Y —* —Y . So with negative values for the computations in figure 4.20, one expects to see arrays of spots which are concentration minima instead ofmaxima. This was verified numerically.The striping propensity of the class of mechanisms with kg = 0 is inherent in thenature of the dynamics of a model and not dependent on the boundary conditions. Calculations with cyclic, no-flux, and dirichlet conditions on domains with various sizes andshapes invariably show this to hold. In addition the presence of gradients in systemparameters do not affect this property. The numerical work shows that the striping tendency of the class is very robust. Figure 4.2t a,b,c show the steady states of a calculationwith equations 4.33 (kg = 0) on a cylindrical domain with no-flux boundaries at top andChapter 4. Numerical Solution of Reaction-Diffusion Equations 79bottom of the cylinder. There are several degenerate striped final states observed depending on the seed given to the random number generator for the concentration fluctuations.Figure 4.21 d,e,f are similar but diffusivities have been reduced by about 50% so thechemical wavelength is shorter. In panel f stripes running the length of the cylinder areseen. At the outset of the project the extreme example of striped pattern formationexhibited in figure 4.21 f was considered impossible for reaction-diffusion mechanisms.On small domains, where diffusive coupling between disparate regions is strong, orclose to the onset of the diffusive instability where there is only a small range of unstablepatterning modes, the morphogen profiles are observed to approach a regular profile veryrapidly. With larger regions (compared to the chemical wavelength) it is found that thepattern initially formed consists of several domains of differing orientations of the stripingdirection. These orientation domains are connected by topological singularities in theorientation field. The nature of the singularities that can occur in this class of patternshas been discussed in the context of the systems of ridges found in fingerprint patterns[118]. The two basic types of singularities for ridge systems are the triradius and the loopas illustrated in figure 4.13. Both of these singularities may be observed in calculationswith the reaction-diffusion system in equation 4.33 (see figure 4.14). The loop is notso common on domains with periodic boundary conditions as there are few free ends ofstripes. Combinations of these two basic types may also be observed by studying thepatterns. The pattern defects have been observed to anneal if the simulation is allowedto run for a large number of iterations. A video sequence [91] of the long time behaviourof reaction-diffusion patterns reveals the vital importance to the annealing process ofrandom concentration fluctuations by which singularities of opposite topological chargemeet and annihilate leaving regularity. Computations in the absence of noise appear toget stuck in disordered states. The factors influencing the annealing have not been fullyexplored. More likely means by which developing organisms may avoid pattern defectsChapter 4. Numerical Solution of Reaction-Diffusion Equations 80Figure 4.13: The two basic types of singularities (or topological defects) for ridge systems,adapted from [118]. A, loop; B, triradius.are biased initial states or spatial inhomogeneity of the dynamics. These are discussedin the next section.The above computations showed it is possible to devise a reaction-diffusion mechanismwhich can produce striped patterns in a very robust fashion. With respect to the questionof how a biological organism might make use of such mechanisms to self-organize a periodic array of repeated parts, there are two immediately obvious unsatisfactory features.First, the striped patterns contain some disorder in the variation of the orientation ofstriping over the pattern-forming field. Domains of locally oriented stripes are connectedby defects in the pattern. The term defect is used in the sense of a topological defect,that is a singularity in the orientation field of the pattern. From the point of view of anorganism trying to develop ordered structures that will be the basis of the architecture ofthe whole body or of essential organs, an alternate definition of the term defect meaningerror or mistake (as it might be used to describe a defective manufactured item) is alsoappropriate. Studies of the long time behaviour of the striped patterns suggest that someChapter 4. Numerical Solution of Reaction-Diffusion Equations 81Figure 4.14: Examples of computed reaction-diffusion patterns displaying the two basictypes of singularities in figure 4.13. Calculated using equations 4.33. A, loop; B, triradius.ABChapter 4. Numerical Solution of Reaction-Diffusion Equations 82that some reaction-diffusion mechanisms are able to correct defects. However, a moresatisfactory way to organize ordered structures is to avoid defects from the outset of pattern formation. A few simulations are presented here showing that the use of monotonicgradients in the initial concentration profiles of the morphogens can accomplish this goalquite easily.A second unsatisfactory feature of the striped patterns discussed so far is that thereremains a degeneracy in the possible orientation of the stripes. Figure 4.21 shows severalcomputations with different initial states. Because of the small size of the patterning fieldthere is a single orientation domain, and absence of defects in the patterns. However,it is clear that several patterns are still possible as steady states of the system withstripes running in different directions. An organism such as Drosophila must be able tocontrol with a high degree of reliability which of a number of degenerate steady statesa particular computation selects, independent of the initial state or stochastic input tothe patterning mechanism. A clue to one solution of the problem comes from the resultthat the final pattern depends on the initial state of the system. For the linear system,two modes with the same growth rate will maintain the same ratio of their amplitudesfor all times in a calculation. For a particular mode to dominate it suffices to give thatmode a bias in the initial state, locating it in the basin of attraction of the desired steadystate. My numerical experiments along these lines have shown that this procedure workswell. The results of two such experiments are shown here. In figure 4.22 the linear Tnringparameters are chosen reasonably close to the primary bifurcation (not so close that thereis only a single positive eigenvalue). A gradient in the X and Y concentrations from topto bottom of a cylinder with no-flux boundaries was used as the initial state. The Fouriercomponents of the initial state consist of plane waves in the direction of the length of thecylinder. The resulting steady state is a set of stripes with k aligned with the gradient.Figure 4.23 shows the same computation but with no added bias in the initial state. TheChapter 4. Numerical Solution of Reaction-Diffusion Equations 83pattern is tending to a single oriented set of stripes, but more slowly (as the defects mustannihilate). The orientation of the final state is not the same.Further away from the bifurcation a gradient in the initial concentrations has similareffects. Figure 4.24 shows a set of wavy stripes which are absent of defects and havethe desired orientation. The origin of the periodic wavering of the stripes is either aresult of the greater range of unstable modes available, or of the increased strength ofthe nonlinearities due to larger pattern amplitude, or both. In the absence of the gradedinitial state, 4.25, one obtains a number of striped domains with pattern defects.When a gradient is used in this way, the symmetry breaking capabilities of reaction-diffusion are no longer being used to the full extent. Rather, in respect of orientation,the reaction-diffusion dynamics is used as a selective amplifier of pre-existing order, aswell as to stabilize a particular wavelength in that order. Biological justification for theuse of gradients is given in chapter 5.One further use of monotonic gradients of reacting substances is in coercing a systemthat would otherwise favour a spotted superposition of modes into making stripes [78]. Aninvestigation of the influence of gradients of reactant concentrations on the propertiesof the brusselator formed a major component of the early part of this project. Thebrusselator with gradients in A and/or B and also in rates of the reactions was the modelsystem used. A large number of computations were undertaken to look for particularcombinations of parameters and gradients that will select stripes. Very few of thesecomputations lead to successful production of striped patterns, and as a result I beganto investigate other ways to produce stripes. One computation is shown in figure 4.26.Other runs produced patterns that are either completely striped, or a mixture of stripesand spots. In many cases a gradient induced no visible striping tendency in the patterns.The conclusion is that this method is not, it itself, sufficiently robust to be a satisfactorystriping mechanism. A possible compromise is the use of gradients on models in whichChapter 4. Numerical Solution of Reaction-Diffusion Equations 84the quadratic nonlinearities are small, to further bias them towards striped patterns.The general conclusion from the computations presented in this chapter is that the mostsatisfactory models for striped pattern formation in development are those in whichstriping is intrinsic to the model and gradients (polarities) are used mostly to orient thestripes properly in the organism.4.6.1 Animation of Pattern Formation SimulationsEarly on in the development of this research program, I decided that it would be worthwhile to be able to animate the dynamics of reaction-diffusion processes in order to gainintuitive insight into their workings. While the current primary interest is in steady statepatterns that are the end products of the patterning process, a feeling for how the systemarrives at this state was found to be valuable in 4eepeuing the understanding of nonlinearTuring models. The exact solution to the Turing equations provides us with only thegrossest picture of dynamics: spatial modes are amplified according to their wavelength,the fastest growing mode swamping the pattern after a short time. The nonlinear modelsexhibit a wide range of phenomena not included in this picture, including movement andannihilation of pattern defects, slowing of pattern growth, mode suppression, the initialonset of pattern.My first efforts to produce an animated sequence resulted in a program for one dimensional pattern formation for IBM compatibles. The program solves discretized reactiondiffusion equations (both implicit and explicit methods have been tried) and the X or1 (or both) profiles are viewed using a simple ‘XOR. animation technique. For a reasonable number of finite-difference nodes (20-40) the programs run quickly enough onan IBM AT with math coprocessor to provide a vivid picture of the behaviour of Turingmechanisms. The formation of stationary patterns as well as oscillating structures haveboth been studied in this fashion. Modified versions of this one-dimensional program alsoChapter 4. Numerical Solution of Reaction-Diffusion Equations 85provided the primary research tool for the early work on a project in the group to modelheart formation in the Axoloti salamander Ambystoma mexicana. The immediate visualfeedback strongly affected the working method on this project, and allowed one of thecollaborators, an experimental biologist1,to quickly gain an intuition for the functioningof the models.Animation of two dimensional simulations seemed out of reach with only a AT compatible. Fortunately, a new scientific visualization facility was made available in theUniversity Computing Services at UBC in the Fall of 1990. The hardware consists ofa Silicon Graphics 4D/25 workstation with Turbo Graphics and a video board. Usingthe Dataview software developed by John Hogg and Helen Salter of the UCS, I wasable to turn my computer ontput into a 15 minute SVHS video that includes several- pattern-forming sequcnces [91]. The video has already proved useful in showing colleagues, especially experimental biologists, what the time evolution of reaction-diffusionpatterns looks like. Figures 4.27 and 4.30 respectively show the early and late stagesof patterning in a computation. The early stages of patterning (up to 10000) show thefast development of short-range correlations in morphogen concentrations (due to localsmoothing arising from diffusion) followed by establishment of wavelength of patterns.The blobs quickly coagulate and thin out, eliminating regions of high curvature, to forma striped pattern. On a longer timescale (by about an order of magnitude) fluctuationsdne to added noise can be seen to drive defects to move towards each other and annihilate. The last defects are eliminated extremely slowly. In the video it is possible tosee the system make several attempts to annihilate the last defects. This underlines theimportance of continuing stochastic input in the dynamics of pattern formation.1Dr. John Armstrong, of the University of Ottawa.Chapter 4. Numerical Solution of Reaction-Diffusion Equations 86Figure 4.15: X and Y morphogen concentrations in a computation with the Prigoginebrusselator. Times are in numbers of iterations (timestep = 0.002). Domain is a unitsquare 33 cells long and 33 cells wide. Parameters used were: A = 2.0, B = 4.17,4.831832 X i0, D = 10D.x,pt=o Yx,t=1000 YChapter 4. Numerical Solution of Reaction -Diffusion Equations 87X,t=3000 Yt=20000 YFigure 4.16: Figure 4.15 Continued.Chapter 4. Numerical Solution of Reaction-Diffusion Equations 88Figure 4.17: X and Y profiles in a computation with the hyperchirality model (see section2.6.2). Conditions are the same as in figure 4.15. Parameters used were: Px = 1.0,Py 0.52125, A = B = 1.0, ky 4.17, = 7.673861, 8.425, and diffusivitiesas in figure 4.15. These parameters were chosen so that the linear behaviour is identicalto that for figure 4.15.x,t=o Yx,t=1000 Ycjq CD rj ‘-1 CD C, 0 CDCj0 0 0 CD 0>< II C)0 0 0II F’)0 0 a 0 -<Chapter 4. Numerical Solution of Reaction-Diffusion Equations 90Figure 4.19: Time evolution of X, for (a)-(d) brusselator and (e)-(f) equations 4.33,with kq = 0. Parameters chosen so that the two systems have the same linearization,and differ oniy in the non-linearities. Computations performed close to the onset of theTuring instability on a grid 57 cells long and 33 cells wide with cyclic boundaries. FromLyons et. al. [89]Chapter 4. Numerical Solution of Reaction-Diffusion Equations 91Figure 4.20: Pattern of X after 30,000 iterations with equations 4.33, showing the effectof increasing the magnitude of the quadratic non-linearity. Computations performed ona unit square 60 cells long and 60 cells wide with cyclic boundary conditions. Parameterswere: k1 = 0.5, k4 = —1.8, k6 = k7 = 1, D = 1.327 x i0, D = 10D, timestep=0.04.kq=O.1 kq=O.5kq=O.777 kq=1.OChapter 4. Numerical Solution of Reaction-Diffusion Equationsad e92Figure 4.21: Striped patterns observed in computations with the reaction-diffusion systemin equation 4.33. Diffusivities in the bottom three are half those in the top three, butparameters were otherwise the same for the six simulations.b CfChapter 4. Numerical Solution of Reaction-Diffusion Equations 93Figure 4.22: Effect of a monotonic gradient in the initial X concentration on patternselection. Calculation with equations 4.33.t=0 t = 4000t = 10000 t=20000Chapter 4. Numerical Solution of Reaction-Diffusion Equations 94t = 10000Figure 4.23: Calculation exactly as in figure 4.22 but without a gradient in the initial Xconcentration.t=0t=20000 t=30000Chapter 4. Numerical Solution of Reaction-Diffusion Equationst = 1000095Figure 4.24: Computation similar to that given in figure 4.22, but paramaters adjustedso that the reaction-diffusion system is further away from the onset of the diffusiveinstability.t=0t=20000 t=30000Chapter 4. Numerical Solution of Reaction- Diffusion Equations 96Figure 4.25: Calculation exactly as in figure 4.24 but without a gradient in the initial Xt=o t= 10000t=20000 t=30000concentration.Chapter 4. Numerical Solution of Reaction-Diffusion Equations 97Figure 4.26: Computations with the Prigogine brusselator (a)-(c) with a linear gradientin reactant A from the top to the bottom of a vertical right circular cylinder; and (d)-(f)with no gradient present.t=0 t=8000 t=1 6000t=0 t=8000 t=1 6000Chapter 4. Numerical Solution of Reaction-Diffusion Equations 98Figure 4.27: Early stages of patterning in a computation with equation 4.33. Theseimages are frames from a video produced with the assistance of the UBC UniversityComputing Services. [91].t=o t=1000t=2000 t=3000tj CD C) 0 p 0 0 CDC) p (r CD p Cf) 0 CS 0• 0 CD p 0II 0) 0 0 0 II o o 0II 0 o o ‘-4- II 0i0 0 0CD CDh!j-4.CD C) 0 0 0 (00 0 0 0II 0 0 0-I II (0 0 o 0LI) 0• 0 -‘-4CD C-’0• IChapter 4. Numerical Solution of Reaction-Diffusion Equations 101t=o t=10000t=20000 t=30000Figure 4.30: Later stages of patterning in the same computation as figure 4.27.cm c.z C) 0 -4. 0 0 cmCb cm C’)0 4— 0 0 Cb C) 0•II 0) 0 0 0 II 0 0 0 0‘-I II 0 0 0 0 II (ii0 0 0 0rj (D Cj C C 0 og CDa 0 a a aII 03 0 0 0 a II (0 0 a 0 0CD 2 CD C.) et 0• 0 CD C.) 0 bChapter 5Biological Examples: Drosophila Segmentation and Ocular DominanceInsects and other arthropods as well as annelids have a body plan that is organized ona pattern of repeated units or segments. In vertebrates, the somites of the backbone areanother instance where a sequential repeat of similar structures provides the foundationof a complex body. To reaction-diffusion modellers, such patterns are strongly suggestiveof Turing structures.The discovery of precisely repeating patterns of proteins and messenger RNA alongthe antero-posterior axis of the Drosophila melanogaster embryo (for example, the pair-rule genes are expressed in seven regularly-spaced bands encircling the egg), raised theexciting possibility that reaction-diffusion might be operating at a crucial, early stage ofdevelopment. These patterns of chemical substances lead directly into the morphologicalsegmentation, and therefore seem to be in the category of chemical prepatterns whichTuring anticipated. For reviews see [1, 65].Drosophila melanogaster is a convenient choice for genetic work, and has become themost popular organism for experiments on the control of biological pattern. Since theearly eighties, there has been an overwhelming number of experimental papers publishedon this topic.A definitive review article on early Drosophila development published in 1987 byMichael Akam [1], suggested that the segmentation gene products known as the pairrule genes might be part of a reaction-diffusion system. In several places, the reviewsuggested reaction-diffusion mechanisms as a likely explanation for the pair-rule gene104Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 105expression patterns. The following excerpt is typical:Initially the transcription of these genes is either uniform, or simply localized,bnt complex periodic patterns rapidly evolve. These patterns characteristically exhibit a stage showing a double segment periodicity. Local order in thisevolving pattern probably depends on some reaction-diffusion mechanismTwo years later the same anthor wrote a commentary in the magazine Nature [2] with thetitle ‘Making Stripes Inelegantly’ stating that the primary pattern forming event in thepair-rule system was probably not reaction-diffusion after all, but static gradient reading.Of the reaction-diffusion mechanisms Akam wrote [2]:One very likely role for these processes is to amplify the initial discontinuitiesin pair-rule gene expression, and to sharpen the boundaries between stripes.These two statements represent a change in the views of many workers on Drosophiladevelopment that took place over a short period of time with little experimental evidencefalsifying the reaction-diffusion hypothesis or supporting the alternatives. The application of reaction-diffusion theory to the control of Drosophila segmentation mnst thereforebe regarded as controversial. In light of this situation one strategy is to analyze commontheoretical problems faced by reaction-diffnsion mechanisms of segmentation in general,with no specific reference to the molecular details of how the mechanisms are implemented. The hypothesis that reaction-diffusion is involved in the early embryogenesis ofDrosphila faces two strong theoretical challenges: the mechanism underlying segmentation must be able to form striped patterns on wide domains in a robnst fashion, and itmnst be able to orient the stripes along the antero-posterior axis of the egg. I discusshow the results of chapter 3 and 4 demonstrate that this may be accomplished using areaction-diffusion mechanism.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 1065.1 Background Material on DrosophilaDrosophila melanogastcr is a holometabolous insect; that is, there is a complete metamorphosis separating the larval and adult stages of the life cycle . A fertilized egg, laidon a suitable substrate such as rotting fruit, undergoes embryonic development for abouta day and then hatches as a worm-like larva which feeds on the surrounding medium.The larva grows, and molts several times as it increases in size. The periods of growthseparating molts are known as the instar stages. After the third instar stage a pupaforms and metamorphoses to produce an adult fruitfly. During metamorphosis most ofthe larval body is consumed by the formation of the adult from the imaginal disks, smallpatches of cells set aside during the development of the larva. The adult appendagesfor mating and reproduction (such as wings, legs, genitalia, and organs of the head)are formed through a process which somewhat resembles the inflation of the balloon-likeimaginal disks. The imaginal disks were widely studied as a model developmental systemin the seventies [20], and a reaction-diffusion theory for some of the disk phenomena wasproposed by Kauffman [69]. Figure 5.33 shows the body of the Drosophila melanogasterlarva and adult fruitfly, which consists of a head region, three thoracic segments, to whichlegs, wings, and halteres are attached, and the abdominal segments, containing digestiveand reproductive organs. Unlike the individual organs, the segmentation pattern persiststhrough metamorphosis.The segmented morphology of the embryo does not become apparent morphologically,until more than five hours after fertilization, well after gastrulation2 in the embryo hastaken place. Recent work, using fluorescent antibodies against a-tubulin, has suggestedthat the distribution of microtubules early in gastrulation is not uniform, but displays1111 contrast to a hemimetabolous insect where the adult structures are formed by modification of theexisting larval ones.2One of the earliest and most fundamental symmetry breaking events in development in which a gutis formed.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 107Figure 5.33: Segmented organization of the Drosophila body plan, at both the larval andadult stages, from Gilbert [32].evidence of metameric organization [15]. However the earliest indications of segmentation in Drosophila are first visible, in spatial patterns of gene activation, much earlier,before the onset of gastrulation. The early embryogenesis of Drosophila is somewhatunusual. The first 13 nuclear divisions are synchronous and take place in a single cytoplasm (syncytium). At about interphase 10 (after the ninth nuclear division) the nucleimove from the center of the syncytium to just inside the peripheral or vitelline membrane surrounding the egg. This stage is called the syncytial blastoderm, in reference tothe single two-dimensional layer of nuclei in diffusive contact with the single cytoplasm.After further nuclear divisions the vitelline membrane moves in enclosing each nucleusto form the cellular blastoderm. Interphase 14 (between the 13th and 14th nuclear divisions, occuring at about hour 3 in development) takes about half an hour. This is longerthan previous cell cycles. It has been established [16] that the developmental fate of eachAbdominalsegmentsLarval AdultChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 108nucleus is determined sometime just around the cellularization of the nuclei. Interphase14 is thus regarded as a critical period during which many of the later events are determined. Gastrulation begins just before separation of the nuclei from the yolk is complete,suggesting that the important spatial patterns must be set up while the embryo is stillin a syncytial state [28].It is during the last few nuclear cycles before cellularization that the first zygoticgenes are transcribed into mRNA and translated into protein. The three classes ofzygotic genes first activated in the embryo are the gap, pair-rule, and segment polaritygenes. Collectively, these are known as the segmentation genes, as mutations affectingthem cause deformity in the metameres of the larva. Zygotic genes are distinguishedfrom the maternal genes. mRNA transcripts of maternal genes are deposited in theegg during oogenesis. The order in which genes become activated in the Drosophilaembryo is roughly maternal, gap, pair-rule, and segment polarity. The expression ofthese genes forms a regulatory hierarchy. For example expression of the pair-rule genesis affected by the maternal, gap, and other pair-rule genes; gap gene expression dependsonly on maternals and other gaps; segment polarity genes may be affected by all of theabove and other segment polarity genes. The segmentation genes were first described ingenetic experiments with mutants displaying abnormal metameric pattern at the larvalstage [109]. Gap mutants such as Hunchback and Knippel were found to lack a number ofcontiguous segments in the embryo. The pair-rule mutants, in which even or odd segmentsare missing from the 14 segment pattern, were found to occur in complementary pairs.The segment polarity mutants describe larvae where part of each segment is deleted andreplaced with a mirror image of the remaining part.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 1095.2 Reaction-Diffusion Models of Drosophila SegmentationThe first suggestion in the literature that segmentation in Drosophila may involve areaction-diffusion mechanism is due to Kauffman et. al. [69]. Kauffman extended hisanalysis of compartmental boundaries in ellipsoidal imaginal disks to pattern formation in the syncytial blastoderm. Kauffman’s work was based on the linear analysis ofreaction-diffusion equations. The model was severely criticized by Bunow et. al. [12]who suggested, on the basis of computer experiments with Kauffman’s equations, thatthe required sequence of modes, forming circumferential bands on the embryo, wouldbe unstable. Bunow’s statements are correct if restricted to the model considered byKauffman. However, as is shown analytically in chapter 3 and numerically in chapter 4,reaction-diffusion models are capable of stabilizing striped modes.Several other models of Drosophila segmentation exist. Meinhardt’s work [94, 96, 97]has contributed mostly to understanding of the role of maternal gradients, and to thehierarchical nature of the gene regulation. Nagorcka [106] presented a detailed descriptionof how a spatial frequency doubling reaction-diffusion mechanism can provide all of thespatial information needed to control expression of all the segmentation genes. A criticismof Nagorcka’s approach is that it relies on an unidentified pair of morphogens whichcontrol all gene expression, so it seems to ignore the regulatory hierarchy suggested byMeinhardt, and which has been borne out by experiment. Another serious criticism,and one which applies also to a paper by Goodwin and Kauffman [34] is that too muchemphasis is is given to the phenomenon of spatial frequency doubling in the constructionof the model. The full frequency doubling sequence has only been observed for thepattern of ftz3 expression. As well, there is only one experimental paper mentioningthis phenomenon [143]. A positive feature of Nagorcka’s work is that he uses a model3ftz=fushi tarazu, a pair-rule gene (Japanese for ‘not enough segments‘.)Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 110with strong intrinsic striping capabilities. Lacalli [78, 80] has considered two approachesto Drosophila segmentation. Using a brusselator he showed that stripes can be formedusing monotonic gradients in reaction rates. He has also used a dual brusselator modelfor two pairs of the pair-rule network. Hunding [61, 63], following Lacalli, uses monotonicgradients to force striping. An unusal feature of Hunding’s work is that his simulationsare on spheroidally shaped domains, reflecting the shape of the egg. This is interesting,but may be a premature fine-tuning of the model. Finally, there is more work underwayin Harrison’s lab. Harrison and Lakowski have examined localization, maternal and gapinput to pair-rule patterning, simulating the effects of some mutations. Lakowski is alsopreparing a phenomenological study of the effect of maternal gradients on localization ofgene expression zones.The most significant theoretical challenge confronting the reaction-diffusion models ofDrosophila segmentation is the problem of forming regular striped patterns on non-narrowdomains. The popular brusselator and Gierer-Meinhardt models both have a strongtendency to make spotted patterns on domains wider than the chemical wavelength. Thesyncytial blastoderm is approximately 420 pm long, and has an approximately circularcross-section with 150 pm [128]. This gives an approximate circumference of 470 pm.Only about 70% of the egglength is striped, pair-rule expression being repressed in themost anterior and posterior regions of the egg. If the patterning domain is cut and rolledfiat, it gives roughly a rectangular domain that would be about 300 pm long and 470pm wide. The wavelength of a seven stripe pattern on such a domain is about 40 pmwhich is comparable to what may be seen in the ft expression pattern (figure 5.34). Thestripe-forming mechanism must be robust to minor perturbations. For example, spottedpatterns of gene expression are not seen.Models with anti-symmetric net production functions are certainly sufficiently robuststriping systems to explain the pair-rule patterns. There is no experimental evidenceChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 111Figure 5.34: Drawing of the pattern of expression of the segmentation gene ftz at inter-phase 14. Dark bands represent areas in which mRNA transcripts of the gene may bedetected.for or against the existence of such a symmetry in the Drosophila pair-rule system. Anatural way to obtain such a symmetry is to use a four-morphogen mechanism such as thehyperchirality (section 2.6.2) or Lacafli’s dual brusselator model (equations 2.90 of section2.6.1) with symmetrica.1 rate constants. Numerical results with equation 4.33) indicatethat the constraint on the non-linearities of a suitable reaction-diffusion mechanism is infact weaker: small quadratic perturbations to the equations do not affect the selection ofstriped patterns.An alternate mechanism for the selection of striped patterns by reaction-diffusion isto use monotonic anteroposterior gradients in reactant concentrations or reaction rates.Lacalli [78] was the first to demonstrate this. His work used the Tyson brusselator, andChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 112gradients in reaction rates. I have extended his results, showing they apply also to gradients in reactant concentrations for the Prigogine brusselator [88] . My experimentationwith this approach showed that it depends on a fine-tuning of the model parameters andgradients. More importantly, it becomes difficnlt to get this approach to work as thedomain width is increased. The conclusion is that this is not a viable way for Drosphilato make stripes. It seems necessary that if a reaction-diffusion process is working in thepair-rule system, then the reactions must have an intrinsic stripe-selecting character.The pair-rule stripes are free from all gross pattern defects (branching, spotting, etc.)and are aligned roughly perpendicular to the antero-posterior axis of the egg. It wasshown in chapter 4 that this can be accomplished using a monotonic gradient in theinitial state of the morphogen. Such a gradient could arise from a gradient in reactantsor reaction rates in a reaction-diffusion mechanism. Such monotonic gradients are knownto exist in developing systems [32]. The best characterized example of such a gradientis the protein product of the maternal gene bicoid in Drosophila [21, 22]. This gradientalmost certainly arises by a simple one component reaction-diffusion process itself. Thebicoid protein concentration has been quantified, and is accurately fitted with a singleexponential [21]. This is typical of the diffusion- first order decay system mentioned inchapter 2. The bicoid gene has been demonstrated to have a strong influence on theexpression of gap and pair-rnle genes [22]. Thus it seems reasonable to conclude that thisgradient could provide the necessary cue for formation of orderly and oriented pair-rulestripes.Bicoid also appears to play a role in regulating the anterior border of the pair-ruleexpression zones [22, 29, 110]. Pair-rule patterns extend over only about 70% of theegg length, in a region in the interior of the embryo. This raises the issue of how thepatterning process is deactivated to localize gene expression. A possible solution to this4Also, M. J. Lyons, unpublished results.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 113problem is provided by work of Herschkowitz-Kaufman [52, 51] on localized dissipativestructures. Herschkowitz-Kaufman showed, using one dimensional simulations, that adissipative structure may be confined to an interior domain of a reactor by using agradient of reactant concentration to take the system into and back out of the patterningregion of parameter space as the domain is traversed [49]. Remarkably, this conceptwas illustrated by localizing a seven peak wave using a hyperbolic cosine gradient ofbrusselator reactant A. Harrison5 has shown how this result may be applied to the pair-rule patterns in Drosophila. I have confirmed this work numerically in one-dimensionalcalculations with the brusselator. The idea is to consider bicoid as either a reactant orrate-constant affecting factor which affects the Turing constants. Bicoid concentrationstoo high or too low shift the reaction-diffusion parameters out of the morphogeneticregion f parameter space. The concentration of bicoid is highest in the anterior of theegg, where no pair-rule activation is seen. De-activation in the posterior of the eggrequires an analogous posterior- high gradient. Such a gene complex is known, but it hasnot been characterized as extensively as bicoid.5.3 Ocular Dominance StripesThe Drosophila pair-rule gene network provided an example of a striped pattern whereregularity of the stripes was of critical importance to the correct development of thebasic body plan. Biological examples of striped patterns not requiring a high degreeof regularity are also numerous. Fish display a wide variety of patterns that includeregular and irregular stripes in various orientations (antero-posterior and dorso-ventral).Zebra coats exhibit the types of topological defects mentioned in chapter 4. Anotherfamiliar example is the patterns of ridges found on the skin of the fingertips. All of these5Rarrison and Lakowski, unpublished results, see also [88].Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 114patterns contain defects which can be described qualitatively as blind alleys, Y and Hbranchings, and loops. Such defects were observed in the two dimensional computationswith reaction-diffusion, where it was found that the particular configuration of defectsof a calculation depends on the stochastic input during the simulation. In this regard,it is interesting that fingerprint patterns have been found to be unique to an individual;even genetically identical individuals (mono-zygotic twins) have very different fingerprintpatterns [86]. At the end of the nineteenth century, a suggestion that fingerprints be usedfor evidence in paternity suits was rejected because of the lack of a heritable componentin the patterns. The genetic control of other patterns such as the zebra stripes, orocular dominance columns is an interesting question but more difficult to study. Theuniversal features shared by the above instances of biological stripes may suggest similardynamical origins. As the computations have shown, reaction-diffusion is certainly acandidate mechanism for any of the examples in this chapter.Ocular dominance stripes (see figure 5.35) are spatial patterns of synaptic connectionsobserved in the visual cortex of some mammals, particularly primates [57]. The stripesbear a strong similarity to zebra stripes. In the following, I will consider an existing classof models [98, 129] for the development of ocular dominance patterns. The nonlinear behaviour of these models is studied using the approximate analytical technique of chapter3.In mammals, visual information from the ganglion cells of the retina passes alongthe optic nerve through the optic chiasm to the lateral geniculate nuclei (LGN), thenonwards to the primary visual cortex (also known as the striate cortex or area 17). Seefigure 5.36 for a schematic illustration of the pathway linking the retina to the primaryvisual cortex. At the optic chiasm there is a crossing over of fibers from the left andright eyes so that each LGN receives afferents from both eye. From there all of the fibresefferent from the left LON synapse at the cortex in the left hemisphere, and similarly forChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 115Figure 5.35: Ocular dominance stripes in layer JVc of the macaque monkey primaryvisual cortex, as they appear in an autoradiograph of the hemisphere contralateral to theeye injected with tritiated proline. Drawn from figure 23 of reference [57].the right LGN. There is a single axon between the retinal ganglion cells and the LGN,and also between the LGN and cortex, so that there is only one synapse separating thecortex from the retina. Most of the terminals synapse primarily in a single anatomicallayer of the striate cortex known as layer IVc.Much of what is known about the cortex arises from observing the activity of singlecells with a microelectrode applied to the exposed brain of an experimental animal, oftena monkey or a cat. The response of cells to controlled stimuli may then be studied.In this way it was found that the visual input from the two eyes is initially segregatedinto distinct non-overlapping but interpenetrating regions of the cortex. These regionscould be mapped out somewhat laboriously using the electrophysiological technique. AnChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 116Figure 5.36: Schematic of the visual pathway in the human brain, adapted from [58].alternate, faster, method of mapping these regions of left or right ocular dominance wasfound and applied to more extensive studies of these patterns. The method involvesinjecting one eye with a radioactive amino- acid, which is transported along the axonsof the optic nerve through the synaptic junctions at the LGN and along axons to thevisual cortex. An autoradiogram of the cortex is then taken, recording which areas ofthe cortex received the labelled amino acid. In macaque monkey the resulting patternof ocular dominance (figure 5.35) was found to be irregular stripes with many branches,but with a uniform spacing over the entire extent of area 17. In the monkey the areasChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 117corresponding to terminals from the ipsilateral (same) or contralateral (opposite) eyewere fonnd to be equivalent, i.e. equal in form and area. Approximately equal numbersof synapses to layer IVc are from the left and right eyes. The patterns recorded using theradioactive tracer were found to correspond to the patterns found using electrophysiology.In other layers of the cortex there are still regions where one eye predominates over theother, but the segregation of left and right is not so complete.The pattern of ocular dominance is found to be absent until the last few weeks beforebirth; layer fl/c consisting initially of an unpatterned mixture of ipsilateral and contralateral terminals. The pattern continues to evolve during the first few weeks followingbirth, and then remains more or less constant. The proper development of the stripesdepends on visual input during this period. In monocular deprivation experiments in themacaque monkey, the visual input to one eye is blocked. Inthe cortex the number ofsynapses from the visually deprived eye is diminished and the area of cortex receivinginput from that eye shrinks. Under these conditions the pattern no longer takes theform of stripes but rather resembles a spotted pattern of ipsilateral ‘islands’ in the sea ofcortex with contralateral input. Binocular deprivation gives the paradoxical result thata striped pattern of ocular dominance is again seen.The other widely studied animal is the cat. With the cat’s visual system it is foundthat the numbers of terminals reaching the cortex from the two eyes are unequal. Approximately 70% of the synapses in layer fl/c are from the contralateral eye. A spottedpattern of ocular dominance is seen. The terminals are not as sharply segregated as inthe monkey. However, there are still distinct regions connected predominately to the leftor right eyes. The ipsilateral-dominated regions are islands in a sea of contralaterallydominated regions (see Figure 5.40)Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 1185.4 Models for Ocular Dominance Pattern FormationTwo widely-known models for ocular dominance patterns [98, 129], take the form ofintegro-differential equations. Swindale’s phenomenological work [129] showed by computation that the general intracortical interactions proposed by Hubel and Wiesel [57]were capable of generating the observed patterns. The model also provided a frameworkfor testing the effects of different forms of intracortical interaction that might result fromsight deprivation experiments. Miller has [98, 99, 100] developed a model with a similarmathematical structure to Swindale’s, but which includes greater physiological detail.The two models take the form of integro-differential equations. As an example I firstconsider Swindale’s equations [129]. The densities of left and right synapses nL(P, t) andnR(r, 1), which are functions of a two dimensional spatial coordinate P on the cortex andtime, develop according to:I—, ,= ffra) J (na(r — r )WRR(r ) + nL(r — r )wLR(r fldr= f(nL)J(nR(P— P’)wRL(P’) + flL(r P’)WLLQ’))dP’. (5.1)The integrals sum competitive and cooperative interactions between like and oppositeterminals over the entire patterned region, weighted by the intracortical interactionsw(P’) (shown in figure 5.37). The function f(n) takes into account the nonlinear effectslimiting the growth of the synaptic densities. It is common to take a simple logistic formf(n) = n(N — n) where N is the maximum density of synapses. Equations 5.1 may bereduced to a single equation for the difference between left and right densities with theassumption that the total density of synapses N = L + fl is constant [129]:=(a*w+K)(N—a)(N+a). (5.2)a=— L, in = 1/4(WRR + WLL), K (N/4) * (WRR — tOLL) and * indicates convolution. Equation 5.2 represents a minimal model for OD stripes which has been shown byChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 119Figure 5.37: Typical weighting functions (kernels) for interactions between likesynapses (w6 = WRR, WLL, solid line) as well as between terminals from opposite eyes(w0 = wLR, WRL., broken line). Pattern formation in this model resnits from short rangeactivation and long range inhibition of like synapses, and vice-versa for the cross interactions (after [129]).computation to possess the essential patterning capabilities exhibited by more complexmodels [129].5.5 Non-linear AnalysisI shall consider models with symmetry between left and right terminals, i.e. WRR = WLLso K = 0. This implies that equation 5.2 is unchanged under the operation a — —a. Ianalyze the simplest possible pattern-forming system having this symmetry [135]:8a 3= a * w — awSweSSSSS4SSSSSSISSISSSISS•0(5.3)The coefficient of the cubic term can be set to unity without loss of generality by changingthe dimensions of the ocular dominance a. Other constants may be absorbed into w.5.3 has homogeneous steady state solution given by a = 0 which corresponds to anChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 120unpatterned state consisting of an equal mixture of left and right terminals. Close tothis state the linear part of the equation predominates:ãa—=a*w. (5.4)Equation 5.4 is exactly soluble by Fourier analysis. For the moment, we treat a finitedomain, but otherwise leave the boundary conditions unspecified. The general solutionis a superposition of plane waves with amplitudes 4k(t)a(r,t) == 4)eWt. (5.5)W(k) is the Fourier tranform of the kernel w(r) which depends only on the magnitudeof Z, and the o(Z) are to be determined from the initial condition.To obtain an approximate solution to the nonlinear equations 5.3 continue to writethe solution as a sum of plane waves 5.5, and derive differential equations for the modeamplitudes. Substituting 5.5 into 5.3, multiplying on the left by 6ic9 and integrating overthe domain using the orthogonality of eigenfunctions:fe’)dA = 6, A = domain area, (5.6)one obtains the mode amplitude equations:= W(k)— (5.7)k ?c “lv “A similar set of equations was derived for reaction-diffusion models in chapter 3. Equation5.7 is equivalent to the original form 5.3, and is intractable. It is possible to understandsome of the behaviour of model 5.3 by studying it with a special form for w(r) whichallows an approximate solution of equations 5.7 by the adiabatic method (or slavingChapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 121principle). The conditions for the approximation to hold are:W(k) > 0 for one value of kW(k) < 0 forkk,W(k) >> IW(k0) . (5.8)If these conditions are satisfied then the relaxation times of the stable modes are muchshorter than that for the unstable modes. The stable mode amplitudes then adjustrapidly to the current values of the In the adiabatic approximation all of the stablemodes (k k0) are then eliminated from the dynamics leaving a set of ODEs for theunstable mode amplitudes, k = k0 at all possible relative orientations gl:= W(k0)-- 0< <K,. (5.9)VThis is essentially the same set of mode amplitude equations found for the symmetricalreaction-diffusion systems discussed in chapter 3. This result partially explains the closesimilarity in patterns produced by the two mechanisms (see figure 5.38).In the macaque monkey, ipsilateral and contralateral afferents to the cortex (i.e.connections from right or left eye to the same or the opposite side of the brain) areapproximately equal in numbers. Also, left and right dominance stripes in layer IVcare similar in form and width (Figure 5.39). The only difference between left and rightdominated bands appears to be their connection to different eyes. This symmetry betweenipsilateral and contralateral inputs is related to the selection of stripes rather than anyother spatial pattern. In monocular deprivation experiments, this symmetry is disrupted,leading to a spotted pattern of ocular dominance in which terminals from the eye receivingvisual information are more numerous than ones from the blindfolded eye. Binoculardeprivation restores the left/right symmetry and stripes are observed when neither eyereceives input during development.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 122II;. .IWd c!+_r (b) t=3000(C) t=5000Figure 5.38: Time evolution of the concentration proffle for one of the morphogens ina reaction-diffusion system with the dynamical symmetry discussed in chapter 3. Thedomain is square with no-flux boundary conditions to mimic the behaviour observed atthe boundary of area 17 of the cortex where ocular dominance stripes run perpendicularto the edge.(d) t=15000Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 123Figure 5.39: Pattern of ocular dominance patches in the visual cortex of a monkey.Adapted from a digitally scanned image of figure 5 of [84].Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 124In the cat, there is normally a quantitative asymmetry between the ipsilateral (about30%) and contralateral (about 70%) inputs to the cortex [83]. A splotchy pattern ofocular dominance is found (figure 5.40). Binocular deprivation gives a pattern quiteunlike stripes. It may be possible to study the effect of varying degrees of asymmetryon the ocular dominance patterns through a comparative survey of various species ofmammals with different ratios of ipsilateral/contralateral afferents.Chapter 5. Biological Examples: Drosophila Segmentation and Ocular Dominance 125Figure 5.40: Pattern of ocular dominance patches in the visual cortex of a cat. Adaptedfrom a digitally scanned image of figure 7b of [3].Chapter 6ConelusionI conclude this thesis with some suggestions for further study.There is no shortage of things to try with non-linear pattern formation models. Asdiscussed in chapter 3 non-linear equations are almost always intractable. The analysispresented in chapter 3 is valid in only a limited region of parameter space, and must beviewed as a first step beyond the linear analysis in understanding the properties of themodels. Other regions of parameter space have been explored using computer calculations, but unless such experimentation is guided by some form of analysis, it is blind. Along term project is to study how patterning modes are coupled away from the onset ofthe instability. At present it is not clear how such a study would proceed. Experimental or computational approaches will improve significantly as increased computation andscientific visualization tools become generally available.With respect to the study of real world instances of Turing structure several commentsare in order. The recent demonstration of Turing structures in vitro [Castets, Tabony],made 1990 an exciting year for reaction-diffusion modellers. Certainly, it would be valuable to attempt detailed and realistic models of the experiments. In the case of Castetsresults such a project has already begun and preliminary results have been published[82]. Several important aspects of the structures have not been explained. The model isonly one-dimensional, and the real patterns are two-dimensional. The knowledge of thereaction kinetics remains incomplete, and will first require experimental elucidation ofthe underlying mechanisms of pattern formation. Further, more realistic representation126Chapter 6. Conclusion 127of the molecular transport than linear Fickian diffusion will probably be reqnired. Modelling, and analysis which could be based on the work presented in this thesis will playan important role in the process of refining our understanding of these chemical patterns.The microtubular dissipative strnctures [131] are even more poorly understood. Furthermore, it seems that there have not been any attempts to model the phenomena. Dr.J. Tabony has written to me stating:We believe that the underlying mechanism for the patterns that we observedis, as you suggest, a reaction diffusion process. ... We also hope by changingconditions. to modify the pattern from stripes to spots. ... I am not aware ofanyone else working on the theory. ... I believe that a model of the type youdescribe provides the underlying explanation for our observations.It is possible to conclude that the striped patterns observed in Tabony and Job’s experiments are a consequence of an underlying symmetry in the dynamics, which mustbe significantly perturbed to observe a spotted pattern. A promising line would be toattempt to develop a realistic model of Tabony and Job’s microtubular solutions, withthe aim of suggesting modifications of their design which would allow control of patternselection.The application of reaction-diffusion to biological problems has had a difficult past,and faces a difficult future. It seems that this is independent of its value as a scientific theory, but rather a sociological consequence of the current state of the researchcommunity in biological development. Most experimental developmental biologists (particularly the geneticists) are unwilling to consider partial differential equations as beingrelevant to their subject. This may be a fear of the mathematics involved, or perhapsa clash between different preconceptions of the nature of the scientific enterprise [49].The situation is unhealthy for both camps, making it harder for theorists to interactChapter 6. Conclusion 128with experiment. The models of experimental geneticists of the patterning in Drosophilaare often inconsistent at best, and at worst fanciful. It is hoped that the situation willimprove in the future, and that dialogue between experiment and theory improves.In the case of ocular dominance models, testable predictions were suggested in chapter6. Synaptic plasticity is currently a topic of active research. Knowledge of the mechanisms by which synapses are weakened and strengthened will allow more accurate modelsof the developmental process.Appendix ATuring’s Two Cell ReactorTuring introduced the concept of diffusion induced instability with a qualitative discussion of the phenomenon in a system consisting of two simple chemical reactors (oridealized ‘cells‘) coupled by linear exchange of matter. The two cell system illustratesthe essential nature of the Turing instability in a simple way, and therefore has pedagogical value. The analysis shares features with other reaction-diffusion systems. Unlikethe one-dimensional ring of cells or continuous tissue, the dynamics of the two cell system can be represented using a two-dimensional phase plane. This provides a visualpicture of the Turing instability which generalizes to the more complex cases. Finally,for certain nonlinear models (specifically the brusselator) there is an exact solution forthe inhomogeneous steady states for these simple systems.A.1 Solution of Turing’s ExampleConsider two chemical reactors (labelled i = 1, 2) with chemical reactions such that thereis a steady state X0,Y0, in the concentrations of two chemical intermediates, X and Y.Let X and Y also stand for departures of the concentrations of substances X and Y fromthe steady state value. Then:= k1X-I- kY= k3X + k4Y (A.1)129Appendix A. Turing’s Two Cell Reactor 130is the most general form for the leading order terms describing the dynamics aronnd thesteady state. The cells are identical as far as chemical reactions, temperatnre, and otherphysical conditions are concerned, bnt may differ in the concentrations of chemical Xand Y in each cell. That is, the dynamics in cells 1 and 2 are identical, although thecells may be in different states.Turing discussed the specific case where parameters take values:k1=5, k2=—6, k3—6, k4=—7 (A.2)The cells are coupled with a linear transport law. That is, the cells exchange X, and Yat ratesDx(X2—Xi),Dy(Y2—1). (A.3)We take:Dx = 0.5,D = 4.5. (A.4)Then the initial value problem takes the form:5 —6 0 0Y1 6—700(A.5)X2 005—6 X2Y2 006—7 Y2in the absence of diffusion, with X(t = 0) = X0 and 11(t = 0) = Yo. In the presence ofdiffusion:4 —6 0.5 0 X16 —16 0 4.5 Y1(A.6)X2 0.5 0 4 —6 X2Y2 0 4.5 6 —16Appendix A. Turing’s Two Cell Reactor 131In the absence of diffusion the solution is asymptotically stable. This may be seen byconsidering the solution to A.5:5 —6 (A.7)6 —7has a single eigenvalue -1 and a single eigenvector:1iY = (A.8)1We look for a solution of the form:1 1 1’et c c2i7+ c2kt (A.9)1 1)Taking:1(A.1O)—1and k is found from:(A—I)i7= kv. (A.11)The general solution to A.5 is:X(t) 1 1 1= et c1 + c2( + 12t (A.12)Y(t) 1 —1 1)The values of e1, c2 are determined from the initial conditions:= (X + 1’)c2 = (X0 — Y0) (A.13)Then:x(t) = e—t(X + 6(X0 — Y0)t) — 0y(t) =et(Yj + 6(X0 — Y0)t) —* 0 , (A.14)Appendix A. Turing’s Two Cell Reactor 132while:X(t) C+Dt= -÷1 (AJ5Y(t) A+DtThe state of each cell approaches (0, 0) after some time, along the vector (1, 1), so thedynamics becomes very simple after an initial transient period.When the cells are coupled the behaviour is very different. The fixed point (0, 0, 0, 0)is unstable to arbitrarily small perturbations in the initial concentrations of the morphogens, and there is an exponential divergence from the initial condition. First, changing to a more convenient set of variables:X = X1+X2,Y=Yx = X1—X2, y=Y1— (AJ6)the equations for the new variables are uncoupled:X 5 —6 XY 6-7 Y4—6 x=. (A.17)6 —16 yIt was shown above that the first of these equations, which describes the time evolution ofthe total concentrations in the two cells, has a solution which asymptotically approachesthe point (0, 0). The total amount of X or Y in the system tends to zero. The secondequation is unstable and the fixed point at the origin is a saddle point of the dynamics.The matrix:4 —6(A.18)6 —16has two eigenvalues A = —14, 2 and corresponding eigenvectors:3 1and (A.19)1 3Appendix A. Turing’s Two Cell Reactor 133The general solution is:c1 e2t + c21e_14t (A.20)1 3where:3 1(A.21)ho 1 3The solution reads:xQ)= (3x0— ho) 62t+(3Y0 — xo) 1e14t (A.22)8 1 8The system possesses a positive eigenvalue and x and y grow with time dependence e2t.After a short time the state is well described byx(t)= (3xo Yo)e2t. (A.23)y(t) 8 1Phase planes are given for the uncoupled and the coupled cases (figure A.41). The essenceof the Turing instabiltiy is that in the presence of diffusion a previously stable systemmay become unstable. In this case a stable node became a saddle point when the cellsare coupled.Appendix A. Turing’s Two Cell Reactor>;-#4>-II‘I-#4II>1= x2-x1134Figure A.41: Phase plane representation of the Turing instability in a two cell reactor.A, absence of morphogen transport. B, Fickian transport of morphogens present.x = X2-X1Appendix BAn Exact Result for Coupled BrusselatorsUsing the symbolic algebra program Math ematica, I was able to show that there is anexact solution for the steady states of two coupled brusselators. states of the two cellbrusselator using the symbolic algebra program Mathematica. After discussions withWilliam Derrick1and Lionel Harrison, I decided to try to simplify the complex expressionsfor the steady state concentrations. The solution to the problem turns out to reduce tothe solution of a quartic polynomial which is a product of two quadratics. This resultprobably does not extend to other nonlinear two cell reaction-diffusion models. Howeverit is of interest, in view of the usefulness of the brusselator model.Consider reactors labelled 1 and 2 coupled by a linear exchange of matter througha semi-permeable membrane. Each cell contains the brusselator reactants at the sameconcentrations for both, i.e. A1 = A2 = A, B1 = B2 = B. The morphogens areexchanged at rates:— X1),Dy(Y — Y1). (B.1)X1,X2, 1<1, Y2 stand for concentrations of X and Y in cells 1 and 2. Then the reactiondiffusion equations to be solved are:= A—(B+l)X1+XYDx(X2X),Y1 = B-XY+Dy(YY,1Dept. of Mathematics, University of Montana135Appendix B. An Exact Result for Coupled Brusselators 136= A—(B+1)X2+XYDx(Xi-),= BX2— XY2 + Dy(Y1 — Y). (11.2)The trivial solution to 11.2 is the spatially uniform one:X1=X2A,Y1=Y2B/A. (11.3)Changing variables to departnres from this solution:x1=X—A,x2—XA, yi=Yi—B/A,y2=Y—B/A. (11.4)The equations then conveniently split into parts linear and nonlinear about the homogeneous steady state (xi = O,x2= O,yi = O,Y2 = 0):(B—l)xi+ A2y1 + 2Ax1y+ + xy1 + Dx(x2 — xi).—Bx1 — A2y1— 2Ax1y— — xyi + Dy(y2—= (B — 1)x2 + A2y + 2Axy+ + xy2 + Dx(xi —= —Bx2 — A2y — 2Axy—— xy2 + Dy(y1— y2). (11.5)The system is symmetric under the interchange of indices 1 and 2 so any solution musthave this symmetry i.e.:xi = X2,y1= 112, (11.6)(as was the case for the homogeneous steady state), or must be accompanied by its pairwhich may be generated by exchanging the indices 1 and 2 in the expressions for thesteady state concentrations.The brusselator has the special property that the departure of the variable X from thehomogeneous state must be nul when integrated over a system with no-flux conditions.This property holds for the present example, and allows elimination of one of the nonlinearequations from 11.5.Appendix B. An Exact Result for Coupled Brusselators 137Since:thi+i=—xi+Dx(x2—xi)+Dy(yyi),X2 + Y2 = + Dx(z1 — x2) + Dy(y1— y2), (B.7)then, at a steady state:0 = (x1 + x2), or x1 = —x2. (B.8)Eliminating the equation for ‘2:o (B — 1)xi + A2y1 + 2Axiyi xxyi — 2Dx1,o = —Bx1 — A2y1 — 2Axiyi —— xy1 + Dy(y2—o = —(B — 1)xi + A2y — 2Axiy + xxy2+ 2Dx1. (B.9)The first two equations may be added and solved for x1:I_______= 1 + 2Dx) (y2 yr). (B.1O)The first and third equations are linear in Yi and y2 respectively, and each contain onlytwo of the three variables. Then these can be re-written:(A2 + 2Ax1 + x)y1 = —((B — 1)xi + — 2Dxx1),(A2— 2Ax1 + x)y2 = ((B — 1)xi —— 2Dxx1). (B.11)The (pathological) solutions x1 = +A implies a special relation between parameters.Barring this, closed form expressions for y and Y2 are:— —xi((B— 1) H- x1 — 2Dx)—(x1 +A)2— xi((B — 1)— x1— 2Dx)B 12Y2— (xi—A)2 . ( . )Appendix B. An Exact Result for Coupled Brusselators 138Then, at the steady state, x1 satisfies:—4 Dy (((B—1)—xi—2Dx) ((B—1)+xi—2Dx) B13X1 1+2Dx)x1 (xi—A)2 + (xi+A)2From the expressions for yi and Y2 that x1 = 0 implies Yl = Y2 = 0 so this solutionis merely the homogeneous steady state. For non-zero x1, dividing through by x1 andsimplifying the result, the equation to be solved is:(Xi + A)2(xi - A)2 = -2Dy(x + A2)+-2BDx1 (B.14)The equation for the steady state values of x1 takes the form:(B.15)with:F = 2Dv_2A2_(f)Q = A2(A2_2Dy(B1x)). (B.16)This is a quadratic in x with solution:X=(A2+iL (B.17)There will be four real roots of the original system if there are two positive roots theabove equation. i.e. x > 0 for both the + and — branch. There will be only 2 roots ifone branch is negative, and no roots if both are negative. The analytical calculation ofthe local stability of these solutions appears to be complicated.Bibliography[1] M. Akam, Development, 101: 1-22, 1987.[2] M. Akam, Nature, 341: 282-283, 1989.[3] P. A. Anderson, J. Olavarria, and R.. C. \7an Sluyters, J. Neuroscience 8: 2183, 1988.[4] \i. I. Arnold, Ordinary Differential Equations, MIT Press, Cambridge MA, 1973.[5] A. Babloyantz, Molecules, Dynamics, and Life- An Introduction to the SelfOrganization of Matter John Wiley and Sons, New York, 1986.[6] J. Bard and I. Lauder, J. Theor. Biol., 45: 501-531, 1974.[7] J. B. L. Bard, J. Theor. Biol., 93: 363-385, 1981.[8] W. E. Boycer and R. C. DiPrima Elementary Differential Equations and BoundaryValue Problems. 3rd ed., John Wiley and Sons, New York, 1977.[9] J. P. Boon and A. Noullez, in On Growth and Form - Fractal and Non-Fractal Patterns in Physics, H. E. Stanley and N. Ostrowsky eds., Martinus-Nijhout, Dordecht,Netherlands, 1986[10] M. Braun Differential Equations and Their Applications. 2nd ed., Springer-Verlag,New York, 1975.[11] N. F. Britten, Reaction-Diffusion Equations and their Applications to Biology, Academic Press, London, 1986.[12] B. Bunow, J. P. Kernevez, C. Joly, and D. Thomas, J. Theor. Biol., 84: 629-649,1980.[13] C. Byrne and E. C. Cox, Develop. BioL, 117: 442-455, 1986.[14] C. Byrne and E. C. Cox, Proc. Natl. Acad. Sci. USA, 84: 4140-4144, 1987.[15] C. Callaini, Development, 107: 35-41, 1989.[16] J. A. Campos-Ortega and V. Hartenstein, The Embryonic Development ofDrosophila melanogaster, Springer-Verlag, Berlin, 1985.[17] 5. B. Carroll, Cell, 60: 9-16, 1990.139Bibliography 140[18] V. Castets, E. Dubs, J. Boissonade, and P. De Kepper, Phys. Rev. Lett., 64: 2953-2956, 1990.[19] E. Conway, D. Hoff, and J. Smoller, Siam J. Appi. Math., 35: 1-16, 1978.[20] F. H. C. Crick and P. A. Lawrence, Science, 189: 340-347, 1975.[21] W. Driever and C. Nüsslein-Volhard, Cell, 54: 83-93, 1988.[22] W. Driever and C. Nflsslein-Volhard, Cell, 54: 9.5-104, 1988.[23] W. Driever and C. Nflsslein-Volhard, Nature, 337: 138-143, 1989.[24] P. Duchateau and D. W. Zachmann, Theory and Problems of Partial DifferentialEquations, Schaum’s Outline Series, McGraw-Hill, New York, 1986.[25] B. A. Edgar, G. M. Odell, and G. Schubiger, Developmental Genetics 10: , 1989.[26] A. Einstein, Investigations on the Theory of the Brownian Movement, Dover Publications, New York, 1956.[27] Th. Erneux andM. Herschkowitz-Kaufman, Biophys. Chem., 3: 345, 1975.[28] V. E. Foe and B. M. Alberts, J. Cell. Sci., 61: 31-70, 1983.[29] H. G. Frohnhöfer and C. Nüsslein-Volhard, Nature, 324: 120-125, 1986.[30] A. Gierer and H. Meinhardt, Kybernetik, 12: 30-39, 1972.[31] A. Gierer, Jahrb. f. Nationalök. u. Stat., 196: 309-331, 1981.[32] 5. F. Gilbert, Developmental Biology, 2nd ed., Sinauer Associates, Sunderland MA,1988.[33] J. L. Gmitro and L. E. Scriven, Sympos. mt. Soc. Cell Biol. 5: 221-255, 1966.[34] B. C. Goodwin and S. A. Kauffman, J. Theor. Biol., 144: 303-319, 1990.[35] H. Haken, Z. Physik, B21: 105-114, 1975.[36] H. Haken and H. Olbrich, J. Math. Biol., 6: 317-331, 1978.[37] H. Haken, Synergetics, Nonequilibrium Phase Transitions and Self-Organization inPhysics, Chemistry, and Biology, 3rd ed., Springer-Verlag, Berlin, 1983.[38] H. Haken, The Science of Structure: Synergetics, Van Nostrand Reinhold, New York,1984.Bibliography 141[39] H. Haken, Information and Self-Organization- A macroscopic Approach to ComplexSystems. Springer-Verlag, Berlin, 1988.[40] L. G. Harrison and T. C. Lacalli, Proc. R. Soc. Lond. B, 202: 361-397, 1978.[41) L. G. Harrison, in Origins of Optical Activity in Nature, D. Walker ed., ElsevierScientific, Amsterdam, 1979.[42] L. G. Harrison, J. Snell, It. Verdi, D. E. Vogt, G. D. Zeiss, and B. It. Green, Protoplasma, 106: 211-221, 1981.[43] L. G. Harrison and N. Hillier, J. Theor. Biol., 114: 177-192, 1985.[44] L. G. Harrison, J. Theor. Biol., 125: 493-515, 1987.[45] L. G. Harrison and M. Kolar, J. Theor. Biol., 130: 493-515, 1988.[46] L. G. Harrison, K. T. Graham, and B. C. Lakowski, Development, 104: 255-262,1988.[47] L. G. Harrison and B. It. Green, Can. J. Chem., 66: 839-851, 1988.[48] L. G. Harrison and M. J. Lyons, Proceedings of the 1991 Les Houches winter workshop: Dynamical Phenomena at Interfaces, Surfaces, and Membranes, in press.[49] L. G. Harrison, Kinetic Theory of Living Pattern Cambridge University Press (tobe pnblished), 1991.[50] D. Hearn and M. P. Baker, Computer Graphics, Prentice-Hall, Engelwood Cliffs NJ,1986.[51] M. Herschkowitz-Kanfman, Bull. Math. Biol., 37: 589-636, 1975.[52] M. Herschkowitz-Kaufman and G. Nicolis, J. Chem. Phys., 56: 1890-1895, 1972.[53] A. Hodges, Alan Turing- The Enigma of Intelligence, Unwin Paperbacks, London,1985.[54] F. C. Hoppensteadt, An Introduction to the Mathematics of Neurons. CambridgeUniversity Press, Cambridge, 1986.[55] It. V/. Hornbeck, Numerical Methods, Quantum Publishers, New York, 1975.[56] K. Howard, Nature, 338: 618-619, 1989.[57] D. H. Hubel and T. N. Wiesel, Proc. It. Soc. Lond. B, 198: 1-59, 1977.Bibliography 142[58] D. H. Hubel and T. N. Wiesel, Scientific American, 241: 150-162, 1979.[59] M. Hülskamp, C. Schröder, C. Pfiefie, H. Jackie, D. Tautz, Nature, 338: 629-638,1989.[60] A. Hunding, in Cell to Cell Signalling: From Experiments to Theoretical Models,A. Goidheter ed., Academic Press, 1989.[61] A. Hunding, S. A. Kauffman, and B. C. Goodwin, J. Theor. Biol., 145: 369-384,1990.[62] A. Hunding and A. K. Saha, pre-print[63] A. Hunding, preprint[64] I. Huntley and R. M. Johnson, Linear and Nonlinear Differential Equations, EllisHorwood Ltd.. Chichester UK, 1983.[65] P. W. Ingham, Nature, 335: 25-34, 1988.[66] V. Irish, Nature, 338: 646-648, 1989.[67] 5. Jakubith, H. H. Rotermund, W. Engel, A. von Oertzen, and G. Ertl, Phys. Rev.Lett. 65: 3013-3016, 1990.[68] T. L. Karr and T. R. Kornberg, Development, 105: 95-103, 1989.[69] 5. A. Kauffman, R. M. Shymko, and K. Trabert, Science, 199: 259, 1978.[70] Johannes Kepler, The Six Cornered Snowflake, first published 1611, english translation, Clarendon Press, Oxford, 1966.[71] L. Edelstein-Keshet, Mathematical Models in Biology, Random House, New York,1988.[72] 5. E. Koonin, Computational Physics, Benjamin Cummings, Menlo Park, 1986.[73] Y. Kuramoto, Prog. Theor. Phys. Suppl., 64: 346.[74] K. Kurata, K. Kismoto, and E. Yanagida, J. Math. Biol., 27: 485-490, 1989.[75] T. C. Lacalli anf L. G. Harrison, J. Theor. Biol. 70: 273-295, 1978.[76] T. C. Lacalli anf L. G. Harrison, J. Theor. Biol. 76: 419-436, 1979.[77] T. C.Lacalli, Phil. Trans. Roy. Soc. Lond. B, 294: 547-588, 1981.Bibliography 143[78] T. C. Lacalli, D. A. Wilkinson, and L. G. Harrison, Development, 102: 1-9, 1988.[79] T. C. Lacalli, J. Theor. Biol., 144: 171-194, 1990.[80] T. C. Lacalli and L. G. Harrison, Seminars in Developmental Biology, 2, 1991.[81] L. Lapidns and G. P. Pinder, NUmerical Solution of Partial Differential Equationsin Science and Engineering, Wiley-Interscience, New York, 1982.[82] I. Lengyel and I. R. Epstein, Science, 251: 650-652, 1991.[83] S. LeVay, M. P. Stryker, and C. J. Shatz, J. Comp. Nenr. 179: 223-244, 1978.[84] S. LeVay, T. N. Wiesel, and D. H. Hubel, J. Comp. Neurol. 191: 1, 1980.[85] C. C. Lin and L. A. Segel, Mathematics Applied to Deterministic Problems in theNatural Sciences, SIAM, Philadelpia, 1988.[86] D. Z. Loesch, Quantitative dermatoglyphics classification, genetics, and pathology,Oxford University Press, Oxford, 1983.[87] M; J. Lyons, A Deuterium NM]? Studyj of Gramicidin A, M. Sc. thesis, Universityof British Columbia, 1985.[88] M. J. Lyons, L. G. Harrison, B. C. Lakowski, and T. C. Lacalli, Can. J. Phys., 68:772-777, 1990.[89] M. J. Lyons and L. G. Harrison, Chem. Phys. Lett. 183: 158-164, 1991.[90] M. J. Lyons and L. G. Harrison, Analysis and Modelling of Neural Systems II, KiuwerAcademic, in press, 1991.[91] M. J. Lyons and L. G. Harrison, Pattern Formation by Reaction-Diffusion, University of British Columbia, 1990, SVHS video produced with the assistance of the UBCUniversity Computing Services.[92] L. J. Manseau and T. Schflpbach, Trends in Genetics, 5: 400-405, 1989.[93] J. Maynard Smith, Mathematical Ideas in Biology, Cambridge University Press,Cambridge, 1968.[94] H. Meinhardt, J. Cell Sci., 23: 117-139, 1977.[95] H. Meinhardt Models of Biological Pattern Formation Academic Press, New York,1982.[96] H. Meinhardt, J. Cell Sci. Snppl., 4: 357-381, 1986.Bibliography 144[97] H. Meinhardt, Development Suppl., 104: 95-110, 1988.[98] K.D. Miller, J.B. Keller, and Michael P. Stryker, Science 245: 605, 1990.[99] K. D. Miller, in Neuroscience and Gonnectionist Theory, M. A. Gluck andD. E. Rumelhart (eds.), Lawrence Erlbanm Associates, Hillsdale NJ, 1990.[100] K. D. Miller and M. P. Stryker, in Connectionist Modeling and Brain Function: TheDeveloping Interface. S. J. Hanson and C. R. Olson (eds.), MIT Press, CambridgeMA, 1990.[101] J. D. Mnrray, J. Theor. Biol., 88: 161-199, 1981.[102] J. D. Mnrray, J. Theor. Biol., 98: 143-163, 1982.[103] J. D. Mnrray, P. K. Maini, and R. T. Tranquillo, Physics Reports, 171: 60-84 1988.[104] J. D. Murray, Scientific American, 258: 80-97, 1988.[105] J. D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1989.[106] B. N. Nagorcka, J. Theor. Biol., 132:Z77-306, 1988.[107] B. N. Nagorcka, J. Theor. Biol., 137: 127-162, 1989.[108] G. Nicolis and I. Prigogine Self Organization in Nonequilibrium Systems, FromDissipative Structures to Order Through Fluctuations, John Wiley and Sons, NewYork, 1977.[109] C. Nflsslein-Volhard and E. Wieschaus, Nature, 287: 795, 1980.[110] C. Nüsslein-Volhard, C. Frohnhöfer, and R. Lehmann, Science, 238: 1675-1681,1987.[111] R. Pool, Science, 251: 627, 1991.[112] G. M. Odell, G. Oster, P. Alberch, and B. Burnside, Develop. Biol., 85: 446-462,1981.[113] G. F. Oster and G. M. Odell, in Mathematical Problems in the Life Sciences, 13,Amer. Math. Soc., Providence, RI, 1980.[114] G. F. Oster, Math. Biosci., 90: 265-286, 1988.[115] G. F. Oster and J. D. Murray, J. Exp. Zool., 251: 186-202, 1989.[116] H. E. Othmer and L. E. Scriven, I. and E. C. Fundamentals, 8: 302-313.Bibliography 145[117] Q. Ouyang and H. L. Swinney, Nature, 352: 610-612, 1991.[118] R. Penrose, Ann. Hum. Genet. Lond., 42: 435-445, 1979.[119] W. H. Press, B. P. Flannery, S. A. Teukoisky, and W. T. Vetterling, NumericalRecipes, The Art Of Scientific Computing, Cambridge University Press, Cambridge,1986.[120] I. Prigogine and G. Nicolis, J. Chem. Phys., 46: 3542-3550, 1967.[121] I. Prigogine and R. Lefever, J. Chem. Phys., 48: 1695-1700, 1968.[122] I. Prigogine, From Being to Becoming, Time and Complexity in the Physical Sciences, W. H. Freeman, New York, 1980.[123] I. Prigogine and I. Stengers, Order Out of Chaos, Man’s New Dialogue with Nature,Bantam Books, Toronto, 1984.[124] J. J. Ramos, Mathematics and Computers in Simulation, 25: 538-548, 1983.[125] E. E. Selkov, Eur. J. Biochem., 4: 79-86, 1968.[126] B. Shorrocks, Drosophila, Ginn and Company, 1972.[127] J. M. W. Slack, From Egg to Embryo, Determinative Events in Early Development,Cambridge Univeristy Press, Cambridge, 1983.[128] B. P. Sonnenblick, in Biology of Drosophila, M. Demerec ed., Hafner PublishingCo., New York, 1950.[129] N. V. Swindale, Proc. R. Soc. Lond. B, 208: 243-264, 1980.[130] G. Strang Introduction to Applied Mathematics, Wellesley-Cambridge Press,Wellesley MA, 1986.[131] J. Tabony and D. Job, Nature, 346: 448-351, 1990.[132] W. Y. Tam, J. A. Yastano, H. L. Swinney, and W. Horsthemke, Phys. Rev. Lett.,61: 2163-2166, 1988.[133] 5. Tanaka, Nec Research and Development, 98: 1-14, 1990.[134] D’Arcy Thompson, On Growth and Form, Abridged edition edited by J. T. Bonner,Cambridge University Press, Cambridge, 1961.[135] J. R. Thompson, Z. Zhang, Wm. Cowan, M. Grant, J. A. Hertz, and M. J. Zuckermann, Physica Scripta, T33: 102, 1990.Bibliography 146[136] A. M. Turing, Phil. Trans. RoyL Soc. Lond. B, 237: 37-72, 19.52.[137] J. J. Tyson and J. C. Light, J. Chem. Phys., 59: 4164-4173, 1973.[138] J. J. Tyson, J. Chem. Phys., 62: 1010-101.5, 1975.[139] J. J. Tyson and J. D. Murray, Development, 106: 421-425, 1989.[140] C. H. Waddington, Principles of Embryology, George Allen and Tlnwin, London,1956.[141] D. Walgraef, G. Dewel, and P. Borckmans, in Advances in Chemical Physics Vol.49, I. Prigogine and S. A. Rice eds., Wiley Interscience, NY, 1982.[142] P. R. Wallace Mathematical Analysis of Physical Problems, Dover Publications,New York, 1972[143] M. P. Weir and T. Kornberg, Nature, 318: 433. 1985.[144] A. T. Winfree, The geometry of biological time, Springer Verlag, New York, 1980.[145] A. T. Winfree, The timing of biological clocks, Scientific American Library,Xis!. H. Freeman, New York, 1987.[146] A. I. Winfree, Nature. 352: 568-569, 1991.[147] D. A. Young, Math. Biosci., 72: 51-58, 1984.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items