UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Reaction-diffusion modelling of somite formation : computed dynamics and bifurcation analysis Orchard, Jeff 1996

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

Item Metadata


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

Full Text

REACTION-DIFFUSION MODELLING OF SOMITE FORMATION: C O M P U T E D DYNAMICS AND BIFURCATION ANALYSIS. By Jeff Orchard . B . Math. (Applied Mathematics) University of Waterloo , 1994 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES D E P A R T M E N T O F M A T H E M A T I C S I N S T I T U T E O F A P P L I E D M A T H E M A T I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A August 1996 © Jeff Orchard, 1996 . In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Hcy-Us e yv\c\4~ ( C S The University of British Columbia Vancouver, Canada Date SeplewAir DE-6 (2/88) Abstract Living organisms exhibit pattern in many forms. The processes underlying this organization, in general, are not well understood. This thesis proposes a candidate model for the dynamic eruption of somites (the precursors for the structures associated with the spinal column) in vertebrate embryos. We apply the Brusselator reaction-diffusion chemical model in a dynamic setting to not only produce the desired pattern, but also to mimic the propagation of the region of somite production (the "cohesive zone") along the embryo. The model can be rendered into a system of partial differential equations, and studied analytically and numerically. A bifurcation analysis is performed on a certain trajectory in the system's parameter space. Further numerical analysis on a more biologically robust version of the problem yields evidence for more than one stable, steady-state solution for certain parameter sets. Movement of the region of somite formation is simulated by adding a feedback mechanism between the peak value of the X morphogen, and a gradient in the source chemical B. Finally, investigation of the effect of adding a diffusion barrier yields some results that may be tested experimentally, creating evidence either for or against the possibility that a Brusselator-like chemical system is responsible for somite production in vertebrate embryos. ii Table of Contents Abstract ii List of Tables v List of Figures vi Acknowledgement vii 1 Introduction 1 1.1 Biological Pattern Formation 1 1.2 Somite Formation 2 1.3 Reaction-Diffusion and the Brusselator 4 2 Computer Simulations of the Brusselator 7 2.1 Description of the Simulation Program 7 2.2 Choosing a Pattern Instigator 9 2.3 Spatial Gradient in bB 11 3 Mathematical Analysis of the Brusselator's Turing Bifurcation 15 3.1 Homogeneous Equilibrium 15 3.2 Linearized System 16 3.3 Mathematical Tools 16 3.3.1 Solution Representation 17 3.3.2 Properties of the Inner Product 18 3.3.3 Adjoint Operator 19 3.4 Parameter Values 20 in 3.5 Linearized Stability of the Homogeneous Equilibrium 22 3.6 Locating the First Bifurcation 24 3.7 Centre Manifold Reduction 24 4 Numerical Analysis of the Brusselator 34 4.1 Numerical Bifurcation Analysis Software Package: A U T 0 8 6 34 4.2 Checking A U T 0 8 6 Output with br. c 37 5 Non-Homogeneous B-Concentration 45 5.1 A U T 0 8 6 Solutions Involving a 6-Gradient 46 5.2 br. c Solutions Involving 6-Gradient 48 5.3 Biological Implications of 6-Gradient 49 6 Dynamic Production of Somites: Cohesive Zone Propagation 51 6.1 Cohesive Zone Boundaries 51 6.1.1 Known Mechanisms 51 6.1.2 Assumed Mechanisms 51 6.2 Changes in the Basic Model 52 6.2.1 Dynamic 6-Concentration 52 6.2.2 Source Placement in br. c 54 6.3 Results of Propagation Model 54 6.4 Mathematical Interpretation of Propagation Model 58 6.5 Diffusion Barriers 59 7 Conclusions 64 Bibliography 66 Appendix A : Algebraic Details 68 iv List of Tables 2.1 Parameter values for equation 1.2 8 3.1 Parameter values for equation 3.3 21 6.2 Dynamic Cohesive Zone Simulation Parameters. A l l values are exact 55 v List of Figures 1.1 Human Embryological Development showing Somites 3 1.2 Scanning Electron Micrograph of an Axolotl Embryo 3 1.3 Zones of Somite Formation 4 2.1 Homogeneous Equilibrium 9 2.2 Non-Homogeneous Equilibrium 10 2.3 Non-Homogeneous Equilibrium induced by a A 11 2.4 Lacalli-Harrison Parameter Space 12 2.5 Equilibrium Solution for Linear 6B-Gradient 14 3.1 Solution Projection Schematic 20 3.2 P versus M when A M = 0 25 3.3 Pitchfork Bifurcation Schematic Diagram 33 4.1 Comparison of Bifurcation Plots: A U T 0 8 6 versus Analysis 38 4.2 br.c Solution Trajectory on Bifurcation Diagram 39 4.3 A U T 0 8 6 and br.c Solution Comparison 40 4.4 Bifurcation Diagram for Homogeneous b (stereogram) 41 4.5 Dispersion Plot for b = 0.435 42 4.6 Dispersion Plots: S and f2 Parameters 43 4.7 A U T 0 8 6 and br.c Solutions Comparison 44 5.1 Non-Symmetric Pitchfork Bifurcation Diagram 47 5.2 Bifurcation Diagram for Non-Homogeneous b (stereogram) 48 5.3 Non-Symmetric Pitchfork Bifurcation Diagram and br. c Solution Trajectory . . 49 vi 6.1 Propagation of Cohesive Zone 57 6.2 Diffusion Barrier with Cohesive Zone Propagation Recovery 61 6.3 Diffusion Barrier with Cohesive Zone Propagation Collapse 63 vii Acknowledgement The work reported in this document is a result of research conducted in conjunction with Dr. Wayne Nagata of the Department of Mathematics at the University of British Columbia, and Dr. Lionel G . Harrison of the Department of Chemistry at the University of British Columbia. Both played a significant role in the development of the ideas in this thesis, as well as aided me in the technical aspect. Conversations with Gregory Lewis contributed to the bifurcation analysis in chapter 3. Finally, Dr. David Holloway helped me understand the methods involved in the simulation program br. c. vii i Chapter 1 Introduction 1.1 Biological Pattern Formation Biological pattern formation includes many instances of small and large numbers of repeated parts: petals of a flower, hair follicles, segments of an insect body, vertebrae, etc. In some instances, all the parts form simultaneously (Acetabularia whorls). In other cases, these parts form sequentially (one-dimensional) or one row at a time with rows appearing in sequence (two-dimensional, chick feather follicles). This thesis is concerned with somite formation in vertebrate embryos, a biological patterning process that has both simultaneous and sequential aspects. There are no universally accepted theories to account for pattern formation in biological development. A fairly general dichotomy in views between physical scientists and biologists has been discussed by Harrison [9]. The concept that pattern is generated by nonlinear dynamics is widespread among mathematicians and physical scientists [8] [9] [15] [17] [18] [23]. In this thesis, one particular type of postulated dynamics, reaction-diffusion theory, is used to set up a model for the phenomenon of somite formation. It is intended that this model should eventually be tested in some detail against experiment for one particular species of salamander, Ambystoma mexicanurn, or the Mexican axolotl, as was done in an earlier project on heart formation [11] [12]. This study has been concerned entirely with a particular reaction-diffusion model, the Brus-selator of Prigogine and Lefever [20]. The strategy of the work has involved two essentially separate investigations with two supervisors: setting-up of the model for somite formation, 1 Chapter 1. Introduction 2 and computer experiments on its behaviour with Dr. L . G. Harrison (Chemistry, and Insti-tute of Applied Mathematics)(chapters 2 and 6); and bifurcation analysis of the Brusselator, with and without pre-existing gradients in an input to the mechanism, with Dr. W . Nagata (Mathematics, and Institute of Applied Mathematics) (chapters 3, 4 and 5). 1.2 Somite Formation We wish to investigate the process of the production of a row of repeated parts in the embryos of all vertebrates. The parts in question are called somites. A somite is an embryological structure composed of a group of cells that cluster and cleave away from other cells, forming a rounded box. Somites are the embryological precursors for structures such as the vertebrae, ribs, some nerve tissue, and many other structures associated with the spinal column. Somites are unusual instances of pattern since their formation is neither fully simultaneous nor fully sequential. At any time, about six somites are forming simultaneously in a "cohesive zone" which moves anteroposteriorly (head to tail) along the paraxial mesoderm to generate 30 to 40 somites. The paraxial mesoderm is a layer of tissue that runs the length of the embryo, alongside the neural tube. Figure 1.1 shows a series of human embryos in sequential stages of development. The somites are visible at each stage shown. Figure 1.2 is a scanning electron micrograph of an axolotl embryo in stage 29 of its development. Chapter 1. Introduction 3 n t Way X atest rtmsii' 23 $ 1 day fc 2* j, ldoy 2S i WNr Figure 1 1: Shown here is a selection of developmental stages of a human embryo (from stage 22 to stage 28). The somites erupt sequentially along the length of the embryo. Taken from Moore [16] Figure 1.2: This scanning electron micrograph of an axolotl embryo (in stage 29) clearly shows the segmentation of the somites (taken from Armstrong and Malacinski [2]). Chapter 1. Introduction 4 For the purposes of this paper, we can think of the paraxial mesoderm as divided into three zones (figure 1.3): a segmented zone, a cohesive zone, and a non-cohesive zone. The segmented zone contains fully-formed somites. The cohesive zone is made up of cells that have gap junctions through which diffusion may occur, forming a continuous diffusion medium for intracellular substances. The non-cohesive zone is filled with cells that are mobile. The cells of the non-cohesive zone do not have gap junctions, so diffusion of substances from one cell to another does not take place. D f \ f Somite V J ° ° ° o o ° o n o o o o o'-Segmented Zone Cohesive Zone o o O o ° 0 ° o ° o — Non-Cohesive Zone Figure 1.3: This schematic diagram of the paraxial mesoderm shows the three zones, identified by stages of somite production. The segmented zone has fully-matured somites. The cohesive zone has partially-formed somites and contains cells that have gap junctions. The non-cohesive zone contains mobile cells that do not have gap junctions. 1.3 Reaction-Diffusion and the Brusselator A number of instances of biological patterning have been investigated with reaction-diffusion models, as first devised by A . M . Turing [23]. At any point in a chemical medium, one can sample (conceptually, at least) the local concentrations of substances. Imagine that there are two chemicals (in a diffusion continuum of some sort) that react with one another. We will call these chemicals "morphogens" if they contribute in some way to pattern formation 1. Perhaps their reaction simply combines the two chemicals to make different ones. More interesting interplay involves somewhat convoluted chemical dynamics such as autocatalysis, or a third chemical adding a step to the reaction sequence. Furthermore, adding diffusion to the system creates the potential for balancing these chemical changes, replenishing the areas of low concentration from nearby areas of high concentration. 1See Harrison [9] (especially chapters 2 and 9) for a discussion of different meanings for "morphogen", and characterization of A and B (in reference to equation 1.1) as "type 1 morphogens" and X and Y as "type 2 morphogens". Chapter 1. Introduction 5 The idea behind reaction-diffusion models for pattern formation is one of prepat tern: A non-homogeneous distribution of chemicals may instigate other activities in a nonhomogeneous manner, and ultimately result in morphogenetic pattern. The Brusselator chemical system is a reaction-diffusion system. It tends to create spatially periodic patterns, and is described by the reaction equations A A X, B + X X Y + D, (1.1) Y + 2D A 3X , X A E , where X, Y, A, B, D and E are the molecules involved. The substances X, Y, A and B are the diffusible morphogens. The cohesive zone seems to be the region of somite formation in axolotl embryos, and hence we focus our attention and application of the Brusselator on that region. In particular, the time evolutions of morphogens X and Y are of interest. Each of the three zones (segmented, cohesive and non-cohesive) is long with respect to its width, and hence can be approximated by a one-dimensional spatial domain. We can therefore represent the spatial concentration profile of a morphogen as a scalar function. For example, we define the concentration of morphogen X as X(s, t) = f the concentration of morphogen X at time t and position s . We will assume that the morphogens X and Y are intracellular, and can diffuse throughout the cells of the cohesive zone. The ends of the cohesive zone form a diffusion barrier for substances X and Y, defining a boundary condition for their concentration. As a result, we can represent the dynamics of our chemical reactions in the cohesive zone with the system of Chapter 1. Introduction 6 partial differential equations, ^ = a A - bBX + cX2Y - dX + £ x 0 , % = bBX-cX*Y + Dv$£, (1.2) ax _ dY_ _ ax = 9yi — o 3s s = 0 as s = 0 as s = 1 as s = 1 A point source of a diffusible chemical is a biologically feasible explanation for the spatial localization of pattern. The nearby high concentration of the resulting spatial gradient produces pattern, while the low concentration in remote regions does not [10]. The concept of monotonic gradients is one that is widely accepted by biologists. We wish to investigate the possibility that the morphogens X and Y form the prepattern for somite formation. Perhaps a high concentration in the X morphogen induces the production and cleavage of a somite. A gradient in either chemical A or chemical B may localize the patterning region. These concepts will be explored in the following chapters. Mathematically, the introduction of spatial dependence in one or more of the parameters a, A, b, B, c and d of equation 1.2 greatly complicates the analysis of the dynamics (see chapter 5). For this reason, simulation will form a large part of our analysis. Chapter 2 Computer Simulations of the Brusselator Computer simulation is used to gain an intuition about the Brusselator's dynamics on a one-dimensional spatial domain. We will use a program called br.c, written in the C programming language by D. Holloway [11], to simulate our system. 2.1 Description of the Simulation Program The program br. c numerically solves the Brusselator system, as posed in equation 1.2, animat-ing the time evolution of the chemical concentrations. Its solution is based on the discretization of the spatial domain into N equally-sized subintervals. Suppose the concentrations of X and Y at spatial element i and at time step j are Xij and Yij respectively. The time derivative can be approximated by dXi X^ — Xij ~df ~ A t ' where A t is the length of the time step. We can also approximate the second spatial derivative term as d^Xij ^ Xi+ij — 2Xij + Xi—ij ds2 ~ ( A s ) 2 ' where A s is the length of the spatial element. Notice that the approximation to the spatial derivatives, when evaluated at one of the domain's endpoints, makes reference to the state of d 2 X - d2Y-• the system outside the domain. To calculate our approximations to dJ3 and gj1 at s = 0 and s = 1, we need to use the boundary conditions, dXpj Xpj—X-ij p. dXpi j XN+I J— XN j n ds ~ As — u > ds ~ As — u i 9 Y 0 j _ Ypj-Y-ij _ n 8 Y N j _ Y N + l j - Y N j _ „ ds ~ As ~ u ' ds ~ As ~ u • (2.1) Chapter 2. Computer Simulations of the Brusselator 8 The exterior points X-ij, X^+ij, Y-ij and Y^+\j can be determined by the boundary condi-tions, and hence eliminated from the system of algebraic equations. For each time step, there are 2(iV + 1) algebraic equations with 2(N + 1) unknowns: {Xij,Yij}^_Q. Variations of these approximations are used to simulate the system. A n implicit/explicit differencing method is used to model the diffusion part, employing the Crank-Nicolson method to solve the resulting tri-diagonal linear system [6]. A more complete description of the solution is given in Holloway [11]. A series of simulations was run in which the spatial domain [0,1] was divided into 100 subintervals. Starting from an initial state, br. c produces an animated display of the time evolution of the solution to the finite-difference approximation. For the parameter values given in table 2.1, figure 2.1 shows the final steady-state solution of the Brusselator on [0,1]. Parameter Group Value aA 0.05 bB 0.4 c 1.8 d 0.05 Dx io- 4 JOy 10~ 3 Table 2.1: Parameter values for equation 1.2. However, if we increase bB to 0.46, we get a different kind of solution, as figure 2.2 depicts. This solution is wave-like, and qualitatively different from a homogeneous solution. Finally, the parameter choice of aA = 0.01 and bB = 0.4 gives us pattern as well (figure 2.3). The spike-like solution was preceded by a long transient (i.e. its convergence to equilibrium was slow). Chapter 2. Computer Simulations of the Brusselator 9 1.5 Cone. 0.5 0 0.2 0.4 0.6 0.8 1 s Figure 2.1: Starting from a random initial state, the ultimate state of the system, for the given parameter values and bB = 0.4, seems to be spatially homogeneous. 2.2 Choosing a Pattern Instigator The somites in an embryo are remarkably regular in their spacing. That is to say, somites appear in regular intervals along the length of the neural tube. Thus, a good candidate model needs to produce concentration peaks that have a regular spatial placement. It also must create the patterns in a timely manner, as embryological development has an agenda. We have seen that certain combinations of levels of the source chemicals A and B can admit evenly-spaced peaks in the resulting equilibrium solutions. Which parameter group, when varied, best produces evenly-spaced peaks in a time-efficient manner? After examining the time evolution of a number of initial states and parameter values, it seems that low values of both aA and bB consistently exhibit longer transients before reaching a steady-state than do cases with high values of both a A and bB. Some light may be shed on this phenomenon by looking at the Lacalli-Harrison parameter space. Linear stability analysis was first done on such two-morphogen models by A. M . Turing [23]. Lacalli reduced the stability of the homogeneous equilibrium solution to a two-parameter space: k[ versus k'4 [13]. For the Chapter 2. Computer Simulations of the Brusselator 10 1.5 X 1 h Cone. 0.5 0 0 0.2 0.4 0.6 0.8 1 s Figure 2.2: Starting from a state near the homogeneous equilibrium, the ultimate state of the system, for the given parameter values and bB = 0.46, seems to be a spatial wave. Brusselator, the parameters k[ and k'4 are defined as Figure 2.4 shows the Lacalli-Harrison parameter space, and the effect of decreasing aA. Ac-cording to Turing's linear analysis, the shaded regions labeled Q, and S contain values of k\ and k'4 which result in instability of the homogeneous equilibrium solution, and hence suggest convergence to a spatial pattern of some sort (since chemical concentrations cannot increase indefinitely). Figure 2.3 shows an equilibrium solution corresponding to parameter values in region £1, while figure 2.2 shows an equilibrium solution corresponding to parameter values in region E . Recall that a point source of a substance produces high concentrations of its chem-ical in the neighbourhood of its source. For our purposes, this suggests that a decrease in the concentration of our pattern-instigating chemical should correspond to a parameter-space tra-jectory that takes us from a region with pattern to a region of patternlessness. Moreover, the d[bB-d aAy/bBc —aA^/c Chapter 2. Computer Simulations of the Brusselator 11 1 0 1 1 i i x — A A y • r i i i i 0 0.2 0.4 0.6 0.8 1 s Figure 2.3: Starting from a state near the homogeneous equilibrium, the ultimate state of the system, for the given parameter values, and aA = 0.01 and bB = 0.4, depicts pattern, but more spike-like than wave-like. region of pattern needs to produce this pattern with a relatively short transient. A decrease in a A can take us from region E , through region f i , to a region in which no particular pattern dominates over another. However, decreasing bB creates a trajectory similar to the one shown in figure 2.4, but in the opposite direction. The trajectory can start in region E (in which transients are short), cross the "Turing Bifurcation", and end in a region of patternlessness. The parameter group bB was chosen to be manipulated as a pattern instigator. In mathe-matical terms, (some scalar multiple of) bB will act as the bifurcation parameter. 2.3 Spatial Gradient in bB The cohesive zone is an area in which somites exhibit varying degrees of maturation. At one end of the cohesive zone, somites are fully mature. At the other end of the cohesive zone, somites are not even visible. If we suppose that, at any position s along the cohesive zone, the concentration of the X morphogen defines the degree to which a somite at that location has matured, then the X solution we want has a wave pattern on one end of the cohesive zone and Chapter 2. Computer Simulations of the Brusselator 12 Figure 2.4: (Lacalli-Harrison parameter space) The shaded regions contain parameter values that, upon linear stability analysis of the homogeneous steady-state solution, admit the growth of patterned solutions, and seem (from computer simulation) to produce similarly-patterned steady-state solutions. Decreasing aA moves the parameter set along the indicated curve in the direction of the arrow, and hence away from the Turing Bifurcation. The slope of the line indicating the Turing Bifurcation is the ratio 5*-. no pattern on the other. On one side of the space, we need our parameter bB to be inside the shaded, pattern-forming region of the Lacalli-Harrison space, and on the other end, we need bB to be on the patternless side of the "Turing Bifurcation" 1. We achieve this trajectory in a biologically realistic way by creating a gradient in the source chemical B, forming a high concentration at one end, and a low concentration at the other. Figure 2.5 shows a solution of the system, with its corresponding B gradient decreasing linearly from left to right. We will return to the analysis of the solutions with a 5-gradient in chapter 6. But first, we 1It should be recognized that the addition of a gradient, in principle, invalidates the linear stability analysis corresponding to the Lacalli-Harrison parameter space, but that this analysis probably remains a good guide in an approximate sense, especially if the gradient is small. Chapter 2. Computer Simulations of the Brusselator want to know exactly what type of bifurcation the "Turing Bifurcation" Chapter 2. Computer Simulations of the Brusselator 14 X B Y 1 h Cone. 0.5 0.2 0.4 0.6 0.8 Figure 2.5: The 5-concentration decreases linearly from left to right. The higher value of B on the left induces pattern, while the lower concentration on the right prevents it. The formation of several peaks simultaneously, but reaching high amplitudes sequentially, in a region with a gradient in a dynamic parameter, was illustrated by Lacalli, Wilkinson and Harrison [14]. Chapter 3 Mathematical Analysis of the Brusselator's Turing Bifurcation Nicolis and Prigogine analysed the Turing bifurcation of the Brusselator in their book [18] in which they calculated the non-homogeneous solutions near the Turing bifurcation as a sum of sinusoids. The following bifurcation analysis has a similar goal to that of Nicolis and Pr i -gogine, but is presented in the context of current dynamical systems practices, containing more conceptual and graphical methods. 3.1 Homogeneous Equilibrium The partial differential equations corresponding to the Brusselator system are shown in equation 1.2. Letting a = a A and b = bB, our system is d2x d 2 Y <^r = a-bX + cX2Y -dX + Dx dx dr dY dr dX w = bX - cX2Y + Dy^, (3.1) dY s=Q \s=Q dX ds _ dY s=l s=l 0 Homogeneous steady-state solutions are found by setting Q£ = ^ = 0 and = = 0. Solving the resulting algebraic system yields the homogeneous equilibrium, Y0 = db d ac With the change of variables u = X—XQ and v = Y — YQ that translates the system's equilibrium (X = Xo and Y = YQ) to the origin (u = 0 and v = 0), the Brusselator system becomes du _ /L J \ i ca2 „, , db„,2 , o«c„„, , „„,2„, , n d2u ^ = (b - d)u + <%-v + fu2 + 2fuv + cu'v + Dx 9 g 2 , v - - 2fuv - cu2v + D y ^ , = 0 . dv dr -bu ca (3-2) du ds dv du dv s=0 ~~ 9 s s=0 ds 5 = 1 ds s=l 15 Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 16 After reseating the time variable, t = dr, the system of equations for the Brusselator becomes (3.3) ^ = (B - l)u + ceo + (5KV? + juv + av?v + D x ^ , dt dv d2i = — du — av — (3KU — 'juv — au v + Dy d s = 0 , du ds dv du dv s=0 ~~ 9 s s=0 ~ d s s=l ds s=l d • where 0 = | , a = « = 7 = 2 f , a = % Dx = % and Dy 3.2 Linearized System We wish to learn more about the equilibrium solutions and their stability. To do this, we linearize the system of partial differential equations 3.3 about its point of equilibrium u = 0, v = 0. Taking only the terms that operate linearly on u and v, we are left with fu = L{(5)u (3.4) where u = and L(B) = -8 a -a + D. d2 yds 3.3 Mathematical Tools Before we go on, we need to define some mathematical constructs that will help us later. Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 17 3.3.1 Solution Representation We seek solutions, u(t, s), of system 3.3, in the space1 X[0,1] . This solution space is a Hilbert space and is defined as, X[0,1] = {u(s) = du du\ / d2u d2u\ du, , du, . : < 0 0 ' { d - s ^ ) < C O > \ d ^ > d ? ) < C O > dS{0) = (3.5) where (it, v) is the inner product denned by (u,v) = 2 / u(s) • v(s) ds . Jo We can use the set of vectors {|.|}u{[; cos nns, cosmrs n=l as a basis for X[0,1], and thus also express our solution space in terms of the coefficients of Fourier cosine series: X[0,1] = { u(s) = + Y^=l u ^ cos nirs 4 2 ) , ^00 (2) n=l 4 L(D n \u. < oo, J2 n 4 u^ n=l < OO (3.6) The boundary conditions are automatically satisfied by the new solution representation. For some purposes, this space can be interpreted as an infinite-dimensional Euclidean space of vectors of the form (uW u(2) u W u{2) u{l) u{2) V The linear differential operator L(0) distributes over the Fourier modes of the solutions in X[0,1] so that we can write L((3)u as (see appendix for details) oo L(P)u = LQ(P)u0 + Y, Ln(P)Un , (3.7) n=l 'Notice that the symbol for the solution space X[0,1] is different than the symbol for the morphogen X (the latter being italicized). Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 18 where we define LQ(B) and Ln{0) as Lo(B) = Ln{0) ^ 0-1 a -0 -a 0 - 1 - n2n2Dx -0 a -a — n2-w2Du (3.8) (3.9) 3.3.2 Properties of the Inner Product One key property of our representation of X[0,1] is that its basis elements are orthogonal with respect to the inner product. For example, assuming m =fi 0, define cos m-ns 0 and consider the inner product "o1' , T-oo (1) 2 + 2_m=l u™ C 0 S n 7 r S «n 2 ) , v-oo (2) 2 + 2^n=l un C 0 S n i r s VL ir n=l 1 0 (1) Un cos nirs (2) un cos nirs (u,eim) = 2 / JO cos ?i7rs cos mirsds 2 2 2 (2) eim. ds \ ) C O S m 7 T S + ^ 71=1 C O S m 7 T S 0 1 0 _ uW _ cos nirs cos m7rs ds "'m • The inner product of an arbitrary vector with this m t h Fourier basis vector returns the amplitude of the m t h Fourier mode. This fact will be useful when we want to look at the dynamics of a bifurcating mode along a particular vector space. Defining our solution space as we do in equation 3.6 decouples the Fourier modes of the Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 19 linearized problem, and our linear system becomes, r i at 0-1 a ' uW ' a (2) L dt J -p —a (2) r d u % i at (0-1) — m 2 7 r 2DX a (2) at -B a I u -a — m2ir2Dy ^ u-m = 1, 2, 3 , . . . an infinite-dimensional system of linear ordinary differential equations. (1) m (2) (3.10) 3.3.3 Adjoint Operator Collapsing a multi-dimensional solution down into one or two dimensions facilitates the obser-vation of specific dynamics of the system. For this reason, it is useful to project solutions in a particular direction. To do this, we need the adjoint operator. The adjoint operator of L is defined as the operator L* such that {Lu, v) = (u, L*v) . For L = L(P), we have (see appendix for details) (L(0)u,v) = (u,L*(0)v) , (3.11) where a L*(0) = provided v satisfies the boundary conditions (P-1) + Ds a2 -0 -ot + D. a2 yap ds s=0 ds s=0 ds s=l ds = 0 . 5 = 1 (3-12) The usefulness of the adjoint operator comes from the Fredholm alternative [22]. Theorem 1 (Fredholm alternative): Let L be a differential operator with a compact resolvent on a separable Hilbert space H. Then, either of the following holds: 1. Lu = 0 has only the zero solution, in which case L*v = 0 has only the zero solution and Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 20 Figure 3.1: The projection of an arbitrary vector u onto NL*, the null space of the adjoint operator. The result is a scalar value q that is proportional to (u,x). Lit = f has precisely one solution for each f in H. 2. The null space, Ni, of L has finite dimension k, in which case NL* also has dimension k and Lu — f has solutions if and only if f is orthogonal to NL* . If one eigenvalue of L((3) is zero, then we have a non-trivial (finite-dimensional) null space: the centre eigenspace. Thus, our application of the Fredholm alternative would be through the second option listed above. In this case, the Fredholm alternative would tell us that the null space of the adjoint operator is orthogonal to the range of L(P). Thus, projecting a vector onto the centre eigenspace (or null space) of L*(P) essentially annihilates all the components that are in the range of L(P), leaving only the component parallel to the centre eigenspace of L*(P). Figure 3.1 gives a schematic of the projection. In the diagram, the projection of an arbitrary vector u G X[0,1] is found. The result is the scalar value q. 3.4 Parameter Values The Brusselator system has a number of parameters. In our bifurcation analysis, we will allow b, and hence P, to change. We must decide on some values for the other parameters. Table 3.1 gives the values of the seven parameters used in equation 3.3. They correspond to the six Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 21 parameter values given table 2.1. The values were taken from D. Holloway [11]. They fall in region S of the Lacalli-Harrison parameter space show in figure 2.4. Parameter Value a 36 P « 8 7 72 K 1 a 36 Dx 0.002 Dy 0.02 Table 3.1: Parameter values for equation 3.3. The values in tables 2.1 and 3.1 are stated without units because we are interested in the capabilities of the model, and not, at this point, in the comparison of our model directly to experimental data. However, we can do a thumbnail calculation to check that some of our parameters are at least of an appropriate order of magnitude. The formula for the wavelength, L , that corresponds to the local maximum of the dispersion curve (discussed in chapter 4) is (see Harrison [9], page 233), L = 2vr -a2c -b + d Dy + Dx I ba2c ~~dF I Dy ~DX Dy-Dx\l dbyDX (3.13) If we suppose that Dy is ten times Dx, then we can substitute Dy = 10DX into equation 3.13. From Armstrong and Graveson [1], we can get the approximate values of 6t = 9290 s for the amount of time it takes for a somite to form, and L = 150 pm for the width of a somite. We can also get measurements for these values from figure 6.1 in chapter 6. Wi th respect to the time unit ut and the length unit u^, the values are 8t = 20 Ut and L = 0.167 u^. Substituting L = 150 pm, as well as the remaining parameter values from table 2.1, except for b = 0.61 (taken from chapter 6), into equation 3.13, we can solve for Dx. Dx = 549.945 ^ m V - 1 Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 22 We can convert from ut time units to seconds by the equation, 20 u t = 9290 s 1 • u t = 464.5 s . Then our diffusivity, in cm 2 per second, is, Dx = 0.6660 x 1 0 - 8 c m V 1 . (3.14) This value is of an appropriate magnitude. Moreover, we can use our unit conversions to convert our value for the parameter Dx to the same units. bx = io" 4 u A r 1 = 0.1744 x 10~ 8 c m 2 s _ 1 The order of magnitude matches the estimate in equation 3.14. 3.5 Linearized Stability of the Homogeneous Equilibrium Suppose we want to find the stability of the equilibrium at the origin. The origin is u(s) = (0,0) T , which is equivalent to ( V ! Q \ U Q \ . . . j = (0,0, 0,0,. . .) in the Euclidean space. To find the stability of this equilibrium, we must find the eigenvalues of the linear operator at that point. / ! = det det (x0I-L(P)\MT A 0 0 0 A 0 6-1 a -6 -a = X20 + X0(-P + l + a) + a = 0 . Thus, Ao =\ (fi - 1 - a ± yJ(-6 + 1 + a)2 - 2a) . Using the parameter values listed in table 3.1 (and setting 6 = 8), our eigenvalues are ap-proximately -1.2996 and -27.7004. Thus, for the given parameter values, the homogeneous steady-state is stable with respect to homogeneous perturbations. Examining further the linear stability of the homogeneous equilibrium, we consider linear Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 23 combinations of (non-homogeneous) perturbations of the form u = eu$ cos MTTS cos MTTS To find the linear stability of the homogeneous equilibrium with respect to a perturbation in any particular mode M, we calculate the corresponding eigenvalues, XM+ and A M - -det (\mI - L(P)M) = det A M 0 0 A M _ B-l- M 2 T T 2 £ > ~P a -a - M 2 v r 2 A , A M + A M (-0 + 1 + a + M2ix2(Dx + £>„)) + ( - / ? + 1 + Af 27r 2£> 1) ( a + M V A , ) + a/? = 0 To conclude that the homogeneous equilibrium is asymptotically stable for the above pa-rameter values, we need to show that the real part of each eigenvalue for the M t h mode is negative and bounded away from zero. In other words, we need to show that, for the parameter values given, the only values of A M that are solutions for A M + A M (-/? + 1 + a + M2ir2(Dx + Dyj) + ( - / ? + 1 + M2TT2DX) (a + M2ir2Dy) + a/3 = 0 are negative and bounded away from zero. Equivalently, we can show that the constant term, + 1 + M 2 T T 2 D X ) (a + M2ir2Dy) + aP , (3 .15) is positive and bounded away from zero. Substituting the parameter values into expression 3.15, we are left with the expression M V - — M V + 36 . (3 .16) 25000 250 V ; For M = 0, expression 3.16 evaluates to 36. However, if we set expression 3 .16 equal to zero and solve for M, we arrive at four complex solutions, all of which are bounded away from the real axis. This fact2 tells us that there is no real value for M that evaluates expression 3 .16 arbitrarily close to zero. Hence, we can conclude that (for the above parameter values) the 2In conjunction with the mathematical theorems associated with persistence of sign and continuity Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 24 homogeneous equilibrium solution is asymptotically stable with respect to any perturbation in X[0,1]. 3.6 Locating the First Bifurcation A local bifurcation occurs when the stability of an equilibrium or periodic solution changes. Since the sign of the real part of the eigenvalues gives us the stability information for the homogeneous equilibrium, an eigenvalue which has a real part of zero indicates a possible bifurcation. Hence, solving the characteristic equation A M + XM(-0 + 1 + a + M2TT2(DX + DY)) + (-/? + 1 + M2TT2DX) ( a + M V A , ) + a0 = 0 for 0 while setting XM = 0 determines the critical value of P for a bifurcation in the M t h mode. 0 = (-P + 1 + M2TT2DX) (a + M V A , ) + aP (1 + M2ir2Dx) (a + M2/K2Dy) PM = P = M27T2Dy If we plot P versus M, the minimum value will give us the value of M that corresponds to the first mode that will bifurcate as we increase p. Figure 3.2 shows that 10 is the integer value of M that gives the lowest value for p. Hence, we expect that the 10 t h mode will bifurcate first. We get our first bifurcation at P = Pw = 8.397702, and the homogeneous equilibrium solution for p < Pio is asymptotically stable. 3.7 Centre Manifold Reduction Now that we know where the first bifurcation takes place, we can analyse it to find out what kind of a bifurcation it is, and what the nearby solutions look like. A n informative entity to derive at this point is the centre eigenspace, Ec. This is the one-dimensional space that corresponds to the eigenvector of the eigenvalue with real part equal to zero. To find the centre eigenvector, we simply solve LM{PM)$ = 0 , Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 25 Figure 3.2: The graph shows the relationship between 0 and M when A M = 0. The lowest point of the graph that corresponds to an interger value of M occurs at M = 10, thus indicating that the 10th mode will bifurcate first as we increase 0. where LM(0) is defined as L*M(P) = (0-1)- M2ir2& a giving us . 1 , M21T2Dy Dy J 0 -a - M2-K2Dn cos Mirs . To study the dynamics of the system at the bifurcation, we need to project the solutions onto the centre eigenspace, Ec. To achieve this, we project along the infinite-dimensional stable eigenspace, Es. The Fredholm alternative tells us that the stable eigenspace is characterized as all those vectors orthogonal to the null eigenvector of the adjoint problem. The centre Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 26 eigenvector of the adjoint operator is <&*, the non-trivial solution to L*m{(3MW = 0 , where L * M is defined in equation 3.12. The solution (up to multiplication by a non-zero scalar) is j + M2n2Dy cos MTTS . It will be useful to select a vector 3?* parallel to in such a way that the inner product (<£, $*) is normalized to unity. To do this, we define <&* = and assert -M2-n2DT. — \ M2TT2DV cos MTTS, 1 M2ir2Dy cos MTTS ) = 1 A f V D y M 2 T T 2 L > x + 1 Thus, ' M2TT2Dy M2TT2Dx + l 1 H a M 2 v r 2 A , M2TT2DV 1 + M27T2Dy 1 cos MTTS . The solutions to our system can be represented as the sum of two vectors taken from complementary subspaces of X[0,1]. We will represent our solution, u, as u(s) = <p(s) + ip(s) where $ e Ec = span{$} and ip G Es = This decomposition will help us later to isolate the dynamics on the centre manifold, since the centre manifold is parallel to Ec at the origin. We have determined that the first bifurcation occurs in mode M = 10 at the parameter value of P = Pio- In order to extract the dynamics at the bifurcation, we need to include the Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 27 non-linear terms of our system: at du&> L dt J -0 a v.™ 1 + - 1 or - u = L(0)u + N(u,0) (3.17) LM = L(0M) = 9 ' x ds'2 a To simplify things, we translate our bifurcation parameter to the origin with the change of variables \x — 0 — 0M- Then we can separate our linear operator into two operators: the constant, singular matrix operator r 0M - 1 + Dx -0M and the varying parameter-dependent matrix 1 0 - 1 0 and hence L(0) = LM + fJ-G. Using the fact that LM</> = 0, our system now takes the form —c?+— if = LM4> + LM$+ l*G4> + HG$+N{$ + $,IJL) (3.18) at at = LMj+uG<£+u.G-ip +N($+->p,u.) . (3.19) Now we can use our inner product to make a projection. To observe the dynamics on the centre manifold near the bifurcation, we project our system onto the centre eigenspace, Ec. We also project our system onto the complementary subspace, Es. We define the projection operator, ^ , by flG = jJL P(u) = ({*,$*} $ . Then P defines a projection of X[0,1] onto Ec, and (I — P) defines a projection onto Es. Under the operations of P and (I — P), our system of equations becomes P (I-P) dt d_ dt <f> = aPG((j> + V) + PN((f> + V, M) , if = LMif+IJ-G<t-iJ.PG(f+ nGip'- uPGtp (3.20) (3-21) Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 28 +N($+$,ii)-PN$+$,ii) , where we have used the fact that P commutes with LM- We augment 3.20-3.21 with d dt H = 0 . Since our centre eigenspace is only one-dimensional, we can parameterize it with one scalar variable. Hence, we let (p(s) = x<&(s). Then, equation 3.20 can be written dx dt $ = nPG{x§ + V) + PN(x§ + ip, u) or dx ~di = (nG(x$ + if) + N(x<S> + i£,u) ,$*) . (3.22) Equation 3.22 contains the dynamics of the system along the vector <fr. Recall that <fr is in Ec, the centre eigenspace. However, equation 3.22 contains ip. Using the invariance of the centre manifold, we can derive an approximation for ip. The centre manifold is locally the graph of, if = H{<j>,p) , = H(x$,u), = h(x,u) . Since u = 0 is a solution for all u, we must have ip = 0, so h(0, u) = 0. Also, the stable x=0 manifold has no component in the direction of the centre eigenspace, so -g^ip = 0. Thus, h(x,0) = O ( |x | 2 + \x\ huxp + h,2ox2 + O (\x\ Iji]2 + \x\2 \p\ + | x | 3 > h^\xp + h*2QX2 + ... h^xu + h^x2 + ... As it turns out, we only need foo to determine the bifurcation. Thus we let u = 0. We substitute h = h(x,0) into equation 3.21, with u = 0, to get dh dx dx dt = LMh + N(x$ + h,0)- PN(x$ + h,0) . Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 29 We can use equation 3.22 to substitute for (2/120Z + ...) ( i V ( s $ + /i,0) , $* ) = LMh + N (x$ + h,0) - (/V* (x$ + h, u) , $ * ) $ (3.23) We will go through equation 3.23 and write each term as a series in x. We start with a;$(i) + x $ ( 2 ) + h(2) + h$X2 + . . . N (x$ + h,o) = N ,0 = N ,0 Recall that for 0 = 0M I 1 N r \ - 1 So (see appendix for details) PMKUW2 + 1UWUW + O-UW2UW) (x2N2 + x3N3 + 0(\x\4)) , (3.24) where Then N2 = / 3 m K $ ( 1 ) 2 + 7$ ( 1 ) $ ( 2 ) , 7V3 = 2 /3 M «* ( 1 ) 4o + 7* ( 1 ) 4o } + 7 * ( 2 ) 4 o ) + a S ^ V ) . + = 2y ' _} (x 2iV2 + x 3/V3 + o(|a;| 4)) • $* ds = 2x2 £ _] • $*JV2 ds + 2xi £ j • $*/V"3 ds + O (\x\4) Note that the Fourier expansions of $(s) and h^s) are defined as N (3.25) $ ( i ) ( s ) $(2)( s ) <p(2) /l20(s) = 4 o ^ ) L(2) cos MTTS + E m = l 4o ,m COSmTTS + E m = l 4 20,m COSmTTS Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 30 Because of the orthogonalities of the Fourier cosine series (see appendix for details), Jo $*N2 ds = 0 . (3.26) Next, we evaluate the integral in the x 3 term of equation 3.25 (see appendix for details): 2x3 ll [A' **N3 ds=x" (Ah^°+Bh™°+chi^2M+m™'2M+E)' (3-27) where *(i) ) A = + ^7<P(2)) ( C = (f3MK^l) + ^ ( 2 ) ) ( V ( 1 ) - <p*(2)) , (3.28) D = > ( 1 ) - ^*(2)) Thus, we now have (TV (x* + h,6), <&*) = x 3 ( A ^ o + + C / i ^ + Dh%M + E) + O ( |x | 4 ) . (3.29) Equation 3.23 now becomes (2K2oa; + . . . ) ( x 3 (Ah% + Bh{2l0 + Ch%M + Dh%M + E)+0 ( | X | 4 ) ) (3.30) = LMh + N ( x$ + h,0)- ( x 3 ( A / I ^ O + + Ch%M + £>4o,2M + E) + 0 ( W ) ) $ Before we do too much work, we should look at equation 3.22 to see what we need to find. After all, equation 3.22 contains the actual dynamics on the centre manifold. We have, d X = (iiG(x$ + j)+N(x$-rj,li) dt P ,$*^ + (N(x$ + h,p),$*) , or dx Ht = 2u f1 Jo $* ds a;$(i) + (^i) -X$W - fed) +x 3 (Ah™Q + Bh% + Ch%M + Dhfl2M + E)+0 ( |x | 2 \U\ + | x | 4 ) (3.31) Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 31 2/x Jo $* ds , The remaining task involves determining the coefficients of the xfi and re3 terms in equation 3.31. Keeping in mind that h = O (|a;|2), we start this by evaluating the integral, rl rl = 2u.j xtpW (ip*(1)-ip*{2)) cos2 MTTS ds + J h^M - ^* ( 2 ) ) cos2 MTTS ds , = W X ) ( V ( 1 ) - </ ( 2 ) ) + O ( |x | 2 ) . (3.32) Expression 3.32 is the xu. term of equation 3.31. A l l the coefficients are known. To completely determine the dynamics on the centre manifold near the bifurcation, we still need to determine ^20 o> ^ 20 o> 2M a n d h$0 2M- We c a n obtain values for these variables by equating the coef-ficients of the terms involving powers of x from equation 3.30. The lowest order terms are of second order, and equating their coefficients gives us (see appendix for details), £ M / I 2 0 + N2 0 PM - 1 + T>x-§^i a -P. M '2 0 L(2) l2 0 = - i ( / W 1 ) 2 + 7^ (V 2 )) (3.33) (1 + cos 2MTTS) If we now equate the coefficients of the different cosines in equation 3.33, we can determine values for the Fourier coefficients of the expansions of and h^2\ Equating the coefficients of the constant terms gives us, -PM a —a 2 A2 0,0 2 A2 0,0 h (1) 2 0,0 l ( 2 ) l2 0,0 = - ^ ( W 1 ) 2 + 7^ ( 1V ( 2 )) = - ( P M K ^ 2 + ^ { 1 V 2 ) ) Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 32 Then we equate the coefficients of the cos 2MTTS terms. pM - 1 - 4M 2 7r 2 £> : a -P. M -a - 4 M2 T T 2 n (i) 2 0.2M (2) 2 0,2M = ~ ( / W 1 ) 2 + 7^{1V2)) 1 - 1 (1) 2 0,2M (2) 20,2M (/W 1 ) 2 + 7^(1V2)) PM - I- 4 M2 T T 2 D X a -PM -a - 4M2ir2Dy We have determined everything we need to evaluate the coefficients in equation 3.30. We will also use the eigenvectors, 1 -1 _ Ac M 2TT 2DV DV cos Mns , <!>* = (, i M V D v M2n2Dx + l l a M27T2Dy - l 1 + M^Dy cos MTTS , to get numeric values for the coefficients in equation 3 .31, which dictates the form of the bifurcation on the centre manifold. It can be written as dx = ciux + c2x3 + O (\x\4 + \x\2 |m|) • dt Let us work out the coefficients c\ and c2. From equation 3.32, we have (see appendix for details) c i = 0 .392309 . (3.34) The pertinent values of the coefficients of h are, = 0 ^2 0,0 h (2) 2 0,0 h{1) "2 0,20 h (2) 2 0,20 = 0 .068052 = - 0 . 2 6 8 9 8 = 0 .030305 Now for c2: C2 Ah, ( 1 ) +Bh% + Ch%, + Dh^ 2 0,0 2 0,20 + E Chapter 3. Mathematical Analysis of the Brusselator's Turing Bifurcation 33 = 1.166704o,0 + 14.1231342o,o + 1.1667o4oj2o + 14.123124o,20 ~ 1-59585 = -0.520567 After all that hard work, the equation for the dynamics on the centre eigenspace is given by, = 0.392309//X - 0.520567x3 + O (|x|4 + \x\2 . (3.35) This is a pitchfork bifurcation. Equation 3.35 gives the leading-order behaviour of the dynamics on the centre manifold. That is, it gives us an idea of how far along the centre eigenspace our equilibrium solution is for a given value of u. Figure 3.3 shows a schematic diagram of the existence and stability of equilibrium solutions for our pitchfork bifurcation. X X < 0 ,-!"""' / » * > ° i M-\ ; X<0 x > 0 Figure 3.3: (Pitchfork bifurcation) The arrows indicate the direction of the evolution of the solution as t increases. For u < 0 (Q < 0w), the trivial solution x = 0 (spatially homogeneous) is asymptotically stable. For u > 0 (0 > 0\Q), the trivial solution is unstable, and there is a pair of nontrivial, asymptotically stable solutions x « ± . / ^ . The nontrivial solutions correspond to a spatial pattern. Chapter 4 Numerical Analysis of the Brusselator Bifurcation analysis is excellent for determining the dynamics of a system in the neighbourhood of an equilibrium solution, indicating the stability of an equilibrium solution, the approximate location of subsequent equilibrium solutions, and the basic shape of the centre manifold in the neighborhood of a bifurcation. However, we also wish to understand more of the global dynamics of the Brusselator system. A numerical package called A U T 0 8 6 is useful for this purpose. 4.1 Numerical Bifurcation Analysis Software Package: A U T O 8 6 A U T 0 8 6 is a software package, written by E . Doedel [7], which tracks equilibria and periodic solutions of systems of differential equations as a specified parameter is varied. This process is called continuation. The user edits a F O R T R A N subroutine that defines the system of differential equations with parameters. A n initial equilibrium solution is also entered into a subroutine. A U T 0 8 6 increments the parameter and uses a nonlinear algebraic solver to converge onto the equilibrium solution that corresponds to the new parameter value. If the solver does not converge, the parameter increment is halved, and the solver is invoked again. This process is continued until either a solution is found, or the parameter increment drops below a user-defined lower bound. At each parameter step, output is generated and the system's eigenvalues are calculated. A bifurcation is indicated when the sign of an eigenvalue changes. The version of the Brusselator system that was implemented into an A U T 0 8 6 function is that from equation 3.1, involving parameters a and b. The parameter b is our bifurcation parameter. A system of ordinary differential equations is obtained from the Brusselator by a 34 Chapter 4. Numerical Analysis of the Brusselator 35 spectral Galerkin approximation. We replace elements of X[0,1] with truncated versions of the Fourier series in equation 3.6. For example, if we take only the first two terms of the Fourier series, then X and Y are approximated by X(s,t) = ^ ^ + X 1(t)cos7rs , Y(s,t) = + Y1(t) cos TTs . Using the above definitions, our system of partial differential equations becomes, 1 dXo dX\ 2 dt + dt • COS ITS a — (b + d) + Xi cos 7rs) ' X +c [ — + X l COS TTS 2 'Y« 0 \ 2 + Y\ COS TTS — 7T DXX\ COS 7TS , ldYndY1 J X ° ^ V — COS TTS — b ( h A l COS ITS 2 dt dt V 2 (4.1) — C + -^1 C 0 S 7 r S ) + ^1 C O S 7 r s ) — TT2DyYi COS 7TS which, in vector form, can be written, (4.2) where ^1+X1(t) COS 7TS ^ + Yl(t)cOS7rS and f(X,s) is the right-hand side of the system in equation 4.1. To separate the four time derivatives into four different equations, we apply the inner products, 1 \ (• 0 0 I • \ • 1 and COS TTS COS TTS to equation 4.2, and simplify the integrals using the orthogonality properties of the cosines. Chapter 4. Numerical Analysis of the Brusselator 36 First, we will focus on the left-hand side of the system (see appendix for details). Likewise, we have, dx dt - X , 1 0 0 1 dX0 dt dYQ dt (4.3) dx COS7TS ) = dXi ~Tt (4.4) dx Tt*' COS TTS ) = dY! dt We still need to apply the inner products to f(X,s), the right-hand side of the system of equations (see appendix for details). 1 0 Evaluating the other inner products, we get, f(X,s), = 2a-(b + d)X0 + c ( ^X2Y0 + XQXM + ^X2Y0 (^f(X,s), = bXQ-c[ ^X2Y0 + XoXrY, + ^X2Y0 ) , (4.5) f(X,8), COS TTS - _ (b + D + TT2Dx) X, + c Q x 2 ^ + ^X0XXY0 + \xlY^j , f(X,s), COS TTS )=bXi- TT DyY\ - c [ -Xfa + " X o X i F o + -jXiY! As a result of applying the inner products, we are left with four ordinary differential equations. dX0 dt = 2a - {b + d)X0 + c (^X$Y0 + XQX^ + \xfY0^ Chapter 4. Numerical Analysis of the Brusselator 37 ) This format is suitable to be entered into a subroutine to be utilized by AUTO86. The following numerical results are for the system in which the Fourier series are truncated after the 2 2 n d mode. The homogeneous equilibrium solution, along with the associated low 6-value of b = 0.417, was input to the appropriate A U T 0 8 6 subroutine. The remaining set of parameter values is given in table 2.1. As A U T 0 8 6 increases 6, it reports a bifurcation at b = 0.419885121. Our analysis predicted a bifurcation at b = 0.419885110, in agreement to seven digits. The final test of the validity of the bifurcation analysis is to compare resulting bifurcation diagrams: the one produced by A U T 0 8 6 , and the one we arrived at analytically in chapter 3. Figure 4.1 shows the two bifurcation plots superimposed for comparison. Since the analytical bifurcation plot is only accurate up to order (\x\ \u\ + \x\3), we only expect the graphs to agree near the bifurcation. 4.2 Checking A U T 0 8 6 Output with br.c We can compare the solutions calculated by A U T 0 8 6 with those by br.c (see section 2.1). The spatial domain [0,1] in br.c is discretized into 100 elements. The Fourier coefficients are periodically saved to a file, and then plotted. Setting b = 0.435, we let br.c run for 800,000 time steps of 0.002 time units, after which the solution calculated by br.c seems stationary, indicating that it is likely close to an equilibrium solution. If we compare the amplitudes of the dominating 10 t h Fourier' cosine modes in the solutions returned by the two methods, the relative error is less than 5%. Figure 4.2 shows the solution trajectory generated by br.c on top of the bifurcation diagram. The solutions themselves are compared in figure 4.3. The plot was generated by taking all 46 Fourier coefficients calculated by A U T 0 8 6 , and plotting the resulting function against the Chapter 4. Numerical Analysis of the Brusselator 38 0.6 O O' V AUT086 O Analysis — -0.6 O. O 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 b Figure 4.1: Superimposed are bifurcation diagrams produced by A U T 0 8 6 and by the bifurcation analysis of chapter 3. The vertical axis represents the amplitude of the 10 t h mode for morphogen X. The curved branches correspond to non-homogeneous solutions. Near the bifurcation, the curvature of the graphs seem to match quite closely near the bifurcation. However, the graphs diverge from one another for parameter values that are far from the bifurcation value. This is expected because the analytical curve is only leading-order accurate. discretized solution produced by br.c. It is interesting to view the other bifurcations, and the role they play in forming the solution set to the Brusselator. For example, the solutions on the branches that were shed from the pitchfork bifurcation that we analysed in chapter 3 are stable. However, A U T 0 8 6 indicates that there are other stable solutions on branches that result from bifurcations at higher parameter values. Thus, multiple stable equilibrium solutions exist for certain parameter values. Figure 4.4 is a stereogram of a portion of the bifurcation diagram involving several bifurcations, one of which is the pitchfork bifurcation that we analysed. The branch of diamonds that extends toward the right-hand-side indicates equilibrium solutions in which the 10 t h mode has amplitude zero, while the 9 t h mode has nonzero amplitude. This suggests that, for appropriate values of b, not only do stable solutions of 5 waves on [0,1] exist, but stable solutions of 4.5 waves on [0,1] also exist. Similarly, increasing b further, the 11 t h mode loses its stability, branches off the Chapter 4. Numerical Analysis of the Brusselator 39 0.5 u (i) 10 -0.5 00000° ooo o o o Stable O Unstable + br.c • I I I I I II I I I I 1 I I I I I I I I I | l I I I I I I I > 0 < >oooo4 ooo 0.415 0.42 0.425 0.43 0.435 b 0.44 0.445 Figure 4.2: The bifurcation diagram and solution trajectory from br.c are superimposed. The solution trajectory output by br.c (for which b = 0.435) starts close to the unstable, homogeneous equilibrium, and evolves towards the stable 10 th-mode dominated solution. homogeneous branch, and bifurcates again, shedding a branch of stable solutions which have 5.5 waves on [0,1]. Note that in figure 4.4, a solution in which the 9 t h and 10th modes have zero amplitude will look like a homogeneous solution since the axes indicate only the 9 t h and 10th modes. • • Some insight may be gained from the dispersion curve, A(m) (shown in figure 4.5). The dispersion curve plots the relationship between a Fourier mode number and the greater of its two eigenvalues. Each Fourier mode has associated with it two eigenvalues (since the original system is two-dimensional). The range of parameter values investigated in this project admits strictly real eigenvalues, so the greater of the two eigenvalues is simple to determine. If, for a particular Fourier mode (i.e. integral value of m), the graph is negative, then both of its eigenvalues are negative and the homogeneous equilibrium is stable with respect to perturbations involving that mode. However, if the graph is positive at a particular Fourier mode, then that mode has at least one positive eigenvalue, and we may conclude that the homogeneous equilibrium is unstable. Moreover, near the homogeneous equilibrium, the mode with the highest eigenvalue Chapter 4. Numerical Analysis of the Brusselator 40 1.5 1 br.c X — AUTO86 X • • •_• br.c Y AUT086^Y 1 h Cone. 0.5 h 0 0 0.2 0.4 0.6 0.8 1 s Figure 4.3: The 10 -mode dominated equilibrium solutions determined by A U T 0 8 6 and br. c for b = 0.4361 are superimposed. The parameter values at which A U T 0 8 6 generates output can not accurately be controled. As a result, the closest value to 0.435 at which output was generated was 0.4361. The br.c simulation, which also used b = 0.4361, was started at a random state close to the homogeneous equilibrium. has the greatest potential for growth. Recall that we chose to vary bB rather than aA as our bifurcation parameter. One of the factors in our decision was the speed at which the system converged to its final steady-state. For some reason, when we used parameter values from region £2 (see figure 2.4), the system converged to its equilibrium much more slowly than if we had used parameter values from region S. We can use the dispersion curves now to see why. Both shaded regions (£2 and £ ) have dispersion curves which have a positive local maximum, and satisfy The difference between the two regions is that region £2 satisfies A(0) > 0, while region £ satisfies A(0) < 0. Figure 4.6 shows typical dispersion curves for parameters from regions £2 and S. Perhaps having an isolated set of unstable modes, bounded away from zero, allows for a smoother convergence to the equilibrium state. Perhaps the fewer the unstable modes, the m lim A(m) < 0 . oo Chapter 4. Numerical Analysis of the Brusselator 41 Figure 4.4: This is a 3-D stereogram of a bifurcation diagram. Hold the page about 20 cm from your face and relax your eyes into a blur. You will see 3 blurry images. Try to focus on the image in the middle. The axis going almost-straight into the page is the 6-axis. The vertical axis is the U^Q axis, and the remaining axis, rising slightly out of the page and to the right, is the axis. Notice which branches are stable, and which are unstable (the stable solutions are labeled with a diamond, while the unstable solutions are labeled with a + sign). less competition there is between modes, and the more quickly the system establishes a balance between these competing modes. Are the solutions suggested by AUT086 's bifurcation diagram accurate for equilibrium solutions far away from the homogeneous equilibrium? The rationale behind considering a finite Fourier series as an approximation to the true solutions is in the insignificance of the non-linear terms for solutions close to the homogeneous equilibrium. For such solutions, the amplitude of the various Fourier modes (other than, perhaps, the zero-th mode) is close to zero, and hence a non-linear combination of such modes would be even smaller. However, when we consider solutions that are not close to the homogeneous equilibrium, we do not have a guarantee that the non-linear cross-terms are small, and hence truncation of the Fourier series may lead us to spurious solutions. Nonetheless, a simulation can build some evidence in favour of the existence of these stable, non-homogeneous solutions. Figure 4.7 shows a solution on the far right-hand branch of stable solutions of figure 4.4. After reaching what seems to be an equilibrium solution, the solution remains 9 t h-mode dominated. Chapter 4. Numerical Analysis of the Brusselator 42 Figure 4.5: For b = 0.435, we see that modes 9 through 11 have positive eigenvalues (the curve is above the zero axis)(curve is below the axis for the 8 t h and 12 t h modes). This tells us that, at the homogeneous equilibrium solution, the 9 t h , 10 t h and 11 t h modes are unstable. The shape and location of the dispersion curve is a continuous function of the parameters. Chapter 4. Numerical Analysis of the Brusselator 43 b) Figure 4.6: Graph a) shows the dispersion curve for the parameter values in table 2.1, with the exception of b = 0.43. These parameter values correspond to a point in region £ in the Lacalli-Harrison parameter space, and hence A(0) < 0. Graph b) shows the dispersion curve for the same parameter values, with the exception of a = 0.01 and b = 0.25. This parameter set falls in region Q of the Lacalli-Harrison parameter space, and hence A(0) > 0. Chapter 4. Numerical Analysis of the Brusselator 44 Cone. 0 0.2 0.4 0.6 0.8 1 s Figure 4.7: Superimposed are plots of the equilibrium solutions calculated by br. c and A U T 0 8 6 for b = 0.4361. See figure 4.3 for comparison. These stable equilibrium solutions are 9 t h-mode dominated, demonstrating the presence of multiple branches of stable equilibria. Choosing a random initial state for the br. c propagation usually causes the system to converge to the 10 t h-mode dominated solution. Instead, the shown solution from A U T 0 8 6 was used as the initial solution for the br. c run. Chapter 5 Non-Homogeneous B-Concentration We have seen that pattern can be obtained with a homogeneous distribution of B. However, as was mentioned in chapter 1, a gradient may produce distributed pattern: creating pattern where the chemical concentration is above a threshold, and no pattern where the concentration is below. In chapter 2, we discovered that a monotonic gradient in the parameter group bB produced pattern corresponding to the cohesive zone during somite formation (see figure 2.5). Here, we further explore the effects and solutions arrived at with a gradient in the concentration of B. Without loss of generality, we can assume that 6 = 1. Since b = bB, we can represent a gradient in B as a gradient in b, implemented as a spatially-varying parameter, b = b(s). From this point on, we refer to the parameter b as if it were a chemical with a spatial distribution b(s) and time dynamics. A point source of b creates a gradient. If we also assume that b decays over time at a constant rate, then the partial differential equation governing its behaviour is, where k is the decay rate and Db is the diffusion rate. The equilibrium solution for b(s) is b{s) = b{0)e V D* . For now, we only need the equilibrium solution for b. We will implement time dynamics on the distribution of b at a later stage. 45 Chapter 5. Non-Homogeneous B-Concentration 46 5.1 A U T O 8 6 Solutions Involving a 6-Gradient Our representation of the equilibrium solution for b(s) has three parameters: 6(0), k, and Db-We can replace these parameters with other, more representative ones. Instead of referring to 6(0), we will define the parameter 6mid to be equal to b(^). We can lump the two parameters k and Db into one parameter called R, which will represent the severity of the exponential drop-off. R is defined as V A. Changing 6mid moves the midpoint of the curve 6(s) up or down, while changing R simply tilts the curve about its midpoint. The gradient is like a see-saw, with 6mj<j acting as the fulcrum, and R affecting the slope. The full function, 6(s), is b(s) = bmide-R(s-l2) (5.2) If R = 0, then b(s) is homogeneous. The greater R is, the more steep the gradient. If R > 0, the equilibrium solution of 6(s) is decreasing on [0,1]. To keep the functional representation in our A U T 0 8 6 modelconsistent, we expand b(s) into a Fourier cosine series, and truncate it as we did for the X and Y morphogens. Hence, oo i = l where = 2 f\(s) Jo cos iizs ds We can determine the coefficients explicitly. 26m;d / R 6 0 = _ ( e a - c a ) k = 2e* — {R — e R (—Rcosin + insmi-iT)) R2 + ih The new definition of 6 = 6(s) can be substituted into the system of differential equations that is given to A U T 0 8 6 , such as that shown in equation 4.1. Note that the additional Fourier modes introduced by b(s) will change the outcome of the inner products that are applied to Chapter 5. Non-Homogeneous B-Concentration 47 achieve a system that is suitable for A U T 0 8 6 input. Spatial variability in b breaks the symmetry of the system's bifurcation diagram. Figure 5.1 shows a bifurcation diagram similar to that in figure 4.1. The branch of homogeneous equilibrium solutions that would have undergone a pitchfork bifurcation avoids this destiny. Instead, it simply flows into the lower branch, while the upper branch loops down, losing its stability at a limit point. 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 u (i) 10 -1 1 1 1 1 1 1 I f Stable O Unstable + - o — -1 1 i i i i i 0.38 0.4 0.42 0.44 0.46 0.48 &mid 0.5 0.52 0.54 0.56 Figure 5.1: Bifurcation diagram of the broken symmetry caused by the spatially non-homogeneous b distribution. This diagram is similar to figure 4.1, except that R m 0.0029 (as opposed to R = 0, as depicted in figure 4.1). The intersection of branches that marked the pitchfork bifurcation in figure 4.1 is broken, leaving two separate branches. The branch of homogeneous solutions that once continued through to the supercritical side of the pitchfork bifurcation is now broken in many places, but still serves to bridge the gaps between pairs of branches of non-homogeneous solutions. Finally, figure 5.2 is a stereogram similar to that of figure 4.4, except with a non-homogeneous (exponential) 6-distribution. Pitchfork bifurcations are no longer present. A l l the branches are non-intersecting. The symmetry is broken. Chapter 5. Non-Homogeneous B-Concentration 48 Figure 5.2: This is a 3-D stereogram of a bifurcation diagram. See caption of figure 4.4 for viewing instructions. Notice that the pitchfork bifurcations shown in figure 4.4 are not present here, leaving separate branches, indicating the symmetry-breaking action of the spatial gradient in b. 5.2 br.c Solutions Involving 6-Gradient How well do the equilibrium solutions of br. c and A U T 0 8 6 agree in the presence of a b-gradient? To check, we perform a series of simulations with br.c for which R = 0.0029, incrementing 6m;d slightly for each run. Once we have found the equilibrium solution for a given value of 6mid, a small increment should move the corresponding equilibrium solution only slightly. We would expect the equilibrium solution from the previous simulation to be in the basin of attraction of its corresponding equilibrium with the slightly incremented 6mid- Hence, we use the previous equilibrium state as the initial state for each simulation. Figure 5.3 compares the bifurcation diagram output by A U T 0 8 6 with an equivalent graph constructed from a number of br.c simulations. Only the stable states are found by br.c. It should be noted that, as 6mid approaches the bifurcation value of the 6-homogeneous problem, convergence to the non-homogeneous equilibrium solutions is very slow. Chapter 5. Non-Homogeneous B-Concentration 49 0.8 Stable O Unstable + br.c (stable) • -0.8 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 Figure 5.3: Shown here is a comparison of bifurcation diagrams produced by br. c and A U T 0 8 6 . Their correlation is good. However, close to the bifurcation point of the corresponding homo-geneous problem (6m;d « &io), the convergence in br.c is slow. 5.3 Biological Implications of 6-Gradient Without a gradient in any of the parameters (i.e. homogeneous b distribution), the bifurcation diagram is symmetrical (as seen in figure 4.1). At appropriate parameter values, for every equilibrium solution that has a peak at s = 0, there exists a corresponding equilibrium solution that has a trough at s = 0. The basins of attraction for these two stable equilibria must also be symmetric. So, for an initial state that is close to the homogeneous (unstable) equilibrium, the final steady state to which it converges is not clear. This ambiguity in convergence behaviour brings doubt to the biological appropriateness of the Brusselator model. As we just witnessed, adding a gradient to b breaks the symmetry of the bifurcation diagram, causing branches that would have undergone pitchfork bifurcations to break into either non-bifurcating branches or into branches with limit points. The symmetry breaking action of the b gradient may play an important role in determining which particular equilibrium solution a dynamic version of this system converges to [15]. Chapter 5. Non-Homogeneous B-Concentration 50 The need to consistently converge to one equilibrium, and not another, might be present in the initialization of somite production. The location of the first peak in X concentration corresponds to the location of the first somite. The difference between starting this process with a peak in X and starting it with a trough is one half of a wave-length, or one half the width of a somite. This offset in positioning may be critical to the development of an embryo. Consider the case in which it is necessary to converge to the equilibrium in which the X concentration exhibits a trough adjacent to the source of b. How can we avoid converging to a solution with a peak in X at the b source? In other words, how can we guarantee that our system converges to a specific (super-critical) branch of the bifurcation diagram in figure 5.1? Start with a sub-critical, graded distribution of b (i.e. b has a gradient, and is below the bifurcation value throughout the cohesive zone). Our equilibrium solution for X is near-homogeneous, and corresponds to a point on the far left-hand side of the bifurcation diagram in figure 5.1: a point on the branch of stable equilibrium solutions. This branch of the bifurcation diagram is continuous, and is connected to the lower branch of the diagram, which corresponds to 10 t h-mode dominated non-homogeneous equilibrium solutions. The fact that the branch of near-homogeneous solutions is connected to a particular branch of non-homogeneous solutions is biologically exploitable. It affords a developing embryo an algorithm for choosing one final steady-state over another. If the X and Y morphogens in the cohesive zone of the embryo are able to converge to their steady-state solutions fast enough, then slowly increasing the production of b will cause the X and Y solutions to stay in (or close to) their steady-state solutions. As the b distribution rises, the X and Y steady-states will continuously flow onto the lower branch of the bifurcation diagram, producing a trough in the X distribution adjacent to the 6-source. Chapter 6 Dynamic Production of Somites: Cohesive Zone Propagation Section 1.2 briefly described the process of somite formation as a dynamic process in which the production of somites takes place in the cohesive zone while the cohesive zone slowly sweeps across the length of the embryo, leaving a trail of mature somites in its wake. In this chapter, we model the motion of the cohesive zone in conjunction with the production of somites. 6.1 Cohesive Zone Boundaries We present a candidate mechanism, involving the b gradient, its source, and its effect on the cohesive zone boundaries, to explain the motion of the cohesive zone. 6.1.1 Known Mechanisms The cells of the cohesive zone are connected by gap junctions that allow intracellular substances to diffuse from one cell to another. A somite is a cluster of these cells. Upon maturation of the somite, the cells break all their gap junctions with other cells except for those cells in the same somite. Hence the cells in the somite detach from the cohesive zone, and no longer contribute to the cohesive zone's intracellular diffusion continuum. 6.1.2 Assumed Mechanisms The action of creating a new somite (by breaking gap junctions) defines the location of the boundary between the segmented zone and the cohesive zone. The source of b (and the interface between the segmented zone and the cohesive zone) is placed at the posterior edge of the newly matured somite, hence shortening the cohesive zone by a somite length. 51 Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 52 We are assuming that the degree of development of a somite at some location is an increasing function of the concentration of the X morphogen. Also, suppose that there is an upper threshold value for the concentration of the X morphogen, denoted Tx, that, when reached, matures the corresponding somite. The other end of the cohesive zone (the cohesive/non-cohesive interface) is defined by a different kind of threshold in the concentration of b. A cell in the non-cohesive zone remains as such until the local concentration of b rises above a certain threshold. Essentially, a cell needs a certain critical concentration of b before it can form gap junctions with its surrounding cells. We will call this critical concentration T/,. Once the gap junctions are established, a decrease in the local concentration of b below T& does not break the gap junctions. For this reason, the boundary between the cohesive and non-cohesive zones can only move in the direction that increases the length of the cohesive zone. Somite formation is the only way to decrease the length of the cohesive zone. In section 1.3, we made the assumption that morphogens X and Y are intracellular, and thus restricted to the diffusion continuum of the cohesive zone. 6.2 Changes in the Basic Model With the introduction of a moving cohesive zone, we must expand the spatial domain of our simulation in order to represent each of the three zones. It is irrelevant how long the spatial domain is, provided it is long enough to accommodate the cohesive zone and a sample of its motion. In our simulations, we define the spatial domain to be 2.5 units long. Initially, the cohesive zone will be approximately one spatial unit in length. 6.2.1 Dynamic 6-Concentration The movement of the cohesive zone forces us to solve the dynamic 6-concentration on the spatial domain, just as we do the X and Y concentrations. Once the source of b moves, the solution for b will undergo a transient. In br.c, the solutions for X and Y are spatially discretized and Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 53 the system is solved using finite difference numerical methods. The same is done here for b. There is a fundamental difference between the way X and Y can be modelled and the way b can be modelled. The morphogens X and Y are intracellular, while b is free to move throughout the three zones in the extracellular fluid. As a result, its boundary conditions are different. The boundaries of the cohesive zone allow no flux of X or Y, but are transparent to b. Thus, b has no boundaries. However, if we wish to simulate the dynamics of our system, we must find a finite rendering of it. A computer cannot represent infinity. Effectively, we have to create boundary conditions for b. One boundary comes from the fact that the 6-solution is symmetric about its source. For the other boundary, the best we can do in this situation is make up a condition that mimics what we would expect with no boundary. The boundary condition of our finite-difference implementation simply depends on how we define the state of the first point outside our domain. Suppose there are n elements in our spatial domain. Then, the boundary condition dictates how we define the state of the (n + l ) t h element. Alternatively, how we define the state of the (n + l ) t h element dictates what the boundary condition should be. Since we expect the solution for b to be near the equilibrium, it seems reasonable to let the first point outside the domain obey the exponential decay that we would expect if the solution actually were in equilibrium. In other words, if our equilibrium solution is then, at some point s* in the spatial domain, b(s*) = b(0)e D b and, thus, The ratio of the equilibrium 6-concentrations between any two adjacent elements is Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 54 Thus, we define b (2.5 + £ ) to be 6 (2.5 + - ) =6(2.5)e' n Since the concentration of b is solved dynamically, we need its dynamic parameters as well as the parameters that define its steady-state. The concentration of b at its source is kept constant. Thus, our dynamic b gradient is represented by three parameters: the concentration at the source (bs = 6(0)), the decay rate (k), and the diffusion rate (Df,). 6.2.2 Source Placement in br.c Determining when the solution of the X morphogen exceeds the threshold, and where to put the new source, needs to be well-defined. For each time step in br. c, a search for a local maximum is started from the source-end of the cohesive zone, and proceeds posteriorly. A local maximum is defined as the spatial element i, where Once the first local maximum of X is found, its value is compared to the somite-maturity threshold, Tx- If the value X{ exceeds the threshold, then the mature somite is added to the segmented zone by moving the segmented/cohesive zone boundary (and hence the source of b). The 6-source (and hence, the posterior end of the mature somite) is placed at the location of the first local minimum in X posterior to the position of threshold breach. A local minimum is defined as the spatial element i, where Xi-2 > Xi-x > Xi and Xi < Xi+\ < Xi+2 -6.3 Results of Propagation Model For each simulation, an initial state of the system must be given. The initial state includes the initial distributions of X, Y and b, and the locations of both cohesive zone boundaries. The simulations are initialized as follows: X i _ 2 < X{-i < Xi and Xi > Xi+\ > Xi+2 • Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 55 1. The b distribution is started at its equilibrium state, and the source is placed at s = 0. 2. For any location s, the initial value of X(s) is ^ and the initial value of Y(s) is ^ r ^ -Gaussian noise is added to both X and Y. 3. The segmented/cohesive zone boundary (and fr-source) is set to s = 0, while the cohesive/non-cohesive zone boundary is set to the place at which the b gradient drops below the thresh-old Tb. The distribution of X and Y in the cohesive zone change as soon as the time propagation starts. Since the non-cohesive zone contains cells that do not have gap junctions, and thus do not allow intercellular diffusion of X and Y (i.e. diffusion between cells), their initial state has only a local dependence. The above "locally-homogeneous" equilibrium states reflect this situation. A series of simulations was performed to observe the effects of the mechanisms incorporated into the model of the dynamic cohesive zone. The list of parameters that initialize and define the simulation are given in table 6.2. Parameter Value a 0.05 bs 0.61 k 2.58175 x 10"4 Db 0.007 c 1.8 d 0.05 bx i o - 4 Dy 10"3 Tb 0.41 TX 1.55 Table 6.2: Dynamic Cohesive Zone Simulation Parameters. A l l values are exact. Figure 6.1 shows a sequence of solutions output by br.c for the above parameter values. Graphs b) and c) exhibit the solutions just before and just after the creation of a mature somite, and the movement of the source of b. Notice that in graph c) the 6-gradient is undefined to the left of the left-most local minimum in X, indicating the location of the source and the boundary Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 56 of the cohesive zone. Graph d) shows the pattern that results from the Brusselator and the added dynamic properties. This behaviour seems robust as long as the upper threshold, Tx, is low enough that the left-most peak in X can reach it. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 57 Figure 6.1: Brusselator cohesive zone propagation. The time intervals are a) t = 0, b) t = 45, c) t = 48.6, and d) t = 192.75. Notice the Gaussian noise added to the initial X and Y distributions. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 58 The motion of both ends of the cohesive zone plays an important part in the spacing of the somites. The bifurcation analysis in chapter 3 assumes that the spatial domain is [0,1]. If the length of the domain changes, the results will change. We can adjust our analysis to fit this scenario by defining our state space, X[0,^], as X[0,£\ = {u(s) # + E £ i ^ 1 } c o s f S # + £ ~ i * 4 2 ) c o s f S oo -9 £2 2 ™i2 , (2) 2 < OO (6.1) i=i i=i and by altering the definition of our inner product appropriately. The extra factor oi,£ changes the mathematical operators in our system (see equation 3.17) in a continuous way.'.',In our bifurcation analysis, the greater eigenvalue of the 10 t h mode has real part zero, while all the other eigenvalues have a negative real part bounded away from zero (for the parameter values given in table 2.1). Hence, there is a neighbourhood about the given parameter values for which the 10 t h mode is still the first to bifurcate as b (or 0) is increased. Thus, for £ in the neighbourhood of unity, we can expect that our bifurcation results will be qualitatively equivalent. The difference, however, is in the wavelength of the bifurcating mode. The wavelength of the m t h mode is ^ . Given small changes in the length of the spatial domain, only small changes in the spacing of the somites will be reflected. 6.4 Mathematical Interpretation of Propagation Model In order for somite production to continue, the states of the system in the cohesive zone have to repeat similar behaviour for each somite that is produced. In other words, after a somite is produced, the resulting state has to evolve through time and arrive at a state that can produce another somite. Then the whole process happens again, and again, and the cohesive zone propagation and sequential somite production is sustained. The object of our propagation model is to produce near-cyclic behaviour in the state of the system. As discussed above, the redefinition of the cohesive zone (by the movement of the cohesive zone boundaries) actually redefines the state space. In what space, then, does the cyclic be-haviour we require exist? If we include the variable £ in our state, we can link a continuum Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 59 of state spaces X [(),•£] together. The result is a larger state space, X , that includes states over cohesive zones of any positive, non-zero length. We define X as x= UXIM • Consider the state immediately prior to the maturing of a somite. This state is a point in X , which we can call Q\. Once the somite has matured, and the cohesive zone is redefined, our system is at a different state which we will call Q2. The key to the propagation of the cohesive zone is that the orbit 1 of state Q2 needs to be sufficiently close to state Q\. If it is, then we can expect that the somite production process will repeat, fixing the system into a near-periodic orbit of X . The propagation of the cohesive zone comes from the dependable flow from one state Q i , back to a region close to Q\. 6.5 Diffusion Barriers Part of the purpose of this project is to create results that can be tested experimentally on real embryos. We have observations of embryos [1] and some indication of the dynamics involved in their production. If our objective is to test whether the somites are formed as a result of a system similar to ours, then we need to create a test-case that isolates its dynamics from other candidate models. One important feature that our simulation uses is the location of the source of b. What will happen if we insert a barrier to the diffusion of chemical 6? Our next series of br. c simulations involves the insertion of a diffusion barrier into the cohesive zone at some point during the propagation. Wi th our model, we can predict what we might expect from a similar embryological experiment. The outcome may give evidence for the validity of our model. The effect of inserting a diffusion barrier is to instantaneously divide the cohesive zone into two diffusion continuums. The side with the source (anterior to the barrier) was not modelled. In this zone, the resulting no-flux boundary would most likely raise the concentration of b, cause a peak in the concentration of X, and finally form a mature somite. However, the new 'The orbit of a state is the set of states that results from its forward propagation. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 60 somite and associated source would still be diffusively separated from the posterior paraxial mesoderm, and cohesive zone propagation would not be imminent. In the section of the cohesive zone that lies posterior to the diffusion barrier (separated from the source of b), cohesive zone propagation may continue or may terminate. Without a source, the concentration of b will eventually vanish due to its diffusion on an infinite domain and decay with time. However, two possible situations still remain. The concentration of b may be high enough, and declining slowly enough, to support the production of another threshold peak in X. If it is, then a new somite will mature, and with it another source of b, allowing the somite production to continue. The system recovers. On the other hand, if a threshold X peak is not produced, the b concentration will continue to decline below the Turing bifurcation level and cause the pattern in X and Y to vanish. The somite production collapses. Figure 6.2 shows a time-sequence of solutions output by br.c. The parameters used are those from table 6.2, except bs = 0.72. The vertical line indicates the location of the diffusion barrier (not inserted in the first graph). Graph b) shows the state of the system the instant the barrier was set in place. Graph c) shows the left-most peak in X approaching Tx, while the level of b is slowly dropping. Finally, the system recovers and graph d) shows somite production resuming. The time evolution in the region to the left of the diffusion barrier is not represented. Instead, the solution is frozen the instant the barrier is inserted. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 61 a) 0 0.5 1 1.5 2 2.5 c) Cone. 0 0.5 1 1.5 2 2.5 d) Cone. 0.5 1 1.5 2 2.5 s Figure 6.2: Recovery of the cohesive zone propagation. The vertical line represents the position of the diffusion barrier, which was inserted at the time of graph b). The time intervals are a) t = 29.4, b) t = 30.6, c) t = 60.0, and d) t — 75.0. Despite being cut off from the source, the propagation of the cohesive zone resumes. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 62 Figure 6.3 shows another time-sequence of solutions output by br.c. The parameters are the same as those used to generate figure 6.2, except bs = 0.71 and Df, = 0.005. Instead of a full recovery and continuation of the cohesive zone propagation, the system collapses as the 6-distribution decays. Graph c) shows the X patterning almost reaching the threshold. As the concentration of b decreases towards the bifurcation value (b = 0.4198851), the amplitude of the corresponding pattern in X and Y gets smaller. In graph d), morphogens X and Y show very little patterning at all. The ultimate state of the system in this setting is homogeneity. Both of these cases are possible with our model. A slight change in the parameters that define the dynamics of b can change the eventual outcome of the system. Also, the location and timing of the diffusion barrier insertion plays a role in the recovery or collapse of the cohesive zone propagation. According to our model, the further from the source the diffusion barrier is placed, the less likely it is that the system will recover. This hypothesis is experimentally testable. Chapter 6. Dynamic Production of Somites: Cohesive Zone Propagation 63 a) 0 0.5 1 1.5 2 2.5 c) I —j i i i i 1 Cone. 0 0.5 1 1.5 2 2.5 d) I i i i i i I Cone. 0 0.5 1 1.5 2 2.5 s Figure 6.3: Collapse of the cohesive zone propagation. The vertical line represents the position of the diffusion barrier, which was inserted at the time of graph b). The time intervals are a) t = 29.4, b) t = 30.6, c) t = 75.6, and d) t = 300.0. The peak in X does not quite reach the threshold, causing the b concentration to collapse. Cohesive zone propagation terminates. Chapter 7 Conclusions This multi-faceted study of the Brusselator chemical system to model somite formation in vertebrate embryos yields a number of important results. The solutions corresponding to parameter values from region £ of the Lacalli-Harrison parameter space tend to converge to their equilibria more quickly than those corresponding to parameter values from region Q. In region £ , the interval over which the dispersion curves is positive is bounded away from zero. This seems to allow for more direct convergence to steady-state. A bifurcation analysis yields that the Turing bifurcation is a pitchfork-type bifurcation. The global bifurcation analysis with a finite Fourier approximation to the Brusselator uncov-ered multiple branches of stable, non-homogeneous equilibrium solutions. It was first thought that these solutions might be an artifact of the finite nature of the solution representation. However, evidence for the existence of these stable solutions in the exact system is gained by numerical simulation. Although these simulations are also finite approximations to the Brusse-lator, confidence is derived from the fact that the method of discretization is different: spatial finite-differencing rather than finite Fourier series. The gradient in source chemical B may play a important role in the outcome of the system. Even a sub-critical gradient in the distribution of B can play a part in determining to which equilibrium the system converges. The motion of the cohesive zone was modelled using feedback between the X morphogen and chemical B. The patterning of X causes the location of the B source to move, while the distribution of B affects the X solution. The simulations that were performed on the system 64 Chapter 7. Conclusions 65 indicate that this type of interplay could possibly explain the sequential production of somites and propagation of the cohesive zone. The behaviour of our model in the cohesive zone seems to settle into a cycle. A newly-matured somite moves the cohesive zone, which in turn causes the maturing of another somite, thus continuing the cyclic behaviour. Moreover, the mere fact that the process of somite production and cohesive zone propagation is observed suggests that this near-periodic behaviour is stable. The introduction of a diffusion barrier into the cohesive zone gives two possible outcomes: somite production collapse, and somite production recovery. Experiments involving diffusion barriers may be conducted in search of this dichotomous phenomenon, the observations of which may give evidence for or against our proposed model for somite production and cohesive zone propagation. B i b l i o g r a p h y [1] Armstrong, J .B. , and Graveson, A . C . (1988). Progressive patterning precedes somite seg-mentation in the Mexican Axolotl (Ambystoma mexicanum). Developmental Biology, 126. [2] Armstrong, J .B. , and Malacinski, G . M . (1989). Developmental Biology of the Axolotl. New York: Oxford University Press. [3] Ascher, U . M . , Ruuth, S.J., and Wetton, B .T .R . (1995). Implicit-explicit methods for time-dependent PDE's. SIAM J. Numer. Anal. 32, no. 3. [4] Burden, R .L . , Faires, J.D., and Reynolds, A . C . (1981). Numerical Analysis. Boston: Prindle, Weber, Schmidt. [5] Crank, J . (1975). The Mathematics of Diffusion. Oxford: Clarendon Press. [6] Crank, J . , and Nicolson, P. (1947). A practical method for numerical evaluation of solutions of partial differential equations of heat-conduction type. Proc. Camb. Phil. Soc. 43, 50. [7] Doedel, E . (1986). AUTO: Software for continuation and bifurcation problems in ordinary differential equations. Pasadena: California Institute of Technology. [8] Edelstein-Keshet, L . (1988). Mathematical Models in Biology. New York: Random House. [9] Harrison, L . G . (1993). Kinetic Theory of Living Pattern. New York and Cambridge: Cam-bridge University Press. [10] Herschkowitz-Kaufman, M . (1975). Bifurcation analysis of nonlinear reaction-diffusion equations. II: Steady-state solutions and comparison with numerical simulation. Bull. Math. Biol. 37, 589-636. [11] Holloway, David M . (1995). Reaction-Diffusion Theory of Localized Structures with Appli-cation to Vertebrate Organogenesis. Ph.D. thesis, University of British Columbia, Canada. [12] Holloway, D . M . , Harrison, L . G . , and Armstrong, J .B. (1994). Computations of Post-Inductive Dynamics in Axolotl Heart Formation. Developmental Dynamics 200, 242-256. [13] Lacalli, T .C . , and Harrison, L . G . (1978). Turing's conditions and the analysis of mopho-genetic models. J. Theor. Biol. 76. [14] Lacalli, T .C . , Wilkinson, D .A. , and Harrison, L . G . (1988). Theoretical aspects of stripe formation in relation to Drosophila segmentation. Development 104, 105-13. [15] Meinhardt, H . (1982). Models of Biological Pattern Formation. London: Academic Press. 66 Bibliography 67 [16] Moore, K . L . (1974). The Developing Human. Philadelphia: W . B . Saunders Company. [17] Murray, J .D. (1989). Mathematical Biology. Berlin: Springer-Verlag. [18] Nicolis, G. , and Progogine, I. (1977). Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations. New York: Wiley. [19] Press, W . H . , Flannery, B.P., Teukolsky, S.A., and Vetterling, W . T . (1990). Numerical Recipes in C. New York, Cambridge University Press. [20] Prigogine, I., and Lefever, R. (1968). Symmetry-breaking instabilities in dissipative sys-tems. II. J. Chem. Phys. 48. 1695-700. [21] Ruuth, S.J., (1995). Implicit-explicit methods for reaction-diffusion problems in pattern formation. J. Math. Biol. 34, 148-176. [22] Stakgold, Ivar (1979). Green's Functions and Boundary Value Problems. New York, John Wiley & Sons. [23] Turing, A . M . (1952). The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B237. A p p e n d i x A Algebra i c Deta i l s Details of derivation of equation 3.7. L{(3)u = L(p) = L(0) = H0) = L(P) = HP) £L 1 2 (2) 2 (1) 2 u ( 2 ) 2 2 (2) u „ ' 2 2 ..(2) + HP) T ^ O O (1) Eoo (2) n = 1 ?in COS n7TS + a -P V^oo (1) Eoo (2) n = 1 t i ^ cos nirs 71 = 1 OO 71 = 1 (/3 — 1 — n2-ir2Dx) Un^ cosn7rs + c r a^ cosmrs —Pu^ cos n 7 r s + (—a — n2DyiT2) un2^ cos mrs P-l- n2n2 A a -a — n2TT2Dy ( i ) •u„ cos n7TS (2) t i n cos mrs oo = Lo(P)u0 + ]P Ln(P)un 71 = 1 Details of derivation of equation 3.11. m m - - , 0 r i ( 0 - l ) u V + D ^ + au(2) (L(P)u,v) - o / i = 2 = 2 Jo ds (/?-l)^) W 2 > + D , ^ ) W 2 > lTi(Dv(l) + au(2)vW _ puWv(2) _ a u ( 2 ) w ( 2 ) ) d s + 2 / D xv ( 1 ) 0 ' ds2  xv^> + Dvv™ ds Jo ds ds2 ds2 68 Appendix A. Algebraic Details 69 Let I = Jo ((ft - l ^ M 1 ) + OLU^VW - Pu^vW - auWyW) ds. Intergrating the second in-tegral by parts and using the boundary conditions for u, we get duWdvW „ duWdvW (L(P)u,v) r fr, ( D d u ( 1 ) r, '^du^Y' I + [ D ' v i ) - a r + D^y-BT) o dv<® Jo = j\^{{P-l)v^ ~Pv^)+u^[ ds ds ds ds -I Dx— — -r-Dy^—ds + D u^-QS2 + U V U Q S 2 m/1) — otv ds aV 2 > ds2 = / Jo = (n,L*(P)v) (P-1) + Dx£ a -a + Dy£ ds Details of derivation of equation 3.24 1 PMK +7 (x$W + h$x2 + ...) (W2> + 4 2 0 ) x 2 + ) +CT (WD + 4 1 0 )o; 2 + . . ( W 2 ) + fc^x2 + ...)) (x2 ( / 3 M ^ ( 1 ) 2 + 7 * ( 1 ) $ ( 2 ) ) +x 3 ( 2 / ? M ^ ( 1 ) 4 1 o ) + 7 $ ( 1 ) 4 o ) + 7 $ ( 2 ) 4 o + V2>) + O ( N 4 ) ) 1 ( x 2 i V 2 + x 3iV3 + o ( |x | 4 ) ) Details of derivation of equation 3.26. Jo $*N2 ds = f1 Jo - f Jo $* ( / 3 m K $ ( 1 ) 2 +7$(1)$(2)) cis = (2) cos MTTS (PMK<P^2 + 7<^1V^) C O S 2 -^ 7 R S ^S Appendix A. Algebraic Details <P' ,(i) .(2) ( & K / ) 2 + 7</?(1V2)) J cos3 MTTS ds Details of derivation of equation 3.27. 2x: Jo $*N3 ds = 2x3 £ (2/3M«*(1)4o) + 7*(1)4o) + 7*(2)4o) + <r*(1)V>) (V ^ _ r (2)) rfs = 2 x 3 ^ ( 2 / ? M ^ ( 1 ) 4 O ) ( « ) + 7¥'(1)42O)(5) + 7^ ( 2 )4 1O )(S) + ^ ( 1 ) V ( 2 ) COS 2 M7TS) (V(1) - </(2)) cos2 MTTS ds = 2xA 2 2 2 2 2 2 N(2) V = X 3 (A4d,0 + -^ 4 0,0 + c4o,2M + ^4o,2M + Details of derivation of equation 3.33. LM^-20 + N2 LM(0M)h2O + ( / ? M ^ ( 1 ) 2 + 7$ ( 1 )$ { 2 )) (i) = 0 = 0 W2 0 + ( / W ( 1 ) 2 + 7</?(V2)) ( / W 1 ) 2 + 7¥> (V 2 )) cos2 MTTS = 0 (1 + COS2MTTS) = 0 PM - 1 + Dx j^z -PM a -a + D. d2 Vds^ l20 l2 0 = 4 ( / W 1 ) 2 + 7 ¥ > ( 1 V 2 ) ) Details of derivation of equation 3.34. (1 + cos 2Mns) Appendix A. Algebraic Details 71 / | M27r2DY M2TT2DX + I\ 1 M2ir2DY a M2TT2D7, a 102TT20.02 102vr20.002 + l \ 1 102TT20.02 = 1 + 36 102TT20.02 36 = 0.392309 Details of derivation of equation 4.3. 1 0 d x T t ' dXn dXi 1 dt + dt COS TTS , dYq dYi 0 dt dt dXn dt 1 ds + 2 [ dXi dt 1 dY„ dt 0 Jo dY: dt 0 fl 1 2/„ 2 /' d X ° d: | ; [' JO dt Jo cos 7TS ds ds + 2 —-— cos 7rs ds lo dt Jo dt dX0 dt Details of derivation of equation 4.4. 1 d x T t ' 0 COS TTS ) = dXn dX-i 1 dt +" dt COS TTS , dY„ dYi 0 dt dt COS TTS Jo 2 dXp dt dYn dt COS TTS ds + 2 Jo dXi dt dYt dt cos2 TTS ds fl dXQ , n f1 dXi 2 = / —— cos TTS ds + 2 / —— cos TTS ds Jo dt Jo dt dXi ~dT ' Details of derivation of equation 4.5. 1 0 = 2 (a - (b + d) (^j- + Xi cosTTS^ J + ' X \ 2 ( Yo C | + X\ COS TTSJ ( — )- Y\ COS TTS ) — TT2DXX\ COS TTS ] ds Appendix A. Algebraic Details 72 = 2a - (b + d)X0 + 2c f1 ds Jo 4 2 + 2 ^ * (c^j-Yi + c X 0 X x y - 7 r 2 D x X i j cos TTS ds + 2 c ^ ^ X o X x F i + I 2 y j cos2 7rs ds + 2c £ XfYi cos 3 TTS ds = 2a - (6 + d ) X 0 + c Qx 2 Yo + XQXIY! + \x2Y^j 


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