A Wave-Pinning Mechanism for Eukaryotic Cell Polarization Based on Rho GTPase Dynamics by Alexandra Jilkine A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2009 c© Alexandra Jilkine 2009 Abstract In response to chemical stimulation, many eukaryotic cells are able to sense the direction of the stimulus and initiate movement. To do so, the cell must break symmetry and develop a front and back in a process known as polarization. During polarization, members of the Rho GTPase family (Cdc42, Rac, and Rho) are recruited to the plasma membrane and localize to form a front and a back of the polarizing cell. These proteins exist in both active forms (on the inner surface of the membrane of the cell), and inactive forms (in the cytosol). In earlier work, I have shown that the property of membrane-cytosol interconversion, together with appropriate feedbacks, endows the Rho proteins with the ability to initiate cell polarity, resulting in a high Cdc42/Rac region, which will become the front, and a high Rho region, which will become the back of the cell. Here I show that this property of polarizability can be explained using a simplified model system comprising of a single active/inactive protein pair with positive feedback. In this model, a travelling wave of GTPase activation is initiated at one end of the domain, moves across the cell, and eventually stops inside the domain, resulting in a stable polar distribution. The key requirements for the mechanism to work include conservation of total amount of protein, a sufficiently large difference in diffusion rates of the two forms, and nonlinear positive feedback that allows for multiple homogeneous steady states to exist. Using singular perturbation theory, I explain the mathematical basis of wave-pinning behaviour, and discuss its biological and mathematical implications. I show that this mechanism for generating a chemical pattern is distinct from Turing pattern forma- ii Abstract tion. I also analyze the transition from a spatially heterogeneous solution to a spatially homogeneous solution as the diffusion coefficient of the active form is increased, and show the existence of other unstable stationary wavefronts. Finally, I argue that this wave-pinning mechanism can account for a number of features of cell polarization such as spatial amplification, maintenance of polarity, and the sensitivity to new stimuli that is typical of polarization of eukaryotic cells. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Co-authorship Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Actin-based Cell Motility . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 Regulation of Cell Motility by Rho GTPases . . . . . . . . . . . 3 1.1.3 Cell Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Model Organisms for Cell Polarization and Motility . . . . . . . 6 1.2 Understanding Cell Polarization in Theoretical Terms . . . . . . . . . . 9 1.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Rho GTPases and the Formulation of the Wave-Pinning Model . . 14 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.1 Guanosine Exchange Cycle of the Rho Proteins . . . . . . . . . 15 iv Table of Contents 2.1.2 A Model for Spatial Segregation of Rho GTPases Based on In- hibitory Crosstalk . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Model Reduction to a Single Rho GTPase . . . . . . . . . . . . . . . . 24 2.2.1 Modelling Assumptions and Simplifications . . . . . . . . . . . . 24 2.2.2 The Well-Mixed System . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.3 The Spatially Distributed System . . . . . . . . . . . . . . . . . 27 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.1 Input Stimuli . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Parameter Values . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.3 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3 Comparison of Cell Polarization Models . . . . . . . . . . . . . . . . . 35 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 A Survey of Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.1 Model Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.2 Turing Type Models . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.3 Adaptation Models . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.4 Wave-Based Models . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.5 Stochastic Models . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.6 Integrating Different Types of Models . . . . . . . . . . . . . . . 51 3.3 A Comparison of Models . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.1 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.2 Time to Polarize . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3.3 Model Comparison Summary . . . . . . . . . . . . . . . . . . . . 64 3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4 Asymptotic Analysis of the Wavepinning Mechanism . . . . . . . . . 67 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 v Table of Contents 4.1.1 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Intuitive Explanation for Wave-Pinning . . . . . . . . . . . . . . . . . . 69 4.2.1 The Single Variable Case (fixed v) . . . . . . . . . . . . . . . . . 69 4.2.2 Wave-Pinning in the Full u, v System: . . . . . . . . . . . . . . . 71 4.3 Matched Asymptotic Expansion in 1D . . . . . . . . . . . . . . . . . . . 73 4.3.1 Scaling the Equations . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.2 Asymptotic Analysis to Leading Order . . . . . . . . . . . . . . 76 4.3.3 Travelling Wave Pins by Depletion of Inactive Protein . . . . . . 80 4.3.4 Asymptotic Analysis to Higher Order . . . . . . . . . . . . . . . 82 4.4 Example with Cubic Kinetic Function . . . . . . . . . . . . . . . . . . . 87 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5 Wave-pinning in 2D: Numerics and Asymptotic Analysis . . . . . . 98 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.2 Asymptotic Analysis in 2D . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.1 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.2 Local Orthogonal Coordinate System . . . . . . . . . . . . . . . 101 5.2.3 Outer Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.2.4 Inner Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2.5 Matching Inner and Outer Solutions . . . . . . . . . . . . . . . . 104 5.2.6 u and v to Order ² . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.2.7 Front Speed to Order ² . . . . . . . . . . . . . . . . . . . . . . . 108 5.3 Numerical Results in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.4 Similarity to the Constrained Allen Cahn Model . . . . . . . . . . . . . 113 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6 Bifurcation Analysis of the Wave-pinning System . . . . . . . . . . . 118 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.2 Model Equations and Stability of the Homogeneous Steady States . . . 119 vi Table of Contents 6.3 Numerical Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.3.1 Pseudo-Arclength Continuation . . . . . . . . . . . . . . . . . . 122 6.3.2 Numerical Continuation on the Full System . . . . . . . . . . . 124 6.4 Reduction of the Steady State System to an Integral System . . . . . . 125 6.4.1 Numerical Continuation on the Time Map System . . . . . . . . 131 6.4.2 Two Parameter Bifurcation Diagrams . . . . . . . . . . . . . . . 134 6.4.3 Full Bifurcation Diagram in the limit D →∞ . . . . . . . . . . 135 6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 Appendix A Numerical Simulation Details . . . . . . . . . . . . . . . . . 158 A.1 Stimuli repertoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 Appendix B Mathematical Background . . . . . . . . . . . . . . . . . . . 160 B.1 Conditions for Diffusion-Driven Instability . . . . . . . . . . . . . . . . 160 B.2 Stability of the Homogeneous States for the Wave-Pinning System . . . 163 B.3 Explicit front speed formula for the cubic kinetics function . . . . . . . 165 B.4 Solvability Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 B.5 Critical Epsilon Calculation for Emergence of Nonhomogeneous Solution 168 B.6 Integral Evaluation for the Time Map System . . . . . . . . . . . . . . 170 vii List of Tables 3.1 Features of polarity explained by various models. . . . . . . . . . . . . . 66 viii List of Figures 1.1 A model for actin filament assembly and disassembly at the leading edge of a cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Downstream effectors of Cdc42, Rac, and Rho. . . . . . . . . . . . . . . 5 1.3 Examples of polarized cells. . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Guanosine exchange cycle of the Rho GTPases . . . . . . . . . . . . . . 17 2.2 Schematic diagram for Rho protein crosstalk based on the scheme in (28) 19 2.3 Geometry of the model in 1D . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Boundary conditions in 1D . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Numerical solutions of system (2.1) . . . . . . . . . . . . . . . . . . . . . 22 2.6 The reaction terms of equation (2.8) as functions of u for a fixed v. . . . 27 2.7 A schematic diagram of stable profiles of u and v across the cell. . . . . 29 2.8 Polarization behaviour in response to six types of distinct stimuli for system (2.12) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.9 Response of the wave-pinning model in response to two pipette experiments. 33 3.1 Two ways of interpreting cell polarity models in one spatial dimension . 37 3.2 Examples of Turing mechanisms . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Perfect adaptation can only arise in two types of networks. . . . . . . . 42 3.4 Examples of adaptation models . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 Hysteresis in a bistable system . . . . . . . . . . . . . . . . . . . . . . . 47 3.6 Polarization behaviour in response to six types of distinct stimuli and/or initial conditions for wave-pinning system (2.12) . . . . . . . . . . . . . 53 ix List of Figures 3.7 Polarization behaviour in response to five types of distinct stimuli and/or initial conditions for ASDM system (3.5) . . . . . . . . . . . . . . . . . . 55 3.8 Polarization behaviour in response to five types of distinct stimuli and/or initial conditions for LEGI system (3.6) . . . . . . . . . . . . . . . . . . 59 3.9 Polarization behaviour in response to five types of distinct stimuli and/or initial conditions for system (3.7) . . . . . . . . . . . . . . . . . . . . . . 61 3.10 Comparison of sensitivity and persistence of the wave-pinning model, the LEGI model and the substrate-depletion model (ASDM) in response to uniform and successively stronger graded stimuli. . . . . . . . . . . . . . 65 4.1 An intuitive explanation of wavepinning based on the Maxwell condition 70 4.2 Wavepinning in 1D for system (4.10) with reaction kinetics (4.12) . . . . 73 4.3 The evolution of the steepest point φ through time and the solution of Eq. 4.70 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.4 The error between the numerical front position φnum obtained from nu- merically solving the full PDE system (4.10) with reaction kinetics (4.12) 90 4.5 The error between the numerical front position φnum, the zero order asymptotic order approximation (Eq. 4.89), and the two term asymp- totic approximation (Eq. 4.89 and 4.96) . . . . . . . . . . . . . . . . . . 95 4.6 The maximum error between the numerical front position φnum, the one term asymtptotic approximation (4.89), and the two term asymptotic approximation (4.89 and 4.96) for a range of ² values. . . . . . . . . . . 96 5.1 For a flat cell (top-down view), region of high Rho GTPase activity is separated from a region of low activity by a planar curve Γ(t). . . . . . 101 5.2 Solutions to (5.2) with cubic reaction term (5.43) on a rectangular do- main for various initial conditions. . . . . . . . . . . . . . . . . . . . . . 110 5.3 Evolution of system (5.2) with a Hill function reaction term (5.44) in response to a stimulus in k0 . . . . . . . . . . . . . . . . . . . . . . . . . 111 x List of Figures 5.4 Evolution of system (5.2) with a Hill function reaction term (5.44) on a rectangular domain from noisy initial conditions in 2D. . . . . . . . . . . 114 5.5 Evolution of system (5.2) with the cubic reaction term (5.43) on a rect- angular domain from noisy initial conditions in 2D. . . . . . . . . . . . . 115 6.1 Change in the nature of solutions to (6.1) as ² is increased. . . . . . . . 120 6.2 An illustration of the supplemental condition (6.4) in pseudo-arclength continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.3 Bifurcation diagrams with for system (6.8) using pseudo-arclength con- tinuation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.4 The solutions of (6.13) viewed in the u-ux phase-plane and the corre- sponding F (u,A,B) given by (6.19). . . . . . . . . . . . . . . . . . . . . 126 6.5 Bifurcation diagram showing the spatially heterogeneous branch in the time map system (6.25) obtained using pseudo-arclength continuation . 133 6.6 Two parameter continuation plots showing the region in which wave- pinning occurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.7 The full bifurcation diagram of the spatially heterogeneous solutions in the time map system (6.34) in the case of D =∞. . . . . . . . . . . . . 135 6.8 Two parameter bifurcation diagram in the time map system (6.34) in the case of D =∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 xi List of Abbreviations Arp2/3 actin-related protein 2/3 GDP guanosine diphospate GTP guanosine triphosphate GTPase GTP hydrolysing protein GEF guanine nucleotide exchange factor GAP GTPase-activating protein GDI guanine nucleotide dissociation inhibitor PI phosphatidylinosytol PI3K phosphatidylinosytol 3-kinase PIP2 phosphatidylinosytol-4,5-bisphosphate PIP3 phosphatidylinosytol-3,4,5-trisphosphate PIP5 kinase phosphatidylinosytol 4-phosphate 5-kinase xii Acknowledgements First and foremost, I would like to thank my supervisor, Leah Edelstein-Keshet, for her persistence, enthusiasm, and constant encouragement throughout the years I’ve been at UBC. I am proud to be her student, and am indebted to her for taking the time to mentor me. I am grateful to Yoichiro Mori for all his help with the more technical aspects of this work and his sage advice. I would like to express my sincere appreciation to my supervisory committee: Dan Coombs, Eric Cytrynbaum, Eirikur Palsson, and Nick Hill. Thank you to Wayne Nagata and Michael Gold for being on my examination committee. I would also like to thank the current and former members of the UBC math biology group for creating a stimulating research environment. Special thanks goes to Alan Lindsay for his help with numerical continuation methods. And, of course, thank you to my friends at the IAM for making grad school a fun experience. Finally, I wish to thank my family for their continual love and support. Funding for this research was provided by a subcontract to Leah Edelstein-Keshet from a NIH grant R01 GM086882 to Anders Carlsson (Washington University, St.Louis, MO), an NSERC CGS fellowship, an NSERC discovery grant to Leah Edelstein-Keshet, and a MITACS team project grant. I would also like to acknowledge the Department of Mathematics and the Institute of Applied Math for providing research facilities. xiii Co-authorship Statement The model of Cdc42-Rac-Rho dynamics presented in Chapter 2 was developed by myself, in conjunction with A.F.M. Marée (Utrecht), and Leah Edelstein-Keshet (UBC) and published in (44). I proposed the reduced wave-pinning model and carried out all the numerical simulations presented in this thesis. Yoichiro Mori (University of Minnesota) played a large role in the analysis of the wave-pinning model, presented in Chapters 4-6. He proposed an explanation for the wave-pinning results based on the Maxwell condition, and carried out the zero order asymptotic calculations in 1D, which were independently verified by me. Parts of that work has previously been published in (76). In that publication I performed all the numerical simulations and participated in the write-up and submission of the paper. The higher order asymptotic calculations in 1 and 2D presented in Chapters 4 and 5 were carried out by me with Yoichiro’s guidance. Yoichiro Mori and Leah Edelstein-Keshet proposed the idea of a “time-map” to study the bifurcation structure of the wave-pinning thesis presented in Chapter 6. I developed the code for the numerical continuation results presented in Chapter 6 ( Yoichiro Mori provided code for integral evaluation used to produce Figure 6.7). My PhD. supervisor, Leah Edelstein-Keshet, has been closely involved in all the work presented in this thesis. She has provided valuable feedback on the results, suggested various computational tests, assisted in solving analytical and computational difficulties, and has tirelessly edited my writing. All the numerical results and figures (except reproductions of experimental data) in this thesis are my own work, and I am entirely responsible for any error or inaccuracies. xiv Chapter 1 Introduction and Motivation 1.1 Introduction 1.1.1 Actin-based Cell Motility The ability to move is important to many eukaryotic organisms, from large animals to single-celled amoeba. Single-celled organisms, such as the cellular slime mold Dic- tyostelium, migrate in search of food. In multicellular organisms, cell motility is vital for embryogenesis, wound repair, angiogenesis and immune system response. In developing embryos, primordial germ cells need to migrate long distances to the sites of future ovaries and testes. In mature individuals, cells like macrophages and neutrophils crawl to sites of infection, while fibroblasts migrate to remodel connective tissue at wound sites. Cancerous cells develop the ability to invade surrounding tissue, becoming ma- lignant. Hence, understanding the mechanisms by which cells migrate is a fundamental research question in cell biology. In tissues, most cells crawl in three dimensions through extracellular matrix. In the laboratory, cell motility is most often studied on a 2D substrate. Cells crawl over solid surfaces by extending a pseudopod (‘false foot’) at the front of the cell, adhering to the underlying substrate, and retracting the tail end of the cell. Pseudopod extension is driven by assembly of the actin cytoskeleton, a highly dynamic network of protein polymers that provides the force necessary for migration (1; 93). By utilizing the energy from ATP hydrolysis, actin monomers (G-actin) polymerize into polar filaments (F- actin) that can add or subtract subunits at their ends. The reaction rate constants differ at the two ends of the polymer, with the so-called ‘barbed’ or plus end polymerizing 1 1.1. Introduction faster than the ‘pointed’ or minus end. This allows for net addition of subunits at the plus end, and net disassembly at the minus end for a range of monomer concentrations. As a result, even though the actin filament may maintain a constant length, there is a constant flux of subunits throughout the polymer (1). At the leading edge of some motile cells, a thin, broad membrane protrusion called the lamellipodium is filled with a dense, branched network of these actin filaments. This network is essentially two-dimensional, as the lamellipodium is very flat (on the order of 0.1-0.2 µm) (53). The cytoskeletal network is oriented with the plus ends pushing on the plasma membrane at the leading edge and the minus ends facing the cell body behind (13). (See Figure 1.1 for a current model of how actin filaments assemble and disassem- ble at the leading edge.) Actin filaments are anchored to adhesion complexes via the plasma membrane, allowing the force that is produced by the cytoskeletal network to change the shape of the membrane and be transmitted to the substrate (75; 91). The study of actin-driven cell motility has seen significant contributions from mathemati- cians and physicists, who quantified theories about the actin remodelling process and used thermodynamics to understand the nature and magnitude of the polymerization force (see, e.g. reviews in (11; 74)). Although actin filaments are too small to be visualized with a light microscope, they can be observed in live cells if labelled with a fluorescent probe, such as green fluorescent protein (GFP) or a fluorescent form of the the toxin phalloidin, which binds to actin filaments. Drugs such as cytochalasin and latruncullin bind to actin filaments and interfere with new assembly from monomers, leading to the disassembly of the actin network. Their application disrupts the cortical actin network on the timescale of seconds, revealing that the cytoskeletal network is undergoing rapid remodelling as the cell is moving. Cells use a multitude of proteins to regulate the assembly of the actin network: some proteins promote assembly, others stabilize filaments, still others depolymerize the filaments. For example, a protein complex called Arp2/3 binds to an actin filament, creating a nucleation site for a new actin filament, allowing a 2 1.1. Introduction branched actin meshwork to emerge (92). The Arp2/3 complex is incorporated into the actin network, and new filaments are capped rapidly, so that activated Arp2/3 complex must be supplied continuously to keep the network growing. Older filaments are depolymerized by the actin depolymerizing factor ADF/cofilin. New barbed ends are also created by uncapping capped filaments and severing existing filaments into multiple smaller ones. While actin polymerization provides a pushing force on the plasma membrane, myosins, a family of motor proteins that can move along actin filaments, provide a contractile force required for the retraction of the trailing edge of the moving cell. Bipolar myosin II filaments bind to actin filaments and cause them to slide past each other, causing the actomyosin bundle to contract (1). Contractile forces seem to be more important in adherent, slowly moving cells such as fibroblasts. Fast, gliding cells such as keratocytes can still translocate without myosin. In keratocytes, myosin activ- ity is important for establishing polarity and motility, but its inhibition does not have much effect on already moving cells (122). 1.1.2 Regulation of Cell Motility by Rho GTPases Rho family GTPases regulate the rates of actin filament polymerization, capping and branching, leading to structural reorganization of the cytoskeleton and cell motility (96; 98). Rho proteins are members of the GTP-binding proteins (G-proteins) superfamily, which are conserved in eukaryotic organisms from yeast and amoeba to mammals. G- proteins are an important class of cell messengers containing a GTP-binding domain that function as molecular switches. The protein must be bound to GTP in order to be active. The proteins hydrolyse GTP, thus returning to inactive form after some time. This hydrolysation of GTP is why G-proteins are also called GTPases. Upon hydrolysing GTP to GDP the protein undergoes a conformational change and becomes inactive. The proteins do not directly interact with actin, but rather relay external stimuli to the cytoskeleton by binding to so-called effector molecules, which then regulate 3 1.1. Introduction (a) (b) Figure 1.1: (a) A model for actin filament assembly and disassembly at the lamel- lipodium, showing how actin binding proteins control the rates of actin filament poly- merization, capping and branching. (b) Signaling pathways from receptors to Arp2/3 complex through Rho GTPases. Images taken from (92). 4 1.1. Introduction the cytoskeleton. Rho proteins are subdivided into three important classes: Cdc42, Rac and Rho. In their active forms Cdc42 and Rac both stimulate a protein family called WASP to acti- vate Arp2/3 and promote actin polymerization by nucleating new barbed ends. Several other pathways downstream of Cdc42/Rac contribute to formation of the branched actin structure found in the lamellipodium. GTP-bound Rho has a different set of targets from Cdc42. Active Rho promotes the phosphorylation of the myosin light chain (MLC) via its downstream target ROCK (Rho kinase). The phosphorylated MLC promotes the interaction of myosin filaments with actin filaments and produces contractile actin bundles such as stress fibers. See Figure 1.2 for details of how Rho GTPases regulate the cytoskeleton. Cdc42 Rac Rho PAKs LIM kinase cofilin WASP ARP2/3 IRSp53 WAVE/SCAR mDIA PIP5 kinase PIP2 ROCK (Rho kinase) gelsolin phophorylated MLC stress fibers/ increased actomyosin contractility Figure 1.2: Downstream effectors of Cdc42, Rac, and Rho compiled from the following review papers: (8; 67; 97; 98). Effectors that have functions unrelated to the cytoskeleton are not shown. Cofilin, Arp2/3, and gelsolin are actin-related proteins. Cofilin cleaves existing filaments, Arp2/3 allows nucle- ation of filaments, and gelsolin cuts filaments and caps their ends. 1.1.3 Cell Polarization Another direction of research in cell motility in which mathematical modeling had an impact concerns the chemical events that occur when cells commit to a direction in which to move. In migrating cells the front and back ends of the cell needs to be struc- turally different. Polarization is the process of building different structures with distinct 5 1.1. Introduction molecular components at each end. Polarization is triggered by external environmental cues, such as chemical gradients (88). There is rapid symmetry breaking in the distri- bution of multiple proteins and lipids: some substances become more concentrated (or activated) in what will become the front of the cell, and others are concentrated at the “back”. Similarly, there is a redistribution of substances from the cytosol to the plasma membrane. A Rho family GTPase called Cdc42, is considered to be a central regulator of polarity in eukaryotic cells. Rho family GTPases connect the assembly and disassembly of the actin cytoskeleton with membrane receptor signalling, transducing external chemotactic stimuli and amplifying internal signals (96; 98). Thus, Rho proteins are key players in polarity establishment, and understanding the role that they play in organizing cell polarity is the motivation for this thesis. In many cell types, polarization is upstream from cytoskeletal reorganization and motility (35; 48; 68; 113). Therefore we do not explicitly include the actin cytoskeleton in the models in this thesis, but note that feedback between the Rho proteins may be mediated through the cytoskeleton. 1.1.4 Model Organisms for Cell Polarization and Motility One classic model organism for cell polarization is the budding yeast Saccharomyces cerevisiae. In the course of the cell cycle, the yeast needs to transition from uniform to polarized growth in order to bud. Cdc42 gathers in an area of elevated concentration, called a polar cap (see (89) for a review). To establish a cap, Cdc42 molecules have to accumulate in a small region of the cell membrane and it is important that this area is defined with high spatial precision, as this site is used as a marker for the growth of a daughter cell during cell division. Polarization can occur in response to an internal signal (old bud scar), or spontaneously in the absence of pre-existing spatial cues. This transition to a polarized state is thought to reflect the action of positive and negative feedback loops that reinforce inequalities in the local concentrations of polarity factors, so that stochastic fluctuations are amplified into a single dominant asymmetry. Many 6 1.1. Introduction (a) (b) (c) (d) Figure 1.3: Examples of polarized cells. (a) A cap of Cdc42 self-assembles in budding yeast to mark the location where the next budding event will take place. (A vacuole inside the cell appears as a white circle) Scale bar= 2 µm. Reprinted with permission from (62). c©Elsevier 2007. (b) HL60 cells (neutrophil-like cell line) treated with a uniform concentration of chemoattractant for 2 min, fixed and stained. The cell develops a polarized morphology, with actin (red) forming a pseudopod at the front and RhoA (green) localizing to the the rear of the cell. Image from (9). (c) A keratocyte visualized with fluorescence and electron microscopy. In the fluorescence image actin is cyan and myosin is red. In the electron micrograph of a similar cell, the branched actin network is revealed. Reprinted with permission from (109). c©University Rockefeller Press. (d) Rac in moving fibroblasts. GFP-Rac shows the total amount of fluorescently labelled RAC in a cell. The FRET signal shows the active fraction of Rac in the cell. Reprinted with permission from (50). c©2000 AAAS. 7 1.1. Introduction proteins leading to yeast polarization appear to be conserved across a broad range of cell types, including migrating mammalian cells. Chemotaxis is defined as cell movement in a direction controlled by a gradient of a diffusible chemical. One common experimental model system for chemotactic cells are white blood cells called neutrophils, which crawl toward a source of bacterial in- fection by chemotaxis, detecting N-formylated peptides, which are only produced by prokaryotic cells. Stimulating human neutrophils with a homogeneous concentration of chemoattractant induces ruffles to form all over the cell surface, which then consolidate into a single pseudopod within minutes (125). Visualizing these cells with fluorescent probes reveals that some molecules localize to the front of the cell (phosphoinositide lipids PIP2 and PIP3, F-actin, and the Rho GTPases Cdc42 and Rac), while others go to the emerging back (RhoA and myosin) (120; 121). There is evidence that these front and back localizing molecules mutually inhibit each other, supporting the idea that cell polarization is a self-organizing mechanism that emerges as a result of feedback loops between various polarity factors (114). Another commonly used organism in the study of chemotaxis and gradient sensing is the amoeba Dictyostelium discoideum. Under starvation conditions the amoeba secretes cyclic AMP (cAMP), attracting other amoebas to migrate toward the chemical source and begin aggregation to produce a multicellular slug. This chemotaxis is dependent on actin remodelling, which is controlled by the Rac subfamily GTPases (Cdc42 and Rho are absent in this model system). However, experiments on latrunculin-treated Dictyostelium cells, where the actin cytoskeleton is poisoned and the cells immobilized, indicate that gradient sensing and polarization can be decoupled from the actin cy- toskeleton (88). As polarization occurs in response to external cues, cells possess a directional sensing system that acts as an internal ‘compass’ to guide them towards the chemoattractant. Dictyostelium is the model organism where the molecular details of how this directional sensing system operates are best understood (119). Other common model systems for cell polarity and migration include: 8 1.2. Understanding Cell Polarization in Theoretical Terms • neurons during neural system development, where growth cones at the tips of developing axons need to travel to synaptic targets and are guided by chemoat- tractants and repellents (58); • fibroblasts, cells that repair connective tissue, in which crosstalk between different Rho GTPases was first discovered (31; 82); • keratocytes, fast moving fish epithelial cells that maintain a consistent shape and smooth gliding motion as they migrate (115). In all these diverse organisms, there are conserved mechanisms that contribute to polarization. In all of them, one or more Rho GTPases play a central role in polarity establishment. Note however that the models for cell motility presented in this thesis are most applicable to neutrophils and neutrophil-like cells, although the same principles may apply in other organisms. 1.2 Understanding Cell Polarization in Theoretical Terms In establishing a front and back, many intracellular components transitions from homo- geneous state to inhomogeneous distributions. Emergence of polarization can thus be viewed as a symmetry breaking phenomenon, making it an attractive subject for math- ematical modelers. To begin to study the mechanisms that enable cells to transition from a homogeneous to a polar state, I first list the features of cell polarization that should be captured by a successful model. 1. When an unpolarized cell is subjected to a transient stimulus in the form of a shallow chemical gradient (from a pipette containing a chemoattractant, for example), it responds by reorganizing into two well-defined regions, front and back. In most cases, this polarized state is maintained even after the transient stimulus is removed (maintenance). 9 1.2. Understanding Cell Polarization in Theoretical Terms 2. Another hallmark of cell polarization is its sensitivity to new incoming signals: when a polarized cell is subject to a sufficiently strong secondary stimulus, it adjusts its polarization accordingly. For example, reversal of orientation is seen in neutrophils (white blood cells) when the stimulus gradient is reversed (2). 3. It has been shown that cells can respond to very shallow gradients (as small as 2-5% difference between front and back). Receptor occupancy and heterotrimeric G-protein activity that transmits the external signal to downstream pathways parallels the concentration of chemoattractant (41). Amplification of the initial signal occurs in the later stages of gradient sensing, resulting in sharp localized responses at the leading edge and rear of a polarized cell. 4. Cells must also possess a mechanism for adaptation in a uniform stimulus when receptor occupancy is held constant. Note that the adaptation mechanism in eukaryotic cells is distinct from that of bacteria, where adaptation is achieved through methylation and demethylation of receptors (20; 49; 73). 5. Polarization can also occur in a uniform pool of chemoattractant. Neutrophils and Dictyostelium cells lacking adenyl cyclase (ACA) break symmetry to form a dis- tinct leading and trailing edge and begin to migrate randomly. This phenomenon is called spontaneous polarization. The random location of this polarization sug- gests that it is stochastic in nature. 6. Migrating polarized cells must also remain sensitive to new stimuli, suggesting that the asymmetry created during polarization should be reversible. Note that some of these phenomena occur under different experimental conditions, and a given model may only describe some of these. For example, maintenance of polarity may be a distinct process from the initial symmetry breaking. Also, there might be a trade-off between robust polarization and sensitivity to new stimuli. It is typically assumed that, like other biological processes, the capacity of cells to polarize must remain robust to environmental and individual cell-dependent parameter changes. 10 1.2. Understanding Cell Polarization in Theoretical Terms How long polarity is maintained after establishment is dependent on both experiment and model organism in question. A cell migrating up a gradient will likely need to change direction, while choosing a location for a new bud site in yeast is probably irreversible. Gradient sensing in F-actin-inhibited Dictyostelium cells is a transient phenomenon. Polarization in keratocyte cell fragments, on the other hand, is sustained long after withdrawal of the stimulus (115). Additional feedback loops are likely responsible for long-term versus transient maintenance of polarity. Successful models for polarization should then be able to address the following ques- tions: 1. How are cells able to sense both steep and shallow gradients within a vast range of concentrations? 2. Why do cells spontaneously polarize in the absence of an external gradient? How can spontaneous polarization be reconciled with adaptation to a uniform stimulus that is also observed? 3. What contributes to the maintenance of cell polarity after it is established? 4. How does the polarized cell remain sensitive to new stimuli, yet robust to envi- ronmental changes? 5. How does the biochemistry of the polarization machinery allow for cell polariza- tion? Are there common conserved properties in such proteins, and, if so, why? 6. What are the minimal requirements needed to establish polarity? In addition to answering these questions, the ultimate goal of a successful model is to understand the mechanism responsible for cell polarization. This allows one to go beyond describing the observed behaviour, and make new predictions for cell behaviour under different experimental conditions. It provides insight into what are essential components, interactions, and dynamics. This is significant in allowing us to understand the universal features of many types of cells and common evolutionary mechanisms. It 11 1.3. Thesis Outline also highlights how that mechanism could be susceptible to mutations, pathogens, or chemical effects. The aim of this thesis is to understand how cells transition from a homogeneous to a polarized state. However, it is difficult to comprehend the large protein interaction (“signal transduction”) network involved in cell polarity, and even more challenging to understand how this entire network of proteins operate in the spatially inhomogeneous environment of the cell. Mathematical analysis of subnetworks (or “modules”) can point to generic phenomena associated with components that make up these subnetworks. This thesis is concerned with the mathematical analysis of one such module, comprised of the Rho GTPase signal transduction network. As I will show that module can undergo polarization by a mechanism we call “wave-pinning”. In my M.Sc. I developed a model of how the interactions of the GTPase family members lead to establishment of polarity (43). That model could reproduce many of the behaviours of polarized cells (44). This Ph.D. thesis is concerned with understand- ing the underlying polarity establishment mechanism in a simplified, mathematically tractable version of that model. This is done by reducing the biochemical details of the model described in (44), while retaining the fundamental properties of the Rho pro- teins necessary for polarity establishment. The purpose of this simplification is to allow a thorough mathematical analysis. The specific mathematical predictions about the mechanism generating the observed cellular responses are then tied back to the biology to generate new questions about cell polarity. 1.3 Thesis Outline The outline of this thesis is as follows: • In Chapter 2, I review the biochemical properties of Rho GTPases and propose a simplified model for cell polarization based on membrane-cytosol cycling of Rho GTPases. I subject this model to a number of tests based on experimental obser- vations and determine if the model can reproduce the behaviours listed above. I 12 1.3. Thesis Outline discuss what biological conclusions can be drawn from that model. • In Chapter 3, I review other mathematical models from the cell polarization lit- erature and compare them to the wave-pinning model in Chapter 2. I subdivide the models into several classes based on proposed mechanism for polarization and discuss what features of polarization can be explained by each class of models. • In Chapter 4, using singular perturbation theory, I explain the mathematical basis of the wave-pinning model from Chapter 2 in 1D. I compare the analytical results to numerical simulations of the PDE system. • In Chapter 5, I extend the asymptotic analysis from Chapter 4 to two spatial dimensions. I discuss the behaviour of the system on various timescales and relate our wave-pinning model to other well-known models for phase separation. • In Chapter 6, I analyze the transition from a spatially heterogeneous solution to a spatially homogeneous solution in the wave-pinning model using numerical continuation methods and phase plane arguments. I discuss my results in terms of bifurcation theory. 13 Chapter 2 Rho GTPases and the Formulation of the Wave-Pinning Model 2.1 Introduction Polarization requires formation and maintenance of different cellular structures at the front and back of the cell. Regulation of spatial localization of Rho GTPase activity is important for cell polarization, as Rho GTPases regulate the type of actin structure that will be built (branched actin network found at the front versus stress fibers that adhere the tail to the substrate, for example). There is strong evidence that different Rho proteins are spatially segregated in stimulated cells. In neutrophils and similar cell lines, active Cdc42 and Rac have been associated with a putative cell front (50; 78), and RhoA with the cell back (114; 121). (See Figure 1.3 for images of localization of Rho GTPases in various cell lines). Sharp spatial segregation between Cdc42 and Rho activity zones is also observed in the egg cells of the frog Xenopus (6). Note, however that in fibroblasts, cells that repair connective tissue, active RhoA is found at the leading edge of the cell, as well as in the tail (it is absent from the middle) (90), suggesting cell specific differences1. 1This difference may be related to the fact that neutrophils are fast moving cells that are loosely adhered to the substrate, while the slower moving fibroblasts have tight adhesion to the substrate. RhoA plays a key role in forming adhesive structures (focal adhesions) that are important for fibroblast locomotion and much less so for neutrophil migration (120) 14 2.1. Introduction In my M.Sc. work I focused on investigating interactions between the three Rho GTPases, Cdc42,Rac, and RhoA, which can give rise to their spatial segregation (43). Here I briefly review the main results of that work, as it directly motivated this Ph.D. research. It was shown that cooperativity in the GTPase interactions, as well as fast diffusion of the inactive forms of the GTPases was essential for obtaining robust spatial bistability (44; 64). The resulting model was a nonlinear system of 6 PDEs (an equation for the active and inactive form for each of Cdc42, Rac, and RhoA). In Section 2.2 we simplify the Rho interaction dynamics for the purpose of doing mathematical analysis on the system, and testing how much of the behaviour of that complicated model can be obtained with just a single GTPase. As this thesis focuses on a single Rho GTP-GDP pair as a basic polarizing unit, and looks at how its biochemistry relates to its role in polarity establishment, I begin by reviewing the biochemistry of Rho GTPases. 2.1.1 Guanosine Exchange Cycle of the Rho Proteins In the cell, Rho proteins undergo two coupled and carefully regulated cycles: cycling between GDP-bound and GTP-bound conformations, as well as cycling between cy- tosolic and membrane-bound forms (see Figure 2.1). These two cycles are regulated by three classes of proteins: GEFs (guanine nucleotide exchange factors), GDIs (GDP dis- sociation inhibitors) and GAPs (GTPase-activating proteins). Upon activation, GDP dissociates from the inactive form. GTP binding to the protein then changes conforma- tion of a region that binds a downstream effector (110). By itself, the dissociation of GDP from the GDP-bound form is an extremely slow reaction, so it is stimulated by various GEFs. These, in turn, are regulated by an upstream signal. In response to stimulation of cell surface receptors by chemical signals, GEFs accumulate at the leading edge of a motile cell (38). A GEF releases bound GDP from a Rho protein and forms a binary complex with the protein. The GEF in the complex is then replaced by GTP. In eukaryotic cells, GEF protein types outnumber Rho protein types by a factor of three, so multiple GEFs are capable of activating one 15 2.1. Introduction Rho protein (104). Why there exist so many apparently functionally redundant GEFs is not known, but this abundance might be the reason why cellular response to the same stimulus can vary from one cell type to the next. It is possible that, through subcellular localization or additional protein-protein interactions, GEFs can determine which downstream pathways are activated by the active protein (104). The second class of regulatory proteins, GAPs, stimulate the rate of GTP hydrolysis of the Rho proteins, causing their inactivation. There are about eighty different proteins, containing a 140 amino acid long RhoGAP-homology domain, encoded in the human genome, though it is not known if all of them function as Rho GAPs in vivo. As with GEFs, it is not yet clear why there are so many proteins that can potentially downregulate Rho activity. It is generally believed that the G-protein has to be both in the GTP-configuration and membrane-bound to interact with its effectors2 and trigger a cellular response, so the cytosolic form of the proteins is assumed to be inactive. When activated, the proteins become attached to membrane sites. However, little is known about the molecular mechanisms that allow the delivery and translocation of Rho proteins to specific sites in cell membranes. Rho GDIs regulate the association of the Rho family proteins with membranes, and have been shown to extract the GDP-bound form of Rho proteins from the membrane to the cytosol, where they are sequestered in an inactive pool (83). The membrane diffusion coefficients of Cdc42, Rac, and Rho are roughly equal, because the three proteins are of a similar size and structure. The diffusion coefficient is approximately 0.1µm2s−1 in the membrane , and 10µm2s−1 in the cytosol, so cytosolic forms diffuse approximately 100 times faster than the membrane-bound forms. (94). Experimental evidence suggests that Rho GTPase cycling described above is directly responsible for the polarization response. In yeast, for example, the accumulation of Cdc42 at sites of polarized growth depends crucially on a gene that codes for a single Cdc42 GEF (89). We and others (44; 64; 86) have hypothesized that exchange between 2Effector is a protein that directly binds to the GTPase when it is the GTP-bound state and mediates downstream signalling 16 2.1. Introduction membrane and cytosolic forms plays a practical role in generating robust cell polarity. In my M.Sc. thesis, I developed a model for the interactions of Cdc42, Rac, and Rho in a polarized cell that requires this membrane-cytosol exchange to achieve polarity. Plasma Membrane upstream signal Rho GDP Rho GDP Rho GTP effector Rho GDI GEF GAP IG δG kon koff Cytosol GTP GDP Pi cytoskeletal effects Figure 2.1: A schematic diagram of the activation and membrane translocation cycles of a typical Rho family GTPase. In accordance with current literature, only the inactive form, Rho-GDP, can come off and on the membrane. To move into the cytosol, Rho- GDP needs to bind to a GDI molecule. Nucleotide exchange can only take place on the membrane. Dissociation of GDP from the GDP-bound form is an extremely slow reaction, so for Rho GTPases to be effective molecular switches, GDP dissociation is stimulated by various GEFs. Although the GDP/GTP exchange reaction is reversible, in reality, it almost never occurs. Instead, GTP is hydrolyzed to GDP, a process that is stimulated by various GAPs. Only the active, GTP-bound, forms of the proteins can bind to their effectors to have downstream effects or crosstalk with other GTPases. 2.1.2 A Model for Spatial Segregation of Rho GTPases Based on Inhibitory Crosstalk I briefly review the Cdc42-Rac-Rho model (44) that motivates a simpler variant dis- cussed in this thesis. In the full model, we keep track of small G-proteins Cdc42, Rac, and Rho that are active (on the membrane) or inactive (both on membrane and in cytosol, bound to GDI). We chose this model to represent the mutual interactions as it is supported by hypotheses in the biological literature about mutual suppression of 17 2.1. Introduction Cdc42 and Rho, together with activation of Rac by Cdc42, and Rho by Rac (28). We do not explicitly account for the GEFs and GAPs that activate/inactivate these proteins. Assuming that all the crosstalk takes place through the GEFs, and that the inactive forms are in a dynamical equilibrium between membrane and substrate, we represent interactions of the small G-proteins by the following model, with C, R, and ρ defined as active GTP-bound Cdc42, Rac and Rho, respectively: ∂C ∂t = fC − δCC +Dm∆C , (2.1a) ∂R ∂t = fR − δRR+Dm∆R , (2.1b) ∂ρ ∂t = fρ − δρρ+Dm∆ρ, (2.1c) where the functions fi were nonlinear terms representing interactions shown in Fig- ure 2.2. We used fC = IC (Ci/Ctot) 1 + (ρ/a1) n , (2.1d) fR = (IR + αC) (Ri/Rtot) , (2.1e) fρ = (Iρ + βR) (ρi/ρtot) 1 + (C/a2) n . (2.1f) Here IC , IR, and Iρ are the basal activation rates in the absence of crosstalk, given that all small G-proteins are in the inactive state. a1 and a2 describe levels of the proteins that lead to a half-maximal mutual inhibition of Rho and Cdc42; and α and β are the activation of Rac by Cdc42, and the activation of Rho by Rac, respectively. δC , δR, and δρ are the basal inactivation rates (mediated by GAPs). Dm describes the diffusion rate of the small G-proteins along the membrane, here assumed to be the same for all three, based on comparable molecular weight and structure. Ci/Ctot is the local concentration of inactive Cdc42, relative to the effective level to be expected in a well-mixed system without any activated Cdc42 (similarly, for Ri/Rtot and ρi/ρtot). In our first analysis, 18 2.1. Introduction we took Gi = Ci, Ri, ρi to be constant parameters Cdc42 Rac Rho + + Figure 2.2: Schematic diagram corresponding to the scheme we propose for the Rho protein crosstalk in equations (2.1a)-(2.1c) based on (28). By simple bookkeeping, each inactive form of the small G-proteins, satisfies an equation of the form ∂Gi ∂t = δG − fG +Dmc∆Gi where G = C,R, ρ. (2.2) Because inactive Rho proteins spend some time on the membrane before GDI extraction Dmc is an “effective” diffusion coefficient that weighs membrane and cytosolic diffusion by relative residence times (44). The value of this diffusion coefficient should thus be in some range between the upper limit established by free diffusion of small G-proteins within the cytosol and the lower limit given by the membrane-bound diffusion, which is 10 to 100-fold lower. Moreover, all concentrations are defined as number of molecules “per unit area” to avoid the rescaling of terms due to unequal compartment volumes. All the variables described above, i.e. C, R, ρ, Ci, Ri, ρi, are expressed in units of [µM] and depict proteins confined to the cell. We assume Neumann boundary conditions along the edge of the cell (see Figure 2.4). In our model, we represented mutual inhibition between Rho and Cdc42 as a bista- bility: either Cdc42 is high (which leads to activation and elevation of Rac) and Rho levels are suppressed, or else, Rho is high and then Cdc42 and Rac are low. Whether bistability actually occurs in the model depends on the Hill coefficient, n, in fC , fρ. When n = 1, bistability is mathematically impossible for all parameter values. For the biologically realistic parameter values (see (44) for values), bistability in the temporal 19 2.1. Introduction back (x=0) Side View front (x=L) Nucleus membrane cytosol lamellipod (a) membrane cytosol Rho-GTP Rho-GDP GDI GAP GEF U V (b) Figure 2.3: Geometry of the model in 1D. (a) A cross-section of a cell (not to scale), in side view, showing membrane (shaded) and cytosol (white). x is position along the cell diameter (0 ≤ x ≤ L). Rho proteins are distributed in the cytosol (inactive, white circles) and across the membrane (active, filled circles). (b) A schematic diagram of the membrane translocation cycle of a typical Rho family GTPase and its approximation by the UV reaction-diffusion system. Straight arrows: exchange between active (u) and inactive (v) protein, curved arrow: positive feedback of u on its activation. Membrane thickness is exaggerated for visibility. Cytosolic (membrane-bound) Rho proteins are represented by open (filled) disks, and assumed to diffuse along the axis of the cell in both membrane and cytosol. The inactive cytosolic form diffuses much faster than the active membrane-bound form and is approximately uniformly distributed. 20 2.1. Introduction front (x=L) back (x=0) Figure 2.4: Schematic representation of a chemically polarized cell. In one spatial dimension we interpret the boundary conditions as the front and back edges of the cell. As proteins cannot leave the cell, it is appropriate to assume no-flux boundary conditions. dynamics only occurs for n ≥ 4. We found that a bistable regime of Rho protein dynamics, however, does not auto- matically lead to persistent polarity in a cell. In general, a stable transition zone does not occur at all in a spatially distributed bistable system with values at ±∞ fixed to the distinct steady states solutions. A travelling wave will, instead, sweep through the domain in a 1D setting (77; 81). Since one or the other equilibrium takes over, no stable coexistence is possible. If Rho GTPases are always on the membrane, and never re- leased into the cytosol, (i.e. if locally the amount of inactive GTPase could be described as Gi = Gtot − G), this type of transient behaviour occurs in our model. Fig. 2.5(a) shows that in the bistable regime, such a travelling wave sweeps through system (2.1), implying that polarity can only be maintained for at most around 5-10min, for realistic diffusion coefficients, and is then lost. These results illustrate the fact that the active forms of the proteins on their own, even with bistability in the ODE formulation do not produce the desired spatial polarization. This transient behaviour changes dramatically when the rapid diffusion of the GDI- bound small G-proteins in the cytosol (two-fold faster than in the membrane) is taken 21 2.1. Introduction space (µm) time (s) 0 298 0 10 3 17.35 Rac space (µm) time (s) 0 298 0 10 0 3.6 Rho space (µm) time (s) 0 298 0 10 0 3.1 Cdc42 (a) space (µm) time (s) 0 298 0 10 0.2 3.0 Rho space (µm) time (s) 0 298 0 10 1.0 9.0 Rac space (µm) time (s) 0 298 0 10 0 2.4 Cdc42 1 (b) Figure 2.5: (a) Space-time plot of system (2.1) on a 1D strip, 10µm wide, with Dm = 0.1µm2s−1. Subpanels show active Cdc42, Rac, and Rho with warmest colours for highest concentrations. (b) Space-time plot of system (2.3) on a 1D strip, 10µm wide, with Dm = 0.1µm2s−1, Dmc = 10µm2s−1. Subpanels show active Cdc42, Rac, and Rho with warmest colours for highest concentrations. All parameters are as given in (44), where they were based on available biochemical data (see also section 2.3.2). 22 2.1. Introduction into account, resulting in the following model: ∂C ∂t = fC(ρ) ( Ci Ctot ) − δCC +Dm∆C , (2.3a) ∂R ∂t = fR(C) ( Ri Rtot ) − δRR+Dm∆R , (2.3b) ∂ρ ∂t = fρ(C,R) ( ρi ρtot ) − δρρ+Dm∆ρ , (2.3c) ∂Ci ∂t = −fC(ρ) ( Ci Ctot ) + δCC +Dmc∆Ci , (2.3d) ∂Ri ∂t = −fR(C) ( Ri Rtot ) + δRR+Dmc∆Ri , (2.3e) ∂ρi ∂t = −fρ(C,R) ( ρi ρtot ) + δρρ+Dmc∆ρi , (2.3f) where the reaction kinetics fC , fR, fρ are as in (2.1). Here we have kept the same dynamics as in equations (2.1a-c), but we now consider inactive forms as dynamic variables with conservation ∫ L 0 G(x, t) +Gi(x, t) dx = constant, (2.4) for G = C,R, ρ and Gi = Ci, Ri, ρi . The rapid cytosolic diffusion means that differences in cytosolic concentrations are very rapidly equilibrated. The exchange between the membrane and the cytosol is also rapid (102). This leads to the results shown in Fig. 2.5(b): with appropriate initial conditions a wave of activation is initiated, moves into the domain, and stops inside. This leads to robust polarization, i.e. a difference in the activation state between the two ends of the cell. A transient gradient in the activation rate IG for one of the Rho proteins, which can be interpreted as a bias in the GEF distribution or activity, were used as “stimuli” to trigger the travelling wave. The activation of each of the small G-proteins, and thus its growth is limited by the availability of the inactive form. Thus, an equilibrium is rapidly established, where 23 2.2. Model Reduction to a Single Rho GTPase expansion of each of the small G-proteins is limited by the availability of the inactive form. The location of the transition front between the high and low steady states is determined by the amount of the inactive form available. Flooding the cell with excess Ci or ρi causes the transition zone to shift, causing the high Cdc42 or high Rho equilib- rium to take over the entire cell in some cases. The fast cytosolic diffusion of Ci, Ri, ρi is essential for stabilizing the travelling wave and establishing robust polarization. The model comprised of equations (2.3) consists of six PDEs, and is challenging to analyze. This motivates my work in this thesis, where I study a simpler variant. In this thesis I will explain how this wave-pinning mechanism works at the level of an individual Rho-protein. 2.2 Model Reduction to a Single Rho GTPase 2.2.1 Modelling Assumptions and Simplifications Based on the previous sections we make the following biologically reasonable assump- tions: 1. The Rho proteins are not synthesized or degraded on the timescale of cell polar- ization (30 seconds- several minutes) but only undergo conversion between active and inactive forms. Therefore, the total amount of protein remains constant over the time scale on which the cell polarizes. 2. The membrane diffusion coefficients of Cdc42, Rac, and Rho are roughly equal, because the three proteins are of a similar size and structure (123). 3. Cytosolic forms diffuse approximately 100 times faster than the membrane-bound forms (94). As a result, the distribution of the cytosolic forms will be approxi- mately homogeneous in space, even when the active forms are polarized on the membrane. 4. The amount of GDI present in the cell is sufficiently large that it is not a limiting 24 2.2. Model Reduction to a Single Rho GTPase factor in extracting GDP-bound Rho protein from the membrane. 5. Each Rho protein has a basal rate of activation (GTP-loading), and a basal turnover or inactivation rate (GTP-hydrolysis). Activation and inactivation of Rho proteins are indirect processes via the GEFs and GAPs, respectively (96). 6. With the appropriate crosstalk, the system of GTPases polarizes, so that, for ex- ample, the active form of Cdc42 is high at the front, and low at the back of a motile cell. We assume that the reaction kinetics for a single GTPase should exhibit bistability once appropriate crosstalk terms are included in the GEF-mediated activation rates. 2.2.2 The Well-Mixed System Consider, first, a well-mixed system consisting of an active form, U , and an inactive form, V , of the signaling protein with interconversion: V kvu À kuv U, (2.5) (See Figure 2.3.) The concentrations of these proteins u, v would satisfy kinetics of the form du dt = Rate of activation− Rate of inactivation, = kvu v − kuv u, ≡ f(u, v) (2.6) where kvu, kuv are rate constants (units of 1/t), possibly dependent on u, and where we have defined the function f(u, v) above. If only this exchange takes place, and no new material is produced or destroyed, then conservation in the well-mixed system implies that dv dt = −du dt . 25 2.2. Model Reduction to a Single Rho GTPase We will refer to kvu v as the rate of activation, kuv u as the rate of inactivation, and f(u, v) as the reaction kinetics. Bistability necessitates that either one or both of kvu, kuv be nonlinear functions of u or v. As v is inactive we assume that one of those quantities depend on u. Here we have chosen one of the simplest assumptions, that there is cooperative positive feedback from the activated form onto its own production (via GEFs), whereas the reverse conversion (mediated by GAPs) takes place at constant basal rate, δ. Expressions for the rates consistent with these assumptions are kvu = k0 + γu2 K2 + u2 , kuv = δ, where k0 is a basal rate, and positive feedback is represented by a Hill function with maximal rate γ and saturation parameter K. The Hill coefficient n = 2 is needed to achieve bistability. The resulting kinetic function f(u, v) = v ( k0 + γu2 K2 + u2 ) − δu, (2.7) has the essential properties we need, and similar kinetics have been used in (71; 87), and in a number of superficially similar models based on Turing pattern formation. (A detailed comparison of these distinct mechanisms will be made in Chapter 3). However, to a large extent, the details of this choice for the kinetic function f(u, v) are arbitrary; only the qualitative aspects of the reaction term, and in particular, its bistability are critical. The well-mixed system has steady states when f(u, v) = 0, i.e., when v ( k0 + γu2 K2 + u2 ) = δu. (2.8) For a fixed value of v, the corresponding values of u at these steady states can be visualized as intersections of a sigmoidal curve (LHS of Eq. 2.8) with a straight line of slope δ (see Fig. 2.6). Up to three such intersections can occur, u−(v), uT(v), and u+(v), the outer two of which are stable. We label the steady states in increasing order, so that u− < uT < u+. The existence of at least two stable values of the activated form of the 26 2.2. Model Reduction to a Single Rho GTPase stable stable saddle u reaction rates uT u+(v) u-(v) Figure 2.6: The reaction terms: LHS (sigmoidal curve) and RHS (straight line) of Equation (2.8) as functions of u for a fixed v. Up to three steady states can occur, u−(v), uT(v), and u+(v), two of which are stable. For the wave-pinning mechanism to work, we will require that there is some range of values of v (vmin < v < vmax) for which the above three steady states exist. protein, u, is strongly suggested by experimental observations: these would correspond to the values of u at the front and at the back of polarized cells. More generally, any kinetic function f(u, v) that admits such steady states in a range of values of v would be suitable for our mechanism to work. Later, we simplify further to cubic kinetics for ease of analysis. 2.2.3 The Spatially Distributed System We now consider a spatial version of a reaction-diffusion system based on the bistable kinetics described above. We first consider a 1D spatial domain, Ω = {x : 0 ≤ x ≤ L}, that forms the axis of a cell. For simplicity, we take this domain to be a strip of small and constant width (w) and thickness (neglecting bulging nucleus of the cell, and neglecting any gradients not in the x direction). We denote by u(x, t) and v(x, t) the concentrations of the active (bound to the cell membrane) and inactive forms (GDP- bound in complex with GDI’s in the cytosol) of the protein, respectively, at position x and time t along this axis. These concentrations are defined as number of molecules per unit “area” to avoid the rescaling of terms due to unequal compartment volumes 27 2.3. Results (for details on appropriate handling of membrane-cytosol conversion please see (44)). In Fig. 2.3, we schematically represent this 1D geometry, including both membrane and cytosolic compartments and show the interconversion between active and inactive Rho GTPases. The concentrations of both u and v that are typical of a resting unpolarized cell and a cell after polarization are shown in Fig. 2.7. Then ∂u ∂t = Du ∂2u ∂x2 + f(u, v), (2.9a) ∂v ∂t = Dv ∂2v ∂x2 − f(u, v). (2.9b) The rate of diffusion of the membrane-bound (active) form is known to be signif- icantly smaller than its inactive cytosolic counterpart (94), so that Du ¿ Dv. We assume no flux at the ends of the domain, so that ∂u ∂x ∣∣∣∣ x=0,L = 0, ∂v ∂x ∣∣∣∣ x=0,L = 0. (2.9c) In our case, a polar pattern would be a chemical distribution that is highest at one end of the cell, and decreases to its lowest value at the opposite end of the cell. Equation (2.9c) implies mass conservation in the spatially distributed system (2.9), i.e. ∫ L 0 (u+ v) dx = Ktotal. (2.10) 2.3 Results We next analyze model (2.9) to see under which conditions the system can exhibit po- larity. To avoid confusion about terminology, we define the following. Bistability means existence of two stable equilibria in the ODE model. Spatial bistability means that the PDE model exhibits zones of stable size, throughout which more or less constant con- centration values are maintained. As was shown for systems (2.1) and (2.3), bistability 28 2.3. Results membrane membrane cytosol cytosol Resting Cell Polarized Cell u v 0 Lx u v 0 Lx Figure 2.7: A schematic diagram of stable profiles of u and v across the cell. Left: spatially homogeneous in resting cells. Right: u is sharply graded, but v is nearly uniform in polarized cells. of the ODE model does not imply spatial bistability in the PDE formulation. Polarity means spatial bistability with exactly two zones, a front and a back. Figure 2.7(b) represents a typical resting (spatially homogeneous solutions) versus polarized state (spatially nonhomogeneous solution). Equation (2.9a) considered by itself (with v as a parameter), supports travelling wave behaviour (depending on parameters), where one of the two stable homogeneous equilibria takes over the entire domain, depending on the initial conditions. This is the same behaviour that was exhibited by the 3 PDE model (2.1) for the active forms of the GTPases only. The full model for GTPase interactions, system (2.3), has a travelling wave solution that quickly pins inside the domain, resulting in a polarized distribution of the active forms. We look for similar behaviour in the reduced model for a single Rho GTPase (2.9) for a variety of input stimuli. 29 2.3. Results 2.3.1 Input Stimuli An external stimulus is known to increase activation of proteins such as Rho GTPases by upregulating the GEFs that convert inactive forms to active forms. Depending on the nature of the stimulus, this increased rate of activation could be space and time dependent. To model the external stimulus, we superimpose an additional transient spatial stimulus-dependent activation function fS . fS(v) = kS(x, t)v, (2.11) where kS(x, t), the increased rate of conversion of V to U due to an external signal has a spatiotemporal dependence. The specific forms of kS we used in simulations are listed in the Appendix. Note that this leads to conversion of v to u resulting in ∂u ∂t = Du ∂2u ∂x2 + f(u, v) + kS(x, t)v, (2.12a) ∂v ∂t = Dv ∂2v ∂x2 − f(u, v)− kS(x, t)v. (2.12b) 2.3.2 Parameter Values For computations, we set parameter values as follows: We take cell diameter to be L = 10µm, and diffusion coefficients in the membrane (Du = 0.1µm2s−1) and in cytosol (Dv = 10µm2s−1) as in (94). Similar values have been measured for membrane-bound Cdc42 in yeast (62). The average membrane lifetime of an activated Rac molecule is 2 s (102), while GAP-stimulated GTP hydrolysis of Rho has been measured as 1.5 s−1 (124). We thus let the inactivation rate be δ = 1 s−1. We assume that GAP hydrolysis must be matched in magnitude by GEF activation, taking γ = 1 s−1. For the baseline GEF activity we take k0 = 0.067 s−1. The parameter K has units of concentration of U . We normalize concentrations so that K = 1 (units of length−1). 30 2.3. Results 2.3.3 Numerical Simulations All numerical simulations were performed using a second-order implicit-explicit finite difference method (101). Crank-Nicholson method was used for the diffusion operator, and Adams-Bashforth method for the reaction. We first demonstrate that the above model reproduces two important features of cell polarization, amplification and maintenance that were discussed in Chapter 1. In the numerical experiments described below, we start with the same total amount of material in the domain. We start with an unpolarized cell in which u is uniformly equal to its lower steady state value u−, corresponding to the spatially uniform value of v. We subject this cell to two types of transient stimuli, localized (klocS ) and graded (k grad S ), as described above and in the Appendix. From the computational results (Figure 2.8(a,b)), we see that both external stimuli klocS and k grad S quickly result in a local increase in the active form u near the cell edge at x = 0. This local rise is spatially amplified by a sharp front of the active form that propagates into the cell. This propagation slows down and eventually halts, a phenomenon we call wave-pinning. The pinned front represents a clear segregation of the cell into front and back. This polarized state is maintained indefinitely, even after kS = 0. The final steady state solution is identical for both klocS and k grad S . We also test whether our model cell exhibits sensitivity to new incoming signals, and whether it can reverse its polarity. In a 1D simulation, we study this as a response to a sequence of two stimuli of the form kgradS (x, t) one following another, but with reversed polarity. In Figure 2.8e): we first give the cell an initial graded stimulus, and then we reverse the gradient of the stimulus. The model cell at first polarizes in the direction of the initial stimulus. After the second stimulus is presented, the front propagates out of the domain, and the cell repolarizes in the opposite direction. We model spontaneous cell polarization by starting with random initial data and no signal u(x, 0) = u−(v0) · 2R, v(x, 0) = v0, kS ≡ 0, 31 2.3. Results 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=20 t=200 u (a) 0 2 4 6 8 10 0 0.5 1 1.5 x t=0 t=20 t=200 (b) 0 2 4 6 8 10 0 0.5 1 1.5 x t=0 t=20 t=200 (c) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=20 t=200 u (d) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=200 t=400 u (e) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=200 t=0 u (f) Figure 2.8: Polarization behaviour in response to six distinct stimuli for the wave- pinning system 2.12 with reaction term (4.4). Detailed forms of stimuli are given in Appendix A.1. Initial conditions (t = 0s), transient (t = 20s), and asymptotic behaviour (t = 200s) are shown. The final position of the wavefront does not depend on the stimuli (dotted lines), but on the total amount of material, which is the same in all cases. (a) A localized stimulus klocS , of amplitude of 0.05 (dotted line, enlarged for clarity) is applied at t = 0 for 20s. (b) Two localized stimuli of same amplitude are applied at t = 0 for 25s. A transient double peak distribution follows. After the stimuli are removed, the system reorients back to a polarized solution. (c) With a smaller diffusion coefficient Du = 0.01µm2s−1 the two peak state becomes metastable. (d) A graded stimulus k grad S of amplitude 0.01 is applied at t = 0 for 20s. (e) A graded stimulus of amplitude 0.02 applied at t = 0 for 20s is followed, at t = 200s, by a similar graded stimulus of reverse orientation for 20s. The first stimulus defines one polarization , but the second stimulus causes the system to reverse polarity. (f) Noise of sufficiently large amplitude, superimposed on the homogeneous initial conditions also leads to polarization. 32 2.3. Results sp ac e time Du=0.1 0 50 100 150 200 0 1 2 3 4 5 6 7 8 9 10 (a) sp ac e time Du=0.01 0 50 100 150 200 0 1 2 3 4 5 6 7 8 9 10 (b) Figure 2.9: Space-time plots showing the response of the wave-pinning model (2.9) to two pipette experiments. Two localized stimuli klocS , of amplitude of 0.05 are applied at t = 0 for 25s and then turned off. (a) For parameter values given in the text, a transient two peak state emerges, with one of the peaks later collapsing. (b) If diffusion coefficient is decreased to Du = 0.01µm2s−1 the two peak stage persists. where 0 < R < 1 is a uniformly distributed random variable. Depending on the initial data, the cell can polarize in either direction (see Figure 2.8(f)). In some runs we also observed solutions with persistent multiple peaks (not shown). Such multiple peaks have been observed experimentally in “two pipette” experiments if the cell simultaneously receives two localized signals at its opposite ends (40; 42). When our model cell is exposed to two localized stimuli, waves enter the domain from opposing directions and appear to set up two opposing “fronts” at the two ends and a “rear” in the middle of the cell (Figure 2.8(b) and Figure 2.9a). After the stimuli are removed, however, the system reorients back to a polarized state with a unique front. Note, however, that this phenomena depends on the size of the domain. For larger domains (or smaller diffusion coefficient Du) the two peak state is retained (Figure 2.8(c) and Figure 2.9b). A peak in the middle of the domain is also possible (not shown). As I later discuss in Chapter 5 such multiple peak solutions are metastable and eventually coalesce to form a right or left polarized state. In simulations with various stimuli of larger slope or magnitude, the response of the cell remains the same (not shown). After some transients, the level of the active form of 33 2.4. Summary protein, and the location of the sharp front dividing front from back is the same. Note also that the position of the pinned wave is stable to small perturbations in the shape of the solution, provided the total amount of the material in the cell is kept constant (not shown). This feature of the model will be explained by our analysis and provides one prediction that distinguishes this mechanism from alternative mechanisms for polarity, which will be discussed in Chapter 3. 2.4 Summary In this chapter I presented the wave-pinning model as a reduction of a six PDE system describing the interactions of Rho GTPases in the cell to a simplified two PDE system representing a single active-inactive Rho GTPase pair. I subjected the system to several types of stimuli (representing a cell exposed to one or multiple pipettes, a cell in a chemoattractant gradient that may disappear or change directions, and a cell placed in a uniform field of chemoattractant) and described the system’s response. These results were all obtained through numerical simulation. I found that the model was sensitive to a bias in these stimuli and polarized accordingly. Polarity was maintained after the stimulus was removed. In Chapter 3, I compare the results obtained here to the predictions of other mod- els for cell polarity based on other proposed mechanisms. In Chapter 4, I will present mathematical analysis of this model that will rigorously explain why wave-pinning oc- curs and why, for a given amount of material, the wave stops at the same location, regardless of the magnitude or type of stimulus. 34 Chapter 3 Comparison of Cell Polarization Models 3.1 Introduction To understand complex emergent phenomena such as cell polarization, an approach taken by some mathematical biologists is to consider the underlying molecular network that gives rise to the phenomena, much like an engineer views a wiring diagram of an electrical circuit (54). In this reductionist view, complex signalling pathways in a molecular network can be viewed as a series of interconnected modules comprised of a few components (proteins and other cell metabolites) that are linked together by chemical reactions. Each module is designed to give rise to a particular behaviour, such as bistability, oscillations, and adaptation to new stimuli (112). When these modules are further coupled by positive and negative feedback loops, more complex behaviour can emerge. By understanding the behaviour of these simple modules in isolation, we can hope to understand the function of the molecular network as a whole. In this chapter, I describe the the various mechanisms used to explain the different features of cell polarity. I survey the models for cell polarization found in the litera- ture, and categorize them according to the types of mathematical behaviour each model exhibits. Finally, I compare these different prototype models to the wave-pinning mech- anism derived in Chapter 2. By subjecting all the proposed mechanisms to the same set of stimuli I can then objectively compare the different models to see if they satisfy the features of cell polarization listed in the introductory chapter. 35 3.2. A Survey of Models 3.2 A Survey of Models 3.2.1 Model Geometry Unless otherwise noted the models of polarization described here are typically reaction diffusion systems considered in 1D. Two common approaches are used. The first, shown in Figure 3.1 (a), is to take a thin slice across the cell. In that case, the domain for the spatial problem is a line, 0 ≤ x ≤ L, where L is a cell diameter. The boundary conditions are at the front and back of the cell. Since the proteins cannot diffuse out of the cell, it is appropriate to have no flux boundary conditions at x = 0, L. In that case, a polar pattern would be a chemical distribution that is highest at one end of the domain, and decreases to its lowest value at the opposite end of the domain3. An alternative, shown in Figure 3.1 (b), is to treat the cell as a disk-like object, and the cell membrane as the perimeter of that disk4. In that case, the domain is an interval in 1 D (0 ≤ x ≤ piL) with periodic boundary conditions. In this case, a polar pattern would be a chemical distribution with a single peak. 3.2.2 Turing Type Models Emergence of polarization can be viewed as a symmetry breaking event: the cell tran- sitions from a spatially homogeneous state into an inhomogeneous one. It is commonly assumed that the unpolarized state is stable to uniform perturbations, but not to spa- tially inhomogeneous perturbations that arise as a result of external chemical gradients or intrinsic random fluctuations (i.e., spontaneous polarization). This supposition is often taken to mean that polarization is due to a diffusion-driven (Turing) instability, that was first described by Alan Turing in 1952 (111). In Turing pattern formation the spatially homogeneous steady state is destabilized by the presence of spatially inho- mogeneous perturbations, and a pattern consisting of stripes or spots with an intrinsic 3For Turing models this requires that only the eigenvalue corresponding to the lowest wave number should be unstable. Excitation of higher modes will lead to a multi-peak pattern, which is undesirable. 4In some of those cases, the interior of the cell is considered to be spatially uniform, and the fast-diffusing component is represented by its spatial average. 36 3.2. A Survey of Models front (x=L) back (x=0) (a) (b) 0 x L 0 x L Figure 3.1: Two ways of interpreting cell polarity models in one spatial dimension: (a) as a slice from the top of the cell and (b) as the boundary of a disk. length scale then emerges (77). The basic idea behind a Turing mechanism rests in having two processes with differ- ent spatial characteristics: one local, promoting activation, the other global, suppressing it. Small fluctuations in protein concentration are amplified by the slowly-diffusing lo- cal activator, which is assumed to be autocatalytic. The global, inhibitory process is necessary to confine the growth of the activator to a localized region and prevent it from taking over the entire domain. The kinetics for the inhibitor must be slow com- pared to the activation, and diffusion must be fast relative to the activator. This idea also appears as so called “lateral inhibition” in other circumstances, e.g., neuronal net- works. In the context of cell biology, the activator is typically assumed to be a large, membrane-bound component, while the inhibitor is typically a small, cytosolic molecule. Alternatively, in such Turing based models the inhibitor can be replaced by a substrate, which is consumed during activator synthesis and stimulates the growth of the activator. Cytoplasmic diffusion is sometimes replaced by averaging in the model. 37 3.2. A Survey of Models In these models the production of the activator is typically assumed to be an auto- catalytic saturating function. Saturation is necessary in order to attain a more realistic stripe pattern (plateau pattern in 1D), rather than high spikes of activator activity seen in the classic activator-inhibitor model. Necessary conditions for diffusion-driven instability in a two component system are found in most classic mathematical biology text books (see, e.g. (77) and Appendix B.1 for detailed derivation). In a two component system, only two types of systems (shown in Fig. 3.2a) can form Turing patterns, given appropriate parameter values: 1. A slow diffusing activator (u) promotes its own activation, and also activates a fast-diffusing component (v), which inhibits both of the reactants. This type of system is called activator-inhibitor system. The patterns generated by these types of systems have the activator and inhibitor in phase, i.e., peaks of one coincide with peaks of the other. 2. A slow diffusing component u promotes its own activation, but during production depletes its fast diffusing substrate v. The substrate, v, promotes the activation of u. This type of behaviour is called a substrate-depletion system. The patterns generated by these types of systems have activator and substrate out of phase. For systems of three or more components, a homogeneous steady state cannot be- come unstable unless the reaction-diffusion system contains a Turing unstable subsystem (103). A criterion for a stable steady state to be destabilized is given in (116). The transition to a polarized state is irreversible in a Turing system, as the homo- geneous steady state is unstable. Regardless of the strength of external stimulus, the system will not evolve back to an unpolarized state. The intrinsic length scale of the Turing pattern can be determined from the dispersion relation, which depends on the reaction kinetics, diffusion coefficients of the components and the size of the domain (77). As discussed earlier, a polar distribution is a pattern where only the lowest wave number is unstable to perturbations. This places a constraint on the parameter values, 38 3.2. A Survey of Models as multiple peak solution can arise due to higher order wave-numbers becoming unsta- ble due to diffusion-driven instability. The timescale on which the polar pattern forms will be discussed in Section 3.3.2. Next, some examples of Turing models with applications to cell polarization are presented. Activator Inhibitor Activator Substrate (a) S F B I (b) Figure 3.2: Examples of Turing mechanisms. (a) A two-component Turing mechanism involves a short range autocatalytic activator and a long range inhibitor. Alternatively, a limiting substrate can play the role of a global inhibitor. (b) A three-component Turing mechanism based on mutual inhibition from (79). An external signal S activates two activators (F and B), which mutually inhibit each other. They, in turn, activate a global inhibitor I. Turing Models for Cell Polarization The first application of a Turing mechanism to biological pattern formation in tissue development were described in the 1970s (27; 72), and analyzed thoroughly for pattern- forming properties (see e.g. (77)). The applicability of a Turing mechanism to cell polarization was championed by Meinhardt (69; 72). Many variants of Turing systems are present in the current cell polarization modelling literature (29; 79; 86; 108) and are discussed below. These models achieve a high degree of signal gradient amplification due to the presence of a positive feedback loop, though parameter space where Turing instability is possible is typically narrow (77). 39 3.2. A Survey of Models The first application of a Turing system to cell polarization was based on the classic Gierer-Meinhardt local activator-global inhibitor model (72). As noted before, in a Turing system, the transition to a polarized state is irreversible. Once the cell polarizes, it will remain polarized. This may be desired under certain conditions (an example is the formation of the Cdc42-cap in budding yeast to mark where the new daughter cell will form), and undesirable in others (a migrating cell such as Dictyostelium may need to rapidly change direction in response to changes in external signal). To alleviate this problem, a second local inhibitor with a longer half life than the activator was added to the activator-inhibitor system (69). The second activator destabilizes the established spatial pattern (eliminating “self-locking”) and enables a second pattern in time to emerge. The biochemical identities of the molecules in this mechanism are still speculative, despite a decade of experimental investigations. A more biochemically realistic Turing model based on the phosphoinositide metabolic cycle was developed by Narang and co-workers (80; 108), where the local activator is given by the membrane phophoinositides, while the global inhibitor is considered to be the cytosolic pool of inositides. To attain more realistic behaviour in these models, an adaptation module (see next section for an explanation) was put upstream of the activator-inhibitor model. Activator-inhibitor and activator-substrate models require the activator to be auto- catalytic, i.e. there is a positive feedback loop promoting activator production. Exper- imental data in neutrophils suggests that the front and back localizing proteins inhibit each other (121). Since double negative feedback is mathematically equivalent to a positive feedback loop (21), it is possible to replace the autocatalytic activator by two mutually inhibitory activators. Narang showed that a Turing instability is impossible in a two component mutually inhibiting system, and derived a minimal three compo- nent Turing model (79). This model consists of two mutually inhibiting local activators that promote the synthesis of a diffusible inhibitor, which, in turn, inhibits both the activators (see Figure 3.2b). 40 3.2. A Survey of Models It is possible for the role of the activator and substrate in a Turing mechanism to be played by the same molecular species with two distinct states. Goryachev et al. (29) presented a model for spontaneous yeast polarization (in the absence of landmark proteins and actin polymerization) that is based on the cycling of Cdc42 between GTP- bound and GDP-bound states. This detailed eight variable model for various forms of Cdc42, its GEF Cdc24, and the Cdc42 effector Bem1 can be reduced to a prototypical two-component Turing activator-substrate system (29). Both phosphoinositide models by Narang (80; 108) and the model by Goryachev (29) have an additional feature in common: the total amount of phosphoinositides in the plasma membrane and the total amount of Cdc42 in membrane and cytosol is conserved. Otsuji et al. (86) consider the effects of mass conservation on several Turing systems. Typically, in diffusion-driven instability, multiple peaks can arise from a homogeneous steady state. As stated earlier, multiple peak solutions are considered undesirable in models of cell polarity. Often, unrealistically high values of diffusion coefficient for the inhibitor are required to prevent multiple peak solutions. Otsuji et al. (86) showed that in their particular mass-conserved system multiple peaks can initially arise, but the depletion of molecules as these peaks grow causes the multi-peak state to become unstable and the smaller peaks collapse. As a result, though transient multi-peak states exist, the stationary state consists of a single stable peak. Although (86) showed that multi-peak states are unstable for a specific Turing model, this “winner-take-all” be- haviour is also seen in models by Narang (80; 108), Goryachev (29), suggesting that mass-conservation may play an important role in polarity, regardless of the detailed molecular network considered by the model. 3.2.3 Adaptation Models Next, we focus on the transient behaviour of cell in response to chemoattractant. To explain the ability of cells to generate a persistent, amplified response to a gradient of chemoattractant, but transiently respond to a uniform stimulus, a series of so-called 41 3.2. A Survey of Models “gradient-sensing” models have been developed. The goal of these models is to present perfect adaptation to spatially uniform signal, meaning that while there may be a tran- sient response to changes in signal strength, the steady-state of the response element will be independent of the signal. These models are motivated by experiments done in latrunculin-treated Dictyostelium cells, earlier described in Chapter 1. Perfect adaptation networks in spatially homogeneous (ODE) settings have been well characterized (see, for example, the “sniffer” mechanism described in (112)). Numerous examples of robust ODE models for perfect adaptation can be found in the literature on chemotaxis of bacteria (5; 73). Recently, a computational study Ma et al. (61) of possible network topologies of three components revealed that perfect adaptation can only be achieved with two types of networks: negative feedback loops and incoherent feed-forward loops5 (see Figure 3.3 for an illustration). (a) Signal A B C Signal A B C (i) (ii) (b) Figure 3.3: (a) Perfect adaptation is a feature of networks to reset their outputs to original levels after a transient response to stimulus. (b) For three component networks, only two types of networks can give rise to adaptation: (i) incoherent feed-forward loops, and (ii) negative feedback loops. Here A and B are intermediates and C is the response element. 5A feed forward loop is a network in which a single input modulates the response via two or more intermediate pathways. In an incoherent feed-forward loop, the intermediate pathways have opposite roles (i.e., activation and inhibition). 42 3.2. A Survey of Models Spatial Models with Perfect Adaptation To explain how cells adapt to uniform increases in chemoattractant levels and undergo persistent signalling in response to a gradient, Levchenko and Iglesias proposed a sim- ple model based on the assumption that the signal activates both a fast acting local activator, A, and a slow acting global inhibitor, I, in direct proportion to an external stimulus S that depends on space and time (55). A third, localized chemical R, creates the response. At steady state the response is determined by the ratio of activator to inhibitor. The idea of the model is that in the case of spatially uniform stimulus S, that ratio does not depend on not extracellular signal intensity, giving rise to perfect adaptation. In the case of a spatial gradient, the assumption that the activator is produced locally, while the inhibitor is globally diffusing leads to more activator being present at the front than the back, and a graded distribution of the response element, R, follows. The LEGI (local excitation, global inhibition) mechanism (see Figure 3.4) is a spatial version of the incoherent feed-forward loop mechanism shown in Figure 3.3a. The LEGI model by itself does not amplify the external gradient, requiring additional mechanisms further downstream to achieve needed amplification (51; 55). Iglesias and coworkers have championed the LEGI model to explain the dynamics of the phospholipid PIP3 (one of the first signaling components to sharply localize in a chemoattractant-stimulated cell) in latrunculin-treated Dictyostelium cells (18; 34; 52; 55; 60). The LEGI mechanism is simple and robust with regards to changes in parameters, but it has been difficult to associate the variables with actual identifiable molecular components. Initially a kinase that synthesizes PIP3, PI3K, was proposed to be the activator, and PIP3’s dephosphotase PTEN was proposed to be the inhibitor (55). However, upon stimulation, PTEN localizes to the back of the cell and PI3K to the front, while in the model both the activator and the inhibitor are activated in the same region. The activator and inhibitor were then postulated to be upstream of PI3K and PTEN. A model with two parallel LEGI mechanisms has been proposed to regulate 43 3.2. A Survey of Models S A I R (a) S A Bm B (b) Figure 3.4: Examples of adaptation models. (a) Local excitation, global inhibition (LEGI) mechanism for gradient sensing. An external signal (S) is assumed to activate two pathways, one involving a local activator, (A) and the other a global inhibitor (I) of the downstream response element (R). (b) Balanced inactivation mechanism for gradient sensing. An external signal activates two pathways, one involving a local membrane-bound activator (A) and the other a cytosolic inhibitor (B). B is converted to a local inhibitor (Bm), which is used up inhibiting A. the binding sites for both PI3K (found at the front of the cell), and PTEN (found at the sides and back of the cell) (60). Despite intensive experimental efforts, the identity of the activator and inhibitor is currently lacking. The LEGI model shows an increasing response to stronger gradients, with steeper gradients exhibiting higher responses. This is seen experimentally in latrunculin-treated Dictyostelium cells. The model also responds appropriately to a reversal in gradient location. However, the LEGI mechanism requires an external gradient in order to remain polarized, returning to baseline response level once the signal is removed. This means that LEGI modules would not be suitable for maintaining robust polarization in cells (such as neutrophils) that retain polarization even after a stimulus is removed. Although the LEGI mechanism accounts very well for the behaviour of unpolarized Dictyostelium cells in the absence of actin cytoskeleton, it is clear that additional mechanisms are needed to explain polarity (39). Another adaptation model that shares some common features with the LEGI model 44 3.2. A Survey of Models is the “balanced inactivation” model (56), that is shown in Figure 3.4. In that model, the external signal leads to the production of a local, membrane-bound activator A and a global, cytosolic inhibitor B at equal rates. Similar to the LEGI system, the fast diffusion of B leads to a higher ratio of A/B at the front than at the back of the cell. In the LEGI model, the response of the system is proportional to this ratio. In the “balanced inactivation” model, B is converted to a membrane-bound inhibitor, Bm that inactivates A, and is used up in the process. The addition of this extra component allows a switch-like response to external gradients, leading to a well-defined front and back. Unlike Turing models, both these gradient-sensing models contain too few nonlin- earities and spatial terms to have inherent pattern-forming features of their own, and are not intended as such. Instead, LEGI and related models are ideal as a readout of persistent macroscopic signals of various amplitudes. They are able to rapidly reverse subcellular organization in response to gradient reversal. They also respond appro- priately to multiple simultaneous chemoattractant sources. Successfully coupling such gradient-sensing mechanisms to modules for polarization is currently the goal of many mathematical modellers in the field (33; 85). 3.2.4 Wave-Based Models While LEGI-like models are monostable, with a unique response to a given stimulus, there is another class of models that have bistable behaviour, informally known as a “toggle switches”. Bistable behaviour is characterized by the presence of two possible steady states for a given input, and is often driven by positive feedback. The particular state at which a bistable system is at a given time is determined by the system’s history. We discuss this property, called hysteresis, below. 45 3.2. A Survey of Models Bistability in an ODE model For bistability to be possible, a system must exhibit some sort of feedback, such as positive feedback or double negative feedback (21). However, that is not a sufficient condition. The set of differential equations governing the system must be nonlinear, and the conditions necessary for the existence of a bistable “switch” are determined by the degree of nonlinearity. A common nonlinear feature of many models with multiple stable steady states is a sigmoidal shaped profile of activation or inactivation kinetics versus concentration. This sigmoidal shape often occurs due to cooperativity. A simple well-known equation that exhibits bistable behaviour is: dA dt = s+ kon An Kn +An − koffA. (3.1) Here A is the level of an autocatalytic chemical that is also activated by an external source s. For n > 1 three steady states are possible, with low and high stable steady states separated by an unstable intermediate state. Consider Figure 3.5, where for scrit1 < s < scrit2 the system is bistable. Suppose system (3.1) is initially at the low A steady state. As the signal strength s increases, the system will stay at low A steady state until s reaches the critical value scrit2, at which the response abruptly increases to the high A value. If we were to now decrease s, the system will stay at the high steady state until s reaches scrit1, when it suddenly drops to the low A steady state. This behaviour is known as hysteresis, and it has been both observed experimentally and implemented in artificial gene networks (26). For more discussion of ODE models that give rise to bistable behaviour see excellent reviews (112) and (21). Travelling waves Generally, adding diffusion to a bistable system creates a travelling wave solution sepa- rating two domains where the system resides in the two stable stationary states. Suppose the system is at one of the stable steady states. If, upon stimulus, the solution on one 46 3.2. A Survey of Models 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.5 1 1.5 2 2.5 3 Response (A) S Scrit2Scrit1 Figure 3.5: Bifurcation diagram corresponding to equation (3.1). For scrit1 < s < scrit2 the system can be either in a high or a low A steady state, depending on the initial conditions of the system.) part of the domain is kicked into the basin of attraction of the other steady state, a travelling wave is triggered, and one of the stable homogeneous states ends up taking over the entire domain. One important feature of such switch-like systems is that in order to trigger the wave from one state to another the external perturbation needs to be larger than some threshold value. This is a feature not seen in Turing systems, where arbitrarily small spatial inhomogeneities destabilize the homogeneous steady state. This is significant if we want to be able to account for resting cells. In a Turing regime, such rest states would not be seen, as they would be destabilized by the slightest noise. The speed of the travelling wave depends on the integral of the reaction kinetics: ∫ u+ u− f(u) du, (3.2) where u− and u+ are the stable steady states of the bistable function f(u) (46). That speed is generally non-zero except at isolated parameter values. Note that a stationary (pinned) wave front will subdivide the cell into multiple domains, where concentration values of the stable steady states are maintained. Spatially distributed versions of bistable systems have been investigated in many biological contexts. For example, bistability in protein-kinase cascades has been noted as a mechanism for promoting travelling waves that allow rapid long-range signalling 47 3.2. A Survey of Models across cells (65). An idea for a wave-based mechanism for directional sensing was proposed in (7). There, a LEGI adaptation module was coupled to an equation for a bistable switch. The two stable steady states (corresponding to the activated and quiescent states at the front and back of the cell) were assumed to be independent of external concentrations, while the unstable steady state was taken to depend on the ratio of activator to inhibitor in the LEGI subsystem. When the activator/inhibitor ratio equaled one a stationary front would emerge. Stationary wavefronts Stable standing waves can arise as a result of a spatially variable diffusion coefficient or a spatially inhomogeneous domain; typically studied in the context of wave propagation in excitable media (36; 57; 95). Our wave-pinning system presented in Chapter 2, provides an example of how a stable stationary wave can be achieved in a reaction- diffusion (RD) system without assumptions of spatial inhomogeneity in the domain or in the parameters. This mechanism does not need a global inhibitor, but, instead, depends on depletion of inactive form that diffuses much faster than the active form. This wave-pinning motif appears in some polarization models proposed by other groups: a model for spontaneous polarization in yeast (87), and in a model for gradient sensing in neutrophils, where another GTPase, Ras, plays the role of a gradient-sensing circuit (84), but the phenomena of wave-pinning was not specifically studied or analyzed in such papers. Reminiscent of some Turing models for polarization discussed earlier (86), the wavepin- ning system is subject to mass conservation of the components that leads to substrate depletion as the wave propagates, and eventually to the stopping of the travelling wave. However, wave-pinning is very different from a Turing instability, as the standing wave will quickly respond to external changes in the chemoattractant concentrations, whereas a Turing system often experiences “locking”. This, and other differences, are discussed further in Section 3.3. 48 3.2. A Survey of Models Note that wave-pinning subsystems can serve as a “symmetry breaking” module in larger models of cell polarity. An idea for coupling an adaptation module with a bistable system is discussed in (7). A logical extension would be to couple an adaptation module to the wave-pinning system. Wave-pinning and Turing mechanisms also need not be mutually exclusive. The interactions between multiple components can result in a model that has wave-pinning behaviour, but can also undergo a Turing instability in some parameter regimes (44). This allows both sensitivity to a small gradient, and an ability to reorient to new stimulus. The presence of additional components can also accelerate the resolution of multiple peaks (14). 3.2.5 Stochastic Models Thus far, we have focused exclusively on deterministic models. However, it must be remembered that cell polarization is inherently a stochastic process. Single molecule studies have indicated that the difference in receptor occupancy between the front and back of a cell is constantly fluctuating. These noise driven fluctuations in the number of receptors are thought to be the means by which cells break symmetry and drive cells to polarize spontaneously. In a deterministic model that is generally dealt with by adding white noise to the initial condition. That may be a reasonable approach when large numbers of molecules are involved. However, for small number of molecules, stochastic fluctuations may drive cellular behaviour to new dynamic states not seen in deterministic models. The stochastic effects of positive feedback and substrate depletion during sponta- neous polarization in yeast was investigated in (4). The deterministic version of their model has a unique spatially homogeneous equilibrium. When the deterministic system is initiated with random noise, it quickly reaches that spatial homogeneous equilibrium. However, when the stochastic model is simulated with a small number of molecules, spontaneous recurrent aggregations of molecules are observed6. This behaviour, due to 6Stochastic descriptions of spatially heterogeneous systems typically divide the reaction volume into a number of subvolumes, chosen small enough to ensure well-mixing within a compartment, and apply 49 3.2. A Survey of Models positive feedback, is not seen in the deterministic model, and as the number of molecules in the simulation increases, the probability of this type of spontaneous polarization de- creases (4). Stochastic simulations of bistable systems with equal diffusion coefficients for all components also show spontaneous separation into spatial domains of opposite phases, which is not seen in ODE or PDE formulations of equivalent systems (19). The occur- rence of this phase separation depends on the shape of the system volume, with tube-like or flat geometries undergoing phase separation more readily. Also, the diffusion coeffi- cients of the components are significant factors in the behaviour of these models. While large diffusivity has a mixing effect on the system, intermediate diffusivity coupled with feedback from the reaction (autocatalysis or double negative feedback) promotes the formation of clusters (19; 24). A mechanism for gradient sensing based on this type of phase separation pro- cess (patch coalescence and merging into a single domain) has also been proposed (15; 24; 25). In their stochastic model, the patches of PIP3 accumulate on the side of the plasma membrane with higher concentration of activated receptors, and a nucle- ation process with smaller patches being absorbed by larger patches eventually takes place. However, in that model the timescale on which polarization factors segregate into separated phases (representing spontaneous polarization) is about 100 minutes, suggest- ing that a mechanism based on patch coarsening is too slow to account for initial stages of polarity. It would be interesting to see how the behaviour of a stochastic version of the deterministic models discussed earlier depends on the number of molecules, and compare the timescale of polarization in the stochastic model to the timescale of polar- ization in the deterministic model when white noise is added to the initial condition. the reaction diffusion master equation. Diffusion between compartments is accounted for as a first- order reaction for the exchange of molecules between subvolumes at rate D/l2, where D is the diffusion constant and l is the side length of the cubic volume. 50 3.3. A Comparison of Models 3.2.6 Integrating Different Types of Models To fully elucidate the mechanisms of cell polarity requires building detailed biochemical models “from the bottom up”. Such models are based on identifiable molecular com- ponents and experimental findings, rather than on theoretical mechanisms of polarity (85). Several such models already exist. Marée et al. (64) and Dawes and Edelstein- Keshet (14) considered the interactions of Rho GTPases, phosphoinositides and the actin cytoskeleton in a variety of motile cells. Onsum and Rao (84) considered many of the same molecules (including phosphoinositides, actin, and the GTPase Ras, but not Rho GTPases) in neutrophil chemotaxis. Goryachev et al. (29) considered Cdc42 and many of its activators and effectors in yeast polarity. However, it is difficult to compare these models as they are often tailored to a specific model organism, and, by necessity, focus on just a subset of polarity regulators and experimental findings about their interactions. By understanding the necessary features for the types of mathemat- ical behaviour discussed in this review, we can categorize them by the types of modules they contain (Turing, adaptation, wave-based patterning) and identify the necessary components that give rise to their global behaviour. 3.3 A Comparison of Models In this section, we compare several models from the literature that represent the different types of models discussed earlier: the local excitation global inhibition (LEGI) adapta- tion model (52; 55), the activator substrate-depletion (ASDM) model (27; 70; 72), and a mass-conserved Turing model (29). I present a detailed comparisons of behaviour of these and our wave-pinning model. 51 3.3. A Comparison of Models 3.3.1 Simulation Results Wave-pinning Model The wave-pinning model was presented in Chapter 2, and consists of equations ∂u ∂t = Du ∂2u ∂x2 + f(u, v), (3.3a) ∂v ∂t = Dv ∂2v ∂x2 − f(u, v), (3.3b) with no-flux boundary conditions. Here u is the active form of the Rho protein, and v is the inactive form. We reproduce the numerical results from Chapter 2 in Figure 3.6, showing the response to six different type of stimuli (details of stimuli are presented in the appendix). For these simulations we used Hill function kinetics f(u, v) = v ( k0 + γu2 K2 + u2 ) − δu, (3.4) and the parameters stated in Chapter 2. For discussion of the simulation results of Figure 3.6 please see Chapter 2. The depletion of v concominant with the production of u resembles the Gierer- Meinhardt substrate-depletion mechanism (27). However, the wavepinning system does not permit a Turing diffusion-driven instability7. One manifestation of the difference of the two types of models is that wave-pinning can act much faster to generate a polar pattern than does a related Turing-type mechanism, and this is discussed further in the next section. Activator depleted substrate model (ASDM) A prototype two-variable Turing substrate depletion variant (courtesy of Prof. H. Mein- hardt) whose various terms most closely resemble those of the wave-pinning model is 7To confirm this, observe that in linearization of system 3.3 about the homogeneous stable steady states u = u−, u+, we obtain coefficients c11 ≡ fu(u, v) < 0, c22 ≡ −fv(u, v) < 0 that fail satisfy a necessary condition for Turing instabilities, namely c11/Du + c22/Dv > 0. 52 3.3. A Comparison of Models 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=20 t=200 u (a) 0 2 4 6 8 10 0 0.5 1 1.5 x t=0 t=20 t=200 (b) 0 2 4 6 8 10 0 0.5 1 1.5 x t=0 t=20 t=200 (c) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=20 t=200 u (d) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=0 t=200 t=400 u (e) 0 1 2 3 4 5 6 7 8 9 100 0.5 1 1.5 x t=200 t=0 u (f) Figure 3.6: Polarization behaviour in response to six types of distinct stimuli and/or initial conditions computed numerically for the wave-pinning system (3.3). See caption to Figure 2.8 for details. 53 3.3. A Comparison of Models given by, ∂a ∂t = b ( sa2 1 + saa2 ) − raa+Da ∂ 2a ∂x2 , (3.5a) ∂b ∂t = bb − b ( sa2 1 + saa2 ) − rbb+Db ∂ 2b ∂x2 . (3.5b) Here a is the activator and b is the substrate. Positive feedback from a to its own pro- duction is governed by a Hill function with an amplitude s and a saturation parameter sa. This production of a in Eq. 3.5a is at the expense of b in Eq. 3.5b. It is assumed that the diffusion of the substrate is greater than the diffusion of the activator (Db À Da). Superficially, this model looks very similar to the wave-pinning system with a Hill re- action term. However, the system 3.5 does not satisfy a conservation principle and consequently cannot admit a wave-pinning mechanism: the decay rates of the activator and substrate, (ra, rb) and constant input of substrate (bb) are inconsistent with simple interconversion of a and b, and the total amount is not conserved. The requirement that a polar pattern can be excited (wavenumber of the form q1 = npi/L for n = 1 on a 1D domain of length L = 10µm with no-flux boundaries) imposes strict constraints on rate constants8. The resulting system was consistent with polar patterns on this 1D domain, but the mode n = 2 was also Turing unstable, and so occasionally a higher mode (with two peaks, one at each edge) would appear. Observe 8Typical parameter values for system 3.5, s = 0.004, sa = 0.1, ba = 0, ra = 0.004, bb = 0.006, rb = 0, Du = 0.002, Dv = 0.4 were suggested by Prof. H. Meinhardt (personal communication), who also noted that rb can be set to zero with no loss of functionality. With these parameters, the spatially uniform steady state of the system, a stable spiral, occurs at u = 1.5, v = 0.81667. It is easy to check that this set of parameters satisfy the conditions for Turing pattern formation on a periodic domain as assumed in most simulations of cell polarity by Meinhardt. To choose a set of parameters for this system that would permit a reasonable comparison of behaviours to be made to the wavepinning model, we scaled the system so that rates of diffusion match those of the active and inactive Rho proteins, as in our wavepinning model (since these parameters can be reasonably determined from known biology.) Scaling time by a factor of 25 leads to the new set of parameters s = 0.1, ra = 0.1, rb = 0 (units s−1), bb = 0.15 (conc/s), sa = 0.1 (conc), Du = 0.05, Dv = 10 (units of µm2/s). The parameter sa is not changed by this rescaling, as it has units of [concentration]−2. This set of parameters is ideally suited to a periodic domain, as it leads to unstable wavenumbers of the form q = npi/L for n = 1, 2, 3, implying that not only the polar pattern, but others with up to 2 peaks could occur. By adjusting the slow diffusion to Du = 0.1µm 2/s, while retaining all other values, we found that the system was even better suited to a no-flux 1D domain. It has identical rates of diffusion as our wavepinning model and is thus suitable for numerical comparisons. Our final choice for parameters was consequently s = 0.1, ra = 0.1, rb = 0, bb = 0.15, sa = 0.1, Du = 0.1, Dv = 10, with units as above. 54 3.3. A Comparison of Models 0 2 4 6 8 100 1 2 3 x t=0 t=40 t=80 t=120 t=200 2.5 1.5 0.5 a (a) 0 2 4 6 8 10 0 0.5 1 1.5 2 2.5 3 3.5 x u t=0 t=500 t=250 (b) 0 2 4 6 8 10 0 1 2 3 x t=0 t=40 t=80 t=120 t=160 t=200 a 0.5 1.5 2.5 (c) 0 2 4 6 8 10 0 1 2 3 x t=400 t=200 t=0 2.5 1.5 0.5 a (d) 0 2 4 6 8 10 0 1 2 3 x t=0 t=80 t=160 t=200 t=400 0.5 1.5 2.5 a (e) Figure 3.7: The ASDM system (Eqs. 3.5) showing responses to the same types of stimuli used for the wave-pinning system. Spatial profiles for the activator (solid lines) at various times (t = 0, 40, . . . , 400s) are shown here to compare the time scale taken for macroscopic polarization to be evident. In (a), it takes well over 40s for the stimulus at the left cell edge to create a noticeable chemical difference between front and back of the cell, and the magnitude of polarization creeps to its full amplitude only by about t = 200s. A similar outcome occurs for a spatial gradient (c). In (b) two stimuli are applied for 20s at opposite ends of the domain and then turned off. Although the model develops two peaks, after the stimulus is removed, one of the peaks wins out and the system eventually evolves to a one peak state. The polarization persists when the stimulus is removed. (d) The polarized pattern becomes “frozen” and does not respond by repolarizing when the gradient is reversed. (e) The model admits a Turing RD instability and polarizes due to (arbitrarily) small perturbations of the homogenous steady state. 55 3.3. A Comparison of Models that rates s and ra are one order of magnitude slower than the corresponding rates γ, δ in the wave-pinning model. This choice was not arbitrary. Rather, it was a requirement for Turing instability of system 3.5. Stimuli were applied in two ways, either by adding klocS , k grad S , described in Chapter 2 to the term in large braces in system 3.5, or by setting s = S(t, x) = klocS , k grad S . Similar results were obtained for both methods. Figure 3.7 illustrates the distinct temporal and spatial behaviours of the ASDM mechanism tested with the same repertoire of stimuli in these ways. Panels (a-d) show that polarity evolves in response to both localized and graded transient stimuli, klocS , k grad S . Qualitative features of the response are distinct from those of wave-pinning : first, in the substrate-depletion model, the relative chemical levels at the two “ends” of the cell (x = 0 vs x = L) change gradually, taking many tens of seconds to create an appreciable “macroscopic” polarization. By t = 40 s, the polarity is still weak. By t = 120 s, the polarity response begins to approach its maximal amplitude. By comparison, in wavepinning, a very fast initial response creates a strong chemical difference between the cell ends (by t = 20 s or earlier) and later the propagating activation front moves further into the cell. Second, panel (d) in Figure 3.7 shows that polarity generated by the ASDM system tends to “freeze” once created9, and fails to respond by reversing when the stimulus gradient is reversed, unlike the wave-pinning model (see Figure 3.6e). Note however in panel 3.7(b) that the system can still sense multiple stimuli, developing multiple peaks. After the stimulus is removed, however, one of the peaks wins out and the system eventually evolves to a one peak state. Finally, comparison of the sensitivity of the two models to spatially noisy initial conditions (Panels e and f in the two figures) shows that ASDM is highly sensitive to signals of arbitrarily low amplitude due to its Turing instability. By comparison, response in the wavepinning model is “triggered” only by stimuli above some threshold. In response to random initial conditions, both models occasionally produce a non-polar pattern, with either more than one peak, or a peak in the center 9Generically, a third intermediate with longer time-constant is assumed to inhibit the activator and unfreeze the system in chemotactic models, e.g. (69). 56 3.3. A Comparison of Models of the domain that is a metastable state in the wave-pinning system. Local excitation global inhibition (LEGI) The LEGI model is given by ∂A ∂t = kaS(t, x)− k−aA, (3.6a) ∂I ∂t = D ∂2I ∂x2 + kIS(t, x)− k−II, (3.6b) ∂R ∂t = kRA(RT −R)− k−RIR, (3.6c) where A is the activator, I is the inhibitor and R is the response element, and S(t, x) is the external stimulus (52; 55). Equations 3.6 have a unique steady state R = A/I kr/kf +A/I . At steady state the response R is then determined by the ratio of activator to inhibitor. In the case of spatially uniform stimulus S, A = ka/k−aS, I = ki/k−iS, the ratio does not depend on extracellular signal intensity, allowing for “perfect adaptation”, as discussed earlier (see Fig. 3.8). No adaptation will occur if the kinetics for the response element are more complicated. In Figure 3.8, we demonstrate prototype responses of the LEGI model to the same repertoire of stimuli used in our wave-pinning tests10. As seen in Figure 3.8(a), a transient stimulus of the left edge of a cell, klocS . causes a transient response in R at t = 20 s, but that response decays rapidly after the signal is turned off. (Note return to homogenous distribution of R by t = 200 s). The same phenomenon is observed in Figure 3.8(c) in response to a transient graded stimulus, kgradS . When presented with two localized stimuli, a two peak response is observed that 10To generate Figure 3.8, we first set S = c0 to obtain the steady state levels, and then used transient stimuli described in the Appendix (e.g., S = c0 + c1x), together with parameter values kR = k−R = 1, ka = k−a = 2, kI = k−I = 1, RT = 1, D = 10. To represent the removal of stimulus we reset the value of S to c0. 57 3.3. A Comparison of Models quickly disappears when the stimulus is removed (Figure 3.8(b)). Note that in both the wave-pinning and the ASDM models, the system remains polarized when the stimuli are removed (whether a single or multi-peak solution remains depends on the size of the domain and the choice of the diffusion coefficient) As a direct result of the LEGI’s passive “readout” property, a model cell rapidly changes polarity in response to gradient reversal, demonstrated in Figure 3.8(d). How- ever, a spatially noisy stimulus (“perturbed rest state”) fails to produce spontaneous polarization, regardless of its amplitude. This is illustrated in Figure 3.8(e), and stems from earlier comments about the absence of a pattern-forming mechanism in the LEGI model. As shown in Fig. 3.10, the LEGI model has an increasing response to stronger gradients, without a threshold. When a stimulus is first turned on, the initial time behaviour before inhibition kicks in is approximated by γ ≈ kRAs (units of seconds). Mass Conserved Turing Model The reduced model two component mass-conserved Turing subsystem in (29) is given by ∂a ∂t = αECa2b+ βEcab− γa+Dm∆a, (3.7a) ∂b ∂t = − (αECa2b+ βEcab− γa)+Dc∆b. (3.7b) Here a is the GTP-bound form of Cdc42 (found on the membrane), b is the GDP- bound form of Cdc42 (both cytoslic and membrane bound) and Ec is the cytoplasmic Cdc24-Bem1 complex given by Ec = E0c 1 + ∫ S f(a) ds . The integral equation for Ec represents the conservation of total Cdc24-Bem1 complex. It is not required for the emergence of Turing instability in the model. However, if Ec is considered to be a constant then all of b will eventually be converted into a form. 58 3.3. A Comparison of Models 0 2 4 6 8 10 x R t=0 t=20 t=200 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (a) 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x R t=0 t=200 t=20 (b) 0 2 4 6 8 10 0 x R t=0 t=200 t=20 0.1 0.2 0.3 0.4 0.5 0.6 0.7 (c) 0 2 4 6 8 10 x R t=0 t=200 t=400 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (d) 0 2 4 6 8 10 x R t=0 t=200 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (e) Figure 3.8: The LEGI model (Eqs. 3.6), showing responses (solid lines) to the same types of stimuli and/or initial conditions (dotted lines) as for the wave-pinning model in Figure 3.6. In each case, starting at t = 0 with the steady state, a transient stimulus klocS , applied locally at the left cell edge in (a) and both edges in (b), or as a gradient (c), kgradS across the cell, leads to a transient response at t = 20 s that disappears by t = 200 s, after the stimulus is turned off. (d) Reversing the gradient leads to a reversal of polarization. (e) Random initial conditions do not lead to any spontaneous polarization, regardless of noise amplitude. 59 3.3. A Comparison of Models Stimuli were applied as for the wave-pinning model, by adding/subtracting klocS or kgradS to the RHS of system 3.7 as described for the wave-pinning model 11 and are shown in Figure 3.9. For all types of stimuli a single unique peak is established. 3.3.2 Time to Polarize Some estimates for the conditions for polarization and the temporal behaviour in the dis- tinct models is obtainable from dimensional considerations, linear stability, and asymp- totic analysis. Here, we concentrate on the specific case of diffusion coefficients relevant to Rho proteins (Da ¿ Db), for a context in which to compare the models. A compar- ison in full generality is beyond our scope, since the essential behaviour of the models is so different. A preliminary remark is that any model based on Turing instability (such as ASDM) is associated with a dispersion relation that describes the growth rate σ(k) of any wavenumber k. At the Turing bifurcation, a single critical wavenumber kcrit becomes unstable, and its corresponding growth rate σ(kcrit) becomes positive. A pattern of the form cos(kcritx) then starts to grow. A tradeoff then exists between tuning parameters to increase the growth rate of this pattern, versus preventing other wavenumbers from also becoming unstable. If the importance of having a specific pattern (e.g. a polar one) is strong, then strict limitations occur on how rapidly that pattern can grow from a small perturbation. Now consider the conditions for polarization in the wave-pinning and ASDM sys- tems. Let representative reaction rate parameters be γw, γs (units=1/time) in these mechanisms and L be the domain size. Then, in order for polarization to occur at all, 11The diffusion coefficients for this model are Dm = 0.0025µm 2/s and Dc = 1 − 10µm2/s, and the domain is considered to be a cross-section of a 6 µm sphere. I could not find a set of parameters that gave a Turing instability for Dm ≈ 0.1 in order to compare to the wave-pinning model. I used Dm = 0.0025µm 2/s and calculated the reaction parameters α, β and γ from the reaction rates given in (29). Because of the much slower diffusion coefficient in this model polarization takes much longer. For simplicity, I took Ec to be a constant. On a long timescale this results in conversion of all material in the model into the active form a. 60 3.3. A Comparison of Models 0 2 4 6 8 10 0 2 4 6 8 10 x u t=200 t=160 t=0 t=40 (a) 0 2 4 6 8 10 0 5 10 15 20 x u t=0 t=50 t=500 t=450 (b) 0 2 4 6 8 10 0 2 4 6 8 10 x u t=0 t=500 t=300 t=400 (c) 0 2 4 6 8 10 0 2 4 6 8 10 x u t=0 t=300 t=400 t=500 (d) 0 2 4 6 8 10 2.5 3 3.5 4 4.5 5 5.5 6 x u t=500 t=400 t=300 t=0 (e) Figure 3.9: Model 3.7, showing responses (solid lines) to the same types of stimuli and/or initial conditions (dotted lines) as for the wave-pinning model in Figure 3.6. In each case,for 50s, a transient stimulus klocS , applied at the left cell edge in (a) and both edges in (b), or as a gradient, kgradS across the cell in (c) and (d) In (d) the gradient is reversed at t = 250s. In all cases a unique peak is established. If the system is simulated for several thousand seconds eventually all material is concentrated in the peak. (e) The model admits a Turing RD instability and polarizes due to (arbitrarily) small perturbations of the homogenous steady state. 61 3.3. A Comparison of Models the parameters must fall approximately within the following ranges : √ Dv γw ∼ L Wave-pinning (3.8)√ Da γs ∼ L Substrate Depletion (3.9) (No such restriction exists in the LEGI model). The first condition is obtained from asymptotics and scaling in the wave-pinning models (the details of which will be pre- sented in Chapter 4); the second is a condition for Turing instability in the case Da ¿ Db. Define chemical reaction times, τw,s, and times for diffusion of the substances, τ a,b D by τw = 1/γw, τs = 1/γs, τaD = L 2/Da, τ b D = L 2/Db. Then, the above conditions can be rewritten as: τ bD ∼ τw Wave-pinning (3.10) τaD ∼ τs Substrate Depletion (3.11) In order to work, both mechanisms are predicated on a balance of a diffusion time scale with a reaction time scale, but note that the diffusive time scale in question are different in the two mechanisms since τ bD ¿ τaD. One implication here is that reaction rates must be rather slow, in keeping with the slow diffusion of the activator in ASDM, or else no polarization would be possible. By comparison, the reaction rates in the wave-pinning could be much faster, as they must only balance the diffusion rate of the rapidly-diffusing inactive form. Suppose now that the above conditions are satisfied. Then, the transition times Tw,s 62 3.3. A Comparison of Models to the polarized state can be expressed approximately as follows: Tw ∼ L√ γwDa = √ τwτaD Wave-pinning (3.12) Ts ∼ 1 γs = τs Substrate Depletion (3.13) The first of these is obtained from considerations of the speed of the wave that moves into the domain in the wave-pinning model. The second estimate is obtained from an approximation for the largest positive eigenvalue (growth rate of the most unstable mode), close to a Turing instability, in the case Da ¿ Db. Thus, in the wave-pinning mechanism, the diffusion coefficient of the inactive form dictates the possibility of the polarization response, whereas the diffusion coefficient of the active form determines the speed with which the polarized state is reached. In the substrate depletion model, the diffusion coefficient of the active form dictates the possibility of polarization whereas the speed of the polarization response is solely determined by the reaction rate. We describe a related timescale for LEGI polarization in the Appendix. From Figure 3.6a, we note yet another timescale of importance in the wave-pinning model, and that is the initial “jump” between the steady states u− and u+. This dynamic is responsible for creating the large difference in cell edges very rapidly in response to a stimulus applied at one edge. This time, TResp is also of the order TResp ∼ 1 γw = τw, Initial jump, Wave-pinning However, as noted above, in wave-pinning, a much faster reaction rate is consistent with polarization, and consequently this response time is also much faster, as discussed in the comparisons of Figure 3.6 and Figure 3.7. For example, for rates of diffusion Du = 0.1, Dv = 10µm2s−1 typical of the active and inactive Rho GTPases, with cell size on the order of L ≈ 10µm, we find the following: reaction times on the order of one to a few seconds that characterize Rho 63 3.4. Summary GTPases, would be consistent with a wave-pinning mechanism ( τw ≈ τ bD ≈ 1 − 10s). However the substrate-depletion mechanism requires that reaction times be on the order of τs ≈ τaD ≈ 1000s, likely too slow to be consistent with Rho GTPase kinetics. Furthermore, if we compare the time to complete polarization in the two mecha- nisms, we find that Tw ≈ √ τwτaD ≈ √ 1 · 1000 ≈ 30s for wave pinning. By comparison, the time to polarization in the substrate depletion mechanism (Ts ≈ τs) would be on the order of some hundreds to thousand seconds (about 20 min), which is slower than seen experimentally. As mentioned above, the wave-pinning mechanism has an even faster response when we consider the fact that an initial jump already defines a clear polarity even before the final polarized state is attained. This occurs on a timescale of τw, i.e. within a few seconds. There is no comparable “fast” response in the substrate-depletion model. These rough estimates illustrate the distinctions between the mechanisms and suggest that the wave-pinning mechanism may be inherently better for polarizing a cell rapidly. 3.3.3 Model Comparison Summary In Figure 3.10, we compare the responses of three models to both persistent (left) and transient (right) stimuli of various gradient strengths. System (3.7) was not consid- ered, as it was developed to explain spontaneous polarity in yeast, and not to detect chemotactic stimuli. Several differences are noteworthy. The wave-pinning model (top) requires gradients larger than some threshold to respond, unlike both ASDM and LEGI models. The LEGI model (middle) responds in a way that increases with stimulus strength, but forgets the stimulus once it is turned off. 3.4 Summary In this Chapter I reviewed other mathematical models from the cell polarization liter- ature and compared them to the wave-pinning model in Chapter 2. I subdivided the models into several classes based on proposed mechanism for polarization. I compared 64 3.4. Summary 0 2 4 6 8 10 0 1 x a 1.5 0.5 (a) 0 2 4 6 8 100 1 x 1.5 0.5 a (b) 0 2 4 6 8 10 x R 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (c) 0 2 4 6 8 100 x R 0.1 0.2 0.3 0.4 0.5 0.6 0.7 (d) 0 2 4 6 8 100 1 2 3 x 0.5 1.5 2.5 a (e) 0 2 4 6 8 100 1 2 3 x 0.5 1.5 2.5 a (f) Figure 3.10: Comparison of sensitivity and persistence of the wave-pinning model (top), the LEGI model (center) and the substrate-depletion model (ASDM) (bottom) in response to uniform and successively stronger graded stimuli. In each panel we show the profile (of the active form a in system 3.3, of the response element R in system 3.6, and of the activator a in system. 3.5) at t = 200 in response to the stimuli S(x) = 0.01(1+m)x wherem = 0, 1, · · · , 9. Left panels: gradient on for entire duration; Right panels: transient gradient, turned off at t = 50. The ASDM (e,f) is most sensitive: it responds to all but the uniform stimulus by polarizing. The wave-pinning system (a,b) exhibits a threshold, and LEGI produces a response that increases with signal. Both WP (b) and ASDM (f) but not LEGI (d) persist after the stimulus is removed. 65 3.4. Summary Table 3.1: Features of polarity explained by various models. Behaviour ASDM Wave-pinning LEGI Maintenance of polarity Yes Yes No Sensitivity Yes Yes Yes Multi-stimuli response Yes Yes Yes Amplification Yes Yes No Adaptation No No Yes Spontaneous Polarization Yes Yes No Reversible Asymmetry No Yes Yes three prototypical model types: the wave-pinning model, the activator substrate deple- tion Turing model (ASDM), and the LEGI model to various stimuli. The behaviours of these models are summarized in Table 3.1, where the features of cell polarization de- scribed in Chapter 1 are listed. I found that no class of models can explain all features of cell polarity. In particular, while Turing and wave-based models are able to maintain polarity once the stimulus is removed, the LEGI model does not. The LEGI model, on the other hand, exhibits adaptation to spatially uniform stimulus, which the other two models do not. The major difference between the wave-pinning and the Turing model is the ability to reverse asymmetry. While the wave-pinning model is able to change direc- tions, and even return back to a homogeneous steady state, with appropriate stimulus, this is not possible in the Turing model. The question of biological relevance of the model is also an important one. The Tur- ing models seem best suited for irreversible polarity, such as the formation of the Cdc42 cap in yeast. LEGI models are able explain the response to stimulus of Dictyostelium cells in the absence of the actin cytoskeleton. The features of the wave-pinning model make it most applicable to spontaneously polarizing cells that can still rapidly respond to a change in stimulus, such as neutrophils. 66 Chapter 4 Asymptotic Analysis of the Wavepinning Mechanism 4.1 Introduction In Chapter 2, I formulated the wave-pinning model, and described the behaviour of the model in the context of cell polarity. In this chapter, I present a detailed mathematical analysis of wave-pinning mechanism under the assumption that Dv À Du, i.e. the inactive form diffuses much faster than the active form. Under this assumption, the wavepinning system can be rescaled to obtain a small parameter ². For ²¿ 1 singular perturbation theory can be used to obtain matched asymptotic expansions for u and v (the active and inactive forms of the GTPase). The wave-front can be viewed as an interior transition layer, where u undergoes a rapid change. Using standard techniques from boundary layer theory (32) allows us to quantitatively predict the location in the domain where the wave will become pinned. 4.1.1 Model Equations Recall that the model for Rho protein concentration in the active, membrane-bound form (u) and the inactive cytosolic form (v) is given by the following system of reaction- diffusion (RD) equations ∂u ∂t = Du∆u+ ηf(u, v), (4.1a) ∂v ∂t = Dv∆v − ηf(u, v), (4.1b) 67 4.1. Introduction where f(u, v) is the rate of interconversion of v to u and η is the speed of the reaction. (Note that η has units of s−1). By previous remarks, rates of diffusion satisfy Du ¿ Dv. The boundary conditions at the ends of the domain are no-flux ∂u ∂x = ∂v ∂x = 0, x = 0, L (4.1c) resulting in mass conservation. That is, the total amount of protein in the domain, Ktotal is constant: ∫ Ω (u+ v)dx = Ktotal. (4.2) For the kinetics function, f(u, v), we require only that there is some range of values of v (vmin < v < vmax) for which f(u, v) = 0 admits three solutions as an equation in u. In order of increasing size, we denote these by u−(v), um(v), and u+(v). Here, u− and u+ are assumed to be stable solutions to the well-stirred (ODE) version of (4.1), for a fixed v, while um(v) is unstable. As a simple example, we use the cubic polynomial f(u, v) = C u(m− u)(u−m− v), (4.3) whose roots are u− = 0, um = m, and u+(v) = m + v, where m > 0 is a constant with dimensions of concentration (defined as above). We define C = 1/m2 > 0 to be a constant carrying units so that dimensions of ηf (concentration/time) match those of ut, vt on the LHS of (4.1). In Chapter 2, I used a more biologically relevant (but less mathematically conve- nient) biochemical kinetics function f(u, v) = v ( k0 + γu2 m2 + u2 ) − δu, (4.4) where k0, γ,m, δ ≥ 0 are constants. However, since the phenomenon of wavepinning does not depend on the exact form of the kinetics function, but only the assumptions 68 4.2. Intuitive Explanation for Wave-Pinning about bistability stated above, the cubic kinetics are chosen for the purpose of analytic tractability. Equation (4.3) can be viewed as a simplified version of (4.4) that retains the necessary features for wavepinning, much like the Fitzhugh-Nagumo model is a reduction of the more biologically realistic Hodgkin-Huxley model. Most of the analysis in this chapter does not depend on details of the function (4.3), and is applicable to other bistable kinetics f(u, v), but using cubic kinetics (4.3) allows one to obtain an exact analytic solution for the speed of the wavefront, which is derived in Section 4.4. 4.2 Intuitive Explanation for Wave-Pinning 4.2.1 The Single Variable Case (fixed v) To understand why wave-pinning occurs, we first briefly recall the following background. Consider a single one-dimensional reaction diffusion equation for u on 0 < x < L, with the above bistability property, and with v taken as parameter: ∂u ∂t = Du ∂2u ∂x2 + f(u, v). (4.5) On an infinite domain, Eq. (4.5) is known to possess a propagating front solution uf (ξ) = uf (x− ct)) with uf → u± as ξ → ±∞, (i.e. rear and front asymptotically approach the homogeneous steady states of f(u, v)) (46; 77). At the rear (respectively, front) of the wave (x→ ±∞) the profile asymptotes to u+ (respectively, u−). Suppose u, obeying Eq. (4.5) on 0 ≤ x ≤ L, develops a propagating front that is sufficiently sharp and far from boundaries; then we can approximate its region of validity by −∞ < x <∞. Let uf (ξ) = uf (x− ct) = u(x, t), where c is the speed of the front, which need not be constant. Then it is known (46; 77) 69 4.2. Intuitive Explanation for Wave-Pinning u v x u 0 Lx Reaction Rates Reaction Rates Reaction Rates u u 0 0 L L x u u v vc (a) Figure 4.1: An intuitive explanation of wavepinning based on the Maxwell condition The left column plots the concentrations of u and v versus position x whereas the right column plots the reaction terms (LHS and RHS of Eq. 2.8) as functions of u. The intersections denote the steady states u−, um, and u+. At the beginning, the value of the integral I(v) (numerator of Eq. 4.6 and difference of shaded areas) is positive, and the wave front moves to the right (top figure). As the wave front sweeps across the domain, v and hence I(v) decreases, and the wave decelerates (middle). The wave eventually comes to a halt when shaded areas are equal (I(vc) = 0, bottom). 70 4.2. Intuitive Explanation for Wave-Pinning that c satisfies c(v) = ∫ u+ u− f(u, v) du∫∞ −∞ (∂uf (ξ)/∂ξ) 2 dξ . (4.6) c depends on v through its dependence on u−, u+, f(u, v) and uf . The denominator of Eq. (4.6) is positive, so the sign of c(v) depends only on the integral in the numerator, which we denote hereafter by I(v). For the reaction term in Eq. (4.4), the quantity I(v) can be represented geometrically. For example, in the top right panel of Fig. 4.1, there are two regions between these intersecting curves (shown in different shades), one clearly having an area of greater magnitude than the other. For our reaction term, there is one critical value v = vc, (in vmin < vc < vmax) at which the areas described above are equal and opposite. At this value of v, the speed c(v) changes sign. If v < vc, c is negative and if v > vc, c is positive. The speed is zero when the following condition (called a Maxwell condition (47; 51)) is satisfied: I(vc) ≡ ∫ u+ u− f(u, vc) du = 0. (4.7) In general, the wave speed c can be zero for isolated parameter values (i.e., v = vc), but zero wave speed solutions are structurally unstable (46; 77). Hence, these are not suitable for robust polarization. 4.2.2 Wave-Pinning in the Full u, v System: The phenomenon of wave-pinning critically depends on the three qualitative features of our model equations (4.1): 1. Mass conservation results from the interconversion kinetics and no-flux bound- ary conditions, as given by Eqs. (4.1) and (4.2). The total amount of protein, Ktotal > 0 is constant even when an external signal is applied, (kS 6= 0) since such signals simply increase the rate of conversion from v to u without affecting the combined total amount of u and v in the domain. 2. The component v diffuses quickly relative to u (since Dv À Du). This means that 71 4.2. Intuitive Explanation for Wave-Pinning v is spatially uniform (to first approximation), though not constant in time. De- pletion of v forms the key reason for stalling the wave. This amounts to replacing v in Eq. (4.1) by the average value. This implies that the inactive form, Rho-GDP acts as a global variable. 3. The assumed bistability of the reaction term f(u, v). That is, there are two stable homogeneous state (u, v) ≡ (u±(v), v) of the reaction diffusion system for a range of v, vmin < v < vmax. Suppose that initially u is at its basal steady state level and v is sufficiently high (u(x, 0) = u−(v), v(x, 0) ≈ constant > vc). A transient stimulus near the edge of the cell at x = 0 triggers a local increase to the level u ≈ u+, via bistability. According to Eq. 4.6, and using the fact that v(0) > vc, this creates a front that starts to propagate in the positive direction, converting a greater fraction of the domain to u ≈ u+(v) (see Figure 4.1). By Eq. 4.2, as the amount of u increases, v must simultaneously decrease, eventually reaching the critical value v = vc at which wave-pinning occurs. The depletion of v to its critical level vc causes the wave to pin. Note that wave-pinning would not occur if there is initially too little inactive form, e.g., v < vc at t = 0. If initially v is too high, if v À vc, wave-pinning would fail as the reduction in v as the front sweeps through the domain is insufficient to reach the level v = vc inside the given domain. For fixed domain size L, the possibility of wave-pinning is determined by the total amount, Ktotal, of u and v. The model predicts that the size of the activated region (i.e. portion of the domain for which u is at its high steady state value) will be an increasing function of the total amount Ktotal. Because the two homogeneous steady states do not lose stability, there is a range of parameters ( i.e., a range of values of Ktotal) where both unpolarized and polarized states coexist. The system asymptotically approaches either one or the other configuration depending on initial conditions and on the applied transient stimulus kS . Calculations of the front speed, and the position where the front pins as a function of the total amount are presented in the next section. 72 4.3. Matched Asymptotic Expansion in 1D 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x t=0 t=1 t=10 (a) Figure 4.2: Wavepinning in 1D for system (4.10) with reaction kinetics (4.12). Param- eters used: ² = 0.1,D = 1µm2/s. The active form u is initialized to be u = u+(v) on 0 < x < L/3 and u = u−(v) on L/3 < x < L (dotted line). This sets up a travelling wave at one end of the domain that moves inward, slows down and freezes by t = 10 s. 4.3 Matched Asymptotic Expansion in 1D 4.3.1 Scaling the Equations Consider a typical scale K0 for the concentration (e.g., K0 = m for (4.3)), a scale T0 for time, and a length scale, L given by the cell diameter. Recall that η is a typical rate of the biochemical reaction of interest (in s−1). Then equations (4.1)can be scaled so that u = K0ũ, v = K0ṽ, x = Lx̃, t = T0t̃, (4.8) where ũ, ṽ, x̃, and t̃ are dimensionless variables. In the dimensionless spatial variable, the domain is Ω′ = {x̃ | 0 ≤ x̃ ≤ 1}. We define: ²2 = Du ηL2 , D = Dv ηL2 , T0 = L√ ηDu . (4.9) Then, based on the assumed unequal rates of diffusion, Du ¿ Dv, ² is a small quantity whereas D = O(1) with respect to ². This is equivalent to saying that √Dv/η ≈ L, i.e. 73 4.3. Matched Asymptotic Expansion in 1D on the timescale of the biochemical reaction, the inactive form can diffuse across the diameter of the cell. Assuming a typical cell diameter of 10 µm, reaction timescale η = 1 s−1, and diffusion coefficientsDu = 0.1µm2s−1 andDv = 10µm2s−1, the dimensionless constants are ² ≈ 0.03 and D ≈ 0.1. Substituting these variables into (4.1) and dropping the ˜ yields the dimensionless equations: ² ∂u ∂t = ²2 ∂2u ∂x2 + f(u, v), (4.10a) ² ∂v ∂t = D ∂2v ∂x2 − f(u, v). (4.10b) The boundary conditions are: ∂u ∂x ∣∣∣∣ x=0 = ∂u ∂x ∣∣∣∣ x=1 = 0, ∂v ∂x ∣∣∣∣ x=0 = ∂v ∂x ∣∣∣∣ x=1 = 0. (4.10c) The total (dimensionless) amount of protein satisfies ∫ 1 0 (u+ v)dx = K. (4.11) where K = Ktotal/K0 (e.g. K = Ktotal/m for (4.3).) Note that after rescaling, the cubic kinetic term (4.3) becomes f(u, v) = u(1− u)(u− 1− v). (4.12) However, as stated earlier, most of the analysis in this section does not depend on the particular form of f(u, v), but only on the fact that f(u, v) is bistable as an equation in u. It is assumed that f(u, v) satisfies the following properties: 1. In the range vmin ≤ v ≤ vmax, the equation f(u, v) = 0 has three roots, u−(v) < um(v) < u+(v). 2. Keeping v fixed, u−(v) and u+(v) are stable fixed points of the equation ut = 74 4.3. Matched Asymptotic Expansion in 1D f(u, v). That is, f(u±(v), v) = 0, and, for stability, ∂f ∂u ∣∣∣∣ u=u±(v) < 0. (4.13) 3. The homogeneous states, (u, v) ≡ (u±(v), v) are stable states of the PDE (4.10a,b), requiring that ( ∂f ∂u − ∂f ∂v )∣∣∣∣ u=u±(v) > 0. (4.14) See Appendix B.2 for a derivation of this condition. 4. There is one value v = vc at which the integral I(v) becomes 0, i.e., I = ∫ u+(vc) u−(vc) f(u, vc) du = 0 (4.15) We assume that I changes sign at vc and that u+(v), um(v), and u−(v) depend on v in such a way that I > 0 for v > vc and I < 0 for v < vc. This condition is satisfied in (4.12) as u+(v) is an increasing value of v (u+(v) and um(v) are independent of v in that case). It is assumed that u and v have the form of a front solution. Let φ(t) denote the (dimensionless) location of this front (let xf (t) be the corresponding dimensional quantity). The diffusion coefficient determines the width of the transition zone. Because ² is small, while D is O(1), u is expected to form a sharp front at φ(t), whereas v would be much smoother. Next, I perform a matched asymptotic calculation to examine such solutions. In the next section I analytically calculate the front speed, and the position where the front will pin, thus reducing a system of two partial differential equations to a single ODE for the propagating front position. 75 4.3. Matched Asymptotic Expansion in 1D 4.3.2 Asymptotic Analysis to Leading Order Outer Solution Consider equations (4.10) in the outer region, that is, away from the front at φ(t). Assume that u has the expansion u = u0+²u1+²2u2 · · · and likewise for v. Substituting these expansions into (4.10)a,b and retaining terms of order ²0 leads to the following equations for the leading order terms, u0 and v0: 0 = f(u0, v0), (4.16a) 0 = D ∂2v0 ∂x2 − f(u0, v0). (4.16b) Equations (4.16) are valid in the region 0 ≤ x < φ(t)−O(²) and φ(t) +O(²) ≤ x < 1, that is, at some small distances to the left/right of the sharp transition zone at the front. Adding (4.16a,b), we find: D ∂2v0 ∂x2 = 0. (4.17) From (4.10c) and (4.17), we conclude that: v0(x, t) = v<(t) 0 ≤ x < φ(t)−O(²), v>(t) φ(t) +O(²) < x ≤ 1. (4.18) where v> and v< (the values of v to the right and to the left of the front) do not depend on x. (Later, it will be shown that v> = v<.) From (4.16a), u0 takes on one of the values u+, u− or um in the outer regions. Since we are seeking a stable front solution, we let: u0(x, t) = u+(v<) 0 ≤ x < φ(t)−O(²), u−(v>) φ(t) +O(²) ≤ x < 1. (4.19) 76 4.3. Matched Asymptotic Expansion in 1D Inner Solution Let w = x− φ(t) (4.20) be the distance from the transition front, i.e. w is a coordinate moving with the front. Introduce a stretched coordinate ξ for the inner layer, in a small strip close to the evolving front: ξ = w ² = x− φ(t) ² . (4.21) Let U(ξ, t) = u((x− φ(t))/², t), (4.22a) V (ξ, t) = v((x− φ(t))/², t), (4.22b) be the inner solution. Then (4.22) is not a traveling front solution in the strict sense, as the speed dφ/dt is not constant. Moreover, the amplitudes of U and V change with time. Hence, we do not assume u(x, t) = U(ξ), but rather u(x, t) = U(ξ, t), and likewise for V . Substitute the new scaling into (4.10) to obtain ² ∂U ∂t − dφ dt ∂U ∂ξ = ∂2U ∂ξ2 + f(U, V ), (4.23a) ² ∂V ∂t − dφ dt ∂V ∂ξ = D²−2 ∂2V ∂ξ2 − f(U, V ). (4.23b) We again assume usual expansions for inner solutions U and V in powers of the small parameter, ². The domain for these equations is −∞ < ξ <∞. Let φ = φ0+ ²φ1+ · · · . Note that dφ/dt is the speed of the front. Then to leading order, we find: −dφ0 dt ∂U0 ∂ξ = ∂2U0 ∂ξ2 + f(U0, V0), (4.24a) ∂2V0 ∂ξ2 = 0. (4.24b) 77 4.3. Matched Asymptotic Expansion in 1D From (4.24b), it follows that V0 = a1(t)ξ + a2(t), (4.25) where a1(t), a2(t) are arbitrary functions of t only, to be determined. Matching Inner and Outer Solutions We must match V0 with the values of v obtained above in the outer regions, i.e., lim ξ→−∞ V0(ξ) = v<, lim ξ→∞ V0(ξ) = v>. (4.26) For these limits to exist V0 must be a constant in the inner layer. We conclude that v> = v< = V0. (4.27) Thus, V0 is spatially uniform throughout the domain and is the same as the outer solution v0. We drop the dependence of v0 on x (or V0 on ξ). We now construct a solution for U0 in the inner layer. Since V0 is spatially constant in the inner layer, equation (4.24a) is an equation in U0 only, where V0 is a parameter (that varies in time). We must solve the boundary value problem (4.24a) with the matching conditions from (4.19) as boundary conditions at ±∞: lim ξ→−∞ U0(ξ) = u+(V0), lim ξ→∞ U0(ξ) = u−(V0). (4.28) Such a heteroclinic solution, unique up to translation, exists for general bistable reaction terms f(U, V ) (46; 77). Denote this solution by Uφ0 (ξ, V0). Multiplying (4.24a) by ∂Uφ0 /∂ξ and integrating from ξ = −∞ to ξ =∞, we obtain: dφ0 dt = c(V0) ≡ ∫ u+(V0) u−(V0) f(s, V0)ds∫∞ −∞ ( ∂Uφ0 (ξ, V0)/∂ξ )2 dξ . (4.29) By previous assumption dφ0/dt changes sign at V0 = vc so that dφ0/dt has the same 78 4.3. Matched Asymptotic Expansion in 1D sign as V0 − vc. Determining the Amount of Inactive Protein at the Interface to Zeroth Order (V0) by a Solvability Condition We next determine the unknown constant V0 by considering the solvability condition for the higher order outer equations (see Section B.4 for details on solvability conditions). Consider equations (4.10a,b) to the next order. Isolating terms of order ² leads to ∂u0 ∂t =fu(u0, v0)u1 + fv(u0, v0)v1 (4.30a) ∂v0 ∂t =D ∂2v1 ∂x2 − fu(u0, v0)u1 − fv(u0, v0)v1. (4.30b) Adding equations (4.30) we obtain ∂ ∂t (u0 + v0) = D ∂2v1 ∂x2 . (4.31) Integrating the above from x = 0 to x = 1 results in ∂ ∂t (∫ 1 0 u0 dx+ v0 ) = 0, (4.32) where we used the boundary conditions (4.10c)and the fact that v0 is a constant through- out the domain. Equation (4.32) represents a solvability condition for equation (4.31). Integrating the above in t, leads to v0 = K − ∫ 1 0 u0 dx (4.33) where K, as before, corresponds to the (dimensionless) total amount of u and v in the domain. The integral of u0 can be approximated by contributions from the two outer 79 4.3. Matched Asymptotic Expansion in 1D regions (to left and right of the front) and the inner region: ∫ 1 0 u0 dx = ∫ φ(t)−O(²) 0 u0 dx+ ∫ 1 φ(t)+O(²) u0 dx+O(²) = u+(v0)φ(t) + u−(v0)(1− φ(t)) +O(²) where we have used (4.19) in the second equality. We ignore terms of O(²), in the leading order approximation. The reaction-diffusion system is then reduced to the following ordinary-differential- algebraic system: dφ0 dt = c(v0), (4.34a) v0 = K − u+(v0)φ0 − u−(v0)(1− φ0). (4.34b) In (4.34b) v0 is the amount of inactive form left over, once a plateau of magnitude u+ and width φ (and a second one of magnitude u− and remaining width 1 − φ) is subtracted from the total amount of material, K. 4.3.3 Travelling Wave Pins by Depletion of Inactive Protein We now show that the front speed, dφ0/dt, and the rate of change dv0/dt, have opposite signs. Differentiating the relation f(u±(v), v) = 0 with respect to v and using (4.14) leads to 0 = ( ∂f ∂u du± dv + ∂f ∂v )∣∣∣∣ u=u±(v) > ( 1 + du± dv ) ∂f ∂u ∣∣∣∣ u=u±(v) . (4.35) Using (4.13) we conclude that: 1 + du± dv > 0. (4.36) Differentiate (4.34b) with respect to t: ( 1 + du+(v0) dv φ0 + du−(v0) dv (1− φ0) ) dv0 dt = −(u+(v0)− u−(v0))dφ0 dt (4.37) 80 4.3. Matched Asymptotic Expansion in 1D Since the front position φ must reside within the domain (0 < x < 1), we have 0 < φ0 < 1. Using this and (5.25), we see that the factor multiplying dv0/dt in (4.37) is positive. Since (u+ − u−) > 0, we conclude from (4.37) that dv0/dt and dφ/dt have opposite signs. Thus, v0 is depleted as the position of the wave, φ moves forward. Suppose v is sufficiently large initially, i.e., v0 > vc at t = 0. Since dφ/dt is positive for v0 > vc, dv0/dt < 0. Thus v0 decreases as the front φ advances further into the domain. The value of v0 asymptotes to vc, at which point the front φ is pinned. The position of the pinned front can then be predicted if the form of the reaction kinetics and the total amount of material in the domain (determined by the initial conditions of the PDE system) is known. This pinned front is stable; if the front is forced to move away from its pinned position, it will relax back to it. Denote the pinned position by φp. Then, vc = K − u+(vc)φp − u−(vc)(1− φp). (4.38) Since vc is a constant that depends only on f(u, v), we can interpret (4.38) as a relation between φp and K. We say that the wave is pinned if the front stalls inside the domain, i.e., provided 0 < φp < 1. This implies that vc + u−(vc) < K < vc + u+(vc) (4.39) is a condition on K for wave pinning to occur. In general, an explicit form for the propagating front solution in the denominator of (4.29) is not obtainable. However, the simplicity of the cubic reaction kinetics (4.12) allows us to evaluate the integrals in (4.29) explicitly as a function of v (77). We obtain an explicit expression for the (dimensionless) front velocity as a function of v, dφ dt = c(v0) = 1√ 2 (u+(v0)− 2um(v0) + u−(v0)) . (4.40) The front becomes stationary when dφ/dt = 0. 81 4.3. Matched Asymptotic Expansion in 1D 4.3.4 Asymptotic Analysis to Higher Order We now obtain approximations to u, v, and φ to the next order in ². This gives a measure of the error of the leading order approximation. The procedure to obtain the second term is similar to finding the first, so for brevity only the major steps are presented here. Outer solution for the first order correction v1 of the inactive protein v The outer equations are given by (4.30). Summing them up we obtain (4.31). Note u0 = u±(v0) and that v0 is spatially constant within the regions x < φ(t) and x > φ(t). From (4.10c), we see that v1 satisfies no flux boundary conditions at x = 0 and x = 1. Integrating equation (4.31) twice and using the boundary conditions we have: v1 = 1 D ( 1 + du+(v0)dv ) dv0 dt x2 2 +B+ = A+x 2 +B+ for x < φ(t) 1 D ( 1 + du−(v0)dv ) dv0 dt (1−x)2 2 +B− = A−(1− x)2 +B− for x > φ(t) (4.41) where A+ = 1 2D ( 1 + du+(v0) dv ) dv0 dt , A− = 1 2D ( 1 + du−(v0) dv ) dv0 dt , (4.42) and B+ and B− are constants to be determined. Substituting in x = ²ξ + φ we rewrite the two term outer expansion in terms of the inner variable ξ to get: v0+²v1 = V0 + ²(A+(²ξ + φ)2 +B+) = V0 + ²(A+φ2 +B+) + · · · 0 ≤ x < φ(t), V0 + ²(A−(1− φ− ²ξ)2 +B−) = V0 + ²(A−(1− φ)2 +B−) + · · · φ(t) < x ≤ 1. (4.43) 82 4.3. Matched Asymptotic Expansion in 1D Matching the inner (V1) and outer (v1) first order terms for inactive protein v In the inner layer, V1 satisfies the following equation, as can be seen from (4.23b) ∂2V1 ∂ξ2 = 0. (4.44) From this we conclude that V1(ξ) has the form: V1(ξ) = A(t)ξ +B(t). (4.45) Rewriting the two term inner expansion in terms of the outer variable x get: V0 + ²V1 = V0 + ²(Aξ +B) = V0 +A(x− φ) + ²B, (4.46) where we used the fact that v0 = V0 is a constant determined previously from the conservation condition. van Dyke’s matching principle states that the power series expansion of the nth order inner solution Y ni (ξ(x)) must match the power series expansion of the nth order outer solution yno (z(ξ)) (45). This means that the inner layer solution expansion in terms of the outer variable and the outer solution expansion in terms of the inner variable must agree. It follows that for expansions (4.43) and (4.46) to be equivalent to first order, we need A = 0 (i.e., V1 must be a constant in ξ in the inner layer), and V1 = B(t) = A+φ2 +B+ = A−(φ− 1)2 +B−. (4.47) We can now rewrite equations (4.41) as v1 = A+(x2 − φ2) + V1 for x < φ(t) A−[(1− x)2 − (1− φ)2] + V1 for x > φ(t) (4.48) 83 4.3. Matched Asymptotic Expansion in 1D The unknown constant V1 will be later determined from the conservation condition (2.10). Determining first order correction for active protein u Let us now consider matching for u1. Equation (4.30a) provides an algebraic relation for u1 in the outer layer: u1 = ∂u0 ∂t − fv(u0, v0)v1 fu(u0, v0) . (4.49) Since we already have an expression for v1, we may substitute this to find the following expression: u1 = 1 fu(u+(v0),v0) du+(v0) dt − fv(u+(v0),v0)fu(u+(v0),v0)(A+x2 +B+) for x < φ(t) 1 fu(u−(v0),v0) du−(v0) dt − fv(u−(v0),v0)fu(u−(v0),v0)(A−(1− x)2 +B−) for x > φ(t). (4.50) Now, consider U1 in the inner layer: ∂2U1 ∂ξ2 + dφ0 dt ∂U1 ∂ξ + fu(U0, V0)U1 = ∂U0 ∂t − dφ1 dt ∂U0 ∂ξ − fv(U0, V0)V1. (4.51) Let us consider the matching conditions to obtain the boundary conditions to be satisfied by U1. Just as in the v1 case, we match the outer and inner solutions using the van Dyke matching principle. Rewrite the two term outer expansion for u: (4.19) and (4.50) in terms of the inner variable to get u0 + ²u1 = u+ + ²G+ + · · · for x < φ(t) u− + ²G− + · · · for x > φ(t), (4.52) where G+ = 1 fu(u+(v0), v0) du+(v0) dt − fv(u+(v0), v0) fu(u+(v0), v0) V1(t), (4.53a) G− = 1 fu(u−(v0), v0) du−(v0) dt − fv(u−(v0), v0) fu(u−(v0), v0) V1(t). (4.53b) 84 4.3. Matched Asymptotic Expansion in 1D Determining the speed of the interface to first order, dφ1dt In order to determine dφ1dt we apply a solvability condition to (4.51). As lim ξ→−∞ U1(ξ) = G+ and lim ξ→∞ U1(ξ) = G− (4.54) equation (4.51) has inhomogeneous boundary conditions. We first consider the related homogeneous problem. Consider the following function Ũ0: Ũ0 = G+ +G− 2 + G+ −G− u+(v0)− u−(v0) ( U0 − u+ + u−2 ) = H1 +H2U0. (4.55) This function satisfies the boundary conditions (4.54) and the equation: ∂2Ũ0 ∂ξ2 + dφ0 dt ∂Ũ0 ∂ξ +H2f(U0, V0) = 0. (4.56) Therefore, the function W = U1 − Ũ0 satisfies the equation: ∂2W ∂ξ2 + dφ0 dt ∂W ∂ξ +fu(U0, V0)W = ∂U0 ∂t −dφ1 dt ∂U0 ∂ξ −fv(U0, V0)V1−fu(U0, V0)Ũ0+H2f(U0, V0) (4.57) where W → 0 for ξ → ±∞. Let: LW = ∂ 2W ∂ξ2 + dφ0 dt ∂W ∂ξ + fu(U0, V0)W. (4.58) The adjoint operator to L is: L∗W = ∂ 2W ∂ξ2 − dφ0 dt ∂W ∂ξ + fu(U0, V0)W. (4.59) Let p(ξ) = ∂U0∂ξ . Take the derivative of the U0 equation (4.24a) with respect to ξ to find: Lp = 0. (4.60) 85 4.3. Matched Asymptotic Expansion in 1D Therefore, p∗(ξ) = p(−ξ) satisfies L∗p∗ = 0. Now ∫ ∞ −∞ LW p∗ − L∗p∗W = ∫ ∞ −∞ p∗ ( ∂2W ∂ξ2 + dφ0 dt ∂W ∂ξ + fu(U0, V0)W ) −W ( ∂2p∗ ∂ξ2 − dφ0 dt ∂p∗ ∂ξ + fu(U0, V0)p∗ ) = [ ∂W ∂ξ p∗ − ∂p ∗ ∂ξ W + dφ0 dt Wp∗ ]∞ −∞ = 0, (4.61) as W = 0 at ξ = ±∞. As L∗p∗ = 0 this implies that ∫ ∞ −∞ LW p∗ = 0. (4.62) This yields the following condition for dφ1/dt: dφ1 dt (∫ ∞ −∞ ∂U0 ∂ξ p∗dξ ) = ∫ ∞ −∞ ( ∂U0 ∂t − fv(U0, V0)V1 − fu(U0, V0)Ũ0 +H2f(U0, V0) ) p∗dξ. (4.63) The same result can be obtained by applying a solvability condition (see Appendix B.4) to (4.58). Determining the unknown constant V1 We have determined u1, v1 and φ1 up to a time dependent constant V1(t). In order to determine this constant, we again use the conservation condition (4.11): K = ∫ 1 0 (u+v)dx = u+(v0)(φ0+²φ1)+u−(v0)(1−(φ0+²φ1))+v0+² ∫ 1 0 (u1+v1)dx+O(²2). (4.64) We know that the leading order term is equal to K. Therefore, we obtain the relation: u+(v0)φ1 − u−(v0)φ1 + ∫ 1 0 (u1 + v1)dx = 0. (4.65) 86 4.4. Example with Cubic Kinetic Function This gives us an algebraic relation between V1(t) and φ1(t), and together with (4.63) gives us the complete description of the approximate solution up order one in ². 4.4 Example with Cubic Kinetic Function Leading order approximation to the front position For kinetics in (4.3) with the concentration scale K0 = m, the dimensionless version of the kinetics is given by f(u, v) = u(1− u)(u− 1− v). (4.66) Thus, we have u− = 0, um = 1, u+ = 1 + v in dimensionless quantities, so that, using (4.40), the leading order equations (4.34) become: dφ0 dt = v0 − 1√ 2 , (4.67) v0 = K − (1 + v0)φ0. (4.68) From (4.68), we find that the wave stops when v0 = 1 ≡ vc. This is also seen by a simple symmetry argument based on the Maxwell condition (4.15), since this value of v produces a cubic with roots that are symmetric about the root um = 1. The corresponding dimension-carrying roots are at u = 0,m, and 2m, i.e. v = m when the wave stalls. Thus, (4.39) reduces to 1 < K < 3, orm < Ktotal/L < 3m, in corresponding dimension-carrying quantities. Solving (4.68) for v0 obtain v0 = K − φ0 1 + φ0 . (4.69) Hence (4.40) becomes dφ0 dt = 1√ 2 ( K − φ 1 + φ − 1 ) . (4.70) The position at which the wave stalls (in dimensionless and dimensioned forms, respec- 87 4.4. Example with Cubic Kinetic Function tively) is given by ∂φ0 ∂t = 0, i.e., (4.71) φp = K − 1 2 , or xp = Ktotal − Lm 2m . (4.72) Comparison with numerical results Numerical solutions to the full PDE system (4.10) with reaction kinetics (4.12) and the ODE (4.70) predicting front position are compared in Fig. 4.3. The exact front position is calculated from the numerical solution of the PDE system by tracking the position φnum at which u = um(v) (u = 1 for reaction kinetics (4.12)). φnum(t0) is used as the initial condition, where t0 ≈ 0 is a time at which the solution to the PDE system has relaxed to assume the form assumed in the asymptotic calculations. Numerical results agree with the asymptotic theory to order ² (see Fig. 4.4). Note that the error between the asymptotic result and the numerical result actually decreases as the wave becomes pinned. The reason for this is discussed in the next section where higher order corrections to the front velocity are derived. Higher order approximation to the front position To get the higher order approximation we have to do a few more calculations. We note that U0 and ∂U0∂ξ may be computed explicitly for the cubic case (see Appendix B.3): U0(ξ) = u+ + u− 2 − u+ − u− 2 tanh ( (u+ − u−)ξ 2 √ 2 ) (4.73) and ∂U0 ∂ξ = − 1 4 √ 2 (u+ − u−)2 1 cosh2 ( (u+−u−)ξ 2 √ 2 ) . (4.74) We now compute relations (4.63) and (4.65). Note in this case that: p∗(ξ) = p(−ξ) = p(ξ) = ∂U0 ∂ξ (ξ). (4.75) 88 4.4. Example with Cubic Kinetic Function 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 x ε =0.05 u(0) v(0) u(10) v(10) (a) 0 2 4 6 8 10 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 po sit io n time Position of front for ε =0.05 Actual front φ0 (b) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x ε =0.1 u(0) v(0) u(10) v(10) (c) 0 2 4 6 8 10 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 po sit io n time Position of front for ε =0.1 Actual front φ0 (d) Figure 4.3: Left: Initial conditions and steady state solutions to system 4.10 with the cubic kinetics function (4.12) are shown for various values of ². The evolution of the steepest point φ through time and the solution of Eq. 4.70 is shown on the right. We see that for small ² there is excellent agreement between the analytic solution and the zeroth order approximation. 89 4.4. Example with Cubic Kinetic Function 0 2 4 6 8 10 −1 −0.5 0 0.5 1 1.5 2 2.5 x 10 −3 ε =0.05 φ num −φ0 (a) 10−3 10−2 10−1 10−4 10−3 10−2 ε Er ro r φ num −φ0 (b) Figure 4.4: (a) The error over time between the numerical front position φnum obtained from numerically solving the full PDE system (4.10) with reaction kinetics (4.12) and the solution of the zero order asymptotic order approximation (Eq. 4.70). ² = 0.05, D = 1. (b) The maximum error between the front position φ and the ODE (4.70) for a variety of ² values. The error is order ². Using fu(u+, v0) = −v0(1 + v0), fv(u+, v0) = v0(1 + v0) (4.76) equation (4.53a) becomes G+ = − 1 v0(1 + v0) dv0 dt + V1(t). (4.77) We see that G− = 0 since fv(u−, v0) = 0. Note that this results in H1 = 0. Therefore, (4.63) becomes: dφ1 dt = (∫∞ −∞ ∂U0 ∂t ∂U0 ∂ξ dξ + ∫ u+ u− (fv(U0, V0)V1 +H2fu(U0, V0)U0 −H2f(U0, V0))dU0 ) (∫∞ −∞ ( ∂U0 ∂ξ )2 dξ ) , (4.78) where H2 = − 1 v0(1 + v0)2 dv0 dt + V1(t) (1 + v0) . (4.79) 90 4.4. Example with Cubic Kinetic Function Consider the integral in the denominator of (4.78): ∫ ∞ −∞ ( ∂U0 ∂ξ )2 dξ = √ 2 16 (u+ − u−)3 ∫ ∞ −∞ 1 cosh4 η dη = √ 2 12 (u+ − u−)3 = √ 2 12 (1 + v0)3. (4.80) Now we need to evaluate the numerator of (4.78). We have ∂U0 ∂t = 1 2 d dt (u+ + u−)− 12 d dt (u+ − u−) tanh ( (u+ − u−)ξ 2 √ 2 ) − ξ 4 √ 2 (u+ − u−) d dt (u+ − u−) 1 cosh2 ( (u+−u−)ξ 2 √ 2 ) . (4.81) The last two functions are odd functions in ξ and since ∂U0∂ξ is an even function, we see that the first term in the numerator of (4.78) reduces to: ∫ ∞ −∞ ∂U0 ∂t ∂U0 ∂ξ dξ = 1 2 d dt (u+ + u−) ∫ ∞ −∞ ∂U0 ∂ξ dξ = u− − u+ 2 d dt (u+ + u−) = −12(1 + v0) dv0 dt . (4.82) The second integral in the numerator of (4.78) yields: V1(t) ( (1 + v0)3 3 − (1 + v0) 2 2 ) +H2 ( −(1 + v0) 4 2 + (2 + v0)(1 + v0)3 3 ) = V1(t) 6 (1 + v0)2(2v0 − 1) + H26 (1 + v0) 3(1− v0). (4.83) Substituting in (4.79) into (4.83) obtain V1(t) 6 (1 + v0)2(2v0 − 1) + ( − 1 v0(1 + v0) dv0 dt + V1(t) ) 1 6 (1 + v0)2(1− v0) = V1(t) 6 v0(1 + v0)2 − 16v0 (1 + v0)(1− v0) dv0 dt . (4.84) Putting everything together, we find the following expression for (4.63). dφ1 dt = √ 2v0 1 + v0 V1(t) + √ 2(1− 2v0) v0(1 + v0)2 dv0 dt . (4.85) 91 4.4. Example with Cubic Kinetic Function Equation (4.85) depends on the zeroth order approximation to the inactive form v0, as well as the first order approximation at the interface, V1. We now determine V1 by computing (4.65). ∫ 1 0 v1dx = V1(t)− 23A+φ 3 0 − 2 3 A−(1− φ0)3 = V1(t)− 13D ( 1 + du+ dv ) dv0 dt φ30 − 1 3D ( 1 + du− dv ) dv0 dt (1− φ0)3 = V1(t)− 13D dv0 dt (2φ30 + (1− φ0)3). (4.86) where A− and A+ were given in (4.41). Likewise, ∫ 1 0 u1dx = ∫ φ0 0 u1dx = G+φ0 − 13A+φ 3 0 = ( − 1 v0(1 + v0) dv0 dt + V1(t) ) φ0 − 23D dv0 dt φ30 (4.87) We therefore have the following algebraic relation between φ1 and V1(t): − ( φ0 v0(1 + v0) + 1 3D (4φ30 + (1− φ0)3) ) dv0 dt + (1 + φ0)V1(t) + φ1(1 + v0) = 0 (4.88) We now have (4.67), (4.68), (4.85) and (4.88) as the equations that describe the solution to order ² in terms of φ0, v0, φ1 and V1. We can eliminate V1(t) and v0 from these equations differential equations for φ0 and φ1. As before, equations (4.67) and (4.68) become: dφ0 dt = K − 1− 2φ0√ 2(1 + φ0) , (4.89) v0 = K − φ0 1 + φ0 . (4.90) Note also that: dv0 dt = −(K + 1)(K − 1− 2φ0)√ 2(1 + φ0)3 . (4.91) 92 4.4. Example with Cubic Kinetic Function We now simplify (4.88). We have: ( (1 + φ0)2φ0 (K + 1)(K − φ0) + 1 3D (4φ30 + (1− φ0)3) ) dv0 dt + (1 + φ0)V1 + K + 1 1 + φ0 φ1 = 0. (4.92) We may therefore solve for V1 to obtain: V1 = − K + 1(1 + φ0)2φ1 − ( (1 + φ0)φ0 (K + 1)(K − φ0) + 1(4φ30 + (1− φ0)3) 3D(1 + φ0) ) dv0 dt . (4.93) In (4.85), we may substitute the expression for v0 in terms of φ0 to get: dφ1 dt = √ 2(K − φ0) K + 1 V1 + √ 2(1 + φ0)2(3φ0 − 2K + 1) (K − φ0)(K + 1)2 dv0 dt . (4.94) Note that when dv0dt = 0, i.e., v0 = 1, equation (4.94) reduces to dφ1 dt = − 1√ 2 ( φ1 1 +K ) , (4.95) which decays to 0. This explains why the error for the zeroth order approximation to the front location, φ0, actually becomes O(²2) as the wave becomes pinned. In general, before the wave becomes pinned we have dφ1 dt =− E(φ0)φ1 + F (φ0), where E(φ0) = √ 2(K − φ0) (1 + φ0)2 , F (φ0) = ( (1 + φ0)φ0 (K + 1)2 + 1(4φ30 + (1− φ0)3)(K − φ0) 3D(K + 1)(1 + φ0) − (1 + φ0) 2(3φ0 − 2K + 1) (K − φ0)(K + 1)2 ) (K + 1)(K − 1− 2φ0) (1 + φ0)3 . (4.96) Numerical verification To check the above theory against numerical experiment, we need to find appropriate initial conditions for the ODEs (4.89) and (4.96). Suppose at t = t0, the solution to the PDE system has relaxed so that the solution assumes the form predicted in the 93 4.5. Discussion asymptotic calculations above. We will use the numerical solution to determine φ0(t0) and φ1(t0), and use them as initial conditions for the ODEs.. We determine the position of the front, which we call φnum to be the point where the numerical solution u is equal to the unstable steady state u = 1. Let: φ0(t0) + ²φ1(t0) ≈ φnum, i.e., φ1 ≈ φ num − φ0(t0) ² (4.97) We need an additional equation to determine φ0(t0) and φ1(t0). We use the numerical value of the total concentration for v and equate it to an analogous expression arising from the asymptotic calculations. Let vnumtotal be the numerically computed total amount of v. Noting that: ∫ 1 0 vdx = v0 + ² ∫ 1 0 v1dx+O(²2) (4.98) we let v0 + ² ∫ 1 0 v1dx ≈ vnumtotal. (4.99) By substituting (4.90) for v0 and (4.86), (4.91) and (4.93) into ∫ v1dx, this gives us an algebraic relationship between φ0 and φ1. We then solve equations (4.97) and (4.99) to find φ0(t0) and φ1(t0). The solutions to (4.89) and (4.96) as a function of time are plotted in Fig. 4.5. Note that as the system goes to steady state, the difference between the two approximations disappears, indicating that at steady state the accuracy of the pinned front position (4.72) obtained from the leading order approximation is O(²2). The maximum error as a function of ² is plotted in Figure 4.6. As ²→ 0 the contribution of the first order correction term becomes marginal. 4.5 Discussion In this Chapter, I applied singular perturbation techniques to the wave-pinning sys- tem in 1 spatial dimension. I derived an algebraic differential equation system for the speed of the transition layer. Similar approaches were previously used to find the speed 94 4.5. Discussion 0 2 4 6 8 10 −1.5 −1 −0.5 0 0.5 1 1.5 2 x 10 −4 Er ro r time ε =0.005 φ num −φ0 φ num −φ0−εφ 1 (a) 0 2 4 6 8 10 −2 −1 0 1 2 3 4 5 6 x 10 −4 Er ro r time ε =0.01 φ num −φ0 φ num −φ0−εφ 1 (b) 0 2 4 6 8 10 −2 0 2 4 6 8 10 x 10 −4 Er ro r time ε =0.02 φ num −φ0 φ num −φ0−εφ 1 (c) 0 2 4 6 8 10 −5 0 5 10 15 20 x 10 −4 Er ro r time ε =0.04 φ num −φ0 φ num −φ0−εφ 1 (d) Figure 4.5: The error over time between the numerical front position φnum and the zero order asymptotic order approximation (Eq. 4.89) (dashed line), and the error between φnum and the two term asymptotic approximation (Eq. 4.89 and 4.96) (solid line). φnum is obtained from numerically solving the full PDE system (4.10). The two term approximation is consistently closer to the true solution than the zero order term. Note that with time the difference between the two approximations disappears. 95 4.5. Discussion 10−3 10−2 10−1 10−5 10−4 10−3 10−2 ε Er ro r φ num −φ0 φ num −φ0−εφ 1 (a) Figure 4.6: The maximum error between the numerical front position φnum, the one term approximation (4.89) and the two term asymptotic approximation (4.89 and 4.96) for a range of ² values. of the travelling pulse in the Fitzhugh Nagumo equations (46), and the speed of the travelling wave for the bistable equation in an inhomogeneous medium (22). An expla- nation for the wavepinning behaviour (based on the Maxwell condition) was presented. The pinning occurs due to depletion of the fast diffusing variable v. As illustrated in Figure 4.1, the wave propagates until v is depleted to the critical level vc, at which point a stable stationary solution emerges. This contrasts with the typical pattern one observes in a bistable system with diffusion, where the travelling front velocity given by (4.34a) depends on the parameters and generally does not vanish. In those cases, the travelling front can be pinned at specific parameter values, but zero velocity fronts are generally structurally unstable, and do not persist if parameters in the model are per- turbed slightly. Front pinning in reaction diffusion systems can also occur due to spatial inhomogeneities in the domain (57; 95). In our coupled bistable system wave-pinning can occur within a range of reaction parameter values on a homogeneous domain. This makes the behaviour much more robust, which as discussed in earlier chapters, is a feature that is important for biological systems. I have shown how the location of the pinned front, given by (4.34), depends on the mean concentration of material and the reaction kinetics, and derived a condition, given by (4.39) that describes the mean concentration of chemical necessary for pinning 96 4.5. Discussion to take place. A second order approximation was used to determine the error of the leading order approximation. The zeroth order approximation was found to agree well with the full numerical solution, with the error approaching O(²2) as the wave became pinned. 97 Chapter 5 Wave-pinning in 2D: Numerics and Asymptotic Analysis 5.1 Introduction In this Chapter I consider the wave-pinning model in two spatial dimensions. Chapters 2-4 presented numerics and analysis of the model in 1D, where we interpreted the do- main as a cross-section across the cell. However, we have to verify that the results will still hold on a more complicated geometry. Biologically, the problem in two dimensions corresponds to studying the distribution that Rho GTPases would attain in a flat cell, neglecting variations in cell thickness. The analysis presented here is done for an arbi- trary shaped cell, as seen from the top down, mathematically represented as a domain Ω ⊂ R2. For simplicity, all numerics were done on a rectangular domain. Many of our observations for a rectangular domain carry over to irregular domains. Numerical simulations show that starting from arbitrary initial data, internal layers are established on a very fast (O(1/²)) timescale, separating the regions in which the active form is at the high steady state, u ≈ u+(v), from the regions at which it is at the low steady state u ≈ u−(v) . On a timescale of tens of seconds (T0 = L√ηDu ≈ 30 s), wave pinning occurs, similar to the one-dimensional case. After these interfaces are pinned, the fronts slowly evolve such that the length of the interface between the high and low steady states is minimized. This phenomenon is known as flow by mean curvature. To complement the numerics, I use matched asymptotic expansions to determine the speed of the wavefront on the T0 timescale, and find that to zeroth order, the results 98 5.2. Asymptotic Analysis in 2D agree with those in one spatial dimension from Chapter 4. A correction to higher order, however, shows that the curvature of the front also affects the front speed. On a longer timescale, I show how the wave-pinning model relates to the Allen-Cahn equation with a mass constraint, a model for phase separation in a binary mixture, which is also known to undergo flow by mean curvature (99). 5.2 Asymptotic Analysis in 2D 5.2.1 Model Equations Consider the wave-pinning model ∂u ∂t = Du∆u+ ηf(u, v), (5.1a) ∂v ∂t = Dv∆v − ηf(u, v), (5.1b) on a region in the plane (Ω ⊂ R2), whose diameter has typical cell size L. Once again, we denote by u(x, t) and v(x, t) the concentrations of the active (bound to the cell membrane) and inactive forms (GDP-bound in complex with GDI’s in the cytosol) of the protein, respectively, at position x = (x1, x2) ∈ Ω and time t. These concentrations are, again, defined as number of molecules per unit “area”. Here Ω represents the two- dimensional top-down view of the cell and ∂Ω is the outer boundary of that 2D region. We impose no flux (Neumann) boundary conditions on this edge, as done previously in the 1D case. We denote the area of the region Ω by the notation |Ω|. The assumptions on the kinetic function f(u, v) are as for the 1D case. As before, equations (5.1) can be rescaled, exploiting the assumed unequal rates of diffusion, Du ¿ Dv, so that the scaled concentrations of active (u(x1, x2, t)) and 99 5.2. Asymptotic Analysis in 2D inactive (v(x1, x2, t)) protein satisfy ² ∂u ∂t = ²2∆u+ f(u, v), (5.2a) ² ∂v ∂t = D∆v − f(u, v), (5.2b) where x = (x1, x2) ∈ Ω is any point in the region, ²¿ 1, D = O(1) and ∆ is given by ∆ = ∂2 ∂x12 + ∂2 ∂x22 . The corresponding boundary conditions are: ∂u ∂n = ∂v ∂n = 0 on ∂Ω. (5.2c) Conservation also continues to hold in this 2D variant: ∫ Ω (u+ v)dx = K. (5.2d) On the timescale τ = t/ε, initial conditions evolve into a moving internal layer located on a curve Γ(t) in Ω (see Figure 5.1). Γ can be considered as a level curve for u(x, t) = um, where um is the unstable steady state. It will be convenient to track the movement of Γ(t) by considering the evolution of the function r(x, t) defined to be the signed distance from x to Γ (i.e., Γ(t) = {r(x, t) = 0}), defined to be positive for u(x, t) > um and negative for u(x, t) < um. By no-flux boundary conditions, Γ is orthogonal to ∂Ω. We assume that Γ divides the plane into two parts: Ω+ where r > 0 and Ω− where r < 0. We will construct an outer solution to (5.2) on Ω± and an inner solution near Γ. More formally, the inner solution is valid on the set Γ²(t) = {r(x, t) < ²} and the outer solution on its complement Ω² = {|r(x, t)| > ²}. As ²→ 0, Ω² → Ω− ∪ Ω+. 100 5.2. Asymptotic Analysis in 2D 5.2.2 Local Orthogonal Coordinate System In the vicinity of Γ, consider a local orthogonal coordinate system (r(x, t), s(x, t)), where s is arclength along the curve Γ, and r is the signed distance from Γ as before. Because Γ is a level curve (so that u and v are constant along Γ), this coordinate system essentially reduces the problem to a quasi-1D setting, where we examine what happens in front of and just behind the interface. On Γ we have |∇r| = 1 and ∆r = κ, where κ is the curvature of Γ, considered positive if Γ is concave as seen from Ω−. We assume that κ = O(1). Define the normal velocity at a point along the curve Γ(t) by the function c(s, t) = −rt. Transformation to local coordinates (s, r) leads to the following relations for the Laplacian and the time derivative: ∆u = ∂2u ∂r2 +∆r ∂u ∂r + ∂u∆s ∂s + ∂2u|∇s|2 ∂s2 = ²−2 ∂2u ∂ρ2 + ²−1κ ∂u ∂ρ + ∂u∆s ∂s + ∂2u|∇s|2 ∂s2 , (5.3a) ∂u ∂t = ∂u ∂t + rt ∂u ∂r + ∂s ∂t ∂u ∂s = ∂u ∂t + ²−1rt ∂u ∂ρ + ∂s ∂t ∂u ∂s , (5.3b) where we used the fact that ∆r = κ. r s G u_(v) u+(v) Figure 5.1: For a flat cell (top-down view), region of high Rho GTPase activity is separated from a region of low activity by some interface. That curve Γ is the boundary between the high and low steady states (u+(v) and u−(v)). Near Γ we consider a local coordinate system (r(x, t), s(x, t)), where r is the signed distance from Γ and s is the arclength along Γ. 101 5.2. Asymptotic Analysis in 2D 5.2.3 Outer Solution Consider equations (5.2) in the outer region, that is, away from the front Γ(t). Assume that u has the expansion u = u0 + ²u1 + ²2u2 · · · and likewise for v. Substituting these expansions into (5.2) and keeping only terms of order ²0 leads to the following equations for u0 and v0: 0 = f(u0, v0), (5.4a) 0 = D∆v0 − f(u0, v0). (5.4b) At first order, i.e. keeping terms of order ² leads to: ∂u0 ∂t = fu(u0, v0)u1 + fv(u0, v0)v1, (5.5a) ∂v0 ∂t = D∆v1 − fu(u0, v0)u1 − fv(u0, v0)v1. (5.5b) Equations (5.4) and (5.5) are valid in the regions Ω±. Adding (5.4a,b), we find: D∆v0 = 0. (5.6) Using no-flux boundary conditions (5.2c) we conclude that: v0(x, t) = v−0 on Ω−, v+0 on Ω+. (5.7) where v−0 and v + 0 do not depend on x. (We will later show that v − 0 = v + 0 .) From (5.4a) and assumptions on f(u, v), u0 takes on one of the values u+, u− or um in the outer regions. Proceeding as in the one dimensional case, we choose the outer solution given by u0(x, t) = u−(v−) on Ω−, u+(v+) on Ω+. (5.8) 102 5.2. Asymptotic Analysis in 2D 5.2.4 Inner Solution Now consider the inner region, close to the moving front at Γ(t). Γ is defined to be the location where the inner solution U = um(V ) to all orders. Introduce a stretched variable: ρ(x, t) = r(x, t) ² . (5.9) Let u(x, t) = U(ρ, s, t) and v(x, t) = V (ρ, s, t) in the inner region. Substitute this into (5.2a,b) to obtain the inner equation for U and V ² ( ∂U ∂t − ²−1c∂U ∂ρ + st ∂U ∂s ) = ²2 ( ²−2 ∂2U ∂ρ2 + ²−1κ ∂U ∂ρ + ∂U∆s ∂s + ∂2U |∇s|2 ∂s2 ) + f(U, V ), (5.10a) ² ( ∂V ∂t − ²−1c∂V ∂ρ + st ∂V ∂s ) = D ( ²−2 ∂2V ∂ρ2 + ²−1κ ∂V ∂ρ + ∂V∆s ∂s + ∂2V |∇s|2 ∂s2 ) − f(U, V ), (5.10b) where we have used the fact that −rt = c(V ) is the normal component of the front velocity. We again assume usual expansions for U , V and the front speed c = c0 + ²c1 + · · · in powers of ². We also expand f(U, V ) = f(U0, V0) + ²fU (U0, V0)U1 + ²fV (U0, V0)V1 +O(²2). (5.11) To leading order, we find: ∂2U0 ∂ρ2 + c0 ∂U0 ∂ρ + f(U0, V0) = 0, (5.12a) ∂2V0 ∂ρ2 = 0. (5.12b) 103 5.2. Asymptotic Analysis in 2D 5.2.5 Matching Inner and Outer Solutions Away from the interface, as ρ → ±∞ the matching conditions for un and Un are determined by matching terms of order ²n: U0(±∞, s, t) = u0(0±, s, t), (5.13a) U1(±∞, s, t) = u1(0±, s, t) + ρ∂u0 ∂r (0±, s, t), (5.13b) where 0+ means limit from the right, and 0− limit from the left. Similar conditions hold for V0 and V1. Note that if the outer solution un is not constant, expansion in terms of the in- ner variable will produce terms of higher order in ², which must be matched to the corresponding term in the inner expansion. For example lim r→0 u0(r) = u0(0) + r ∂u0(0) ∂r + · · · = u0(0) + ²ρ∂u0(0) ∂r + · · · . As the outer solutions u0 and v0 do not depend on r, it follows that (5.13b) simplifies to U1(±∞, s, t) = u1(0±, s, t), (5.14) V1(±∞, s, t) = v1(0±, s, t). (5.15) From (5.12b), it follows that V0 = a1(t)ρ+ a2(t). (5.16) We must match this solution with the values of v obtained above in the outer regions, i.e., lim ρ→−∞V0(ρ) = v − 0 , limρ→∞V0(ρ) = v + 0 . (5.17) In order for the above limits to exist, it follows that V0 is constant in the inner layer. 104 5.2. Asymptotic Analysis in 2D We conclude that v−0 = v + 0 = V0. (5.18) and thus, V0 is spatially uniform throughout the domain. We drop the dependence of v0 on ρ (or V0 on ξ). Hence, to first order, the inner solution is identical to the inner solution for the one-dimensional case. We now construct a solution for U0 in the inner layer. Since V0 is spatially constant in the inner layer, equation (5.12a) can be considered as an equation in U0, where V0 is a parameter (that varies in time). We must solve the boundary value problem (5.12a) with the following matching conditions as boundary conditions at ±∞: lim ξ→−∞ U0(ξ) = u−0 (V0), lim ξ→∞ U0(ξ) = u+0 (V0). (5.19) Let U0 = Ψ(ρ, V0) be the travelling wave solution to (5.12a) with boundary conditions (5.19). It can then be shown that, as in the 1D case, c0(V0) = ∫ u+(V0) u−(V0) f(u, V0) du∫∞ −∞ (∂Ψ(ρ, V0)/∂ρ) 2 dρ . (5.20) Note that by assumption 4 on f(u, v) (see Chapter 4), c0 changes sign at v0 = vc and c0 = 0 at v0 = vc. It can also be shown, as in 1D, that v0 = K |Ω| − ∫ Ω u0 dx |Ω| (5.21) for the constant K, which corresponds to the total amount of u and v. The integral of u0 can be evaluated by dividing it into contributions from the two outer regions (Ω±) and the inner region. ∫ |Ω| u0 dx = ∫ Ω− u0 dx+ ∫ Ω+ u0 dx+O(²) = u−(v0)|Ω−|+ u+(v0)|Ω+|+O(²) (5.22) 105 5.2. Asymptotic Analysis in 2D where |Ω±| denote the areas of the corresponding regions and where we have used (5.8) in the second equality. Hence to first order, |Ω|v0 = K − u+(v0)|Ω+| − u−(v0)(|Ω| − |Ω+|. (5.23) Differentiating (5.23) with respect to t and regrouping terms, obtain ( |Ω|+ du+(v0) dv0 |Ω+|+ du−(v0) dv0 (|Ω| − |Ω+|) ) dv0 dt = −(u+(v0)−u−(v0))d|Ω+| dt . (5.24) Note that the the normal velocity of the front position Γ, c = −rt, is related to the rate of change of the area behind the front d|D+|/dt through the relation d|Ω+| dt = ∫ Γ rt · ∇r|∇r| dr. For v0 = vc, c0(vc) = −rt = 0, so the front Γ will become stationary. As, in the 1D case, using the assumptions on f(u, v) it can be shown that 1 + du± dv > 0. (5.25) Using the fact that 0 < |Ω+| < |Ω| and (5.25), we see that the factor multiplying dv0/dt in (5.24) is positive. Since (u+ − u−) > 0, we conclude from (5.24) that dv0/dt and d|Ω+|/dt have opposite signs. So as the wavefront propagates (i.e., |Ω+|, the area behind the front increases) v0 is depleted until v = vc and the front position becomes stationary. The opposite signs of dv0/dt and d|Ω+|/dt means that this “locked” front is stable. 5.2.6 u and v to Order ² So far, the analysis has been fully analogous to the 1D case. However, once the front has evolved so that c0 = 0, higher order terms in the expansion c = c0+²c1+· · · become dominant. Here, we determine corrections to the front velocity c of order ², 106 5.2. Asymptotic Analysis in 2D We collect terms of order ² in (5.10a), obtaining: ∂2U1 ∂ρ2 + c0 ∂U1 ∂ρ +fU (U0, V0)U1 = (−c1−κ)∂U0 ∂ρ −fV (U0, V0)V1+ ∂U0 ∂t + ∂s ∂t ∂U0 ∂s . (5.26) Similarly, collecting terms of order ² in (5.10b) leads to ∂2V1 ∂ρ2 + κ ∂V0 ∂ρ = 0, (5.27) which can be simplified to ∂2V1 ∂ρ2 = 0, (5.28) as V0 is independent of ρ. From this, we conclude that V1(ρ) has the form: V1(ρ) = A(t)ρ+B(t) (5.29) We determine the constants by matching V1 to the outer solution v1. As in the V0 equation (5.12b), for this matching to be possible, V1 must be constant in ρ. Adding (5.5a,b) we obtain ∆v1 = 1 D ∂u0 ∂t + ∂v0 ∂t on Ω, (5.30) with no flux boundary conditions. The solution is given by v1(x, t) = ∫ Ω− 1 D ( ∂u−(v0) ∂t + ∂v0 ∂t ) G− dx0, on Ω−∫ Ω+ 1 D ( ∂u+(v0) ∂t + ∂v0 ∂t ) G+ dx0, on Ω+, (5.31) 107 5.2. Asymptotic Analysis in 2D where G±(x;x0) are the Green’s functions defined by ∆G± = 1 |Ω±| − δ(x− x0) on Ω±, (5.32a) ∂G± ∂n = 0 on ∂Ω ∩ ∂Ω± (5.32b) G± = V1 on ∂Ω± ∪ ∂Ω (5.32c)∫ Ω± G± dx = 0. (5.32d) The form of the Green’s function G depends on the domains Ω+, Ω−, and Ω. ∂ denotes the boundary of the domain. From equation (5.5a) it follows that u1(x, t) = ∂u0 ∂t − fv(u0, v0)v1 fu(u0, v0) . (5.33) Substituting (5.31) for v1 into (5.2.6) will lead to an expression for u1. 5.2.7 Front Speed to Order ² Next we will use a solvability condition on (5.26) to determine c1. (See Appendix B.4 for background information on solvability conditions.) Consider the operator D̂ ≡ ∂ ∂t + ( ∂s ∂t ) ∂ ∂s . It can be shown by elementary differential geometry (see Chapter 1 in (22) for a deriva- tion) that D̂F (ρ, x, t) = ∂F ∂t + c ∂F ∂r , (5.34) where c = −rt is the normal velocity of the front as before. Note that the last two terms on the RHS of (5.26) corresponds to D̂U0. As ρ does not depend on t it follows that D̂U0(ρ, V0) = ∂U0 ∂V0 ∂V0 ∂t + c ∂U0 ∂V0 ∂V0 ∂r = ∂U0 ∂V0 ∂V0 ∂t , as ∂V0 ∂r = 0. (5.35) 108 5.2. Asymptotic Analysis in 2D Hence, (5.26) can be rewritten as LU1 ≡ ∂ 2U1 ∂ρ2 +c0 ∂U1 ∂ρ +fU (U0, V0)U1 = (−c1−κ)∂U0 ∂ρ −fV (U0, V0)V1+ ∂U0 ∂V0 ∂V0 ∂t , (5.36) where U0(ρ, V0) is the travelling wave solution to (5.12a) as before. Following exactly the same solvability procedure on (5.36) as in Section 4.3.2 of Chapter 4, we obtain (−c1 − κ)α(t) + β(t) = 0, (5.37) where α = ∫ ∞ −∞ ∂U0 ∂ρ p∗dρ, (5.38a) β = ( ∂U0 ∂t − fv(U0, V0)V1 − fu(U0, V0)Ũ0 +H2f(U0, V0) ) p∗dρ, (5.38b) and p∗ is given by p∗ = ∂U0 ∂ρ (−ρ). (5.39) Note that α and β are identical to the LHS and RHS terms of equation (4.63) derived in Chapter 4 to describe the evolution of the first order correction to front velocity. As in Chapter 4, Ũ0 is given by: Ũ0 = G+ +G− 2 + G+ −G− u+(v0)− u−(v0) ( U0 − u+ + u−2 ) = H1 +H2U0, (5.40) where G+ = 1 fu(u+(v0), v0) ( du+(v0) dt − fv(u+(v0), v0)V1(t) ) , (5.41a) G− = 1 fu(u−(v0), v0) ( du−(v0) dt − fv(u−(v0), v0)V1(t) ) . (5.41b) 109 5.3. Numerical Results in 2D It follows that, to order ², the front speed is given by c = c0(V0)− ²κ+ ²β(t) α(t) . (5.42) The first term in the expansion of c is identical to the speed of the front determined in 1D to order one. The term ²κ represents slow movement of the front due to curvature12 of the interface Γ. The form of the second term of order ² is identical to the O(²) correction to the wave-speed in the one-dimensional case derived in Section 4.3.4. Note that once the system has evolved to c0(V0) = 0, movement of the front due to curvature becomes dominant. This is discussed in the next section. In the case where f is a function of u alone, (5.42) reduces to the result of (100) for a single bistable equation. 5.3 Numerical Results in 2D x y t=2.0 t=0.5 t=0 1 1 (a) x y t=0 t=0.5 t=0.5 t=0 t=2.0 t=1.0 1 1 (b) x y t=0 t=15 t=15 t=10 t=5 t=0 t=10 t=5 t=20 1 1 (c) Figure 5.2: Solutions to (5.2) with cubic reaction term (5.43) on a rectangular domain for various initial conditions. Level curves corresponding to the position of the front (i.e. the curve Γ) are shown. ² = 0.05. (a) The wave front propagates outward at a decreasing speed, eventually pinning. (b) The cell is initialized with two Gaussian peaks of high u value. The peaks merge and the front then evolves to minimize curvature. (c) The cell is initialized with two step-like peaks of high u value. A transient two- front state is seen, but one of the peaks later collapses, resulting in a single stable peak solution. To go beyond the regime in which asymptotic results are valid, and to check those results where they apply, I here discuss numerical simulations on a rectangular domain Ω. The model consisting of equations (5.2) was simulated numerically on a 2D grid in 12Note that there is an implicit assumption in the asymptotics that κ = O(1). 110 5.3. Numerical Results in 2D 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x y t=0.05 t=2.0 t=0.4 (a) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (b) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (c) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 2 (d) Figure 5.3: (a) Solutions to (5.2) with a Hill function reaction term (5.44) on a rectan- gular domain in response to a graded stimulus in k0. Noisy initial conditions are used. Level curves corresponding to the position of the front (i.e. the curve Γ) are shown at the indicated times. ² = 0.05. (b-d) Evolution of system (5.2) with a Hill function reaction term (5.44) in response to a reversal of stimulus. (b) Initially a gradient in k0 results in a high steady state at x = 0. (c) When the stimulus is reversed the wave-front travels in the opposite direction. (d)A new high steady state is reestablished at x = 1. 111 5.3. Numerical Results in 2D Matlab using an IMEX (implicit-explicit) method13. The system was investigated for reaction terms: f(u, v) = u(1− u)(u− 1− v), (5.43) and for the more biologically realistic f(u, v) = v ( k0 + γu2 m2 + u2 ) − δu, (5.44) for various initial conditions. As in Chapter 2, we obtain the following biologically relevant results that extend to 2D that are shown in Figure 5.3: 1. The homogeneous steady state distribution v, u = u±(v) is stable to small pertur- bations in either variable14. 2. If the perturbation is above some threshold (which could itself be relatively small), a wave response is triggered. Perturbations we tested included a transient and/or localized input to the equation (5.2a) for u or some suitable choice of non-uniform initial conditions that are close (but not too close) to the stable steady state. 3. Imposing a graded initial condition in u or in the activation parameter k0 in (5.44), e.g. a linear ramp across the domain, results in wave front propagation, slowing down, and pinning of the wave provided conditions previously derived (on parameters and total amount of material) are satisfied. (See Figure 5.3 a) 4. The non-homogeneous distribution can be reversed by applying a (suitably large) input, whether graded or step-like, at the opposite side of the domain. (See Figure 5.3 b-d) The biological relevance of the above points have already been discussed in Chapter 2. 13Crank-Nicholson method is used for the diffusion operator and a second order Runge-Kutta scheme is used for the nonlinear reaction term. See (101) for details about implicit-explicit methods. 14As stated in Chapter 3, this contrasts with many other mechanisms based on Turing diffusion- driven pattern formation, where the homogeneous steady state is unstable to small amplitude spatially heterogeneous perturbations. 112 5.4. Similarity to the Constrained Allen Cahn Model For small values of ², results of numerical simulations are shown using level curves in Figure 5.2. For values of ² larger than some threshold, whose value depends on the total amount of material, no pinning occurs and the perturbation resolves into a homogeneous steady state solution (not shown). The evolution of the front in response to multiple stimuli is shown in Figure 5.2b,c. Here, we observe either merging of multiple peaks (Figure 5.2b) or a metastable multiple peak pattern followed by a sudden collapse of one of the peaks, resulting in a stable solution with a single peak (Figure 5.2c). The evolution of the front in response to noisy initial conditions is shown in Figure 5.4 and Figure 5.3 for different types of reaction terms. For small values of epsilon (² ≈ 0.01 for |ω| = 1) multiple peak patterns are quickly established (Figure 5.4a,b and Figure 5.3a). We observe that on a longer timescale these peaks merge to minimize the curvature of the wave front (see Figure 5.4and Figure 5.3). Note that the length of the interface Γ is also being minimized. The exact timescale on which the peaks merge depends on the type of reaction kinetics that are used. 5.4 Similarity to the Constrained Allen Cahn Model This flow by mean curvature behaviour seen in the numerical simulations is not sur- prising if one considers the results for some related systems in the literature. As stated previously, a single RD equation of the type ∂u ∂t = Du ∂2u ∂x2 + f(u; v), (5.45) with v as a parameter, has fronts that move along its normal with constant velocity determined by (5.20). Here, the reaction kinetics are assumed to satisfy f(u; v) =W ′(u), where W (u) is the double well-potential and v is a parameter, However, it is known for equation (5.45) that if the Maxwell condition is satisfied (i.e., the numerator in (5.20) is zero and the travelling wave in 1D is stationary), the front’s normal velocity is ²κ, where κ is the front’s mean curvature (3; 100). A front will then either shrink to a point 113 5.4. Similarity to the Constrained Allen Cahn Model 0 0.5 1 0 0.5 1 0.2 0.4 0.6 0.8 1 (a) 0 0.5 1 0 0.5 1 0.2 0.4 0.6 0.8 1 (b) 0 0.5 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 (c) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (d) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (e) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (f) Figure 5.4: Evolution of system (5.2) with a Hill function reaction term (5.44) on a rectangular domain from noisy initial conditions in 2D. Snapshots of u are shown at t = 0.01, t = 0.2, t = 0.4, t = 0.6, t = 0.8, and t = 1.0. As time progresses multiple peaks coalesce in a coarsening process. ² = 0.01. 114 5.4. Similarity to the Constrained Allen Cahn Model 0 0.5 1 0 0.5 1 0.4 0.6 0.8 1 (a) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (b) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (c) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (d) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (e) 0 0.5 1 0 0.5 1 0 0.5 1 1.5 (f) Figure 5.5: Evolution of system (5.2) with the cubic reaction term (5.43) on a rectan- gular domain from noisy initial conditions in 2D. Snapshots of u are shown at t = 0.1, t = 2, t = 4, t = 6, t = 8, and t = 10. As time progresses multiple peaks coalesce in a coarsening process. ² = 0.01. 115 5.4. Similarity to the Constrained Allen Cahn Model in finite time, or tend to a to a locally shortest diameter of the domain Ω. For example, in a dumbbell-shaped region, Ω, curve shortening by mean curvature may lead to a stationary front that is orthogonal to the neck of the domain. In general, for a single reaction diffusion equation (5.45), a nonconstant stable steady state can result only if the domain Ω is convex (66). In our case, the Maxwell condition is satisfied when the interfaces become pinned. Mass conservation due to (5.2d) prevents the interfaces from shrinking to a point. Once the interfaces are pinned on the t timescale, on the timescale T = ²t our system is closely related to a mass-constrained Allen-Cahn equation (99), consisting of one bistable reac- tion diffusion equation with a global integral constraint, used to model phase separation: ² ∂u ∂t = ²2∆u+ f(u)− 1|Ω| ∫ Ω f(u) dx, (5.46) where f(u) = W ′(u) for bistable potential W . Because the total amount of u in the domain is conserved: ∫ Ω u(x, t) dx =M, system (5.46) support structurally stable zero velocity front solutions, that are also observed in the wavepinning system. Our wave-pinning analysis above rested upon the reduction of our reaction diffusion system (5.1) to a single equation where the equation for the inactive form was reduced to an integral constraint. For this reason, the two models show qualitatively similar behavior (coarsening) on a long time scale. In the constrained Allen Cahn, equation (5.46), the normal velocity, c, of the inter- face Γ also exhibits flow by mean curvature on the O(²) timescale: c ∼ ²(κ− 1|Γ| ∫ Γ κ ds), (5.47) where |Γ| is the length of the interface. A single, closed interface will evolve to a circle enclosing the same area. Whereas, in the unconstrained Allen Cahn equation (5.45), this circle will shrink to a point, for the mass conserved system, the circle will exponentially 116 5.5. Discussion slowly drift towards the closest point on the boundary ∂Ω (117). When multiple interfaces are present in the constrained Allen Cahn system, large areas of high u will grow while smaller areas of high u shrink, preserving the total area enclosed by the interfaces. We see the same type of coarsening behaviour in the wavepinning system (5.2). However, the wave-pinning model has conservation of the total amount of u and v, rather than the total area enclosed by the interfaces, resulting in some differences between the two models. For example, unlike the wavepinning model, the globally constrained Allen Cahn model does not support stable homogeneous solutions. 5.5 Discussion The study of phase separation of two or more media has been explored mathematically for decades, (3; 10; 25) and is an ongoing area of interest in applied mathematics. Ap- plications have primarily included physico-chemical systems, but here we have extended methods used in this field to address a problem from cell biology. Here we extended the asymptotic analysis done in Chapter 4 to a 2D domain. Not surprisingly, we found that to leading order, the asymptotic approximation to the 2D solution behaves as it did in the 1D case close to the front. An integral condition analogous to the Maxwell condition is obtained. The curvature of the front locally affects the wave speed only in the terms of order ², as demonstrated by Eq. (5.42). Analogous results have been shown for wavefronts in the Fitzhugh-Nagumo and Belousov-Zhabotinskii models (23). Our numerical results verify that dynamics tend to minimize the curvature. However, we believe that flow by mean curvature we observe is too slow to play a major role in cell polarization. Domain shape can also introduce curvature effects. Simulations by Marée (63) suggest that in two dimensions this effect of curvature minimization might be important for accelerating polarization in cells with dynamic, rather than static, shapes. 117 Chapter 6 Bifurcation Analysis of the Wave-pinning System 6.1 Introduction In this chapter we consider the steady-state behaviour of the wave-pinning model when the parameter ² is no longer small, that is, when the diffusion of the inactive form of the protein is not too different from that of the active form. Biochemically, this can be achieved by interfering with membrane-cytosol cycling that is regulated by RhoGDI. As a result, the proportions of fast and slow diffusing subpopulations of inactive forms will change, affecting the effective diffusion coefficient of the inactive form. To motivate this study, we simulated the wave-pinning behaviour of the 1D time- dependent system (4.10) with cubic reaction kinetics to steady state, increased the value of the parameter ² by a small increment, and continued the simulation repeatedly. For small ², this resulted in a gradual change in the wave shape and stall position. Beyond some critical ² value, the stationary wave disappears and is no longer sustained. Results of this simulation are shown in Figure 6.1 (a-b). Such qualitative change in the nature of the solution is called a bifurcation. Figure 6.1c plots the difference between the max- imum and the minimum of the solution, and suggests that there is an unstable branch of solutions joining the two stable branches (stationary wave solutions and spatially homogeneous solutions) shown in the diagram. In this regime, the asymptotic analysis described in the previous sections is no longer valid, as the assumption ² ¿ 1 no longer holds. Hence, we study the transition 118 6.2. Model Equations and Stability of the Homogeneous Steady States from a spatially heterogeneous to a spatially homogeneous solution, that occurs as the parameter ² is increased, using numerical continuation methods. Some analytical insights can also be obtained from phase plane analysis. In this chapter, we use both approaches to study the bifurcation structure of the wave-pinning system. 6.2 Model Equations and Stability of the Homogeneous Steady States As before, the wave-pinning system is given by ² ∂u ∂t = ²2 ∂2u ∂x2 + f(u, v), (6.1a) ² ∂v ∂t = D ∂2v ∂x2 − f(u, v). (6.1b) The system we consider is the steady state version of (6.1), namely 0 = ²2uxx + f(u, v), (6.2a) 0 = Dvxx − f(u, v). (6.2b) As before, the boundary conditions are no flux: ∂u ∂x ∣∣∣∣ x=0 = ∂u ∂x ∣∣∣∣ x=1 = 0, ∂v ∂x ∣∣∣∣ x=0 = ∂v ∂x ∣∣∣∣ x=1 = 0. (6.2c) Recall that in this dimensionless system, the length of the domain has been scaled to 1. All the numerics in this chapter were carried out with cubic reaction kinetics f(u, v) = u(m− u)(u−m− v), (6.2d) though analysis only requires that the reaction terms be bistable. The homogeneous steady states of this system are u− = 0, um = m and u+(v) = m + v. (Note that the dimensionless model discussed in earlier chapters corresponds to the case m = 1.) It 119 6.2. Model Equations and Stability of the Homogeneous Steady States sp ac e time ε =0.145 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 (a) sp ac e time ε =0.146 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 (b) 0.13 0.135 0.14 0.145 0.15 0 0.5 1 1.5 ε m a x− m in Bifurcation diagram (c) Figure 6.1: Change in the nature of solutions to (6.1) as ² is increased. (a,b) Space-time plots of the solution to (6.1) with cubic reaction (6.2d). Parameters used: m = 0.8, K = 1.49, D = 1. (a) For ²=0.145 a wavefront solution is established. (b) For ²=0.146, with the same initial condition, the wave-front solution collapses, and a homogeneous solution emerges. (c) Bifurcation diagram obtained by solving the time-dependent system (6.1) with cubic reaction kinetics given by (6.2d) for ²i to steady state and then using this steady state solution as an initial condition for the PDE for ²i+1. This approach allows us to observe the stable solution branches only. The difference between the maximum and minimum of the solution is plotted. 120 6.3. Numerical Continuation can easily be shown that, for a given v, u = u± are saddle points of system (6.2a), while um is a center. It is shown in Appendix B.2 that u = u± are linearly stable homogeneous solutions of the PDE system (6.1), while the homogeneous steady state u = um is an unstable steady state of the PDE system. That is, u = u± are stable to perturbations of the form ũ ṽ = α1 α2 cos kx eσt. (6.3) Note that these types of perturbations satisfy the no-flux boundary conditions and preserve the mass constraint, provided k = npi/L = npi (as domain length L = 1). We’re restricted to considering perturbations in the u, v that preserve that mass constraint in order that (6.1) be satisfied. 6.3 Numerical Continuation To study the transition from a spatially heterogeneous branch of solutions to a spatially homogeneous branch of (6.1) as ² is varied, I used numerical continuation to follow the spatially heterogeneous branch of solution as ² is increased. Several numerical software packages for bifurcation analysis of dynamical systems are publicly available. Since the software is written for ODE systems, PDE system are typically discretized using the method of lines. We tested two of the most successful packages: AUTO, developed by E.J. Doedel, and MATCONT, continuation software in Matlab developed by Y.A. Kuznetsov. Neither of these packages was able to successfully compute solution branches for system (6.1). The reason for this failure is that the numerical continuation algorithms require the ODE system to be non-singular (i.e., the matrix corresponding to the discretized version should be invertible) at non-bifurcation points. It is easily seen that system (6.1) does not satisfy that requirement. Thus, another condition, such as the conservation condition (2.10), must be appended to the system. However, appending the conservation condition to system (6.1) created numerical instability in 121 6.3. Numerical Continuation the resulting ODE system, so we were unable to use both AUTO and MATCONT. We also reformulated system (6.2) as a boundary value problem15, but were unable to use the Matlab BVP solver, bvp4c, to converge to the solution. For this reason, we designed our own numerical continuation code, implementing a pseudo-arclength continuation algorithm (105), details of which are described below. 6.3.1 Pseudo-Arclength Continuation The main idea of numerical continuation is to compute a sequence of points (ui, ²i), i = 1, 2, · · · which approximate the solution curve of R(u; ²) = 0 to some specified tolerance. This typically involves a predictor-corrector algorithm. Once a good initial guess is determined, Newton’s method is commonly used as a corrector. In zero order continuation for some accepted solution (ui; ²i) the next point on the solution curve is determined by setting ²i+1 = ²i +∆², and (ui+1) is determined using (ui) as a predic- tor. This technique gives no indication of the sensitivity of solution to the bifurcation parameter and will fail near a fold point, where the branch of solutions turns around. An improvement to this is first order continuation, which uses the tangent vector of the solution with respect to the bifurcation parameter in the initial guess. R(u; ²) = 0 =⇒ J(u; ²)∂u ∂² + ∂R ∂² = 0. Solving for ∂u∂² , and setting the initial guess for ui+1 to be ui + ∆²i ∂u ∂² alleviates this problem. Note that at the fold point the Jacobian J(u; ²) is singular, so Newton’s method will fail there. An extra scalar condition has to be appended to the system R(u; ²) = 0. In pseudo-arclength continuation, system R(u; ²) is augmented with an equation for the 15System (6.2) can be turned into a boundary value problem (BVP) by adding two additional ODES ²′ = 0 and M ′ = u + v subject to the boundary conditions M(0) = 0, M(1) = K. Together with system (6.2) this results in 6 first order equations and 6 boundary conditions. 122 6.3. Numerical Continuation n=0 (ui,ei) (ui+1,ei+1) (ui+Deif,ei+Dei) R(u,e)=0 Figure 6.2: An illustration of the supplemental condition (6.4) in pesudo-arclength continuation. The hyperplane given by equation (6.4) is perpendicular to the tangent vector (u̇i, ²̇i). The intersection of this hyperplane with the bifurcation curve will be non zero provided the curvature of the branch R(u, ²) = 0 and the distance ∆s between the hyperplane and (ui, ²i) are not too large. plane, whose normal is the tangent vector to the bifurcation curve: n ≡ u̇i(u− ui) + ²̇i(²− ²i)−∆s = 0, (6.4) where (u̇i, ²̇i) is the tangent vector at (ui, ²i), and ∆s is the distance of the plane from (ui, ²i) (overdot represents differentiation with respect to arclength s). (See Figure 6.2 for an illustration.) The Jacobian of this new augmented system is J(u; ²) ∂R∂² u̇i ²̇i, which is non-singular at fold points. For an accepted point (ui, ²i) we iterate to (ui+1, ²i+1) using the following algorithm:un+1 ²n+1 = un ²n + ∆un ∆²n , (6.5) 123 6.3. Numerical Continuation where Jn ∂Rn∂² u̇i ²̇i ∆un ∆²n = − Rn u̇i(u− ui) + ²̇i(²− ²i)−∆s . (6.6) Now, u̇i and ²̇i can be determined from the following considerations: R(u; ²) = 0 =⇒ Ju̇i + ∂R ∂² ²̇i = 0, and ‖u̇i‖2 + ‖²̇i‖2 = 1. (6.7) Choose u̇i = a ∂u ∂² , and ²̇i = a, where J ∂u ∂² = −∂R ∂² and a = ± 1√ 1 + ‖∂u∂² ‖2 . The sign of a is chosen to preserve the direction in which the branch is traversed. 6.3.2 Numerical Continuation on the Full System We spatially discretize system (6.2a-6.2b) with the 3 point approximation for the Lapla- cian L (modified to account for no-flux BCs), which results in a system of 2N equa- tions with a spatial step of ∆x = 1/N . The spatially discretized system is given by R(u,v) = 0, where R(u,v) = ²2Lui + f(ui, vi), DLvi − f(ui, vi) (6.8) for i = 1 · · ·N . Note, however, that system (6.8) is singular, so it needs to be augmented with the discretized version of the conservation condition ∆x i=N∑ i=1 (ui + vi) = K. (6.9) 124 6.4. Reduction of the Steady State System to an Integral System System (6.8) is modified by replacing u1 in the first equation of (6.8) with u1 = K ∆x − i=N∑ i=2 vi − i=N∑ i=2 ui, (6.10) which is obtained from the conservation condition. This modified system was then solved using pseudoarclength continuation. Several ways of appending the conservation condition (6.9) were tested and gave the same result. We observed a fold (saddle-node) bifurcation at the expected location (compare Figure 6.3 to Figure 6.1c), but were only able to continue along a part of the unstable branch before the continuation converged to the homogeneous branch of the solution. The reasons for this are discussed later in the chapter. 0.1425 0.143 0.1435 0.144 0.1445 0.145 0.1455 0.146 1.395 1.4 1.405 1.41 1.415 1.42 1.425 1.43 1.435 1.44 ε Um ax −U m in Bifurcation diagram (a) 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 1.405 1.41 1.415 1.42 1.425 1.43 1.435 1.44 1.445 D Um ax −U m in Bifurcation diagram (b) Figure 6.3: Bifurcation diagrams for system (6.8) using pseudo-arclength continuation. Compare to Figure 6.1c. (a) Bifurcation diagram with respect to ². D = 1,m = 0.8,K = 1.49. (b) Bifurcation diagram with respect to D. ² = 0.14, m = 0.8, K = 1.49. 6.4 Reduction of the Steady State System to an Integral System By eliminating v we can rewrite system (6.2) as an ODE system for u and ux, where x is a “time-like” variable. This allows us to use phase plane methods (30) to study the nature of the stationary solutions to the PDE system in question. This technique is well- 125 6.4. Reduction of the Steady State System to an Integral System known for the bistable case (30), but leads to a somewhat more complicated analysis for our system. Here, we will only be concerned with the existance of (monotonic) front solutions and homogeneous steady state solutions. u- u+ um u0 u1 (a) -4 -3 -2 -1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 u- u+ u ux u0 u1um (b) (c) −0.5 0 0.5 1 1.5 2 2.5 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 u F( u,A ,B ) (d) -8 -6 -4 -2 0 2 4 6 8 0 0.5 1 1.5 2 2.5 ux u u- um u+ (e) (f) Figure 6.4: The solutions of (6.13) viewed in the u-ux phase-plane and the correspond- ing F (u,A,B) given by (6.19). Equations (6.1) were solved numerically to steady state. Values of A and B were computed from (6.12) and (6.17), and the corresponding po- tential function F (u,A,B) and the u-ux phase-plane are plotted. Top: System (6.2) is in the wave-pinning regime for ² = 0.13, D = 1, m = 0.8. (a) F (u,A,B) has four zeros. (b) The phase plane corresponding to the potential in (a) given by Eq. (6.13). The spatially heterogeneous solution is the orbit going from u0 to u1. This trajectory satisfies both (6.21) and (6.24). (c) The trajectory from u0 to u1 viewed as a function of x. Bottom: When ² is increased to ² = 0.147, system (6.2) is no longer in the wave- pinning regime. (d) The corresponding potential function F (u,A,B) is recalculated for the new values of A and B. F (u,A,B) has the two lowest zeros coalesce into a double root, which is a local minimum. (e) The phase plane corresponding to the potential in (b). There exists no trajectory satisfying (6.21) and (6.24). (f) The steady state solution to (6.1) is the spatially homogeneous solution u− = 0. We reduce our system by the following steps: Add (6.2a,b) to obtain ²2uxx +Dvxx = 0 (6.11) 126 6.4. Reduction of the Steady State System to an Integral System Integrate this equation twice and use the boundary conditions (6.2c) to obtain ²2u+Dv = A. (6.12) Substituting in v = A/D − (²2/D)u into (6.2a), we rewrite system (6.2) as a single equation in u and A: 0 = ²2uxx + f̃(u,A), (6.13) where f̃(u,A) = f(u, v(A)). By mass conservation, the constant of integration, A, is given by A = DK − (D − ²2) ∫ 1 0 u dx. (6.14) Equation (6.13) can be written as a system of two first order ODEs: ux = w, (6.15a) wx = −f̃(u,A). (6.15b) This means that we can study solutions in the u-w phase plane (also denoted u-ux phase plane), and construct trajectories such as those in Figure 6.4b,e. Note that A is not an arbitrary constant. Rather, it is determined by the total mass (fixed by the initial conditions) and the nonlocal term ∫ 1 0 u dx. Given the same total mass, different steady states of u have different A values, and the phase plane corresponding to (6.13) will shift for each steady state solution. Hence, the appropriate value of A can only be determined once we know the full solution. Suppose, for example, that the value of A corresponding to some steady state solu- 127 6.4. Reduction of the Steady State System to an Integral System tion is known. Multiplying both sides of (6.13) by ux we obtain 0 = d dx ( 1 2 ²2u2x + F̃ (u,A) ) , where F̃ (u,A) = ∫ f̃(u,A) du. (6.16) Integrating with respect to u and using the no-flux boundary conditions we obtain ²2u2x + F̃ (u,A) = B, (6.17) where B is an integration constant. Let F (u,A,B) = B − F̃ (u,A). This function is analogous to a potential energy whose level curves in the uux plane are the solution curves we desire. For a given steady state solution, we can plot this potential (i,e., obtain Figure 6.4). By BCs (6.2c) we consider only trajectories with u(0) ≡ u0, ux(1) = 0; and u(1) ≡ u1, ux(1) = 0. (6.18) As we are only interested in the monotonic fronts, such trajectories are “half-loops” that start and end on the u axis. Multiple loops correspond to multiple peak solutions. Note that as ux = 0 at x = 0, 1, it follows from (6.17) that B = F̃ (u0, A) = F̃ (u1, A). Hence, we can evaluate B once we determine the value of u at either endpoint of the domain. Since f̃(u,A) is bistable, F (u,A,B) has the form of a double well potential (see Figure 6.4). In the specific case of the cubic reaction term (6.2d), we compute that F (u,A,B) = B + αu4 + βu3 + γu2 = B + α 2 u4 − 2 3 (mα+ γ(A))u3 +mγ(A)u2, (6.19) where α = 1 + ²2 D , γ(A) = m+ A D (6.20) Not all curves described by (6.17) that start and end on the ux = 0 axis are feasible solutions of (6.2). First, we must satisfy the constraint that the domain length is 1. 128 6.4. Reduction of the Steady State System to an Integral System This can be written in the form 1 = ∫ 1 0 dx = ∫ u1 u0 du |ux| = ² ∫ u1 u0 du√ F (u,A,B) , (6.21) where u0,1 are the values of u at the ends of the domain x = 0, 1. We have used, from (6.17) that |ux| = √ F (u,A,B)/², (6.22) and the fact that the solution curve u is a monotone function of x. A second constraint that is imposed by our model is conservation of total protein, ∫ 1 0 (u+ v) dx = K. . This leads to the second integral equation that F (u,A,B) must satisfy, namely: K = ∫ 1 0 (u+ v)dx = ∫ u1 u0 ( u+ A D − ² 2 D u ) du |ux| = ( 1− ² 2 D ) ² ∫ u1 u0 u du√ F (u,A,B) + A D . (6.23) Using (6.12) this second constraint can be rewritten as K −A/D 1− ²2/D = ² ∫ u1 u0 u du√ F (u,A,B) . (6.24) Once A is given, we can reduce system (6.2) to a conservative system in the u− ux phase plane. Curves of constant F (u,A,B) that start and end on the u axis, and satisfy constraints (6.21) and (6.24), correspond to solutions of system (6.2). However, different steady states have distinct A and B values, and have different versions of the phase plane diagram shown in Figure 6.4. From the numerics, we observe that at a homogeneous steady state of (6.1) there is a coalescence of two roots of the potential function F (u,A,B) (see Figure 6.4d). When the system has a spatially heterogeneous solution, F (u,A,B) has four zeros. At the high u, low v steady state, the two highest zeros merge into a double root at the local 129 6.4. Reduction of the Steady State System to an Integral System minimum u+. At the low u, high v steady state, the two lowest zeros become a double root at the local minimum u−. When system (6.2) is at the unstable homogenous steady state um, the two intermediate zeros of F (u,A,B), u0 and u1, coalesce into a double root at the local maximum um. We believe that this coalescence of the zeros of F (u,A,B) near the homogeneous steady states is the reason that our numerical continuation on the discretized system described in Section 6.3.2 was unable to completely capture the unstable branch. As ² is decreased, u0 moves closer and closer to the fixed point u−, so the resulting trajectory spends most of its “time” near u− (i.e., the front solution is close to 0 for most of the domain). When u0 gets below a certain threshold, Newton’s method for the pseudo- arclength continuation converges to the homogeneous solution u−. In summary, the dynamical system and the two integral constraints to be satisfied contains two unknown constants, A and B. We have reduced the problem to finding values of A and B such that there exists a heterogeneous trajectory satisfying ²2u2x = F (u,A,B), (6.25a) T (u0) ≡ ² ∫ u1 u0 du√ F (u,A,B) = 1, (6.25b) M(u0) ≡ ( 1− ² 2 D ) ² ∫ u1 u0 u du√ F (u,A,B) + A D = K. (6.25c) The constraint (6.25b) is analogous to the “time-map” outlined in (30), i.e. this integral is similar to “time” (length of domain) needed to transverse an orbit segment from u0 to u1. Smoller and Wasserman (106) showed that for cubic reaction kinetics, the “time- map” T (u0) for the scalar Neumann problem uxx+f(u) = 0 is monotone, so there exists at most one nonconstant solution. Furthermore, this nonconstant solution is always unstable on convex domains (66). Our system is more complicated by the presence of a second integral constraint that comes from mass conservation. This complication allows the emergence of a bifurcation structure that is absent in the case of a single bistable RD equation. 130 6.4. Reduction of the Steady State System to an Integral System Since each monotonic heterogeneous steady state solution of system (6.1) is char- acterized by a unique (A,B) pair, we have reduced the system (6.2) to a system of two integral relations for two unknowns A and B. These integrals are singular at the endpoints u0 and u1, so a transformation has to be made to eliminate the singularity for numerical evaluation (see Appendix B.6). Note that the homogeneous steady states cannot be described by a time-map, since at these steady states the time map integrals become undefined. 6.4.1 Numerical Continuation on the Time Map System We also used pseudoarclength continuation with an adaptive step in arclength ds on a reduced ‘time-map’ system of two integral equations: R(A,B; ²) = ² ∫ u1 u0 du√ F (u,A,B) − 1 (1− ²2D ) ² ∫ u1 u0 u du√ F (u,A,B) +A/D −K. (6.26) This means that numerical continuation can be done on a system of two integral equa- tions rather than on a spatial discretization of system (6.2), which, with the 3 point approximation for the Laplacian, results in a system of 200 equations when using a spatial step of ∆x = 0.01. We obtained an initial numerical approximation for the integration constants A and B by solving the time-dependent problem (6.1) to steady state and then evaluating expressions (6.12) and (6.17) (to avoid numerically approximating the derivative, it is best to evaluate B at x = 0 or 1). Using these estimates for A and B we validate expressions (6.25b) and (6.25c) numerically, obtaining values very close to 1 for (6.25b) and K for (6.25c) (the total amount is fixed by the initial conditions for u and v). For a given initial value of ²0, an initial guess A0 and B0 are determined by solving the full time dependent system to steady state, and calculating A0 and B0 as described 131 6.4. Reduction of the Steady State System to an Integral System above. Newton’s method An+1 Bn+1 = An Bn − R(An, Bn; ²)J(An, Bn; ²) (6.27) is then used to refine (A,B) to some specified tolerance. Pseudoarclength continuation was then carried out with this initial guess. For Newton’s method, we have to calculate the Jacobian: J = ∂ ∂A ( ² ∫ u+ u− du√ F (u,A,B) − 1 ) ∂ ∂B ( ² ∫ u+ u− du√ F (u,A,B) − 1 ) ∂ ∂A ( (1− ²2D ) ² ∫ u+ u− u du√ F (u,A,B) + AD −K ) ∂ ∂B ( (1− ²2D ) ² ∫ u+ u− u du√ F (u,A,B) + AD −K ) . (6.28) We numerically estimate the Jacobian by using the approximation ∂F ∂A ≈ F (A+ h,B)− F (A,B) h (6.29) and similar approximations for the other partial derivatives including ∂R/∂². For details on integral calculation see Appendix B.6. The resulting solution curves are shown in Fig. 6.5. The bifurcation diagram agrees with the results obtained from solving the full time-dependent system to steady state (Fig. 6.1c). We also obtain comparable results to numerical continuation on a spatial discretization of system (6.2) with the 3 point approximation for the Laplacian (see Figure 6.3). Again, we were only able to compute a partial bifurcation diagram. This occurs because for small ², the variable B ≈ 0 making it difficult to find the endpoint of the integral u0. This issue arises because at B = 0 the two roots around u− = 0 coalesce, and it becomes impossible to accurately determine u0 (see Figure 6.4). Ways of rescaling the integral are discussed in Appendix B.6. 132 6.4. Reduction of the Steady State System to an Integral System 0.13 0.132 0.134 0.136 0.138 0.14 0.142 0.144 0.146 0.85 0.9 0.95 1 1.05 1.1 ε A (a) 0.13 0.132 0.134 0.136 0.138 0.14 0.142 0.144 0.146 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 x 10−4 B (b) 0.13 0.132 0.134 0.136 0.138 0.14 0.142 0.144 0.146 1.36 1.38 1.4 1.42 1.44 1.46 1.48 1.5 Um ax −U m in ε (c) Figure 6.5: Bifurcation diagram showing the spatially heterogeneous branch in the time map system (6.25) obtained using pseudo-arclength continuation. A fold point occurs at ²=0.145. Compare to Figure 6.1c. (a) A = ²2u + Dv is plotted. Parameters used: D = 1, m = 0.8, A = 0.848, K = 1.49. (b)B = ²2u2x + ∫ f̃(u,A) du is plotted. (c) The difference between the endpoints of the time map integral (u1 − u0) is plotted. 133 6.4. Reduction of the Steady State System to an Integral System 6.4.2 Two Parameter Bifurcation Diagrams For a range of ² values we also computed bifurcation diagrams with respect to the dimensionless parameter D = Dv/ηL2 and the total amount of material in the domain, K. The resulting bifurcation diagrams are plotted in Figure 6.6. As expected, we observe that when ² is held fixed and D is decreased the wavepinning system transitions from a stable spatially heterogeneous solution to spatially homogeneous one in a saddle node bifurcation with an unstable spatially heterogeneous solution (see Figure 6.3b and Figure 6.6b). Treating K as a bifurcation parameter, we observe a 2 parameter cusp bifurcation (see Figure 6.6a). Note that there is a range of K values for which a spatially hetero- geneous solution exists. For ² ¿ 1 this range corresponds to m < Ktotal/L < 3m, as predicted from the asymptotic analysis. As ² is increased this range becomes smaller, with a stationary wavefront solution persisting only at the isolated value Ktotal/L = 2m for large ² (see Fig. 6.6). K ε 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 (a) D ε 0.02 0.04 0.06 0.08 0.1 0.12 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (b) Figure 6.6: Two parameter continuation plots showing the region in which wave-pinning occurs (red). Blue indicates a region with spatially homogeneous solutions only. (a) Plot of K versus ² obtained by zero-order continuation on system (6.8) with respect to ² for each value of K. The asymptotic prediction for occurrence of wavepinning is 0.8 < K < 2.4 (for ² ¿ 1). Other parameters used: D = 1, m = 0.8. (b) Plot of D versus ². Other parameters used: K = 1.6, m = 0.8. 134 6.4. Reduction of the Steady State System to an Integral System 0 0.1 0.2 0.3 0.4 0.5 0m 0.5m 1m 1.5m 2m 2.5m 3m ε v εcrit 0.2 0.215 0.225 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16 ε v K=1.95m 0.2 0.202 0.2035 0.9 0.92 0.94 0.96 0.98 1 1.02 ε v K=2.05m 1 2 K=3.0m K=2.0m K=1.5m K=2.0m K=1.5m stable stable unstable unstable K=2.5m unstable unstable K=2.0m Figure 6.7: The full bifurcation diagram of the spatially heterogeneous solutions in the time map system (6.34) in the case of D = ∞. F (u, v,B) is given by (6.32). ²crit =√ v pi2/m+m/D is the value at which the spatially heterogeneous branches of solutions merges with the unstable homogeneous steady state. Figure produced using code for integral evaluation courtesy of Yoichiro Mori. We show as inserts, two examples for 0.2 < ² < 0.22 of the stable branches corresponding to the heterogeneous solutions in the time map system (6.34) for D = ∞ obtained using pseudo-arclength continuation for K = 1.95m (v > 1) and K = 2.05m (v < 1). Note that v = A/D when D → ∞. Thus, these branches correspond to curves just above, at, or below v = 1m in the full bifurcation diagram 6.4.3 Full Bifurcation Diagram in the limit D →∞ Using (6.12) we can rewrite the second variable v as v = A D − ² 2 D u. (6.30) Note that as D →∞ (i.e., the v variable is assumed to be well-mixed) then v = A D . (6.31) 135 6.4. Reduction of the Steady State System to an Integral System We can then rewrite (6.19) as F (u, v,B) = B + 1 2 u4 − 2 3 (2m+ v)u3 +m (m+ v)u2, (6.32) and the integral constraints (6.25) as 1 = ² ∫ u1 u0 du√ F (u, v,B) , (6.33a) K = ² ∫ u1 u0 u du√ F (u, v,B) + v. (6.33b) System (6.33) can be rewritten as a single equation ∫ u1 u0 du√ F (u, v,B) − ∫ u1 u0 u du√ F (u,v,B) v = 0. (6.34) For a given K and v, (6.34) can be solved for B using the bisection method. From this value of B a corresponding ² value can then be determined using ² = 1∫ u1 u0 du√ F (u,v,B) . (6.35) (See Appendix B.6 for details of integral evaluation, which are non-trivial.) The result- ing diagram of spatially heterogeneous solutions are plotted in Figure 6.7. Note that as mentioned earlier, the time map approach only describes monotonically increasing or decreasing solutions. Hence, multiple peak solutions are not shown in this diagram. The branches emanating from the branch v = 1.0m, correspond to the stable stationary wave front solutions, and are seen to undergo a fold bifurcation as ² increases. Plotting K versus the ² values for which a fold point occurs (Figure 6.8), we obtain a figure similar to Figure 6.6a, which was attained with pseudo-arclength continuation. The other branches seen in the full diagram are all unstable spatially heterogeneous solu- tions that we could not obtain using the pseudo-arclength continuation methods. At 136 6.5. Discussion 0 0.05 0.1 0.15 0.2 0.25 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 K ε Figure 6.8: Two parameter bifurcation diagram in the time map system (6.34) in the case of D = ∞. Plot of K versus ² is obtained by plotting the location of the fold points shown in Figure 6.7 for each value of K. Wave-pinning regime occurs in the area between the curves. Red circle indicates the maximum ² value for which wave-pinning occurs (at K = 2m). ²crit = √ v pi2/m+m/D these spatially heterogeneous branches of solutions merges with the the unstable homogeneous steady state (see Appendix B.5 for details). How much of this global bifurcation structure is retained in the case D < ∞ is currently under investigation. 6.5 Discussion In this chapter we considered the bifurcation structure of the wave-pinning system. We focused on the case when ² is no longer small, that is, when the diffusion of the inactive form of the protein is not too different from the active form. This could result if the inactive forms spend a greater fraction of time on the membrane, i.e., in the case of defective GDI. With pseudo-arclength continuation methods, we were able to find a second, spatially heterogeneous frozen wavefront solution to the wave-pinning system. 137 6.5. Discussion This unstable spatially heterogeneous solution can be viewed as a threshold solution separating the stable heterogeneous and the homogeneous steady states. Hence the level of stimulus required to switch the wave-pinning system from a spatially homogeneous to a spatially heterogeneous solution regime increases with ². This has implications for possible values of the diffusion coefficient of the active form Du, the size of the domain L, and the rate of the reaction kinetics η under which wave-pinning is possible. As ² is increased, the threshold solution approaches the stable heterogeneous solution, and a saddle node bifurcation occurs at the ² value where a transition between the heterogeneous and the homogeneous steady state occurs in the time-dependent system. The bifurcation we observe is reminiscent of travelling wave blocking for a bistable equation on a spatially inhomogeneous domain. In that system when a travelling wave passes over a “gap” region with decreased excitability, a stable and an unstable station- ary wave solutions are created in a saddle node bifurcation at a critical gap size (57). Similar results hold for systems of bistable reaction-diffusion equations (36). However, we have demonstrated that pinning can occur on a spatially homogeneous domain if the conditions given in section 4.2 are satisfied. 138 Chapter 7 Conclusion In this thesis, I have explained how a single pair of proteins in active and inactive forms, cycling between plasma membrane and cytosol, together with conservation and cooper- ative positive feedback could suffice to provide a basic cellular polarization mechanism. Part of my analysis points to the functional importance of such membrane-cytosol exchange and its role in separating the rates of diffusion of the intermediates. Candidates for proteins that have such features are the Rho-GTPases with their GEF-mediated activation and GAP-mediated inactivation. Cytosolic Rho-GTPases, the inactive forms, diffuse up to two orders of magnitude faster than the active membrane-bound forms. Thus, the inactive forms inherently act as “global information carriers”. One of the contributions of this thesis is the proposal that not only well-known Tur- ing pattern forming mechanisms, but also wave-propagation based aspects of reaction- diffusion systems, specifically the freezing and pinning of those waves, can lead to cell polarity. A second contribution is the explicit formulation of the necessary qualitative properties that the reaction kinetics must satisfy for the wave-pinning mechanism to work. Unlike some models for cell polarity with global or long-ranged inhibition discussed in Chapter 3, our model does not require a global inhibitor. Instead, it is based on a depletion mechanism, in which simple interconversion between two species guarantees a conservation principle. Wave pinning depends on three properties: 1. Mass conservation, stemming from no-flux boundary conditions and simple inter- conversion of the two forms. 2. The fast diffusion of the inactive form v relative to u (i.e., Dv À Du). This means 139 Chapter 7. Conclusion that the level of v across the cell is spatially uniform, though not constant in time. Depletion of v forms the key reason for stalling the wave. Thus, overexpressing or changing the amount of Rho GTPases in the cell, should affect the cell’s ability to polarize. 3. The assumed bistability of the reaction term f(u, v). That is, there are two stable homogeneous states corresponding to a high level of the active form and a low level of the active form everywhere in the cell. Based on these assumptions I was able to provide an explanation for the pinning be- haviour based on the well-known Maxwell condition that is illustrated in Figure 4.1. In Chapter 4, using matched asymptotic expansions, I derived an ODE for the location of the moving wave front. I showed how the location of the pinned front depends on the total amount of material and the reaction kinetics, and derived a condition on the total amount of material necessary for pinning to take place. In Chapter 5, I looked at the wave-pinning model in 2 dimensions, and was able to relate its behaviour to the mass- constrained Allen-Cahn equation, a well-known model for phase separation. In Chapter 6 I considered the wavepinning system in a regime where the matched asymptotics are no longer valid. I was able to show that in addition to the spatially homogeneous so- lutions and the stationary wave front solution that we observed through numerics, the wavepinning system also supports several unstable spatially heterogeneous solutions. Our model captures a number of essential properties of cellular polarization that were described in Chapter 1. First, our “model cell” responds to a variety of stimuli (described in Chapter 2) by spatially amplifying the stimulus, and attaining a polarized state. This occurs for both localized as well as graded signals. Second, the polarization is maintained even after the signal is removed. Third, the cell remains sensitive to further stimuli, and reverses its polarity in response to a sufficiently strong second signal of opposite orientation. Finally, the model cell polarizes in response to spatially random (“noisy”) signals under some circumstances. These behaviours are observed in numerous cell types upon exposure to chemoattractant, from Dictyostelium discoideum 140 Chapter 7. Conclusion to neutrophils. Several types of experimental findings lend support to application of this model to Rho GTPases. Some biochemical evidence of a positive feedback loop in Cdc42 activation (59) supports the Hill function term in our kinetics. Other evidence points to phosphorylation of RhoGDI as a possible mechanism to control which proportion of Rho GTPase is in the active form (16; 17). Experimental evidence in yeast suggests that GTPase cycling is directly responsible for the polarization response, and abolishing it leads to loss of polarity (37; 118). A mutation that interferes with the cycling of Cdc42 between its GTP and GDP bound states of Cdc42 in yeast results in mutants with multiple polarization sites (12). Furthermore, other proteins signalling to Cdc42 such as Rsr1p (Bud1p) can exist in both membrane-bound and cytosolic fractions, and deletion of these proteins leads to a travelling patch of active Cdc42 rather than a stable cap formation (87). The mechanism described here may prove to have greater implications in cell biology than just Rho GTPases: any bistable kinetics, together with reduced mobility of one of two interchangeable states (e.g. due to binding to the cytoskeleton, membrane, or other cellular structures) could, in principle, lead to similar phenomena. Other interconvert- ible protein pairs include the GTPase Rap1. Disruption of Rap1’s GTP/GDP cycle in melanoma cells leads to abolishment of polarity (Mike Gold, personal communication.) It would be also interesting to alter the GTPase cycle in motile cells such as neutrophils and keratocytes and see its effect on polarity. The model suggests a number of experimental tests that could be used to check its validity. For example, we predict that incremental overexpression of the Rho protein in the cell would incrementally shift the position of the transition zone. Beyond some level, further overexpression should abolish cellular polarity altogether, as the total amount of protein in all forms no longer satisfies the condition required for wave-pinning inside the domain. The prediction that overexpression of a Rho GTPase leads to abolishment of polarity has already been observed with Cdc42 in yeast (4). 141 Chapter 7. Conclusion We have also focused more directly on the effect of changing the small parameter ², which stems from the difference in the diffusion of active membrane-bound and inactive cytosolic protein. We showed that a transition from spatially heterogeneous to a spa- tially uniform solution occurs as the diffusion coefficient of the active form approaches that of the inactive form. Assuming a typical cell diameter of 10 µm, reaction timescale η = 1 s−1, and diffusion coefficients Du = 0.1µm2s−1 and Dv = 10µm2s−1, the dimen- sionless constants are ² ≈ 0.03 and D ≈ 0.1, which is far away from the bifurcation regime. However, we should start to see polarization failure under the following kinds of scenarios: 1. Increasing the diffusion coefficient of the active, membrane-bound form to Du = 1µm2s−1 would bring the system into the bifurcation regime. 2. Decreasing the cell size to L ≈ 3µm would similarly bring the system into the bifurcation regime. Since small keratocyte cell fragments devoid of organelles are capable of polarization (115), Cell fragments could be made successively smaller to test the effect of domain size and determine if there is a critical size necessary for polarization in that experimental system. 3. Decreasing the rate of the interconversion of the two GTPase forms η to 0.1 s−1 would also make the bifurcation apparent. We predict that slowing down the membrane cytosol exchange would also lead to problems with maintaining cell polarization. In normal cells under usual conditions, rates of diffusion of proteins in the cell interior and on its membrane are indeed vastly different (by 100-fold). However, interestingly, the diffusion constants Du, Dv are not necessarily fixed values. Indeed, inactive Rho GTPases have to be extracted from the membrane into the cytosol by specialized pro- teins, i.e., GDIs. The abundance and affinity of GDIs (which could be genetically or biochemically controlled) affects the fraction of time that inactive protein v spends in cytosol versus membrane, and hence modulates the magnitude of its effective diffusion 142 Chapter 7. Conclusion coefficient. Our results suggest that GDI regulation could be an important mechanism for regulating polarity by affecting GTPase cytosol membrane-cycling and activity (16)). A defective GDI mutant, where the inactive proteins spend more time on the plasma membrane, would be particularly useful to validate our model. 143 Bibliography [1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter. Molecular Biology of the Cell. 5th ed. Garland Science, 2008. [2] E. Albrecht and H. R. Petty. Cellular memory: Neutrophil orientation reverses during temporally decreasing chemoattractant concentrations. PNAS, 95(9):5039– 5044, 1998. [3] S. Allen and J. Cahn. A Microscope Theory for Antiphase Boundary Motion and its Application to Antiphase Domain Coarsening. Acta Metall., 27(6):1085–1095, 1979. [4] S. Altschuler, S. Angenent, Y. Wang, and L. Wu. On the spontaneous emergence of cell polarity. Nature, 454:886–9, 2008. [5] N. Barkai and S. Leibler. Robustness in simple biochemical networks. Science, 387:913–917, JUNE 1997. [6] W. M. Bement, H. A. Benink, and G. von Dassow. A microtubule-dependent zone of active RhoA during cleavage plane specification. J Cell Biol., 170(1):91–101, JUL 2005. [7] C. D. Beta, G. Amselem, and E. Bodenschatz. A bistable mechanism for direc- tional sensing. New J. Phys., 10:083015, 2008. [8] A. L. Bishop and A. Hall. Rho GTPases and their effector proteins. Biochem. J., 348:241–255, June 2000. 144 Chapter 7. Bibliography [9] H. Bourne. Neutrophil polarity. http://www.cmpharm.ucsf.edu/bourne/lab_ science/Polarity.html. [10] J. Cahn and J. Hilliard. Free Energy of a Nonuniform System. I. Interfacial Free Energy. The Journal of Chemical Physics, 28:258, 1958. [11] A. E. Carlsson and D. Sept. Mathematical modeling of cell migration. Methods Cell Biol., 84:911–937, 2008. [12] J. P. Caviston, S. E. Tcheperegine, and E. Bi. Singularity in budding: A role for the evolutionarily conserved small GTPase Cdc42p. PNAS, 99(19):12185–12190, 2002. [13] L. P. Cramer, T. J. Mitchison, and J. A. Theriot. Actin-dependent motile forces and cell motility. Curr. Opin. Cell Biol., 6(1):82–86, Feb. 1994. [14] A. Dawes and L. Edelstein-Keshet. Phosphoinositides and rho proteins spatially regulate actin polymerization to initiate and maintain directed movement in a 1d model of a motile cell. Biophys J, 92:1–25, 2007. [15] A. de Candia, A. Gamba, F. Cavalli, A. Coniglio, S. D. Talia, F. Bussolino, and G. Serini. A simulation environment for directional sensing as a phase separation process. Sci. STKE, 2007(378):pl1 – pl1, 2007. [16] C. DerMardirossian, G. Rocklin, J.-Y. Seo, and G. M. Bokoch. Phosphorylation of RhoGDI by Src Regulates Rho GTPase Binding and Cytosol-Membrane Cycling. Mol. Biol. Cell, 17(11):4760–4768, 2006. [17] C. DerMardirossian, A. Schnelzer, and G. M. Bokoch. Phosphorylation of RhoGDI by Pak1 Mediates Dissociation of Rac GTPase. Molecular Cell, 15:117–127, 2004. [18] P. Devreotes and C. Janetopoulos. Eukaryotic chemotaxis: Distinctions between directional sensing and polarization. J. Biol. Chem., 278(23):20445–20448, JUN 2003. 145 Chapter 7. Bibliography [19] J. Elf and M. Ehrenberg. Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. Syst. Biol., 1(2):230–236, DEC 2004. [20] R. G. Endres and N. S. Wingreen. Precise adaptation in bacterial chemotaxis through assistance neighborhoods. Proc Natl Acad Sci U S A., 103(35):13040– 13044, AUG 2006. [21] J. E. Ferrell, Jr. Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability. Curr. Opin. Cell Biol., 14(2):140–148, Apr. 2002. [22] P. Fife. Mathematical aspects of reacting and diffusing systems, volume 28 of Lect. notes in biomath. Springer, New York, 1979. [23] P. Fife. Dynamics of internal layers and diffusive interfaces. SIAM, 1988. [24] A. Gamba, A. de Candia, S. Di Talia, A. Coniglio, F. Bussolino, and G. Serini. Diffusion-limited phase separation in eukaryotic chemotaxis. Proc. Natl. Acad. Sci. USA, 102(47):16927–16932, NOV 2005. [25] A. Gamba, I. Kolokolov, V. Lebedev, and G. Orentzi. Patch coalescence as a mechanism for eukaryotic directional sensing. Phys. Rev. Lett., 99(15):158101–1– 158101–4, 2007. [26] T. Gardner, C. Cantor, and J. Collins. Construction of a genetic toggle switch in Escherichia coli. Nature, 403(6767):339–42, 2000. [27] A. Gierer and M. H. A theory of biological pattern formation. Kybernetik, 12:30– 39, 1972. [28] E. Giniger. How do Rho family GTPases direct axon growth and guidance? a proposal relating signaling pathways to growth cone mechanics. Differentiation, 70(8):385–396, OCT 2002. 146 Chapter 7. Bibliography [29] A. B. Goryachev and A. V. Pokhilko. Dynamics of cdc42 network embodies a turing-type mechanism of yeast cell polarity. FEBS letters, 582(10):1437–1443, 2008. [30] P. Gridnrod. The theory and applications of reaction-diffusion equations: patterns and waves. 2nd ed. Oxford University Press, 1996. [31] A. Hall. Rho GTPases and the actin cytoskeleton. Science, 279:509–514, 1998. [32] M. Holmes. Introduction to Perturbation Methods. Springer, 1998. [33] P. A. Iglesias and P. Devreotes. Navigating through models of chemotaxis. Current Opinion in Cell Biology, 20(1):35 – 40, 2008. Cell structure and dynamics. [34] P. A. Iglesias and A. Levchenko. Modeling the cell’s guidance system. Sci. STKE, 2002(148):RE12, Sept. 2002. [35] M. Iijima, Y. E. Huang, and P. Devreotes. Temporal and spatial regulation of chemotaxis. Developmental Cell, 3(4):469–478, OCT 2002. [36] H. Ikeda and M. Mimura. Wave-blocking phenomenona in bistable reaction- diffusion systems. SIAM J. Appl. Math, 49(2):515–538, APR 1989. [37] J. Irazoqui, A. Gladfelter, and D. Lew. Scaffold-mediated symmetry breaking by cdc42p. Nature Cell Biol., 5(12):1062–70, DEC 2003. [38] R. E. Itoh, K. Kurokawa, Y. Ohba, H. Yoshizaki, N. Mochizuki, and M. Mat- suda. Activation of Rac and Cdc42 video imaged by fluorescent resonance energy transfer-based single-molecule probes in the membrane of living cells. Mol. Cell. Biol., 22(18):6582–6591, Sept. 2002. [39] C. Janetopoulos and F. R. A. Directional sensing during chemotaxis. FEBS Letters, 582:2075–2085, 2008. 147 Chapter 7. Bibliography [40] C. Janetopoulos, J. Borleis, F. Vazquez, M. Iijima, and P. Devreotes. Temporal and spatial regulation of phosphoinositide signaling mediates cytokinesis. Dev. Cell, 8(4):467–477, 2005. [41] C. Janetopoulos, T. Jin, and P. N. Devreotes. Receptor-mediated activation of heterotrimeric g-proteins in living cells. Science, 291:2408–2411, 2001. [42] C. Janetopoulos, L. Ma, P. N. Devreotes, and P. A. Iglesias. Chemoattractant- induced phosphatidylinositol 3,4,5-trisphosphate accumulation is spatially ampli- fied and adapts, independent of the actin cytoskeleton. Proc. Natl. Acad. Sci. USA, 101(24):8951–8956, JUN 2004. [43] A. Jilkine. Mathematical study of Rho GTPases in motile cells. Master’s thesis, University of British Columbia, 2005. [44] A. Jilkine, A. F. M. Marée, and L. Edelstein-Keshet. Mathematical model for spatial segregation of the Rho-family GTPases based on inhibitory crosstalk. Bull. Math. Biol., April 25, Epub, 2007. [45] J. Keener. Principles of Applied Mathematics. Perseus Books, 2000. [46] J. Keener and J. Sneyd. Mathematical Physiology. Springer, 1998. [47] D. A. Kessler and H. Levine. Stability of traveling waves in the Belousov- Zhabotinskii reaction. Phys. Rev. A, 41(10):5418–5430, May 1990. [48] A. Kimmel and C. Parent. The signal to move: D-discoideum go orienteering. Science, 300(5625):1525–1527, JUN 2003. [49] B. Knox, P. Devreotes, A. Goldbeter, and S. L.A. A molecular mechanism for sensory adaptation based on ligand-induced receptor modification. Proc Natl Acad Sci U S A., 83:2345–2349, 1986. [50] V. Kraynov, C. Chamberlain, G. Bokoch, M. Schwartz, S. Slabaugh, and 148 Chapter 7. Bibliography K. Hahn. Localized Rac activation dynamics visualized in living cells. Science, 290(5490):333–337, OCT 2000. [51] J. Krishnan and P. A. Iglesias. Receptor-Mediated and Intrinsic Polarization and Their Interaction in Chemotaxing Cells. Biophys. J., 92(3):816–830, 2007. [52] B. Kutscher, P. Devreotes, and P. Iglesias. Local excitation, global inhibition mechanism for gradient sensing: An interactive applet. Science’s STKE, 219:pI3, 2004. [53] M. Ladwein and K. Rottner. On the Rho’d: The regulation of membrane protru- sion by Rho-GTPases. FEBS Letters, 582:2066–2074. [54] Y. Lazebnik. Can a biologist fix a radio? Or what I learned while studying apopotosis. Cancer Cell, 2(179):179–182, 2002. [55] A. Levchenko and P. Iglesias. Models of eukaryotic gradient sensing: Application to chemotaxis of amoebae and neutrophils. Biophys. J., 82(1):50–63, JAN 2002. [56] H. Levine, D. A. Kessler, and W.-J. Rappel. Directional sensing in eukary- otic chemotaxis: A balanced inactivation model. Proc. Natl. Acad. Sci. USA, 103(26):9761–9766, 2006. [57] T. J. Lewis and J. P. Keener. Wave-block in excitable media due to regions of depressed excitability. SIAM J. Appl. Math, 61(1):293–316, 2000. [58] Z. Li, C. Aizenman, and H. Cline. Regulation of Rho GTPases by crosstalk and neuronal activity in vivo. Neuron, 33(5):741–750, FEB 2002. [59] Q. Lin, W. Yang, D. Baird, Q. Feng, and R. A. Cerione. Identification of a DOCK180-related Guanine Nucleotide Exchange Factor That Is Capable of Me- diating a Positive Feedback Activation of Cdc42. J. Biol. Chem., 281(46):35253– 35262, 2006. 149 Chapter 7. Bibliography [60] L. Ma, C. Janetopoulos, L. Yang, P. Devreotes, and P. Iglesias. Two complemen- tary, local excitation, global inhibition mechanisms acting in parallel can explain the chemoattractant-induced regulation of PI(3,4,5)P-3 response in Dictyostelium cells. Biophys. J., 87(6):3764–3774, DEC 2004. [61] W. Ma, A. Trusina, H. El-Samad, W. A. Lim, and C. Tang. Defining network topologies that can achieve biochemical adaptation. Cell, 138(4):760 – 773, 2009. [62] E. Marco, R. Wedlich-Soldner, R. Li, S. Altschuler, and W. L.F. Endocytosis optimizes the dynamic localization of membrane proteins that regulate cortical polarity. Cell, 129:411–22, 2007. [63] A. Marée, V. Grieneisen, and L. Edelstein-Keshet. Phosphoinositides mediate competing lamellipods: insights from a 2d computational model of cell motility. In preparation. [64] A. Marée, A. Jilkine, A. Dawes, V. Grieneisen, and L. Edelstein-Keshet. Polarisa- tion and movement of keratocytes: a multiscale modelling approach. Bull. Math. Biol., 68(5):1169–1211, 2006. [65] N. Markevich, M. Tsyganov, J. Hoek, and B. Kholodenko. Long-range signaling by phosphoprotein waves arising from bistability in protein kinase cascades. Mol Syst Biol., 2:61–69, 2006. [66] H. Matano. Asymptotic behavior and stability of solutions of semilinear diffusion equations. Publ. Res. Inst. Math. Sci., 15(2):401–454, 1979. [67] T. Matozaki, H. Nakanishi, and Y. Takai. Small G-protein networks - their crosstalk and signal cascades. Cellular Signalling, 12(8):515–524, AUG 2000. [68] R. Meili and R. Firtel. Two poles and a compass. Cell, 114(2):153–156, JUL 2003. [69] H. Meinhardt. Orientation of chemotactic cells and growth cones: models and mechanisms. J. Cell Sci., 112(17):2867–2874, SEP 1999. 150 Chapter 7. Bibliography [70] H. Meinhardt. The Algorithmic Beauty of Sea Shells, 3rd Ed. Springer, Heidelberg, New York, 2003. [71] H. Meinhardt and P. A. J. de Boer. Pattern formation in Escherichia coli: A model for the pole-to-pole oscillations of Min proteins and the localization of the division site. Proceedings of the National Academy of Sciences, 98(25):14202–14207, 2001. [72] H. Meinhardt and A. Gierer. Application of a theory of biological pattern forma- tion based on lateral inhibition. J. Cell Sci., 15:321–346, 1974. [73] B. A. Mello and Y. Tu. Perfect and near-perfect adaptation in a model of bacterial chemotaxis. Biophys J., 84(5):2943–2956, MAY 2003. [74] A. Mogilner. Mathematics of cell motility: have we got its number? J. Math. Biol., 58:105–134, 2009. [75] A. Mogilner and G. Oster. Cell motility driven by actin polymerization. Biophys. J., 71(6):3030–3045, DEC 1996. [76] Y. Mori, A. Jilkine, and L. Edelstein-Keshet. Wave-pinning and cell polarity from a bistable reaction-diffusion system. Biophysical Journal, 94(9):3684–97, 2008. [77] J. Murray. Mathematical Biology, Second Edition. Springer, 1993. [78] P. Nalbant, L. Hodgson, V. Kraynov, A. Toutchkine, and K. Hahn. Activation of endogenous Cdc42 visualized in living cells. Science, 305(5690):1615–1619, SEP 2004. [79] A. Narang. Spontaneous polarization in eukaryotic gradient sensing: A mathe- matical model based on mutual inhibition of frontness and backness pathways. J. Theor. Biol., 240(4):538–553, JUN 2006. [80] A. Narang, K. Subramanian, and D. Lauffenburger. A mathematical model for chemoattractant gradient sensing based on receptor-regulated membrane phos- 151 Chapter 7. Bibliography pholipid signaling dynamics. Annals Of Biomedical Engineering, 29(8):677–691, AUG 2001. [81] Y. Nishiura. Singular limit approach to stability and bifurcation for bistable reaction diffusion-systems. Rocky Mt. J. Math., 21(2):727–767, 1991. [82] C. Nobes and A. Hall. Rho, Rac and Cdc42 GTPases regulate the assembly of multimolecular focal complexes associated with actin stress fibers, lamellipodia, and filopodia. Cell, 81(1):53–62, APR 1995. [83] B. Olofsson. Rho guanine dissociation inhibitors: pivotal molecules in cellular signalling. Cell. Signal., 11(8):545–554, Aug. 1999. [84] M. Onsum and C. V. Rao. A mathematical model for neutrophil gradient sensing and polarization. PLoS Comput Biol., 3(3):e36, MAR 2007. [85] M. Onsum and C. V. Rao. Calling heads from tails: the role of mathematical modeling in understanding cell polarization. Curr. Opin. Cell Biol., 21(1):74–81, 2009. [86] M. Otsuji, S. Ishihara, C. Co, K. Kaibuchi, A. Mochizuki, and S. Kuroda. A mass conserved reaction-diffusion system captures properties of cell polarity. PLoS Comput Biol., 3(6):e108, 2007. [87] E. M. Ozbudak, A. Becskei, and van Oudenaarden A. A system of counteracting feedback loops regulates cdc42p activity during spontaneous cell polarization. Dev. Cell, 9(4):565–71, OCT 2005. [88] C. Parent and P. Devreotes. A cell’s sense of direction. Science, 284(5415):765– 770, APR 1999. [89] H.-O. Park and E. Bi. Central Roles of Small GTPases in the Development of Cell Polarity in Yeast and Beyond. Microbiol. Mol. Biol. Rev., 71(1):48–96, 2007. 152 Chapter 7. Bibliography [90] O. Pertz, L. Hodgson, R. Klemke, and K. Hahn. Spatiotemporal dynamics of rhoa activity in migrating cells. Nature, 440(7087):1069–1072, 2006. [91] C. Peskin, G. Odell, and G. Oster. Cellular motions and thermal fluctuations: the brownian ratchet. Biophysical Journal, 65(1):316 – 324, 1993. [92] T. D. Pollard, L. Blanchoin, and R. D. Mullins. Molecular mechanisms controlling actin filament dynamics in nonmuscle cells. Annual Review of Biophysics and Biomolecular Structure, 29(1):545–576, 2000. [93] T. D. Pollard and W. C. Earnshaw. Cell Biology. Elsevier, 2008. [94] M. Postma, L. Bosgraaf, H. Loovers, and P. Van Haastert. Chemotaxis: signalling modules join hands at front and tail. Embo Reports, 5(1):35–40, JAN 2004. [95] A. Prat and Y.-X. Li. Stability of front solutions in inhomogeneous media. Physica D: Nonlinear Phenomena, 186(1-2):50 – 68, 2003. [96] M. Raftopoulou and A. Hall. Cell migration: Rho GTPases lead the way. Dev. Biol., 265(1):23–32, JAN 2004. [97] A. Ridley. Rho family proteins: coordinating cell responses. Trends Cell Biol., 11(12):471–477, DEC 2001. [98] A. Ridley. Rho GTPases and cell migration. J. Cell Sci., 114(15):2713–2722, AUG 2001. [99] J. Rubinstein and P. Sternberg. Nonlocal reaction–diffusion equations and nucle- ation. IMA Journal of Applied Mathematics, 48(3):249, 1992. [100] J. Rubinstein, P. Sternberg, and J. Keller. Fast reaction, slow diffusion, and curve shortening. SIAM Journal on Applied Mathematics, 49(1):116–133, 1989. [101] S. Ruuth. Implicit-explicit methods for reaction-diffusion problems in pattern formation. J Math Biol, 34:148–176, 1995. 153 Chapter 7. Bibliography [102] Y. Sako, K. Hibino, T. Miyauchi, Y. Miyamoto, M. Ueda, and T. Yanagida. Single- molecule imaging of signaling molecules in living cells. Single Mol., 1(2):159–163, 2000. [103] R. A. Satnoianu, M. Menzinger, and P. K. Maini. Turing instabilities in general system. Journal of mathematical biology, 4. [104] A. Schmidt and A. Hall. Guanine nucleotide exchange factors for Rho GTPases: turning on the switch. Genes & Development, 16(13):1587–1609, JUL 2002. [105] R. Seydel. From Equilibrium to Chaos: Practical Bifurcation and Stability Anal- ysis. Elsevier, 1998. [106] J. Smoller and A. Wasserman. Global bifurcation of steady state solutions. Journal of Differential Equations, 39:269–290, 1981. [107] I. Stakgold. Green’s functions and boundary value problems. John Wiley and Sons, 1979. [108] K. Subramanian and A. Narang. A mechanistic model for eukaryotic gradient sensing: Spontaneous and induced phosphoinositide polarization. J. Theor. Biol., 231(1):49–67, NOV 2004. [109] T. M. Svitkina, A. B. Verkhovsky, K. M. McQuade, and G. G. Borisy. Analysis of the Actin-Myosin II System in Fish Epidermal Keratocytes: Mechanism of Cell Body Translocation. J. Cell Biol., 139(2):397–415, 1997. [110] Y. Takai, T. Sasaki, and T. Matozaki. Small GTP-binding proteins. Physiological Reviews, 81(1):153–208, JAN 2001. [111] A. Turing. The Chemical Basis of Morphogenesis. Philosophical Transactions of the Royal Society (B), 237:37–72, 1952. [112] J. Tyson, K. Chen, and B. Novak. Sniffers, buzzers, toggles and blinkers: dynamics 154 Chapter 7. Bibliography of regulatory and signaling pathways in the cell. Curr Opin Cell Biol, 15(2):221– 231, 2003. [113] P. Van Haastert and P. Devreotes. Chemotaxis: signalling the way forward. Nature Reviews Molecular Cell Biology, 6(8):626–634, 2004. [114] A. Van Keymeulen, K. Wong, Z. A. Knight, C. Govaerts, K. M. Hahn, K. M. Shokat, and H. R. Bourne. To stabilize neutrophil polarity, PIP3 and Cdc42 augment RhoA activity at the back as well as signals at the front. J. Cell Biol., 174(3):437–445, 2006. [115] A. B. Verkhovsky, T. M. Svitkina, and G. G. Borisy. Self-polarization and direc- tional motility of cytoplasm. Curr. Biol., 9(1):11–20, Jan. 1999. [116] L. Wang and M. Y. Li. Diffusion-driven instability in reaction-diffusion systems. Journal of Mathematical Analysis and Applications, 254(1):138 – 153, 2001. [117] M. Ward. Metastable Bubble Solutions for the Allen-Cahn Equation with Mass Conservation. SIAM Journal on Applied Mathematics, 56(5):1247–1279, 1996. [118] R. Wedlich-Soldner, S. C. Wai, T. Schmidt, and R. Li. Robust cell polarity is a dynamic state established by coupling transport and GTPase signaling. J. Cell Biol., 166(6):889–900, 2004. [119] S. S. Willard and P. N. Devreotes. Signaling pathways mediating chemotaxis in the social amoeba, dictyostelium discoideum. European Journal of Cell Biology, 85(9-10):897 – 904, 2006. A special issue dedicated to Guenter Gerisch on the occasion of his 75th birthday. [120] K. Wong, O. Pertz, K. Hahn, and H. Bourne. Neutrophil polarization: Spatiotem- poral dynamics of RhoA activity support a self-organizing mechanism. Proc. Natl. Acad. Sci. USA, 103(10):3639–3644, MAR 2006. 155 Chapter 7. Bibliography [121] J. Xu, F. Wang, A. Van Keymeulen, P. Herzmark, A. Straight, K. Kelly, Y. Takuwa, N. Sugimoto, T. Mitchison, and H. Bourne. Divergent signals and cytoskeletal assemblies regulate self-organizing polarity in neutrophils. Cell, 114(2):201–214, JUL 2003. [122] P. T. Yam, C. A. Wilson, L. Ji, B. Hebert, E. L. Barnhart, N. A. Dye, P. W. Wiseman, G. Danuser, and J. A. Theriot. Actin myosin network reorganization breaks symmetry at the cell rear to spontaneously initiate polarized cell motility. J. Cell Biol., 178(7):1207–1221, 2007. [123] M. Zerial and L. Huber, editors. Guidebook to the small GTPases. Oxford Uni- versity Press., 1995. [124] B. Zhang and Y. Zheng. Regulation of RhoA GTP hydrolysis by the GTPase- activating proteins p190, p50RhoGAP, Bcr, and 3BP-1. Biochem., 37(15):5249– 5257, APR 1998. [125] S. H. Zigmond, M. Joyce, C. Yang, K. Brown, M. Huang, and M. Pring. Mecha- nism of Cdc42-induced actin polymerization in neutrophil extracts. J. Cell Biol., 142(4):1001–1012, Aug. 1998. 156 Appendices 157 Appendix A Numerical Simulation Details A.1 Stimuli repertoire Adding stimulus to the wave-pinning equations results in the following system (2.12): ∂u ∂t = Du ∂2u ∂x2 + f(u, v) + kS(x, t)v, (A.1a) ∂v ∂t = Dv ∂2v ∂x2 − f(u, v)− kS(x, t)v. (A.1b) The specific form of kS used in simulations of (2.12) as well are listed below: Figure 2.8(a): Transient localized stimulus: We stimulated 10% of the domain, using klocS = s(t)(1 + cospix), 0 ≤ x ≤ L/10, 0 otherwise. where the time dependence was as follows: s(t) = S 2 0 ≤ t ≤ t1, S 4 (1 + cos ( pi (t−t1)(t2−t1) ) , t1 ≤ t ≤ t2, 0, otherwise. S = 0.05s−1, t1 = 20s, t2 = 25s were used. Figure 2.8(b): Transient graded stimulus: kgradS = s(t)(10− x), 0 ≤ x ≤ L, 158 A.1. Stimuli repertoire where s(t) = S, 0 ≤ t ≤ t1, S ( 1− t−t1t2−t1 ) , t1 ≤ t ≤ t2, 0, otherwise. (A.2) S = 0.07s−1, t1 = 20s, t2 = 25s were used. Figure 2.8(c): Reversal of graded stimulus:We applied same stimulus kgradS as in (b) followed by kgradS = s(t)x, 0 ≤ x ≤ L where s(t) = S, 200 ≤ t ≤ t1, S ( 1− t−t1t2−t1 ) , t1 ≤ t ≤ t2, 0, otherwise. (A.3) S = 0.07s−1, t1 = 120s, t2 = 125s were used. Figure 2.8(d): Random initial conditions: Here kS ≡ 0. We used nonuniform initial condi- tions of the form v(x, 0) = v0 , u(x, 0) = u−(v0) · 2R, with a uniformly distributed random number R = UNIF (0, 1) so that the noise has zero mean about the ho- mogeneous steady state. Using u−(v0) = 0.268v0 = 2.0, lead to some outcomes with polarization (as in Fig. 2.8d), some with long-lived multiple peaks, and some with decay back to homogenous distribution. 159 Appendix B Mathematical Background B.1 Conditions for Diffusion-Driven Instability A RD system exhibits diffusion-driven (also known as Turing) instability if the ho- mogeneous steady state is stable to small perturbations in the absence of diffusion, but unstable to small spatial perturbations when diffusion is present. For simplicity, consider a two component system: ∂u ∂t = ∆u+ f(u, v) on Ω, (B.1a) ∂v ∂t = d∆v + g(u, v) on Ω, (B.1b) un = 0, vn = 0 on ∂Ω, (B.1c) where n is the unit outward normal to ∂Ω. System (B.1) has been scaled so that d is the ratio of the diffusion coefficients of u and v with the assumption that d 6= 1. We assume that in the absence of diffusion the homogeneous steady state (u0, v0) given by solution of f(u, v) = 0, g(u, v) = 0 (B.2) is stable. That is, letting w = u− u0 v − v0 (B.3) we obtain from the ODE form of (B.1) wt = γAw, (B.4) 160 B.1. Conditions for Diffusion-Driven Instability where A is the Jacobian given by A = fu fv gu gv evaluated at (u0, v0). For linear stability of (u0, v0) it is required that trace(A) = fu + gv < 0, |A| = fugv − fvgu > 0. (B.5) Linearizing the full system (B.1) about (u0, v0) we obtain wt = Aw +D∆w, (B.6) where D = 1 0 0 d . Look for solutions of the form w = ∑ k cke λtWk, (B.7) whereWk is the eigenfunction corresponding to wave-number k of the eigenvalue prob- lem ∆W + k2W = 0 on D, Wn = 0 on ∂Ω. (B.8) Substituting this form of w into (B.6), we obtain the following equation for each wavenumber k: λWk = AWk −Dk2Wk. (B.9) Nontrivial solutions forWk are determined by the roots of the characteristic polynomial |λI −A+Dk2| = 0, 161 B.1. Conditions for Diffusion-Driven Instability which evaluates to λ2 + λ[k2(1 + d)− (fu + gv)] + [dk4 − (dfu + gv)k2 + γ2|A|]. (B.10) By assumption that the steady state is stable in the absence of diffusion we have that Reλ(k2 = 0) = 0, which give rise to conditions (B.5). For the steady state to be unstable to spatial perturbations we require that Reλ(k2 > 0) > 0. Now, the coefficient of λ is always greater than zero, as k2(1 + d) > 0 and fu + gv < 0 by (B.5). So, the only way that Reλ(k2 > 0) > 0 is if H(k2) = dk4 − (dfu + gv)k2 + |A| < 0 (B.11) for some k2. As the determinant of A, |A| > 0 by (B.5), it follows that the only way H(k2) < 0 is if dfu+gv > 0. As fu+gv < 0 it follows that fu and gv must have opposite signs and that d 6= 1. Differentiating H with respect to k2 we find that the minimum value of H occurs at Hmin = |A| − (df)u+ gv) 2 4d at k2 = dfu + gv 2d . (B.12) It follows that the condition for H(k2) < 0 is given by (df)u+ gv)2 4d > |A| =⇒ (dfu + gv)2 − 4d(fugv − f)vgu) > 0. (B.13) At the bifurcation point H(k2) = 0, and the critical wave-number k2c = dfu+gv 2d becomes unstable16. Close to the bifurcation point, this wave-number is the fastest-growing mode in the emerging spatially heterogeneous solution. At the onset of Turing instability the speed at which this mode grows is approximated by expλ(k2)t, where λ(k2) is the largest eigenvalue of (B.10). In summary, the conditions for diffusive-driven instability of steady state (u0, v0) in 16On a finite domain there is a discrete set of wavenumbers for the eigenvalue problem (B.8), so the wavenumber closest to kc will be the first to become unstable. 162 B.2. Stability of the Homogeneous States for the Wave-Pinning System a two component system are fu + gv < 0, fugv − fvgu > 0, (B.14a) dfu + gv > 0, (dfu + gv)2 − 4d(fugv − f)vgu) > 0, (B.14b) where all derivatives are evaluated at (u0, v0). Note that fu and gv must have opposite signs, so it follows by condition (3) that fvgu < 0. In the case that fv < 0, gu > 0 we have an activator-inhibitor system. In the case that fv > 0, gu < 0 the Turing system is of the substrate depletion type. B.2 Stability of the Homogeneous States for the Wave-Pinning System Here, we show that steady states u− and u+ of the wave-pinning system (4.10) are linearly stable homogeneous solutions of the PDE system. Let (u0, v0) be a homogeneous solution to (4.10), and (ũ, ṽ) be perturbations about that homogeneous steady state. Then from (4.10) we obtain the following linear system ² ∂ũ ∂t = ²2 ∂2ũ ∂x2 + fuũ+ fvṽ, (B.15a) ² ∂ṽ ∂t = D ∂2ṽ ∂x2 − fuũ− fvṽ (B.15b) on x = [0, L] with no-flux boundary conditions. Consider perturbations of the form ũ ṽ = α1 α2 cos kx eσt (B.16) where k = npi/L, for L = 1, n = 1, 2, , cdots. These types of perturbations satisfy the no-flux boundary conditions and preserve the mass constraint (4.11). Substituting 163 B.2. Stability of the Homogeneous States for the Wave-Pinning System (B.16) into (B.15) and simplifying we obtain the following system of linear equations: α1(²σ + ²2k2 − fu) + α2(−fv) = 0, (B.17a) α1(fu) + α2(²σ +Dk2 + fv) = 0. (B.17b) For a non-trivial solution we require that the determinant of the system be equal to 0: ∣∣∣∣∣∣∣ ²σ + ²2k2 − fu −fv fu ²σ +Dk2 + fv ∣∣∣∣∣∣∣ = 0. (B.18) Solving for the growth rate of the perturbations we obtain σ = 1 2² ( −τ ± √ τ2 − 4∆ ) = 1 2² ( −[(D + ²2)k2 + fv − fu]± √ ((D + ²2)k2 + fv − fu)2 − 4(D²2 + ²2fv −Dfu) ) . (B.19) For Re(σ) < 0 we require that τ > 0 for all k2 > 0, or equivalently that ( ∂f ∂u − ∂f ∂v )∣∣∣∣ u=u±(v) > 0, (B.20) which is assumption (4.14) on f(u, v). In particular, for the cubic reaction kinetics (4.12) f(u, v) = u(1− u)(u− 1− v), (B.21) we have fu = (1− u)(u− 1− v) + u[−(u− 1− v) + (1− u)], fv = −u(1− u) (B.22) 164 B.3. Explicit front speed formula for the cubic kinetics function which results in fu|u=u−(v) = −(1 + v), fv|u=u−(v) = 0, (B.23) fu|u=um(v) = v, fv|u=um(v) = 0, (B.24) fu|u=u+(v) = −v(1 + v), fv|u=u−(v) = v(1 + v). (B.25) It follows that for u = u−(v) = 0, σ = −Dk 2 ² , −² 2k2 + 1 + v ² , (B.26) for u = um(v) = 1, σ = −Dk 2 ² , v − ²2k2 ² . (B.27) Note that for epsilon2k2 > v, the growth rate σ in (B.27) becomes positive. For u = u+(v) = 1 + v, σ = − 1 2² ( (D + ²2)k2 + 2v(1 + v)± √ (²2 −D)2k4 + 4v2(1 + v)2 ) . (B.28) It follows that u = u±(v) are linearly stable homogeneous states of system (4.10), while u = um(v) is linearly unstable. B.3 Explicit front speed formula for the cubic kinetics function Consider a single RD equation (4.5) in 1D for u(x, t) with v taken as parameter. In general, (and specifically for reaction terms of type (4.4)), an explicit form for the propagating front solution U0 in the denominator of (4.29) is not obtainable. However, the simplicity of the cubic reaction kinetics (4.3) allows us to obtain a closed form 165 B.3. Explicit front speed formula for the cubic kinetics function solution to (4.5), and thus evaluate (4.29) explicitly. We consider f cubic(u, v) = λ(u− u1)(u2 − u)(u− u3) (B.29) where λ is a positive constant and u1 < u2 < u3. The reaction diffusion equation (4.5) is then given by ∂u ∂t = D ∂2u ∂x2 + λ(u− u1)(u2 − u)(u− u3). (B.30) Assume that (B.30) has wavefront solutions of the form u(x, t) = U(z), where z = x− ct, U(−∞) = u3, U(∞) = u1. Substituting into (B.30) gives DU ′′ + cU ′ + λ(U − u1)(u2 − U)(U − u3) = 0. (B.31) Following (77) assume that U satisfies U ′ = a(U − u1)(U − u3), (B.32) where a is to be determined. Substituting (B.32) into (B.31) and simplifying obtain (U − u1)(U − u3) [ (2Da2 − λ)U − (Da2(u1 + u3)− ca− λu2) ] = 0. (B.33) It follows that 2Da2 − λ = 0 and Da2(u1 + u3)− ca− λu2 = 0, (B.34) which gives a = √ λ 2D and c = √ λD 2 (u1 − 2u2 + u3). (B.35) 166 B.4. Solvability Conditions Note that substituting (B.32) into the integral formula (4.6) gives the same result for the wave speed c. It follows from (B.32) that U(z) = u3 + u1 2 − u3 − u1 2 tanh ( (u3 − u1)z 2 √ 2 ) . (B.36) B.4 Solvability Conditions Consider the Poisson equation subject to no-flux boundary conditions: ∆u = f on Ω (B.37a) ∂un = 0 on ∂Ω. (B.37b) Integrating over Ω and using the boundary conditions it can be seen that equa- tion (B.37a) has a solution if ∫ Ω f dx = 0. (B.37c) This is both a necessary and sufficient condition (107) that has a simple physical inter- pretation. Equation (B.37a) describes the steady state of a heat equation subject to a heat source f . Since heat cannot escape the domain due to no-flux boundary conditions, the only way that a steady state can exist is for net generation of heat in the domain to be zero, which is given by (B.37c). In general, for a linear operator A on a closed, bounded domain Ω, we need a solvability condition on the forcing term f to ensure that Au = f on Ω (B.38) has a solution. The adjoint of A is defined to be the operatorA∗ for which< Au, v >=< u,A∗v > for all u in the domain of A and v in the domain of A∗. Typically, for 167 B.5. Critical Epsilon Calculation for Emergence of Nonhomogeneous Solution Au ∈ L2[a, b], the inner product is given by < u, v >= ∫ b a u(x)v(x)w(x) dx, where w is a weight function. In general, integration by parts will yield an expression of the form < Au, v >= J(u, v)+ < u,A∗v > where J(u, v) represents the boundary conditions on the operator A. The Fredholm Alternative (45; 107) states that for a differential operator A with boundary operator B, a solution to Au = f, J(u, v) = g exists if and only if < u,A∗v >=< Au, v > −J(u, v) =< f, v > −J(u, v) = 0 for all v satisfying A∗v = 0, B∗v = 0. It follows that a solution to (B.38) can be found if and only if f is orthogonal to the null space of the adjoint operator. B.5 Critical Epsilon Calculation for Emergence of Nonhomogeneous Solution Consider system (6.2) together with no flux boundary conditions. Recall that in this dimensionless system, the length of the domain has been scaled to 1. We think of (6.2) as an ODE system for u and v, where x is the time variable. Add (6.2) to obtain ²2uxx +Dvxx = 0 (B.39) Integrating twice and using the no flux boundary conditions we obtain ²2u+Dv = A, ⇒ v = A D − ² 2 D u. (B.40) 168 B.5. Critical Epsilon Calculation for Emergence of Nonhomogeneous Solution Substituting the form of v from (B.40) into (6.2a), we rewrite the steady state system as a single equation 0 = ²2uxx + f̃(u,A), (B.41) where f̃(u,A) = f(u, v(A)). We can linearize equation (B.41) around the unstable steady state (uuss) and find the critical value of ² required for emergence of nonhomogeneous solution from the un- stable homogeneous steady state. This value indicates an ² value at which the unstable spatially heterogeneous branch of solutions merges with the the unstable homogeneous steady state. The linearized system 0 = uxx + 1 ²2 f̃ ′(uuss)u (B.42a) ux(0) = ux(1) = 0 (B.42b) has a first non-trivial solution emerge when 1 ²2crit f̃ ′(uuss) = pi2. (B.43) For example for f(u, v) = u(m − u)(u − m − v), with the unstable steady state uuss = m, it follows that f̃(u,A) = u(m− u)(u−m− A D + ²2 D u, (B.44) with f̃ ′(m) = m ( A D − ² 2 D m ) . (B.45) This results in ²crit = √ A/D pi2/m+m/D . (B.46) Since the value of A depends on the solution u the relation stated in this form is not 169 B.6. Integral Evaluation for the Time Map System very useful. Note, however, that using the conservation of total amount we can rewrite A/D as a function of u and the total amount K: A D = K − ( 1− ² 2 D )∫ 1 0 u dx. (B.47) Note that u = uuss at ²crit, so for equation (B.44) A D = K − ( 1− ² 2 D ) m ≈ K −m (B.48) for ²2 ¿ D. Equation (B.46) then becomes ²crit = √ K −m pi2/m+m/D . (B.49) For K = 2m as D →∞ ²crit → mpi . B.6 Integral Evaluation for the Time Map System Note that the endpoints u0 and u1 of the integral ∫ u1 u0 du√ F (u,A,B) , (B.50) appearing in (6.25) satisfy F (u,A,B) = 0 since ux(0) = ux(1) = 0. This means that the two integral relations (6.25b) and (6.25c) have a square root singularity at the endpoints. These types of singularities can still be evaluated numerically in the newer versions of Matlab (using quadgk function or adaptive Gauss-Kronrod quadrature). To evaluate the integral using more traditional methods such as adaptive Simpson quadrature, we need to rescale the integral in order to eliminate the singularity. This is done in three steps. First, rewrite F (u,A,B) = αu4 + βu3 + γu2 +B as F (u,A,B) = (u− u0)(u− u1) ( αu2 + (β + α(u0 + u1))u+ B u0u1 ) 170 B.6. Integral Evaluation for the Time Map System where a = α 4 (u1 − u0)2, (B.51a) b = α(u21 − u20) + β 2 (u1 − u0), (B.51b) c = α 4 (u0 + u1)2 + (β + α(u0 + u1)) ( u0 + u1 2 ) + B u0u1 . (B.51c) Then we rescale the integral using the transformation p = 2 u1 − u0u− u0 + u1 u1 − u0 (B.52) to get ∫ u1 u0 du√ F (u,A,B) = ∫ 1 −1 dp√ (p+ 1)(p− 1)(ap2 + bp+ c) . (B.53) Finally, to eliminate the singularity at the endpoints of the integral, we make the following variable transformation: p = 1 2 (3s− s3), (B.54) giving us I3 ≡ ∫ u1 u0 du√ F (u,A,B) = 3 ∫ 1 −1 ds√ (s− 2)(s+ 2) (a4 (3s− s3)2 + b2(3s− s3) + c) . (B.55) Similarly, I4 = ∫ u1 u0 u du√ F (u,A,B) = 3(u1 − u0) 4 ∫ 1 −1 (3s− s2) ds√ (s− 2)(s+ 2) (a4 (3s− s3)2 + b2(3s− s3) + c) + (u0 + u1) 2 I3. (B.56) To verify that the integral transformation is correct both integrals (6.25b) and (B.55) were evaluated in Matlab (using quadgk function for (6.25b) and quad function for 171 B.6. Integral Evaluation for the Time Map System (B.55)) using the values of A and B calculated from the time-dependent problem as described above, obtaining the correct result. Similarly, equations (6.25c) and (B.56) give the same result. Another problem with evaluating (6.25b) and (6.25c), is that for B ≈ 0 it is difficult to accurately determine the roots of F (u,A,B). In particular, at B = 0 the two roots around 0 coalesce, and it becomes impossible to determine u0 (see Figure 6.4). This is alleviated by splitting the integral into two parts: ∫ u1 u0 du√ F (u,A,B) = ∫ η u0 du√ F (u,A,B) + ∫ u1 η du√ F (u,A,B) , (B.57) and doing a different integral transformation for the first integral. In particular, close to u ≈ 0 we can rewrite F as F = (u− δ2)(αu2 + βu+ γ). Then ∫ η u0 du√ F (u,A,B) becomes ∫ η δ du√ (u− δ2)(αu2 + βu+ γ) . (B.58) We rescale the integral, letting u = δ cosh z to get ∫ cosh−1(η/δ) 0 dz√ αδ2 cosh2 z + βδ cosh z + γ . (B.59) Letting λ = cosh−1 (η δ ) ⇔ δ = η coshλ (B.60) we obtain ∫ λ 0 dz√ αη2 ( cosh z coshλ )2 + βη ( cosh z coshλ ) + γ . (B.61) 172
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A wave-pinning mechanism for eukaryotic cell polarization...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A wave-pinning mechanism for eukaryotic cell polarization based on Rho GTPase dynamics Jilkine, Alexandra 2009
pdf
Page Metadata
Item Metadata
Title | A wave-pinning mechanism for eukaryotic cell polarization based on Rho GTPase dynamics |
Creator |
Jilkine, Alexandra |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | In response to chemical stimulation, many eukaryotic cells are able to sense the direction of the stimulus and initiate movement. To do so, the cell must break symmetry and develop a front and back in a process known as polarization. During polarization, members of the Rho GTPase family (Cdc42, Rac, and Rho) are recruited to the plasma membrane and localize to form a front and a back of the polarizing cell. These proteins exist in both active forms (on the inner surface of the membrane of the cell), and inactive forms (in the cytosol). In earlier work, I have shown that the property of membrane-cytosol interconversion, together with appropriate feedbacks, endows the Rho proteins with the ability to initiate cell polarity, resulting in a high Cdc42/Rac region, which will become the front, and a high Rho region, which will become the back of the cell. Here I show that this property of polarizability can be explained using a simplified model system comprising of a single active/inactive protein pair with positive feedback. In this model, a travelling wave of GTPase activation is initiated at one end of the domain, moves across the cell, and eventually stops inside the domain, resulting in a stable polar distribution. The key requirements for the mechanism to work include conservation of total amount of protein, a sufficiently large difference in diffusion rates of the two forms, and nonlinear positive feedback that allows for multiple homogeneous steady states to exist. Using singular perturbation theory, I explain the mathematical basis of wave-pinning behaviour, and discuss its biological and mathematical implications. I show that this mechanism for generating a chemical pattern is distinct from Turing pattern formation. I also analyze the transition from a spatially heterogeneous solution to a spatially homogeneous solution as the diffusion coefficient of the active form is increased, and show the existence of other unstable stationary wavefronts. Finally, I argue that this wave-pinning mechanism can account for a number of features of cell polarization such as spatial amplification, maintenance of polarity, and the sensitivity to new stimuli that is typical of polarization of eukaryotic cells. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0068766 |
URI | http://hdl.handle.net/2429/17453 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_jilkine_alexandra.pdf [ 6.55MB ]
- Metadata
- JSON: 24-1.0068766.json
- JSON-LD: 24-1.0068766-ld.json
- RDF/XML (Pretty): 24-1.0068766-rdf.xml
- RDF/JSON: 24-1.0068766-rdf.json
- Turtle: 24-1.0068766-turtle.txt
- N-Triples: 24-1.0068766-rdf-ntriples.txt
- Original Record: 24-1.0068766-source.json
- Full Text
- 24-1.0068766-fulltext.txt
- Citation
- 24-1.0068766.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0068766/manifest