Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Reaction-diffusion theory of localized structures with application to vertebrate organogenesis Holloway, David M. 1995

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

Item Metadata


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

Full Text

REACTION-DIFFUSION THEORY OF LOCALIZED STRUCTURES WITH APPLICATION TO VERTEBRATE ORGANOGENESIS By David M. Holloway B.A. (Chemistry) The University of Puget Sound, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF CHEMISTRY We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 1995 © David M. Holloway, 1995 In presenting this thesis in partial fulfillment 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. (Signature) Department of CflgfrliVry The University of British Columbia Vancouver, Canada Date Ar:i u. im. Abstract Biological pattern formation is addressed within the framework of Turing's (1952) chemical theory. His (linear) reaction-diffusion equations for chemical concentration have spatially inhomogeneous solutions. These stationary chemical waves are a possible mechanism for the spatial patterning of cellular differentiation in embryogenesis. Chemical mechanisms which give Turing dynamics have nonlinear rate equations. Much of the dynamics can be predicted from linear analysis, but solution is generally only possible numerically. Numerical methods for reaction-diffusion equations, especially implicit techniques, are reviewed. This work is specifically concerned with the localization of pattern. Reaction-diffusion models tend to form repeated peaks. Under some circumstances, however, patterns of localized peaks surrounded by large areas of low concentration can be formed. The requirements for this are developed and explored through comparison of two models, the Brusselator and the Gierer-Meinhardt. It is found that both can give localized pattern if the ratio of the component diffusivities is high. This leads to partly-ordered pattern, linking Turing theory with the inhibitory field theory of Wigglesworth (1940). However, the Gierer-Meinhardt gives less ordered pattern in two dimensions and is more susceptible to gradients in precursor concentration than the Brusselator. The different dynamics in the two models, underlying the different behaviours, are discussed. ii These principles of pattern localization in reaction-diffusion models are applied to heart formation in the axolotl, Ambystoma mexicanum. A form of the Gierer-Meinhardt model is used (but derived from a chemical mechanism), due to its greater tendency to form localized pattern. Heart formation proceeds by induction of the mesoderm tissue by the underlying endoderm, followed by more localized cardiomyocyte differentiation in the mesoderm. This sequence is modelled by treating induction as a process which establishes the distribution of a precursor in the reaction-diffusion mechanism postulated to be responsible for cardiomyocyte differentiation. Normal development, and numerous transplant and explant experiments with wild-type and cardiac-lethal mutant (c) tissues are successfully modelled. Agreement with experiment is achieved in terms of spatial patterns and timing of events. Predictions for future experiments are given, as well as some suggestions on the probable chemical nature of the morphogens, with some relationship to signal transduction mechanisms. iii Table of Contents Abstract ii List of Tables viii List of Figures ix Acknowledgement xii 1 Introduction 1 Part One: Theoretical Aspects 9 2 Turing Theory and Thermodynamics 10 2.1 Kinetic Theory of Chemical Pattern Formation 10 2.2 Thermodynamics of Chemical Pattern Formation 20 3 Numerical Methods 24 3.1 Introduction 24 3.2 Numerical Solution of Reaction-Diffusion Equations 25 3.2.1 Explicit Methods 27 3.2.2 Implicit Methods 32 Fully Implicit Method 33 Crank-Nicolson Method 34 Peaceman-Rachford Method 35 iv 3.3 Random Input 42 3.4 Computing Facilities 47 4 Localization of Reaction-Diffusion Pattern 49 4.1 Introduction 49 4.2 Diffusion Terms 51 4.2.1 Evaluation of Pattern Order and Peak Sharpness 53 4.2.2 Models and Methods 55 4.2.3 Results 58 4.2.4 Discussion 69 4.3 Precursor Gradients 72 4.3.1 Methods 73 4.3.2 Results 75 Self-effects 75 Cross-effects 78 Exponential gradients 81 Gradient growth 82 4.4 Reaction Terms 83 4.5 Conclusion 90 5 A Chemical Basis For the Gierer-Meinhardt Model 94 5.1 Babloyantz-Hiernaux Mechanism 94 5.2 The Enzyme Kinetic Scheme 98 Part Two: Biological Application 102 6 Experimental Work on Axolotl Heart Formation 104 6.1 Introduction 104 6.2 Induction 104 6.3 The cardiac lethal mutation in the axolotl 112 6.4 Experimental characterization of heart-forming chemicals 118 6.5 Conclusion 121 7 Modelling Induction 124 7.1 Introduction 124 7.2 The inductive process 124 7.3 Inducer distribution 129 7.4 Saturation 139 7.5 Region of specification 140 7.6 The product of induction 141 7.7 Conclusion 144 8 Modelling Heart Formation 147 8.1 Introduction 147 8.2 Chemical Mechanism 147 8.3 Results 150 8.3.1 One-Dimensional Computations 151 vi Wild-Type Heart Development 151 Stage 20 Transplants 154 Stage 20 Explants 158 Stage 29 Transplants 160 Early in vivo studies 162 Exogenous stage 35 rescues 166 8.3.2 Two-Dimensional Computations 168 8.4 Discussion 177 8.4.1 Predictions for experiment 178 8.4.2 Molecular possibilities 180 8.4.3 Reaction-diffusion vs. mechanochemistry 182 8.4.4 Reaction-diffusion vs. inhibitory field theory 183 8.5 Conclusion 184 9 Conclusion 185 Bibliography 191 vii List of Tables 4.1 Brusselator parameters 57 4.2 Gierer-Meinhardt parameters 58 4.3 Clark and Evans R parameters 67 4.4 Average distance to 1/4 X peak height 69 7.1 Time to beat for in vitro induction 127 7.2 Position of leading edge of mesoderm 130 7.3 Extent of specified region (SR) in dorso-ventral direction 141 8.1 Parameter values for equations (8.1) 151 8.2 Time to beat (in days) for stage 20 transplants 156 8.3 Time to beat (in days) for stage 20 explants 158 viii List of Figures 2.1 Chemical pattern formation by reaction-diffusion 2.2 Dispersion plot in region D of Lacalli parameter space 2.3 Lacalli parameter space 2.4 Dispersion diagrams in different regions of Lacalli parameter space 3.1 The numerical instability in the explicit technique 3.2 Comparison of computations with different time step 3.3 Comparison of computations with different spatial step 3.4 Further comparison of the effects of the time step 4.1 Position in Lacalli space of the four parameter sets used 4.2 Two-dimensional Brusselator patterns with increasing n 4.3 Radial distribution functions for the Brusselator patterns 4.4 Two-dimensional Gierer-Meinhardt patterns with increasing n 4.5 Radial distribution functions for the Gierer-Meinhardt patterns 4.6 Dispersion plots for the Brusselator and Gierer-Meinhardt with increasing n 4.7 Clark and Evans R parameters for Brusselator and Gierer-Meinhardt patterns 4.8 Gradient effects in the linear decay terms of the Brusselator and Gierer-Meinhardt 4.9 c gradient effects in the Gierer-Meinhardt 4.10 Growth and halt of pattern in the Brusselator and Gierer-Meinhardt 5.1 Comparison of the Babloyantz-Hiernaux mechanism with the Gierer-Meinhardt 6.1 Cross-section through the axolotl embryo at the mid-thoracic level 6.2 Mesoderm advance in neurula stage axolotl embryo ix 6.3 Side views of the post-neurula axolotl embryo, showing specified region 6.4 Cross-sections oiAmbystoma punctatum embryo, showing Copenhaver's surgery 7.1 Side view of stage 14 axolotl, showing different regions of endoderm 7.2 Cross-section through stage 13 axolotl embryo 7.3 Stage 14 mesoderm explants with corresponding computations 7.4 Modelled concentration distribution of the endodermal inducer 7.5 Angular coordinates in the endoderm 7.6 Model for mesoderm movement over the endoderm with fixed dorsal point 7.7 Model for mesoderm movement over the endoderm with fixed dorsal region 7.8 Computed concentration of the product of induction 7.9 Model for the effect of mesoderm explantation on concentration of the inducer product 8.1 Computation of normal heart development 8.2 Computation of stage 20 bilateral mesoderm transplant with c/c host 8.3 Computation of stage 20 unilateral mesoderm transplants 8.4 Computation of stage 20 mesoderm explants 8.5 Computation of stage 29 mesoderm transplants 8.6 Computation of cardia bifida 8.7 Computation of Copenhaver (1926) experiments 8.8 Computation of stage 35 rescue of c/c heart 8.9 Orientation of two-dimensional computation on the axolotl 8.10 Two-dimensional computation of stage 20 unilateral transplant 8.11 Rescue by R in two dimensions 8.12 Two-dimensional computation of normal heart development x 8.13 Two-dimensional computation of stage 20 mesoderm explants XI Acknowledgment I would foremost like to thank Prof. Lionel Harrison for being both a supportive and a critical supervisor. His clear vision of the interplay of experiment and theory has inspired my work, as has his ability to evaluate the logical underpinnings of scientific works. I would like to thank Dr. Michael Lyons, who passed on much of his graduate wisdom and who was instrumental in getting me off the ground in this research. Thanks to Jacques Dumais for the many philosophical diversions we have shared, which have helped to develop my perspective. Drs. John Armstrong, Steve Smith, and Heather Easton have been very generous with their axolotl data and their thoughts on the heart model. The close collaboration with John Armstrong has led to a very detailed model, and has allowed me more insight into experiment than is usual in a theoretical project. I also thank Helen Salter, Paul Hilchey and Susan Mair of the University Computing Services, all of who have helped to guide me through the computer maze. I would also like to acknowledge the contributions my family have made, most especially by their examples for pursuing one's passions. In particular, I would like remember my grandmother, Maria Neill. And final thanks to Cynthia, with who I am honoured to share my life. xii Chapter One Introduction The origin of form has long been a central question in human thought. Early Greek philosophy, arising in the 6th century B.C., dealt with explanations of cosmological order. By the early 5th century B.C. thinkers such as Empedocles had begun to address the questions of structure and order in biology (Kirk et al., 1983). The events of biological development have continued to be a chief source of speculation concerning form. The basic question of how the egg develops into an adult organism is far from being answered. The myriad levels of organized pattern-forming events constituting this progression are only dimly perceived. Though there are many philosophical aspects to the questions of development, the approach in this thesis will be scientific. In particular, the theory pioneered by Turing (1952) is the major topic of this thesis. A scientific attempt to answer some of the questions of biological development rests on two pillars: experiment and theory. A detailed experimental description of the phenomena must exist before a theory can be attempted. The difficulty of obtaining data in biology (especially of the many well-hidden events in embryology) requires great skill and time of the experimental biologist. This alone creates a tendency towards experimentalism in biology. As well, the plethora and complexity of the data obtained resists theoretical generalizations. For these reasons, biology tends to be a field characterized by experimental description rather than theoretical unification. Physical scientists have a different perspective: the precision of data lends itself to generalization into theory. In general, one can observe simple physical systems obeying deterministic laws. The dynamics of deterministic systems are described by differential Chapter 1. Introduction 2 equations. The rich theory behind differential equations allows for connections to be made between disparate experiments, as well as providing strong predictions for future experiment. For this reason, physical sciences tend to have a larger proportion of theory to experiment than is the case in biology. The differential equations which underlie Newtonian physics are conservative and describe systems with orbital (non-asymptotic) stability. If the parameters of the equation are changed slightly, the whole system moves to a new trajectory. Perturbations of the system remain and do not die out. Description of biological systems must incorporate different stability characteristics. Biological systems are stable in a fluctuating environment, and can evolve. Differential equations that incorporate dissipation can have asymptotic stability which insures movement towards a specific steady state despite fluctuations in the environment. Changes in system parameters lead to changes in the asymptotically stable steady state, not to simple shifts in trajectory. The current state of experimental developmental biology is quite advanced. Morphological characterization of the major events in embryonic development was completed by the mid-twentieth century. With the advent of modern genetics and molecular biology in the 1960's, the tools have been available to determine the molecular species responsible for the morphological events. Hierarchies of genes and the interactions of their products have now been characterized for many of the major events in development. The picture which is developing is incredibly detailed and complex. In the past decade experimental methods have allowed for the in vivo visualization of a few of the spatial chemical patterns which precede tissue differentiation1. At this level, experimental biology is coming to the point of describing the chemical events underlying morphology. However, the biochemical approach lrrhe existence of these chemical patterns was not an expectation of most developmental biologists, but had been predicted much earlier by Turing. Chapter 1. Introduction 3 has done little to explain chemical pattern formation. At present, unification of disparate experimental results lies in the identification of recurrent motifs in molecular structure. A physicochemical theory of pattern formation was first set forth by Turing (1952). The Turing theory consists of a system of differential equations incorporating the processes of reaction and diffusion in a chemical system. These equations describe the growth of chemical concentration pattern from a homogeneous steady state. The specific (nonlinear) chemical models which are described by Turing dynamics generate concentration patterns which have the asymptotic stability necessary for biological description. Although the equations have a wider applicability, they were applied to the problem of biological morphogenesis in Turing's original paper. This application represents a bridge between the rich fields of physical theory and experimental biology. With the current state of genetics and molecular biology it is becoming possible to discuss the reaction-diffusion dynamics of chemical pattern formation with respect to specific molecules, or classes of molecules. Given the complexity of biochemistry, however, there is still a great deal of work needed to link the kinetics of biological pattern formation to particular molecular species. The simplicity of the Turing equations serves as a starting point for a conceptual understanding of these processes. In the present work, I apply the theory of chemical pattern formation to the problem of vertebrate organogenesis. A primary question in organogenesis is how a single organ differentiates out of a larger area of identical tissue. Why do embryos form single organs such as the heart instead of multiple repeats? As first formulated, Turing theory describes patterns of repeated parts (the patterns are basically sinusoidal). Some additions must be made to the theory to describe localized patterns, such as those typified by vertebrate organogenesis. Chapter Four is a discussion of the properties of reaction-diffusion systems which can be used to describe localized pattern. It is found that pre-existing gradients can localize Turing pattern. If a concentration gradient is present in the system, the Turing Chapter 1. Introduction 4 dynamics tend to amplify it. At the relatively late developmental stage of organ formation numerous gradients are present in the embryo. Organogenesis occurs well after the polarities of the embryo have been established, when the major axes are well defined by chemical gradients (Ding and Lipshitz, 1993). This work is divided into two Parts. In the first, I develop the theoretical aspects necessary for application of reaction-diffusion theory to vertebrate organogenesis. Chapter Two is a short explanation of the Turing theory and an introduction to the effects of nonlinearities. As well, the thermodynamic approach to far from equilibrium systems is discussed. Chapter Three details the numerical methods used in solution of the reaction-diffusion equations. Chapter Four examines the properties of reaction-diffusion equations which can be used to produce localized pattern. Chapter Five describes a chemical mechanism devised to give reaction-diffusion equations that tend to produce localized solutions. In the second part of the thesis, I apply the theory developed to vertebrate organogenesis. In particular, I describe a reaction-diffusion model of early stages of heart formation in the Mexican salamander Ambystoma mexicanum, or axolotl. In this work we have collaborated with experimental biologists at the University of Ottawa. Because of the detailed experimental data we have received, the model has been developed into a quantitatively detailed description of the heart forming events in the axolotl. Chapter Six is a review of the experimental work to date on heart formation in the axolotl, with some comparison to other vertebrates. Chapter Seven is a discussion of how the dynamics of inter-tissue signalling (induction), which underlie heart differentiation, are incorporated into the model. Chapter Eight is a description of the reaction-diffusion system used to describe heart differentiation, the results of the modelling and the conclusions that can be made from these results. Chapter 1. Introduction 5 The theoretical basis for this thesis is the chemical theory of pattern formation as developed by A. M. Turing (1952). This was the first complete formalism for a kinetic theory of pattern formation. It embodies dynamic principles which are more general than Turing's specific use of them in reaction-diffusion, and which have later appeared in other kinds of kinetic theory (Harrison 1993, Ch. 8). Turing's system consists of two mutually interacting chemicals which diffuse. These dynamics give rise to stable spatial concentration patterns in the two chemicals. In his original paper, Turing developed the mathematical formalism and then applied it to problems in biological morphogenesis. His theory remained in relative obscurity for some years, but has achieved a wider readership more recently (Harrison, 1994). His theory has been elaborated with specific nonlinear models which simplify to the general linear Turing equations. Most notable of these contributions are the Brusselator model of Prigogine (1967a; also Prigogine and Lefever, 1968) and the model of Gierer and Meinhardt (1972). Although the processes of reaction and diffusion had been coupled in previous work (Rashevsky, 1940), Turing's paper offers a complete formalism. He explicitly discussed the local interactions which could give stable macroscopic pattern. As well, Turing's realization that these processes needed to only amplify existing modes in the thermal noise was original. For these reasons, Turing's work stands as a landmark in the theory of pattern formation. Turing's theory involves the dynamics of two chemical intermediates, X and Y. These intermediates (morphogens) require sources maintaining constant concentrations of precursors (e.g., A and B in the Brusselator, see Chapter Four), hence the system is open. There are several important aspects to the precursors, or reactants, in any chemical mechanism which has Turing dynamics. The first is that they provide a communication via material flow between the pattern forming system and the rest of the universe. This communication allows the entropy reducing events of pattern formation to occur without Chapter 1. Introduction 6 violation of the second law of thermodynamics (see section 2.2). Second, the precursors may themselves be products of previous developmental events. In modelling embryology, then, it is possible to account for a whole hierarchy of developmental processes within the context of a single model. Finally, the earlier events may lead to spatial patterning of the precursors, for instance in simple monotonic gradients. These additional spatial dependencies change the dynamics enough to make analysis far more difficult. However, the altered dynamics can lead to new types of pattern for the intermediates. Hunding (1987; 1989) has described reaction-diffusion mechanisms with pre-existing gradients as "Turing systems of the second kind". Although mathematical analysis is difficult and of limited scope (presently), these patterns can be studied through numerical analysis. Significant parts of this thesis (section 4.3, Chapter Seven) are concerned with the generation of localized structures by control of reaction-diffusion dynamics through precursor gradients. An understanding of localized chemical pattern is necessary for the current application to vertebrate organogenesis. Kinetic models which give the Turing dynamics, such as the Brusselator or Gierer-Meinhardt, involve nonlinear chemical reactions. Specifically, self-enhancement of the activator requires a reaction order greater than unity. As well, other reactions in the scheme can be of order greater than unity. Linearizing these models about the homogeneous steady state results in linear and higher order terms. The linear dynamics in the models can provide accurate information about the nature and growth rate of the pattern in its early stages. There are, however, effects from the nonlinear terms, without which the Turing dynamics would be of limited scope in trying to match theory and experimental details precisely. Nonlinearities provide feedback in the model which allows for arrest of pattern growth. The linear dynamics lead to exponential growth of modes, with no halt. This description is valid for initial growth of pattern near the homogeneous steady state, but breaks down away from the homogeneous state when the patterned state begins to dominate. In this regime, the Chapter 1. Introduction 7 pattern amplitudes are of such a scale that the nonlinear terms are of comparable magnitude to the linear terms. When this occurs, the nonlinearities serve to halt the growth of pattern. Nonlinear models also allow for mode competition. In the linear theory, the modes are independent, each undergoing exponential growth with a specific growth rate. However, with the addition of nonlinear dynamics, there is interaction between modes. The final pattern tends to be dominated by the fastest growing modes, but is a mixture of modes which linear analysis is unable to predict. For instance, in two-dimensional pattern there are cases in which the nonlinear interactions between modes can lead to striped pattern where linear analysis would predict spotted pattern. In these cases, nonlinear dynamics lead to instability (with respect to the patterned state) for all but one of the modes with positive linear growth rate (Lyons and Harrison, 1991; 1992). This is a good example of how different nonlinearities (stemming from different chemical mechanisms) can give different final pattern with properties that cannot be described with linear theory. Description of these nonlinear effects must be through numerical analysis or nonlinear analysis. In general, nonlinear analysis is difficult and of limited scope. There are still a number of questions that can be pursued only through numerical techniques. Chapter Four will be a numerical study of some differences between two reaction-diffusion models. Thermodynamic characterization of systems which are out of equilibrium has been developed to a certain extent over the past sixty years. Chemical pattern formation, and processes which derive order out of random initial conditions in general seem to contradict the basic tenets of the second law of thermodynamics. In fact, biological pattern formation drove Schrodinger to devise the concept of 'negative entropy' (Pauling, 1987). Such additions to thermodynamics are not necessary, however. The second law is usually stated in terms of the Clausius relation dS^dQ/T (1.1) Chapter 1. Introduction 8 where the inequality holds for irreversible processes. This shows that the entropy change is always positive for irreversible, adiabatic processes (or for free isothermal expansion of an ideal gas). Since dQ = 0 for any isolated system, this leads to entropy increase in the universe, as long as any irreversible processes are occurring. However, positive entropy change is not required for non-isolated systems. Therefore, the entropy can decrease for a given system (such as a Turing system), but always at the expense of increased entropy in the surroundings. The study of non-equilibrium thermodynamics attempts a description of the entropy flow for non-isolated systems and the derivation of thermodynamic stability criteria for non-equilibrium steady states. This will be discussed in more detail in section 2.2. It is in principle possible to formulate reaction-diffusion theory of pattern formation in terms of nonequilibrium thermodynamics. However, the thermodynamic approach tends to be more cumbersome than the kinetic formulation, as well as being more removed from the current investigations in biological development. The kinetic approach is more promising for eventual connection with biochemical and genetic approaches to development. The focus of this thesis will be kinetic. The emphasis in this work is on stability criteria and the spatial patterns which emerge far from equilibrium. The reaction-diffusion theory of Turing as developed in the next chapter provides simple stability analysis. As well, Turing theory provides information on the growth of instabilities and the wavelengths of the patterned instabilities. It is a simple and elegant approach to the questions addressed in this work. PART ONE THEORETICAL ASPECTS bink' epeco; TOTE u.ev yap ev T]ij§n0r] povov etvat 8K JtXeovcov, TOT8 5' av Siecjnj JtXeov' e | evas etvat. ...Tnt pev Ytyvovxat xe Kat ou otptatv epjteSos atcav; r|t 6e StaXXaaaovta Stapjr.epe.s oi)5apa Xr\yei, Tocurnt 5' oaev eaatv aKtvnTot Kara ICUKXOV. [I will tell a twofold tale; for at one time they [i.e. the roots] grew to be one alone out of the many, at yet another time they escaped to be many out of the one. .. .in this way, they both come into being and yet do not have an eternal lifetime; insofar as they never cease their continual interchange, thus far are they always permanent in the cycle.] -Empedocles, early 5th century B.C. (in Kirk et al. 1983) This part of the thesis will introduce the theory and methods used in characterizing a reaction-diffusion theory of local structures. Chapter Two will introduce the formalism of the Turing theory necessary for a discussion of the reaction-diffusion patterns in this thesis. Chapter Three will discuss the numerical methods used in solution of the specific nonlinear models that will be worked with. In Chapter Four, I will present results which help to define the characteristics of reaction-diffusion models which can be manipulated to create localized pattern. In Chapter Five a chemical mechanism will be developed which leads to a reaction-diffusion model that tends to form localized pattern. This model will be used in the subsequent application to vertebrate organogenesis described in Part Two. 9 Chapter Two Turing Theory and Thermodynamics In this chapter I present a short introduction to the formalism of Turing theory. Although summaries of the Turing theory are published elsewhere (Harrison, 1993; Lyons, 1991), it will be necessary to refer to specific details of the theory throughout this thesis. Therefore, I offer a short summary containing the essential aspects necessary for discussion of my results. As well, a section on the thermodynamics of far from equilibrium processes is included. Although the focus of this thesis is kinetic, it is important to recognize the information provided by the thermodynamic approach. Specifically, through a detailed entropy balance, the interactions with the surroundings necessary to maintain pattern formation in the system can be elucidated. 2.1 Kinetic Theory of Chemical Pattern Formation To begin the derivation of the Turing formalism, it is instructive to consider a one-morphogen system (Harrison 1993, pp. 184-190). If the morphogen X is autocatalytic, undergoes first-order decay and diffuses, then the equation governing its dynamics is as follows: dX/dt = -kdX + kfX2+ Dx d2X/ds2 (2.1) where s is the spatial variable, Dx is the diffusivity and kd and k, are rate constants. We are interested in the stability of the homogeneous steady state. That is, given a uniform concentration of X, and the assumption of a steady state (dX/dt = 0), can this state be maintained? The (non-trivial) X concentration at the steady state is X0-kd/kf (2.2) Chapter 2. Turing Theory and Thermodynamics 11 since the diffusion term drops out in a homogeneous system. If equation (2.1) is linearized about the steady state, one obtains BU/dt = kdU+Dxd2U/ds2 (2.3) where U = X-X0 is the deviation from the steady state. If one considers that U is homogeneous, that it has no spatial dependence, then equation (2.3) becomes dU/dt = kdU (2.4) whose solution is U = Cekit (2.5) where C is a constant, k^ is positive since it is a real rate constant (also, negative kj would imply negative XQ by equation 2.2, which is physically meaningless). Therefore, homogeneous perturbations grow exponentially. This is due to the autocatalysis in the mechanism. To study the effects of perturbations with spatial dependence, we assume that U is sinusoidal. This is consistent with a picture of U as a Fourier component of thermal fluctuations. For U = sin cos or U = cos cos, where co = 2nfk is the wavevector (co is used for the wavevector instead of the usual k, since the traditional chemical usage of ^s for rate constants is adopted in this work), d2U/ds2 - -co2£/ (2.6). Substituting into equation (2.3) gives dU/dt = (kd-<a2Dx)U. (2.7) The solution of this equation is U m Ceat (2.8) where o = kd- oi2Dx is the eigenvalue and C is a constant. The steady state is stable with respect to sinusoidal perturbations if a < 0 , unstable if a > 0 . a is negative at large wavevector (short wavelength) and patterned perturbations die out. However, at small wavevector (long wavelength), a is positive and sinusoidal pattern can grow. There is not a Chapter 2. Turing Theory and Thermodynamics 12 good discrimination of pattern, however, as the eigenvalue is strictly increasing for decreasing wavevector. The eigenvalue reaches a maximum at infinite wavelength (wavevector zero). The final pattern is therefore homogeneous, but growing. Thus, the single-component reaction-diffusion system is unable to amplify and maintain patterned perturbations. With the addition of a second chemical species, however, spatially dependent perturbations can have maximum growth rate. Such a system of two chemical intermediates is described by Turing dynamics. Turing's system has two chemicals, each of which has a self-interaction and a cross-effect with the other chemical. There are numerous mechanisms which can be devised to describe these interactions. Turing theory, however, assumes that the essential dynamics of the chemical interactions can be captured in the terms linear in the deviations from the homogeneous steady state. These four linear terms are mathematically the same for all two-morphogen mechanisms. It is important to make a strong distinction between the actual concentrations of the intermediates, denoted by X and Y, and their deviations from the homogeneous steady state, U = X-X0 and V = Y-Y0. The linear Turing equations are stated in terms of the deviations U and V. The linearized rate equations contain three types of kinetic terms. These kinetics are: 1) autocatalysis, or self-enhancement, for U; 2) cross-catalysis of Vby U; and 3) inhibition of U by V1. These interactions can come about through numerous actual reactions between X and Y. Therefore it is important to distinguish between the deviations and the true concentrations. For example, Meinhardt (1989) has shown that an inhibition of an inhibition leads to the autocatalysis in point 1) above. Simple bimolecular autocatalysis of Xis not necessary to produce Turing pattern. However, because bimolecular autocatalysis and direct inhibition are the simplest ways to think of the Turing dynamics, X is frequently called the activator and Y the inhibitor. When the three types of interactions are Pattern-forming dynamics also exist if effects 2) and 3) are both reversed. Chapter 2. Turing Theory and Thermodynamics 13 JJ=X-X0 ^ X . XQ, Y0,orU=V = 0 .. V= Y- Y, o Figure 2-1: Chemical pattern formation by reaction-diffusion. The horizontal line in each frame represents the homogeneous steady state concentrations XQ and YQ. Concentration is plotted in the vertical direction, the spatial dimension in the horizontal direction, a) This illustration starts with a small positive perturbation of X from its steady state value XQ (positive U). b) This is amplified through autocatalysis. The U peak cross-catalyzes production of V, and V diffuses faster than U (their diffusivities are required to be unequal in the Turing theory), c) Through inhibition of U, the spread of V causes regions of negative U laterally from the peak, d) The cross-effect of U on V generates negative V peaks where U is negative. This illustration suggests progress towards a patterned state. Mathematical analysis shows that this state exists, with spacing dependent on the rate and diffusivity parameters in the model. Turing models can also form pattern with out-of-phase U and V peaks if the catalysis and inhibition of the cross-effects are reversed. Redrawn (Harrison, 1993), following Maynard Smith (1968). Chapter 2. Turing Theory and Thermodynamics 14 coupled with Fickian diffusion (with DY * Dx), pattern formation can occur (Figure 2-1). This leads to the Turing equations: dU/dt m kfJ + k2V + Dx d2U/ds2 (2.9a) dV/dt = kJJ + k4V + DY d2V/ds2. (2.9b) If [/and Fare homogeneous perturbations, then equations (2.9) become dU/dt = k1U+k2V (2.10a) dV/dt = k3U + k4V. (2.10b) The eigenvalues of this system are ^ • ^ f ^ (2.11) where ± denotes the two eigenvalues, Tr = kx + k4 and det = kft4 - k2k3 (Lyons 1991, p. 17). For the system to be stable to homogeneous perturbations a* < 0, so Tr < 0 and det > 0. For the case of sinusoidal perturbations, the eigenvalues are still given by equation (2.11), but Tr = ky - (o2Dx + k4- oi2DY (2.12) and det = DxDy(a4 - {k^Dy + k4Dx)(a2 + (kxk4 - fcjfc,). (2.13) The Tr is always negative if homogeneous perturbations are assumed to be stable. However, the det can become negative for certain values of co2, allowing for instability with respect to sinusoidal perturbations. Equation (2.13) shows that the det is a quadratic function in co2. It is a parabola opening upwards (Edelstein-Keshet 1988, p. 513). The minimum of this parabola becomes negative for (k^Dy + kADxf - ADxDY(kJc4 -k2k3)>0. (2.14) Equation (2.14) gives the condition for the Turing bifurcation. If the inequality of equation (2.14) is met, then sinusoidal perturbations can grow. Figure 2-2 maps equation (2.11) for a particular set of wavevectors. This shows the selection of wavevectors afforded by Turing Chapter 2. Turing Theory and Thermodynamics 15 "9 0.004 0 0.004 Figure 2-2: Relation (2.11) is plotted versus wavelength (X = 2JI/G>), for the choice of parameters ^ = 0.03; k4 = -0.04; - k^ - 0.05; Dx - 0.02; DY - 0.2 (Harrison et al, 1984). Only a finite set of wavelengths have positive eigenvalue. The wavelength in this range, fitting boundary conditions and having maximum growth rate k (a+ in equation 2.15), will determine the spacing of the final pattern. This type of plot is also referred to as a dispersion diagram. From Harrison (1993). Chapter 2. Turing Theory and Thermodynamics 16 dynamics (its "band-pass" quality). Generally, the wavevector which fits the boundary conditions (Lacalli and Harrison, 1978a) and has the largest eigenvalue within this band-pass will grow to dominate the system. The eventual macroscopic pattern will have this wavevector. This wavevector can be found by calculating the maximum in the dispersion relation between the eigenvalue and the wavevector (equation 2.11). This gives iK-K)±{Dx+DY)}-k*yDD co2max = \ / x Y (2.15) where the positive sign is used for Dx>Dy and the negative sign is used for Dy>Dx- The wavelength of the pattern, the so-called "chemical wavelength", can be determined by substituting equation (2.15) into the relation comax = 2nfkm. The growth rate of this pattern is found by substituting comax into equation (2.11). A parameter space has been developed by Lacalli and Harrison (1979) which gives a clear understanding of the different qualitative behaviours of the Turing system (Figure 2-3). The six parameters in the Turing model, the rate constants, k\-k^, and the diffusivities, Dx and Dy, can be reduced to three. The self-effects, k\ and £4, can be divided by the cross effects, &2 and £3, to form the reduced rate parameters k^ = k^j^kjc^ and kA = kAj^-kJc^. The essential features of the diffusion terms are retained in their ratio, n = DY/DX. With these three reduced parameters, a two-dimensional k\ , A4 space is defined for a given value of n. Boundaries in this space indicate the conditions placed upon pattern formation in the Turing theory. Derivations of these conditions can be found in Harrison (1993, pp. 236-243). Line 3 is the Turing bifurcation defined by equation (2.14). The slope of this line is equal to n, the diffusivity ratio. Points in the space to the right of this line give positive pattern growth. Points above line b will have a maximum eigenvalue at a Xm less than infinity. Line c defines the boundary between solutions of the Turing equations with real and complex eigenvalues. Points below line c give monotonic sinusoidal pattern growth, while points above the line give sinusoidal pattern with oscillating amplitude. Above line e, Xm is Chapter 2. Turing Theory and Thermodynamics Figure 2-3: The parameter space of Lacalli and Harrison (1979). The dynamics of the Turing system are completely determined by the reduced parameters k\, £4, and n. For a given n, the dynamics can be plotted in this two-dimensional space. The lines, denoted by lower case letters, which divide the space are defined in the text. The regions, denoted by upper case letters, are discussed in Figure 2-4. From Harrison (1993). Chapter 2. Turing Theory and Thermodynamics 18 complex and oscillations grow faster than non-oscillating pattern. Line d corresponds to the Hopf bifurcation in homogeneous systems. It is the boundary between solutions of the Turing equations with complex eigenvalues that have positive real part and negative real part. Finally, the parabola f defines systems for which the eigenvalues are positive for a finite range of wavelengths. These boundaries define four main regions of the parameter space. These regions are bounded by lines a, b and e, outside of which stable pattern formation is impossible. In region A, there are complex eigenvalues with positive real part. However, Xm corresponds to a real eigenvalue. Region A is therefore of possible morphogenetic significance, as non-oscillating pattern will grow fastest (Figure 2-4c). In region B there are no complex eigenvalues, but there are positive eigenvalues at infinite wavelength (Figure 2-4a). In this region, therefore, there is not good selection of pattern. The selection can be adequate for simple patterns (Lacalli and Harrison, 1978a), but ordered patterns of many repeated parts do not arise for these parameters. Regions C and D display the best mode selection, as the eigenvalues for these parameters are positive only along a finite interval of wavelengths. Region C has complex eigenvalues, but they have negative real part (Figure 2-4d). Therefore, parameters in this region give damped oscillations, but are dominated by stable pattern. Region D has only real eigenvalues (Figure 2-4e), and therefore stable pattern grows monotonically. All modelling in this thesis is done with parameters in regions C and D. Figure 2-4: Dispersion diagrams for different regions of the parameter space, c) corresponds to region A, a) corresponds to region B, d) corresponds to region C, and e) corresponds to region D. b) is to the left of the Turing bifurcation (line a in Figure 2-3). The pattern dynamics in each region are discussed in the text. From Harrison (1993). Chapter 2. Turing Theory and Thermodynamics Figure 2-4 Chapter 2. Turing Theory and Thermodynamics 20 2.2 Thermodynamics of Chemical Pattern Formation Thermodynamics is most often used to provide information about state variables at equilibrium. However, work has been conducted since the late nineteenth century to try to characterize the value of state functions away from equilibrium. Two regimes can be defined for non-equilibrium thermodynamics. The first of these is linear non-equilibrium thermodynamics. This theory has been well developed since the time of Onsager (1931). The second regime is that of far-from-equilibrium, or nonlinear, thermodynamics. Though much work has been done in this area over the past forty years, there is still much controversy as to the approach necessary to define thermodynamics for these cases. I will discuss the general theory for the linear and nonlinear cases and the information which can be gained from thermodynamics concerning far from equilibrium pattern formation. The thermodynamics of processes close to equilibrium (linear non-equilibrium thermodynamics) are relatively well characterized. In this regime, it can be assumed that fluxes (e.g. material flow, reaction rate) are linear functions of their forces (e.g. density gradient, affinity). This can be expressed as J.-YL*** (2-16) where Ja is the flux, X^ are the forces and L^ are linear coefficients (the phenomenological coefficients). Equation (2.16) is known as the linear phenomenological equation. Phenomenological coefficients with identical indices (a=b) relate "conjugated" (Wisniewski et al. 1976, p. 29) fluxes and forces (e.g. material flux and the density gradient), while coefficients with a*b relate cross-effects on a flux from a non-conjugated force (e.g. affinity on material flux). The coefficients are given by , (c*b). (2.17) Lab rdJ ^ KdXbJXr Chapter 2. Turing Theory and Thermodynamics 21 In the late nineteenth century, it was observed that there were "reciprocal relations" in diverse phenomena (Lavenda 1978, p. 26), which can be stated as Lab~Lba. (2.18) This means that for two forces, a and b, a change in force a affects flux b to the same extent that a change in force b affects flux a. The reciprocal relations (2.18) were derived from theoretical considerations by Onsager (1931). Onsager assumed microscopic reversibility; a certain time-scale for fluctuations from equilibrium (longer than several collisions and shorter than a complete decay to equilibrium); and that these fluctuations obey the macroscopic relaxation laws. His derivation verifies equation (2.18) in cases for which the fluxes are time derivatives of the extensive state variables (Casimir, 1945). The local entropy production, s, can be related to the fluxes and forces by S'-^JJXJ. (2.19) This quantity is positive if irreversible processes are occurring in the system. Substituting equation (2.16) into equation (2.19) gives s-JJVi*,. (2.20) « J For the case of two forces it can be directly shown that minimization of the local entropy production with respect to the forces coincides with the steady state. This steady state is equilibrium if the forces can vary freely, but is not equilibrium if one of the forces is held fixed (Reichl 1980, p. 521). That the steady state is the state of minimized entropy production was first shown by Prigogine (1946, 1967b). Lavenda (1978, Ch. 5) has shown that this is really a state of least dissipation of energy, which coincides with minimized entropy production in the linear regime, but does not coincide far from equilibrium. The linear result can be generalized to n forces and fluxes. Chapter 2. Turing Theory and Thermodynamics 22 When equation (2.20) is expanded about the steady state (Reichl 1980, p. 522), a criterion for stability can be derived. The second term in the expansion is the excess entropy production, P, (due to the fluctuation from the steady state) P - i j r f V j ^ A X i A X } (2.21) where V is volume and the sum is over the forces which are not fixed. P is always positive because of the properties of L. However, its time derivative is always negative. The excess entropy production is, therefore, a Liapounov function which establishes the stability of the steady state in the linear regime. Non-equilibrium thermodynamics has been pushed into the far from equilibrium regime (Glansdorff and Prigogine 1971; Nicolis 1971; Reichl 1980, Ch. 17). First, the entropy is expanded about the steady state. This gives a second order term which represents the entropy associated with fluctuations from the steady state, as in the linear case (equation 2.21). For the far from equilibrium case, however, the matrix L cannot be said to be positive definite. Therefore, the second order term (which can be equated with the excess entropy production) cannot be proven to be positive. An assumption of local equilibrium is made to give a definite sign to the second order term and support further analysis in the far from equilibrium regime. It is argued that the second order term (which is the curvature of the entropy) is required to be negative if equilibrium holds locally. However, this requires the entropy to be a maximum at the steady state. The entropy is only truly a maximum at equilibrium in the case of an isolated system. The assumption of local equilibrium requires each locality to be isolated, which is not in accordance with the attempt to derive global stability conditions. For a non-isolated locality, the entropy may not be at a maximum, in which case the second derivative may not be negative. Further analysis in the far from equilibrium regime depends on showing that the time derivative of the second order term is positive, thus establishing this term as a Liapounov function defining the stability of the steady state. Lavenda (1978, Ch. 4) discusses the difficulties in the local equilibrium Chapter 2. Turing Theory and Thermodynamics 23 assumption for the far from equilibrium case, and demonstrates that the stability conditions placed on the second order term are postulates of the theory, not derived conditions. Thermodynamic stability analysis in the far from equilibrium regime is, therefore, full of difficulties. Lavenda (1978, Ch. 7) argues that the only meaningful stability analysis for far from equilibrium systems must come from kinetic arguments. The knowledge gained from thermodynamic balance equations complements that gained from kinetic analysis. Stability criteria derived from thermodynamics must always agree with kinetic stability analysis. But by linking thermodynamic knowledge with kinetics, greater understanding of the forces at work in far from equilibrium systems can be gained. Therefore, thermodynamics is still important in the study of far from equilibrium systems, but it cannot provide necessary and sufficient conditions for stability. Having discussed the thermodynamic aspects of far from equilibrium pattern formation, I will devote the rest of the thesis to kinetic considerations. The introduction in this chapter to the kinetics should be sufficient for an understanding of the rest of the thesis. For a slightly different approach to the derivation of growth rate and wavelength of pattern, see Harrison (1993). The linear algebraic approach used here can also be found in Edelstein-Keshet (1988) and Lyons (1991). As the specific models used in this thesis are nonlinear, Turing analysis is used to discuss the general features of the pattern formation. Solutions to the models are actually found through numerical analysis. Chapter Three is a discussion of the numerical methods used in this work. Chapter Three Numerical Methods 3.1 Introduction The work in this thesis is largely of a numerical nature. In working with reaction-diffusion theory, two general approaches can be made. The first is mathematical analysis. This can in principle provide fully general results, as compared to the specific results obtained by numerical analysis. The previous chapter contained the linear analysis first developed by Turing. These techniques are relatively easy to apply and will be used throughout the thesis. However, for any analytical information on the effect of nonlinearities, much more advanced techniques are necessary. Analytical techniques from dynamical systems, or bifurcation, theory have been applied to standing wave solutions of reaction-diffusion equations with some success (e.g. Lyons and Harrison 1991; Ermentrout 1991; Takagi 1982, 1986). These techniques are most effective near the bifurcation point, that is in the weakly nonlinear regime. Most techniques require approximations which limit their applicability. At present, nonlinear analysis is unable to describe much of the behaviour of reaction-diffusion systems in the strongly nonlinear regime. However, it is the nonlinear steady states that are of most interest in the study of biological pattern formation. While the large scale dynamics of reaction-diffusion systems can be described by linear analysis and some nonlinear effects can be described through nonlinear analysis, the great majority of the details of final pattern are left in the realm of numerical investigation. In this thesis I am interested in several aspects of the nonlinear effects in reaction-diffusion systems. My main theoretical objective is to establish the conditions for pattern localization in reaction-diffusion systems. In the next chapter, I will examine these conditions by Chapter 3. Numerical Methods 25 studying the effects of the diffusion terms, the reaction terms and the presence of precursor gradients. The first two parts, studying the reaction terms and the diffusion terms, involve investigation of very long-time, nonlinear pattern. At present, this type of investigation is possible only through computation. Study of the pattern effects of precursor gradients is beyond the scope of current nonlinear analytical techniques, and must also be studied numerically. Chapter Five involves the development of a chemical mechanism for the Gierer-Meinhardt model. This mechanism is found by computation to give stable pattern, while the earlier mechanism of Babloyantz and Hiernaux (1975) does not. Part Two of this thesis involves an application of the theory developed in Part One. In this modelling of axolotl heart formation, new requirements on the numerical technique become apparent. Specifically, in order to model some of the experimental data, it is found that the diffusivity of one of the species in the model needs to be quite fast. The simplest numerical techniques, which are used in the next chapter, are unstable for solution of equations with fast diffusivities. Therefore, more sophisticated techniques are used for the biological application. This chapter will provide an explanation of the numerical techniques used and a discussion of their specific application in the thesis. As well, the modelling of stochastic influences will be discussed. Finally, a short summary of the computing facilities used will be presented. 3.2 Numerical Solution of Reaction-Diffusion Equations Although the numerical solution of reaction-diffusion equations is very important, given their analytical intractability, there are few numerical techniques which are designed specifically for these systems. The basic nature of these equations is that reaction terms and diffusion terms must be solved simultaneously, with consideration of the accuracy and Chapter 3. Numerical Methods 26 stability characteristics of both types of terms. In practice, the reaction terms and the diffusion terms are treated separately, frequently by different techniques. In almost all numerical solutions of reaction-diffusion equations, the Euler technique is used for the reaction terms. The error associated with this method at each time step is relatively high (truncation after the first order term in a Taylor series expansion, see footnote 1 below). However, the method is fast and simple. Higher accuracy methods, such as the Runge-Kutta (truncation after the fourth order term of the Taylor series expansion; see Boyce and DiPrima 1986, p. 421), can be used to lower the error from the reaction term. I will describe the Euler technique for the reaction terms in some detail below. It was used in all computations in this work. Most numerical solutions of reaction-diffusion equations differ in the treatment of the diffusion terms. Several different techniques were used in this work and will be described below. All numerical techniques start with the discretization of the differential system of interest. This involves specification of a set of nodes which approximates the continuous system. For solution of the equations in one spatial dimension, the nodes are points on a line. In this thesis, only the case of uniformly spaced nodes will be addressed. In two dimensions, the set of nodes can be arranged in various geometries. The mesh can be hexagonal, square, etc. There are no computations of reaction-diffusion equations in more than two dimensions in this thesis. For a reaction-diffusion system discretization of the equations is relatively straightforward. Reaction terms are entirely local, so any mesh is acceptable. Diffusion occurs in accordance with Fick's law, so there is material exchange down any local concentration gradient. In the computation, neighbouring pairs of points are compared and material is exchanged depending on the direction and magnitude of the concentration gradient. Chapter 3. Numerical Methods 27 Restrictions on the mesh geometry in two dimensions arise for specific techniques of solution. Grids in which each node is connected to four other nodes are used in this thesis due to the fact that the technique (Peaceman-Rachford) used in Part Two flows in two specific directions. This is possible to implement on a rectangular grid. In practice, I have used square grids, due to their ease of implementation. Non-rectangular grids (e.g. hexagonal) can be used with certain methods, but the discretization is generally more complex than for the rectangular case. 3.2.1 Explicit Methods The simplest technique for numerical solution of differential equations is the Euler method. Two-morphogen reaction-diffusion equations involve solving for the values of the two morphogens, X and Y, at discrete time steps. An assumption is made that the reaction and diffusion terms can be treated separately. This means that the change in time of the morphogen concentration due to reaction is computed, then the change in time due to diffusion is computed, and finally these results are added together to get the overall change in time of the morphogen concentration. For reaction, the Euler technique simply means that Z"+1 = X" + At(dXn/dt) (3.1) where At is the discrete step in time and n is the number index of the time step (n and n+1 are superscripts, not exponents). The equation is only shown for the X morphogen, the Y equation is identical in form. This method is also known as the tangent line method (Boyce and DiPrima 1986, p. 399), since the value of X at the next time step is obtained by multiplying the first derivative of X at the current time step (the tangent line) by the time step and adding this to the value of X at the current time step. The value of the tangent is found from the rate equations of the reaction-diffusion scheme, since the rate equations Chapter 3. Numerical Methods 28 involve time derivatives of the morphogen concentrations. The error associated with the method at each time step is proportional to (At)2. While this seems a relatively high error for many applications, it is on the same order as the error in many of the techniques for diffusion. To limit the accumulated error when using this technique, small time steps are taken. The Euler technique can also be applied to the diffusion terms. For diffusion, this technique is also known as forward differencing, or explicit. This terminology arises from the fact that given the X concentration at every point in the system at the current time step, the X concentration for every point can be calculated for the next time step. Only nearest neighbour interactions will be considered in the discretized model of diffusion in this work. Thus, for the one-dimensional case each point exchanges material with two other points on the line, while on a two-dimensional rectangular grid each point exchanges with four other points. In one dimension, the diffusion equation is f-4# (3.2) dt ds where D is the diffusivity and s is the spatial variable (again, the equation is only shown for the X morphogen). This can be discretized1 to give •x;+1-2x; + xj_x Xf =XJ + D At (3.3) (As)2 where j is the spatial node index and As is the discrete spatial step between nodes. This equation can be understood heuristically in that the time derivative is as in equation (3.1) and the spatial second derivative is represented by material exchange between the ;th point and its neighbours. Discretization can be achieved formally by expanding X in Taylor series around the current point: one series in space expanded around a negative perturbation; one series in space expanded around a positive perturbation; and one series in time expanded around a positive perturbation. The two spatial series are added. The two series (one in space, one in time) are then truncated at the tenn following the derivative of interest (first derivative of time, second derivative of space). This gives expressions for the derivatives in terms of the lower order terms of the Taylor series. The terms not included in the expressions for the derivatives constitute the local (at every iteration) truncation error. For an explicit derivation, see Lyons (1991, pp. 65-66). Chapter 3. Numerical Methods 29 This scheme has several advantages. For many problems it can give fast solution, because there are relatively few arithmetical operations per time step. However, this scheme is only first order accurate with respect to the time step (error proportional to its square), so small time steps must be taken. Frequently, though, a balance can be found between the small time step and the simplicity of the computation to make this a favourable technique. Certainly the conceptual simplicity of the technique is convenient. Furthermore, boundary conditions are easily worked into this scheme. It is straightforward to specify the type of exchange with the j-1 point at the left hand boundary or the j+1 point at the right hand boundary. This can be readily generalized to two dimensions. As well as easily incorporating various boundary conditions, the explicit discretization can be extended to a variety of grid geometries. Because it is straightforward to incorporate varied boundary conditions in the explicit method, this is used in the next chapter. In Chapter Four, solutions of the Brusselator and Gierer-Meinhardt models are computed with periodic boundary conditions. This is easy to do with the explicit scheme, but quite involved for the implicit schemes used elsewhere in this work (Strikwerda 1989, p. 80). Besides the accuracy limits on the time step, the explicit scheme has a numerical instability. Stability analysis for differencing schemes was developed by von Neumann (see Press et al., 1988, section 17.1, or Strikwerda 1989, section 2.2). One expresses the solutions of the difference equations as the Fourier elements XI = |» e-*(^) (3.4) where \ is a complex number known as the amplification factor and GO is the wavenumber (as in the previous chapter), with index k. The time dependence of the solutions are integer powers of | . Therefore, the differencing scheme is stable for Sj <, 1. Substituting equation (3.4) into equation (3.3) gives the requirement that Chapter 3. Numerical Methods 30 7W*1 (3'5) for the explicit method to be stable. This means that the method cannot be used in the case of fast diffusivities or large time steps for a given spatial step. While stability with an explicit technique can typically be gained by lowering the time step, this was not a good option for the biological computations in Part Two. At the time steps used in that modelling, the large diffusivity of the fastest diffusing substance violated inequality (3.5). Maintaining use of the explicit technique would have meant reducing the time step roughly fourfold. However, the diffusivity of this substance was defined in a fixed relation to other variables in the model, namely the rate of tissue movement (see beginning of section 7.3) and the diffusivity of other species (e.g. the diffusivity of the slowest substance was 1000 times slower than that of the fastest diffusing substance). Slowing the whole computation by roughly four times would result in a significant increase in computational time. Therefore, the implicit methods of the next subsection were used for the biological modelling. The explicit scheme was used for the computations of Chapter Four, however, due to the use of periodic boundary conditions. For these computations, small time steps were taken to insure low accumulated error. The methods discussed in the next subsection are more accurate, and were used for the computations in Chapter Five for this reason. Inequality (3.5) can also be illustrated on a more local level which gives a better physical picture of the instability. This illustration is due to L.G. Harrison (personal communication). Imagine three adjacent points at which the concentration is defined (Figure 3-la). The endpoints Xt and X3 are fixed concentrations. The midpoint X2 is allowed to vary according to Fick's Law. The asymptotic (fc=°°) concentration of X2 will be the average X2m = (X3 + Xx)/2. An instability results if the time step carries the value of X2 beyond X2X. This can be expressed as (A*)max = X*. ~X2 = £ ± £ - X 2 = Z 3" 2f 2 + Z l . (3.6) Chapter 3. Numerical Methods 31 c o c V o c o O Distance ^2« c » u c o O Figure 3-1: A physical picture of the instability in the explicit technique for solving the diffusion equation, a) Consider three points on a line, each with different concentration. X3 andZ l 5 the endpoints, are of fixed concentration, while X2 varies according to Fick's law. X2<X) is the equilibrium value ofZ2. I fZ2 exceeds A^o,, a numerical instability results, b) is a plot of time, t, versus concentration of X2. Time is on the horizontal axis, and the concentration of X2 on the vertical axis. The origin represents the current (nth) time step. The curve represents the true exponential relaxation of X2 to the equilibrium value X2oo. Since the numerical technique approximates the value of X2 at the next (n+1) time step by the slope dX2/dt, too large a step forward in time, to point b, extrapolates the value of X2 past X2<x> resulting in the instability. At is too large in this case. If smaller time steps are taken, for example to point a, then the approximation is both more accurate and keeps X2 less than A o^o- From L.G. Harrison (personal communication). Chapter 3. Numerical Methods 32 The requirement for stability is expressed as (AX)max 2» AX. Using equation (3.3) for AX, this gives X3-2X2+X1 (X3-2X2+X,) 2 (As)2 K ' The concentrations cancel from each side, giving I DAt 2* (As) i*77ZT (3-8) which easily rearranges to expression (3.5). Figure 3-lb shows additional graphical representation of this problem. The value of X at future time in equation (3.3) is approximated by the slope at the present value of X (as in equation 3.1). Analytical solution of the diffusion equation for this 3-point system shows an exponential relaxation in time of X2 to the value X2«,. However, if the time step is too large, the approximation from the slope carries X2 beyond this asymptotic value and instability results. Only by using a small time step can the slope be re-evaluated frequently enough that the interpolated value of X2 remains close to the exponential relaxation of the exact solution. 3.2.2 Implicit Methods Due to the numerical instability (3.5) and the second order error with respect to the time step, explicit diffusion was only used for the computations in Chapter Four. Incorporation of the periodic boundary conditions used therein is much simpler with the explicit technique. For the biological modelling, techniques without time step instabilities were found to be faster. These techniques are also second order accurate in the time step, and were hence used for the computations of Chapter Five. Chapter 3. Numerical Methods 33 Fully Implicit Method The diffusion equation (3.2) can be spatially differenced in the same manner as equation (3.3), but differenced in time for future values (the Taylor expansion is around a negative perturbation in time). This gives the fully implicit scheme A)+l ~ZAj + Aj-1 yn+l Vn T Vrt+1 1 vn+1 , \rn+l At (As)2 (3.9) Substitution of equation (3.4) into this equation gives the stability condition 1 . 2/W.A5 1 +4a sin ' <:1 (3.10) where a = DAt/(As)2 . This differencing scheme is therefore stable for all time and space steps. This stability can also be shown by Harrison's local method. For the case of an implicit scheme, the concentration atX2 (Figure 3-la) is given at the next time step by 1 + 2-(As)2 This is just a rearrangement of equation (3.9) evaluated at the three local points. In the limit that At -* 00, this expression becomes X2 = ^ L ( 3 . i 2 ) which is exactly the definition of the asymptotic value, XlX , given above. The implicit method is again shown to be stable with respect to any size time step. However, the method is only first order accurate with respect to the time step. Higher order accurate methods are discussed below. Chapter 3. Numerical Methods 34 Solving implicit difference equations for Xn+1 involves the simultaneous solution of a system of linear equations. Thus, at each iteration a matrix must be solved to find the Xn+1. Luckily, the matrix is tridiagonal: it only has non-zero elements on the diagonal and at plus or minus one column index from the diagonal. In other words, equation (3.9) only involves elements of Xn+1 at the spatial indices j+l,j, and;'-l. Tridiagonal matrices are relatively fast to solve, as all the values of the matrix need not be stored. The method used for matrix solution will be discussed below. Crank-Nicolson Method To achieve better accuracy while retaining the stability characteristics of the implicit scheme, the explicit and implicit differencing methods can be averaged. This gives xf -x] D \(x£-ixf + xffl + (x;+1 -2X] + xnhl) At (3.13) (As)2 This method was developed by Crank and Nicolson (1947). It has second order accuracy with respect to the time step and is unconditionally stable. The von Neumann stability criterion is 1 - 2a sin . , v ' © i As 1 + 2a sin =sl (3.14) where a is as in equation (3.10). Increasing At or D always decreases the left hand side of inequality of (3.14). Thus, this method is used in the computations in Part Two in which very fast diffusivities are required to model the experimental data. As well, the second order accuracy with regard to the time step makes it a suitable method for the computations in Chapter Five. Equation (3.13) can be written as AX"+1 = BX" (3.15) Chapter 3. Numerical Methods 35 where A and B are tridiagonal matrices. Following the algorithm of Burden et al. (1981, p. 528; with boundary conditions modified from Dirichlet to no-flux), the matrices A and B are solved to give Xn+1. The method of LU decomposition is used (Press et al., 1988, section 2.3). Matrix A can be broken into a lower triangular matrix (L) and an upper triangular matrix (U). In the first step of the algorithm, the values of an intermediate vector Z are found (Z is defined by the relation LZ = BX"). These values are found from the defined value of Z at the left hand boundary (since there is only one element of L at the left hand boundary) and proceeding by forward substitution through the rest of the vector. The second step of the algorithm is to solve for Xn+1 according to U • X"+1 = Z . The algorithm starts with the defined value of X"+1 at the right hand boundary (there is only one element of U at the right hand boundary) and computes the rest of the vector through back substitutions. Due to the nature of the matrix solution, implicit schemes are not suitable for ill-defined boundaries. Also, since the algorithm specifically starts with the boundary values, use with periodic boundary conditions is complicated. Periodic boundary conditions can be used by an application of the Sherman-Morrison method for computing matrix inverses (Strikwerda 1989, p. 80). This involves use of additional vectors to specify the boundary conditions and solution of the equations for these vectors. By contrast, the diffusion equation can be solved at a local level (at each node) for the explicit method and implementation of periodic boundary conditions is straightforward. Peaceman-Rachford Method The Crank-Nicolson method can be extended to two spatial dimensions by including material exchange in the second direction. This is quite easy for a square grid with the spatial step equal in both directions. However, the matrices for this difference equation are sparse, but no longer tridiagonal. Therefore, each iteration requires more than the square (considering the number of nodes) of arithmetic operations compared to the one dimensional case. The Chapter 3. Numerical Methods 36 solution is inherently more complex (in terms of arithmetic and storage) in two dimensions than one. A method developed by Peaceman and Rachford (1955) gets around this difficulty by splitting the time step. For one half-iteration, the diffusion is computed in the x-direction, and for the other half-iteration the diffusion is computed in the y-direction. For this reason it is also called an Alternating Direction Implicit (ADI) technique. The total diffusion is still computed over the whole iteration. This scheme can be shown as *x? - (2+p)*rr+Ktj2 - -*r,-i+(2 - p n - K,* (3-i6a) X£t - (2 + p ) ^ 1 + X^ = -X^f + (2 - p)X%y2 - X£>2 (3.16b) where p = (As)2/.DA* is the reciprocal of the a used above (As is the same in both directions). The solution of each of equations (3.16) only involves the solution of tridiagonal matrices, thus avoiding the extra arithmetic and storage of the two-dimensional Crank-Nicolson technique. The Peaceman-Rachford technique is used in all two-dimensional computations in Part Two. Lawson and Morris (1977) published a further splitting of the Peaceman-Rachford technique which should save 30-40% of the arithmetic with a cost of extra storage of a matrix of the dimensions of the system (same number of nodes as X). Their scheme (applied to the Peaceman-Rachford technique), as published, is: Wfi1 = -X»hl + (2 - p ) ^ . - X»j+1 (3.17a) Kif - (2 + p)*r/1/2 + K™ = K-1 (3-17b) Xl% - (2 + p ) ^ ; 1 + X1%, = 2X$m - Wr? (3.17c) for the first iteration, p is as in the previous equations. W"+1 is the new intermediate matrix n+l required by this technique for solution of X" . The next iteration is given by Wfi2 = 2X?f - 2X£m + Wtf1 (3.18a) Kl'j2 - (2 + p)X"fn + x£j2 - WiT (3-18b) Xgt - (2 + p)X?f + Xlf+l = 2X£3/2 - W£2 (3.18c) Chapter 3. Numerical Methods 37 Concentrations at greater times are computed through iteration of equations (3.18). However, this scheme does not work as written. The technique gives negative concentrations and will not compute the diffusion equation. Step (3.17a) merely defines a matrix equal to the RHS of equation (3.16a). Step (3.17b) sets this equal to the LHS of equation (3.16a). The left-hand sides of equations (3.17c) and (3.16b) are equal. The RHS of equation (3.17c) is derived by recognizing that the RHS of equation (3.16b) is the negative of W"?1 (the LHS of 3.16a) minus a factor of X"J112. However, the factor of 2 reported in Lawson and Morris' paper is incorrect. Equation (3.17c) should read Z* - (2 + p ) ^ 1 + X& = -29X»?2 - Wtf (3.19) Equation (3.18c) will be identical in form, but evaluated at times greater than n+1. Equation (3.18a) corresponds to equation (3.17a), but evaluated at later time steps. Equation (3.18a) can be expressed in terms of the matrices computed in equations (3.17). Since equation (3.18c) corresponds to equation (3.17c), the form of equation (3.18a) will be correct for all time steps greater than n+1. The negative of the LHS of equation (3.17c), minus a factor of X"f, has the same form (but evaluated at the next time step) as the RHS of equation (3.17a). Again, the factor of X"*1 must be corrected from 2 to -2p. Performing these operations on the RHS of equation (3.19) (corrected equation 3.17c) gives a corrected expression for equation (3.18a): Wff = -2pJ^J1 + 2p^ ; 1 / 2 + Wtf. (3.20) When corrections (3.19) and (3.20) are implemented, the technique does solve the diffusion equation. In several tests against the straight Peaceman-Rachford technique, I found the Lawson and Morris scheme always slower. The tests consisted of starting with one half of a square at high (A=35) concentration and the other half of the square at low concentration (X=l). Diffusion between these two regions was computed for a set number of iterations. It was Chapter 3. Numerical Methods 38 found in all cases that the Lawson and Morris scheme was slightly slower than the Peaceman-Rachford. As well, the corrected scheme is far more difficult to programme than the Peaceman-Rachford (as it has three steps and an additional stored matrix compared to the two step Peaceman-Rachford). Thus, there is no good reason to use the Lawson and Morris scheme. For the reaction terms, the simple Euler (explicit) method was felt to be sufficient for all computations. The simplicity results in less arithmetic at each iteration, which speeds computation. The second order error in the time step is acceptable, as shown in the tests in Figures 3-2, 3-3 and 3-4. For diffusion, the Crank-Nicolson is not much more complicated than the fully implicit method. Therefore, the Crank-Nicolson was implemented for its accuracy characteristics. Using a higher order method for the diffusion increased accuracy for the whole method: use of a further higher order method in the reaction was felt to be unnecessary, given the results of the tests in Figures 3-2, 3-3, and 3-4. Overall then, Euler solution of the reaction terms was coupled with a variety of schemes for the diffusion terms. Explicit diffusion was used for the computations in Chapter Four, while the Crank-Nicolson (one dimension) and Peaceman-Rachford (two dimensions) techniques were used for all other computations. In Part Two, implicit techniques were used in order to conduct fast computations with high diffusivities. There has been some concern expressed recently about the coupling of an explicit reaction solver with an implicit diffusion solver in the same overall scheme (Ascher et al., 1995; Ruuth, submitted). Ruuth (submitted) has found that a first order accurate explicit reaction with fully implicit diffusion tends to give lower wavenumber solutions than the exact solution. His results are for single wavenumber errors at large wavenumber solutions (e.g., 12 peaks computed instead of the exact 13). He does not address computations of reaction-Chapter 3. Numerical Methods 39 delt=0.075 X delt=0.3 75 c o ca L. * - > c a? o c o O Distance Figure 3-2: Comparison of the effects 6*f the time step for the two most commonly used values in this thesis. The heavy dashed line is the solution at 33,000 iterations for Af=0.075. The lighter dashed line is the solution for Af=0.375 at the same time. The spatial step, As, is 1/80 for both cases. Chapter 3. Numerical Methods 40 diffusion patterns with pre-existing gradients, nor does he address computations with random input at every iteration. He has also found that a Crank-Nicolson method with second order accurate reaction tends to give misleading results. In this case, higher frequency modes tend to linger too long. This is due to a misrepresentation of nonlinear terms by the numerical method. He found that this can lead to computed results which are a reflection of the exact solution. Errors of this sort can be hard to detect, as decrease in the time and space steps can still lead to similar erroneous results. For the first problem above, Ruuth (submitted) recommends use of higher order schemes. For the problems with the straight Crank-Nicolson, he recommends use of modified schemes which split the spatial differencing over more time steps. Both of these suggestions require more arithmetic per iteration, and hence longer computation for the same time step, compared to the simpler methods. As well, the problems which he describes are somewhat far from the concerns of this thesis. The only chapter in which high wavenumber patterns will be addressed is Chapter Four, in which solutions of different chemical models will be compared with each other, but not with exact solutions. Furthermore, I have used a fully explicit technique for these computations, due to the boundary conditions. The only inaccuracy in such a method is the well-known second order error with respect to the time step. At small time step, these errors are negligible. In the computations in which I have used implicit diffusion, reaction-diffusion patterns from pre-existing gradients are studied. In these cases, there is extreme forcing of the pattern from the gradient. It would take error greater than that found by Ruuth (submitted) to give erroneous results with these computations. In my initial work, I ran computations at a number of time and space steps. A general rule of thumb in numerical anaylsis is to run a computation at a given time step and then to halve that time step and see if the results agree within a small error. This can be done similarly for Chapter 3. Numerical Methods 41 c o L_ c o c o o 12 10 8 6 4 2 -0 mtiinuiiHimmiii i i i i i i i i i i i i i i i i i i i i i inrnT ™™"'« IIIHIlHmHmmMIMHmi***. Distance Figure 3-3: The effects of changing the spatial step are shown. The solid line shows As=l/80, the value used in most computations in this thesis. The dashed line shows the computation to the same time with the spatial step halved: As=l/160. The time step is Af=0.075 in both computations. Chapter 3. Numerical Methods 42 the spatial steps. All time and space steps used in my computations have been tested in this way and found to give accurate results. For instance, in Part Two two different time steps, Afc=0.375 and A*=0.075 are used. For Euler reaction with Crank-Nicolson diffusion, there is no significant difference in the results at these time steps (Figure 3-2). As well, there is no significant difference in computed results between a spatial step, As, of 1/80 (which I have used throughout Part Two) and one of 1/160 (Figure 3-3). At the suggestion of Ruuth, I tested a At of 0.0075, one-tenth of my usual value. At the smaller spatial step of 1/160, there is no appreciable difference from the results at Af=0.075 (Figure 3-4). These tests show the Euler reaction/Crank-Nicolson diffusion method to be accurate for the computations that I have used it for. These computations differ from Ruuth's in that a pre-existing gradient is used and thermal noise is simulated at every iteration. 3.3 Random Input In addition to numerical solution of the differential equations, random input was added at every time step to simulate the thermal fluctuations present in any chemical system. The modes contained in the simulated fluctuations are amplified according to the Turing dynamics. In this way, the equations solved are really stochastic differential equations. The concentrations are calculated from the deterministic reaction and diffusion terms, and then receive random kicks. The concentration dynamics are similar to the velocity dynamics in the Langevin equation for Brownian motion. The random number generator supplied in the C language libraries was used to fill an array with random numbers between 0 and 215. An element of this array was then chosen with another random number. In this way, a sequence of calls for random numbers will be effectively free of sequential correlations (Press et al., 1988, Ch. 7). New random numbers are generated for every time step. The random numbers obtained by the above procedure are Chapter 3. Numerical Methods 43 delt=0.0075 X delt=0.075 c o ' • ^ CD C <D O C o O 12 10 h 8 6 4 2 0 i i i i i i i i i i i i i i i i i i i i i i i i i M i i i i i i i m i i i i i i i r r r • n T m r r i n i i i i i i i i i i i i i i Distance Figure 3-4: Further testing of the time step Afc=0.075. The solid line is the solution at 33,000 iterations with the usual time step of Af=0.075. The dashed line is the solution at 330,000 iterations for Ar=0.0075. The spatial step is As=l/160 for both computations. Chapter 3. Numerical Methods 44 scaled to values between -1 and 1 and then multiplied by a factor so that their magnitude will be approximately 1% of the initial X concentration. This is an approximation of the noise level in a Poisson distribution, assuming that the initial X concentration is in the 10"9 mol/1 range. This level of noise correlates with the connection made between the spatial step in the computations and real length in Part Two (section 7.3) with regard to the biological modelling. There, the spatial step is identified as 0.025 mm. This corresponds to a volume of 2xl0"8 ml. Multiplying by the assumed concentration of 10"9 mol/1, this gives 2xl0'20 moles per node in the computation. This is roughly 12,000 molecules. If the concentration at the node obeys a Poisson distribution, then the standard deviation is on the order of ^/12,000 s l l O . Therefore, fluctuations on the order of 1% of the initial concentration are allowed at each node. For further discussion of the relation between computational and real units, see Harrison et al. (1984). The biological data modelled in Part Two contains stochastic aspects not directly related to the thermal noise. Two main quantities are of interest in the axolotl heart experiments. These are the percent embryos that develop heartbeat, and the time to develop heartbeat. The ratio of these two quantities defines the Heart Determination Coefficient (HDC) of Jacobson and Duncan (1968). This expression is used as a scale of embryo viability for a given experiment. From the standpoint of modelling, however, the use of the HDC confuses the stochastic effects present in the data. The percent beating hearts term contains variability in a multitude of factors. For instance, in normal heart development it includes variability in the rate of relative tissue movements. For experiments done in culture, variability in tissue geometry and tissue-to-tissue contact become incorporated in the HDC (see Figure 7-9). The magnitude and distribution of these stochastic factors is generally not well characterized. Therefore, it is impractical to attempt to model these factors through use of the random number generator. Chapter 3. Numerical Methods 45 The time to heartbeat term is more deterministic and hence more suited to application of the reaction-diffusion model. The modelling in Part Two achieves quantitative agreement with this aspect of the experimental data. There is still a stochastic aspect to this term: for all embryos which do attain heartbeat, there is still variation in the time to beat. Attempts to model this are difficult because numerous factors could influence the distribution of times to beat. Two major causes could be variation in the tissue velocities (again) and thermal fluctuations in the morphogen concentrations. Variation in the tissue velocities affects the amount of morphogen produced by a given stage of development. Thus, retarded tissue movement can lead to late beating (or none, as discussed above), while advanced tissue movement leads to early beating. Approximate error limits on tissue velocities in normal development are known. However, incorporating this into the model is complicated by the effect of thermal fluctuations in morphogen levels on the time to beat. For a given noise level, different seeds to the random number generator give different times for a morphogen to reach a given concentration (peak firing). The seed affects the time of peak firing, since different seeds give different random kicks to the crucial nodes at which the peak is forming. The site of peak formation is always the same in the biological modelling, due to the use of precursor gradients in the reaction-diffusion model. The values of the morphogen concentrations are critical in the region of the gradient maximum, and random fluctuations in those values can significantly change the time of peak firing. This is somewhat counterbalanced by the fact that the variability of firing times found by variation of the random seed is diminished in reaction-diffusion systems with growing precursor gradients. This will be discussed in more detail in section 4.3.2. Given just these two causes of variation in times to beat, it becomes apparent how difficult a separation of their effects becomes. The variation due to thermal effects in the model does not have a one-to-one relationship with the variation in the experimental parameter. The Chapter 3. Numerical Methods 46 experimental variation also reflects variation in tissue velocity. It is already a difficult task to model the effects of the tissue movement on the model concentrations (Chapter Seven). To include stochastic tissue movement would add a new level of complexity to this task. As there is no way of directly correlating these effects to the experimental fluctuations, the necessary impetus to incorporate stochastic effects into the tissue velocities is absent. The approach in Part Two has instead been to model the mean values of the experimental times to beat. The same random seed was used in all computations. In this way, a given computation always gives the same time to beat. This enables the comparison of computations of different experiments on a relative basis. The model is used deterministically to fit the experimental data, including the relations between the times to beat for the different experiments. Different configurations of the morphogens, corresponding to different random seeds, would give longer or shorter times to beat, but retain the relations between the different experiments. Quantitative agreement is reached between computation and experiment in the relative times to beat for all experiments. Due to the complexity of the large scale stochastic effects in the biological data, the only fluctuations directly addressed in this thesis are thermal fluctuations. This random input to the morphogen concentrations is necessary for Turing pattern formation, as discussed in the previous chapter. In this section, the levels of random input have been correlated with morphogen concentrations in the nanomolar range. This range is not crucial to the results of the modelling, but represents a reasonable estimate of initial morphogen concentrations involved in embryogenesis. Chapter 3. Numerical Methods 47 3.4 Computing Facilities All computations were conducted on a Silicon Graphics Iris Indigo workstation with a MIPS R3000A processor running IRIX. All programmes were compiled with the C language compiler available from Silicon Graphics. Graphs of position versus concentration for one-dimensional computations were written with C-language functions contained in the Silicon Graphics Graphics Library. In two dimensions, results were visualized either as surface plots, with concentration represented on the vertical (z-) axis, or as colour maps, with concentration represented as a colour (or grey) scale. Dataview, a package written at the University of British Columbia, was used for the colourmapping. This package also provides easy interactive graphics during the running of a computation. This two dimensional animation is reasonable only with at least the speed and storage capabilities of a small workstation such as the Iris Indigo. The interactive graphics of the one-dimensional or Dataview package in two dimensions was a key to the completion of this project. Since this was a chiefly numerical work, the fast processing of numerous parameter values was necessary. This would have been incredibly time-consuming without interactive graphics in both one and two dimensional computations. Some two dimensional output is presented as surface plots, which allows for the simultaneous visualization of two concentration profiles. For this, Explorer, a product of Silicon Graphics, was used. This package is not amenable to animation during computation time. It is used after completion of computation on saved data sets. Therefore, the computation was animated with Dataview and important stages of the data saved. This was then processed through Explorer to provide surface plots. Both two dimensional packages provide colour output. However, these have been translated into greyscale or black wireframe for use in this thesis. Only one colour plate has been used in this work (Figure 8-12). Chapter 3. Numerical Methods 48 This concludes discussion of the methods used in this thesis. The remaining two chapters of Part One present results on different aspects of reaction-diffusion pattern formation. Chapter Four Localization of Reaction-Diffusion Pattern 4.1 Introduction The term localized pattern is used to denote a pattern element which stands alone and without periodic repeats. In biological development, there are many examples of both localized and periodic pattern. Many large organs (e.g. heart, liver) form singly and without repeats, in contrast to such patterns as the somitic division of the dorsal mesoderm in vertebrates, or segmentation in insects1. All of these patterns are tightly controlled: any error in formation of patterns, of either type, can be lethal to the developing embryo, or lead to large-scale defects in the adult. There are numerous types of pattern forming mechanisms which can describe the events of development (see Harrison, 1993). However, within the context of this thesis, I will concentrate on reaction-diffusion dynamics. As computations have shown, reaction-diffusion can produce both types of pattern: localized morphogen peaks surrounded by large areas of low concentration, and regular patterns of repeated peaks (see Lacalli, 1981). These different types of pattern have been produced by different nonlinear reaction-diffusion models. In the case of repeated pattern, there are nonlinear models whose behaviour approximates very closely to the sinusoidal solutions of the linear Turing system. However, it is difficult to describe the localized pattern produced by other nonlinear models in terms of the linear theory. In Part Two, I develop a reaction-diffusion 1 The examples from the animal kingdom represent the bias of this thesis. There are certainly myriad examples of important pattern-forming events of both types in development in other kingdoms. Chapter 4. Localization of Reaction-Diffusion Pattern 50 model of vertebrate heart development. This project requires an understanding of the characteristics of reaction-diffusion equations which lead to localized pattern. In this chapter, I investigate computationally the effect of the various components of two different reaction-diffusion models on final pattern. There has been recent interest in the area of long range inhibition (leading to pattern localization) in reaction-diffusion systems. Computations of the simple Gray and Scott (1983) model have shown a tendency to form large areas of inhibition around local activator maxima (Pearson, 1993). This has now been observed for a chemical system in the laboratory (Lee et al., 1994). However, these patterns involve travelling waves. Localization occurs because the inhibitor is very fast diffusing with respect to the activator. Initially, activator and inhibitor peaks form in phase, but the fast inhibitor diffusion allows a wave of inhibition to escape from the front of the travelling wave. When this occurs, wave propagation is stopped or strongly altered. In such a system, a pattern of localized peaks does occur, but without much stability. A continual growth and decay of pattern is seen in these systems. In the current work, I am interested in stable pattern. This is not provided by the propagation and annihilation of travelling peaks. However, the principle of fast inhibitor diffusion to create localized structures can be applied to models and parameters which form stationary pattern. In this chapter, I explore the characteristics of the Gierer and Meinhardt (1972) model, developed from the concept of lateral inhibition in cybernetics, but applied to biological development. These equations have been used chiefly to model localized patterns. Chapter 4. Localization of Reaction-Diffusion Pattern 51 In characterizing the properties of the Gierer-Meinhardt model, I will compare and contrast it with the Brusselator (Prigogine 1967a; Prigogine and Lefever 1968). These two models are representative of two major classes of reaction-diffusion model (Meinhardt 1982; Harrison 1993, section 8.2.2): activator-inhibitor (Gierer-Meinhardt) and depletion (Brusselator). The different nonlinearities in the two models lead to differences in pattern formation, and has led to use of the Gierer-Meinhardt for modelling localized structures, while the Brusselator is more often used for repeated patterns. First, I will investigate the pattern effects of the diffusion terms. Changing the diffusivity ratio DY/DX = n can have pattern effects similar to those of changing the nonlinear kinetic terms of particular models. It is therefore important to study the effects of the diffusion terms before discussing the characteristics of particular models. Second, I will investigate the effects of precursor gradients. These effects are found to be model specific, and can be related to the properties of the reaction terms in the models. The effects of the reaction terms will be discussed last, with emphasis on the intrinsic properties of the two models which lead to their differences. 4.2 Diffusion Terms One of the necessary conditions for pattern formation in reaction-diffusion systems is that the diffusivities of the two morphogens be different. This requirement can be seen explicitly in equation (2.15). For Dx=Dy, the wavenumber becomes infinite. Computations of reaction-diffusion equations have typically used n (the ratio of DYID^) values of 45 (Gierer and Meinhardt 1972; Lacalli et al., 1988), and lower (e.g. Lacalli and Harrison 1979; Murray Chapter 4. Localization of Reaction-Diffusion Pattern 52 1989). Although a variety of n's have been used, a systematic investigation of the effects of changing this ratio has not been published. I have computed patterns with the Brusselator and Gierer-Meinhardt models at a series of diffusivity ratios for constant linear wavelength. In one-dimensional computations, roughly sinusoidal patterns occur at low n for both models. The patterns become increasingly irregular in spacing as the diffusivity ratio is raised. As well, the peaks become sharper in comparison to the troughs as n is raised. The same phenomena are more striking for patterns in two dimensions, and I therefore concentrate on these. I have used two methods for quantifying degree of order: the R-parameter of Clark and Evans (1954) and the radial distribution function. For evaluation of peak sharpness, a simple programme was written to find the distance at which the concentration fell below some fixed fraction of the peak concentration. Peak sharpness is a local measure of pattern localization at the level of a single peak. Measurement of pattern order can demonstrate a more system-wide localization of pattern. However, this measure is more indirect, in that ordered patterns are associated with simple combinations of modes which lead to repeated pattern, while disordered patterns result from more complex combinations of modes. These complex combinations can lead to localized waveforms. Chapter 4. Localization of Reaction-Diffusion Pattern 53 4.2.1 Evaluation of Pattern Order and Peak Sharpness The Clark and Evans (1954) R-parameter gives a measure of how far a pattern deviates from being completely random. It is based on nearest neighbour distances. The average nearest neighbour distance in a completely random arrangement is given by where p is the overall area density of structures. The average nearest neighbour distance in the observed pattern is defined by r -2: where the r's are all the nearest neighbour distances in the pattern and N is the total number of structures. The R-parameter is R = r°>*/ . (4.3) / ran R=l for a random distribution of structures. For a more clustered packing of structures than random, R<1 (Lacalli and Harrison, 1978b). Ordered structures without clustering give higher R (>1). The maximum value, R=2.1491, is achieved in a hexagonal array; for a square array, R=2. The R-parameter is simple to calculate and serves as a quick measure of order in a pattern, provided that visual inspection has fairly clearly eliminated orderly clustering. The radial distribution function g(r) is a commonly used measure of order in fluids (McQuarrie, 1976). It has also been used to measure the order in biological patterns Chapter 4. Localization of Reaction-Diffusion Pattern 54 (Markovics et al., 1974). It is the ratio of the average density of structures at a distance r from an arbitrary structure to the average overall density: g(r)=l for a completely random distribution. G(r) gives much more structural information than the R-parameter. For very fine subdivision of the r scale, g(r) versus r is a line spectrum with lines at nearest neighbour, second neighbour, etc. separations. A plot of g(r) versus r contains information not only of order in the pattern itself, but of the type of order. Roughly random patterns will generate plots in which g(r)=l at all r. Hexagonal arrays give plots in which g(r) has peaks at r=l, and at >/3/, 21, etc, where / is the nearest neighbour distance. G(r) goes asymptotically to 1 as r -» oo : the plot looks like one of damped oscillations. To calculate g(r), first, local distribution functions &(r) = - ^ r (4.4) are calculated for each (ith) structure. These are calculated for a discrete set of distances, each representing an annulus of area A/r). nj(r) is the number of structures in the annulus, and p is again the overall density. These are averaged to give the radial distribution function at each r S(r)--nrflAi(r)gi(r)- (4.5) A(r) is defined by A(r) = JjA;(r). (4.6) The results can then be displayed as a plot of g(r) versus r. Chapter 4. Localization of Reaction-Diffusion Pattern 55 A programme was written for evaluation of one-dimensional peak sharpness. This serves as one measure of pattern localization. The sharper the peak, the faster the morphogen concentration falls off with distance from the peak. For each peak, the distance from the peak to the first point below 1/4 of peak height was found2. The distance was measured separately for each side of the peak, to avoid complications arising from the asymmetry of peaks due to discretization. These distances were then all averaged to give a measure of the mean peak sharpness. 4.2.2 Models and Methods Two parameter sets were used with each model. The rate equations for the Brusselator are given by AY — =aA- bBX + cX2Y -dX+ DXV2X (4.7a) dt — = bBX - cX2Y + DYV2Y (4.7b) dt where V2 is used for the Laplacian operator. The mechanism leading to this model will be presented and discussed in the next section. One set of parameters corresponded to a point in region C of the Lacalli parameter space (Figure 4-1). For the other parameter set, I used the values from Harrison and Kolar (1988), which corresponds to a point in region D of the 2 The peak height is measured as the difference between the peak value and the average trough value. The actual concentration at 1/4 peak height is therefore obtained by dividing the peak height by 4 and adding the average trough value. Chapter 4. Localization of Reaction-Diffusion Pattern 56 Figure 4-1: Location in the linear space (Lacalli and Harrison, 1979) of the parameter sets used in this chapter. BRl and BR2 are the two sets used in the Brusselator, GMl and GM2 are the two sets used in the Gierer-Meinhardt. w=10 in this diagram. Chapter 4. Localization of Reaction-Diffusion Pattern 57 Lacalli space (Figure 4-1). These parameter sets will henceforth be referred to as BRl and BR2, respectively, and are summarized in Table 4-1. Table 4-1: Brusselator Parameters Set BRl BR2 a 0.01 0.05 A 1 1 b 0.026 0.5 B 1 1 c 0.01 1.8 d 0.007 0.05 At 0.2 0.01 Different time steps are used to compensate for different linear growth rates (maximum eigenvalues) in the two parameter sets. The Gierer-Meinhardt equations are given by dX „ cSX2 — = aS + dt Y BY -\xX+DxV2X — = c'S'X2-vY + DYV2Y. dt Y (4.8a) (4.8b) The specific kinetics of this model will be discussed in section 4.4. For these computations, I used one parameter set close to that of Figure lc in Gierer and Meinhardt (1972) and another parameter set from Lacalli (1981). These will be referred to as GM1 and GM2, respectively. Both parameter sets are in region C of the Lacalli space3. The parameter sets are summarized in Table 4-2, and plotted in Figure 4-1. It can be shown that only a narrow range of the decay constants, [i and v, will put the Gierer-Meinhardt model into region D of the linear parameter space. Very noisy pattern resulted in the computations run with these constraints. Chapter 4. Localization of Reaction-Diffusion Pattern 58 Table 4-2: Gierer-Meinhardt Parameters Set GM1 GM2 a 0.0006 0.0006 S 0.4 0.4 c 0.05 0.05 V 0.009 0.21 c' 0.065 0.025 S' 0.4 0.4 V 0.01 0.27 At 0.125 0.025 Again, the time steps differ in order to compensate for the different maximum eigenvalues of each parameter set. For both models, the equations were computed on a 60 node by 60 node square domain. In order to minimize boundary effects on pattern formation, periodic boundary conditions were used. Due to its ease of implementation for these boundary conditions, the Euler method was used to solve the diffusion terms, as well as the reaction terms. The time step was kept small, so as to stay well below instability (3.5). 4.2.3 Results Computations were conducted with the above parameters in the two models, for increasing n value. Linear wavelength was kept constant at 0.125 times system length by changing/^ and Dx at constant n in order to facilitate comparison of patterns. All computations were iterated until the pattern had shown no significant change in about half the computing time. This required from 20,000 to 560,000 iterations, depending on the specific parameters used. All parameter sets show a trend of decreasing order with increasing n. Clark and Evans R-parameters and radial distribution functions were calculated for each pattern. Figure 4-2 Chapter 4. Localization of Reaction-Diffusion Pattern 59 shows computations with BR2 at the n values of 15, 30 and 100. At an n of 15 (Figure 4-2, top), the X pattern is very close to hexagonal. Increasing n to 30 (Figure 4-2, middle) destroys some of the regularity in the pattern, but a rough hexagonal arrangement can still be seen. When n is increased to 100 (Figure 4-2, bottom), the pattern becomes nearly random. This decrease in order can be seen in the radial distribution functions for these computations. Figure 4-3 shows plots of g(r) versus r for the patterns in Figure 4-2. The r interval used was 1 unit. To this accuracy, the peak at r=8 corresponds to the linear wavelength X=7.5. In a perfect hexagonal lattice, the sequence of peaks should be, to 1-unit accuracy, 8, 13 ,15 ,20 ,23, etc. Only the first, third and fifth of these values appear as peaks in Figures 4-3a and b. The second and fourth neighbour values are missing. A probable reason for this can be seen by viewing Figure 4-2 (top computation) from close to the plane of the paper. The peaks tend to line up in parallel lines with proper spacing at /, 21, 31, etc., but to have adjacent lines incorrectly staggered for hexagonal packing, which disturbs distances which are diagonals in the rhombic representation of the lattice. These parallel lines of peaks suggest incipient stripe formation. Because of the type of nonlinearities in the model, however, a pure striped pattern does not form (Lyons and Harrison 1991, 1992). In Figure 4-3b, the regular pattern can still be seen in the plot of the radial distribution function for n=30. However, the plot is not Figure 4-2: Patterns computed with the Brusselator model with parameters BR2 (see Table 4-1). The two spatial dimensions are represented in the x- and y-directions, with the concentration of the X morphogen shown as a 256 value greyscale. White represents high X concentration, while black represents low concentration. Each computation was conducted on a 60 by 60 grid with periodic boundary conditions. Top) n=15; middle) n=30; bottom) n=100. Chapter 4. Localization of Reaction-Diffusion Pattern 60 Figure 4-2 Chapter 4. Localization of Reaction-Diffusion Pattern 61 Figure 4-3: Radial distribution functions for each of the computations in Figure 4-2. a) g(r) plotted against r (in grid points) for the computation in Figure 4-2 (top). The width of the annuli used in the calculation of g(r) was 1 grid point. The plot is truncated at a distance of 42, since the distribution of peaks was found to be random in all computations at longer distances. The three peaks at integer multiples of the linear wavelength (7.5 grid points) are indicative of the regular pattern observed in Figure 4-2 (top), though not of hexagonal pattern, as discussed in the text, b) Radial distribution function for the pattern in Figure 4-2 (middle). The amplitude of g(r) is somewhat diminished in comparison to a), indicating a decrease in order, c) Radial distribution function for the pattern in Figure 4-2 (bottom). This plot is characteristic of a random distribution of peaks beyond an inhibitory field distance of about 5 grid points (approximately 2k/3). Chapter 4. Localization of Reaction-Diffusion Pattern 62 as sharply defined as that in Figure 4-3a. The peaks in the plot are lower than at n=15, indicating a lower level of order in the pattern. Finally, the radial distribution function for n=100 (Figure 4-3c) is characteristic of a random pattern, except for an "inhibitory field" of radius about 5 (approximately 2k/3). This represents a minimum separation of peaks. The inhibitory field concept originated in the work of Wigglesworth (1940) on patterns of bristles on insect integuments. Claxton (1964) followed up the concept in relation to hair follicles on Australian sheep, and did a computation by hand, obtaining R=1.757. Lacalli and Harrison (1978b) discussed inhibitory fields with time-dependent radii in relation to pore patterns in the cell walls of desmids. The concept is well known to biologists, who commonly refer back directly to Wigglesworth (e.g. Sater and Jacobson, 1990). The closeness of the relationship between Turing models and inhibitory field models, especially at high inhibitor diffusivity, has not been generally recognized. This relation was mentioned in Holloway et al. (1994). The trend of decreasing order with increasing n can also be seen with the Gierer-Meinhardt model. Figure 4-4 shows computations with GM2. The n values shown are 15, 30, and 60. The radial distribution functions are shown in Figure 4-5. As with the Brusselator, these show a progression from nearly hexagonal pattern to nearly random arrangement. The Figure 4-4: Patterns computed with the Gierer-Meinhardt model, using parameters GM2 (see Table 4-3). The spatial dimensions and concentration of the X morphogen are plotted as in Figure 4-2. The computational methods used are also identical to Figure 4-2. Top) n=15; middle) n=30; bottom) n=60. Chapter 4. Localization of Reaction-Diffusion Pattern 63 Figure 4-4 Chapter 4. Localization of Reaction-Diffusion Pattern 64 I i I o I 1 i I I .1 L-l .1 I I Figure 4-5: Radial distribution functions for the patterns shown in Figure 4 -* a) shows g(r) plotted against r for the pattern in Figure 4-4 (top). The plot indicates a relatively well-ordered pattern, in agreement with visual inspection of the computation, b) Radial distribution function for Figure 4-4 (middle). This plot shows very little order in the pattern, although there is a certain peaking of g(r) at the first multiple of the linear wavelength, c) This plot, corresponding to Figure 4-4 (bottom), shows even less order than b). This plot is indicative of a random distribution of peaks, with an inhibitory field effect. Chapter 4. Localization of Reaction-Diffusion Pattern 65 Gierer-Meinhardt tends to become random at lower n than the Brusselator. This is the reason why n=60 is shown instead of n=100. This difference between the models will be discussed in section 4.4. Although the contrast between the two models shows that kinetic nonlinearities are important, linear analysis can be of some aid in explaining the effect of changing n. In Lacalli's parameter space (Figure 2-3), the slope of line a, the Turing bifurcation, is the diffusivity ratio n. Any set of kinetic parameters is represented by a fixed point inside this space (Figure 4-1). As the diffusivity ratio is increased, the Turing bifurcation is pulled farther and farther away from the parameter point inside the space. As a result, more modes have positive eigenvalue (for a different demonstration of this, see Murray 1989, Fig. 14-10). The dispersion plots for the Brusselator parameters used in Figure 4-2 are shown in Figure 4-6a. The result is that the final waveform, which is a mixture of these modes, becomes more complex. Figure 4-6b shows this effect for the Gierer-Meinhardt parameters of Figure 4-4. The remainder of the computations conducted are summarized in Table 4-3. This lists the Clark and Evans R-parameters for the four parameter sets run (Brusselator and Gierer-Meinhardt) at the various diffusivity ratios. The same information is shown graphically in Figure 4-7. Chapter 4. Localization of Reaction-Diffusion Pattern ——• n*15 - - - n=30 66 — " • n»100 9 3 « > C o a iH mil I I i milium ininii innnnnnnmin inii intnn in mi mini 3000 60O0 9000 12000 15000 n*15 - - - n-30 — — n=60 I 0.20 0.10 0.00 -0.10 -O "Pd '""""""" 3000 6 0 0 0 9 0 0 0 Wavevector squared 12000 Figure 4-6: Dispersion plots for the computations in Figures 4-2 and 4-4. The eigenvalue o+ is plotted against the square of the wavevector, co2 = co2 +co2. Different lines are used for different diffusivity ratios. Each set of curves illustrates that the range of wavevector values with positive eigenvalue, and hence the number of pattern modes with positive growth rate, increases with n. This is one contributing factor to increasing pattern disorder. However, this is only linear analysis, and contains no information about the different mode interactions which result from different nonlinear terms in the two models, a) BR2: all eigenvalues are real; b) GM2: eigenvalues are complex at low wavevector, to the left of the plotted curves. Chapter 4. Localization of Reaction-Diffusion Pattern 67 Table 4-3: Clark and Evans R-Parameters4 n 10 15 30 45 60 100 BR1 1.72 1.71 1.65 1.64 1.60 1.60 BR2 1.84 1.88 1.82 1.78 1.82 1.64 GM1 1.82 1.71 1.58 1.50 1.52 1.54 GM2 1.78 1.80 1.64 1.56 1.53 1.53 The R parameter is useful in that it allows the degree of order to be represented quantitatively by a single number. However, the simplicity of this representation is balanced by some significant limitations on the information conveyed by R. Specifically, the lower limit of R drops below 1 if there is clustering (Lacalli and Harrison, 1978b), but rises above 1 if there is an inhibitory field effect: a sort of "anti-clustering". In the present work, the highest n results (n=100 for the Briisselator and n=60 for the Gierer-Meinhardt) are best described as random, except for the minimum separation imposed by the inhibitory field. This is clearly seen both in the patterns and in the radial distribution functions. The R parameter is decreased (towards 1) by randomness, but increased by the inhibitory field. It thus reaches a lowest value of about 1.5, which does not give information about the two opposing effects. 4 The R-parameters for the computations in Figures 4-2 and 4-4 are underlined. Chapter 4. Localization of Reaction-Diffusion Pattern 68 0> E co CO a i BR1 hex. BR2 square - GM1 - Claxton •" GM2 - random 2.30 2.10 1.90 1.70 1.50 1.30 1.10 0.90 j _ j _ J L J_ J L J_ J L _l_ J_ J L 10 15 20 25 30 35 40 45 50 55 60 65 70 76 80 85 90 65 100 n Figure 4-7: Plots of the Clark and Evans (1954) R parameter versus diffusivity ratio n, for all four parameter sets. Values of R for hexagonal pattern, square pattern, Claxton's (1964) computed inhibitory field pattern, and random pattern are marked with horizontal lines. The curves show numerically the decrease in order with increasing n, but illustrate also a limitation of the utility of the R parameter. The R values stay above 1.5, when there is an inhibitory field effect, event though the pattern has a very random appearance. This is a graphical representation of the data in Table 4-3. Chapter 4. Localization of Reaction-Diffusion Pattern 69 The average peak sharpness at n=15 and n=100 was calculated for one-dimensional computations with the four parameter sets. The results are summarized in Table 4-4. All parameters sets showed a shorter average distance to 1/4 X height at higher n. There is no overall difference between the two models. BR2 gives the broadest peaks, while BR1 gives the sharpest peaks. Both Gierer-Meinhardt parameter sets fall between the Brusselator sets in terms of peak sharpness. Table 4-4: Average Distance to 1/4 X Peak Height (in grid points) n 15 100 BR1 2.000 1.000 BR2 2.875 1.800 GM1 2.500 1.100 GM2 2.375 1.125 This sharpening of peaks is also observed in two-dimensional pattern, as is readily seen in Figures 4-2 and 4-4. It can be understood intuitively in that the range of the activator becomes very short with respect to the overall pattern wavelength (a function of both diffusivities). In these computations, the pattern wavelength is kept constant as the relative range of the activator is shortened by raising n. The peak is effectively becoming more localized within the measure of the linear wavelength. 4.2.4 Discussion The intent of this chapter is to assess the relative importance of quantitative parameter values and qualitative types of nonlinearity for the formation of localized pattern surrounded by long-range lateral inhibition. In this section, it has been shown that large values of the diffusivity ratio DY/DX= n lead to the formation of a kind of partly-ordered pattern with large Chapter 4. Localization of Reaction-Diffusion Pattern 70 regions in which no structures (peaks) form. This shows a basic type of lateral inhibition. The two models studied are not equally sensitive to n: the Gierer-Meinhardt model is more easily induced to generate disordered pattern than the Brusselator. In Figure 4-4 (bottom), lateral inhibition is quite definite. The pattern contains some spaces nearly two wavelengths in diameter in which no peaks form. Ordered and partly-ordered two-dimensional patterns have commonly been approached in quite different ways, even within the realms of mathematics or computation. The Wigglesworth (1940) concept of the inhibitory field leads to computations such as those of Claxton (1964) in which pattern points are initiated sequentially, one at a time, in random positions except for an inhibitory circle (created instantaneously upon point initiation) surrounding each established point. No new points can be initiated within these inhibited regions. This process of pattern formation does not contain the concept of a mode, or superposition of modes, to which all points belong. The latter concept is the basis of computations on Turing (1952) type reaction-diffusion models with mode interaction afforded by various nonlinearities stemming from specific chemical models. The present work has shown that the latter type of computation can yield patterns of the Wigglesworth type when inhibitor diffusion is dominant. The results suggest a recognition that most of the principal elements of the inhibitory field and the Turing concepts are the same: self-enhancing activation; generation in the activated region of a substance which spreads rapidly; and the action of that substance as an inhibitor of the activation process. There are two elements of the Turing theory missing from the inhibitory field model. First, Chapter 4. Localization of Reaction-Diffusion Pattern 71 the diffusion of the activator, which must be similar to, but slower than, that of the inhibitor. This is essential for the formation of orderly pattern. Second is the time dependence of the size of the inhibitory field. This is implied from the physical basis of the inhibitory field in outward diffusion of a substance from the activated point. Inhibitory fields with diffusional time-dependence were considered by Lacalli and Harrison (1978b). These were shown to be capable of producing quite disorderly pattern (R=1.35). The present work suggests that the high diffusivity ratio limit of two-morphogen reaction-diffusion theory may be the same as inhibitory field theory. Whether one must use the full two-morphogen mathematical formulation in discussing a pattern, or whether one may use the conceptually simpler inhibitory field idea depends very much on the degree of order in the patterns being considered. Low degree of order is amenable to the simpler discussion. As well, it has been shown that the activator (X) peaks become sharper as the diffusivity ratio, n, is raised. Computations were run at high and low n with the four parameter sets used in this section (Table 4-4). These demonstrate a type of pattern localization, in that the activator peaks become a smaller portion of the chemical wavelength as n is raised. For this measure, pattern localization is shown not to be model specific: neither model inherently forms sharper peaks than the other. Chapter 4. Localization of Reaction-Diffusion Pattern 72 4.3 Precursor Gradients In the previous section, it was shown that perturbation of the non-model specific diffusion terms can give different results with different reaction-diffusion models. These differences are due to the specific reaction terms in each model. Before addressing the specific reaction dynamics which lead to these differences, I would like to further explore the differences between the models in response to other types of perturbation. In this section, I will present results on pattern localization by spatial gradients in morphogen precursors. Although the previous work in this thesis has involved reaction-diffusion dynamics from homogeneous initial conditions, a patterned initial state can be easily incorporated. Hunding (1987, 1989) has made the distinction between reaction-diffusion models with no pre-pattern ("Turing models of the first kind") and models which include pre-pattern ("Turing models of the second kind"). For Turing models of the first kind, linearization about the homogeneous steady state is straightforward, and has provided the basis for nonlinear analytical techniques (Lyons 1991; Lyons and Harrison 1991, 1992). For Turing models of the second kind, linearization and nonlinear analysis have not yet been formulated. Linear properties, such as the chemical wavelength, can only be defined at each point in the system. This is similar to the local equilibrium assumption of far-from-equilibrium thermodynamics discussed in section 2.2. This local linear analysis is helpful in several ways, but cannot be regarded as a precise description of the system. Chapter 4. Localization of Reaction-Diffusion Pattern 73 Precursor gradients have been used in reaction-diffusion models for over twenty years. The Gierer and Meinhardt (1972) model was introduced with precursor gradients. This was effectively used in their first application, localization of the head structure in the coelenterate Hydra, and in much subsequent work. Gradients have also been used in the Brusselator (Herschkowitz-Kaufman, 1975) to limit pattern formation to some fraction of the system. In that example, the gradient carries the system in and out of the pattern-forming region of the linear parameter space (across the Turing bifurcation), effectively limiting the region wherein peak formation can occur. Such localization has been used more recently to model segmentation patterns in the fruit fly, Drosophila melanogaster (Lyons et al. 1990; Harrison 1993). As well, precursor gradients have been shown to cause the Brusselator to form striped patterns in two dimensions (Lacalli et al., 1988). 4.3.1 Methods In this work, I am interested in studying the effect of gradients on model dynamics. In order to do this more clearly, all computations are conducted within the pattern forming region of the Lacalli parameter space. This ensures that the action of long range gradient effects upon the model dynamics is not simply turning off Turing type pattern formation. While the latter use of gradients has proven effective in modelling (as in Herschkowitz-Kaufman, 1975), it gives an incomplete picture of gradient effects on reaction-diffusion dynamics5. For simplicity, I will only discuss one-dimensional gradients in this section. 5 In the first demonstration of long-range inhibition (Gierer and Meinhardt 1972, fig. lc) the gradient nowhere takes the model to the wrong side of the Turing bifurcation. Chapter 4. Localization of Reaction-Diffusion Pattern 74 Many different shapes of gradient can be used for these "Turing models of the second kind". The linear gradient is an obvious starting point. Also commonly used, and known to exist in biological systems (Driever and Niisslein-Volhard, 1988), is the exponential gradient. Both of these will be used in the present comparison of gradient effects in the Gierer-Meinhardt and Brusselator models. In Chapter Seven, I develop a different shape of gradient, based on the experimental data available for axolotl heart induction. If a chemical is released from a localized source and destroyed at a separate localized sink, the resulting steady-state concentration profile is linear (see Edelstein-Keshet 1988, p. 409). This type of gradient is easily included in the precursor concentration in a reaction-diffusion programme. Because of this, linear gradients are used in many of the computations in this section. Exponential gradients can be established by release of a chemical from a point source, with first-order decay throughout the system (see Harrison 1993, section 2.1). Exponential gradients provide a chemically established distance measure: the distance at which the concentration is lie of the source concentration is given by jD/k, where k is the rate constant for the first-order decay and D is the diffusivity. From the early years of this century it has been apparent, through experiment, that gradients control important events in biological development (Boveri 1904, 1910; Child, 1915). For reviews, see Huxley and de Beer (1934, Ch. 8) and Sander (1994). In recent years, some of the chemicals forming these gradients have been elucidated. For example, gradients of Chapter 4. Localization of Reaction-Diffusion Pattern 75 retinoic acid have been implicated in numerous patterning events in many organisms. The earliest retinoic acid gradient characterized was that associated with limb generation in the chick (Tickle et al. 1982; Eichele and Thaller, 1987). For a review of some of the biological gradients characterized to date, see Ding and Lipshitz (1993). In Part Two, I review the evidence for a gradient of inductive capacity in the formation of the axolotl heart. The role of a precursor gradient in the localization of this structure in the biological system corresponds well to the role of precursor gradients in reaction-diffusion models. Computations were conducted with the BR2 and GM2 parameter sets in one spatial dimension. No-flux boundary conditions were used since outside constraint on peak formation might result from Dirichlet conditions and periodic conditions cannot be used with a monotonic gradient. 4.3.2 Results Self-effects The Gierer-Meinhardt and Brusselator are first compared in terms of parameters that are shared by the two models. Such parameters are the zeroth order rate constants for production of X, designated a in both models, and the first order decay constants for X, called u, in the Gierer-Meinhdardt and d in the Brusselator. Other constants are either not shared by the two models, or are involved in cross-catalysis, which is dynamically quite different between the two models. Chapter 4. Localization of Reaction-Diffusion Pattern 76 Gradients in a tend to give equivocal results and do not allow for good comparison between the models. In the Gierer-Meinhardt, increased a leads to ill-defined, noisy pattern. There is some decrease in peak number from that predicted from the linear wavelength, but the pattern is too noisy to be useful in biological modelling. With the Brusselator (both BR1 and BR2 parameter sets) a cycle of peak number was found. That is, the initial pattern had peak number close to that of the linear wavelength. As the computation progressed, peaks were eliminated. At later time peaks formed again, giving pattern close to the linear wavelength. These computations show that a gradients do not give the smooth, stable pattern needed for biological modelling. However, this points out a difference between the gradients used in this thesis and those which carry the system back and forth across the Turing bifurcation. Herschkowitz-Kaufman's (1975) localization of Brusselator pattern was achieved with an a gradient, and the pattern produced was smooth and regular. Such differences in response to a gradients could stem from different gradient slopes. In attempting to eliminate peaks, I used a gradients with a ratio of 6 between the maximum and minimum values. Herschkowitz-Kaufman used exponential gradients, with a ratio of about 1.6. Computation with \i (Gierer-Meinhardt) and d (Brusselator) gradients gives results that appear more significant. Figure 4-8 shows the effect of these gradients upon GM2 and BR2. These computations were run at n=15 to minimize the disordering effect of high diffusivity ratio. A linear wavenumber of 3 was used. In both models, the steepest gradient possible Chapter 4. Localization of Reaction-Diffusion Pattern 77 Y xS X 2 0 Y xS d X20 • • 1 1 1 1 1 • * i ' 1 1 1 1 • 1 1 * *%L • • • • Y XB mil x5 • • • Y X5 •mi xS Distance Distance Figure 4-8: The effect of a gradient in the linear decay term is shown for the Brusselator and Gierer-Meinhardt models. The linear wavenumber is 3 for all computations, and n=15. a) the Brusselator (BR2) with no gradient, d=0.05. b) BR2 parameters with a gradient of d=0.1075 to d=0.033. c) the Gierer-Meinhardt (GM2) with no gradient, u.=0.21. d) GM2 parameters with a gradient of ^=0.27 to u.=0.11. Chapter 4. Localization of Reaction-Diffusion Pattern 78 (while remaining in the pattern forming region of linear parameter space) eliminates a half peak. The Gierer-Meinhardt pattern becomes slightly more localized towards a dominant peak at the low end of the gradient than does the Brusselator pattern. However, the pattern remains one of multiple peaks in both models. More robust localization methods are needed for biological modelling of single structures. Cross-effects The next step is to investigate gradient effects in the cross-catalytic terms. Since the dynamics of cross-catalysis are different in the two models, direct comparison of gradient effects for a specific parameter is not possible. In each model, there are two terms to investigate. These are: c (with S held constant) and c' (with 5" held constant) in the Gierer-Meinhardt; and b and c in the Brusselator. The Gierer-Meinhardt v has no direct analogue in the Brusselator, and gradients in this parameter have not been investigated, v has a strong effect on the position of the system in parameter space, as does u., and would therefore be predicted to have similar gradient effects (see the discussion on the next page linking linear behaviour with gradient effects for c). In the Brusselator, gradients in b and c do not produce substantially different effects from gradients in d. For both the BR1 and BR2 parameters, the steepest gradients possible (while remaining in the pattern forming region) give only a half peak deletion from the linear wavelength, as in Figure 4-8 for the d gradient. With the BR2 parameters, a gradient in b (steepest possible) actually gives an increase of a half peak from the linear wavelength. Chapter 4. Localization of Reaction-Diffusion Pattern 79 The Gierer-Meinhardt responds somewhat differently to gradients in the cross-catalytic terms. Significant peak deletion from the linear wavelength is found for gradients in both c and c'. Figure 4-9 shows c gradients with the GM2 parameters. The computation with no gradient is shown in Figure 4-8c. With a gradient from c=0.075 to c=0.025, a half peak is deleted (Figure 4-9a). If the slope is nearly doubled, so that the gradient goes from c=0.1 to c=0.000125, then another full peak is deleted (Figure 4-9b). This is a much stronger response than any seen with the Brusselator. Some of the difference between the models with respect to c gradients can be seen through comparison of the resulting trajectories in parameter space. In the Brusselator, k\ and kf4 are proportional to the square root of c. In the Gierer-Meinhardt, k\ and A^ 4 contain expressions in which c cancels out at high values. The result is that a moderate c gradient in the Brusselator will produce a trajectory which traverses the whole of the pattern forming region, while in the Gierer-Meinhardt the trajectory asymptotically approaches a point in the pattern forming region for higher and higher values of c. There is no limit, from these considerations, upon the slope of the c gradient in the Gierer-Meinhardt. Practically, it is found that no further peak deletion can be accomplished with larger slope. In the case of a linear wavenumber of 3, linear c gradients could not cause more than one and a half peaks to be deleted. This was found with both the GM2 parameters and a slight variation of the original Gierer-Meinhardt parameters. Chapter 4. Localization of Reaction-Diffusion Pattern 80 Y x5 c x20 I I I J Distance Figure 4-9: The effects of a c gradient in the Gierer-Meinhardt model (GM2 parameters) are shown, a) linear gradient of c=0.075 to c=0.025. b) linear gradient in c with nearly double the slope shown in a): c=0.1 to c=0.000125. c) an exponential gradient in c, with a maximum ofc=0.1. Chapter 4. Localization of Reaction-Diffusion Pattern 81 Gradients in c' in the Gierer-Meinhardt also have strong effects. With the GM2 parameters, a gradient of c' from 0.1 to 0.0001 leads to deletion of two and a half peaks for a linear wavenumber of eight. In comparison, a c gradient of the same magnitude leads to deletion of two peaks. In this case, the c' gradient has stronger effects. However, c gradients were used in Part Two due to easier explanation as putative chemistry of the biological system. For axolotl heart formation, the process of inter-tissue induction that precedes differentiation is modelled as increase in magnitude of a gradient in c, which is a concentration gradient of an enzyme. To keep the same effects with a c' gradient would necessitate a decrease in c' over time. This would imply a patterned depletion of c', which is possible, but not as conceptually simple as the patterned growth of c due to tissue movement and inductive signalling developed in Chapter Seven. Exponential gradients Pattern localization through precursor gradients has been shown in the Gierer-Meinhardt, but to what extent is this dependent on gradient shape? For comparison, computations have been run with exponential gradients. In general, exponential gradients tend to localize pattern to a greater extent than do linear gradients. For moderate decay constants, the precursor concentration is lower at a given position than is the case for a linear gradient with the same maximum (and which stays positive throughout the system). Figure 4-9c shows a computation with the GM2 parameters. The exponential gives a lower c concentration at the position of the second peak in Figure 4-9b. This is sufficient to repress peak formation, resulting in the half peak pattern shown. Chapter 4. Localization of Reaction-Diffusion Pattern 82 Gradient growth In the biological modelling of Part Two, the effects of inter-tissue induction increase over time. Induction is modelled with a c gradient in a Gierer-Meinhardt type model. The increase is modelled as a growth of c magnitude in time. This growth of the c gradient has been found to be a powerful localizer of morphogen pattern. The gradient from c=0.1 to c=0.000125 (as in Figure 4-9b) was generated by incrementing all values of c by 2xl0"5 per iteration for 5000 iterations. Instead of the one and a half peaks obtained in Figure 4-9b with a static gradient, the growing gradient resulted in a half peak pattern (similar to Figure 4-9c). Growing gradients also show slope effects: generation of a gradient of c=0.075 to c=0.025 by the same increments resulted in a one and a half peak pattern, instead of the two and half peak pattern shown in Figure 4-9a for the static gradient. For the axolotl modelling in Part Two, gradient growth has a significant effect on the final pattern. If a static gradient is used, the single peak pattern is very broad and noisy. Gradient growth serves to limit the area of tissue differentiation, represented by the morphogen peaks, from the area of specified tissue, represented by the gradient. Achieving this difference in areas is very important in fitting the experimental data. Gradient growth may be a general method of pattern localization in embryogenesis. A telltale indication of gradient growth could be tissue movement during an inductive event, as in axolotl heart formation. Gradient growth serves to localize pattern through establishing a temporal advantage for certain peaks over others. For instance, in the half peak pattern obtained by growth of a c=0.1 to c=0.000125 gradient, the half peak is established before any other parts of the system have sufficient c for (rapid) peak formation. An inhibitory field grows outwards from Chapter 4. Localization of Reaction-Diffusion Pattern 83 this established half peak as Y diffuses. This is sufficient to repress formation of nearby peaks, even though nearby regions eventually reach high enough c levels for rapid peak formation. For these patterns, temporal primacy gives spatial primacy. In the Brusselator, c gradient decay (lower c gives slower destruction of Y) results in no change in pattern from the static case. This difference will be addressed in the next section. Gradient growth also smooths the effect of the random seed on peak firing time (see section 33). Test computations were run with a Gierer-Meinhardt model with a c gradient in which the time at which the X peak passed the value of 5 was recorded for four different random seeds. This gave a range of times from 1317 iterations to 3311 iterations. When the same gradient was generated by incremental growth (same rate as above), the range of times went from 1228 iterations to 1238 iterations. This decrease in variability of firing times with growing gradients increases the confidence with which predictions can be made from the axolotl model in Part Two. 4.4 Reaction Terms In the previous two sections, I have discussed two major ways in which reaction-diffusion pattern can be localized. In addition to presenting general results on the effects of precursor gradients and the diffusivity ratio, I have shown a difference in the way the two models respond to these perturbations. In this section, I would like to discuss the reaction dynamics which distinguish the two models and how these dynamics affect the specific pattern forming properties of each model. Chapter 4. Localization of Reaction-Diffusion Pattern 84 A well-known distinction between the Brusselator and the Gierer-Meinhardt is that the former is a depletion model, while the latter is an activator-inhibitor model. This chemical mechanistic classification corresponds to the signs of the Turing constants in the linearized systems. Given the Jacobian of the Turing system at the homogeneous steady state: k3 K4 two possibilities exist for the signs of the matrix elements. These are + -' + -which characterizes an activator-inhibitor system, and + + (4.14) (4.15) (4.16) characterizing a depletion system (see Edelstein-Keshet, 1988, p. 295). In a depletion model, V (U and V are used, as in Chapter 2, for the dependent variables in the Turing equations) has a positive effect on U (k2>0). This can be realized in a chemical system in which Y is necessary in the autocatalytic formation of X. As Y is used up (reflected in k3<0), X autocatalysis slows and peak growth is halted. Depletion models form U and V patterns out of phase with one another. By contrast, in an activator-inhibitor system, V has a negative effect on U (k2<0). In a chemical system, this can be achieved by direct inhibition by Y of the formation of X. In turn, X catalyzes Y production (reflected in &3>0). Therefore, as X increases, it more effectively produces its own inhibitor, and peak growth is terminated. The activator-inhibitor model forms U and V peaks in phase with one another. Chapter 4. Localization of Reaction-Diffusion Pattern 85 Can the dynamics characterizing the two classes of reaction-diffusion model account for the different effects in the Brusselator and the Gierer-Meinhardt presented earlier? I think that the two models serve to exemplify these two major types of reaction-diffusion model, and that by study of the dynamics of the nonlinear models, some general understanding can be gained. It has been pointed out that the depletion dynamics in the Brusselator are analogous to the kinetics of radical chain reactions (Harrison 1993, Ch. 9). The Brusselator equations stem from the mechanism: A^-^X (4.17a) B + X—b-^Y + D (4.17b) Y + 2X—C—*3X (4.17c) X—*-+E. (4.17d) Steps (4.17b) and (4.17c) correspond to the propagation steps in a radical chain reaction, where the intermediates are repeatedly formed and destroyed. The major difference of the Brusselator from the kinetics of radical chain reactions is that step (4.17c) is autocatalytic in X. Propagation steps such as (4.17b) and (4.17c), in which X and Y are alternately formed and destroyed, do not change the sum X+Y. In the rate equations (4.7), the terms -bBX and +cX2Y in the X equation are balanced by the same terms with opposite signs in the Y equation. Harrison (1993, Ch. 9) has pointed out that in a system with no-flux boundaries this conservation of X+Y implies that peaks above the spatially homogeneous steady-state must be of equal area to the troughs below the steady-state. This restricts the capacity of the Brusselator to form highly localized peaks surrounded by long-range inhibition. This is reflected in the generally sinusoidal, ordered pattern observed for low to moderate n, though localized pattern can be forced by high n (section 4.2), in which case the morphogen Chapter 4. Localization of Reaction-Diffusion Pattern 86 concentrations in the inhibited areas are very close to the homogeneous steady-state values. By contrast, the Gierer-Meinhardt equations (4.8) have no such balance of terms between the two rate equations and correspondingly no balance between the formation and destruction of the two morphogens in putative chemical reaction schemes (equations 5.4). Peaks and troughs, measured as deviations from the homogeneous steady-state need not be equal in area. This allows the model to form localized pattern more easily than the Brusselator, as seen in the computations of section 4.2. Figure 4-10 is a schematic description of pattern growth and termination in the two models. Figures 4-10a, b, and c are a time series of a Brusselator computation. Due to the balanced cross-effects in the model, X growth is always closely balanced by Y depletion. This balance is maintained throughout the system. As Y is depleted, X growth and cross-catalysis slow, and the system comes to a patterned steady-state. In the Gierer-Meinhardt, the process is more closely described as peak 'firing' (to use Gierer and Meinhardt's terminology). In the initial stages (Figure 4-10d), the growth of X is greater than the growth of Y (especially for nonzero o, or cS[Y>c'S'). The X peak escapes somewhat from the inhibitory effects of Y and becomes large. With large X, cross-catalysis becomes important and Y growth overtakes that Chapter 4. Localization of Reaction-Diffusion Pattern 87 Figure 4-10: Schematic depiction of pattern growth and termination in the Brusselator and Gierer-Meinhardt models. The arrows represent the relative rates of peak growth. The X morphogen is shown with a solid line, the Y morphogen with a dashed line. For the Brusselator: a) departure from the homogeneous steady-state; b) growth of pattern amplitude; c) the patterned steady-state. For the Gierer-Meinhardt: d) departure from the homogeneous steady-state, with X much faster growing than Y; e) growth of pattern amplitude, with X growth approximately equal to Y growth; f) slowing of pattern growth, with Y growing faster than X The morphogens eventually reach a patterned steady-state after a period of oscillations in amplitude. Chapter 4. Localization of Reaction-Diffusion Pattern 88 of X (Figure 4-10e). This leads to the slowing (Figure 4-10f) and eventual reversal of X growth. The steady-state is attained after several oscillations in peak amplitude. The oscillations are characteristic of region C (complex eigenvalues) of linear parameter space, in which the Gierer-Meinhardt model is run (see footnote 3). With such peak 'firing' dynamics, the final pattern is determined more by the peak regions than is the case for the Brusselator. This difference in symmetrical versus asymmetrical cross-effects can be used to interpret the results of the previous two sections. I suggest that the higher degree of order in two-dimensional Brusselator pattern than in the Gierer-Meinhardt, with the same n, is due to the greater degree of symmetry in the kinetics of the Brusselator. Because the X and Y concentrations stay balanced throughout the system, the pattern tends to amplify system-wide modes more than the Gierer-Meinhardt. For the peak dominated picture of the Gierer-Meinhardt, X escapes Y inhibition at the peak during firing, but stays inhibited in non-peak regions. These dynamics tend to establish peaks where random fluctuations are high, but not to emphasize amplification of the system-wide modes. Note that these pictures of the two models emphasize their differences: both models amplify sinusoidal modes at low n and are better described by inhibitory fields at high n. The susceptibility of the Gierer-Meinhardt model to precursor gradients can also be interpreted in terms of its asymmetrical cross-effects. In the Brusselator, gradient effects were small in all terms. This was partially due to the limitations placed on gradient slopes by the requirement of keeping the whole system in the pattern forming region of the linear parameter space. As well, though, the balance of cross-effects throughout the system could Chapter 4. Localization of Reaction-Diffusion Pattern 89 counteract the asymmetrizing effects of an applied gradient. These small gradient effects were also seen for gradients in the zeroth order term, a, and in the decay term, u, a self-effect, in the Gierer-Meinhardt. However, strong effects were observed for gradients in the cross-effects, c and c'. One reason for this is that there is no limit on the magnitude of the gradient slope for these parameters. It is also possible, though, that gradients in these asymmetrical cross-effects tend more strongly to form asymmetrical pattern than do gradients in self-effects. For the case of c gradient growth (or decay, in the case of the Brusselator) the first peak established in the Gierer-Meinhardt computation dominates the pattern at later times. With the Brusselator, there is enough balance of X and Y growth throughout the system that new peaks are not suppressed by the first peak. Some further comments on the uniqueness of the two models studied should be made. As discussed above, the Brusselator mechanism is closely related to mechanisms for radical chain reactions. How important are the specific details of the initiation (4.17a) and termination (4.17d) steps in the Brusselator? This can be answered in part by the Schnackenberg (1979) mechanism: X —!± A (4.18a) B^^Y (4.18b) 2X + Y—^3X. (4.18c) (4.18c) is identical to (4.17c), retaining the essential autocatalysis and balance of terms. (4.18b) does not require X for cross-catalysis of Y, and so the balanced bBX term in the X and Y equations of the Brusselator is absent in the Schnackenberg model. The zeroth order Chapter 4. Localization of Reaction-Diffusion Pattern 90 production and first order decay of X has been rolled into the equilibrium (4.18a). The Schnackenberg model is therefore a simplified and pared down version of the Brusselator, with only one of the balanced terms retained. Though it is unknown whether the Schnackenberg has all the pattern forming properties of the Brusselator that have been explored in this chapter, it does successfully form pattern and has been used in modelling (Murray 1989; E. Dulos, personal communication). On the other hand, the Gierer-Meinhardt model is quite sensitive to changes in model terms. Alteration of the c'S'X2 term, the cross-catalysis of Y by X, destroys the ability of the Gierer-Meinhardt model to halt pattern growth. I tested models with both c'S'X and c'S'X2IY cross-catalysis. Both models resulted in unstopped pattern growth. This shows that the peak 'firing' dynamics of the Gierer-Meinhardt model are very sensitive to alterations of the cross-catalysis. Their model was very carefully crafted to achieve termination of pattern growth. 4.5 Conclusion In this chapter, I have presented results on several different methods of pattern localization in reaction-diffusion models. The effect of the diffusivity ratio was explored in section 4.2. This is caused by variation in the Fickian diffusion terms, which are identical in all reaction-diffusion models. Therefore, the pattern localization observed at high n holds for all models. This localization was measured through X peak sharpness, and indirectly through pattern Chapter 4. Localization of Reaction-Diffusion Pattern 91 order. No differences were observed for X peak sharpness between the two models. However, the models did respond differently in terms of pattern order. This is due to differences in the reaction terms of the specific models. In section 4.4, I discussed the two major classes of reaction-diffusion models. Their different linear dynamics imply nonlinear interactions for slowing and termination of peak growth. The Brusselator and Gierer-Meinhardt models each fall into a different dynamic class, and this can serve as a starting point for discussion of their different pattern forming properties. It is suggested that the balanced (due to the form of the depletion) propagation steps of the Brusselator lead to well-regulated pattern growth. The more peak driven dynamics of the Gierer-Meinhardt do not regulate growth as well, and tend to form more random pattern. Section 4.3 presented results on the susceptibility of the two models to precursor gradients. It was found that the effects were minor for the Brusselator and for self-effects in the Gierer-Meinhardt. However, strong gradient effects were observed for gradients in the cross-catalytic parameters of the Gierer-Meinhardt. These differences certainly have some basis in linear dynamics, but could also result from the different types of dynamic balance in the two models. In light of the results in this chapter, it is interesting to re-evaluate patterns obtained by previous workers. In Murray's (1980, 1981) computations of animal coat patterns with the Thomas (1975) model, diffusivity ratios of 7 and 10 were used. These low values agree well with the adaptability of pattern to initial and boundary conditions observed in his work. The early morphogen peaks do not become strongly established and dominate later pattern, as has Chapter 4. Localization of Reaction-Diffusion Pattern 92 been observed in the present work at high n. In another case, one of the first examples of pattern localization in a reaction-diffusion system was Figure lc of Gierer and Meinhardt's original (1972) paper. That figure, as well as subsequent modelling of localized biological structures with similar parameters, has served to characterize the Gierer-Meinhardt model as one of localized pattern. However, inspection of the parameters used for Figure lc in Gierer and Meinhardt (1972) shows that non-model specific methods of pattern localization were used, as well as the model specific precursor gradient. Using an n of 45 and a linear wavelength of 0.6 system length gave a long range of the inhibitor and low wavenumber, emphasizing pattern localization. However, the Brusselator (BR2) will not give this half peak pattern with similar parameters. It is the model-specific gradient effect which gives the half peak pattern for the Gierer-Meinhardt. The localization properties developed in this chapter will be used in the application to heart formation in Part Two. Specifically, a modified Gierer-Meinhardt type model (developed in Chapter Five) is used in conjunction with a growing precursor gradient. The diffusivity ratio effect is not used: low n is used to fulfill some of the requirements imposed by experimental data. With the model and gradient localization, single peak patterns are generated, despite a linear wavelength of approximately 0.27 times system length6. This ability to reduce the actual wavenumber from the linear prediction gives a robustness to reaction-diffusion models of biological pattern formation. With this ability, identical patterns can be generated for a range of morphogen diffusivities (which change the linear wavelength). In the case of The gradient used in Part Two takes the middle half of system into the pattern-forming region of parameter space. Therefore, this gradient is something like that used by Herschkowitz-Kaufman (1975). Such constraint of the pattern to a portion of the system follows from the experimental data for heart induction. Chapter 4. Localization of Reaction-Diffusion Pattern 93 axolotl heart formation, enough experimental data exists to place restraints on the morphogen diffusivities, as will be discussed in Chapter Eight. For preliminary modelling of a biological phenomenon, though, this versatility is very useful. Chapter Five A Chemical Basis for the Gierer-Meinhardt Model The results of Chapter Four show that the Gierer-Meinhardt model inherently tends to form localized patterns of morphogen concentration. Although it is possible to form localized patterns with the Brusselator model, the Gierer-Meinhardt does this much more robustly. For this reason, the Gierer-Meinhardt model is used in the application of reaction-diffusion theory to the problem of vertebrate heart formation in Part Two. There is an inherent problem for chemists and biochemists in using the Gierer-Meinhardt model. The original model cannot be derived from a chemical mechanism. Gierer and Meinhardt's emphasis was on the dynamics which lead to lateral inhibition in cybernetics. Though they addressed the problem of localized structures in biological development, their paper was published in the journal Kybernetik. They found that the representation of reaction-diffusion dynamics in their rate equations can produce lateral inhibition. This stems from asymmetries of exchange between the morphogens, X and Y, as discussed in the previous chapter. For this asymmetry, they conceived of a direct inhibition of X by Y, with an unequal cross-catalysis of Y by X. The inhibition enters the rate equations as the X^fY autocatalytic/inhibitory term. While it produces the desired lateral inhibition, this term cannot be derived from chemical considerations. It contains the essence of direct chemical inhibition, but its infinite value at Y=0 shows that it cannot stem from any real chemical mechanism. One of the aims of modelling axolotl heart formation is to help future experimenters find the dynamics and the chemicals which constitute the process. In order to do this, it is crucial that the dynamic modelling connect to a possible chemical mechanism. Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 95 5.1 Babloyantz-Hiernaux Mechanism Babloyantz and Hiernaux (1975) published a mechanism with rate equations which approximate to the Gierer-Meinhardt reaction terms. This mechanism actually contains three chemicals with varying concentrations: the two usual Turing intermediates, and a third, P, which allows for exchange between the morphogens. To be consistent with the usage in Part Two, I will use A and / for the morphogens (for Activator and Inhibitor) instead of the previous X and Y. Their scheme is 5—^-* A (5.1a) 2A + S'^-^I + 2A (5.1b) A —B—» destruction (5.1c) /—i:—> destruction (5. Id) I + P 7 ^ A (5.1e) S + A + P—^P+2A. (5.1f) If equation (5.1e) is at equilibrium, then P = k2AlklI. This can be substituted into the rate equation for reaction (5.If), with k^ = kz = 1. If rate equations are written for A and / from the whole scheme and diffusion terms added, the Gierer-Meinhardt equations are produced: ¥L = aS + ^ ^ - - MA + DAV2A (5.2a) dt I A K — = c'SA2 - v/ + Z),V2/. (5.2b) dt V2 is used for the Laplacian operator. The original Babloyantz and Hiernaux scheme allows for different precursors in steps (5.1a) and (5.If). However, in order to derive the Gierer-Meinhardt rate equations (5.2), these precursors must be the same, or related by some constant which can be rolled into A:3 and a. I have treated the precursors in both steps as the same (5). The analysis in the Babloyantz and Hiernaux paper is limited to a study of the Gierer-Meinhardt equations. They did not study the dynamics of the full scheme (5.1). Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 96 In order to fully evaluate the Babloyantz-Hiernaux mechanism, the dynamics of the three species A, I, and P must be solved explicitly. I have found in computations that the dynamics of (5.1) are actually quite different from the dynamics of the Gierer-Meinhardt model. Writing rate equations for (5.1) from mass action gives — = aS + k,PI+k3SAP - k2A - \xA + DAV2A (5.3a) dt — = c'SA2 + k2A - kxPI -vl + D7V2/ (5.3b) dt dP k2A (5.3c) dt Diffusion terms have been added for A and J, but not for P. Whenever dP/dt=0, equations (5.3) become the Gierer-Meinhardt equations. All steady states of (5.3) are therefore steady-state solutions of (5.2); these include the initial spatially homogeneous state and the final patterned state. Babloyantz and Hiernaux made the assumption that dP/dt=0 throughout pattern formation. However, if dP/dt departs from zero between the initial and final state, my computations show that the expected Gierer-Meinhardt pattern may not be readily achieved. For instance in one computation (Figure 5-1), the Gierer-Meinhardt equations (5.2) gave a 4.5 peak pattern in A and /, but a computation using the corresponding parameters in equations (5.3) gave large amplitude pattern only in P. In this case, A and / oscillated and then collapsed to very low values. Both computations are shown at 150,000 iterations. This is far past the point at which the Gierer-Meinhardt pattern quit changing (about 20,000 iterations). While it is uncertain why equations (5.3) lead to the pattern in Figure 5-lb, especially the collapse of A and / concentrations, it is certain that these equations do not produce Gierer-Meinhardt pattern on the same time scale as equations (5.2). The steady states can readily be shown as identical in the two models, but the Babloyantz-Hiernaux scheme cannot be expected to give the proper pattern on the same time scale. The transient dynamics are not the same. Therefore, their scheme does not serve as a mechanistic basis for the Gierer-Meinhardt equations. Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 97 i X5 • - - A/I /5 x5 — - - P /5 b 6 I I 4 -2 -0 Distance Figure 5-1: Comparison of the Babloyantz and Hieraaux (1975) scheme with the original Gierer-Meinhardt model. These patterns are shown at 150,000 iterations, which is far past when the Gierer-Meinhardt model achieves steady state (about 20,000 iterations), a) computed Gierer-Meinhardt pattern with the two morphogens, A and /, shown, as well as their ratio. Parameters are: Af=0.005; 0=0.0006; n=0.21; v=0.27; k3 (c in previous chapter)=0.05; c'=0.025; 5=5=0.4; Z)=lxl0'4; JD/=2xl0"3. b) If dP/dt is not assumed to be constantly zero, the Babloyantz-Hiernaux mechanism does not approximate the Gierer-Meinhardt model. The A and / peaks oscillate and then drop to very low concentrations (approximately one-thousandth of the initial values). P remains large. For this computation, &X=A:2=1. These rate constants are much faster than the others in the model (see above). Even so, dP/dt * 0 during peak firing, and the transient pattern does not match the Gierer-Meinhardt pattern in a). Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 98 5.2 The Enzyme Kinetic Scheme Since the Babloyantz-Hiernaux mechanism does not reproduce the Gierer-Meinhardt dynamics, and also for a better correlation with possible biochemistry, we were led to devise a different chemical scheme (Holloway et al., 1994). In this approach, we followed Gierer and Meinhardt's (1972, p. 36) suggestion that the A2/I term arises from direct inhibition of autocatalysis by the inhibitor molecule. For this, we have proposed a mechanism incorporating competitive inhibition in a Michaelis-Menten scheme which produces the activator A. The Michaelis-Menten scheme involves the binding of an enzyme E to a substrate S. The enzyme can also be bound by an inhibitor /, which precludes substrate binding. In simple Michaelis-Menten kinetics, the ES complex decays to give the product, A. This, however, does not give autocatalysis of A. To get this, A must aid in the creation of the ES complex. We have incorporated steps of allosteric activation of the enzyme by the activator A such that it is an EA2S complex which decays to give the product A. Combined with competitive inhibition of the enzyme by the inhibitor /, this leads to a reaction term much like the A2/I found in the Gierer-Meinhardt equations. To give the rest of the Gierer-Meinhardt equations, we have used simple steps such as equations (5.1a-d). In the biological application of Part Two, there is evidence for a fast-diffusing rescuer of heart function. This entity is stable for several days in experimental conditions in culture. We have incorporated this molecule, R, into the scheme as an additional allosteric regulator necessary for activator production1. Thus, the overall model is an enzyme kinetic scheme which involves competitive inhibition and allosteric activation of the enzyme by the product A (giving a form of autocatalysis) and Such multiple allosteric regulation is well-known for the membrane-spanning proteins of the signal transduction pathway, one of our suggestions for the nature of enzyme E; see Chapter Eight. This type of allosteric regulation has been developed theoretically for some time, for instance see Segel (1975). Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 99 the rescuer R. The Turing morphogen pair are the activator A and the inhibitor /. The enzyme E and rescuer R are necessary for activator production. These species are treated as variables in Part Two. The concentration distribution of E will be used to reflect the effects of the underlying tissue on the heart forming tissue. The distribution of the rescuer R will be used to model rescue experiments of mutant embryos. S and S' are substrates for A and / production, respectively. S and S' are taken to be constant. The mechanism is as follows: E + I^-^-*EI (5.4a) E + 2A*-S^*EA2 (5.4b) EA2 + R «-£-»• EAfi (5.4c) S + EA2R —^-* EA^^R + A (5.4d) 2A + S'—C—*I + 2A (5.4e) S—*->A (5.4f) / —:L-* destruction (5.4g) A —•*—»• destruction. (5.4h) It is assumed that the bimolecular activation pathway is preferred over any other side reactions for enzyme catalyzed activator production, such as unimolecular or zero order activation. This mechanism is intended to capture the essence of the dynamics. It may do so whether or not it correctly represents the types of chemical substances involved. These equations may be abbreviations for a biochemical reaction scheme with more steps. In order to derive the following reaction-diffusion equations from the above mechanism, it is assumed that all equilibrium constants, except for K±, are small. That is, the enzyme is not saturated with A or R, but does bind significantly with /. This means that only a constant term and the term from competitive inhibition are significant in the denominator of the velocity expression for enzyme catalyzed A production (the autocatalytic/inhibitory term of the reaction-diffusion equations). Writing kinetic equations from mass action considerations and adding Fickian diffusion, the following reaction-diffusion equations are derived: Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 100 dA _ kSERA2 . _ _ 2 . / r r x — ° ^ + 1 — — ^ + ^ V 2 A (5.5a) or 1 + Aji — = c'S'A2 - vl + £>,V2/. (5.5b) dt where £=^2 ^3 ^4 a n ^ ^ *s t n e t o t a l enzyme concentration. These equations are quite similar to the Gierer-Meinhardt equations. But, since these equations were derived from a chemical mechanism with inhibition, the denominator of the autocatalytic/inhibitory term in equation (5.5a) is 1+KiI instead of just /. This eliminates the unrealistic property of the Gierer-Meinhardt equations in which the production of A goes to infinity if/ goes to zero. Equations (5.5) are used for all modelling in Part Two. These equations do not exhibit the inaccuracies or instabilities found with the Babloyantz-Hiernaux mechanism. The derivation of the autocatalytic/inhibitory term in equation (5.5a) from steps 5.4a-d in the mechanism maintains the dynamics of the Gierer-Meinhardt equations. By establishing this chemical mechanism, insight into the reactions which give rise to localized pattern can be gained. In the previous chapter, the dynamical conditions for localized pattern were explored. In the present chapter, these conditions are linked with a specific chemical mechanism. This provides the connection between dynamics and chemicals necessary to make predictions for future experimental investigation. For instance, in Part Two diffusivities are predicted for A, I and R from modelling of experimental data with equations (5.5). Experimental discovery of such diffusivities for chemicals involved in heart formation could serve as a test of the model. If such chemicals were isolated, further work could involve direct investigation of the kinetics in mechanism (5.4) and testing of the details of the present model. If experiments do get to such a stage it is likely that many of the steps in mechanism (5.4) would need to be altered in order to fit the data. This would not be a refutation of the present model, however. This mechanism is proposed as a starting point for making the connection between the dynamics and chemistry. As shown in Part Two, the evidence is already strong for the presence of reaction-diffusion dynamics in axolotl heart Chapter 5. A Chemical Basis for the Gierer-Meinhardt Model 101 formation. By examining the characteristics of a specific model, this evidence can be examined in a quantitative way. It is to be hoped that the model will continue to be refined by the experiments to which it can in turn offer predictive guidance. PART TWO BIOLOGICAL APPLICATION "...the arrangement of the words is even more important than the words themselves. And in the same way life is a pattern of chemical processes." - J.B.S. Haldane (1949, pp. 58-62) In Part Two of this work, the preceding theoretical considerations are applied to a biological system. These chapters represent a collaboration with John Armstrong of the University of Ottawa and his co-workers Steven Smith and Heather Easton. The Armstrong laboratory has been conducting experiments on amphibian heart formation for many years. Our group was approached to conduct detailed modelling of experimental data which indicated that a reaction-diffusion mechanism might be present in heart formation. Chapter Six summarizes the experimental data on heart formation in amphibians. Chapters Seven and Eight contain the results of the numerical modelling of heart formation. Many biological systems rely on the ability to accurately produce a localized structure for normal development. Prime examples of this can be found in vertebrate organogenesis. In these cases, it is of extreme importance for the developing embryo that single organs form, with no periodic repeats. In our collaboration with the Armstrong lab, we have applied the reaction-diffusion theory of localized structures to the problem of heart development. With the theory we seek to explain how heart tissue can differentiate in a spatially controlled way out of a much larger field of identical tissue. Amphibian heart development involves both the 103 morphogenesis of a tube which later forms the chambered heart, and the cellular differentiation of the myocardium which leads to beating tissue. In normal development, the morphogenetic region and the differentiated region must eventually coincide. While the application of the theory has implications for vertebrate heart formation in general, the detailed modelling is conducted to quantitatively match the experimental data from the salamander Ambystoma mexicanum (Mexican salamander; axolotl; ajolote). Chapter Six Experimental Work on Axolotl Heart Formation 6.1 Introduction This chapter will summarize the experimental work on vertebrate heart development, with emphasis on amphibians. First, some of the concepts of classical heart induction will be introduced. This will lead into the evidence for post-inductive events which are necessary for heart differentiation. The existing data for the chemicals which might be involved in these events will be reviewed. Finally, evidence for a reaction-diffusion mechanism for post-inductive differentiation will be presented. 6.2 Induction The post-gastrula (neurulation stage) embryo consists of three major tissue layers: the endoderm, the mesoderm and the ectoderm (Figure 6-1). The endoderm is innermost and gives rise to such tissues as the gut and liver. The mesoderm gives rise to muscles and parts of the kidney and heart. The ectoderm is outermost and gives rise to neural tissue and epidermis, among other tissues. The present study is concerned with the endoderm and mesoderm in the axolotl between developmental stage 14 (58 hours of development: all times from Bordzilovskaya et al., 1979; 1989), at the start of neurulation, and developmental stage 35 (120 hours of development), at the commencement of heartbeat. From stage 14 to stage 29 (97 hours of development) the embryo is shrinking in diameter and elongating such that the mesodermal sheets advance ventrally over the underlying endoderm (Figures 6-1, 6-2) without increasing in dorsoventral dimension. At stage 29, the mesoderm primordia meet Chapter 6. Experimental Work on Axolotl Heart Formation 105 Figure 6-1: A cross-section through the axolotl embryo at an anterior thoracic level. The three tissue layers of endoderm, mesoderm and ectoderm are shown at three developmental stages. A) stage 20; B) stage 24; C) stage 28. The shaded areas of the mesoderm represent the extent of the region of specification at that stage. Dorsoventral lengths of the specified region are: 0.28mm at stage 20; 0.38mm at stage 24; and 0.50mm at stage 28. v: ventral midline. From Easton et al (1994). Chapter 6. Experimental Work on Axolotl Heart Formation 106 and begin to fuse. Between stages 29 and 35, a single tube heart forms from the two primordia. The heart chambers form after commencement of heartbeat in normal development. The early tube heart beats peristaltically to drive blood circulation. It was discovered quite early that mesoderm has a strong potential to form heart, which is not easily disturbed. By mechanically preventing mesoderm union at the ventral midline, two hearts can be formed (known as cardia bifida). This was first demonstrated in chicks by Graper (1907). Ekman (1921) and Copenhaver (1926) showed that removal of lateral portions of the mesoderm did not inhibit heart formation. They also both observed that two hearts formed if a side of mesoderm was split lengthwise. Copenhaver (1926) further showed that removal of the anterior portion of one mesoderm side did not inhibit heart formation. A central mechanism for initiation of cell differentiation is inter-tissue signalling. Such signalling was found by Jacobson (1960; 1961) for heart formation in the salamander Taricha torosa. He found that removing the endoderm from an embryo suppressed heart formation in the mesoderm if the embryo was young enough. If the embryo was older than developmental stages 24-26 (mid-tailbud), then heart formation was not suppressed. Jacobson and Duncan (1968) found this effect in vitro, as well. Mesoderm primordia explanted at an early stage (mid-neurula) displayed a functional heartbeat only if cultured with pharyngeal (anterior) endoderm. This action of the endoderm to cause heart formation in the mesoderm was termed induction. Induction is a gradual process in which the mesoderm becomes increasingly fated to form heart as time in contact with endoderm progresses. Endoderm was found to be most inductive in the neurula stages for T. torosa. In the axolotl, inductive capacity in the endoderm is relatively constant through neurulation and tail-bud stages (Smith and Armstrong, 1990). Induction can be delayed by the presence of inhibitory tissues, such as the neural plate and neural fold inhibition reported in T. torosa (Jacobson and Duncan, 1968). Although such inhibitory actions have not been observed in Chapter 6. Experimental Work on Axolotl Heart Formation 107 Figure 6-2: Views of the developing axolotl embryo with the ectoderm peeled away, a-d) side views (anterior on the right) of stage 14, 16, 18, and 20 embryos, respectively, showing the advance of the mesoderm over the underlying endoderm. e) is a ventral view at stage 29, the stage at which the mesoderm primordia meet. NF: neural fold; NP: neural plate; HFM: heart forming mesoderm; PE: pharyngeal endoderm; NT: neural tube; MM: mandibular mesoderm; MA: mandibular arches; GB: gill bulges. From Smith (1990). Chapter 6. Experimental Work on Axolotl Heart Formation 108 other amphibian species, it is important to note that the timing of induction can vary widely between species. Induction is complete by the mid-tailbud stages 24-26 in T. torosa. In the axolotl, induction occurs during neurulation (Smith and Armstrong, 1990; Muslin and Williams, 1991). Smith and Armstrong (1990) found that mesoderm at the start of neurulation (stage 14) was unable to form beating tissue if explanted. Mesoderm explanted at the end of neurulation (stage 20, 70 hours of development) was as capable of forming beating tissue as if it had been left in the embryo. Therefore, in the axolotl, induction is taken to start at stage 14 and to be complete by stage 20. By contrast with these two salamanders, both urodeles, specification of cardiac mesoderm occurs much earlier for the South African clawed frog Xenopus laevis, an anuran. Sater and Jacobson (1989) found that mesoderm explanted at the end of gastrulation (the very early stages in which mesoderm is first formed) was capable of forming beating heart tissue. Even extirpation of the superficial pharyngeal endoderm at the beginning of gastrulation did not affect the ability of the mesoderm to form beating tissue upon explanation. Similar results have recently been reported in quail heart development (Antin et al, 1994). Again, the earliest available mesoderm is capable of forming beating tissue upon explantation. In this case, it was found that endoderm has an important role in providing a substrate for heart development, but it is not necessary as in urodeles. Interestingly, it has also been recently shown that endoderm is necessary for cardiac mesoderm differentiation in chicks (Sugi and Lough, 1994). Therefore, in the species studied, there seem to be two distinct pathways for heart development. The first, which is used by the urodeles and chicks, requires an inducing signal from the underlying endoderm to the mesoderm in order for cardiac differentiation to occur. This induction occurs after gastrulation. The second mechanism, used by Xenopus and quails, involves a much earlier specification of cardiac mesoderm. These species also go through the other stages of heart development more rapidly. Since many vertebrates, including humans, also progress rapidly through heart development, it has been postulated Chapter 6. Experimental Work on Axolotl Heart Formation 109 that the faster pathway, with early cardiac specification, is more typical for higher vertebrates (Antin et al, 1994). However, the slower pathway allows for distinct study of cardiac induction and of post-inductive events. The fast pathway is a much more difficult system in which to separate the events of heart formation. This thesis is concerned specifically with the problem of axolotl heart formation, in which the inductive and post-inductive events occur over a matter of several days, and are therefore quite distinct. The following discussion will centre on species that undergo a distinct induction. Induction appears to have a definite spatial distribution, as well as the temporal one. Fullilove (1970) reported a spatial gradient in inductive capacity of the endoderm in T. torosa. Pharyngeal (anterior) endoderm was found to have high inductive activity, whereas posterior endoderm displayed no inductive capacity. In the axolotl, Smith and Armstrong (1990) found a more localized area of inductive capacity than the previous work had indicated. They found inductive capacity to be localized in the ventral half of the pharyngeal endoderm. Culturing un-induced (stage 14) mesoderms with different portions of endoderm they report a Heart Determination Coefficient (HDC; the percent of cultures beating multiplied by 3, the time needed for normal onset of heartbeat, divided by the days to beat; Jacobson and Duncan, 1968), of 63 for mid-ventral pharyngeal endoderm, an HDC of 38 for mid-lateral endoderm and an HDC of approximately 19 for all other regions tested (including dorsal anterior endoderm and posterior endoderm). The more detailed analysis of Smith and Armstrong (1990) in combination with the general results of Fullilove (1970) indicate that there is a gradient of inductive capacity in the endoderm during amphibian heart development. All these data, however, relate both to the distribution of the inducing signal and the kinetics of its transfer to the mesoderm. Modelling of these two factors is described in Chapter Seven. Chapter 6. Experimental Work on Axolotl Heart Formation 110 There are indications that induction alone cannot be responsible for the final pattern of heart differentiation in the embryo. Specification maps have been obtained both for Xenopus (Sater and Jacobson, 1990) and for the axolotl (Easton et al, 1994). Mesoderm is said to be specified for heart formation if it achieves heartbeat after explantation (Slack, 1983). The specified region is that mesoderm which has been induced. By explanting small portions of the mesoderm and testing for heartbeat, a spatial map can be made. In both Xenopus and the axolotl, maps were made at different developmental stages in order to track the dynamics of the induced area. Figures 6-1 and 6-3 show the specified region in the axolotl at three post-neurula (tailbud) stages. In three species, the axolotl, Xenopus and T. torosa (Jacobson, 1960), it has been reported that the induced region is larger than the subsequent area of heart tissue differentiation. In the axolotl, the specified region is approximately twice as large as the differentiated region. Given the disparity between the induced region and the final heart region, the question becomes one of finding the mechanism which establishes the size of the heart-forming region. In normal development, only the ventral portions of the induced mesoderm are allowed to form heart. The more lateral portions become inhibited with respect to formation of heart tissue. In Xenopus, Sater and Jacobson (1990) found that the specified region decreased after stage 20, the end of neurulation. However, this is also the stage at which the mesoderm primordia meet in Xenopus. In the axolotl, the mesoderm is still advancing over the underlying inductive endoderm until stage 29. It follows that the specified region increases in the axolotl between the observed stages 20, 24 and 28. The difference in timing is important between these species, since axolotl heart differentiation requires a mechanism to form the smaller differentiated region in spite of a growing region of specification. Xenopus heart development could involve the same mechanism as the axolotl (the differentiated region is still smaller than the specified region), but in Xenopus the specified region is decreasing after neurulation. The decrease in the specified region observed in Xenopus cannot account for the delineation of the differentiated region in the Chapter 6. Experimental Work on Axolotl Heart Formation Figure 6-3: Side views of the axolotl embryo at stages 20, 24 and 28. The pharyngeal cavity is shown with the broken line. The specified region in the mesoderm is shaded and outlined by the solid line. The ventral portion of the mesoderm is the bottom edge of the solid line. From Easton et al (1994). Chapter 6. Experimental Work on Axolotl Heart Formation 112 axolotl. A further mechanism is required for heart tissue differentiation. Although these in vitro maps of the specified region give very detailed data for the mesoderm which has been induced, similar in vivo assays were conducted much earlier. Copenhaver (1926) removed ventral portions of mesoderm sides from stage 27 Ambystoma punctatum embryos, and determined a spatial limit to the mesoderm that could be removed and have the embryo still initiate heartbeat (Figure 6-4). He termed the mesoderm ventral to this limit the "heart-field mesoderm". Copenhaver found that if mesoderm migration was blocked by filling the wound with non-heart forming mesoderm, that a smaller region could be extirpated and still prevent heartbeat. This smaller region corresponds to the region of specification determined through in vitro methods. The portions of Copenhaver's heart-field which are dorsal to the region of specification represent that mesoderm which migrates over the inductive endoderm after surgery and subsequently forms heart tissue . 6.3 The cardiac lethal mutation in the axolotl In the axolotl, a mutant has been isolated which displays the phenotype of no heartbeat (Humphrey, 1968; 1972). This cardiac lethal (c) mutation displays itself at stage 35, at the normal onset of heartbeat. Embryos homozygous for the mutation (c/c) have morphologically normal tube hearts (Humphrey, 1972; Lemanski, 1973; Fransen and Lemanski, 1988), but there is no onset of heartbeat. The embryos live for some time after stage 35, but eventually die. C/c embryos swim normally, which indicates that skeletal muscle is functional. Further characterization of the mutant determined that the cardiac muscle proteins do not organize into sarcomeres. This was seen by traditional microscopy (Lemanski and Fuldner, 1977; Hill and Lemanski, 1979), and more recently by confocal microscopy of tropomyosin Chapter 6. Experimental Work on Axolotl Heart Formation 113 C Presumptive Presumptive Endocardium Myocardium Figure 6-4: Ambystoma punctatum embryo at Harrison (1969) stage 27, from Copenhaver (1926). C) shows a cross-section through the ventral portion of the embryo at the level of the mesoderm primordia (level shown in B). If the segment ZW is removed and the space filled with non-heart forming tissue, no heartbeat arises. If ZW is removed and the space not filled, then heartbeat arises. If XY is removed, there is no onset of heartbeat, even if the space is left unfilled. ZW corresponds to the specified region recently determined by in vitro methods (Easton et al., 1994). Chapter 6. Experimental Work on Axolotl Heart Formation 114 chelate fluorescence (LaFrance et al., 1994). Confocal microscopy reveals organized tropomyosin distribution in sarcomeres for wild-type hearts, and amorphous distributions in c/c hearts. It appears that the c/c hearts contain the usual muscle proteins at normal or near-normal levels. The proteins that have been characterized are: actin, myosin, and a-actinin (Lemanski and Fuldner, 1977; Lemanski et al, 1975; Lemanski, 1978; Starr et al, 1989); troponin-T (Fuldner et al, 1984); desmin (Shen and Lemanski, 1989); vimentin and vinculin (Shen and Lemanski, 1986). Starr et al (1989) found decreased levels of tropomyosin in post-heartbeat c/c axolotls (stage 39, 98 hours after onset of heartbeat in normal development, and stage 41, 143 hours after normal onset of heartbeat). Although chemical concentrations can change drastically over time in development, these data could be reflecting pre-heartbeat stage concentrations of tropomyosin. If so, these data could indicate that tropomyosin is directly affected by the c mutation. Tropomyosin chain folding has been identified as one of the first steps in the myofibrillar assembly process (Epstein, 1990). If tropomyosin levels were reduced during heart forming stages, sarcomere assembly would be diminished. LaFrance et al (1994) corroborated the low levels of tropomyosin in post-heartbeat stage embryos with confocal microscopy. Humphrey (1972) postulated that the c mutation could be either of two things. First, it could be an inductive failure, that is the mutation is expressed in the endoderm as a faulty inducing signal. Second, the mutation could be located in the mesoderm, that is the presumptive heart tissue might not respond properly to the inductive signal. Humphrey postulated that this mesodermal failure could arise from an excessively inhibitory environment to heart differentiation in the mesoderm of cfc embryos. The debate around this question has led to much clarification of the processes of heart induction and post-inductive heart differentiation. Chapter 6. Experimental Work on Axolotl Heart Formation 115 The early work on the c mutant indicated that an inductive failure might be occurring. Lemanski et al. (1977) offered microscopic evidence that c/c endoderm contained more microvilli than normal endoderm. Since normal tissue gains more microvilli as development progresses, it was postulated that the c mutant embryos were progressing too fast through the heart formation stages. More direct evidence was reported by Lemanski et al. (1979) by demonstrating that stage 35 c/c hearts could be rescued in culture by wild-type stage 30 endoderm. Hill and Lemanski (1979) found that explanted stage 35-40 wild-type hearts continued to beat, whereas the same stage c/c hearts never began to beat. These results were interpreted as showing that there was no mesodermal failure, as c/c hearts did not beat even when removed from the inhibitory environment implied by Humphrey's (1972) mesodermal failure hypothesis. Kulikowski and Manasek (1978) provided experimental evidence which directly contradicted the results of Lemanski et al. (1979). They found that the explanation of c/c (stage 35 or later) hearts into culture led to heartbeat in as short a time as 30 minutes. It was suggested (Justus, 1978) that these contradictory results could be due to two different strains of the c mutant (Justus, 1978). However, it is also quite possible that experimental conditions were inconsistent between the two groups. Smith and Armstrong (1990) found no heartbeat in explanted stage 20 c/c mesoderms (one culture contained a few beating cells). The rescue of late stage c/c embryos by endoderm was further refined by Davis and Lemanski (1983; 1987). They found that an RNA extract of wild-type endoderm-conditioned medium had the ability to rescue the c mutant hearts. Digestion with ribonuclease (RNase) eliminated the rescuing ability of the medium. LaFrance et al. (1994) argue that this series of experiments demonstrate that the c gene acts by causing abnormalities in the endoderm. While induction had been characterized in T. torosa by Jacobson and Duncan (1968), it was unknown what the timing of axolotl induction was until quite recently. In the tests for inductive failure, it was assumed that the axolotl developed at a similar rate to T. torosa. Two new studies (Sater and Jacobson, 1989; Smith and Armstrong, 1990) show that heart Chapter 6. Experimental Work on Axolotl Heart Formation 116 specification occurs much earlier in other amphibians. With the discovery that axolotl induction is complete by stage 20 (Smith and Armstrong, 1990), any claims that the stage 35 c/c heart rescues are linked to induction become suspect. These rescues performed some 50 hours after the completion of induction are most probably a manifestation of a post-inductive differentiation. In fact, Smith and Armstrong tested the RNA preparations on un-induced hearts (stage 14, both wild-type and c/c) and found no inductive capacity. However, they did corroborate the stage 35 rescues of c/c hearts with the same RNA preparations. In addition, no specific inductive effect was found from culturing endoderm tissue homogenates with un-induced mesoderms or by culturing un-induced mesoderms in endoderm conditioned media (Smith, 1990). Smith and Armstrong (1991) offer more direct evidence against the inductive failure hypothesis. They found c mutant endoderm to be fully inductive. They cultured un-induced (stage 14) wild-type heart mesoderm with mid-ventral pharyngeal endoderm from both c/c and wild-type embryos. In both cases, heartbeat resulted, with no difference in the mean time to beat. To be able to directly study these early developmental stages, Smith and Armstrong used a new technique for identifying c/c embryos. Since there are two sides to the mesoderm prior to stage 29, only one side is explanted for any given experiment. The donor embryo is then closed and allowed to heal. At stage 35, heartbeat commences in the donor if it is of wild-type phenotype (wt/wt or wt/c), and does not if it is a c/c embryo. This technique bypasses the limitations of waiting for the embryo to display its phenotype before conducting the experiment. Finally, a number of different experiments have shown that juxtaposition of wild-type and c/c mesoderm can rescue the mutant tissue, leading to a beating heart. The earliest transplant experiments were reported by Humphrey (1972). In these, the heart-forming mesoderm (the most ventral portion, approximately 25% of the mesoderm) was taken out of a donor embryo Chapter 6. Experimental Work on Axolotl Heart Formation 117 at stage 29 and placed in the same position in a host. Experiments were done with both c/c hosts and wild-type donors, and with wild-type hosts and c/c donors. In the case of c/c donor mesoderm into wild-type host, beating hearts were observed. Therefore, the c/c mesoderm is rescued by the wild-type host. Since the endoderm is unaffected by the c mutation, Humphrey's transplant is evidence for a mesodermal rescue of the mutant tissue. In the reverse experiment, with wild-type heart-forming mesoderm placed into a c/c mutant host, no heartbeat was observed. The mutant host somehow swamps out the viability of the wild-type tissue. Smith and Armstrong (1993) carried out the Humphrey transplants with wider tissue portions. They found that wide (approximately 70% of the mesoderm) wild-type transplants into c/c hosts resulted in beating hearts. This was opposite from Humphrey's results with the smaller portion of donor tissue. Therefore, there is some size effect in these transplants. Presumably the wild-type tissue in these cases is a large enough proportion of the post-surgery mesoderm that it does not become swamped by the c/c tissue. A model of post-inductive dynamics must account for this size effect. Bilateral transplants were conducted on embryos at the very end of induction (stage 20; Smith and Armstrong, 1993). Approximately 70% of each side of the mesoderm of known wild-type embryos was placed into wild-type and c/c hosts. As in the wide stage 29 transplants, the c/c host mesoderm was not enough to swamp the activity of the wild-type mesoderm, and heartbeats resulted. Heartbeat was observed in similar proportions for the wild-type and c/c mutant hosts. Unilateral transplants were also conducted at stage 20 (Smith and Armstrong, 1993). Single sides of mesoderm (approximately 70% portions of total mesoderm side) were explanted from the hosts and replaced by donor tissue from the same region of mesoderm. This Chapter 6. Experimental Work on Axolotl Heart Formation 118 unilateral transplantation allowed for the use of both wild-type and c/c donors. Donor embryos were identified by determining whether the remaining mesoderm side commenced heartbeat or not. Host phenotypes were identified by observing whether or not the explanted mesoderm side beat. It was observed that both c/c donors into wild-type hosts and wild-type donors into c/c hosts developed beating hearts. Of these experiments, the most severe test upon the viability of the wild-type mesoderm is when it is placed into the c/c host. The proportion of wild-type to c/c tissue is approximately 1:3. However, the transplants do rescue in this case. Finally, Smith and Armstrong (1991) demonstrated that mesoderm alone is indeed sufficient for rescue. Single mesoderm sides out of wild-type and c/c embryos were juxtaposed in vitro. Surgery was again conducted at stage 20. The c/c tissue formed a beat within 6.3 days after surgery, on average. Two beating foci formed in these experiments, with the original wild-type tissue beating at 4.0 days on average. This difference in timing represents some type of delayed rescue. A simple mechanism to account for this would be the diffusion of a rescuing molecule between the tissues. As well, the tissue size effects observed in the stage 29 transplants suggest a diffusible rescuer. When the proportion of wild-type to mutant tissue is too low, the rescuer presumably diffuses to too low a level to affect a rescue. 6.4 Experimental Characterization of Heart-Forming Chemicals In chicks, an RNA extract from heart tissue was shown to increase beating tissue in explants of presumptive heart mesoderm (Butros, 1965). RNA extracts from liver (Butros, 1965), brain, kidney and thymus (Niu and Deshpande, 1973) did not enhance heartbeat in presumptive heart tissue explants. It was subsequently found that the active species in the heart RNA extracts was a 7S polyadenylated RNA (Niu and Deshpande, 1973; Deshpande and Siddiqui, 1977; Deshpande et al., 1977). More recent work on this RNA, however, has Chapter 6. Experimental Work on Axolotl Heart Formation 119 shown it to be 99% homologous with a portion of the chick mitochondrial ATPase 6 gene (Desjardins et al., 1989). This means that the 7S RNA could be a degradation product of mitochondrial RNA that is produced during the RNA isolation (Smith, 1990). It is also worth noting that the 7S RNA was isolated from mature heart tissue, which is not an inductive tissue in normal heart development. However, the fact that the heart RNA extracts do contribute to enhanced heartbeat in mesoderm explants suggests that there is some active species present along with the (apparently) mitochondrial 7S RNA in the overall extract. The enhanced heartbeat stimulated by the active species could stem from an association of the active species either with the inductive pathway or with post-inductive differentiation events. As heartbeat exists without the application of the active species, it is more likely that this species is associated with post-inductive events. Much of the evidence presented for the inductive failure hypothesis in c mutant axolotls contained information on the chemical nature of the species involved (Davis and Lemanski, 1987; LaFrance et al., 1994). While it has subsequently been shown that these rescues cannot pertain to inductive processes, the characterization of a post-inductive rescuer is very important. Smith and Armstrong (1991; 1993) have provided evidence for a mesodermal post-inductive rescuer. Combined with the stage 35 rescue of c/c hearts with RNA extract (Davis and Lemanski, 1987; Smith and Armstrong, 1991), there is strong evidence for a post-inductive organizer (or organizers) of cardiac myofibrils which is absent in the c mutant. LaFrance et al. (1994) present evidence for a rescuer in both mesoderm and endoderm. Combination endoderm/mesoderm portions were extracted for RNA to give higher yields in the RNA isolation process. Pure mesoderm extracts were also studied in order to separate mesoderm effects from endoderm effects. Pure mesoderm, pure endoderm and mesoderm/endoderm RNA extracts were shown to rescue stage 35 c/c embryos. Digestion with RNase destroyed the rescuing ability of the mesoderm/endoderm extracts. Lemanski and Davis (1987) found the same in pure endoderm. Treatment of the Chapter 6. Experimental Work on Axolotl Heart Formation 120 mesoderm/endoderm extracts with Proteinase K decreased the rescuing ability, but did not destroy it. In conjunction with the RNase experiments, this argues against the presence of protein rescuer(s), and for the presence of RNA rescuer(s) in both mesoderm and endoderm. However, it is also conceivable that Proteinase K treatment destroys a protein rescuer, but that the protein is quickly replenished from its RNA. A dose response was found in that the rescuing ability of the extracts varied with amount of RNA in the extracts. It is unknown whether the mesodermal rescuer reported by LaFrance et al. (1994) is the same as that indicated by Smith and Armstrong (1991; 1993), or whether this is the same rescuer as that found in the endoderm. Work is in progress in Lemanski's lab to further characterize the RNA extracts from endoderm. Erginel-Unaltuna et al. (1993) have created a complementary DNA library from the RNA extracts and isolated 20 clones. Of these, 13 show homology with ribosomal RNAs. An antibody was made against a protein from one of the reading portions found in the RNA. Immunofluorescent studies show spatial distribution of a 200 kD protein which starts at stage 24 and persists into adult tissue. It is unknown whether this protein is actually involved in heart formation, but is certainly a step towards characterizing the active species in the rescuing extracts. In concurrent studies (LaFrance et al., 1993), the RNA extract was separated on a sucrose column. Fractions from the upper end of the sucrose gradient, mostly small RNAs, had the highest rescuing activity. Various peptide growth factors have also been tested for effect on axolotl heart formation. Muslin and Williams (1991) found a significant increase in heart determination for stage 14 wild-type mesoderm cultured with Transforming Growth Factor-pl (TGF-|31), or Platelet Derived Growth Factor (PDGF). Mesoderms cultured with TGF-pi were almost as likely to form beating tissue as mesoderms cultured with endoderm. No significant increase in heart determination was found by culturing stage 14 mesoderm with Epidermal Growth Factor Chapter 6. Experimental Work on Axolotl Heart Formation 121 (EGF), Fibroblast Growth Factor (FGF, actually found to have an inhibitory effect on later stage mesoderms), insulin, or retinoic acid. Mangiacapra et al. (1993) found that stage 35 c/c hearts could not be rescued by application of TGF-fJl or Activin A. They also found that TGF-J31, Activin A and TGF-(35 could not induce stage 14 c/c mesoderm, but that these growth factors could cause myocardial differentiation in stage 14 wild-type mesoderms. These results contrast with the lack of heartbeat stimulation found for endodermal RNA extracts cultured with stage 14 wild-type mesoderms (Smith and Armstrong, 1990). The current evidence indicates that small RNAs are responsible for the late stage rescue of c/c hearts and are hence associated with post-inductive heart differentiation. Peptide growth factors (TGF-pi and PDGF) can greatly increase the likelihood of heartbeat in stage 14 mesoderm explants, and are hence implicated in the induction process. 6.5 Conclusion At this time, the conclusive evidence offered by Smith and Armstrong (1990; 1991; 1993) supports the hypothesis that c mutant mesoderm cannot respond properly to a fully functional inductive signal from the endoderm. The inductive failure hypothesis is still supported by some (LaFrance et al., 1994) based partly on mis-representation of Smith and Armstrong's results and partly on disagreement over the concept of induction. To conclude this chapter, a picture of axolotl heart formation is offered which attempts to synthesize the above data on induction and post-inductive events. It appears that endodermal induction is fully functional in c mutant axolotls (Smith and Armstrong, 1991). Induction occurs during neurulation, starting at stage 14 and being complete by stage 20. Between stage 14 and stage 29, the mesoderm advances over the underlying inductive endoderm, resulting in a growing region of specification (induced mesoderm) during this period (measured at stages 20, 24, and 28; Easton et al., 1994). The Chapter 6. Experimental Work on Axolotl Heart Formation 122 long period between the end of induction and the onset of heartbeat in the axolotl has allowed for detailed study of post-inductive processes. Evidence has been presented for the post-inductive differentiation which must occur before onset of heartbeat. First, the area of specification is larger than the area of final heart tissue differentiation (Jacobson, 1960; Sater and Jacobson, 1990; Easton et al., 1994). Second, stage 35 c/c hearts can be rescued and form heartbeat (Lemanski et al., 1979; Smith and Armstrong, 1990). Finally, juxtaposition of wild-type mesoderm with c/c mesoderm in post-inductive stages leads to rescue of the mutant tissue (Humphrey, 1972; Smith and Armstrong, 1991; 1993). The c mutation involves a species necessary for this post-inductive differentiation. Specifically, the c mutant embryos lack a species necessary for post-inductive sarcomere assembly. Application of this species, or one of its precursors, to induced c mutant mesoderm leads to rescue. Humphrey (1972) originally postulated that a mesodermal failure (as opposed to an inductive failure) could arise from production of a substance inhibitory to normal cardiomyocyte formation. In our work, we treat the mutation as a lack of a necessary species. Both postulates are in agreement with the present data, but would involve different parameters in a dynamical model of the phenomena. As the area of differentiation is smaller than the area of specification, some post-inductive mechanism must be present which partitions heart tissue from induced tissue. Armstrong (1989; also Smith and Armstrong, 1991; 1993) has postulated reaction-diffusion dynamics which achieve this spatial partition. Chapters Seven and Eight summarize detailed numerical modelling based upon this hypothesis. It is a hierarchical model, with induction feeding into the post-inductive differentiation process. Induction is treated as providing a precursor necessary for the patterned chemical intermediates of the reaction-diffusion mechanism which lead to final heart differentiation. The species with which stage 35 c/c hearts can be rescued is treated as a separate precursor in the reaction-diffusion mechanism1. Through the The chemical analyses discussed in the previous section can be summed up within this interpretation. Late stage (induced) c/c hearts exposed to the rescuing RNAs are being exposed to the rescuer (R in the next two Chapter 6. Experimental Work on Axolotl Heart Formation 123 detailed modelling, we have developed a framework within which to interpret the existing data on axolotl heart formation. It is our hope that this framework can lead to insight into heart formation which can guide future experiments. As well, we make predictions for experimental tests of the theory. The following two chapters discuss the specific data modelled and the numerical results obtained. chapters). Early stage (uninduced) wild-type mesoderm can be induced by certain peptide growth factors. The inducer (E in the next two chapters) is provided for these hearts, which already have rescuer. Both precursors are present for the mechanism of differentiation. When uninduced c/c mesoderms are exposed to these inducing growth factors, no beat is observed because the rescuer is not present in the mutant tissue. Both precursors are not present and the mechanism of differentiation fails. Chapter Seven Modelling Induction 7.1 Introduction In the next two chapters, a dynamical model of heart formation is developed. An attempt has been made to quantitatively match the experimental data with computational results. In establishing a theoretical model for heart formation in the axolotl, a distinction is made between two stages in the dynamics. The first stage is induction. For this, the data on the region of specification and the mesoderm velocities between stages 13 and 29 has been incorporated. The second stage of the dynamics is post-inductive, in which reaction-diffusion pattern arises within the established specified region. This chapter will discuss the approach to modelling the inductive dynamics. Chapter Eight will address the computation of the second stage of the model, the reaction-diffusion pattern. The present chapter will discuss the dynamics of mesoderm movement, assumptions about the induction process, and an interpretation of the data concerning levels of inductive capacity in the pharyngeal endoderm in terms of these considerations. The computations will also be matched with data from Easton et al. (1994) concerning the region of specification. 7.2 The Inductive Process Induction is the action of the endoderm to cause heart formation in the mesoderm (section 6.2). We take this action to be chemical in nature, as induction can occur both in vivo (Jacobson, 1960; 1961) and in vitro (Jacobson and Duncan, 1968). This makes mechanical action highly unlikely, as the inter-tissue forces are much different in the two cases. Also, inductive action can be reproduced by the application of various peptide growth factors Chapter 7. Modelling Induction 125 (Muslin and Williams, 1991). This is highly indicative of a chemical basis for induction. Given a chemical action of the endoderm on the mesoderm, one can postulate the simple reaction: inducer{endo) <s> product(meso). (7.1) However, since the product of induction tends to remain in the mesoderm (explanted, induced mesoderms, at stage 20, form heartbeat, see Smith and Armstrong 1990), the reverse reaction must be negligible in comparison with the forward reaction. Hence, our simple reaction becomes: inducerleM)-* product(mao). (7.2) This irreversible chemical reaction implies that the concentration of the product in the mesoderm will increase as long as the mesoderm is in contact with the endoderm (assuming other product decay reactions are slow). The concentration of the product in the mesoderm is an integral over time of the concentration of the inducer in the endoderm. This can be expressed as follows: product,^ = fenducer{endo)dt. (7.3) The inducer concentration is assumed to be constant. No attempt has been made to account for a decay in inducer concentration as the reaction proceeds. The inducer is assumed to be replenished continuously and be at a stable steady state concentration. Computation of product concentration would be extremely difficult if both inducer concentration decay and mesoderm movement were modelled. Also, it has been reported that there is no significant difference in the inductive capacity of endoderm between stages 14 and 30 (Smith and Armstrong, 1990). This argues against a fixed starting concentration of inducer which then decays, and for maintenance of a constant concentration in a non-equilibrium steady state. Computation of equation (7.3) results in a concentration gradient for the product of induction. This gradient is quite different from the exponential gradient described in Chapter Four. The exponential gradient arises through dynamics which involve chemical decay for the gradient molecule, whereas equation (7.2) assumes a lack of these processes. Chapter 7. Modelling Induction 126 Computation of equation (7.3) requires values for the inducer concentration. Values were used that reproduced the data from Smith and Armstrong (1990) for an inductive capacity gradient in the endoderm. The levels were originally reported in terms of the HDC of Jacobson and Duncan (1968; see section 6.2). The HDC is: HDC = 3* % beating cultures mean No. of days to beat 3 represents the number of days to beat in normal development. Therefore, if all of the cultures beat and the mean time to beat is the same as in normal development, the HDC is 100%. The % beating term represents the stochastic aspect of the explant experiments. Fluctuations in these experiments arise from a number of sources. One source is thermal fluctuation in chemical concentration. This is modelled as described in Chapter Three. Further sources of fluctuation include variation in mesoderm position with respect to the endoderm for a given time, and variations in geometry and inter-tissue contact in explant experiments. Since we have no model for the size and distribution of fluctuations in these cases, we do not attempt to match the HDC data of Smith and Armstrong (1990) directly. Instead, we model the time to beat for those cultures which do beat. The time to beat data corresponding to the original HDC data are summarized in Table 7-1. Chapter 7. Modelling Induction 127 Table 7-1: Time to Beat for In Vitro Induction (John Armstrong, personal communication) Experiment Days to Beat control (st. 14 meso) 9.5±0.3 st. 14 meso+mid-lateral endo 6.9±0.4 st. 14 meso+m id-ventral endo 5.2±0.3 The portions of endoderm used in the cultures are shown in Figure 7-1. Schematic diagrams of the experiments are shown in Figures 7-3 A, C and E. In the original publication (Smith and Armstrong, 1990) it was reported that the control mesoderms had the same HDC as those cultured with posterior (yolky) or anterior dorsal (roof) endoderm. It is assumed that this relation to holds for the time to beat parameter. In the model, heartbeat is a product of the tissue differentiation due to post-inductive reaction-diffusion dynamics. The product of induction in the mesoderm is treated as a precursor (reactant) in the reaction-diffusion scheme. The concentration of this precursor has a great deal of effect on the timing of peak firing in the reaction-diffusion computation (corresponding to the onset of heartbeat in the embryo). For a given experiment (mesoderm cultured with mid-ventral endoderm, mid-lateral endoderm, or no endoderm), the precursor concentration is found that will give peak firing corresponding to the experimental time to beat. Therefore, the precursor concentration is tightly constrained for a given set of reaction-diffusion parameters. Since the precursor concentration in the mesoderm is an integration of the inducer concentration in the endoderm, a definite inducer concentration distribution is produced by fitting the three different times to beat of the above in vitro induction experiments (Table 7-1). In other words, the average concentration of inducer in the three different tissue portions is a tightly controlled parameter in the modelling. Since these different concentrations come from distinct portions of the embryo, there is a definite spatial Chapter 7. Modelling Induction 128 A mv Figure 7-1: Side view of a stage 14 axolotl embryo. A: anterior, P: posterior, mv: mid-ventral endoderm; ml: mid-lateral endoderm; r: roof endoderm; y: yolky endoderm. From Smith and Armstrong (1990). Chapter 7. Modelling Induction 129 concentration distribution of inducer that follows from the model. These considerations will be combined with in vivo data in the next section into a computed distribution for the inducer concentration. 7.3 Inducer Distribution In this section, I start the quantitative modelling of the timing of induction. Therefore, a short explanation of the computational time-scale is necessary. As discussed in Chapter Three, the dynamic equations are solved forward in time by approximating the infinitesimal derivatives with finite time steps. Each time step in which the equations are solved is an iteration. An arbitrary decision as to the relation between the time step and real time is initially made. For most computations, I use 12,000 iterations = 1 day. The time step used with this ratio is A£=0.075. For some computations the ratio is decreased by a factor of five, so that 2400 iterations = 1 day. In these computations, the time step needs to be increased by a factor of five to Afc=0.375 so that the equations will effectively react and diffuse more during the longer iteration. Once the relation between computational and real time is established, the computational mesoderm velocities must correspond to the experimentally determined movement of the mesoderm. The kinetic parameters in the reaction-diffusion equations which model the experimental times to beat are linked to real time, since they are related to the real time of the mesoderm movement. If the time step is made five times larger, then the mesoderm moves five times faster and reaction and diffusion happen five times faster. Through the relation between the real and computational time-scales, the computed kinetic parameters can be converted into real units. By combining the levels of inducer necessary to model the data in the previous section with data for the mesoderm position at given times, a final inducer distribution can be established. Chapter 7. Modelling Induction 130 The dorsal-ventral distance from the leading edge of the mesoderm to the ventral midline is given for the stages in which mesoderm movement is occurring in Table 7-2. Table 7-2: Position of Leading Edge of Mesoderm (John Armstrong, personal communication) Developmental Stage 13 14 20 24 27 28 29 Distance from Ventral Midline 1.66 mm (extrapolated) 1.52 mm 1.00 mm 0.52 mm 0.25 mm 0.07 mm 0.00 mm Hours of Development1 55 58 70 80 86 92 97 The mesoderm is a constant length of 1.00 mm during these stages. No significant growth or cell division occurs (John Armstrong, personal communication). The movement of the mesoderm over the endoderm is a product of a shrinking in diameter and lengthening of the embryo. This change in shape is apparently due to cell sorting rather than active cell division. The fact that the mesoderm remains a constant length throughout heart formation simplifies the model, since the positions of all parts of the mesoderm are linked to the above data for the position of the leading edge of the mesoderm. Computationally, the constant 1.00 mm length of the mesoderm is represented as a constant 40 points on a line for which the values of the model concentrations are computed. Since the space between grid points corresponds to 25 jxm, the concentration at each point is effectively calculated in a volume of 2xl0-8 ml. ^rom Bordzilovskaya et al. (1979; 1989). Chapter 7. Modelling Induction 131 MESODERM NP | Hf ENDODERM -ECTODERM Figure 7-2: Cross-section of a stage 13 embryo at the level of the mesoderm primordia. NP: neural plate; NF: neural fold; NC: eventual site of notochord. The three layers of the endoderm, mesoderm and ectoderm are shown. LE: leading edge of mesoderm. VML: ventral midline. L. G. Harrison sketch, extrapolated from stage 17 micrograph. Chapter 7. Modelling Induction 132 After gastrulation, at stage 13, the mesoderm is just starting to move ventrally over the endoderm (Figure 7-2). Since stage 13, and the majority of stage 14, explanted mesoderms do not develop a heartbeat, this indicates that these mesoderms have not yet had sufficient contact with inductive endoderm. The endoderm which underlies the mesoderm at stage 14 (that endoderm which is at least 1.52 mm from the ventral midline, see Table 7-2) is taken to have inducer concentrations too low to cause induction. A small percentage of explanted stage 14 mesoderms do develop a heartbeat. We have accounted for this variability in the model by identifying it with the variability in mesoderm position at a given stage. The position of the leading edge of the mesoderm with respect to the ventral midline can vary by as much as ±0.1 mm between embryos at a given stage (John Armstrong, personal communication). In the model, we assume that it is the spatially advanced stage 14 mesoderms which can develop a heartbeat. Therefore, the stage 14 explants (the controls, which beat in 9.5±0.3 days) were modelled by allowing the computational mesoderm to advance beyond the average position at stage 14. It was found that by letting the computational mesoderm advance 0.08 mm past the average position, the experimental time to beat was duplicated. This means that the leading edge of the advanced mesoderms was 1.44 mm from the ventral midline instead of 1.52 mm. In the model, any mesoderm that advances between 0.08 mm and 0.10 mm from the average stage 14 position will achieve heartbeat. This allows for a certain population of the explants to achieve heartbeat in the model, which is in agreement with the 20% of controls which achieve heartbeat experimentally (Smith and Armstrong, 1990). The level of inducer in the endoderm subjacent to that mesoderm position must then be enough to cause firing of the reaction-diffusion peaks at the time corresponding to 9.5 days (Figure 7-3 B). The endoderm gradient was further mapped by modelling the data for the mid-lateral and mid-ventral endoderm cultured with stage 14 mesoderm (Figures 7-3 D and F, respectively). These computations differ from the control in that there is now growth in the concentration Chapter 7. Modelling Induction 133 Figure 7-3: The experiments of Smith and Armstrong (1990) are shown with the corresponding computations. A) shows the explantation of mesoderm at stage 14 (sketch of stage 17 embryo shown). B) shows the concentrations of A and / (xO.l), the Turing morphogens (see Chapter Five), and the product of induction, E (xl5), versus distance at 9.5 days, when heartbeat is observed experimentally. The maximum in the E gradient is 0.023. C) is the co-culturing of stage 14 mesoderm with mid-lateral stage 14 endoderm. D) shows the levels of A, I (xO.l), E (product of induction, xl5) at 6.9 days, corresponding to the experimentally observed onset of heartbeat. The maximum in £ goes from 0.0179 at stage 14 to 0.028 by 6.9 days after surgery. E) is the culturing of stage 14 mesoderm with mid-ventral stage 14 endoderm. F) shows the computed levels of A, I, and E (concentration adjustments as above) at 5.2 days, the experimental onset of heartbeat. The maximum in E goes to 0.031 by 5.2 days post-surgery. Points c and g on the embryo cross-sections correspond to points c and g on the computations. The E values quoted in this legend include further considerations of endoderm to mesoderm transfer considered below (in this section and in sections 7.4, 7.5 and 7.6). From Holloway et al. (1994). Chapter 7. Modelling Induction 134 mid-ventral In 8.31x10A mid-lateral b e \mMMMmmmttwtmmMUMauuumumaKmaMMmMMMmMauwm V / \ ' 0.0 mm Distance 1.44 mm 1.66 mm *> 90° e„ 0° 2.66 mm 146° Figure 7-4: Model for concentration distribution of the endodermal inducer. The positions of the mid-ventral and mid-lateral endoderm explants are shown. Other portions (roof and yolky) correspond to areas of zero or below threshold levels of the inducer. The maximum concentration of inducer is 8.31x10"'* in the concentration scale of the model. The angular co-ordinates used in computing positions of the mesoderm with respect to the endoderm are shown at several locations. Distances to the ventral midline at stage 13 (the beginning of the computation), in millimetres, are shown for those angles. Chapter 7. Modelling Induction 135 of the product of induction. This reflects the fact that the mesoderm remains in contact with the endoderm throughout the computation. The inducer concentration necessary to cause average stage 14 mesoderms cultured with mid-ventral endoderm to beat is the same as that for the advanced stage 14 mesoderms (controls). The endoderm-cultured mesoderms beat earlier because the product of induction is allowed to grow. Computations of stage 14 mesoderm cultured with mid-lateral endoderm require one half of the level of inducer necessary for the mid-ventral cultures. Putting this together with the requirement for no firing of reaction-diffusion peaks in stage 13 and regular stage 14 explants leads to an inducer distribution in the endoderm as shown in Figure 7-4. The concentration steps up sharply between the positions of 1.66 mm and 1.44 mm from the ventral midline, and is then flat for all positions closer to the ventral midline. Compare this profile to the endoderm portions shown in Figure 7-1. Mid-ventral endoderm is a portion of the flat part of the profile. Mid-lateral endoderm is a portion centred over the step in the profile (note the position of the leading edge of the mesoderm in Figure 7-1). Since these experiments include extended portions of endoderm, the inducer level is taken as the average over the tissue sample. All other portions shown in Figure 7-1 are treated as having zero inducer concentration. The fact that the relative positions of the endoderm and mesoderm are now allowed to vary in the model must be incorporated into equation (7.3). Position in the mesoderm, omeso, Can be related to an angular position in the endoderm, 6endo (see Figure 7-5). The assumption that the inducer profile is unchanging in time with respect to the angular coordinate is equivalent to the assumption that embryo narrowing and elongation occurs by cell movements which are exclusively anteroposterior. If angle 0° corresponds to the ventral midline, with 180° as the notochord, then the dorsal boundary of the specified region can be set at 90°. It is found that the dorsal margin of the mesoderm (always 1.0 mm linearly from the leading edge) stays at a relatively constant angle of 146° through the tailbud stages for which specified region Chapter 7. Modelling Induction 136 20 24 28 Figure 7-5: Angular co-ordinates in the endoderm. 0° corresponds to the ventral midline, with 180° corresponding to the dorsal-most point on the embryo. Angles corresponding to mesoderm positions are shown at three stages: 20, 24 and 28. The angle corresponding to the dorsal margin of the mesoderm is 143° at stage 20, 146° at stage 24, and 150° at stage 28, with an average of 146°. The dorsal extent of the specified region is chosen to always be at 90°. This places the angle corresponding to the leading edge of the mesoderm at 65° for stage 20, 52° at stage 24, and 15° at stage 28. Chapter 7. Modelling Induction 137 data is available. Angular velocities, v, can be calculated from the angle of the leading edge of the mesoderm at given times. The dorsal margin of the mesoderm remains fixed, while the rest of the mesoderm is allowed to move over the endoderm in the model. This means that the position of 1.0 mm in the mesoderm always corresponds to the endoderm angle of 146°. The angular position corresponding to the leading edge of the mesoderm (0.0 mm) can be expressed by the following equation: <U, = e 0 -v (7.5) where 90 is the angular position of the leading edge at stage 13 (the beginning of the computation) and t is the time elapsed. The positions for all points in the mesoderm can then be assigned angles according to: 8— = (B0 - « ) + S_[146° - (60 - vt)] (7.6) where Smeso is expressed as a proportion of position over total length (value between 0 and 1). This relation allows the mesoderm positions to be continuously reassigned according to the position of the leading edge (Figure 7-6). Since the inducer in the endoderm is a linear gradient between 60 (1.66 mm from ventral midline at stage 13) and 90° (dorsal boundary of the specified region), the expression for the inducer concentration in this region becomes: inducer(endo) = m(dendo) + yiot = m{(80 - vt) + Smeso[146° - (60 - v*)]} + yM (7.7) where m is the slope of the line and yim is its value at the ventral midline. Substituting this into equation (7.3) gives: product(meso) =^ (m{(9 0 - vt) + Sneso[U6° - (8 0 - vt)} + yixA)dt (7.8) where t = iterations x At in the computation. Integrating this expression gives: Product^ = m[B0(l-Smeso)+ 146° Smeso]t + yiJ + m(S^>t2 ( ? 9 ) The product concentrations for points in the mesoderm corresponding to angles between 90 and 90° in the endoderm are calculated from this expression. Mesodermal product concentrations for positions corresponding to angles less than 90° are calculated from a constant level of inducer (8.31xl0"4). The mesoderm velocities are calculated for the Chapter 7. Modelling Induction 138 B -1.00 mm-In T -M h 0.0 mm [\ 0. ! ! 0.28 50 mm mm 1.00 mm VML 65° 90° ^ A 146° (2.00 mm from VML) \ 0.0 mm 4 In -tr 1 1.00 mm i\ 0.50 mm \ .0° 9 0 ° \ 146° (1.00mm from VML) VML Figure 7-6: Mesoderm movement over the endoderm with a fixed dorsal point. This movement is computed from equation (7.6). The angle of 146° in the endoderm always corresponds to 1.0 mm in the mesoderm. 90° in the endoderm always corresponds to the dorsal boundary of the specified region in the mesoderm. A) stage 14+1.9 hours, the point at which the leading edge of the mesoderm is first specified in normal development; B) stage 20; C) stage 29. The mesoderm distances, both within the mesoderm and with respect to the ventral midline are shown. Some corresponding endoderm angles are shown. The distance from the 146° point in the endoderm to the ventral midline is shown for each stage. 6Q is the angle corresponding to the position of the leading edge of the mesoderm at stage 13. Chapter 7. Modelling Induction 139 intervals available in Table 7-2. These velocities are used in the appropriate stage of the computation in equation (7.9). With equation (7.9), the growth of the product of induction in the mesoderm can be computed from stage 13 to stage 29. For a computation of normal, in vivo development it turns out that the reaction-diffusion peaks do not fire in time to correspond to the onset of heartbeat. This is because the inducer level necessary to fire the reaction-diffusion peaks in 5.2 days (62,400 iterations) for the stage 14 mesoderm+mid-ventral endoderm simulation is lower than the inducer level necessary to fire the reaction-diffusion peaks in the 62 hours (31,000 iterations) of normal development between stages 14 and 35. To account for this, a decreased efficiency of induction is assumed for the in vitro experiments as compared to the in vivo cases. Proper timing of peak firing for both the in vitro experiments and normal development is achieved when the ratio of inductive efficiency between in vivo and in vitro is 40:1. This decrease in efficiency which is necessary for the modelling corresponds to the different physical environments of the in vitro and in vivo cases. Inter-tissue contact is presumably much tighter and evenly distributed in vivo than in vitro. 7.4 Saturation An additional consideration must be made for in vivo induction. The level of inducer encountered by the mesoderm is constant from shortly after stage 14 until stage 29. The product of induction in the mesoderm grows by accumulation (integration). Given the t-squared growth from the integration, the product of induction reaches high levels by the later stages of development. This becomes a problem when modelling stage 29 transplant experiments (Humphrey, 1972). To properly model the experimental data, the product of induction must increase in concentration in early stages, but its growth must plateau off at later stages. If the product concentration is too low at early stages, then early stage (14 and Chapter 7. Modelling Induction 140 20) explant and transplant experiments cannot be modelled. If the concentration is too high at stage 29, the Humphrey transplants will not work. To account for this non-linear growth, we assume a saturation of the induction product in the mesoderm. This saturation is implemented by dividing equation (7.9) by 1 + k(product,meso))2, where k is just a rate constant for the saturation. This slows down the induction as the product concentration increases: the product is self-inhibitory. A product-squared saturation is used to slow growth more effectively at higher product concentrations. Product growth is effectively halted at later stages of development (corresponding to tailbud stages). Several physical causes for this saturation can be postulated. One possibility is that the product molecules undergo a further reaction with the inducer. Another possibility is that the saturation effect is the completion of differentiation of specified cells to express the phenotype represented by the product of induction. Saturation could also arise from the filling of membrane receptor sites. This saturation would not depend on the product concentration squared, however, unless the product is a binding molecule which needs two receptor sites. Inducer concentration decay does not give a good account of this saturation effect. If inducer decay were incorporated into the model, the leading edges of the mesoderm would still be encountering new (initial) levels of the inducer. The product of induction would continue to grow significantly if no saturation effect were incorporated. Furthermore, the evidence for constant inductive capacity in the endoderm over several days (Smith and Armstrong, 1990) argues against a significant decay in the inducer concentration. 7.5 Region of Specification It was hoped that the inclusion of simple chemical induction, the mesoderm velocities and the saturation of the product in the model, as described above, would result in accurate computation of the specified region (area of high induction product) at each stage. However, Chapter 7. Modelling Induction 141 it was found that these considerations alone led to somewhat large computational specified regions (see Table 7-3). Table 7-3: Extent of Specified Region (SR) in Dorso-ventral Direction Stage 20 24 28 SR (expt)2 0.28 mm 0.38 mm 0.50 mm SR (comp): fixed dorsal pt. in meso. 0.30 mm 0.45 mm 0.55 mm SR (comp): dorsal 25% of meso. fixed 0.30 mm 0.43 mm 0.48 mm To decrease the size of the specified region in later stages, an additional consideration is incorporated into the model. Instead of fixing only the most dorsal point of the mesoderm with respect to the endoderm, a region of the dorsal mesoderm is constrained. The closest agreement with the measured areas of specification is when this constrained region is the dorsal 25% of the mesoderm (Table 7-3; Figure 7-7). This means that mesoderm positions corresponding to endoderm between the angles 146° (dorsal limit of mesoderm) and O.250O+1O9.5° (from equation 7.6 with Smeso=0.75 and t=0) do not move. 7.6 The Product of Induction With these dynamics incorporated into the model, we can compute concentration distributions for the product of induction that agree with the experimental data. Computed distributions are shown for four developmental stages in Figure 7-8. The stages shown are 14; 20; 24; and 29. The effects of the saturation can be seen in the flat profile of the product. The lack of product in the dorsal portions of the mesoderm stems from the requirement that the angles 90° to 146° remain unspecified and from the constraint on movement between the mesoderm and the endoderm in the dorsal 25% of the mesoderm. The regions defined by positive induction product concentration in Figure 7-8 correspond to the region of 2From Eastern et al. (1994). Chapter 7. Modelling Induction 142 —0.52 mm — • 0.0 mm 0.38 mm h i 00 mm i , / 1 7 0.75 mm I. s0° V 52° (a) 90° \ ^ \ 146° (point b, 1.52 mm from v) O.250o+109.5° (point e) Figure 7-7: Mesoderm movement over the endoderm with a fixed dorsal region. The fixed region extends from 1.0 mm to 0.75 mm in the mesoderm. A) stage 20; B) stage 24; C) stage 28. No relative movement is allowed between endoderm and mesoderm along segments eb or df. The angles e and f are calculated from equation (7.6) with Smeso=0.75 and t=0. All mesoderm movement is for segments ae and cf. v is the ventral mid-line. D) shows endoderm angles with corresponding mesoderm positions at stage 24. Chapter 7. Modelling Induction 143 St 14 • • • • St 2 0 • • • st 24 " st 29 c o c V o c o o 4 3 2 0 0.0 mm / I 0.50 mm 0.30 mm 0.43 mm Distance 1.00 mm Figure 7-8: Concentration distribution of the mesodermal product of induction at four stages of development. Stages 14; 20; 24; and 29 are shown. This figure shows the effect of the saturation of product growth at later stages of development. The lack of product in the dorsal halves of the computation follows from the discussions of section 7.5. The maximum concentration of £ is 0.057 at stage 20 and 0.085 at stage 29. Chapter 7. Modelling Induction 144 specification in Easton et al. (1994). The concentration gradient of the product of induction is identical for all in vivo computations. This includes all of the transplant experiments and normal development. However, the profile is changed for the in vitro computations. When the epithelial sheet of the mesoderm is explanted into culture, mechanical forces cause it to curl (Figure 7-9 A). The mesoderm loses the flat geometry and the close contact with the endoderm which it has in vivo. To reflect this rolling, the profile of the induction product (assumed to be non-diffusing, or very slow-diffusing) is altered. The maximum in the profile is pulled away from the ventral midline, reflecting that if the most product was at the leading edge of the mesoderm in vivo, curling of the mesoderm would place that former leading edge at a non-extremal portion of the explant (Figures 7-9 B, C). The system length is shortened in these computations to reflect two factors. First, the explant is only 70% of the original mesoderm. This makes the explant length 0.70 mm. Second, I have tried to account for the projection onto one dimension of the rolled tissue, which leads to a further shortening. A length of 0.55 mm is used for the explanted mesoderms. This length is used for both the stage 20 explants shown in Figure 7-9 and for the stage 14 in vitro studies described in section 73 . 7.7 Conclusion The inducer concentration distribution in the endoderm (Figure 7-4) has been generated by assuming a chemical action of the inducer and determining what levels of inducer are necessary to fire the reaction-diffusion peaks in the experimentally determined times (Smith and Armstrong, 1990) for each portion of endoderm tested. This inducer distribution could be a possible test of this model of induction. If the inducer could be quantitatively visualized in the future (after identification of the inducer, of course), then the model predicts that a relatively flat concentration profile would be found throughout the ventral portion of the Chapter 7. Modelling Induction 145 Distance 0.55 mm Product of Induction in Mesoderm In Vitro at 20 ! B 0.0 mm 0.55 mm explant 0.0 mm Distance 0.55 mm Figure 7-9: Model for the effect of explantation and tissue rolling upon the concentration profile of the induction product. A) the explant process, showing the size of tissue explanted and the tissue curling in vitro. B) The concentration profile of the product just prior to surgery at stage 20. C) The concentration profile of the product after explantation. Chapter 7. Modelling Induction 146 pharyngeal endoderm. Also, the hypothesis that in vivo induction is 40 times more efficient than in vitro induction is a possible point for experimental verification of the induction model. By considering general aspects of the nature of induction and the dynamics of mesoderm movement over the endoderm, a model has been developed which generates computed spatial distributions of the product of induction that agree quantitatively with the experimentally determined region of specification. As well, the concentrations of the induction product at each stage account for the temporal data for onset of heartbeat in the experiments conducted. This concludes the discussion of inductive dynamics and how they have been incorporated into the model. The next chapter will discuss the second stage of the heart-forming process, the post-inductive differentiation which leads to viable heart tissue. The dynamics discussed in this chapter and the resulting chemical distributions will serve as precursors to the post-inductive reaction-diffusion dynamics. The reaction-diffusion aspect of the model has already been alluded to in this chapter to explain how the inducer profile and inductive process were modelled. Chapter Eight will discuss the reaction-diffusion patterns and their correspondence with the data on normal heart development and the c mutant transplant and explant experiments. Chapter Eight Modelling Heart Formation 8.1 Introduction This chapter concludes Part Two of this thesis. In it, the computations of reaction-diffusion dynamics to model final differentiation of heart tissue in the axolotl will be presented. The results from one-dimensional computations will be presented for normal development and for the experiments described in Chapter Six. Two-dimensional computations will be discussed, and several examples will be presented. Finally, the results will be discussed in terms of predictions for future experiments and the possible chemical nature of the species in the model. The applicability of reaction-diffusion theory for the case of heart formation will be evaluated in comparison with other pattern formation theories. 8.2 Chemical Mechanism The modified Gierer-Meinhardt equations developed in section 5.2 are used for modelling the post-inductive dynamics of axolotl heart formation. These equations are: dA/dt = aS+ 4g(jf - uA + DA V2A (8. la) 61/dt = c'S'A2 - vl + D,S72I. (8.1b) These equations differ from the Gierer-Meinhardt equations in that they are derived from a chemical mechanism. In section 5.2, the enzyme mechanism leading to equations (8.1) was presented. It is a scheme of eight steps which involves the bimolecular activation of an enzyme E by the activator A, followed by the binding of a substrate R which enables the conversion of substrate S into a molecule of A. The inhibitor / is a competitive inhibitor of Chapter 8. Modelling Heart Formation 148 the enzyme. In addition,A bimolecularly catalyzes the production of/ from a substrate 5". A is made in a zeroth order process from substrate S and both A and / decay in a linear fashion. Substrates S and S" are assumed to be in excess concentration at all times and are treated as constant. The rates a, k, \i, c', and v are assumed constant, as is the equilibrium constant Ku The species E and/J, however, are functions of both time and space. The inductive dynamics discussed in Chapter Seven are incorporated into the concentration distribution of the enzyme E, which is postulated as the product of induction. The E concentration is shown at various stages in Figure 7-8. In Gierer and Meinhardt's original paper (1972), a gradient was placed in S. This put a gradient in two of the terms in equation (8.1a). The E gradient used in the present work allows for more direct control of the autocatalytic/inhibitory term. By incorporating the inductive dynamics as a parameter in equation (8.1a), the two stages of the dynamics, inductive and post-inductive, have been linked in a hierarchical manner. In light of the experimental evidence presented in section 6.4, the inducer, or a chemical which mimics the action of the inducer, could be a peptide growth factor such as PDGF or TGF-pi. In this case, the growth factor would lead to expression of the enzyme E, perhaps by a signal transduction pathway to switch on nuclear production of the enzyme. The enzyme could be either dissolved in the cytosol, or be a membrane receptor itself, such as an autophosphorylating tyrosine kinase receptor. E is non-diffusing on the scale of the model, that is intercellularly. We identify R with the diffusible rescuer described in section 6.3. R is a further activator of E, since the complex EA^R produces more A. We postulate that R is present at normal levels in wild-type mesoderm, but is absent in the c/c mutant mesoderm. This is the source of the mutant phenotype in the model. R is diffusible in the model, allowing for the rescue of c/c mesoderm by contact with wild type mesoderm and by exogenous application of endoderm Chapter 8. Modelling Heart Formation 149 extracts. The spatial distribution of R is calculated from a Fickian diffusion equation with no kinetic terms. Exogenous application of small RNAs (section 6.4) has been identified with rescue of stage 35 c/c mutant hearts. These RNAs could either contain the species corresponding to R in the model, or they could mimic the action of R. If R is strictly an extracellular molecule that diffuses through the extracellular matrix, then a receptor at the cell surface is required for R. This receptor could either be E, or the receptor could bind an intracellular E (the receptor could even affect E in some sort of cascade of chemical signals). Alternately, there could be some mechanism fori? to enter the cell. This could mean thati? is lipid-soluble (unlikely for an RNA), or that there is some passive channel for R through the membrane: a gap junction (probably too small for an RNA1) or more specific channel. The species A and / are the Turing morphogens, which achieve a patterned concentration distribution dictated by equations (8.1). A and / are likely to be more closely related chemically to their substrates S and S' than to either E or R. These substrates, and A and /, are likely to be intracellular. A or / is a necessary component of organized sarcomeres, either as a signaller for the initiation of sarcomere assembly, or as a structural component of the sarcomeres. In a reaction-diffusion mechanism with a true kinetic inhibitor, the spatial patterns of A and / are in phase with each other. That is, where there is the most activator there is also the most inhibitor. This contrasts with depletion models in which A and / are out of phase (Chapter Four). Therefore, either A or / could be the necessary patterned organizer of sarcomere assembly. We have fine-tuned our model to treat A as the entity responsible for final differentiation, but / could also be readily assigned this role. In the model, any portion of the system which has an A concentration above unity is considered to be differentiated to form heart. This is the threshold necessary for sarcomere assembly. The modelling shows the A and / diffusivities to be much slower than R. This fits with a picture of A and / as peptides which diffuse intercellularly either through gap junctions (for small A ori) or through some 1 Gap junctions allow passage of molecules less than 2000 daltons (Caveney, 1985). Chapter 8. Modelling Heart Formation 150 other passive channels. Chemically, A could possibly be tropomyosin, which is found at reduced levels in c/c mutant hearts (section 6.4). If this is the case, then A could not diffuse through gap junctions (see footnote 1), but would have to diffuse by some other type of passive membrane transport. The details of how the diffusivities are determined computationally and further implications for the chemical nature of the species in the model are discussed below. 8.3 Results Equations (8.1) are solved numerically (see Chapter Three) from stage 13 (55 hours of development) until the time of heartbeat onset (e.g. stage 35, or 120 hours, for normal development). The dynamics of E and R are incorporated into the reaction-diffusion equations (8.1) at each iteration. Random noise is also added to the A and / concentrations at each iteration to simulate thermal fluctuations. In the modelling, there is a consistent delay between the establishment of morphogen peaks ("firing the gradient" in the terminology of Gierer and Meinhardt, 1972) and the subsequent completion of myocardial differentiation to onset of beating. For all simulations, this delay is a constant 48 hours after the A peak has fired. The A peak is considered to have fired when it has reached a value of unity on the arbitrary concentration scale. This is about 5 times the initial homogeneous steady state value of A concentration, and less than 1/10 of its final peak value. The delay of 48 hours in the model is postulated as the time needed for the events of final heart differentiation and onset of heartbeat once the activator is present in above threshold concentration. Chapter 8. Modelling Heart Formation 151 8.3.1 One-Dimensional Computations Wild-Type Heart Development Once the concentration of E is fixed for each position and stage, as discussed in Chapter Seven, the other reaction parameters in the model can be fixed by modelling wild-type development. The parameters outside of the autocatalytic/inhibitory term in equation (8.1a) are basically the same as those used in the original Gierer-Meinhardt model (1972, fig. lc). The value of R is constrained by E in that their product (the coefficient for the autocatalytic term) must reconstruct the experimental data. E and R are the only dynamic coefficients in the model. Since the dynamics of £ are fixed by modelling the inductive dynamics, the time to beat can be reconstructed by proper choice of R. The maximum values of E given in Chapter Seven are for an R value of 0.36. The values of all other reaction parameters in the model are summarized in Table 8-1. Table 8-1: Parameter values for equations (8.1) Parameter a S k Kj li V c' S' Value 0.0006 1.0 1.0 1.0 0.0035 0.0045 0.0234 0.4 Peak width and the time required for peak firing are functions of both the reaction terms and the diffusion coefficients. Therefore, for given values of the reaction parameters, diffusivities were found to fit the differentiated region in the embryo and time to onset of heartbeat. The dorsal-ventral length of the differentiated region in the embryo was roughly Chapter 8. Modelling Heart Formation 152 calculated to be 0.50 mm, based on the observed diameter of 0.16 mm for the early heart tube. The diffusivities of A and / (for a given diffusivity ratio) are therefore constrained in order to produce a peak width corresponding to 0.50 mm, given the E gradient produced by induction. The necessary diffusivities are Dj^=2.0xl0"6 and Z)/=3.0xl0"5 in computational units. This converts to roughly D^=8.3xl0"10 cm2/s and Dj=1.25xl0'8 cm2/s using the fact that the 1 iteration that is 0.075 in computational units corresponds to 7.2 seconds (based on the relation of 12,000 iterations per day) and converting the distance variable of the computation to centimetres. The diffusivity ratio of 15 is small enough to allow convergence of peaks, as in the Humphrey stage 29 transplants below, and large enough to have significant growth rate, or maximum eigenvalue (Chapter Two). Only one parameter, Z3/j, remains to be determined. It is determined by modelling the stage 20 explants, which will be discussed below. DR is not a significant parameter in the wild-type modelling, as all wild-type tissue has an identical concentration of R. To account for wild-type heart development and the experiments with the cardiac mutant, only a very tight range of the reaction and diffusion parameters is possible. Computations of the concentrations of A, /, E and R, using the above parameters, are shown for wild-type heart development in Figure 8-1. Equations (8.1) were solved at 80 nodes, representing the 2.00 mm of length in the mesoderm (1.00 mm for each side). The gap shown in the centre for Figures 8-1 A and B represent the fact that the mesoderm sides are not in contact. The centre of the computation is a no-flux boundary to diffusion until stage 29. The dorsal boundaries of the computational mesoderm are no-flux boundaries for all time, representing the separation of the lateral mesoderm from the somites. The boundaries for the region of specification are always met by the limits of the enzyme gradient. The boundaries of the differentiated region defined by the A level do not change greatly over time, but the boundaries of the specified region, represented by the enzyme concentration, Chapter 8. Modelling Heart Formation 153 12 10 a 6 4 h 2 0 12 10 8 6 4 2 "I I I I t I I I I I I I I I I I 1 t I I I t I I I I t I I I I I I I I I M I I I I I I I I I I I I I I I I I | I I I I I I e a c f d B o "i 11111111111111111111 ut\m b e ' ' ' !i * ' # ' ' ' " 1111 v w i j.J.jf — m > - 1 . * a c f d 12 2 10 i f H 1 a o a S 6 _ -1**1 »/ 1 •f ' * # ' "i 111111111111111111 U ^ M J b e 1 •" 1 u ac • i • [ • l • imi i IA I 11 I I 1111111 i i 111 f • I I I d Distance n u n R Figure 8-1: Computations of normal heart development in a wild-type embryo showing the A, I (xO.l), E (xl5) and R (x5) concentrations. A) stage 20; B) stage 29, just before the leading edges of the mesoderm meet at the ventral midline; and C) stage 35, at the onset of heartbeat. Distances (from a to the left, and from c to the right) are dorsal-ventral distances from the leading edges of the mesoderm. Lettering a, b, c, d is the same as in Figure 7-7. Distances ab and cd remain 1.00 mm from stage 20 to stage 35. In A and B, a gap is left in the middle of the diagram to indicate that leading edges a and c have not yet touched. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 154 do. This agrees with a picture of constant heart size arising from a growing specified region. The constancy of the size of the differentiated region is largely due to the growth of the enzyme gradient, as discussed in section 4.3.2. The diffusivities of A and / are large enough that they will form a peak the width of the specified region if the enzyme gradient is static. This was tested with a static stage 29 enzyme gradient, computed from local steady state concentrations of the morphogens. Figure 8-1 shows that proper modelling of the enzyme gradient coupled with the right choice of rate and diffusion coefficients will account for the wild-type pattern. Next, it is shown that these parameters are also sufficient for the in vivo and in vitro experiments on the c mutant. All of these experiments involve juxtaposition of wild-type and c/c mesoderm. Experimentally, the contact of these tissues is sufficient to initiate a rescue of the mutant mesoderm. We model this rescue as the action of a diffusible substance, R, in equations (8.1). Diffusion of R is allowed in the computations when separated tissues come into contact, either through surgery or through the meeting of the mesoderm primordia at stage 29. The diffusivity of R is the last parameter to fit to experimental results. Stage 20 Transplants Figure 8-2 shows a stage 20 bilateral transplant (section 6.3) and the computation corresponding to it. Figure 8-2 A shows the experiment, in which both sides of mesoderm from a known wild-type donor are placed into a host from a wild-type/c mutant spawning (Smith and Armstrong, 1993). Of these, 25% of the hosts are c/c and the other 75% are wild-type in phenotype. The donor mesoderm is placed in the same position in the host. Approximately 70% of each side of the mesoderm is removed. The surgery heals quickly and the embryo develops a heartbeat, regardless of the host phenotype. Computations of this experiment, with a c/c host, are shown in Figures 8-2 B and C. For c/c hosts, the heartbeat is Chapter 8. Modelling Heart Formation 155 A 8 • I I I I I I I I I I I M I I I I I I I I I I I I I I I I I I I I I I I I I m i l l R Figure 8-2: The stage 20 bilateral transplant of Smith and Armstrong (1993). A) shows the bilateral transplant of wild-type precardiac mesoderm into a c/c host. B) shows the computed levels of A, /, E, and R (concentration adjustments as in Figure 8-1) at stage 20 when the surgery takes place. C) shows the concentrations of these species at stage 35, the onset of heartbeat. Distances are as for Figure 8-1. In B, the gaps on each side represent the no-flux boundaries present between the tissues just prior to implantation at stage 20. The central gap indicates that the mesoderm primordia have not yet touched. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 156 somewhat delayed from normal development (Table 8-2). Table 8-2: Time to Beat (days) for Stage 20 Transplants Transplant Wild-type Development4 Bilateral: wt into c/c Unilateral: c/c into wt Unilateral: wt into c/c Experimental2 2.15 3.0 3.2±0.4 3.2±0.4 Computational3 2.10 2.14 2.10 2.14 This establishes that c/c hosts cannot suppress induced wild-type mesoderms when the transplants are of this size. In the case of the bilateral transplants (and unilateral transplants with c/c hosts, see below) the average R level in the wild-type tissue drops to 70% of its original value upon contact with the c/c tissue. Sufficient levels remain to fire the morphogen peaks, A and I, but with a one hour delay from the wild-type time. In this case, the model qualitatively fits the experimental data (Table 8-2). In the unilateral transplants with wild-type hosts, the wild-type side fires morphogen peaks at the normal time. The model predicts no difference in firing between this case and normal development. In order to conduct transplant experiments with both c/c and wild-type hosts, unilateral transplants were conducted (section 6.3). Approximately 70% of one side of the mesoderm was removed from one animal and placed in the same position in a host animal. In both cases, that is c/c tissue transplanted into a wild-type host (Figure 8-3 A) and wild-type tissue transplanted into a c/c host (Figure 8-3 E), heartbeat was observed. This demonstrates that the presence of the induced wild-type mesoderm was enough to rescue the cardiac lethal tissue. Computation with the diffusible rescuer, R, re-creates the experimental phenomena. 2From S.C. Smith, personal communication. All times are for experiments run at 18-20°C; tissues were checked daily for beat. 3The computational time to beat is the time at peak firing + 48 hours. 4A11 experimental values are given as the time between stage 20 and the onset of heartbeat. 5From Bordzilovskaya et al (1979,1989), at 20°C. Chapter 8. Modelling Heart Formation 157 B 2 ' H Ml i in i m n iu t i I ill Dist&ace n u n R Figure 8-3: The stage 20 unilateral transplants of Smith and Armstrong (1993). A) shows the unilateral transplant of cic mesoderm into a wild-type host (host shown on right). B), C), and D) show the concentrations of A, /, E, and R (concentration adjustments as in Figure 8-1) at stages 20, 29, and 35, respectively. Distances are as in Figure 8-1. Gaps indicate no-flux boundaries. E) the unilateral transplant of wild-type mesoderm into a cic host (host shown on right). F), G), and H) show the levels of A, /, E, and R (concentrations as above) at stages 20, 29, and 35, respectively. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 158 This is most dramatic in the wild-type transplant into the c/c host (Figures 8-3 F-H), as this represents the smallest proportion of wild-type R levels to c/c R levels (which are zero). The final level of R after diffusion into the mutant tissue (35% of the original wild-type concentration) is still sufficient to fire the morphogens, A and /. The model can easily produce a peak of A and/for the other experiment (Figures 8-3 B-D). In this case, though, it is important that the diffusivities of A and / be fast enough to suppress peak formation along the dorsal edges of the specified region as R diffuses into the dorsal portion of the c/c donor tissue. In terms of the discussion in Chapters Two and Four, the wavelength must be long enough to suppress secondary peaks, or the diffusivity ratio, n, must be high enough to create areas of inhibition. Stage 20 Explants In vitro experiments analogous to the unilateral transplants were conducted (section 6.3; Smith and Armstrong, 1991). Stage 20 mesoderms were explanted from wild-type animals and from c/c embryos in the same sizes as the above experiments. These two mesoderms were then placed in contact with each other in culture (Figure 8-4 A). Due to mechanical forces in the epithelioid sheets, these pieces roll up. The effect of this rolling is modelled as discussed in section 7.6 (Figure 7-9). The wild-type tissue rescues the c/c tissue and two beating foci are observed. The times to beat are given in Table 8-3. Table 8-3: Time to Beat (in days) Tissue wild-type c/c Experimental6 4.0±0.5 6.3±0.4 for Stage 20 Explants Computational7 4.1 6.2 6From Smith and Armstrong (1991). All times are for experiments run at 18-20°C. Tissues were checked daily for beat. 7The computational time to beat is the time at peak firing + 48 hours. Chapter 8. Modelling Heart Formation 159 Figure 8-4: Computations of the stage 20 explants reported in Smith and Armstrong (1991). A) the explantation and juxtaposition of mesoderm from wild-type and c/c embryos. B) the concentrations of A, / (x0.25), E (x3) and R (x5) at stage 20, at the time of surgery (one iteration before contact between the tissues). The distance axis represents 1.1 mm, indicating that two explanted mesoderms are present. Length of explanted mesoderm is as in Figure 7-9. C) concentrations at 3.0 days, the time that heartbeat is experimentally observed for the original wild-type tissue, minus one day. D) concentrations at 5.3 days, when heartbeat is observed in the original c/c tissue, minus one day. The concentrations are shown prior to the experimental time of heartbeat in order to emphasize the different times of peak firing that correspond to the different times to heartbeat for the two tissues. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 160 We have achieved quantitative agreement with these data by modelling the rescue process as the action of the diffusible rescuer, R. Our modelling has indicated that R must diffuse quite fast. We calculate an approximate diffusivity of 8.3xl0"7 cm2/s for R. In the model we have postulated that A and / do not diffuse between the tissues in this case. The post-operational dynamics are shown in Figures 8-4 C and D. Stage 29 Transplants Humphrey (1972), conducted transplant experiments on stage 29 embryos (section 63). He took central pieces of mesoderm out of donor embryos and transplanted these into the same position in host embryos. Humphrey conducted both wild-type transplants into cardiac lethal hosts (Figure 8-5 A), and c/c transplants into wild-type hosts (Figure 8-5 G). The c/c transplants into wild-type hosts were rescued and led to beating hearts. Wild-type transplants into c/c hosts, however, did not rescue. This seems counter to the results of Smith and Armstrong with the earlier stage transplants and explants. Smith and Armstrong repeated Humphrey's experiments and found that the results were dependent on the size of transplanted mesoderm (Smith and Armstrong, 1993). Humphrey used approximately 25% of the mesoderm in his transplants. Smith and Armstrong found that if 70% of the mesoderm was used for the transplants (consistent with the size of tissue for the stage 20 transplants), then rescue did occur for wild-type transplants into c/c hosts (Figure 8-5 D). The difference between the Humphrey and the Smith and Armstrong experiments was computed purely as a difference in transplant size. This reproduces the experimental results. Figures 8-5 B and C show a computation of the Humphrey experiment. The A and / peaks have been suppressed by stage 35, indicating that no beat will occur. Computation of the Smith and Armstrong experiment with a wider transplant is shown in Figures 8-5 E and F. Strong A and / peaks are maintained through stage 35, indicating that a heartbeat will be established. In the context of the model, it is easy to see how the transplant size could have this effect. In the Humphrey experiment, the rescuer (R) concentration in the wild-type tissue is at normal Chapter 8. Modelling Heart Formation 161 Figure 8-5: Computations of stage 29 transplants as reported in Humphrey (1972) and Smith and Armstrong (1993). A) The Humphrey (1972) transplant of ventral mesoderm from a wild-type embryo into a c/c host. B) and C) concentrations of A, /, E and R (concentration adjustments as in Figure 8-1) at stages 29 and 35, respectively. Distances are as in Figure 8-1. D) the transplant of 70% of the mesoderm from a wild-type embryo into a c/c host (Smith and Armstrong, 1993). E) and F) levels of A, /, E and R (concentrations as above) at stages 29 and 35, respectively. G) the Humphrey (1972) transplant of ventral mesoderm from a c/c embryo into a wild-type host. H) and I) computations of the A, I, E and R (concentrations as above) levels at stages 29 and 35, respectively. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 162 levels, but the transplant is only 25% of the total mesoderm area. Upon surgery, the substrate rapidly diffuses to the average concentration between the wild-type and c/c tissues. This average value is weighted by the 75% c/c tissue, and is hence too low to sustain the A and / peaks. In the Smith and Armstrong experiment, the transplant is 70% of the total mesoderm area. The average level to which the substrate diffuses after surgery is hence much higher. It is sufficient to sustain the A and / peaks. Finally, Figure 8-5 G shows the reverse Humphrey experiment in which a c/c transplant was made into a wild-type host. It is observed in this case that a heartbeat does occur. It is shown in the computations that the mutant tissue in this case is small enough to be rescued by the wild-type host (Figures 8-5 H and I). The rescuer levels are at wild-type levels throughout 75% of the mesoderm at the time of operation. The rescuer concentration, therefore, only drops slightly from the wild-type level as it diffuses into the transplanted tissue. Both sides of the surgical cut fire morphogen peaks. This is due to residual A present in the wild-type host. Although two peaks are fired, the resulting A levels are above the threshold of unity necessary for sarcomere assembly throughout the heart forming region. The model predicts the single beating heart observed in this rescue experiment. It is important in this computation that the diffusivity of A is fast enough to allow A to flow into the ventral part of the embryo. In other terms, the diffusivity ratio, n, must be small enough so that the ventral portion of the embryo is not inhibited by the range of/. Early in vivo studies In this section, I will discuss computations of the early in vivo studies, such as cardia bifida, extirpation of mesoderm sides, and Copenhaver's (1926) work delineating the heart-field mesoderm (section 6.2). Cardia bifida is created experimentally by blocking the union of the mesoderm sides with a mechanical barrier. In the model, retention of the no-flux boundary at the ventral midline until onset of heartbeat re-creates this barrier. This results in the maintenance of two separate peaks past the onset of heartbeat (Figure 8-6). This corresponds Chapter 8. Modelling Heart Formation 163 c 0 +* o h . *-c 9 O c 0 O 12 10 8 6 4 2 0 Distance X0.1 E x 1 5 R X5 Figure 8-6: Cardia bifida computation at stage 35. Experimentally, a mechanical barrier prevents meeting of the leading edges of the mesoderm. Computationally, the no-flux boundary is maintained between the two halves of mesoderm past stage 29. Two separate beating foci are produced in the computation, which reflects the experimental result. The no-flux boundary is represented by the 1 cm gap between the two halves of the computation. Chapter 8. Modelling Heart Formation 164 to two separate areas of mesoderm differentiating into heart tissue. The fact that the reaction-diffusion peaks fire well before onset of heartbeat makes the model robust to perturbations at late stages (stages later than peak firing). Likewise, extirpation of parts or whole sides of mesoderm still leads to heartbeat. The mesoderm sides are basically independent in the model, reflecting the experimental evidence. In modelling normal development, removal of the no-flux boundary at stage 29 has the effect of allowing diffusion between the two sides. This is a relatively minor effect, as both sides have fully fired morphogen peaks prior to stage 29. Copenhaver's (1926) experiments delineating the heart field in Ambystoma punctatwn can be modelled by removing the enzyme (E) gradient. This gradient corresponds to the region of specification found by in vitro methods. If the portion of the computational system corresponding to the induced region is removed, then there is no enzyme left to catalyze the self-enhancement of the activator, A, and no heart tissue will differentiate (Figure 8-7 A). This corresponds to the experiment where non-heart forming mesoderm was used to fill the wound and there was no further mesoderm movement. Heartbeat was found if mesoderm movement was allowed. The corresponding computation is shown in Figure 8-7 B. Here, the E gradient is removed at stage 27, representing surgical removal of that mesoderm, and mesoderm movement is allowed until stage 29. In the axolotl, the mesoderm advances 0.25 mm between stage 27 and stage 29. By inspection of Figure 6-4, the mesoderm travels about one-half of the length of the specified region between stages 27 and 29. Assuming that the specified region in A. punctatum is of similar size to that in the axolotl, this is a distance of approximately 0.23 mm. This is roughly in agreement with the distance travelled in the axolotl during the same stages. These stages are much later in A. punctatum, and they are much closer to each other in time. Stage 27 is at approximately 120 hours, and stage 29 is at 125 hours (Harrison, 1969). Also, surgery was performed at various late tailbud stages, not always at stage 27. A complete analogy to the axolotl should not be made, therefore. The Chapter 8. Modelling Heart Formation 165 i xO.I E • ' • I I I R x15 x6 1 8 O 1 -Distance B Distance Figure 8-7: A) Copenhaver mesoderm extirpation with non-heart forming tissue placed into wound, shown at stage 35. All concentrations are set to zero in the middle half of the computational system at stage 27. Diffusion is allowed into the non-heart forming tissue in this computation. Since no mesoderm movement is allowed, no new mesoderm is exposed to inductive endoderm, and no new enzyme is produced in the mesoderm. There is no morphogen peak firing in this computation, which agrees with the lack of heartbeat observed experimentally. B) the ventral space is left unfilled after extirpation of the mesoderm. Therefore, mesoderm movement over inductive endoderm is allowed and new enzyme is produced in the mesoderm. Diffusion is allowed into the ventral space in this computation. A single heart is predicted by the model if the leading edges of the mesoderm meet, two hearts if the leading edges do not meet. In the embryo shown in Figure 6-4, the mesoderm primordia do not meet. The computation shows this case, in which the mesoderm velocity is not enough to bring the primordia together after surgery. Copenhaver (1926) observed double heart formation in the majority of cases. Chapter 8. Modelling Heart Formation 166 computations in Figure 8-7 are meant to model the experiments in a qualitative sense, not quantitatively. The model predicts the formation of one heart if the mesoderm sides fuse at the ventral midline, and two hearts if they do not meet. Since the region outside of the specified region is shorter than the specified region in Figure 6-4, the mesoderm primordia probably do not have a chance to meet prior to heartbeat. Copenhaver (1926) reported that in eleven out of sixteen embryos with beating hearts, double hearts did indeed form. Exogenous Stage 35 Rescues Late stage rescues of c/c embryos have been reported by several groups. First, Lemanski et al. (1979) found that stage 35 c/c hearts could be rescued by contact with stage 30 wild-type endoderm. The same result was obtained by culturing stage 35 c/c hearts with an RNA extract from wild-type endoderm (Davis and Lemanski, 1987). Smith and Armstrong (1990) corroborated these results for rescue of stage 35 c/c hearts. They tested RNA preparations from both stage 14 and stage 35 endoderm. Since these preparations had no inductive effect on stage 14 mesoderm, the RNA rescuer is identified with the rescuer, R, in the model, not the product of induction, E. These late stage rescues can be modelled by application of R to a c/c system at the time corresponding to stage 35. As the specific geometries of the tissue cultures differ and are complex, R levels were instantaneously raised from c/c to wild-type levels at stage 35 in the computations. This allows for general modelling of the different experiments (endoderm, endoderm-conditioned media, and RNA extract) in one computation. Because the diffusivity of R is relatively fast, this instantaneous rise of R in the mesoderm is not a large approximation. Computation of these late stage rescues are shown in Figures 8-8 A and B. Chapter 8. Modelling Heart Formation 167 xO.1 E x15 111111 R X5 ! Distance B ] Distance Figure 8-8: Stage 35 rescue of c/c heart. R concentration is brought to wild-type levels at stage 35. This simulates experiments wherein contact with wild-type endoderm, culture with (wild-type) endoderm conditioned media, or culture with an RNA extract of (wild-type) endoderm at stage 35 resulted in heartbeat for c/c mesoderm. Since the specific geometries of the experiments are unknown, this computation is intended to model the general effect. Actual diffusion of the rescuer, R, into the mesoderm from an external source would take longer than modelled. However, R has a very fast diffusivity, so these results are generally accurate. A) Computation at one hour prior to stage 35. B) Computation at one day after stage 35 rescue. The width of the morphogen peaks is not intended to correspond quantitatively to the width of the beating tissue, as the specific tissue geometry for these in vitro experiments is unknown. Detailed modelling of the spatial scale was not conducted. Chapter 8. Modelling Heart Formation 168 8.3.2 Two-Dimensional Computations The one-dimensional computations were conducted first because they are fast and provide an easy graphical interface. These computations were used to develop the model and find the parameters which are consistent with experimental observations. However, the mesoderm is a two-dimensional sheet of cells and this fact could have important implications for the dynamics. Therefore, all of the one-dimensional computations described above were repeated in two dimensions. This enabled study of the effect of anteroposterior extension in the mesoderm on the reaction-diffusion dynamics. As the mesoderm is a thin sheet of tissue (several cells thick), it was felt that the essential dynamics should be captured in a two dimensional computation. The very small third dimension should not affect the results. Therefore, all computations were done in two dimensions, but not in three. Due to the fast diffusivity of R, the Peaceman-Rachford method (Chapter Three) was used to give reasonable computation times. This method is most easily implemented on a square grid, with equal spatial step in the two directions. There are several orientations which can be fit onto the shape of the mesoderm (Figure 8-9). If the square is made to fit over the most Figure 8-9: Side views of axolotl embryos showing possible orientations of a square grid with respect to the two-dimensional shape of the mesoderm. The specified region is shaded. In A) (stage 20) and B) (stage 24) the square grid covers the most mesoderm area possible. This orientation represents the leading edge of the mesoderm with a corner of the square. Computation of diffusion across the ventral midline after stage 29 would be difficult and inaccurate with this grid orientation. In C) (stage 20) and D) (stage 24) an edge of the square grid represents the leading edge of the mesoderm. This orientation allows more realistic modelling of the diffusional dynamics after stage 29. However, not as much mesoderm area is covered by the square. Line segment Id ' shows the orientation and extent of the one-dimensional computation in each case. Chapter 8. Modelling Heart Formation 169 B Figure 8-9 Chapter 8. Modelling Heart Formation 170 mesoderm possible (Figures 8-9 A and B), then the leading edge of the mesoderm corresponds to a corner of the square, with an adjacent edge corresponding to the anterior edge of the mesoderm. This leads to difficulties when the leading edges of the mesoderm come into contact at stage 29. In the computation, the diffusion across the ventral midline would have to occur through the corner points of the squares only. This would be biologically incorrect, since the leading edges of the mesoderm do have extent and do allow diffusion across a region of tissue after stage 29. Also, diffusion between single points would make the timing of the rescues quite different in two dimensions from the analogous cases in one dimension. The other option for the orientation of the square is shown in Figures 8-9 C and D. This allows for diffusion along an edge between the two sides of mesoderm . It does not cover the mesoderm quite as well as the option in Figures 8-9 A and B, but it is preferable due to the better account of the region of heart formation, and as a result, of the diffusive process. Also, the discrepancy between mesoderm coverage for the two orientations decreases from stage 20 (Figures 8-9 A, C) to stage 24 (Figures 8-9 B, D). Using a square grid to represent the mesoderm, with diffusion allowed along an edge after stage 29, leads to a simple extension of the one-dimensional computations. Again, the mesoderm does not change size or shape greatly during the events of heart formation, so only a static domain is considered. The rows of the two-dimensional computations (in the x-direction, dorsal-ventral in the embryo) have the same input parameters (R value, E gradient, etc.) as the corresponding one-dimensional computation (Figure 8-10). In other words, each row is a copy of the input parameters of the one-dimensional case. Each row is represented by 80 nodes, as in the one-dimensional computations. The second dimension is represented by columns in the y-direction (anteroposterior in the embryo). These columns contain identical values of the input parameters (Figure 8-10). Each column is represented by 40 nodes. This follows from an assumption of a square mesoderm side. Each computation corresponding to the 1.00 mm by 1.00 mm square mesoderm is conducted on a grid of 40 by Chapter 8. Modelling Heart Formation 171 : : . ' : • : • • ' • : ' • '-:'•. • : • : • • : ' : ' . : Figure 8-10: Computations of the stage 20 unilateral transplant of c/c mesoderm into wild-type host in two-dimensions (compare to Figure 8-3 B). Greyscale represents concentration: white is the highest concentration, black is the lowest. Top: the concentration pattern of enzyme, E, at stage 20. Any row is a copy of the one-dimensional E gradient. Any column contains 40 copies of a given value of E. Bottom: the concentration pattern of the rescuer, R, at stage 20. Again, any row is a copy of the one-dimensional case and the R values in any column are identical to one another. Chapter 8. Modelling Heart Formation 172 40 nodes. The total computation, representing both sides of the mesoderm, is conducted on a grid of 80 by 40 nodes. The columns link the rows by diffusion of A, / and R. There are, a priori, three possible outcomes: a) Each row is a repeat of the one-dimensional computation, with no patterning in the y-direction. b) There is additional patterning in the y-direction (repeated A and / peaks) but the minima between the peaks are above the threshold necessary for sarcomere assembly. Possibilities a and b imply a single beating heart, c) There is patterning in the y-direction with low enough minima to produce multiple hearts in an anteroposterior column. This would be a failure of the model. All computations had either outcome a or b. The patterning in the y-direction took longer to develop than the patterning in the x-direction. Therefore, it was seen chiefly in the computations of in vitro explants, which covered longer time periods. Overall, no quantitative difference is found for the agreement between experiment and theory for the two-dimensional computations as compared to the one-dimensional computations. The time scales for peak firing in the two-dimensional model still agree within experimental error with the observed times for onset of heartbeat. The timing of the rescues in the two-dimensional case is very close to that of the one-dimensional computations because R is diffusing along a linear front. Since the value of R remains the same in any column, there is no net diffusion in the y-direction. The R concentration in any row of the two-dimensional computation remains very much the same as the one-dimensional computation for all times (Figure 8-11). Some of the resulting reaction-diffusion patterns are shown in Figures 8-12 and 8-13. Figure 8-12 shows a two-dimensional computation of wild-type heart development, corresponding to the one-dimensional computation of Figure 8-1. The result at stage 35 closely (but not entirely) resembles 40 repeats of Figure 8-1 C stacked side-by-side in the y-direction with no Chapter 8. Modelling Heart Formation 173 Figure 8-11: The R concentration is shown for the unilateral transplant of c/c mesoderm into wild-type host (compare Figure 8-3 C) just after stage 29, as the diffusive rescue is occurring. Greyscale is as in the previous figure: white is high concentration, black is low concentration. Chapter 8. Modelling Heart Formation 174 Figure 8-12: Two-dimensional computation of wild-type heart development (compare Figure 8-1). The spatial dimensions are plotted in the x- (dorsal-ventral) and y-(anteroposterior) directions, while the A (white) and E (green, xl5) levels are plotted in the z-direction. The computation is shown at stage 20 (top), stage 29 (middle), and stage 35 (bottom). The moire' patterns at left and right are an artefact of the computer graphics, but serve as shading for regions where A=E=0 (undifferentiated and unspecified regions). This computation shows the difference between the specified (positive E) and differentiated (A>1) regions. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 175 interaction between the repeats. Figure 8-13 shows the final stage of a two-dimensional computation of the stage 20 explant juxtaposing wild-type mesoderm with c/c mesoderm. Compare this to the one-dimensional case in Figure 8-4 D. In the two-dimensional computation, the morphogen patterning in the y-direction can be seen. This is an example of how a precursor gradient can affect the morphogen pattern. In the x-direction, the gradient in E constrains the morphogen pattern into the single peak corresponding to the single heart. In the y-direction, there is no precursor gradient and the morphogen pattern has multiple peaks. The pattern is freer to conform to its chemical wavelength in the absence of the gradient. The pattern of 11/2 peaks in 1 mm has a slightly longer wavelength (approximately 0.67 mm) than that predicted from linear theory (approximately 0.56 mm for the value of E along the maxima of the gradients). Perhaps some gradient effects are still present. However, the relatively close agreement between the two wavelengths does show the predictive power of Turing theory for nonlinear models. The A concentration is still above unity throughout the peak area, so the model predictions are still the same as in the one-dimensional case. The patterning in the y-direction does not predict multiple areas of heart differentiation. Since this secondary pattern depends on homogeneous concentrations of precursors, it would be destroyed by application of an anteroposterior gradient in one of these species. Anteroposterior gradients do exist in vertebrate heart formation, as reported recently for atrial/ventricular differentiation in zebrafish (Stainier and Fishman 1992,1993) and in chicks Figure 8-13: Two-dimensional computation of the stage 20 in vitro rescue of c/c mesoderm by wild-type mesoderm (compare Figure 8-4 D). The anteroposterior axis is plotted in the y-direction, while the dorsal-ventral direction is plotted in the x-direction. In the greyscale representation, black is the lowest concentration and white is the highest concentration. The top, middle and bottom show the concentrations of A, / and E, respectively, at 6.3 days past stage 20, the time at which both tissues are observed to beat experimentally. This computation shows the secondary pattern that arises in the y-direction in the absence of any precursor gradient in that direction. From Holloway et al. (1994). Chapter 8. Modelling Heart Formation 176 sUsss Figure 8-13 Chapter 8. Modelling Heart Formation 177 (Yutzey et al., 1994). Therefore, this secondary patterning might not exist in the biological system. 8.4 Discussion Armstrong (1989) claimed the experimental work on axolotl heart formation as the first concrete example of reaction-diffusion dynamics in vertebrate development. This claim was based on the fact that there were obviously important post-induction dynamics present and that these dynamics required a level of control and self-organization that could be provided by a reaction-diffusion scheme. This has now been supported by a quantitatively detailed application of a reaction-diffusion model. By including the inductive dynamics in the E parameter of a Gierer-Meinhardt type model, a two stage hierarchically linked system was developed. Figure 8-12 summarizes the chemical patterns generated in the two stages: induction gives the pattern of E which defines the specified region; while the post-inductive dynamics give the pattern of A, which defines the smaller region of myocardial differentiation. By combining models for the inductive and post-inductive events, we have been able to reproduce the experimental observations to date, and to make several predictions (sections 8.4.1 and 8.4.2). The close collaboration between experimentalists and theorists in this work has allowed a more rigorous test of the reaction-diffusion hypothesis than is frequently possible. The requirements of modelling led to more quantitatively detailed experimental work on some aspects of the phenomena than if the modelling had not been in progress. These include: a) the position of the leading edge of the mesoderm relative to the endoderm at a series of stages and b) the boundaries of the specified region of the mesoderm. Chapter 8. Modelling Heart Formation 178 Quantitative agreement was reached between computational values and experimental values for both the boundaries of the specified and differentiated regions, as well as the times to beat, within experimental error. As detailed in the results for the stage 29 transplant experiments, the experimental results for the rescue of c/c mesoderm by wild-type mesoderm was reproduced from transplant size considerations alone. This fact is further indication that a diffusible rescuer is present. 8.4.1 Predictions for Experiment In addition to the possibilities for experimental testing of the model of induction discussed in section 7.7, predictions can be made from the whole model, with regards to both induction and heart differentiation. These are as follows. 1) There is a possibility of producing multiple hearts. In the absence of gradient effects, and with a wavenumber greater than one, the reaction-diffusion dynamics will generate multiple concentration peaks. With low to moderate diffusivity ratios (as in the present case of n=15) the spacing between peaks will be regular, and generally will not greatly exceed the peak width. In other words, the reaction-diffusion system will obey the linear Turing dynamics. If these conditions are met, it might be possible to generate multiple hearts in vivo, e.g. laterally located repeats of the normal mid-ventral heart. To do this experimentally, the induced gradient would have to be either completely flattened in vivo, or the mesoderm pulled over the inductive endoderm faster than in normal development. It will be observed that the E gradient in the computations is already quite flat, and yet it constrains the morphogen distribution to a single peak. This constraint is due not only to the concentration gradient in E, but also to the fact that the E gradient grows throughout development, as discussed in Chapter Four. At early stages, positive E concentration is constrained to a small Chapter 8. Modelling Heart Formation 179 region. Consequently, the morphogens are able to form only in that small region. At later stages, when the E concentration is positive in a region larger than the chemical wavelength, the already established, central morphogen peak inhibits lateral morphogen peak formation within the range of inhibition defined by the diffusivity of/. Multiple hearts have been observed in nature, as in a chick with 7 relatively normal hearts observed by Verocay (1905). In this case, the model would predict that the inductive gradient had been flattened, or that the mesoderm had progressed over the endoderm at a much earlier stage than in normal development. By progressing over the endoderm at an early stage, a portion of mesoderm larger than the chemical wavelength could be induced (relatively) simultaneously. This would negate the peak formation constraints that are postulated for normal heart development, namely precursor gradients and growth of those gradients. 2) The Gierer-Meinhardt type of reaction-diffusion model was adopted for this work because of its sensitivity to suppression of multiple peaks by application of an input gradient. In this type of model, the peaks of activator and inhibitor coincide. It is therefore expected that both A and / are to be found in greatest concentration in the heart-forming region in vivo. However, the morphogens are the most labile of the species in the model, and discovery of the precursors, E and R, is predicted to be easier. 3) For the unilateral in vivo rescues of c/c mesoderm by wild-type mesoderm, the computed area of differentiation is always skewed towards the original wild-type tissue (Figure 8-3). This is due to the fact that the activator in the wild-type tissue takes some time to diffuse into the mutant tissue. While experimental verification of this would be difficult, we predict that more wild-type tissue is incorporated into the hearts of these unilateral transplants. Chapter 8. Modelling Heart Formation 180 8.4.2 Molecular Possibilities Only a limited set of model parameters will reproduce all of the experimental results. On this basis some more speculative predictions can be made as to the possible chemical nature of the substances in the model. First, it is necessary that the diffusivities of the three diffusible species (A, I and R) be tightly controlled. To account for the proper heart size, it is necessary for the diffusivities of A and / to be relatively slow (Dj^=S3xlO'10 cm2/s and Z)y=1.25xl0'8 cm2/s). These diffusivities are too slow for simple diffusion through gap junctions (see Safranyos and Caveney, 1984). Either A and / are large molecules that diffuse by passive transport through carrier proteins, or they are smaller molecules that are sterically or electrostatically hindered in passing through gap junctions. Since the diffusivity of A is 15 times smaller than the diffusivity of /, they could be quite different molecules diffusing via different pathways. R must diffuse much faster (D/p8.3xl0~7 cm2/s) in order to account for the experimental rescue time in the stage 20 explants. Such a fast diffusivity could lead to two pictures of the chemical nature ofi?. First, R could interact with the other substances in the model intracellularly. This would require R to be a relatively small molecule. It would probably be lipid-soluble, such as a steroid hormone or a retinoid. Since we have modelled diffusion in the stage 20 explants as blocked for A and / (blocked on this time scale; a decrease in D^ and Dj would be sufficient), but permitted for R, the transport mechanism must be different for R than for A and /. For instance, a lipid-soluble molecule could still diffuse between two in vitro tissues that had not formed functional gap junctions. A lipid-soluble R would naturally remain in tissues which were suspended in aqeous solution. As discussed at the beginning of this chapter, lipid-solubility is unlikely for the RNA rescuer postulated by Lemanski's group (section 6.4). It is possible, however, that this externally applied RNA could affect the production of R in c/c tissue through some signalling pathway. Other scenarios can be Chapter 8. Modelling Heart Formation 181 postulated in which an intracellular R diffuses via gap junctions and A and / diffuse via carrier proteins, or vice versa. A second possibility is fori? to be a larger extracellular molecule. In this case, R would need a membrane receptor in order to affect events in the cell. This receptor could be E, which could be an autophosphorylating tyrosine kinase, as discussed at the beginning of the chapter. E would then bind A and I upon extracellular binding of R. If E itself was not the receptor, but was completely intracellular, it could still be signalled by the R receptor in some sort of reaction cascade. Signal transduction pathways are promising systems for exhibiting the dynamics of the model, provided that the pathways involve reversible binding of R, A and /, in agreement with the proposed mechanism (equations 5.4). Many of the known facts of signal transduction are those of irreversible binding, because the stable products are capable of isolation and biochemical characterization. If R is extracellular, some additional mechanism, such as a reversible binding to components of the extracellular matrix would be necessary to prevent a great deal of loss to the aqueous culture medium. The model will accomodate some loss, however. Investigation of the possible mechanisms for transport of A, I andR could be initiated by blocking different types of membrane transport (e.g. gap junctions) in the wild-type-c/c rescue experiments and observing whether the rescues still occurred. The possible relationship between the model and signal transduction suggests a different possible source of the bimolecularity essential to the nonlinear dynamics. Rather than two molecules of A binding to one of E, bimolecularity could arise from the interaction of two EA complexes by intermolecular autophosphorylation. This is a known feature of many membrane-spanning receptors with intramolecular tyrosine kinase domains, such as the epidermal growth factor (EGF) or platelet-derived growth factor (PDGF) receptor (Panayotou and Waterfield 1993, review). Such a source of nonlinear dynamics was Chapter 8. Modelling Heart Formation 182 suggested by Lisman (1985) and used by Harrison et al. (1988) in connection with morphogenetic events at the cell surface of the alga Acetabularia mediterranea. If this dimerization of receptors is the source of bimolecularity in the reaction-diffusion mechanism, a factor of E2 would appear in equation 8.1a, instead of simply E. This means that the levels of E developed in this work to account for the experimental data would actually be the values of E2. Twice the E molecules would have to be present to get the same effect, so an inducer molecule would have to produce two E molecules. 8.4.3 Reaction-Diffusion vs. Mechanochemistry Kinetic theory of pattern formation is not restricted to reaction-diffusion. It is only one of several classes (Harrison 1993, Ch. 8). Another important class is mechanochemical theory, in which mechanical forces are explicitly accounted for in the modelling. If theory is to be incorporated as a major area of interest in developmental biology, then future experiments need to be designed to distinguish between reaction-diffusion and mechanochemistry. Only through a body of experimental data will the actions of the different types of dynamics in development, whether chemical, mechanical, or some combination, be elucidated. In the present study, it is clear that there are stresses in the mesoderm. But in explant experiments these forces are drastically changed, without the post-operative restoration to near normal forces that presumably occurs in transplant experiments. The fact that morphogenetic interactions are only mildly changed, and that our model will readily reproduce the results of both types of experiment, appears to be strong evidence in favour of reaction-diffusion and against mechanochemistry in this instance. In addition to the differing mechanical conditions in vivo and in vitro, a putative mechanical defect in c/c tissue would have to be capable of swamping the activity of wild-type tissue, as in the Humphrey transplants (Figures 8-5 A-C), as well as accounting for the tissue size Chapter 8. Modelling Heart Formation 183 differences seen in the rest of the transplant experiments (especially between Figures 8-5 A and D). Such activity by a mechanical mechanism seems highly unlikely. The experimental evidence for c rescue fits much more easily within a diffusion model. 8.4.4 Reaction-Diffusion vs. Inhibitory Field Theory Sater and Jacobson (1990) discussed the concept of inhibitory fields (Wigglesworth, 1940) with respect to the problem of loss of specification in Xenopus heart development. Wigglesworth's idea is that cells newly committed to differentiate and initiate a new organ will produce a surrounding inhibited region in which cells are prevented from becoming committed. This requires outward diffusion of an inhibitor, or inward diffusion of a substance needed for differentiation. The concept is within the reaction-diffusion category in the most general sense (Chapter Four). Inhibitory field models remain useful in relation to partly-ordered patterns (Lacalli and Harrison, 1978b). But for either well-ordered patterns or precise placement of single organs, a reaction-diffusion mechanism must have additional features. These, as first shown by Turing (1952), are that the activator must also diffuse (more slowly than the inhibitor) and that the inhibitor must interact with the activator. In terms of the continuum that reaction-diffusion theory and inhibitory field lie upon, based on relative diffusivities (Chapter Four), the present model is well within the realm of Turing dynamics, as it requires a diffusivity ratio of 15. The dynamics of an inhibitory field do not begin to dominate until higher diffusivity ratios. The degree of control of final pattern (reproducibility of heart size and position) for the experiments modelled requires the robustness of Turing dynamics. Chapter 8. Modelling Heart Formation 184 8.5 Conclusion In this chapter and the previous two, the localized reaction-diffusion theory developed in Part One has been applied to a problem of vertebrate organogenesis. The biological evidence for heart induction by the endoderm and post-inductive dynamics was presented. Within the constraints of this experimental evidence, a model was developed that incorporates both induction and final heart differentiation. Only a small range of the model parameters accounts for the experimental data. For this reason, predictions can be made as to the chemical nature of the species in the model. Also, predictions for future experiments can be made. It is hoped that the collaborative nature of this work is indicative of future efforts in developmental biology. It is only through the wealth of detailed information that the model has been constrained to such a narrow range of parameters that predictions can be made. Much of this information was unpublished data from Armstrong's lab in Ottawa. Attempting a model on published data alone would be a far less rewarding experience, as the type and amount of data necessary to build up a physicochemical model is not commonly published in biological papers. Detailed modelling of the phenomena of heart formation has shown that reaction-diffusion theory is more promising than either the inhibitory field theory or mechanochemistry, in this case. In applying kinetic theories to biological processes, it is important to distinguish between the possible kinds of underlying dynamics. This analysis of axolotl heart development has provided a physicochemical framework through which to analyze the experimental data. This framework accounts for more of the heart forming phenomena than models which depend significantly on mechanical forces, or which do not have reproducible steady state solutions, such as the inhibitory field theory. Chapter Nine Conclusion This thesis has been presented in two Parts: theory and application. However, the overall aim has been to address a single question: how can one, through chemical dynamics, account for localization of structures in biological systems? Given this basic assumption, that biological development is determined by chemical processes, what are the processes which give the observed progression of tissue differentiation in a spatially controlled manner? This question has been addressed within the general framework of Turing's theory of reaction and diffusion of two chemical intermediates Original results on this topic have been presented in Chapters Four, Five, Seven and Eight. In Chapter Four, I explored mechanisms whereby reaction-diffusion pattern can become localized. Chief among the findings are: a) high diffusivity ratio gives localized, partly-ordered pattern in the Brusselator and Gierer-Meinhardt models, forming a link between the theories of Turing and Wigglesworth; b) the Gierer-Meinhardt model is more susceptible to localization through precursor gradients than the Brusselator; and c) gradient growth gives greater localization than can be achieved with static gradients. Overall, the Gierer-Meinhardt was found to have a greater tendency to form localized pattern than the Brusselator, consistent with earlier work on these models by others. Since the original Gierer-Meinhardt model was not directly based on a chemical mechanism, one was developed in Chapter Five. Chapter 9. Conclusion 186 The theory was then tested against the evidence from a case of localized pattern formation in biology: heart development in the axolotl. Modelling involved consideration of the events of induction by the endoderm, as well as the process of myocardial cell differentiation in the mesoderm. The model for induction was presented in Chapter Seven. Data on mesoderm velocities, inductive capacity and the region of specification were incorporated, leading to formulation of a gradient in E, the putative enzyme responsible for morphogen production in the proposed chemical mechanism. This E gradient served as prepatterned input to the reaction-diffusion model of heart differentiation in Chapter Eight. Agreement with the experimental evidence was achieved, in terms of the times to heartbeat as well as the regions of tissue differentiation, for various combinations of wild-type and cardiac lethal (c) mutant tissue in vivo and in vitro. Based on the limitations placed on model parameters to reach such agreement, predictions for future experiments were made. These were presented in section 8.4.1. They can serve not only to guide future experiments, one of the functions of theory, but more specifically to guide experiments which will test the theory. The results of Chapter Eight show the versatility of the reaction-diffusion approach for this application. The disadvantages of other approachs were addressed in sections 8.4.3 and 8.4.4. Because of this, it is likely that future experiments will serve to refine the model, perhaps greatly, but that its essential reaction-diffusion nature will remain. Further modelling is possible at the moment, for example in the use of more realistic domain Chapter 9. Conclusion 187 shapes in the two-dimensional computations. However, I think that further modelling would be premature until there is experimental response to the present work. Such testing of the model predictions should lead to refinement of the model. This could be in many aspects, such as the details of the kinetics in the proposed mechanism, or the assumptions upon which the modelling of induction is based. These assumptions include: the geometry of the in vitro mesoderms and the effect of this on the E gradient; the non-diffusion of activator (A) and inhibitor (I) between tissues in vitro; the lack of movement in the dorsal 25% of the mesoderm; and the saturation of E. If the kinetics of heart forming chemicals can be determined in the future, then the chemical mechanism could be altered accordingly, and the interpretation of the results in Chapter Eight would become more precise. Determination of the character of even one of the constituents of the model could affect the interpretation in Chapter Eight, i.e. the discovery of whether or not the rescuer (R) is an extracellular or intracellular messenger. Future experiments may also help to unravel some of the aspects of noise or variability in the current data. Only thermal noise has been simulated in this work. Further stochastic aspects of the experimental data were discussed in Chapters Three and Seven. The Heart Determination Coefficient used to grade viability in axolotls has two major components: time to beat and percent of beating hearts. Both of these terms contain a large array of possible stochastic causes. For example, the time to beat parameter depends on mesoderm velocity, among other things, which can be expected to vary with such general effects as temperature or tissue geometry. Percent beating could also be affected by such things as Chapter 9. Conclusion 188 variability in surgical damage. Though it would be quite difficult to control these completely in an experiment, any possible elimination of noise, or at least isolation of separate noise effects, would help in testing the deterministic model. Perhaps one of the weaknesses of the present model is the sensitivity of the time-to-beat to the concentration of E. If the time-to-beat parameter can be better controlled, and the E gradient manipulated experimentally, it might become clear whether or not the present model of the inductive effects is realisitic. The best opportunity for further theoretical work, without more experiment, lies in expanding on the work of Chapter Four. This includes development of new methods for analysis of peak shape and evaluation of pattern order; comparison of the present results to those with other reaction-diffusion models; study of clustering and anti-clustering in the patterns; and investigation of what implications these partly-ordered patterns have on the general question of noise-reduction and reproducibility in reaction-diffusion patterns. All of these are potentially rich areas of study. Overall, I think this work shows the advantages of close collaboration between experimentalists and theorists. The model could never have become as detailed as it is, with its resulting predictions, without direct contact with experimentalists who wished to help to develop a dynamical theory. Much published biological evidence, commonly dealing with characterization of specific chemicals, or studies of their general effects, does not bear directly on the dynamics of a given process. The types of information necessary for developing a dynamical model have been readily available in this collaboration. For Chapter 9. Conclusion 189 example, the accurate measurements of tissue velocities and of the extent of spatial patterns has been indispensable. Such detail is necessary, I think, to convince experimentalists that dynamical theories, such as reaction-diffusion, are important in explaining developmental events. Most experimentalists study specific classes of molecules, and a common question of reaction-diffusion is whether or not a given molecule or class of molecules could be a Turing morphogen. Models become convincing when they can give some sort of answer to this question. Only through such communication will the dynamics of developmental processes be measured, and an understanding of them be attainable. What other pattern-forming events in vertebrates presently hold promise for similar interaction between theory and experiment? Study of the development of any single organ has potential, but additional factors are important: a) an experimental group currently active on the phenomenon that is amenable to collaboration with theoreticians; and b) initial data that suggest induction is not the whole mechanism. In the case of the axolotl, the long period between completion of induction and onset of heartbeat afforded a window in the developmental sequence for study of post-inductive events. This window is not long enough in other amphibians, such as Xenopus laevis or Taricha torosa, to easily study post-inductive differentiation. There are also cases of repeated structures in vertebrates which would be worth pursuing with reaction-diffusion theory. One example found in mammals is the interspecies variability in number of mammary glands (e.g. mice versus cattle versus humans, see Gilbert 1988). This would be worth a literature search to ascertain the available body of data. Another case, Chapter 9. Conclusion 190 pertinent to all vertebrates, is the sequential formation of somites, in which the dorsal mesoderm becomes segmented. The precise control of the final number of somites within a species might possibly proceed by a reaction-diffusion mechanism. Bibliography Antin, P.B., Taylor, R.G., and Yatskievych, T. (1994) Precardiac mesoderm is specified during gastrulation in quail. Dev. Dyn. 200, 144-154. Armstrong, J.B. (1989) A Turing model to explain heart development. AxolotlNewsletter 18, 23-25. Ascher, U.M., Ruuth, S.J., and Wetton, B.T.R. (1995). Implicit-explicit methods for time-dependent PDE's. SIAMJ. of Numerical Analysis 32, no. 3 (in press). Babyolantz, A., and Hiernaux, J. (1975) Models for cell differentiation and generation of polarity in diffusion-governed morphogenetic fields. Bull. Math. Biol. 37, 637-657. Bordzilovskaya, N.P. and Dettlaff, T.A. (1979) Table of stages of the normal development of axolotl embryos and the prognostication of timing of successive developmental stages at various temperatures. Axolotl Newsletter 7, 2-22. Bordzilovskaya, N.P., Dettlaff, T.A., Duhon, S.T., and Malacinski, G.M. (1989). Developmental stage series of axolotl embryos. In Developmental Biology of the Axolotl, Armstrong, J.B. and Malacinski, G.M. (eds). New York: Oxford University Press, pp. 201-219. Boveri, T. (1904). Ergebnisse tiber die Konstitution der chromatischen Substanz des Zellkerns. Jena: Gustav Fischer. Boveri, T. (1910). Die Potenzen derAscara-Blastomeren bei abgeanderter Furchung, zugleich ein Beitrag zur Frage qualitativ-ungleicher Chromosomen-Teilung. In Festschrift fur Richard Hertwig, vol. 3, Jena: Gustav Fischer, p.131. Boyce, W.E., and DiPrima, R.C. (1986). Elementary Differential Equations and Boundary Value Problems. 4th edition, New York: John Wiley and Sons. Burden, R.L., Faires, J.D., and Reynolds, A.C. (1981). Numerical Analysis. 2nd edition, Boston: Prindle, Weber, Schmidt. Butros, J. (1965) Action of heart and liver RNA on the differentiation of segments of chick blastoderms. J. Embryol. Exp. Morphol. 13,119-128. Casimir, H.B.G. (1945). On Onsager's principle of microscopic reversibility. Rev. mod. Phys. 17, 343-350. Bibliography 192 Caveney, S. (1985). The role of gap junctions in development. Ann. Rev. Physiol. 47: 319-35. Child, CM. (1915). Individuality in Organisms, Chicago University Press. Clark, P.J., and Evans, F.C. (1954). Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35, 445-453. Claxton, J.H. (1964). The determination of patterns with special reference to that of the central primary skin follicles in sheep. J. Theor. Biol. 7, 302. Copenhaver, W.M. (1926) Experiments on the development of the heart of Amblystoma punctatum. J. Exp. Zool. 43, 321-371. Crank, J., and Nicolson, P. (1947). A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Proc. Camb. Phil. Soc. 43, 50. Davis, L.A. and Lemanski, L.F. (1983) Inductive properties of a factor produced by endoderm. J. Cell Biol. 97, 58a (abstract). Davis, L.A., and Lemanski, L.F. (1987) Induction of myofibrillogenesis in cardiac lethal mutant axolotl hearts rescued by RNA derived from normal endoderm. Development 99, 145-154. Deshpande, A.K., Jakowlew, S.B., Arnold, H.-H., Crawford, P.A., and Siddiqui, M.A.Q. (1977) A novel RNA affecting embryonic gene functions in early chick blastoderm./. Biol. Chem. 252, 6521-6527. Deshpande, A.K, and Siddiqui, M.A.Q. (1977) A reexamination of heart muscle differentiation in the post nodal piece of chick blastoderm mediated by exogenous RNA. Dev. Biol. 58, 230-247. Desjardins, P., L'Abbe, D., Lang, B.F., and Morais, R. (1989) Putative chicken "muscle-specific 7S RNA" is related to the mitochondrial ATPase 6 gene. J. Molec. Biol. 207, 625-629. Ding, D., and Lipshitz, H.D. (1993). Localized RNAs and their functions. BioEssays 15, 651-658. Driever, W., and Niisslein-Volhard, C. (1988). A gradient of bicoid protein in Drosophila embryos. The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner. Cell 54, 83-93, 95-104. Easton, H.S., Armstrong, J.B., and Smith, S.C. (1994) Heart specification in the Mexican axolotl (Ambystoma mexicanum). Dev. Dyn. 200, 313-320. Bibliography 193 Edelstein-Keshet, L. (1988). Mathematical Models in Biology. New York: Random House. Eichele, G., and Thaller, C. (1987). Characterization of concentration gradients of a morphogenetically active retinoid in the chick limb bud. J. Cell Biol. 105, 1917-1923. Ekman, G. (1921) Experimented Beitrage zur Entwicklung des Bombinator-Herzens. Ofversikt af Finska Vetenskaps-Societetens Forhandlingar 63,1-37. Epstein, H.F. (1990) Myofibril assembly. Bioessays 12, 447-448. Erginel-Unaltuna, N., Dube, D.K., and Lemanski, L.F. (1993) Immunohistochemical studies of a unique protein from axolotl. Mol. Biol. Cell 4, Supplement, 143a (abstract). Ermentrout, B. (1991). Stripes or spots? Nonlinear effects in bifurcation of reaction-diffusion equations on the square. Proc. R. Soc. Lond. A 434, 413-417. Fransen, M.E., and Lemanski, L.F. (1988) Myocardial cell relationships during morphogenesis in normal and cardiac lethal mutant axolotls, Ambystoma mexicanum. Amer. J. Anat. 183, 245-257. Fuldner, R.A., Lim, S.-S., Greaser, M.L., and Lemanski, L.F. (1984) Accumulation and localization of troponin-T in developing hearts of Ambystoma mexicanum. J. Embryol. Exp. Morphol. 84, 1-17. Fullilove, S.L. (1970) Heart induction: Distribution of active factors in newt endoderm. J. Exp. Zool. 175, 323-326. Gierer, A , and Meinhardt, H. (1972) A theory of biological pattern formation. Kybernetik 22,30-39. Gilbert, S.F. (1988). Developmental Biology. 2nd edition, Sunderland, Mass. USA: Sinauer Associates. Glansdorff, P., and Prigogine, I. (1971). Thermodynamic Theory of Structure, Stability, and Fluctuations. New York: Wiley-Interscience. Graper, L. (1907) Untersuchungen uber die Herzbildung der Vogel. Roux'Arch. Entw. Mech. 24, 375-410. Gray, P., and Scott, S.K. (1983). Autocatalytic reactions in the isothermal, continuous stirred tank reactor. Isolas and other forms of multistability. Chem. Eng. Sci. 38, 29-43. Haldane, J.B.S. (1949). What Is Life? London: Alcuin Press. Bibliography 194 Harrison, L.G. (1993) Kinetic Theory of Living Pattern. New York and Cambridge: Cambridge University Press. Harrison, L.G. (1994). Kinetic theory of living pattern. Endeavour 18, 130-136. Harrison, L.G., Snell, J., and Verdi, R. (1984). Turing's model and pattern adjustment after temperature shock, with application to Acetabularia whorls. J. Theor. Biol. 106, 59-78. Harrison, L.G., Graham, K.T., and Lakowski, B.C. (1988) Calcium localization during Acetabularia whorl formation: evidence supporting a two-stage hierarchical mechanism. Development 104, 255-262. Harrison, L.G. and Kolar, M. (1988) Coupling between reaction-diffusion prepattern and expressed morphogenesis, applied to desmids and dasyclads. J. Theor. Biol. 130, 493-515. Harrison, R.G. (1969) Harrison stages and description of the normal development of the spotted salamander, Ambystoma punctatum (Linn.). In Organization and Development of the Embryo.Wilens, S. (ed.). New Haven, Ct., USA: Yale University Press, pp 44-66. Herschkowitz-Kaufman, M. (1975). Bifurcation analysis of nonlinear reaction-diffusion equations II. Steady-state solutions and comparison with numerical simulations. Bull. Math. Biol. 37, 589-636. Hill, C.S., and Lemanski, L.F. (1979) Morphological studies on cardiac lethal mutant salamander hearts in organ cultures. J. Exp. Zool. 209, 1-20. Holloway, D.M., Harrison, L.G., and Armstrong, J.B. (1994) Computations of post-inductive dynamics in axolotl heart formation. Dev. Dyn. 200, 242-256. Humphrey, R.R. (1968). A genetically determined absence of heart function in embryos of the Mexican axolotl {Ambystoma mexicanum). Anat. Rec. 160, 475 (abstract). Humphrey, R.R. (1972) Genetic and experimental studies on a mutant gene (c) determining absence of heart action in embryos of the Mexican axolotl {Ambystoma mexicanum). Dev. Biol. 27, 365-375. Hunding, A. (1987). Bifurcations in Turing systems of the second kind may explain blastula cleavage plane orientation. / . Math. Biol. 25, 109-122. Hunding, A. (1989). Turing patterns of the second kind simulated on supercomputers in three curvilinear coordinates and time. In Cell to Cell Signalling: From Experiments to Theoretical Models, Goldbeter, A. (ed.), London: Academic Press, pp. 229-236. Bibliography 195 Huxley, J.S., and De Beer, G.R. (1934). Elements of Experimental Embryology. Cambridge University Press. Jacobson, A.G. (1960) Influences of ectoderm and endoderm on heart differentiation in the newt. Dev. Biol. 2, 138-154. Jacobson, A.G. (1961) Heart determination in the newt. J. Exp. Zool. 146,139-151. Jacobson, A.G., and Duncan, J.T. (1968) Heart induction in salamanders. J. Exp. Zool. 167, 79-103. Justus, J.T. (1978) The cardiac mutant: an overview. Amer. Zool. 18, 321-326. Kirk, G.S., Raven, J.E., and Schofield, M. (1983). The Presocratic Philosophers. 2 ed. Cambridge University Press. Kulikowski, R.R., and Manasek, F.J. (1978) The cardiac lethal mutant ofAmbystoma mexicanum: A re-examination. Amer. Zool. 18, 349-358. Lacalli, T.C. (1981). Dissipative structures and morphogenetic pattern in unicellular algae. Phil. Trans. R. Soc. Lond. B 294, 547-588. Lacalli, T.C, and Harrison, L.G. (1978a). The regulatory capacity of Turing's model for morphogenesis, with application to slime moulds. / . Theor. Biol. 70, 273-295. Lacalli, T.C, and Harrison, L.G. (1978b). Development of ordered arrays of cell wall pores in Desmids: a nucleation model. J. Theor. Biol. 74, 109-138. Lacalli, T.C, and Harrison, L.G. (1979). Turing's conditions and the analysis of morphogenetic models. J. Theor. Biol. 76, 419-436. Lacalli, T.C, Wilkinson, D.A., and Harrison, L.G. (1988). Theoretical aspects of stripe formation in relation to Drosophila segmentation. Development 104, 105-113. LaFrance, S.M., Dube, D.K., Nakatsugawa, N., Erginel-Unaltuna, N., Capone, R., Lemanski, S.F., and Lemanski, L.F. (1993) Partial characterization of an RNA which promotes myofibrillogenesis in cardiac mutant axolotl hearts. Mol. Biol. Cell 4, Supplement, 143a (abstract). LaFrance, S.M., Fransen, M.E., Erginel-Unaltuna, N., Dube, D.K., Robertson, D.R., Stefanu, C , Ray, T.K., and Lemanski, L.F. (1994) RNA from normal anterior endoderm/mesoderm-conditioned medium stimulates myofibrillogenesis in developing mutant axolotl hearts. Cell, and Mol. Biol. Res. 39, 547-560. Lavenda, B.H. (1978). Thermodynamics of Irreversible Processes. 1993 Publication: Mineola , N.Y. USA: Dover. Bibliography 196 Lawson, J.D., and Morris, J.Ll. (1977). A note on the efficient implementation of splitting methods in two space variables. BIT17, 492-493. Lee, K-J., McCormick, W.D., Pearson, J.E., and Swinney, H.L. (1994). Experimental observation of self-replicating spots in a reaction-diffusion system. Nature 369,215-218. Lemanski, L.F. (1973) Morphology of developing heart in cardiac lethal mutant Mexican axolotls (Ambystoma mexicanum). Dev. Biol. 33, 312-333. Lemanski, L.F. (1978) Morphological, biochemical and immunohistochemical studies on heart development in cardiac mutant axolotls, Ambystoma mexicanum. Amer. Zool. 18, 327-348. Lemanski, L.F., and Fuldner, R.A. (1977) Immunofluorescent studies for myosin, tropomyosin, and a-actinin in developing hearts of normal and cardiac lethal mutant salamanders {Ambystoma mexicanum). J. Cell Biol. 75, 327a (abstract). Lemanski, L.F., Joseph, X., and Iyengar, M.R., (1975) Quantitation by radioimmunoassay of absolute amounts of myosin in embryonic hearts of cardiac lethal mutant axolots {Ambystoma mexicanum). J. Cell Biol. 67, 239a (abstract). Lemanski, L.F., Marx, B.S., and Hill, C.S. (1977) Evidence for abnormal heart induction in cardiac-mutant salamanders {Ambystoma mexicanum). Science 196, 894-896. Lemanski, L.F., Paulson, D.J., and Hill, C.S. (1979) Normal anterior endoderm corrects the heart defect in cardiac mutant salamanders {Ambystoma mexicanum). Science 204, 860-862. Lisman, J.E. (1985). A mechanism for memory storage insensitive to molecular turnover: a bistable autophosphorylating kinase. Proc. Nat. Acad. Sci. U.SA. 82, 3055-3057. Lyons, M.J. (1991). Mathematical Models for Biological Pattern Formation in Two Dimensions. Ph.D. thesis, University of British Columbia, Canada. Lyons, M.J., Harrison, L.G., Lakowski, B.C., and Lacalli, T.C. (1990). Reaction diffusion modelling of biological pattern formation: application to the embryogenesis ofDrosophila melanogaster. Can. J. Physics 68, 772-777. Lyons, M.J., and Harrison, L.G. (1991). A class of reaction-diffusion mechanisms which preferentially select striped patterns. Chem. Phys. Lett. 183, 158-164. Lyons, M.J., and Harrison, L.G. (1992). Stripe selection: an intrinsic property of some pattern-forming models with nonlinear dynamics. Dev. Dyn. 195, 201-215. Bibliography 197 Mangiacapra, F., Fransen, M.E., and Lemanski, L.F. (1993) Activin A and transforming growth factor-p stumulate heart formation in axolotls but do not rescue cardiac lethal mutants. Mol. Biol. Cell 4, Supplement, 143a (abstract). Markovics, J., Glass, L., and Maul, G.G. (1974). Pore patterns on nuclear membranes. Exptl. CellRes. 85, 443-451. Maynard Smith, J. (1968). Mathematical Ideas in Biology. Cambridge University Press. McQuarrie, D.A. (1976). Statistical Mechanics. New York: Harper Collins. Meinhardt, H. (1982). Models of Biological Pattern Formation. London: Academic Press. Meinhardt, H. (1989). Models for positional signalling with application to the dorso-ventral patterning of insects and segregation into different cell types. Development [Suppl. J, 169-180. Murray, J.D. (1980). A pattern formation mechanism and its application to mammalian coat markings. In Lecture Notes in Biomathematics 39. New York: Springer, pp. 360-399. Murray, J.D. (1981). On pattern formation mechanisms for lepidopteran wing patterns and mammalian coat markings. Phil. Trans. Roy. Soc. Lond. B 295, 473-496. Murray, J.D. (1989). Mathematical Biology. New York: Springer-Verlag. Muslin, A.J., and Williams, L.T. (1991) Well-defined growth factors promote cardiac development in axolotl mesodermal explants. Development 112, 1095-1101. Nicolis, G. (1971). Stability and dissipative structures in open systems far from equilibrium. Adv. chem. Phys. 19, 209-324. Niu, M.C., and Deshpande, A.K. (1973) The development of tubular heart in RNA-treated post-nodal pieces of chick blastoderm. / . Embryol. Exp. Morphol. 29, 485-501. Onsager, L. (1931). Reciprocal relations in irreversible processes I. Phys. Rev. 37, 405-426; Reciprocal relations in irreversible processes II. Phys. Rev. 38, 2265-2279. Panayotou, G., and Waterfield, M. (1993) The assembly of signalling complexes by receptor tyrosine kinases. BioEssays 15, 171-177. Pauling, L. (1987). Schrodinger's contribution to chemistry and biology. In Schrodinger: Centenary celebration of a polymath, Kilmister, C.W. (ed). Cambridge University Press. Bibliography 198 Peaceman, D.W., and Rachford, H.H. Jr. (1955) The numerical solution of parabolic and elliptic differential equations. J. Soc. Indust. Appl. Math. 3, no. 1, 365-375. Pearson, J.E. (1993). Complex patterns in a simple system. Science 261, 189-192. Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T. (1988). Numerical Recipes, The Art of Scientific Computing. Cambridge University Press. Prigogine, I. (1946). Moderation and irreversible transformations of open systems. Bull. Acad. Roy. Belg. CI. Sci. 31, 600-606. Prigogine, I. (1967a). Dissipative structures in chemical systems. In Fast Reactions and Primary Processes in Chemical Kinetics, Claesson, S. (ed). New York: Interscience, pp. 371-82. Prigogine, I. (1967b). Introduction to Thermodynamics of Irreversible Processes, 3 ed. New York: Wiley-Interscience. Prigogine, I., and Lefever, R. (1968). Symmetry-breaking instabilities in dissipative systems. II. / . Chem. Phys. 48, 1695-1700. Rashevsky, N. (1940). An approach to the mathematical biophysics of biological self-regulation and of cell polarity. Bull. Math. Biophys. 2, 15-25. Reichl, L.E. (1980). A Modern Course in Statistical Physics. Austin, Tx. USA: University of Texas Press. Ruuth, S.J. (submitted). Implicit-explicit methods for reaction-diffusion problems in pattern formation. J. Math. Biol. Safranyos, R.G.A., and Caveney, S. (1984). Rates of diffusion of fluorescent molecules via cell-to-cell membrane channels in a developing tissue. J. Cell. Biol. 100: 736-747. Sander, K. (1994). Of gradients and genes: developmental concepts of Theodor Boveri and his students. Roux'sArch. Dev. Biol. 203, 295-297. Sater, A.K., and Jacobson, A.G. (1989) The specification of heart mesoderm occurs during gastrulation in Xenopus laevis. Development 105, 821-830. Sater, A.K. and Jacobson, A.G. (1990) The restriction of the heart morphogenetic field in Xenopus laevis. Dev. Biol. 140, 328-336. Schnackenberg, J. (1979). Simple chemical reaction systems with limit cycle behaviour. J. Theor. Biol. 81, 389-400. Segel, I.H. (1975). Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady State Enzyme Kinetics. New York: Wiley-Interscience. Bibliography 199 Shen, P.-S., and Lemanski, L.F. (1986) Immunofluorescent studies for desmin, vimentin and vinculin in developing hearts of normal and cardiac mutant Mexican axolotls, Ambystoma mexicanum. J. Cell Biol. 103, 123a. Shen, P.-S., and Lemanski, L.F. (1989) Immunofluorescent, immunogold, and electrophoretic studies for desmin in embryonic hearts of normal and cardiac mutant Mexican axolotls, Ambystoma mexicanum. J. Morphol. 201, 243-252. Slack, J.M.W. (1983) From Egg to Embryo. Determinative Events in Early Development. Cambridge: Cambridge University Press. Smith, S.C. (1990) Control of heart development in the Mexican axolotl (Ambystoma mexicanum). Ph.D. thesis, University of Ottawa, Canada. Smith, S.C. and Armstrong, J.B. (1990) Heart induction in wild-type and cardiac mutant axolotls (Ambystoma mexicanum). J. Exp. Zool. 254, 48-54. Smith, S.C, and Armstrong, J.B. (1991) Heart development in normal and cardiac-lethal mutant axolotls: a model for the control of vertebrate cardiogenesis. Differentiation 47, 129-34. Smith, S.C, and Armstrong, J.B. (1993) Reaction-diffusion control of heart development: Evidence for activation and inhibition in precardiac mesoderm. Dev. Biol. 160, 535-542. Stainier, D.Y.R., and Fishman, M.C. (1992) Patterning the zebrafish heart tube: acquisition of anteroposterior polarity. Dev. Biol. 153, 91-101. Stainier, D.Y.R., and Fishman, M.C. (1993) Cardiac morphogenesis in the zebrafish, patterning the heart tube along the anteroposterior axis. In Molecular Basis of Morphogenesis, Bernfield, M. (ed). New York: Wiley-Liss, pp. 79-91. Starr, CM., Diaz, J.G., and Lemanski, L.F. (1989) Analysis of actin and tropomyosin in hearts of cardiac mutant axolotls by two-dimensional gel electrophoresis, western blots, and immunofluorescent microscopy./. Morphol. 200,1-10. Strikwerda, J.C (1989). Finite Difference Schemes and Partial Differential Equations. Pacific Grove, Ca. USA: Wadsworth and Brooks/Cole. Sugi, Y., and Lough, J. (1994) Anterior endoderm is a specific effector of terminal cardiac myocyte differentiation of cells from the embryonic heart forming region. Dev. Dyn. 200,155-162. Takagi, I. (1982). A priori estimates for stationary solutions of an activator-inhibitor model due to Gierer and Meinhardt. Tohoku MathematicalJournal, 113-132. Takagi, I. (1986). Point-condensation for a reaction-diffusion system. J. Diff. Eqn. 61: 208-249. Bibliography 200 Thomas, D. (1975). Artificial enzyme membranes, transport, memory, and oscillatory phenomena. In Analysis and Control of Immobilized Enzyme Systems, Thomas, D., and Kernevez, J.-P. (eds.). New York: Springer, pp. 115-150. Tickle, C , Alberts, B.M., Wolpert, L., and Lee, J. (1982). Local application of retinoic acid to the limb bud mimics the action of the polarizing region. Nature 296, 564-565. Turing, A.M. (1952) The chemical basis of morphogenesis. Phil. Trans. Roy. Soc. Lond. B 237,37-72. Verocay (1905) Multiplicitas cordis (Heptocardia) bei einem Huhn. Verhandl. der Deutschen pathol. Ges., Erg. Heft, 16, 192-198 Wigglesworth, V.B. (1940) Local and general factors in the development of "pattern" in Rhodniusprolixus (Hemiptera). J. Exp. Biol. 17, 180-200. Wisniewski, S., Staniszewski, B., and Szymanik, R. (1976). Thermodynamics of Nonequilibrium Processes. Hingham, Mass. USA: D. Reidel. Yutzey, K.E., Rhee, J.T., and Bader, D. (1994) Expression of the atrial-specific myosin heavy chain AMHC1 and the establishment of anteroposterior polarity in the developing chicken heart. Development 120, 871-883. 


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