Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Comparing cell polarization models using local perturbation analysis DuTôt, Meghan 2014

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

Item Metadata


24-ubc_2014_spring_dutot_meghan.pdf [ 2.18MB ]
JSON: 24-1.0166907.json
JSON-LD: 24-1.0166907-ld.json
RDF/XML (Pretty): 24-1.0166907-rdf.xml
RDF/JSON: 24-1.0166907-rdf.json
Turtle: 24-1.0166907-turtle.txt
N-Triples: 24-1.0166907-rdf-ntriples.txt
Original Record: 24-1.0166907-source.json
Full Text

Full Text

Comparing cell polarization models using localperturbation analysisbyMeghan DuToˆtB.Sc., University of British Columbia, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)The University Of British Columbia(Vancouver)April 2014c©Meghan DuToˆt, 2014AbstractPatterns are ubiquitous in nature, but the underlying mechanisms giving rise tothem are not always understood. One such pattern is the inhomogeneous domainof a cell that is established during polarization. Cell polarization is a way in whichcells respond to and interact with their environment. For example, white bloodcells locate and destroy bacteria, and yeast cells create buds for reproduction. Sig-nalling proteins such as GTPases are redistributed throughout the cell and, throughdownstream effects, rearrangement of the actin cytoskeleton follows. This redis-tribution can occur in response to a stimulus, such a chemoattractant, or it may bespontaneous. Because many biological details are unknown, mathematical modelsare developed to recreate features of cell polarization and determine the minimalmodules or characteristics for these features. Polarization models are often sim-ple, conceptual reaction-diffusion equations for one or more signalling molecules.But comparing these models is often difficult, and there are many models in theliterature for different cell types or behaviour.A new method of nonlinear stability analysis, called local perturbation analy-sis (LPA), was developed by Stan Mare´e, and later Bill Holmes, to take advantageof models with substantial diffusion disparities. This method recapitulates the dy-namics of a pulse applied to a reaction-diffusion system using a system of ordinarydifferential equations. Bifurcation analysis of these equations is relatively easy,and LPA detects pattern formation through threshold and Turing dynamics and pro-vides bifurcation maps of these regimes for any parameter. LPA is well-suited tocell polarization models, because the signalling proteins we model often have bothfast-diffusing inactive and slow-diffusing active forms. In this thesis, I introduceLPA and its use through a wave pinning model and its extension, a model for actiniiwaves. I then review and analyze five additional cell polarization models usingcombinations of LPA, simulations, and Turing analysis. In many cases, I discov-ered new dynamics of the models. LPA helps us to map patterning regimes and theirrobustness to changes in parameters, and provides a new avenue for us to comparemany current models for cell polarization.iiiPrefaceThis dissertation is the original intellectual property of the author, M. Dutot. Asoutlined below, versions of the material have been published, but none of the textor figures of the dissertation are taken directly from these articles.Some material in this thesis has been published in:1. Edelstein-Keshet, L., Holmes, W.R., Zajac, M., and Dutot, M. (2013). Fromsimple to detailed models for cell polarization. Philosophical Transactionsof the Royal Society B: Biological Sciences, 368(1629): doi:10.1098/rstb.2013.0003.2. Mata, M.A., Dutot, M., Edelstein-Keshet, L., Holmes, W.R. (2013). A modelfor intracellular actin waves explored by nonlinear local perturbation analy-sis. Journal of Theoretical Biology, 334:149-161.My contribution to these works was: participating in preliminary discussions to de-sign the research with members of the Keshet group; reading and presenting resultsof polarization models; using software in Matlab to run bifurcation analysis on LPAsystems; creating and running simulation programs (using Matlab) to produce PDEsolutions; creating LPA bifurcation diagrams and publication-quality figures; ana-lyzing pulse width and speed in a model for actin waves (analysis omitted from thisthesis); providing summaries of my own work for these publications; and readingand commenting on drafts of the manuscripts.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Biological background . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Actin’ up: key components in an intracellular pathway . . . . . . 52.1.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Polarization in different cell types . . . . . . . . . . . . . . . . . 92.2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Mathematical modelling background . . . . . . . . . . . . . . . . . 163.1 Introduction to modelling with reaction-diffusion equations . . . . 16v3.2 Patterns in reaction-diffusion systems . . . . . . . . . . . . . . . 193.2.1 Turing patterns . . . . . . . . . . . . . . . . . . . . . . . 193.2.2 Travelling waves . . . . . . . . . . . . . . . . . . . . . . 203.2.3 Local excitation / global inhibition . . . . . . . . . . . . . 213.3 Models for cell polarization . . . . . . . . . . . . . . . . . . . . . 223.3.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 273.4 Connecting models and biology . . . . . . . . . . . . . . . . . . 283.4.1 Abstracting features of cell polarization . . . . . . . . . . 283.4.2 Polarization modelling questions . . . . . . . . . . . . . . 314 Local perturbation analysis . . . . . . . . . . . . . . . . . . . . . . . 324.1 Current methods of analysis for reaction-diffusion equations . . . 324.1.1 Numerical solutions . . . . . . . . . . . . . . . . . . . . 334.1.2 Turing analysis . . . . . . . . . . . . . . . . . . . . . . . 344.1.3 Nonlinear stability analysis . . . . . . . . . . . . . . . . . 354.2 Local perturbation analysis . . . . . . . . . . . . . . . . . . . . . 354.2.1 Mass conservative systems . . . . . . . . . . . . . . . . . 374.2.2 Model signatures . . . . . . . . . . . . . . . . . . . . . . 384.3 Benefits and caveats . . . . . . . . . . . . . . . . . . . . . . . . . 395 Demonstrating local perturbation analysis through examples . . . . 415.1 Wave Pinning model . . . . . . . . . . . . . . . . . . . . . . . . 425.1.1 Model equations . . . . . . . . . . . . . . . . . . . . . . 425.1.2 Preliminary investigations . . . . . . . . . . . . . . . . . 435.1.3 The LPA system . . . . . . . . . . . . . . . . . . . . . . . 445.1.4 LPA bifurcation analysis . . . . . . . . . . . . . . . . . . 465.1.5 PDE simulations . . . . . . . . . . . . . . . . . . . . . . 485.2 Actin Waves model . . . . . . . . . . . . . . . . . . . . . . . . . 495.2.1 Model introduction . . . . . . . . . . . . . . . . . . . . . 495.2.2 Preliminary investigations . . . . . . . . . . . . . . . . . 535.2.3 The LPA system . . . . . . . . . . . . . . . . . . . . . . . 535.2.4 LPA bifurcation analysis . . . . . . . . . . . . . . . . . . 565.2.5 PDE simulations and connection to LPA . . . . . . . . . . 59vi5.3 Wave pinning versus actin waves . . . . . . . . . . . . . . . . . . 605.4 Additional application: periodic cell migration . . . . . . . . . . . 626 Analyzing cell polarization models . . . . . . . . . . . . . . . . . . . 656.1 Goryachev Model . . . . . . . . . . . . . . . . . . . . . . . . . . 676.1.1 LPA system and bifurcation diagrams . . . . . . . . . . . 706.1.2 Turing analysis . . . . . . . . . . . . . . . . . . . . . . . 706.1.3 PDE simulations and connection to LPA . . . . . . . . . . 746.2 Otsuji Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 766.2.1 LPA system and bifurcation diagrams . . . . . . . . . . . 796.2.2 Turing analysis . . . . . . . . . . . . . . . . . . . . . . . 806.2.3 PDE simulations and connection to LPA . . . . . . . . . . 826.3 Meinhardt Activator-Substrate Depletion Model . . . . . . . . . . 836.3.1 LPA system and bifurcation diagrams . . . . . . . . . . . 876.3.2 PDE simulations and connection to LPA . . . . . . . . . . 896.4 Pseudopod Centered Model . . . . . . . . . . . . . . . . . . . . . 906.4.1 LPA system and bifurcation diagrams . . . . . . . . . . . 956.4.2 PDE simulations and connection to LPA . . . . . . . . . . 966.5 Balanced Inactivation Model . . . . . . . . . . . . . . . . . . . . 996.5.1 Model modifications . . . . . . . . . . . . . . . . . . . . 1016.5.2 LPA system and bifurcation diagrams . . . . . . . . . . . 1046.5.3 PDE simulations and connection to LPA . . . . . . . . . . 1076.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1097 Results and comparisons . . . . . . . . . . . . . . . . . . . . . . . . 1107.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . 1107.1.1 Wave Pinning Model . . . . . . . . . . . . . . . . . . . . 1107.1.2 Goryachev Model . . . . . . . . . . . . . . . . . . . . . . 1117.1.3 Otsuji Model . . . . . . . . . . . . . . . . . . . . . . . . 1127.1.4 Meinhardt Activator-Substrate Depletion Model . . . . . . 1127.1.5 Pseudopod Centered Model . . . . . . . . . . . . . . . . 1137.1.6 Balanced Inactivation Model . . . . . . . . . . . . . . . . 1147.2 Comparing models according to patterning regimes . . . . . . . . 115vii7.2.1 Turing instabilities . . . . . . . . . . . . . . . . . . . . . 1157.2.2 Threshold regimes (upward and downward) . . . . . . . . 1167.2.3 Stable regions . . . . . . . . . . . . . . . . . . . . . . . . 1177.3 Comparing models according to qualitative features . . . . . . . . 1187.4 Wave Pinning versus Meinhardt’s Activator-Substrate DepletionModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1207.5 Connecting patterning regimes to simulation stimuli . . . . . . . . 1217.5.1 Turing regimes . . . . . . . . . . . . . . . . . . . . . . . 1217.5.2 Threshold regimes . . . . . . . . . . . . . . . . . . . . . 1227.6 Comparing simple and detailed models . . . . . . . . . . . . . . . 1228 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1248.1 Conclusions concerning LPA . . . . . . . . . . . . . . . . . . . . 1258.1.1 LPA quickly reveals information about pattern dynamics . 1258.1.2 Model modifications and assumptions . . . . . . . . . . . 1268.1.3 Other methods provide additional information . . . . . . . 1268.2 LPA and cell polarization models . . . . . . . . . . . . . . . . . . 1288.2.1 Answering polarization modelling questions . . . . . . . . 1288.2.2 Classifying models . . . . . . . . . . . . . . . . . . . . . 1288.3 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . 130Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132A URL references . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141B Crank-Nicolson method . . . . . . . . . . . . . . . . . . . . . . . . . 142C Turing Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146D Derivation of the Thomas algorithm . . . . . . . . . . . . . . . . . . 151D.1 Forward sweep . . . . . . . . . . . . . . . . . . . . . . . . . . . 152D.2 Backward sweep . . . . . . . . . . . . . . . . . . . . . . . . . . 153viiiList of TablesTable 2.1 Summary of polarization in different cell types . . . . . . . . . 14Table 3.1 Summary of polarization models . . . . . . . . . . . . . . . . 28Table 3.2 Qualitative polarization features for cell types and models . . . 30Table 5.1 Wave Pinning model parameters . . . . . . . . . . . . . . . . 43Table 5.2 Actin Wave model parameters . . . . . . . . . . . . . . . . . . 51Table 6.1 Goryachev model parameters . . . . . . . . . . . . . . . . . . 68Table 6.2 Otsuji model parameters . . . . . . . . . . . . . . . . . . . . . 77Table 6.3 Meinhardt’s Activator-Substrate Depletion model parameters . 85Table 6.4 Pseudopod Centered model parameters . . . . . . . . . . . . . 93Table 6.5 Balanced Inactivation model parameters . . . . . . . . . . . . 100Table 7.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . 111ixList of FiguresFigure 1.1 Simple diagram of cell polarization . . . . . . . . . . . . . . 2Figure 2.1 Schematic diagram of the actin network . . . . . . . . . . . . 7Figure 2.2 Rho GTPase dynamics . . . . . . . . . . . . . . . . . . . . . 8Figure 2.3 Micrographs of six different cell types . . . . . . . . . . . . . 10Figure 3.1 Reaction diagrams . . . . . . . . . . . . . . . . . . . . . . . 22Figure 4.1 Visualization of the LPA system . . . . . . . . . . . . . . . . 36Figure 5.1 Initial investigation of the Wave Pinning model . . . . . . . . 45Figure 5.2 Wave Pinning LPA bifurcation diagrams . . . . . . . . . . . . 47Figure 5.3 Wave Pinning kymographs . . . . . . . . . . . . . . . . . . . 49Figure 5.4 Reaction diagram of Actin Wave model . . . . . . . . . . . . 52Figure 5.5 Initial investigation of the Actin Wave model . . . . . . . . . 54Figure 5.6 Actin Wave LPA schematic . . . . . . . . . . . . . . . . . . . 55Figure 5.7 Actin Wave LPA bifurcation and eigenvalues . . . . . . . . . . 57Figure 5.8 Additional Actin Wave LPA bifurcations . . . . . . . . . . . . 59Figure 5.9 Actin Wave PDE simulations . . . . . . . . . . . . . . . . . . 61Figure 5.10 Wave Pinning LPA Bifurcation diagram for total concentration 64Figure 6.1 Initial investigation of the Goryachev model . . . . . . . . . . 69Figure 6.2 Goryachev LPA bifurcation diagrams . . . . . . . . . . . . . . 71Figure 6.3 Maximal Turing eigenvalues for Goryachev . . . . . . . . . . 74Figure 6.4 Goryachev kymographs . . . . . . . . . . . . . . . . . . . . . 75xFigure 6.5 Initial investigation of the Otsuji model . . . . . . . . . . . . 78Figure 6.6 Otsuji LPA bifurcation diagram . . . . . . . . . . . . . . . . . 80Figure 6.7 Maximal Turing eigenvalues for the Otsuji model . . . . . . . 83Figure 6.8 Otsuji kymographs . . . . . . . . . . . . . . . . . . . . . . . 84Figure 6.9 Initial investigation of Meinhardt’s Activator-Substrate Deple-tion model . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 6.10 Meinhardt’s Activator-Substrate Depletion model LPA bifurca-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Figure 6.11 Meinhardt’s Activator-Substrate Depletion model kymographs(ra) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Figure 6.12 Meinhardt’s Activator-Substrate Depletion model kymographs(sa) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 6.13 Initial investigation of the Pseudopod Centered model . . . . 94Figure 6.14 Pseudopod Centered LPA bifurcation diagrams . . . . . . . . 96Figure 6.15 Pseudopod Centered kymographs for s . . . . . . . . . . . . . 97Figure 6.16 Pseudopod Centered kymographs for rc . . . . . . . . . . . . 98Figure 6.17 Balanced Inactivation schematic diagrams . . . . . . . . . . . 102Figure 6.18 Initial investigation of the Balanced Inactivation model . . . . 105Figure 6.19 Balanced Inactivation LPA bifurcation diagrams . . . . . . . . 106Figure 6.20 Balanced Inactivation simulations (stimulus) . . . . . . . . . 107Figure 6.21 Balanced Inactivation simulations (parameters) . . . . . . . . 108xiGlossaryLPA local perturbation analysisRDE reaction-diffusion equationPDE partial differential equationODE ordinary differential equationASDM Activator-Substrate Depletion ModelHSS homogeneous steady stateNPF nucleation promoting factorLEGI local excitation / global inhibitionxiiAcknowledgmentsI would like to express appreciation to my supervisor Dr. Leah Edelstein-Keshet forcontinuing to encourage and support my research, and for providing valuable com-ments and insight into my work. I would also like to thank Dr. William Holmes forhis guidance and suggestions, as well as supervision over parts of this work. Spe-cial thanks to my mathletes Laura Liao, May Anne Mata, and Hildur Knu´tsdo´ttirfor generating great discussions during and after group meetings, beers, and fillingthe only bar in Naramata to sing karaoke.I would not have made it through writing and editing of this thesis withoutnumerous encouraging coffee shop sessions with my partner in crime ♥, and I ex-press special appreciation to Lost & Found cafe´ for their free coffee refills. I cannotexpress how grateful I am to my parents and Kaitlin for their support through allof my academic endeavours. I truly would not have made it this far without theirlove.This research was funded through a NSERC discovery grant to Leah Edelstein-Keshet, NSERC postgraduate scholarship to myself (Meghan DuToˆt), and partlysupported by National Institutes of Health under grant no. R01 GM086882 toAnders Carlsson.xiiiTo my cat, Lenny, who gave me the valuable contribution“23ghejrnjfrjjjjjjjjjjkkk,” which was subsequently removed during editing.xivChapter 1IntroductionStimulus-response pathways occur everywhere in nature; it is how living organ-isms, both macroscopic and microscopic, sense and respond to their environment.Take for example the response that we undergo when we see a particularly attrac-tive individual. The attractive person is the stimulus, and the stimulus is receivedas light collected on the retina and interpreted by the brain. After receiving thisstimulus, neurochemicals such as dopamine, norepinephrine, and serotonin are re-leased, making us feel giddy and nervous. We may respond calmly, or by blushing,sweating, and stumbling before deciding that our time is better spent on compu-tational biology. For another pathway on a different sensory system, consider thescent of biryani wafting down the street. The pungent combination of slow-cookedrice, aromatic spices, vegetables and mutton leads to the release of acetylcholine,resulting in salivation.The above descriptions are two examples of stimulus-response pathways thatoccur in humans, who are multicellular organisms with complex sensory organs.Now imagine a single cell. While cells clearly do not have the sensory organs ofhumans, we can draw analogies between the two. The cell still needs to sense its en-vironment, and it does so using receptors located on its surface. Certain molecules(the stimulus) bind to these cell surface receptors and trigger the induction of a sig-nalling pathway, similar to our neurochemical pathways. The signal is transmittedthrough the cell, eliciting a physiological response [67].An important physiological response for cells is polarization. During polariza-1Stimulus binds toreceptor PolarizationoccursFigure 1.1: This figure (not to scale) shows a simple depiction of cell po-larization. The cell begins in a homogeneous state, represented here bythe solid light purple colour. A chemoattractant (red dot) binds to a rep-resentative receptor on the surface of the cell, triggering a complicatedsignal transduction pathway inside the cell. This pathway eventuallyleads to polarization, in which certain proteins are found in high con-centration at either the front (direction of gradient) or rear of the cell. Inreality, there may be millions of chemoattractant molecules in the envi-ronment, and the cell responds in the direction of highest concentration.tion, proteins and lipids are redistributed within the cell to form distinct regions[7]. It is the first step to several important processes such as cellular division (cy-tokinesis, bacterial division, budding in yeast), and cell migration. For example,chemotaxis is directed cell migration along a chemical gradient; chemoattractantsare the chemicals that form this gradient and serve to attract the cells. Some leuko-cytes (white blood cells) migrate along chemoattractant gradients to find bacteriain our bodies. See Figure 1.1 for a simplified depiction of cell polarization in re-sponse to a chemoattractant. During migration, the front of the cell expresses thestructural protein actin and the back is rich in myosin proteins. In cells undergoingchemotaxis, a chemoattractant is the stimulus that induces polarization [56, 67].A cell undergoes both chemical and mechanical polarization, and both are ac-complished by the redistribution of either signalling or structural proteins, respec-tively. Redistribution of signalling proteins (chemical polarization) typically leadsto mechanical polarization through downstream effects to the structural proteins.Chemical polarization occurs much faster than changes to cell shape. Interest-ingly, polarization can be induced by different types of stimuli, or it can proceed2spontaneously [7]. The mathematical models are mostly partial differential equa-tions (PDES) that represent changes of protein distributions in both time and space.It is well known that PDES can be challenging to analyze.Even though cell polarization has been widely studied, the exact mechanismsby which it occurs are still unknown [53]. A fundamental question that modellersaim to answer is how this redistribution actually occurs. How does a cell innatelyorganize certain internal components into an asymmetric pattern? Which pieces ofthe signalling pathways are responsible and how do they re-organize? How com-parable are different cell types? Mathematical models attempting to answer someof these questions are abundant, and there is a wealth of competing models in theliterature [9]. Predictions made by these models are now being tested, and newexperiments are leading to discoveries of the mechanisms, modules, and pathwaysthat mediate polarization [24, 53]. It is important to understand these existing mod-els as best we can, in order to continue connecting mathematical and experimentalresults. This thesis provides a mathematical method that helps with the analysisof a specific class of PDE models, reaction-diffusion equations (RDES), and withanswering some important questions concerning polarization modelling.A review article by Jilkine and Edelstein-Keshet [29] uses numerical simu-lations to compare four different cell polarization models. Using a single set ofparameter values (either those given in the original publication, or slightly modi-fied as described in the article), the authors subject the models to different types ofstimuli. These stimuli are applied as perturbations to one equation of the model,and include single or multiple localized stimuli, a persistent graded stimulus (whichvaries in strength, or reverses direction), and noisy initial conditions. They eveninvestigate changes to the cell size. The models are compared according to thequalitative features of their simulations. This analysis is a step towards connectingsome of the existing models, but it is time consuming to carry out such explorationsfor many models and only reveals information for one parameter set at a time.In this thesis, I use a new method called local perturbation analysis (LPA) toaid in investigating a set of polarization models. LPA takes advantage of the factthat many cell polarization models, as we will see later, exhibit a disparity in dif-fusion rates between their components. Mathematical details will be presented inChapter 4, but ultimately we will remove the spatial component and use a system3of ordinary differential equations (ODES) to approximate the dynamics of a pulseapplied to the full reaction-diffusion system. ODES are much simpler to analyzeusing readily available bifurcation programs. Through the nature of LPA, regionsexhibiting threshold behaviour, stability, and even Turing instabilities are mappedout for any chosen parameter. This method allows us to quickly see the patterningregimes present, along with their robustness.1.1 Thesis outlineI will begin by introducing the biological and mathematical details that are relevantto the work in this thesis. In Chapter 2, I will present key components of a pathwayfor actin, a structural protein with many regulators, and survey how polarizationis realized in six different cell types. Chapter 3 will feature an introduction tomathematical modelling with RDES in the context of cell polarization. I will alsouse this chapter to provide descriptions of prominent cell polarization models inthe literature, which I will subsequently analyze using LPA.There are two goals that I wish to accomplish with this thesis. The first is tooutline and describe the method of LPA. In Chapter 4, I will present the LPA methodin a theoretical context. I will use two models, one of which is an extension ofthe other, to demonstrate the use and highlight the features of LPA in Chapter 5.The second goal is to analyze and compare a set of cell polarization models usinginformation gained via LPA and PDE simulations, while exploring how the benefitsand caveats of LPA help or hinder the analysis. In Chapter 6, I will present themathematical details of the models and their analysis using LPA, Turing analysis (intwo cases), and PDE simulations. Interpretations of the analysis and comparisonsof the models will be discussed in Chapter 7. In Chapter 8 I will document myconclusions for LPA and the polarization models.4Chapter 2Biological backgroundStudies in cellular biology involve complicated pathways of signalling proteins andlipids that perform different functions as messengers within the cell. Because thepathways are so complicated, and discussing many of the components is beyond thescope of this thesis, I will present a simplified view of one network: the structuralprotein actin and a few of its regulators. External stimuli trigger the activation ofsignalling pathways that result in actin redistribution and serve to mechanicallypolarize the cell. Various physiological responses result, and I will examine howpolarization is realized in six different cell types.2.1 Actin’ up: key components in an intracellularpathwayA network of filaments called the cytoskeleton is responsible for giving shape andstructure to the cell. It is composed of different kinds of filaments, but I will focuson actin filaments for this thesis. Unlike our static skeletons, the cytoskeleton isdynamic because its components can be broken down (depolymerized) and assem-bled (polymerized) in different parts of the cell as they are needed [17]. As such, italso functions in cell motility, contraction during cell division, and transportationof vesicles and organelles within the cell. Actin has two forms: a monomeric formcalled globular, or G-actin, and a filamentous form called F-actin. G-actin poly-merizes into F-actin, and filaments are depolymerized back into monomers. This5dynamic remodelling occurs continuously, and is typically directed by signallingproteins in response to stimuli that the cell receives from the environment.Figure 2.1 is a schematic depiction of the actin network whose details aredescribed below. The context under which we will discuss the actin network ischemotaxis. Receiving a stimulus initiates signal transduction pathways leading tothe growth of actin filaments in a specific location. Elongation of these filamentscreates an outward force at the cell’s leading edge, resulting in forward protru-sion of the cell [58]. To complete the forward motion, a network of actomyosinfilaments at the rear of the cell creates contraction, pulling the back of the cellforward. Actomyosin is a combination of actin filaments and the motor proteinmyosin, which binds to two actin filaments and uses energy to slide the filamentstogether. Protrusion at the front coupled with contraction at the rear drives the cellforward, resulting in migration [40].Polymerization is the process by which G-actin monomers combine to formactin filaments, and nucleation is the event that initiates polymerization. Nucleationis often initiated by other proteins, because spontaneous formation of the actinnucleus by the monomers themselves occurs very slowly. The nucleating proteincomplex called Arp2/3 binds to the side of a pre-existing filament of actin, andprovides a nucleation core for G-actin to bind [79]. Polymerization results anda new actin branch grows out from the original filament. Because the G-actinmonomers are asymmetrical, growth of F-actin is unidirectional. The filament ispreferentially disassembled at the minus end, and G-actin monomers are recycledto the fast growing leading edge, or plus end [31]. Capping proteins can stabilizeF-actin by binding to the plus end and preventing further polymerization.Rearrangement of actin is a physiological response the cell undergoes after re-ceiving a stimulus. Just like salivation prepares one for digestion, rearranging thestructure of the cytoskeleton prepares the cell to achieve a certain goal such as di-vision or migration [12, 18]. And just as salivation is triggered by neurochemicals,a signal transduction pathway is required to initiate actin rearrangement. The Rhofamily of GTPases is an important set of proteins responsible for the regulation ofseveral aspects of actin dynamics. Proteins such as Rac, RhoA, and Cdc42 – allmembers of the Rho GTPase family – target proteins such as Arp2/3 to stimulatenucleation of actin [37, 60].6Figure 2.1: A schematic depiction of the actin network. Actin nucleationand dissociation form a continuous circle with many different regulatorsmodifying the speed of each step. The following descriptions are addi-tional components not outlined in the text, see the text for details of thepathway. PIP2, which is recruited along with activated Rho GTPases,aids in activation of the WASp/Scar proteins, which then activate theArp2/3 complex. The nucleated actin monomers undergo ATP hydroly-sis to become ADP-actin, and ADF/Cofilin disassemble these older re-gions of the actin filament into monomers. Profilin aids in the turnoverof actin monomers. c©2006 Jeffrey R Kuhn, adapted with permission(originally appeared in [58]).GTPases often have two forms, giving them the ability to act as switches thatare either on or off [3, 4]. The active form binds GTP (guanosine triphosphate).Hydrolysis of GTP into GDP (guanosine diphosphate) releases energy to drive areaction. A GTPase bound to GDP is in the inactive conformation. In the Rhofamily, transitions between these two conformations are regulated by an additionalset of proteins, called GEFs and GAPs. GEFs, or guanine nucleotide exchangefactors, are catalysts for the exchange of GDP for GTP. This reaction puts theRho GTPase into its active conformation, and so GEFs are activators. Conversely,7Figure 2.2: Schematic diagram of Rho GTPase cycling between active andinactive forms. The activated form has an exposed hydrophobic tailwhich is inserted into the cell membrane. While bound to the mem-brane, Rho GTPases diffuse slowly. The inactive form is found withinthe cytosol where diffusion is quicker. The GEF and GAP proteins serveto regulate the activation and inactivation, respectively, of the Rho GT-Pases.GAPs (GTPase-activating proteins) induce the hydrolysis of GTP into GDP, andare therefore inactivators of Rho GTPases [10].The Rho GTPase undergoes conformational changes as it switches betweenthe active and inactive forms. As a result, the two forms have different locationswithin the cell. In its active conformation, the GTPase has an exposed hydrophobictail which anchors the GTPase into the cell membrane. Here, the active GTPasediffuses slowly on the membrane but does not move into the interior of the cell.When a cell is polarized, there is a high concentration of the active GTPase locatedat a certain patch on the membrane. The inactive form diffuses within the aqueousinterior of the cell, known as the cytosol. Note that proteins in the cytosol diffusemuch faster than those bound to the cell membrane. See Figure 2.2 for a schematicdepiction of Rho GTPase cycling.2.1.1 SummaryLet us summarize this information from the top down, and look again to Figure 2.1for a visual depiction of the network. The first step to cell polarization is typicallya stimulus. However, some cell types will randomly initiate polarization. When8specific molecules bind to cell surface receptors, a signal transduction pathway isinitiated. Many proteins are involved, and we have mentioned only a few key com-ponents for this response. Downstream in the pathway, the activities of GEFs andGAPs are modified to create distinct regions of high and low activity of the RhoGTPases (RhoA, Rac, or Cdc42 for example). Rho GTPases affect the activityof the protein complex Arp2/3 through the downstream action of different targetproteins, and Arp2/3 induces the nucleation of actin. G-actin nucleation is upreg-ulated in certain regions, and F-actin filaments form. This rearrangement of theactin cytoskeleton helps the cell respond to the stimulus it received.2.2 Polarization in different cell typesNow that we have become familiar with a pathway leading to cytoskeletal reorgani-zation, let us examine some examples of polarization in actual cells. The followingsystems are mentioned specifically because they are well studied, and are depictedin the mathematical models discussed next in Chapter 3. Figure 2.3 contains mi-crographs of each of these cells for visual reference. They form a diverse set ofsystems that simultaneously highlight the similarities and the differences in po-larization among different cell types. Keep in mind the types of stimuli, internalsignals and regulators, and the resulting responses as points of comparison. Mostof the cells considered in this thesis are eukaryotic cells, which are larger cellswith a defined nucleus and other complex structures. A distinguishing feature ofeukaryotic cells is their ability to sense spatial gradients as low as 1-2% across theirdiameter [24, 29]. In contrast, bacterial cells are small and can only sense temporalchanges in chemical gradients.Neutrophils are white blood cells that are part of the innate immune system, theinitial response to infection. Their main functions are to ingest bacteria and re-lease antimicrobial proteins to the surrounding environment. Chemotaxis is a keyphysiological response that neutrophils must undergo in order to locate bacteriaand sites of inflammation. They do so by migrating up a chemoattractant gradient,in this case formylated peptides which are by-products released by bacteria [77].Neutrophils also migrate randomly in a uniform chemoattractant gradient [2].9a) b) c) d) e) f) Figure 2.3: Micrographs of the six cell types that are described in the text: (a)neutrophil (larger cell), a type of white blood cell, chasing a bacterium(small black dot); (b) the slime mold Dictyostelium discoideum [19];(c) keratocyte, a skin cell from fish scales [30]; (d) neuron [72]; (e)Saccharomyces cerevisiae, a type of yeast [44]; and (f) the bacteriumEscherichia coli. See Appendix A for citations of (a) and (f). Note thatthese cells are not shown to scale in comparison with each other. A celldiameter is typically: (a) 8-10 µm, (b) 10-20 µm, (c) 30-40 µm, (d)4-100 µm, (e) 5 µm, and (f) 2 µm.Neutrophils respond with high sensitivity to gradients of chemoattractant throughasymmetric activation of different signalling pathways. This response forms mu-tually exclusive front and back regions [82]. The frontal zone becomes protrud-ing pseudopods pointing toward the source of infection. Exclusive to this regionis increased activation of Rac and production of F-actin [82]. The back regionforms a contracting appendage called a uropod. Activated Rho, Rho-dependentkinases (ROCK, regulators of F-actin organization and contraction) and myosinII (promotes contraction of actin filaments) are localized here [82]. This spatialdistribution of front and back responses occurs within minutes [80].10Dictyostelium discoideum is a social, soil-dwelling amoeba. D. discoideum un-dergoes several stages in its lifecycle, beginning as unicellular amoeba, coalescinginto a multicellular slug, and then into a fruiting body which releases reproductivespores. During the unicellular amoeboid stage D. discoideum find and consumebacteria using folic acid produced by the bacteria as a chemoattractant. If no foodis present, starvation triggers the cells to produce cyclic AMP (cAMP) to attractneighbouring cells. The unicellular amoeba stick together and aggregate to formthe next stage in their life cycle.D. discoideum move along chemical gradients via cellular extensions calledpseudopods (Latin for “false foot”). There are two mechanisms by which pseu-dopods form. New pseudopods may be created or, more commonly, they mayform when a pre-existing pseudopod splits. A dominant pseudopod, biased by thedirection of a chemoattractant, takes over and the other pseudopod retracts [1, 26].Cells polarize in response to folic acid or cAMP, but may also spontaneously polar-ize in the absence of external stimuli [39]. In the latter case, D. discoideum move ina random direction via zig-zagged movements, using this motion to increase theirchances of finding food [34]. It is suggested that pseudopod formation is inherentto the cell and is biased by, but does not require, external signals, making sponta-neous motion possible. Polarization patterns are mediated by F-actin (front), PI3K(front), and PTEN (back) [39]. The latter two proteins affect the level of PIP3 (seeFigure 2.1) at the leading edge versus the rear of the cell.Keratocytes are epithelial (skin) cells that are isolated from fish scales. In con-trast to neutrophils and D. discoideum which migrate using pseudopods, kerato-cytes have a large flat structure at the cell front called a lamellipod (Latin for “sheet-like foot”) [83]. The structure and shape of the lamellipod is preserved as the cellmoves, and results in a rapid gliding motion. Another feature contrasting the previ-ous two migratory cells, keratocytes are not known to respond to chemoattractantsor other chemical stimuli. Instead, polarization is initiated through a mechanicalstimulus and is maintained long after that stimulus has been removed [83]. Kera-tocytes also respond to electrical stimuli.In the absence of a stimulus, stationary keratocytes are capable of initiating mi-gration by polarizing spontaneously after approximately 30 minutes. It is also pos-11sible to create fragments of the lamellipod which contain no nuclei, microtubules(another type of cytoskeletal filament), or organelles. After applying a mechanicalstimulus, such as a stream of fluid from a micropipette, these fragments still exhibitrapid and persistent migration, indicating that the necessary components for polar-ization are already present in the cellular fluid [73]. Symmetry breaking is thoughtto initiate at the rear of the cell through the actomyosin network, as opposed toprotruding structures at the front [83].Neurons, or nerve cells, are animal cells that transmit information via a seriesof chemical and electrical signals. The neuron is composed of a cell body, manysmall dendrites which receive stimuli from other neurons, and a large structurecalled an axon which transmits the signal. Proper development and placement ofthese structures are crucial for the organism, and asymmetric distribution of cel-lular regulators plays a key role. Polarization is induced through environmentalcues, but spontaneous polarization in the absence of a signal is also possible. Intra-cellular regulators include the Rho GTPases among others. Once established, thecytoskeleton maintains polarity for the lifetime of the neuron [68].The neuron begins its development as a round cell with actin-rich growth conesfound on the ends of small, finger-like appendages that eventually become the axonand dendrites. The growth cones support extension of the appendages by sensingand responding to spatial signals. While neurons are morphologically very differ-ent from motile cells, zooming in to a growth cone reveals a picture very similarto neutrophils or keratocytes [71]. Their movement is dictated by a lamellipodand filopodia (hair-like projections) which are guided by environmental stimuli. Asymmetry breaking event eventually occurs to distinguish the axon and dendrites,and different cues in the environment serve as chemoattractants or repellents forthese navigating structures. For example, a chemoattractant called semaphorin 3Aattracts dendrites, but is a repellent for the axon [68]. These cues ensure the newneural structures are oriented correctly with the surrounding neurons to properlytransmit signals.12Saccharomyces cerevisiae is an extensively studied type of budding yeast usedboth in research as a model organism, and in food chemistry through winemaking,baking, and brewing. S. cerevisiae cells polarize for mating and budding, a type ofasexual reproduction where a bud forms and pinches off the parent cell to becomea new daughter cell. Two overlapping polarization pathways have been suggested:one involving the actin cytoskeleton and the GTPase Cdc42, and another involvingthe Cdc42 GTPase signalling network [64]. Polarization still occurs in experimentswhen actin or Bem1 (a protein within the Cdc42 signalling network) are inhibited,but fails if both are inhibited [76].During mating, pheromones are the stimuli for Cdc42 polarization, and areanalogous to chemoattractants [57]. The yeast senses a pheromone gradient, andextends the tip of a mating projection in the direction of the source. During bud-ding, alignment of the polarization axis is a result of positional landmarks on thecell surface. These landmarks guide the subsequent bud formation. It has beenshown that cells lacking these landmarks are still able to form buds normally, butdo so in a random direction [27, 55].Escherichia coli is a rod shaped bacterium which comprises part of the gut mi-croflora population of many animals. I mentioned previously that most of the mod-els studied in this thesis are based on eukaryotic cells; a model for E. coli is theexception. In this case, I am considering the role polarization plays in the cellulardivision of E. coli, and not as a response to chemical gradients. Bacteria are toosmall to sense a spatial gradient across their diameter and as such only respond totemporal changes in chemical concentrations.During cell division, the division plane must be accurately positioned in orderto ensure two equal sized daughter cells are produced. In bacterial cells, the Zring, composed of the protein FtsZ, determines the location for the division plane[36, 61]. Forming an asymmetric distribution (polarization) of FtsZ is the first stepin cell division, and occurs after the genome has replicated. The Min system inE. coli is composed of the proteins MinC (inhibits FtsZ polymerization), MinD,and MinE (both regulate MinC). These proteins are responsible for localizing theZ ring to the center of the cell [42]. The Min system oscillates from pole to poleof the cell, causing FtsZ to localize at the center. Contraction of the Z ring leads to13Cell Type Stimulus IntracellularComponentsResponseNeutrophils Formylated peptides(from bacteria)Rac, F-actin, Rho,ROCK, myosin IIChemotaxis viapseudopodsDictyosteliumdiscoideumFolic acid (bacteria),cAMP (starved cells)PI3K, PTEN, F-actin Chemotaxis viapseudopodsKeratocytes Mechanical changes Actomyosin network Motility vialamellipodNeurons External chemicalgradients (can bedendrite/axonspecific), adhesionmoleculesRho GTPases, PI3K,cofilinChemotaxis vialamellipod andfilopodiaSaccharomycescerevisiaePheromone gradient (1) Actin and Cdc42;(2) Cdc42 signallingnetworkMatingprojection,buddingEscherichiacoliEnd of genomereplicationFtsZ, Z-ring, MinSystemCell divisionTable 2.1: Summary of polarization in the six different cell types. Thecolumns describe the stimulus, some intracellular components involvedin the polarization response to that stimulus, and the morphological re-sponse. See the text in Section 2.2 for additional details. Note that thesedescriptions are not exhaustive, as the intracellular pathways are verycomplex.division of the cell.2.2.1 SummaryTable 2.1 provides a summary of the six different cell types discussed above. Ineach case, the cell receives a type of stimulus, has intracellular componentswhich polarize or help to regulate polarization, and elicits a response.The actual intracellular pathways involved are much more complex than I havedescribed here. It is often known how components and regulators of the actin14network act locally to influence the cytoskeleton, but the mechanisms by whichthese signalling molecules and regulators coordinate to organize the actin networkis not fully understood. We can use qualitative features that are observed acrossdifferent cell types as points of comparison. When creating a mathematical model,the same qualitative comparisons can be made between simulations and biology,keeping in mind what type of stimulus the cell receives and the morphologicalresponse it elicits.15Chapter 3Mathematical modellingbackgroundThe goals of mathematical modelling in cellular biology are diverse. Some modelsseek to replicate the complexity of a biological system, while others isolate onlythe key components and features. The level of realism in a model ranges frompurely conceptual to fully mechanistic; models may focus on a few interactingcomponents or include an entire network [48]. Mathematical modelling can pro-vide predictions for experimental outcomes, or reveal the underlying mechanismsin biological systems that are already well observed. With the broad scope of mod-elling styles available for computational biology, both the type of biological systemand the information one might hope to gain from research will influence the type ofmodel that is used or developed. The models I analyze for this thesis are diverse intheir approaches. I do not extend beyond three-dimensional models (i.e., systemswith three equations) for my analysis, but some of the models are simplified fromdetailed models with many components.3.1 Introduction to modelling with reaction-diffusionequationsDifferential equations represent how certain dependent variables (i.e., concentra-tion) change as independent variables (i.e., time, space) change; ordinary differen-16tial equations (ODES) involve one independent variable, while partial differentialequations (PDES) involve two or more independent variables. Differential equa-tions can be combined into systems, which then take into account multiple depen-dent variables and how they affect each other. In this thesis I will consider modelsthat form a subset of PDES: reaction-diffusion equations (RDES). A key featureof biological systems is their ability to inherently develop patterns, and RDES arecapable of a wide range of behaviours (some of which are discussed later in Sec-tion 3.2). RDES model the change in concentration of a particular quantity overtime and space, and include reaction kinetics and the spatial effects of diffusion[65].Diffusion is a method of passive transport of chemical species, meaning it doesnot require cellular energy. It can also be thought of as a “smoothing” process;individual molecules still move randomly, but on average species move from areasof high concentration to lower concentration. Diffusion leads to a uniform distri-bution provided there are no other influences. The PDE for diffusion in 1D is:∂u(x, t)∂ t =∂∂x[D(u,x)∂u(x, t)∂x]. (3.1)The variable u(x, t) represents the density or concentration of the quantity undergo-ing diffusion, at a particular time t and location x. The diffusion coefficient D(u,x)is the rate of diffusion, dependent on the density u and location x. Throughout thisthesis, I will assume D(u,x) is constant, i.e., D(u,x) ≡ D. This assumption leadsto the simplified equation:∂u∂ t = D∂ 2u∂x2 . (3.2)I will use Equation 3.2 to model the diffusion of cellular signalling components.The diffusion coefficient D represents the rate at which diffusion occurs fora particular chemical species. Depending on their size or the medium, chemicalspecies diffuse at different rates. In cells, there is a distinct disparity in diffusionrates between species diffusing in the cytosol or on the cell membrane. This dis-parity is due to the high viscosity and large number of obstacles encountered onthe cell membrane, as compared with the cytosol.17In Equation 3.2, a reaction component can be included by appending a func-tion f (u). This function, also termed the reaction kinetics, models the creation orconsumption of a chemical species in a given reaction. With a system of more thanone RDE, interconversion terms may represent conversion from one substance intoanother.The general form for a system of two RDES is given by:∂u∂ t = Du∂ 2u∂x2 + f (u,v), (3.3a)∂v∂ t = Dv∂ 2v∂x2 +g(u,v). (3.3b)There are two chemical species, u and v, with distinct diffusion rates (Du and Dv)and reaction kinetics ( f (u,v) and g(u,v)). Specific models will be explored thatbuild on the structure of Equation 3.3.In the context of cell polarization, RDES are a natural way to model the sig-nalling components [29, 65]. Take for example the Rho GTPases. They are smallmolecules which can be modelled by Equation 3.3. Let us define u as the activeform and v as the inactive form. The former diffuses slowly on the cell membranewhile the latter diffuses quickly within the cytosol, so the diffusion coefficients willexhibit the property that Du Dv.The GTPases can also react within their environment, represented mathemati-cally by the terms f (u,v) or g(u,v). Within the “reaction” part exists interconver-sion between the active form and the inactive form of the GTPase. The activationrate represents conversion from v into u, and therefore is dependent on the concen-tration of v. In an autocatalytic reaction, positive feedback exists from the activatorto itself. Therefore, the activation rate will also increase with the concentration ofu. Our activation term might have the form kA(u) · v, and would be positive in theequation for u but negative in the equation for v. The inactivation rate representsconversion from u into v, and will be dependent on the concentration of u. Thisterm could have the form kI · u. Biologically, GEFs and GAPs catalyze activationand inactivation of the GTPase, and in some papers modellers represent them as18parameters influencing the reaction rates. Production and decay of the GTPasecould also be included in the reaction kinetics. The fact that polarization occursfast relative to production and decay leads some modellers to ignore the latter pro-cesses.3.2 Patterns in reaction-diffusion systemsFrom microscopic organization of proteins in cells to animal coat patterns visibleon the macroscopic scale, patterns are ubiquitous in nature. Pattern formation inchemical systems is the self organization of two or more species into a nonuniformspatial distribution. The changes in concentration of these species are typicallydescribed through the use of RDES. In this section, I will explore the major typesof pattern forming mechanisms behind cell polarization models. These are Turingpatterns, waves, and the local excitation / global inhibition (LEGI) mechanism.Interestingly, none of these mechanisms on their own are sufficient to reproduceall features of cell polarity. The reaction diagrams in Figure 3.1 show examples ofthe three mechanisms.3.2.1 Turing patternsThe idea behind Turing models existed before it was applied to cell polarization.In fact, Turing patterns were suggested as a mechanism for morphogenesis, whichis pattern formation on the multicellular scale. Alan Turing proposed that undercertain conditions, a system of chemical interactions that is stable can becomeunstable to small perturbations in the presence of diffusion [70]. The diffusionrates must be different, and there are additional constraints on the parameters anddiffusion coefficients (see Section 4.1 and Appendix C for more detail). In a twocomponent system, an additional restriction for Turing patterns applies: the systemmust follow either an activator-inhibitor or a substrate-depletion scheme.Gierer and Meinhardt [14] were the first to simulate pattern formation by shortrange activation with long range inhibition. Note that the authors were initiallyunaware of Turing patterns and came upon this idea independently. Two speciesmake up this activator-inhibitor model: a slowly diffusing, membrane-bound ac-tivator with positive feedback, and an inhibitor. The global inhibitor, capable of19long-range diffusion since it moves through the cytosol, serves to inhibit the ac-tions of the activator. Localized peaks of activation, i.e., polarization, arise fromthe interplay between these two species (see Figure 3.1(a)). Substrate depletionhas been proposed as an alternate mechanism for generating Turing patterns in cellpolarization. Here, the concept of the inhibitor is replaced with that of a cytosolicsubstrate. The substrate is depleted as it becomes activated and associates with themembrane, leading to a reduced rate of activation [53].In general, Turing models replicate spontaneous polarization and display a dis-tinct internal amplification of the gradient. In a Turing regime, the modelled cellis sensitized, since even minute perturbations lead to polarization. Turing modelscannot account for a stable resting state in which a cell remains unpolarized. Thepatterns generated through Turing instabilities are usually frozen, which allowsthese models to maintain polarization once a stimulus is removed. As a caveat,the pattern does not respond to additional stimulation in different directions. Cer-tain models or parameter regimes may lead to the existence of multiple peaks ofactivation.3.2.2 Travelling wavesExcitable and wave-based models also existed prior to cell polarization modelling,but this approach was not applied until several years after Gierer and Meinhardt’smodel. The idea that polarization occurs through waves of activity was sparkedby experimental data. Waves or pulses of activation have been experimentally ob-served in different cell systems, such as actin polymerization [47, 74, 75, 78] or itsregulators [5, 13, 81]. Travelling wave solutions are appealing, as a suitably largeperturbation is required to destabilize a rest state, and the information travels fast.However, while the domain is polarized when the wave is present, in most systemsthe wave sweeps through the domain and leaves behind a uniform distribution [9].For a polar distribution to be stable, the wave or pulse must stop somewherein the domain, termed “wave pinning”. While stopping the wave would usuallyrequire an unjustified inhomogeneity in the domain, a model by Mori et al. [49]successfully explains how wave pinning could result from a depletion effect in auniform domain. For wave pinning to occur, a model requires several features:20bistability, different diffusion rates, and mass conservation. Nonlinear reactionterms lead to the existence of two stable steady states (bistability). The model be-gins at a resting state with a uniform distribution at the lower steady state. Thewave pinning mechanism works by initiating a wave that depletes its limited sup-ply of fuel and has to stop, leaving the cell with a stable polar distribution. Asufficiently large perturbation is required to initiate the wave, whose height is de-termined by the higher steady state. The final pattern is stable, even in the absenceof a perturbation. See Figure 3.1(b) for the wave pinning reaction diagram.Wave pinning models feature a strong degree of amplification, reach the polar-ized state quickly, and maintain the polar distribution. While Turing-type modelsare always sensitized to even the smallest perturbations, wave pinning models ac-count for resting states since a sufficiently large stimulus is required for the acti-vator to reach the second steady state. This fact also means these models do notspontaneously polarize.3.2.3 Local excitation / global inhibitionAn interesting feature of many chemotacting cells is the ability to adapt to dif-ferent levels of stimuli. When the stimulus increases homogeneously over space,the signalling pathways respond with transient polarization. Even when a highbackground concentration of chemoattractant is present, the cell only maintainspolarization when a gradient is present [24, 32]. Gradient sensing is the persistentpolarization response which occurs when a chemoattractant gradient is present. Po-larization featuring both adaptation and gradient sensing is what the LEGI modelsought to replicate. The first LEGI model was created by Levchenko and Iglesias[32].In contrast to models featuring Turing and wave patterns, LEGI models are notcreated with inherent patterning mechanisms in mind. They lack the necessarynonlinearities and spatial terms for forming patterns on their own [49]. Instead,a membrane bound activator and a cytosolic inhibitor are produced in proportionto a stimulus. Since the activator diffuses slowly, it accumulates at the front ofthe cell in the direction of the gradient, while the inhibitor takes on a uniformconcentration. A polarized distribution is established through these differences21uvuvRSuvS Sa) b) c) Figure 3.1: Examples of reaction diagrams leading to polarization by eachof the three mechanisms: (a) Turing instability, (b) wave pinning, and(c) local excitation / global inhibition (LEGI). The reacting species areshown in purple (gradient = membrane bound and slow diffusing, solid= cytosol and fast diffusing), and the stimulus is depicted in red. De-pending on the model, u is the activator, v is an inhibitor or the inactiveform of u, and R is a response concentrations by a response element, as long as the stimulus is present [53](see Figure 3.1(c)). When the stimulus is removed, the model does not maintainpolarization.LEGI models in general can account for features such as adaptation and sensi-tivity, and are robust to changes in parameters [32]. Since they lack inherent patternforming abilities, LEGI models on their own do not maintain polarity and cannotrandomly polarize, but could be coupled to other modules which would accountfor these features [25, 29]. LEGI models can even be coupled with each other toqualitatively account for observed directional sensing behaviours in cells [38].3.3 Models for cell polarizationThere are numerous models for cell polarization in the literature; here I survey aselection of six polarization models. All these models feature disparities in diffu-sion rates among their components. This section serves to introduce the differentmodels considered in this thesis, to provide the context in which the models were22developed and to describe some results observed in the original publications. Notethat for each of these models, the authors applied a stimulus to initiate pattern for-mation. The stimulus can be transient or persistent, and it is either superimposedon the “active form” equation or incorporated into the model as a spatially vary-ing parameter. Mathematical details of the models and stimuli will be discussed insubsequent chapters.Wave Pinning Model The Wave Pinning model is a two component RD systemintroduced by Mori et al. [49] for GTPase-based polarization. Originally, a sixcomponent model was developed for active and inactive pairs of three Rho GT-Pases: Cdc42, Rac, and Rho. In the simpler minimal model, only one pair isconsidered. Since polarization occurs on a short enough timescale that proteinsynthesis or decay can be considered negligible, the developers chose a model inwhich mass is conserved. It is based on interconversion between the two differentforms (active, inactive), with positive feedback leading to its own activation. Themodel also features bistability, in which there are two stable steady states. Thismodel is meant to describe a set of Rho GTPases in a motile eukaryotic cell, suchas a neutrophil, but is not specific to one cell type.Several different types of stimuli were applied to the active form; these includea transient localized stimulus, a persistent gradient, and large amplitude noise. Un-der certain conditions, the model responds with a stable polar distribution of high“active” concentration at one end and low concentration at the other. After the stim-ulus is applied, the combination of bistability, unequal diffusion rates, and conser-vation result in a wave front that propagates into the domain and stalls (hence theterm “wave pinning”). The polarization response is stable even after the stimulusis removed, and remains sensitive and can repolarize to stimuli in other directions,provided they are sufficiently large. It was noted that multiple activation peaksmay form when the model is stimulated in two directions simultaneously. Thepeaks will eventually resolve into one peak, but depending on the distance betweenthe peaks, that process generally occurs too slowly to be biochemically relevant[49]. Additionally, the model does not account for cells that continually send outpseudopods (which correspond to continually growing and decaying peaks of acti-vation).23Otsuji Model The model presented by Otsuji et al. [54] is a conceptual, simpli-fied model following a more complicated model for Rho GTPases developed in[62]. In the original paper, the authors developed a model of six reaction-diffusionequations for the membrane bound and cytosolic forms of Rac, Cdc42, and RhoA.It was found that this model exhibited the following features of cell polarization:a switch-like accumulation of Cdc42 and Rac at the front and RhoA at the back;emergence of a single stable peak if multiple peaks are present; and the ability ofthe polarized peak to sense a new stimulus. Similar features were also observedin models by Subramanian and Narang [50, 66], even though these models weredeveloped for different molecular species.The authors hypothesized that the two different models contained similar un-derlying features, and found that both models are mass conserved, RD systemswith diffusion-driven instability (Turing patterns). To investigate whether thesefeatures are sufficient for cell polarization responses, a simpler conceptual modelwas developed. This conceptual model retained cell polarization features similar tothose observed in the more complicated Rho GTPase model. The authors proposedthat a reaction diffusion system featuring mass conservation and diffusion driveninstability may be a principle of cell polarity in migrating cells.Goryachev Model The two component Goryachev system models the polariza-tion inherent in yeast cell budding [15]. The first phase in budding involves theformation of a cluster of Cdc42 (polarization) at the bud scar, and is independentof actin. Actin is recruited to continue the process of budding in the second phase.This model focuses on the first phase, and was initially comprised of eight variablesfor the cycling of Cdc42 between its membrane-bound active form and cytosolicinactive form, its effector Bem1, and the GEF Cdc24. In order to elucidate funda-mental features of the system, the model was reduced to a two component systemfor the active/inactive forms of Cdc42. Both the complete and the simplified ver-sions of the Goryachev model feature conservation, based on the principle that budformation happens on too short a timescale for production or decay of any of thecomponents. This model was developed as a detailed biochemical model: it takesinto account known molecular components and experimental findings.A Turing instability is the core pattern forming mechanism of the Goryachev24model. The reduced model can be abstracted as an activator substrate-depletionmodel, where active Cdc42 functions as the activator and it depletes inactive Cdc42,the substrate. Positive feedback to the activation of Cdc42 proceeds through twopathways which are incorporated into the model. The authors find that this re-duced model is sufficient for the emergence of polarization in the form of a singlepeak of activation. The response is robust to noise in the stimulus, and changesin parameter values. This model marks an interesting case in which a Turing-typemechanism emerges from known, experimentally determined molecular networks.The authors propose that Rho GTPases continually expend energy to maintain a po-larized structure and combat diffusion, so membrane cycling is required for GTPhydrolysis. Competition for conserved resources leads to the formation of a uniquepeak.Meinhardt’s Activator-Substrate Depletion Model (ASDM) This model, presentedby Meinhardt and de Boer in [46], is unique from the other models consideredhere as it was developed for bacteria, rather than eukaryotic cells. Despite thedifferences in biological motivation, the RDES for this model are very similar tothe others, as we will see in Chapter 6. The Meinhardt model seeks to simulatethe accurate placement of the Z ring in Escherichia coli cells. The first step inpositioning the ring involves polarization of FtsZ. The model then simulates poleto pole oscillations of the MinD and MinE proteins, leading to localization of the Zring (represented here by membrane-bound FtsZ) in the center of the domain. Thefull model includes a total of six equations for membrane bound and free formsof the proteins FtsZ, MinD, and MinE. However, since FtsZ is responsible for theinitial polarization event in this model, in subsequent analysis I will focus only onRDES for FtsZ.This model is based on an activator-substrate depletion mechanism, and ex-hibits Turing-type instabilities. The cytosolic form of FtsZ acts as a substratewhich is depleted as FtsZ localizes to the membrane. Note that this model doesnot employ conservation, and it is possible for multiple peaks of active FtsZ toform and persist. In the context of the full model however, any peaks that do formare eventually resolved into a single peak in the center of the domain by MinDand MinE. Simulations by the authors reveal that the full system is self-organizing,25since patterning can be generated through random fluctuations alone.Pseudopod Centered Model This model was originally developed as a system ofODES with a discrete spatial component by Meinhardt [45], and was intended tomodel polarization in growth cones of neurons. Neilson et al. [52] then extendedthe model to a PDE system and applied it to pseudopod generation. I will refer tothis model as the “Pseudopod Centered” model, since it models chemotaxis cen-tered around the dynamic formation of pseudopods. The model represents an in-trinsic patterning system that orients the direction of pseudopod formation. Ratherthan assuming pseudopods are created in response to external asymmetries, thismodel assumes that pseudopods are continually extended and retracted. The direc-tion of the cell is modulated through the splitting of pseudopods that are already inexistence. An external gradient serves to bias which pseudopods become dominantand which are retracted.To achieve the effect of pseudopod formation and retraction, three componentsare present in the model: a local activator, a local inhibitor, and a global inhibitor.Together, the local activator and global inhibitor form the basis for Turing instabil-ities, leading to spontaneous peaks of activation. The local inhibitor destabilizesthese peaks, leading to splitting and retraction of old peaks while new peaks ofactivation form in their place. Ultimately the RDE solutions resemble the dynamicformation, splitting, and retraction of pseudopods. The Turing pattern that formsallows for more than one activation peak at a time, which is biologically relevantin this case since a cell may have multiple pseudopods. A superimposed stimulusbiases the location of dominant activation peaks, but random peaks will form in theabsence of a stimulus. When this biochemical model is linked to deformation of acell edge in 2D, simulations show a cell that performs U-turns in shallow gradients,and will form de novo pseudopods in response to changes in a steep chemoattrac-tant gradient, similarly to D. discoideum cells.Balanced Inactivation Model The Balanced Inactivation model was developedfor directional sensing seen in cells like D. discoideum amoeba, fibroblasts, andneutrophils, that respond to external chemical gradients with high sensitivity. Levine26et al. [33] base the development of their model on the following biological obser-vation: several components quickly and uniformly localize to the cell membrane inresponse to a stimulus, and shortly after, these components disassociate from theback of the cell [28, 63]. Due to this type of response, it was noted that subcellularreorganization is not merely an amplification of the external chemical gradient butappears as a switch. The Balanced Inactivation model is a LEGI-type model incor-porating this switch-like behaviour in terms of abstract components: a membrane-bound activator and an inhibitor with membrane-bound and cytosolic forms.Simulations reveal an initial localization of activator at the membrane, with aloss of localization at the back. Diffusion determines the timescale of activatorresponse at the back of the cell. The model responds to shallow gradients, po-larization can be reversed, and multiple points of stimulation lead to responses atthe sources. The switch-like behaviour is realized by the equal production of theactivator and inhibitor chemical species, the diffusion of the cytosolic form of theinhibitor throughout the cell, and the inhibitory effects of membrane-bound in-hibitor on the activator. A common criticism of LEGI type models is the lack ofbiological evidence for the inhibitor. The authors suggest that a heterotrimeric G-protein may play the role of all three chemical species described in this model, butso far the inhibitor has not been identified.3.3.1 SummaryA summary of the six models is provided in Table 3.1. It includes the type of cellthe model is attempting to capture, and the signalling components represented inthe equations. The model may be specific (the Meinhardt ASDM model targetsE.coli) or general (motile eukaryotic cells). Some models were developed for aparticular cell type, but can be extended to other cell types as well. An example isthe Pseudopod Centered model. It was designed for nerve cells and later applied toa chemotaxis model for D. discoideum. As the models range from mechanistic toconceptual, the signalling components and reaction kinetics vary from being bio-chemically based on experimentation (Goryachev) to abstract, with the biologicalcounterpart yet to be determined (Balanced Inactivation).Currently and with the exception of Balanced Inactivation, a Turing regime is27Model Cell type/s Signalling Components PatterningRegimesWave Pinning Motile eukaryoticcellsActive/inactive pair of RhoGTPasesTuring,ThresholdOtsuji Motile eukaryoticcellsActive/inactive pair of RhoGTPases (conceptual)TuringGoryachev Yeast cell budding Cdc42 network TuringMeinhardt ASDM E.coli cell division FtsZ, Min system TuringPseudopodCenteredNerve cells,D. discoideumAbstractactivator/inhibitors tripletTuringBalancedInactivationMotile eukaryoticcellsActivator/Inhibitor pair(conceptual, potentiallyheterotrimeric G-protein)NoneTable 3.1: A summary of the six cell polarization models presented in Sec-tion 3.3. The model may not be specific to one cell type, and the sig-nalling components represented by the RDES range from conceptual tospecific. The patterning regimes listed were discovered by the originalauthors of the model.featured in each model. The Wave Pinning model is the only set of RDES here fea-turing a threshold response. No known inherent patterning exists in the BalancedInactivation model, as it is an extension of the LEGI system. The Wave Pinning,Otsuji and Goryachev models all feature conservation, but the other three modelsdo not. Several models target eukaryotic cells undergoing chemotaxis. Are theymodelling the same aspect of the signalling pathway, and how do they compareagainst each other?3.4 Connecting models and biology3.4.1 Abstracting features of cell polarizationJilkine and Edelstein-Keshet [29] outline the following list of features for polariz-28ing cells in a review article on cell polarization models:1. Amplification2. Sensitivity3. Maintenance4. Spontaneous polarization5. Adaptation6. Single axis of polarity7. Pseudopod-splitting.Not every cell exhibits each feature, but different combinations of these featuresare observed across many different cell types. I have also highlighted these fea-tures for the different pattern forming mechanisms above in Section 3.2, and I willfrequently refer to these features throughout this thesis.During the polarization response, the external gradient is amplified internallythrough signalling pathways. The external stimulus gradient can be as shallow as a1-2% difference across the diameter of the cell, and amplification leads to macro-scopic asymmetry of signalling components within the cell. Some cells exhibitsensitivity to new stimuli. That is, if a cell is already polarized, it can sense stimuliin different directions and repolarize accordingly. This feature could be detrimen-tal in certain cell types like neurons; once polarization occurs it does not reorient.In cells featuring maintenance, polarity is maintained even if the initial stimulus isremoved. Some cells initiate polarization spontaneously with no external stimu-lus, while others require the stimulus to reach a certain threshold before polarizing.A cell could also have a resting state, in which it is stable to small levels of stimulusin the environment. Adaptation ensures that a cell responds to an actual gradient,and not a uniform change of a chemoattractant. When the concentration of stimu-lus increases uniformly, a cell might initially polarize in response to the temporalchange, but return to a non-polarized state shortly after. In cells that feature asingle axis of polarity, only one polarization front or zone is observed. Multiplestimuli do not result in multiple zones of activation. Pseudopod-splitting is only29Cell Type PolarizationFeaturesModel PolarizationFeaturesNeutrophils 1 2 3 4 6 Wave Pinning 1 2 3 4Dictyosteliumdiscoideum1 2 3 4 5 7 Otsuji 1 3 4 6Keratocytes 2 3 4 6 Goryachev 1 3 4 6Neurons 1 2 3 5 7 Meinhardt ASDM 1 3 6Saccharomycescerevisiae3 4 6 Pseudopod Centered 1 2 4 7Escherichia coli 3 4 6 Balanced Inactivation 1 2 5Table 3.2: A list of the qualitative polarization features that are observedin the cell types and models (presented in Section 2.2 and Section 3.3,respectively). The models and cell types are not intended to be comparedacross rows. Using qualitative features such as these, we can comparethe models to each other and to the cell types they represent.relevant in cells that move via pseudopods. Rather than forming new pseudopods,reorientation occurs through splitting of pre-existing pseudopods. The pseudopodsplits into a dominant side that extends while the other retracts.Table 3.2 lists the qualitative features observed for each of the cell types dis-cussed in Section 2.2 and the models presented in Section 3.3. These qualitativefeatures are useful points of comparison between mathematical models and biol-ogy, or between the models themselves. It is helpful to answer the question ofwhich polarization features each model accounts for, and if these features coincidewith experimental observations in the cell type being depicted.While these are commonly observed features, the cellular processes respon-sible for them may still be largely unknown. For example, neutrophils, D. dis-coideum, neurons, and keratocytes all exhibit (2) sensitivity to new stimuli and (3)maintenance of polarization. It is unknown how these cells are able to maintain po-larization, but still have the ability to sense stimuli and reorient their polarizationaxis accordingly.303.4.2 Polarization modelling questionsThe signalling network surrounding cell polarization is complicated. We want toidentify key or universal underlying mechanisms which would produce essentialfeatures of the polarization response. We can then examine the minimum modulesthat produce polarization patterns. Depending on the cell type being modelled,the desired polarization features may differ. Several questions arise surroundingthese models and many of them are linked to the qualitative features I have justdescribed.• Some cells polarize spontaneously while others require a stimulus to surpassa threshold value. Which models produce these responses?• Is maintenance of polarity observed in simulations?• How robust is the particular polarization response (pattern) to changes inparameter values? Cells and their environments are not identical, and so thepattern should be achievable over a robust set of parameters.• Is there a distinction in the underlying mechanisms for cells that move vialamellipods versus those that move via pseudopods? Recall that pseudopodsform on their own and are continually extending and retracting. Direction ofmovement is based on a dominant pseudopod rather than a constant “front”as in lamellipodial motion.• Does the model exhibit a unique axis of polarity or are multiple zones ofactivation observed?If a set of reaction diffusion equations produce the sought after responses of po-larization, they may indicate a corresponding module within the cell signallingpathway which is responsible for organizing the polarization response.31Chapter 4Local perturbation analysisIn general, partial differential equations must be analyzed on a case-by-case basis;there is no single method that can be applied to all partial differential equations(PDES). Historically, pattern formation in reaction-diffusion equations (RDES) hasbeen analyzed through numerical solutions, linear and nonlinear stability analysis.This chapter introduces a new method, called the local perturbation analysis (LPA),for analysing a class of RDES with slow and fast diffusing components. LPA ap-proximates the dynamics of a local pulse in a system of RDES using ordinary dif-ferential equations (ODES). One can then use commonly available bifurcation soft-ware packages to investigate regions in which the pulse grows (leading to patternformation) or decays (no patterning).For clarity, recall that RDES are a class PDES. The generic PDE system I amreferring to is a set of RDES of the form given in Equation Current methods of analysis for reaction-diffusionequationsIn contrast to ODES, there are few general methods for analyzing PDES, makingclassification among models of a similar theme difficult. Methods are often specificto classes of PDES, but even comparison between models within a particular classcan be challenging. A small amendment to one term (for example, modifying thereaction term in a RDE) may result in significant changes to the dynamics. These32changes may be predictable with further simulations and analysis, or may requirean entirely new method of analysis. Here we outline some existing methods foranalyzing RDES.4.1.1 Numerical solutionsOften the first method to investigate PDES is to examine their solutions analyticallyor numerically. If a PDE is linear, there exist methodical ways of obtaining ananalytical solution. Since most biologically inspired PDES are nonlinear, they areusually not analytically solvable. Numerical solutions can be used to visualize theultimate pattern that forms in each region.Numerical solvers are implemented in well-established ways, however notethat PDES do exist which require more tricks to solve numerically, such as thosewith sharp solutions or shocks [59]. Such PDES are beyond the scope of this thesis.The finite difference methods are a set of numerical methods that discretize thecontinuous domain of the equation to provide a numerical approximation for thePDE. For example, the finite difference method I use in this thesis is the Crank-Nicolson scheme:un+1i −uni∆t=D2(∆x)2((un+1i+1 −2un+1i +un+1i−1 )+(uni+1−2uni +uni−1))+ f (uni ), (4.1)where unk is the approximate solution at time t = n∆t and location x = k∆x. Here,D is the diffusion coefficient and f (u) is the reaction term. Refer to Appendix Bfor additional details on the Crank-Nicolson scheme.An important consideration for biological models is the robustness of a particu-lar type of solution or pattern. Numerical solutions can help to determine the effectsof different parameters on the system. However, once a numerical solver is imple-mented, adjusting parameter values to observe their effects on the PDE solution istime consuming even if the number of parameters is small. Complexity increasesfor models with large parameter sets. In Section 4.2, I will discuss a method whichcan be used to create maps of pattern forming regimes. Using this map to directparameter choices, numerical solutions become a powerful visualization tool.334.1.2 Turing analysisLinear stability analysis is a method for testing how the homogeneous steady state(HSS) of a RDE behaves when small, periodic perturbations are applied to the initialconditions. In the 1950’s, Alan Turing made an interesting observation concerningRDES [70]. If the HSS of a model is stable in the absence of diffusion, the addi-tion of diffusion may result in destabilization, leading to pattern formation. Thisphenomenon is referred to as a Turing instability. Turing instabilities only occur insystems of two or more RDES with different diffusion rates; additional constraintson the values for parameters and diffusion coefficients are also required. For a Tur-ing instability to occur in a 2x2 system, one of two schemes are required: (1) localactivation with global inhibition, or (2) substrate depletion. These two schemeswere described in terms of cell polarization models in Section 3.2.A Turing analysis maps the parameter space in which diffusion-driven instabil-ities occur. To briefly outline the method, we begin with a HSS which is stable tospatial perturbations in the absence of diffusion, uSS. After introducing diffusion,a small perturbation of the form:u′(x, t) = αeσtcos(qx), (4.2)is applied to the HSS. Here, α is the amplitude of the perturbation at t = 0 andis assumed to be small, σ is the eigenvalue of the linearized system, and q is thewavenumber. We want to determine which values of the wavenumber q result inan eigenvalue with Re(σ) > 0, that is, which wavenumbers will result in unstableperturbations. This question leads to a specific set of conditions for instability de-pendent on the parameters and diffusion coefficients of the original RDE system.We now have the parameter space and wavenumbers relevant for Turing instabili-ties. For a detailed derivation of Turing instabilities, see Appendix C.Note that the results of a Turing analysis are not valid once perturbations growtoo large, and so do not provide information on a long timescale. The ultimatepattern that forms must be determined through numerical simulations. Addition-ally, depending on the complexity of the reaction terms, it may be difficult to carryout a fully analytic Turing analysis, especially in cases where it is not possible tosolve for the HSS. Scaling and dimensional analysis can aid in Turing analysis and34exploration through simulations by reducing the number of parameters.4.1.3 Nonlinear stability analysisDepending on the type of PDE, methods of nonlinear stability analysis may revealmore information about the stability of a system. For example, asymptotic analysisconsiders the limiting behaviour of a system. It explores the effects as parametersor terms become very large or small, approaching the limiting cases of infinity orzero. Nonlinear methods are designed to treat PDES on a case-by-case basis, andcan be technically onerous to implement.4.2 Local perturbation analysisLPA was developed by Mare´e and Greineisen for the class of RDES in which thespecies involved exhibit a disparity in diffusion rates [16]. It was developed furtherin [21] and [22]. Recall the general RDES given by Equation 3.3, in which thereare two quantities u and v. Without loss of generality, for this discussion considerDu Dv. Here it is assumed that u is slow diffusing and v is fast diffusing. LPAtakes advantage of this fact to approximate the behaviour of pulse dynamics inthe PDE system with a system of ODES. With a large enough disparity betweendiffusion rates, consider the limit Du→ 0 (u is infinitely slow diffusing) and Dv→∞ (v is infinitely fast diffusing).Consider applying a perturbation to the system in the form of an arbitrarily tall,but very narrow, pulse. By assumption, the area of the pulse is essentially zero. Inthe limiting case where u does not diffuse, this perturbation will not spread. Speciesu will exhibit a local behaviour close to the perturbation and a global behaviour ev-erywhere else. These differences in behaviour are denoted ul and ug, respectively.Since v diffuses infinitely fast in the limiting case, any inhomogeneities causedby the perturbation will immediately be smoothed out. As v is always spatiallyhomogenous, only the global quantity, denoted vg, needs to be considered. Theperturbation results in distinctive differences between ul and ug, but vg remains un-changed. See Figure 4.1(a) for a visualization of this description. My discussionfocuses on a system of two variables, but LPA can easily be extended to higherdimensional systems (I will provide an example of a 3D system in Section 5.2).35DomainConcentrationPerturbationulugvga) b) ulugvgFigure 4.1: (a) Visualization of the LPA system. Solid lines indicate the ap-proximation by the LPA system, where Du→ 0 and Dv→∞. The dashedlines are a more realistic portrayal of the system, in which the effects ofdiffusion are evident. (b) Interaction between the LPA variables. Theglobal variables influence each other, while ul is only influenced by vg.Due to the lack of diffusion, the two forms of u do not interact. Since ulhas negligible area, it does not influence vg.Figure 4.1(b) is schematic which shows how vg, ug, and ul interact. In theasymptotic limit of Du = 0, the variable ul will only interact with itself and vg,because the lack of diffusion does not allow for communication to ug. By the samereasoning, ug depends only on itself and vg. The fast diffusing variable vg interactswith itself and ug, but not ul since the perturbation has a negligible area. In general,the global variables influence each other as they would in the well-mixed system.The local variables are influenced by the variables that are purely global and otherlocal variables. Pulse dynamics in the original system of PDES can be approximatedby the following system of ODES, denoted the LP system, which models the local36and global variables separately:duldt(x, t) = f (ul,vg), (4.3a)dugdt(x, t) = f (ug,vg), (4.3b)dvgdt(x, t) = f (ug,vg). (4.3c)Regions where the pulse grows, indicated by ul , can now be investigated throughbifurcation software packages and other ODE techniques. Approximations by theLPA system hold only at early stages of the pulse growth. In the original PDES, dif-fusion leads to interaction between the local pulse and the global variables, causingthe pulse to spread and approximations from the LPA to fail. Numerical simulationsare required to determine long-term pattern formation, but the LPA helps map pa-rameter space, providing a more directed search for patterning.While linear stability methods like Turing analysis consider very small pertur-bations, LPA considers large but narrow perturbations. This fact allows the LPAmethod to detect threshold behaviour that may be present in some systems. In ad-dition, the LPA method is able to detect regions in parameter space correspondingto Turing instabilities. Since LPA takes limiting assumptions about the spatial termsin the PDE, we do not retain the spatial information. Hence a full Turing analysisis required to determine which wavenumbers destabilize the system. Also, whilethe LPA technique will provide a good indication of pattern forming regions, it can-not guarantee the long-term behaviour of the system. It is useful to consider fullnumerical solutions to the original system of PDES in order to test the long-termdynamics of each parameter region. Given the relative ease of analyzing ODEScompared to PDES, the LPA method provides researchers with a useful and easy touse tool for investigating parameter regions.4.2.1 Mass conservative systemsA simplification to the LP system can be made when the quantities u and v areconserved. For a biological example of conservation, suppose u and v represent the37active and inactive forms of a particular protein, respectively. On a short timescale,production and decay of this protein are negligible. Then u and v can convertbetween the two forms, but the total mass remains constant.Consider the conserved amount of mass T and the length of the domain L.Then in the LP system:ug + vg =1L∫(ug + vg)dx≈1L∫(ug +ul + vg)dx =TL, (4.4)since the width of the pulse is negligible (∫uldx ≈ 0). Any instances of vg can bereplaced by T/L−ug. The LP system from Equation 4.3 becomes:duldt(x, t) = f (ul,T/L−ug) (4.5a)dugdt(x, t) = f (ug,T/L−ug). (4.5b)This simplification reduces the dimensionality of the LP system by removing theequation for vg. It is important to note that for many bifurcation programs, thisreduction is a necessary step as it avoids a singular Jacobian.4.2.2 Model signaturesEdelstein-Keshet et al. [9] refer to the LPA bifurcation diagrams as “signatures” forthe models. Each model has a unique LPA signature, which may vary dependingon which parameter is chosen. The signature indicates the pattern regimes thatare present, and how robust those pattern regimes are. We can compare modelsbased on how similar the LPA signatures are; if two models have comparable LPAsignatures, then they operate by similar mechanisms. Equivalently, if two modelsoperate through distinct mechanisms, then their LPA signatures will be differenteven if the models appear superficially similar.384.3 Benefits and caveatsThe local perturbation analysis is a useful tool for analyzing RDES which exhibit adisparity in diffusion rates. The following is a summary of the benefits this methodprovides:• LPA recovers information about the HSS of a model (number of steady states,their stability, and their values).• LPA detects threshold regions and predicts the height of a perturbation thatwill destabilize the HSS in these regions.• The LPA method detects regions of Turing instabilities. A Turing analysis istherefore not required to merely detect the presence of Turing instabilities.Information from LPA can help guide a Turing analysis if one is required.• LPA detects limit cycles in a model.• The LPA method is relatively simple to implement, and bifurcation packagesfor ODES are widely available and easy to use.• Information on the size of a patterning regime (i.e., robustness) with respectto any parameter can easily be determined from an LPA bifurcation diagram.The following list is a summary of the caveats:• The LPA method only works for this specific class of PDES, that is, RDESwith no advection and with a disparity in diffusion rates.• LPA only provides information on the initiation of a pattern: whether a pulse-like perturbation grows or decays. PDE simulations with the original set ofequations are required to determine the long-term dynamics.• Only a partial Turing analysis is recovered, as no information is gained onspecific wavenumbers or other spatial aspects of the pattern that arise.• LPA assumes extreme cases for diffusion, so any information specific to thediffusion coefficients is not provided. Since Turing regimes are dependenton the diffusions speeds, the boundaries may change during simulations.39• As the LPA is an approximation, the boundaries for the patterning regimesor locations for threshold behaviour are only approximate, and will not pre-cisely coincide with the PDE simulations. The extent to which the boundariesand threshold values differ depend on the PDES.40Chapter 5Demonstrating local perturbationanalysis through examplesThis chapter serves to demonstrate the use of the local perturbation analysis (LPA)method through two models. I will use the Wave Pinning model by Mori et al.[49] to outline and demonstrate the basic use of LPA. This model was introducedin Section 3.3. The Actin Waves model by Holmes et al. [20] is an extension of theWave Pinning model. I will use this model to highlight additional features of LPA,such as the connection between the patterning regimes and the LPA eigenvalues.The Wave Pinning model is a 2D system, and so applying LPA using the guidelinesfrom Chapter 4 will be straightforward. Applying LPA to the Actin Waves modelis slightly more complicated, as it is a 3D system.The LPA bifurcation diagrams are produced using MatCont, a numerical con-tinuation program designed for MATLAB (MathWorks) [8]. I used the Crank-Nicolson method with no-flux boundary conditions implemented in MATLAB (seeAppendix B) to solve the reaction-diffusion equations (RDES). The models werestimulated by introducing an arbitrarily large pulse with a narrow width to the ini-tial conditions. 5% of the 1D domain was stimulated, and the size of the pulsedepended on the patterning regime: threshold regions require a specifically-sizedpulse, while a very small perturbation to the homogeneous steady state (HSS) issufficient in a Turing regime.415.1 Wave Pinning model5.1.1 Model equationsThe Wave Pinning model is a reaction-diffusion system representing the active andinactive forms of a single GTPase. The equations involve an inactive form v whichdiffuses quickly, and a slowly diffusing active form u. The RDE system, presentedin [49], is as follows:∂u∂ t = Du∂ 2u∂x2 + v(k0 +γu2K2 +u2)−δu, (5.1a)∂v∂ t = Dv∂ 2v∂x2 − v(k0 +γu2K2 +u2)+δu. (5.1b)The parameter values are given in Table 5.1. Notice that the reaction terms in theWave Pinning model have the form g(u,v) = − f (u,v). There is only interconver-sion between the forms u and v, no synthesis or decay of the proteins; the WavePinning model is conservative.Description of model terms and parametersThe term δu represents the linear conversion from the active form u into the in-active form v, where δ is the inactivation rate. Activation (i.e., conversion fromof v into u) is given by: v(k0 +γu2K2+u2). We will consider this term in two parts.First is the basal activation of u, given by the rate parameter k0. Second is thepositive feedback to activation, modelled by a second degree Hill function. Thepositive feedback term is very low at resting levels of u, when u K, but the basalactivation rate retains a constant level of activation. In the presence of a stimulus,positive feedback increases to become the dominant term for activation. Param-eters γ and K are the maximal rate of autocatalysis (positive feedback) and thesaturation parameter, respectively.42Parameter Value DescriptionDu 0.1 Diffusion of active form (µm2s−1)Dv 10 Diffusion of inactive form (µm2s−1)k0 0.067 Basal activation of u (s−1)γ 1 Maximal rate of conversion (s−1)K 1 Saturation parameter (µM)δ 1 Inactivation rate (s−1)T 2.27 Conserved concentration (µM)L 10 Cell length (µm)Table 5.1: Parameter values for the Wave Pinning model in Equation 5.1. Val-ues were obtained from [49].The wave pinning mechanismIn certain parameter regimes of the Wave Pinning model, a stimulus applied to theactive form u will lead to the formation of a polarized domain. Autocatalysis, orpositive feedback of u to itself, along with bistability, contribute to the formationof a wave front which travels into the domain. Species v is converted into u as thewave front moves forward, and the amount of v is depleted due to conservation.When the inactive form v is sufficiently depleted, the wave front decelerates andeventually stalls. The stalled wave results in a macroscopic difference between thelevel of activation at the front and back of the cell, indicating polarization. Notethat the level of v is essentially homogeneous, and it is the active form u which hasa higher concentration at the front of the cell and a lower concentration at the back.5.1.2 Preliminary investigationsI conducted a preliminary analysis of the Wave Pinning model using the originalparameters (Table 5.1) and stimulus from [49]. The stimulus is given by:43kstim ={s(t)(1+ cospix) if x≤ 1,0 if x≥ 1,(5.2)where:s(t) =S 0≤ t ≤ t1,S4(1+ cos(pi t−t1t2−t1))t1 ≤ t ≤ t2,0 t ≥ t2.Here, S = 0.5, t1 = 20s, and t2 = 25s. The stimulus is incorporated into the WavePinning model by appending kstim to the model equations: Equation 5.2 is added toEquation 5.1(a), and subtracted from Equation 5.1(b) to conserve mass. After thetime indicated by t2, the stimulus has decayed and no longer affects the model.Figure 5.1(a,b) are 1D simulations of Equation 5.1, presented as kymographsand line plots. A kymograph shows spatiotemporal dynamics of a given variableas intensity plots. The x-axis represents time and the y-axis is the spatial domain,with the front of the cell oriented at y = 0. These diagrams show the wave pin-ning response to the above stimulus, where the wave front stalls and stabilizes ina polar pattern. Bifurcation plots for the well-mixed system (i.e., reaction kinet-ics) in Figure 5.1(c,d) are conducted with respect to the parameters δ and k0. Inboth diagrams there is a unique HSS and it is universally stable, indicating that anypotential pattern formation does not extend from the reaction kinetics alone.5.1.3 The LPA systemI will outline a procedure for using LPA to find patterning regimes in the WavePinning model. As noted previously, Du Dv and so u has both local and globalforms, and v is purely global. We can model the initial dynamics of a pulse inEquation 5.1 by using the following system of three ordinary differential equations44TimePosition of u  100 200 300051012TimePosition of v  100 200 30005101.851.950 5 10012DomainConcentration of ut = 0t = 5, 50, 200Stimulus00.51Stimulus (added)0 1 2012δConcentration of u0 0.1 0.2012k0Concentration of ua) b) c) d) tFigure 5.1: Preliminary analysis of the Wave Pinning model in Equation 5.1using the default parameter values given in Table 5.1. (a) Kymographof solutions u (top) and v (bottom). Note that u sustains an inhomoge-neous (polarized) distribution, whereas v is nearly uniform across thedomain. (b) Line plot of solutions for u at t = 0, 5, 50 and 200 seconds(time progresses in the direction of the arrow). The stimulus (right axis)is shown in red. Bifurcation diagrams of the well-mixed system withrespect to parameters (c) δ and (d) k0. These diagrams show that thereis a unique and stable HSS for this set of parameters.(ODES):duldt= vg(k0 +γu2lK2 +u2l)−δul, (5.3a)dugdt= vg(k0 +γu2gK2 +u2g)−δug, (5.3b)dvgdt=−vg(k0 +γu2gK2 +u2g)+δug. (5.3c)45The spatial element of the equations is removed, and we no longer consider diffu-sion. The variables involved in the LP ODES are: ul (local u, represents the pulsedynamics), ug (global u, represents the global dynamics but not the pulse), and vg(global v). Notice that the interactions between the LP variables follow the samescheme as presented in Equation 4.3.Recall that the Wave Pinning model is conservative, since in the original systemEquation 5.1, the reaction kinetics have the form f (u,v) = −g(u,v). Define Tas the total conserved concentration of u and v. The value of T is provided inTable 5.1. We can make the substitution vg = T − ug, reducing our system inEquation 5.3 from three ODES to two:duldt= (T −ug)(k0 +γu2lK2 +u2l)−δul, (5.4a)dugdt= (T −ug)(k0 +γu2gK2 +u2g)−δug. (5.4b)5.1.4 LPA bifurcation analysisFrom Equation 5.4, bifurcation diagrams with respect to parameters δ and k0 (thesame parameters from Figure 5.1) were constructed, and are shown in Figure 5.2.Dash-dot lines represent regions of instability, and solid lines are stable regions.The y-axis represents ul to ensure pulse dynamics are shown in the LPA bifurcationdiagrams. Notice that, with the minor exception of a new unstable region, theblack branches are the same as the well-mixed bifurcations from Figure 5.1(c) and(d). I will refer to these well-mixed branches as “global” branches. The pulsedynamics are represented by the local branches shown in red. Clearly the additionof diffusion results in interesting new behaviour.Several distinct regions in Figure 5.2 relate to three general patterning regimes:stable, Turing, and threshold. The regions are defined by examining vertical cross-sections of the LPA bifurcation diagram. A stable regime will contain a stableglobal branch with no local branches. The Turing regime appears as an unstable460 0.5 1 1.5 200.511.52δConcentration of ulBP BP 0 0.05 0.1 0.15 0.200.511.52k0Concentration of ulBP BP a) b) I II III IV V Figure 5.2: LPA bifurcation diagrams for the Wave Pinning model with re-spect to parameters (a) δ and (b) k0. Legend: Global branch (black);Local branch (red); Solid line (stable); Dash-dot line (unstable); BP =branch point. There are three types of regions in each diagram: purelystable, threshold response, and Turing regime. The boundaries betweenpatterning regimes were explicitly drawn in (a) and labelled as: I sta-ble, II downward threshold, III Turing, IV upward threshold, and V sta-ble. In (b), the regimes from left to right are: upward threshold, Turing,downward threshold, and stable, but are not distinguished by lines. Notethat an upward (downward) threshold regime requires a perturbation ofa certain size that is above (below) the HSS. More detail is presented inthe branch. Note that the presence of local branches in the Turing regime isinconsequential since, once the pulse starts to diffuse, any predictions by LPA breakdown. In a threshold regime, vertical cross-sections contain three branches: awell-mixed branch followed by two local branches. Here, the global branch isstable, the interior local branch is unstable, and the exterior local branch is stable.The unstable local branch represents the threshold size of a perturbation that willdestabilize the HSS. If the local branches are located above the global branch, aperturbation in the upward direction (i.e., a value larger than the HSS) is requiredto initiate polarization. The local branches can also be located below the globalbranch, and a downward perturbation is then required. I will refer to these differentthreshold regimes as “upward threshold” and “downward threshold” respectively.47The patterning regimes in Figure 5.2(a) are labelled as follows: I stable, IIdownward threshold, III Turing, IV upward threshold, and V stable. Figure 5.2(b)does not have explicitly drawn boundaries between the regimes, but we can see thatthey occur in the following order: upward threshold, Turing, downward threshold,and stable. In addition to the types of patterning regimes, LPA also shows theparameter ranges over which they occur. As both bifurcation diagrams exhibit thesame types of responses, but in a different order, I will only discuss Figure 5.2(a)in more detail.The first regime is purely stable, and occurs when δ < 0.36. In this regime,there are no perturbations that will lead to pattern formation of the model. InRegion II, when 0.36 < δ < 0.81, patterns are formed when a threshold valueis surpassed. The threshold is given by the unstable local (red dash-dot) branch,and the lower steady state occurs at the stable local (red solid) branch. Here adownward perturbation is required, which corresponds to a depletion of activatoru. In Region III, 0.81 < δ < 0.98, the black dash-dot branch indicates a Turingregime, in which patterns will form from small amplitude oscillations. Anotherthreshold region occurs in Region IV over the range 0.98 < δ < 1.21. In thisregion an upward perturbation, or an increased amount of activator u, is required.In Region V, when δ > 1.21, the system is again stable.5.1.5 PDE simulationsNow that we have started to map the parameter space using the LPA bifurcations,simulations of the original RDES in Equation 5.1 can be conducted for each regime.Note that our map is for a limited set of parameter values, since the LPA bifurcationdiagrams are continuations of just one parameter at a time. In Figure 5.3, the fivedifferent kymographs correspond to the five different regions outlined above andshown in Figure 5.2(a). For these simulations, I used a stimulus that is applied as aperturbation to the initial conditions of u and does not depend on t. I applied a pulseto 5% of the far left edge (i.e., y = 0 in the kymograph) whose size is determinedby the type of patterning regime, as described in the paragraph below.In both the stable regimes, no patterning results regardless of the stimulus size(Figure 5.3(a,e)). A very small perturbation (uss + 0.05, where uss is the steady48100 200 3000510Time100 200 3000510Time100 200 3000510100 200 3000510x100 200 3000510Timex012a) b) c) d) e) Figure 5.3: Kymographs of the Wave Pinning model in Equation 5.1, con-firming the patterning regimes in Figure 5.2(a). The parameter δ takesthe values: (a) 0.2 (stable), (b) 0.7 (threshold, down perturbation) (c)0.9 (Turing) (d) 1.1 (threshold, upward perturbation) and (e) 1.3 (sta-ble). Other parameter values are given in Table 5.1.state) results in a polarization response in the Turing regime (Figure 5.3(c)). In thetwo threshold regimes, either a large downward or upward perturbation is required.In the downward threshold regime (Figure 5.3(b)), the perturbation has a value of0.1, and must be applied to 20% of the domain for the desired pattern to form. Here,I increased the width because an increase in magnitude would lead to a negatively-valued perturbation. A value of 2 is sufficient for patterning in Figure 5.3(d), whereδ is in the upward threshold regime.5.2 Actin Waves model5.2.1 Model introductionIn the second example, I apply LPA to a model developed by Holmes et al. for actinwaves. A nondimensionalized version of the model is presented in [43], alongwith connections to LPA and investigations of pattern transitions. The Actin Wavesmodel was inspired by biological experiments showing evidence of in vivo wave-like patterns of actin [47, 74, 75, 78]. Actin monomers diffuse within the cytosolof most cells, and nucleate to form linear polymers (F-actin), as discussed in Sec-tion 2.1. This partial differential equation (PDE) model is based on the biologicalinteractions between F-actin and signalling molecules that promote the nucleation49of F-actin, here termed nucleation promoting factors (NPFS). The GTPases arecandidates for the NPF in the Actin Waves model.Model equationsThe nondimensionalized model is given as follows:∂A∂ t = DA∂ 2A∂x2 + f (A, I,F), (5.5a)∂ I∂ t = DI∂ 2I∂x2 − f (A, I,F), (5.5b)∂F∂ t = εh(A,F), (5.5c)where:f (A, I,F) =(k0 + γA3A30 +A3)I−(s1 + s2F1+F)A, (5.5d)h(A,F) = knA− ksF. (5.5e)A table of parameter values is given in Table 5.2. In this model, A and I rep-resent the active and inactive forms, respectively, of the NPF. Since it is assumedthat A is membrane-bound and I is located in the cytosol, DA DI . F-actin is rep-resented by the variable F . Notice that F-actin does not diffuse. Even though thecytoskeleton is dynamic, once the actin filaments form, they are bound to adhesioncomplexes and do not diffuse within the cell. A schematic depiction of the ActinWaves model reaction terms is presented in Figure 5.4.Description of model terms and parametersLet us focus on Equation 5.5a and Equation 5.5b first. These two equations aremodified from the Wave Pinning model for GTPases that I discussed in Section 5.1.The function f (A, I,F) represents interconversion terms, allowing the NPF to changebetween its two states. Production and decay of the NPF are assumed to be neg-ligible. The first term of f (A, I,F) represents activation, or conversion of I intoA, and depends on the amount of I present. Activation occurs at a basal rate (k0),50Parameter Value DescriptionDA 0.00033 Diffusion of active formDI 0.033 Diffusion of inactive formk0 0.1-0.35 (0.05) Basal activation rate of NPFγ 1 Magnitude of NPF positive feedbackA0 0.4 Level of A for positive feedback to occurs1 0.7 (0.5) Basal inactivation rate of NPFs2 0.7 (0.6) F-actin mediated inactivation rate of NPFkn 2 F-actin nucleation rate by NPFks 0.25 F-actin disassembly rateε 0.1 F-actin dynamics timescaleTnp f 1 Total amount of NPFL 1 Length of the cellTable 5.2: Parameter values for the nondimensionalized Actin Wave model,taken from [43]. Unless otherwise stated, the bracketed values are onlyused in the preliminary investigations.but is also influenced by autocatalysis (the nonlinear term). Inactivation, conver-sion from A to I, is represented by the second term and depends on the value ofA. Again, there is a basal rate of inactivation s1, and here a nonlinear term rep-resents increased inactivation when F-actin is present. Since F is present in theinactivation term, negative feedback exists from F to A.Now consider the equation for actin filaments, Equation 5.5c. The activeform A is found on the membrane of the cell and leads to nucleation of the actinmonomers to create F-actin, indicated by the term knA. When the NPF is in theinactive state I, it does not affect F . Depolymerization of F-actin is representedby ksF . The cell receives a stimulus which triggers a signal transduction pathway,51A FIFigure 5.4: Schematic Actin Wave model reaction diagram. The slowly dif-fusing variables A and F are shaded, I is quickly diffusing and solid.There is interconversion between A and I, positive feedback from A toitself (autocatalysis), and negative feedback from F to A as F is pro-duced. The dashed arrow represents depolymerization of F .leading to a high concentration of A in the part of the cell closest to that stimulus.In response to this activation, actin filaments nucleate from actin monomers. Areorganization of the cytoskeleton follows, and allows the cell to respond to its en-vironment, for example by migrating along a chemical gradient. Actin monomersare available in excess, and we assume that they remain at a constant level that doesnot limit the speed of reactions. The F-actin dynamics are assumed to occur on atimescale slower than the NPF, approximated through the presence of ε 1, whichis small in comparison to the reaction kinetics timescale of the NPF. Parameter val-ues are justified in reference [20].Pulse mechanismUnder certain parameter regimes, when a stimulus is applied to the Actin Wavesmodel, the wave pinning mechanism is initiated. Bistability and autocatalysis ofthe active form A contribute to the production of a wave front, which travels into thedomain. In the Wave Pinning model, mass conservation between (analogous) A andthe inactive form I resulted in a stalled wave. However, negative feedback from theactin filaments F to A is present in the Actin Waves model. As the concentrationof F increases, negative feedback pulls part of the wave of A back down to thelower steady state, creating a pulse of A rather than a wave. That is, some of A isconverted back to I. Since mass conservation no longer limits activation of A atthe wavefront, the pulse continues to move forward, while the trailing end follows.52Therefore, the pulse mechanism can lead to self-organized dynamic patterns.5.2.2 Preliminary investigationsFigure 5.5 presents preliminary investigation of the Actin Waves model. In part (a),kymographs show PDE solutions for each variable A, I, and F . The dynamic patternwas initiated by introducing a large pulse (A0 = 0.5) to 5% of the 1D domain, atthe far left (i.e., y = 0) edge. Note that the stimulus is only applied to the initialconditions of the model (refer to Section 5.1.2), and does not depend on time.Parameter values for this simulation give rise to reflecting pulses, but other typesof patterns including boundary localized patterns and pulse trains, are possible[20, 43]. We can also examine the reaction kinetics, given by the functions inEquation 5.5e. The bifurcation diagram in part (b) shows no evidence of patternformation with respect to the basal activation parameter k0; the entire branch isstable. This observation indicates that in the absence of diffusion there are nointeresting dynamics.5.2.3 The LPA systemThe Actin Waves model Equation 5.5 is a reaction-diffusion system in which thechemical species diffuse and react on different timescales (DA DI,DF = 0). Wecan make the approximation that I is fast diffusing, and our slow diffusing vari-ables are A and F . Using the method outlined in Section 4.2, pulse dynamics inEquation 5.5 can be approximated by a system of five ODES. The spatial element isremoved, so we no longer consider diffusion, and the variables involved become:Ag, Al , Ig, Fg, and Fl . The inactive form I is a fast diffusing variable, and does notrequire a local term. These modifications lead to the following LP system:53A  200 400 600 800 100000.510.20.6I  200 400 600 800 100000.510.70.76F  200 400 600 800 100000.5124PositionTime0 0.1 0.2 0.300. of Aa) b) Figure 5.5: This figure shows preliminary investigation of the Actin Wavesmodel in Equation 5.5, using standard parameter values from [43] (thebracketed values in Table 5.2). (a) Kymographs for the PDE simulationsof A, I, and F . The x-axis represents time and the y-axis is the spatialdomain, with the front of the cell oriented at y = 0. These parametervalues result in a pattern of reflecting waves, which travel from one endof the cell to the other. (b) Bifurcation diagram of the reaction kineticsfrom Equation 5.5e with respect to the basal activation parameter k0.The active NPF A is represented on the y-axis. There is a unique HSS forthese parameter values, and it is universally stable.dAgdt= f (Ag, Ig,Fg), (5.6a)dAldt= f (Al, Ig,Fl), (5.6b)dIgdt=− f (Ag, Ig,Fg), (5.6c)dFgdt= εh(Ag,Fg), (5.6d)dFldt= εh(Al,Fl). (5.6e)Refer to Figure 5.6 for a schematic depiction of the interactions between variables54AlAgIgFlFgFigure 5.6: Interactions between the Actin Wave LPA variables from Equa-tion 5.6. Global variables are solid, local variables are shaded. Thereare interactions between every global variable, and between every localvariable. The local variables Al and Fl are influenced by Ig. Since A andF are slow variables in the original system Equation 5.5, they lack dif-fusion in the LP system. Therefore the local variables interact with eachother at the perturbation, but do no interact with the “slowly diffusing”global variables Ag or Fg. The perturbation has a small area, so the localvariables do not influence the purely global Ig.of the LP system.Substituting these variables into the reaction kinetics gives the following sys-tem of ODES:dAgdt=(k0 + γA3gA30 +A3g)Ig−(s1 + s2Fg1+Fg)Ag, (5.7a)dAldt=(k0 + γA3lA30 +A3l)Ig−(s1 + s2Fl1+Fl)Al, (5.7b)dIgdt=−(k0 + γA3gA30 +A3g)Ig +(s1 + s2Fg1+Fg)Ag, (5.7c)dFgdt= ε(knAg− ksFg), (5.7d)dFldt= ε(knAl− ksFl), (5.7e)55In the original model, conservation of mass exists between A and I. Therefore,as outlined in Section 4.2, we can reduce the LP system to four equations by intro-ducing Ig(t) = T −Ag(t), where T in this case represents the total amount of NPF.This substitution gives the following reduced set of LPA equations:dAgdt= f (Ag,T −Ag(t),Fg), (5.8a)dAldt= f (Al,T −Ag(t),Fl), (5.8b)dFgdt= εh(Ag,Fg), (5.8c)dFldt= εh(Al,Fl). (5.8d)5.2.4 LPA bifurcation analysisBifurcation diagrams with respect to different parameters can now be constructedusing the LP system given by Equation 5.8 and parameters from Table 5.2. Fig-ure 5.7(a) is a bifurcation diagram for the basal activation rate of the NPF, k0, andFigure 5.7(b) depicts eigenvalues of the global branch. In part (b), notice thatthere are both global and local eigenvalues, which arise from the global or localvariables. The global variables and the variables from the well-mixed system areequivalent; the spatial component is realized through the local variables. FromFigure 5.7(b), we can see that the instability on the global branch is a result of thelocal eigenvalues, i.e., the spatial features.The black curve in Figure 5.7(a) represents the global part of the LP system,and is shared with the well-mixed system shown in Figure 5.5. It occurs when thereis no perturbation, and Al = Ag, Fl = Fg. The y-axis is the concentration of Al , thelocal part of the active NPF. By using Al , we can examine the local state (red), andrecover the well-mixed state. Solid versus dash-dot lines represent stable or unsta-ble solutions, respectively. We can now make inferences about whether patternswill form and their potential initiation given different values of k0. Eigenvaluesfor the well-mixed branch are provided in Figure 5.7(b); they were computed by560 0.1 0.2 0.300.20.40.6k0Concentration of AlBP BP H H 0.15 0.2 0.25 0.3−λmax Re( λmaxLP )Im(λmaxLP )λmaxglobalBP BP H H I II III IV V I II III IV Va) b) Figure 5.7: (a) LPA bifurcation diagram and (b) maximal eigenvalues forEquation 5.8 versus the basal activation rate k0. Bifurcations (• BP= branch point,  H = Hopf bifurcation) and pattern regimes can becompared on both diagrams. Parameter values are given in Table 5.2.(a) Black = global branch, red = local branch, solid = stable, dot dash= unstable. There are four different regions to notice: threshold be-haviour (I), oscillatory Turing instability (II and IV), stationary Turinginstability (III), and a stable regime (V). (b) The value of the maximumeigenvalue λmax on the global branch among the global variables (Agand Fg, black) or the local variables (Al and Fl , red). The imaginary partof λmax is included for the local variables.MATCONT during the LPA bifurcation analysis. Here I show the maximal eigen-value among the global variables (Ag and Fg, black), and among the local variables(Al and Fl , red). Corresponding pattern regimes are labelled on both diagrams forcomparison.Region I occurs between k0 = 0 and the first H (Hopf) bifurcation. The globalbranch is stable, the lower local branch is unstable, and the upper local branch isstable. The eigenvalues also show stability in the global branch, as Re(λmaxG ) andRe(λmaxLP ) are negative. This region indicates a threshold response to a perturbation:a narrow perturbation with a height exceeding the threshold level (given by theunstable local branch) will be pushed to the upper stable steady state defined bythe stable local branch. Note that this description only defines the initiation of apattern with threshold dynamics. Once the perturbation spreads, the LPA breaks57down and no longer approximates the pulse dynamics in the PDE system. Theultimate pattern that forms must be determined through long-term simulations ofthe full PDE model. This fact is also true for the subsequently described regions.Regions II and IV occur between H and BP (branch point) bifurcations. Inthese regions, the global branch is unstable. Examining the eigenvalues for thisregion reveals that the Re(λmaxG ) < 0, indicating that the well-mixed system is sta-ble. In the case of the local eigenvalues however, Re(λmaxLP )> 0 and Im(λmaxLP )> 0.This region represents an oscillatory Turing instability, and it is through the localvariables that this instability occurs.Region III, located between the two BP bifurcations, is unstable. Examin-ing the eigenvalues of the LP system, we notice again that Re(λmaxG ) < 0 andRe(λmaxLP ) > 0. Now however, we have that Im(λmaxLP ) = 0. This region then rep-resents a stationary Turing instability. In this and the previously defined region,pattern formation results from a small perturbation which spreads. The stability ofthe local branch is inconsequential, as the LPA breaks down once the perturbationis no longer local.Region V occurs for values of k0 beyond the second H bifurcation. In thisregion, the global branch is stable and there are no local branches present, indicat-ing an absence of any pattern forming dynamics. Again the eigenvalues indicatestability in the global branch, as Re(λmaxG ) < 0 and Re(λmaxLP ) < 0.In summary, the Al-k0 bifurcation diagram reveals four different regions:• Region I: threshold dynamics,• Regions II and IV: oscillatory Turing instability,• Region III: stationary Turing instability,• Region V: no pattern formation.We also see the robustness of each region. Patterns will only form provided k0 issmall enough (i.e., smaller than the second H bifurcation). This general patternforming region is broken down further, given the type of stimulus that initiates thepatterning. Very small regions are not robust to changes in the bifurcation parame-ter, and may not be biologically relevant as a result. Conversely, large regions willbe robust against changes in the bifurcation parameter.580 0.5 1Concentration of Al0 0.5 2Concentration of AlBP BP BP BP H H H H a) b) Figure 5.8: Bifurcations for the local active NPF variable Al against (a) thebasal NPF inactivation parameter s1 and (b) the F-actin mediated NPFinactivation parameter s2. Legend: black = global branch, red = localbranch, solid = stable, dot dash = unstable, • BP = branch point,  H =Hopf bifurcation. Parameter values are given in Table 5.2.Additional bifurcation diagrams are shown in Figure 5.8; each diagram exhibitsregions that are qualitatively the same as one of the four described in the abovelist. The two figures show bifurcations for parameters (b) s1 and (c) s2, with Alon the y-axis. Notice that the bifurcation parameters s1 and s2 both correspondto inactivation of the NPF, compared to k0 which is an activation parameter. Thisdifference results in diagrams that horizontally mirror Figure PDE simulations and connection to LPAIn order to observe the long-term patterns that form, PDE simulations using theoriginal model in Equation 5.5 were conducted. All diagrams were produced usingthe parameters given in Table 5.2 unless otherwise specified. In simulations forthreshold regimes, the perturbation should be large enough to push the level of Apast the threshold described by the local branch in Figure 5.7. It is important toremember that the LPA is an approximation, and so inferences about the size ofthe perturbation or pattern initiation on the boundaries may not be exact. I adjustthe height of the pulse so as to trigger the pattern in the full PDE simulations. Thesolutions shown in the kymograph are for A, the active NPF.59The simulations presented in Figure 5.9 move through the k0 parameter do-main, from small to large values, to illustrate each of the regions described above.Parameter values for (a) and (b), with k0 = 0.05,0.12, are in the threshold regime,and only a large perturbation results in pattern formation. As k0 increases, the re-quired size for perturbation decreases from 2.5 (at k0 = 0.05) to 0.4 (at k0 = 0.12).Parts (c)-(h) feature Turing instabilities, and only require a small perturbation(Ass + 0.05) to initiate the wave. As predicted, k0 = 0.35, shown in Figure 5.9(i),is stable to any perturbations.Notice the transition of patterns as k0 increases. At first, a pulse of active NPFcrosses through the domain and reflects off the boundaries. As k0 increases, thewidth and frequency of the pulses increase. In part (f), the initial pattern changesbut resolves back into reflecting pulses at longer times. Part (g) features a com-pletely new pattern: standing waves. The amplitude of the pattern decreases dras-tically in part (h), and diminishes altogether in the stable regime, part (i). It isinteresting to see that, even though (d)k0 = 0.1814 and (g)k0 = 0.25 are locatedat the branch points, the patterns are very different. The same observation can bemade about the oscillatory Turing regimes II and IV in parts (c)k0 = 0.171 and(h)k0 = 0.26. Also, while the threshold and Turing regimes initiate patterns in verydifferent ways, the results from (a)-(f) are all reflecting pulses. These observationsemphasize that LPA only predicts how the pattern is initiated, not the long-termresult.5.3 Wave pinning versus actin wavesThe first model I analyzed in this section was the Wave Pinning model, the 2Dsystem in Equation 5.1 representing the two forms of a single GTPase. IntroducingF-actin and actin-dependent negative feedback to the Wave Pinning model lead tothe Actin Waves model, the 3D system in Equation 5.5. The Wave Pinning modelis essentially nested in the Actin Waves model. If we examine the kymographsfor the Wave Pinning and Actin Waves models in Figure 5.1 and Figure 5.9, re-spectively, the simulations appear very different. Introducing F-actin results in thepossibility for different kinds of dynamic patterns. Visually, we can see how ex-tending the Wave Pinning model into the Actin Waves model affects the ultimate60400 80001x400 80001400 80001Timex400 80001Time400 80001Time400 80001400 80001400 80001400 80001xa) b) c) d) e) f) g) i) h) 0 0.5Figure 5.9: PDE simulations for the Actin Wave model in Equation 5.5 withparameters from Table 5.2. Values of k0 were chosen to investigate re-gions in Figure 5.7. The kymographs show the results for A. Simula-tions for Turing regimes all used the same small stimulus, ASS +0.005,perturbations in the threshold regimes are the minimum height whichgives rise to patterns. (a) Threshold regime k0 = 0.05, perturbation= 2.5; (b) threshold regime k0 = 0.12, perturbation = 0.4; (c) oscilla-tory Turing regime k0 = 0.171; (d) at the first branch point k0 = 0.1814;(e) Turing regime k0 = 0.21; (f) Turing regime k0 = 0.24; (g) at the sec-ond branch point k0 = 0.25; (h) oscillatory Turing regime k0 = 0.26; (i)stable regime k0 = 0.35, perturbation = 3.61pattern. However, from the kymographs alone we do not know if the core mecha-nism responsible for initiating these patterns has changed.Edelstein-Keshet et al. [9] compare these two models using the LPA signaturesfor the common parameter k0 (in this thesis: Figure 5.2(b) for Wave Pinning andFigure 5.7(a) for Actin Waves). Extending from Wave Pinning to Actin Waves doesnot change the core bifurcation diagram appreciably, but it does add the additionallayer of Hopf bifurcations. With a few quantitative differences, the well-mixedbranches and local branches have the same shape. Both LPA bifurcation diagramscontain similar Turing, threshold, and stable regimes. This observation indicatesthat the core mechanism behind pattern initiation does not change with the additionof F-actin dynamics, and in the Actin Waves model, the NPF equations contain thekey motif for patterning.Using the Wave Pinning and Actin Waves models as examples shows how wecan transition from simple to detailed models by adding or changing feedback andintroducing new components. Simulations allow us to observe the changes in dy-namics between simple and extended models, but they do not provide a completepicture. By comparing the bifurcation signatures, LPA can be used to observe thechanges in pattern formation dynamics that occur when simple models are ex-tended or modified into more detailed models.5.4 Additional application: periodic cell migrationIn this section, I will show qualitatively how LPA could be used to explain the spon-taneous emergence of periodic cell migration in a paper by Camley et al. [6]. Theauthors develop a 2D model in response to observations of cell motility on a mi-cropatterned substrate. In experiments, surfaces with micropatterns are able to in-duce polarization, directed motion, and periodic motion in crawling cells [11, 69].The computational model is composed of several modules. The stresses inducedthrough actin protrusion and myosin contraction are coupled to modules for cellshape, substrate adhesion, and flow of the cytoskeleton. The Wave Pinning RDESin Equation 5.1 model an actin promoter, and serve to induce polarization. The au-thors also model the substrate with a micropatterned strip in the center, and assumethat adhesions between the cells and substrate only occur on the micropattern.62Simulations of the computational model gave rise to dynamic results consis-tent with biological observations. In particular, they observed periodic migrationwhich emerged without underlying oscillatory dynamics in any of the modules.The simulated, polarized cell began migrating on the micropatterned strip. Forcesincurred on the cell’s shape during migration caused the cell to shrink, which led todepolarization and a corresponding halt in migration. After migration ceased, thecell suddenly increased in size and repolarized in the opposite direction. Migrationresumed in this new direction, leading to contraction of the cell, and the cycle re-peated. Interestingly, the authors also developed a 1D version which produced thesame results.Mori et al. [49] show analytically that the length of a cell must exceed a certainvalue in order for polarization to occur. When the cell contracts, it depolarizes.When the cell expands again, the actin promoter requires a perturbation for repo-larization. This perturbation, Camley et al. [6] suggest, occurs when an increasedamount of myosin at the trailing edge causes the leading edge to expand faster,resulting in a depletion of active actin promoter at the front. When this depletedvalue reaches a certain threshold, the cell repolarizes with a high concentration ofactive actin promoter in the opposite direction.Since polarization arises from the Wave Pinning model, I decided to examine abifurcation diagram for the Wave Pinning LP system in Equation 5.4. A qualitativeanalysis of Figure 5.10, whose bifurcation parameter T is the total concentrationof actin promoter, is consistent with the authors’ suggested causes of periodic mi-gration. Parameter values were taken from the supplementary information of [6].The authors indicate that the model begins in the wave pinning regime (II). An in-crease of the total concentration T results when the cell contracts, possibly endingin region (V) where the cell has a stable regime. As the cell expands, T decreasesto region (IV). Recall that at this stage the leading edge is expanding faster, corre-sponding to depleting active actin promoter at the front. In region (IV), a downwardperturbation (i.e., depletion) of sufficient size leads to polarization (see Figure 5.3).But since depletion drives the active actin promoter to the lower steady state at thefront, the higher steady state is now at the “back” of the cell. Migration proceedsin this new direction, and as the cycle repeats, T moves between region (IV) and(V).631.5 2 2.5 300.40.81.2TConcentration of AlBP BP I II III IV VFigure 5.10: An LPA bifurcation diagram for the Wave Pinning model inEquation 5.4. The bifurcation parameter is the total concentrationof the actin promoter T . Black bars indicate the different patterningregimes, and the red arrow indicates the possible regions in whichspontaneous periodic migration occurs for the computational modelin [6]. As the cell shrinks, T increases, and as the cell expands, T de-creases. For parameter values, see the supplementary information in[6].64Chapter 6Analyzing cell polarizationmodelsThe following sections contain a collection of models obtained from literature oncell polarization. These models were introduced in Chapter 3; my goal now is toanalyze them. Analysis for each model will be conducted similarly to the modelsin the previous chapter. Throughout this section, I will apply local perturbationanalysis (LPA) to each model, perform bifurcation analysis using the LP system tovisualize the pattern forming regions, and use numerical partial differential equa-tion (PDE) simulations to test these regions.As the analysis for the next sections is repetitive and we want to be able tocompare each system, what follows here is some housekeeping for consistency.Recall the general form for a system of two reaction-diffusion equations (RDES)presented in Chapter 3:∂u∂ t = Du∂ 2u∂x2 + f (u,v), (6.1a)∂v∂ t = Dv∂ 2v∂x2 +g(u,v). (6.1b)In general, the variable u represents an activator, or the active form of a signallingprotein while v corresponds to an inhibitor or the inactive form. The domain is65a slice across the cell diameter, given by 0 ≤ x ≤ L, and I use no flux boundaryconditions:∂u,v∂x∣∣∣∣x=0= 0 =∂u,v∂x∣∣∣∣x=L, (6.2)Unless otherwise noted, in all cases, Du  Dv. All polarization models will bepresented in this form or with intuitive modifications (i.e., in the case of a threedimensional system of equations). I use the Crank-Nicolson method implementedin MATLAB for 1D PDE simulations, and MatCont [8] for the bifurcation analysis.Analysis for all models begins with simulations using the type of stimulus fromthe original publication, and a default set of parameter values from either the orig-inal publication or the cell polarization review article by Jilkine and Edelstein-Keshet [29]. Bifurcation diagrams of the reaction kinetics are also included. Thesediagrams indicate any patterning regimes that occur purely within the reaction ki-netics. Line plots display the system as it approaches polarization and the originalstimulus that was applied.When confirming the regimes mapped out by the LPA bifurcation diagrams(using numerical simulations), a consistent form for the stimulus is applied amongall the models. In general, the stimulus is applied to the leftmost 5% of a domainof length L (in the kymographs, the stimulation occurs at the bottom) as follows:u0 ={stimulus if x < 0.1L,uSS if x > 0.1L,(6.3)where:stimulus ={uSS +0.005 for Turing regimes,ε for threshold regimes.The value of ε depends on the patterning regime, uSS is the homogeneous steadystate (HSS) of u for the given parameter values, and L is the length of the 1D cell.In a Turing regime, the model only needs a small amount of noise to generate apattern, so the stimulus is just a very small perturbation applied to the HSS. Forthreshold regimes, LPA bifurcation diagrams indicate a value which the perturba-tion must surpass in order to stimulate a response. The stimulus for a threshold66regime is given by ε , which is assigned the minimum value necessary to generatepolarization.During the simulations for some downward threshold regimes, even very smallperturbation values are not sufficient to destabilize the HSS. To avoid situationswhere ε ≤ 0, the width of the stimulus is increased from 5%. Note the biologicaldistinction between a wider pulse versus a larger pulse. A wider pulse refers to agreater area of the cell being stimulated, while a larger pulse refers to a stimulus ofgreater intensity.Bifurcation diagrams and simulations were conducted for all parameters. Rep-resentative diagrams or those with interesting/non-repetitive results are shown. Inorder to ensure the ultimate pattern that forms is not transient, PDE simulationswere conducted for longer durations than shown. Several parameter values andstimulus sizes were tested in each pattern forming regime.6.1 Goryachev ModelThe Goryachev model was initially developed as a detailed biochemical modelfor yeast, and exhibits Turing-type behaviour. The simplified model includes twocomponents: u represents active Cdc42, and v represents the total concentration ofinactive Cdc42. The system by Goryachev and Pokhilko [15] is given as follows:∂u∂ t = Du∂ 2u∂x2 +αEcu2v+βEcuv− γu, (6.4a)∂v∂ t = Dv∂ 2v∂x2 −αEcu2v−βEcuv+ γu, (6.4b)Ec =E0c1+∫S f (u)ds=E0c1+ζ . (6.4c)The final equation represents a total conserved amount of the Cdc24-Bem1complex. It is present in the autocatalysis terms as this complex responds in apositive feedback loop to increase activation of Cdc42 from the inactive cytoplas-mic form. In the Jilkine review article, Ec is taken to be a constant. The authorshypothesize that only a single peak of activation forms due to the conservation67Parameter Value DescriptionDu 0.0025 Diffusion of active form (µm2s−1)Dv 10 Diffusion of inactive form (µm2s−1)ζ 100 Scaling parameter (unitless)α 1 Conversion to activator via pathway one (µM−3s−1)β 2 Conversion to activator via pathway two (µM−2s−1)Ec 1/3 Conservation of total CDC42-Bem1 complex (µM)γ 0.01733 Conversion to inactive form (s−1)T 5 Conserved concentration (µM)L 10 Cell length (µm)Table 6.1: Parameter values for the Goryachev model in Equation 6.4. Valuesare calculated from the supplementary material in [15]. I chose values forDv, ζ , and Ec according to [29].feature. The two terms αEcu2v and βEcuv represent the two different pathways bywhich activation of Cdc42 occurs in budding yeast cells. Inactivation is representedby the term γu.Initial investigation of the Goryachev RDE model is shown in Figure 6.1 usingparameter values given in Table 6.1. A polarization response is elicited from smallamplitude perturbations around the steady state uSS. A narrow but tall peak formson the edge of the domain. Note the timescale: progression into patterning occursslowly and a stable polarization front occurs at about t = 700 seconds due to theslower [with respect to the Wave Pinning model] diffusion coefficient.In the line plot of Figure 6.1(c) it appears as though two peaks begin to form oneither edge of the domain. The leftmost peak grows faster, utilizing the availableresource v, and the right peak recedes (see Figure 6.1(d)). The universally stablebifurcation diagram, constructed with respect to the rate parameter α for the firstpathway, indicates no inherent patterning present in the reaction kinetics.68TimePosition of u  200 600 100005101030TimePosition of v  200 600 100005100.20.60 1 2 32345αConcentration of u0 5 103.63.844.24.44.6DomainConcentration of ut = 50, 100, 150, 200Stimulus (t = 0)0 5 1002040DomainConcentration of ut = 0t = 300, 500, 1000Stimulus45Stimulus (IC)a) b) c) d) tttFigure 6.1: Initial investigation of the Goryachev model with parameter val-ues from Table 6.1. (a) Kymograph of solutions for u (top) and v (bot-tom). (b) The bifurcation diagram for α is stable. (c) Line plots for theinitial stage of polarization, t = 50, 100, 150, 200 sec (times progressin direction of arrow). (d) Line plots for the final, stable polarizationpattern, t = 0, 300, 500, 1000 sec. While two peaks appear to forminitially, they resolve into one activation peak. The stimulus (red, rightaxis) is given by small amplitude noise about the steady state uSS.696.1.1 LPA system and bifurcation diagramsUsing LPA and the fact that Du Dv, Equation 6.4 can be reduced to a system ofthree ordinary differential equations (ODES) as follows:duldt= αEcu2l vg +βEculvg− γul, (6.5a)dugdt= αEcu2gvg +βEcugvg− γug, (6.5b)dvgdt=−αEcu2gvg−βEcugvg + γug. (6.5c)Since the Goryachev model also features conservation (vg = T −ug), I eliminate vgto obtain:duldt= αEcu2l (T −ug)+βEcul(T −ug)− γul, (6.6a)dugdt= αEcu2g(T −ug)+βEcug(T −ug)− γug. (6.6b)Equation 6.6 is the LP system that is used for constructing the bifurcation di-agrams in Figure 6.2. LPA bifurcation diagrams in Figure 6.2 are shown for theparameters (a) α (activation via pathway one) and (b) T (the conserved total con-centration of Cdc42, u+ v). Both feature only a well-mixed branch which is uni-versally unstable, indicating a Turing regime. While the LPA bifurcations do notindicate the number of stable peaks, from the results in the original publication weexpect only a single peak to be stable.6.1.2 Turing analysisAfter an initial attempt at a Turing analysis using Equation 6.4, I decided to firstnondimensionalize the model using the scalings:700 1 2 32345αConcentration of ul0 1 2 3 4 501234TConcentration of ula) b) Figure 6.2: LPA bifurcation diagrams for the LP Goryachev model in Equa-tion 6.6 with respect to parameters (a) α (activation via pathway one)and (b) T (the total concentration of Cdc42). Both bifurcation diagramshave only universally unstable, global branches (black, dash-dot lines).t =tˆγ , x =√Duγ xˆ, u =γβEcuˆ, v =γβEcvˆ,This produces the following nondimensional equations:∂ uˆ∂ tˆ =∂ 2uˆ∂ xˆ2 + ku2v+uv−u, (6.7a)∂ vˆ∂ tˆ = D∂ 2vˆ∂ xˆ2 − ku2v−uv+u, (6.7b)whereD =DvDu, k =αγβ 2Ec.With the standard parameters from Table 6.1, we have the default values of D =4000 and k = 1.31, and the identity u+ v = T becomes:T¯ = uˆ+ vˆ =(βEcγ)T = 1.90.71We will drop the “hats” for subsequent analysis.Our first step in the Turing analysis is to find the HSS values, (uss,vss) usingthe reaction terms in Equation 6.7. That is, we want f (uss,vss) = g(uss,vss) = 0.Using the fact that u+ v = T¯ , I find that the only relevant HSS is:vss =T¯ k+1−√(T¯ k+1)2−4k2k,uss = T − vss = T −T¯ k+1−√(T¯ k+1)2−4k2k.The Jacobian has the form:J =[fu fvgu gv]∣∣∣∣∣ss=[2kuv+ v−1 ku2 +u−2kuv− v+1 −ku2−u]∣∣∣∣∣ss.Evaluating J at the steady state (uss,vss) yields complicated equations which I donot list here, but notice that a = −c and b = −d. To ensure the HSS is stable, werequire that:1. Det(J) > 0,2. Tr(J) < 0.Since the Jacobian is singular, Det(J) = 0. For condition (2), Tr(J) < 0 whenT¯ > 1, which corresponds to T = 2.6 in the dimensional system, and k > 0. Wecan see from the LPA bifurcation diagram for T (Figure 6.2) that when T = 2.6 < 0,uss < 0. Since we always require uss > 0 (and therefore T = 2.6 > 0) and k > 0for biological relevance, we may say that Tr(J) < 0 for the entire domain we areconsidering.I introduce diffusion to the Jacobian, according to Equation 6.7, to test fordiffusion-driven instabilities:JD =[a−q2 cb d−Dq2],where a,b,c, and d are the entries of J. For Turing patterns to occur, we mustsatisfy the following two conditions:721. a+ dD > 02. ad−bc < D4(a+ dD)2.Condition (2) is always true, since ad− bc = Det(J) = 0, and the right hand sideof the inequality is positive. Solving the first condition will give values of theparameters k and T¯ for which q2 > 0. When solving for condition (1), we have twocases: solve for T¯ in terms of k, or solve for k in terms of T¯ .Case 1: Solving for T¯ gives the following range:1 < T¯ <−200+4001√10k200k,where 1/4000 < k < 4000 to ensure the right hand side of the inequality is greaterthan 1. Using the standard parameter value of k = 1.31, the range for T¯ becomes1 < T¯ < 54.5 which corresponds to 2.6 < T < 141.8 in the dimensional model.Case 2: When condition (1) is solved in terms of k, the parameter falls betweenthe two values:−8000T¯ +16008001±4001√−16000T¯ +160080018000T¯ 2.After substituting the standard parameter value T¯ = 1.90, the range becomes k ∈(0.0003,1102.4). This range decreases as T¯ increases, up to the limit T¯ = 1000 toensure the square root stays real. As T¯ decreases, the range increases.If we pick a specific values for T and k, there are multiple unstable modes, butonly at extreme values of the parameters.Connection to LPA Note that as we push the diffusion coefficients to their limits,the ranges from Case 1 and 2 become:1 < T¯ < ∞,0 < k < ∞,730 50 100 15000. Turing eigenvalue Du = 0.0025,   Dv = 10Du = 0.00025,   Dv = 100Figure 6.3: Maximal Turing eigenvalues as a function of T¯ , with the dimen-sionless parameter k = 1.31. Two different sets of diffusion coefficientsare plotted: the default values and a set with increased disparity (seelegend for values). As the diffusion rates approach the limiting val-ues, the range in which the maximal eigenvalue is positive increases to1 < T¯ <∞, matching the range for Turing instabilities predicted by LPA.The maximal Turing eigenvalues for Case 1 are plotted in Figure 6.3 with two dif-ferent sets of diffusion coefficients. Figure 6.3 shows how the range of values forT which admit a Turing instability increases as the disparity in diffusion rates in-crease. Under the given diffusion rates in Table 6.1, the Turing analysis predictsstability when T or k are increased. But, as the diffusion rates approach the lim-iting cases of zero and infinite diffusion, the Turing analysis begins to match thepredictions made by the LPA bifurcation diagrams in Figure 6.2. Note that T¯ = 1corresponds to T = 2.6, matching the value in Figure 6.2(b). In the LP bifurcations,our predictions are made based on extreme diffusion rates, so we do not expect tosee any stable regions.6.1.3 PDE simulations and connection to LPAThe LPA bifurcation diagrams and Turing analysis show the predominant patternforming mechanism to be a Turing instability. The Turing analysis also predicts a741000 30000510x2610200 600 100005102040a) b) c) d) 200 400 600 800 10000510Timex 1030200 400 600 800 10000510Time2040Figure 6.4: Kymographs showing different responses of the Goryachevmodel in Equation 6.4 for (a) α = 0.05, (b) α = 5, (c) α = 30, and(d) β = 0.07. The patterning regime of this model is almost universallyTuring, but note the difference in polarization times. Parts (c) and (d)show how even extreme parameter values affect the final pattern verylittle.stable region, but because this region occurs at such extreme parameter values, Iomitted the investigations from this thesis. PDE simulations were carried out forEquation 6.4, varying the parameter α , as shown in Figure 6.4. By taking α to bevery small, we can test the edge of the regime. With α = 0.05, a small stimulusstill results in pattern formation. However the length of time before a stable patternis achieved is greatly increased. In fact, a stable polarization response is not yetreached in Figure 6.4(a) at t = 4000.Varying α does not qualitatively change the polarization peak; it remains thesame height and width when the final pattern is achieved. It does, however, changethe length of time for this final pattern to be reached. When α is taken to belarger, as in Figure 6.4(b,c), time for polarization is reduced. The Turing analysispredicted stability at very extreme parameter values, and we see in Figure 6.4(c,d)that even with α = 30 or β = 0.07, the polarization pattern changes very little.Results (not shown) are similar for γ . The qualitative features of the peak do notchange, but the timescale for polarization shifts. Increasing the inactivation rate γtoo much leads to the inability of the Goryachev model to polarize; the active form75u is inactivated too quickly for any peaks to grow.6.2 Otsuji ModelThe Otsuji model is a two component conceptual model for polarization. Similarlyto the Wave Pinning model, the authors originally developed a six component RDsystem for active/inactive pairs of Rac, Cdc42, and RhoA. The simple conceptualmodel was created to test the ability of certain underlying features, namely massconservation and diffusion driven instability, to create stable activation peaks.In the publication [54], the authors focused on the following system:∂u∂ t = Du∂ 2u∂x2 +a1[v−u+ v(a2s(x)(u+ v)+1)2], (6.10a)∂v∂ t = Dv∂ 2v∂x2 −a1[v−u+ v(a2s(x)(u+ v)+1)2]. (6.10b)Here, the variable u represents a slowly diffusing active form and v representsthe quickly diffusing inactive form. The Otsuji model features mass conservationsince there is only interconversion between the two different forms. The kineticsterms can be rearranged into the more familiar form of u′ = (activation rate)·u -(inactivation rate)·v as follows:f (u,v) = a1[(u+ v(a2s(x)(u+ v)+1)2)u+(1(a2s(u+ v)+1)2−1)v]. (6.11)The parameter a1 can be modified to adjust the timescale of simulations. Astimulus is represented within the RDES, rather than in the initial conditions, bythe spatially dependent function s(x). For LPA and later simulations, we will makethe assumption that s(x)≡ s = 1. That is, we will assume s(x) is a constant, sincethe spatial component is removed in LPA. Note that in this conceptual model, v isassumed to affect both the activation and inactivation rates. Since v is inactive, thisassumption is not biological.Results for an initial investigation of the Otsuji model are shown in Figure 6.5.76Parameter Value DescriptionDu 0.1 Diffusion of active form (µm2s−1)Dv 10 Diffusion of inactive form (µm2s−1)a1 25 Timescale parameter (s−1)a2 0.7 Reaction parameter (µM−1)s 1.0 Intensity of stimulation (unitless)T 2 Conserved concentration (µM)L 10 Cell length (µm)Table 6.2: Parameter values for the Otsuji model in Equation 6.10. Valueswere taken from [29], where they were scaled from the original publi-cation in order to match the rates of diffusion with those in the WavePinning model. The stimulus, originally a spatially dependent functions(x), is reduced to a constant s.The kymograph and line plots use the parameter values from Table 6.2, and theoriginal form of the stimulus:s(x) = 1+0.1cos(2pi(x/L+0.5)), (6.12)where L is the length of the domain. The stimulus quickly resolves into a peak ofactivation in the center of the domain.With respect to the reaction kinetics, where the spatial aspect is not considered,the parameters s and a2 both have the same effect on the model. Therefore, ex-amining them independently yields bifurcation diagrams that are qualitatively thesame. In Figure 6.5(b), I consider the bifurcation parameter for the reaction kinet-ics as the combination a2s, and s is given a constant value. The reaction kineticsare purely stable.770 5 10048DomainConcentration of ut = 0t = 5, 10, 50Stimulus0.81.2Stimulus (Parameter s)TimePosition of u  100 200051026TimePosition of v  100 20005100.250.30.350 1 2012a2SConcentration of ua) b) c) tFigure 6.5: An initial look the Otsuji model from Equation 6.10 using pa-rameter values given in Table 6.2. (a) Kymographs for u (top) and v(bottom), with the stimulus quickly resolving into a peak in the centreof the domain. (b) Bifurcation diagram of the reaction kinetics with re-spect to a2s. (c) Line plot showing solutions of u for t = 0, 5, 10, 50sec (times progress in direction of arrow). For (a) and (c), the stimu-lus s (red, right axis) is given by a cosine function (see text), in (b) s isconsidered to be a constant.786.2.1 LPA system and bifurcation diagramsAs in the previous systems, the activator u has a local and a global form, and v hasonly a global form since Du Dv. The LP system is given by:duldt= a1[vg−ul + vg(a2s(ul + vg)+1)2], (6.13a)dugdt= a1[vg−ug + vg(a2s(ug + vg)+1)2], (6.13b)dvgdt=−a1[vg−ug + vg(a2s(ug + vg)+1)2]. (6.13c)Through conservation and subsequent simplifications with vg = T −ug, I have thefollowing system:duldt= a1[(T −ug)−ul +(T −ug)(a2s(ul +(T −ug))+1)2], (6.14a)dugdt= a1[(T −ug)−T(a2sT +1)2]. (6.14b)The LPA bifurcation diagram in Figure 6.6 is constructed with respect to theparameters a2s. There are two regions present: one in which the well-mixed branchis stable with an unstable local branch, and the other which occurs after crossingthe branch point (a2s = 0.5), leading to a change in stability for both branches. Thesecond region can be interpreted as a Turing instability, and small amounts of noiseshould give rise to spatial patterning. This region is where the default parametervalues are located, and where much of the analysis by the original authors wascarried out.The first region, where a2s < 0.5, is indicative of threshold behaviour. Ac-cording to the diagram, this region is stable to small perturbations but unstable tolarge perturbations. A perturbation extending past the unstable local branch shouldcause u to increase, apparently without bound. As one of the properties of the Ot-suji model is mass conservation, I expect u to reach a finite value. PDE simulations790 1 2012a2SConcentration of ulBP Figure 6.6: LPA bifurcation diagram for the LP Otsuji model in Equation 6.14,with the bifurcation parameter a2s. Legend: Well mixed branch (black);Local branch (red); Solid line (stable); Dash-dot line (unstable); BP =branch point. There are two patterning regions: threshold (a2s < 0.5)and Turing (a2s > 0.5).are needed to determine the long term pattern that forms, but the LPA bifurcationreveals a new region of pattern formation.6.2.2 Turing analysisOur first step in the Turing analysis is to find the HSS from the reaction terms inEquation 6.10. Recalling that u+ v = T , we have:vss =T(a2sT +1)2,uss = T − vss = T −T(a2sT +1)2.To simplify the analysis, let k = a2sT .The Jacobian evaluated at the steady state is:J =[fu fvgu gv]∣∣∣∣∣ss=k−1(k+1)3k(k2+3k+4)(k+1)3− k−1(k+1)3 −k(k2+3k+4)(k+1)3 .80Notice that our matrix entries have the form a =−c and b =−d. From the Turinganalysis description in Appendix C, there are two conditions we must satisfy inorder to ensure our HSS stable. However, we can see that the eigenvalues λ for Jhave the numeric values 0 and -1 by the following analysis.|J−λ I|= 0 =∣∣∣∣∣[a−λ bc d−λ]∣∣∣∣∣(6.16a)= (a−λ )(d−λ )− cb= (a−λ )(d−λ )−ad= λ (λ − (a+d)),where:a+d =k−1(k+1)3−k(k2 +3k+4)(k+1)3=−(k+1)3(k+1)3=−1.Therefore, we can conclude that the HSS is always stable. Now I introduce diffu-sion to the Jacobian to test for diffusion-driven instabilities:JD =[a−Duq2 cb d−Dvq2],where a,b,c, and d are the entries of J. For Turing patterns to occur, we mustsatisfy the following two conditions:1. aDu +dDv> 02. ad−bc < DuDv4(aDu+ dDv)2.As in the Goryachev model, condition (2) is always true since ad−bc = Det(J) =0. Solving the first condition will give values of the parameter k for which q2 > 0.81Substituting values for the diffusion rates Du and Dv according to Table 6.2and solving for condition (1) gives the interval k ∈ (1.09,7.74). Any values fork outside of this range, given these particular diffusion coefficients, will not havea positive value for q2. If we pick a specific value for k, plots of the maximaleigenvalue versus q reveal that multiple modes are unstable, and multiple peakscan form. This observation is confirmed in Figure 6.8(b,c), where two peaks arepresent at the beginning of the simulation.Connection to LPA Notice that if we push the diffusion rates to extreme values,the above range approaches k ∈ (1,∞). Turing eigenvalues showing the approachto this limit are plotted in Figure 6.7. Therefore, in the LPA bifurcations we donot expect to observe an upper limit for stability, but we do expect a lower limitcorresponding to k = 1. To match the regions of instability predicted by the Turinganalysis and by the LPA bifurcation in Figure 6.6, recall that k = a2sT and ourbifurcation parameter is the combination a2s, with T = 2. At k = 1, we have thevalue a2s = 0.5 when T = 2, which is exactly where the LP system predicts achange in stability.6.2.3 PDE simulations and connection to LPAThe PDE simulations in Figure 6.8 show the long term patterns for parameter val-ues in each of the two regions. I continue to consider s as a parameter with constantvalue as opposed to being spatially variable. Simulations for a value of s = 0.67confirm the region is stable to small perturbations (not shown). When a perturba-tion of ε = 8 is applied, a clearly defined polarization pattern results. This result isshown in Figure 6.8(a). Further reducing the value of s or a2 requires much largerperturbations for pattern formation, as was suggested by the LPA bifurcation dia-gram. In addition, due to conservation in the Otsuji model, the peak of activationreaches a finite value which remains stable rather than increasing without bound.The next set of kymographs, in Figure 6.8(b,c) show PDE simulations in theTuring regime. By manipulating the values for s and a2, while still remainingin the Turing regime, patterns with two or even three peaks can be initiated. Asmentioned by the authors of the model in [54], multiple peaks are unstable and820.9 1 1.1−505101520 x 10−3kMaximal Turing eigenvaluea) b) 0 10 20 3000. Turing eigenvalueDu = 0.1,   Dv = 10Du = 0.001,   D v = 1000Figure 6.7: Maximal Turing eigenvalues as a function of k (recall that k =a2sT ) for the Otsuji model. Two different sets of diffusion coefficientsare plotted: the default values and a set with increased disparity (seelegend for values). As the diffusion coefficients approach the limitingvalues, the range in which the maximal eigenvalue is positive increasesto 1 < k < ∞, matching the range for Turing instabilities predicted byLPA. Part (a) shows a larger view of the domain k, while part (b) focuseson the convergence to k = 1.disappear on a short timescale. This property is consistent among our simulations.As evidenced by Figure 6.8(c), the length of time for the second peak to disappearis dependent on the parameter settings. The parameter a2 has increased slightlyfrom part (b), which increases the lifetime of the second peak. For biologicalrelevance, the combination of a2s must be restricted to a region in which the secondpeak disappears quickly. The Turing analysis predicts that a2sT > 7.74 will be astable region. In part (d), when a2 and s are increased to satisfy the inequality,stability occurs.6.3 Meinhardt Activator-Substrate Depletion ModelThe type of cell for the next model deviates from what we have previously con-sidered. For the Wave Pinning, Goryachev, and Otsuji models, our cells were eu-karyotic. Meinhardt and de Boer [46] introduce the Activator-Substrate Depletion83100 200 3000510Timex 24100 200 3000510510a) b) x  100 300 500051026c) 500 1000 15000510Timed) 22.12.2Figure 6.8: Kymographs of the Otsuji model in Equation 6.10 correspondingto the patterning regimes mapped in Figure 6.6. Parameter values arefrom Table 6.2, except (a) s = 0.67 (threshold regime, upward perturba-tion) (b) a2 = 0.7 (Turing regime) and (c) a2 = 1.2 (Turing regime) (d)a2 = 1.7,s = 2.5, which corresponds to k = a2sT = 8.5 (Stable regime).Part (a) confirms that the peak of activation reaches a stable value, ratherthan increasing without bound. In the Turing regime, multiple peaksmay form initially, but they always resolve into one peak. The lengthof time for the other peaks to diminish is dependent on the parame-ters. The final regime confirms that when the parameter combinationk = a2sT > 7.74, the model becomes stable.Model (ASDM) for the bacteria Escherichia coli. Recall from Chapter 3 that bacte-rial cells are too small to sense spatial gradients, so this is not a chemotaxis model.Instead, the RDES model the correct placement of the division ring within a divid-ing bacterial cell. The following RD equations:∂u∂ t = Du∂ 2u∂x2 + svu21+ sau2− rau, (6.18a)∂v∂ t = Dv∂ 2v∂x2 +bb− svu21+ sau2− rbv, (6.18b)correspond to the membrane bound and free forms (u and v respectively) of theFtsZ protein. Note that we still have Du Dv, since FtsZ is located either in the84Parameter Time (scaled by 25) DescriptionDu 0.01 Diffusion of active form (µm2s−1)Dv 10 Diffusion of inactive form (µm2s−1)s 0.1 Conversion rate (µM−2s−1)sa 0.1 Half saturation concentration (µM−2)ra 0.1 Decay of active form (s−1)bb 0.15 Production of inactive form (µM · s−1)rb 0.01 Decay of inactive form (s−1)L 10 Cell length (µm)Table 6.3: Parameter values for Meinhardt’s activator substrate depletionmodel in Equation 6.18. The original values for the model can be foundin [46]. The values given in this table were scaled (see the supplemen-tary material in [49]) so that the diffusion rates coincided with those inthe Wave Pinning model.cytoplasm or on the membrane, so LPA will still apply. The full system involvesequations for membrane-associated and free forms of MinD and MinE. Since theinitial symmetry breaking event occurs through FtsZ, it is these equations alonethat I will consider as the polarization model.This system follows an activator-substrate depletion mechanism. Species ufunctions as an activator, which effectively depletes the substrate v. Parametersra and rb represent decay of u and v, and the production rate of v is given by bb.Through the inclusion of production and decay terms, this model does not satisfyconservation. Positive feedback to activation is given by the term sv u21+sau2. It issuggested that the parameter rb can be given a value of 0, but here I have decidedto give it a small value of 0.01.A preliminary investigation of Meinhardt’s ASDM is shown in Figure 6.9, usingthe parameters from Table 6.3. The stimulus from the original publication is givenas small amplitude noise around the steady state uSS. In parts (a) and (b), the850 5 10024DomainConcentration of uTimePosition of u  100 300 500051024TimePosition of v  100 300 50005100.70.80.9t = 0t = 100, 150, 250Stimulus1.451.55Stimulus (IC)0 0.5 101234s aConcentration of A0 0.1 0.2012raConcentration of AH a) b) c) d) tFigure 6.9: Initial analysis of Equation 6.18, Meinhard’s ASDM with param-eters from Table 6.3. (a) Kymograph of numerical PDE solutions for u(top) and v (bottom). (b) Line plot of solutions for u at t = 0, 100, 150,and 200 sec (times progress in direction of arrow). The stimulus (red,right axis) is given by small amplitude noise around the steady state uSS.The bottom two plots are bifurcation diagrams of the reaction kineticsfor (c) ra (instability due to Hopf at ra = 0.13) and (d) sa (stable).kymographs and line plot feature a stable polarization pattern at the edge of thedomain. In the bifurcation diagram for ra, we see that a Hopf bifurcation actuallygives rise to an instability in the reaction kinetics. The kinetics bifurcation for sa ispurely stable.866.3.1 LPA system and bifurcation diagramsApplying LPA to the Meinhardt RDES yields the following system of equations:duldt= svgu2l1+ sau2l− raul, (6.19a)dugdt= svgu2g1+ sau2g− raug, (6.19b)dvgdt= bb− svgu2g1+ sau2g− rbvg. (6.19c)Again we have that u takes on both local and global forms as ul and ug, and v ispurely global. Note here that, unlike the polarization models discussed until now,the Meinhardt model does not satisfy conservation. Thus we must use the full threedimensional system for our LPA bifurcation analysis.The location of the default parameter values in the LP bifurcation diagramsconfirms that pattern formation in the original model is indeed due to a Turinginstability. There are two different regions in the bifurcation diagram with respectto parameter sa. The first region, to the left of sa = 0.45, is characterized by anunstable well-mixed branch and a stable local branch. To the right of sa = 0.45,both branches switch stability. This form suggests the possibility for thresholdbehaviour in the full PDE simulations. A downward perturbation of u past thethreshold defined by the local branch may give rise to polarization.The bifurcation diagrams for parameters ra and bb contain three possible re-gions of different behaviour. Two are similar to those described above, and a thirdregion is defined after crossing a Hopf bifurcation. The first region is a thresh-old regime where bb > 0.32 or ra < 0.046. From either 0.10 < bb < 0.32 or0.046 < ra < 0.13 is a Turing regime. Then, for bb < 0.10 or ra > 0.13, we cross aHopf bifurcation, which is present on both the global and the local branches. In thecase of both parameters, the well-mixed curve is unstable on either side of the Hopfpoint. The LPA diagrams do not provide further information on how these regionsmay differ in pattern formation, and full PDE simulations will be used to reveal anydifferences. The diagrams with respect to parameters s and rb (not shown) do not870 0.2 0.4 0.6051015bbConcentration of Al H LP BP H 0 0.05 0.1 0.15 0.2051015raConcentration of AlH LP BP H 0 0.5 101234s aConcentration of AlBP a) b) c) Figure 6.10: LPA bifurcation diagrams for the Meinhardt ASDM in Equa-tion 6.19. Diagrams are constructed for parameters (a) bb, (b) ra, and(c) sa. Legend: Well mixed branch (black); Local branch (red); Solidline (stable); Dash-dot line (unstable); BP = branch point; H = Hopfbifurcation. Regions involved are: downward threshold regimes, Tur-ing regimes, and a threshold regime with an instability due to a Hopfbifurcation in the reaction kinetics.8850 100 150 2000510Timex50 100 150 2000510Time 50 100 150 2000510Time2 4 6a) b) c) Figure 6.11: PDE simulations of u for Equation 6.18 testing the regionsmapped in Figure 6.10(b). For this set, three values of ra are testedin two different regions: (a) ra = 0.03 (threshold region, stable re-sponse), (b) ra = 0.15, and (c) ra = 0.23 (both threshold, unstabledue to Hopf). Despite large downward perturbations, polarization wasnot maintained in the threshold region (a). Parts (a) and (c) featurethe second threshold regime with Hopf instabilities. Notice that thelower state is stable, but the pulse oscillates. Polarization in the Turingregime is not shown here (but see Figure 6.9(a)).reveal any other possibilities for behaviour and do not contain local branches.The well-mixed branches in the LPA diagrams of bb and ra fold underneaththemselves and run close to or along ul = 0. This behaviour suggests an additionallow steady state value of active FtsZ. Each of these diagrams also contain limitpoints, and it was expected that the well-mixed branches would become stable aftercrossing these points. The presence of the Hopf bifurcations prevent this change instability from appearing in the bifurcation diagrams.6.3.2 PDE simulations and connection to LPAWhile most of the authors’ simulations on this model were for a regime in whichpolarization was initiated through a Turing instability, the LPA bifurcations revealpossibilities for different behaviour. Long-term PDE simulations for the RD systemin Equation 6.18 were conducted to match regions in both Figure 6.10(b) and (c)(that is, parameters ra and sa). See the kymographs in Figure 6.11 and Figure 6.12,respectively. As expected from the bifurcation diagrams, simulations for the pa-rameter bb were similar to those of ra and are not shown.The parameter ra is the decay rate of u. Its bifurcation diagram features threepatterning regimes, two of which are tested in Figure 6.11. The first of these,89part (a), is a threshold regime, requiring a downward perturbation. Despite testingdownward perturbations with large width (up to 30% of the domain) and large mag-nitude, polarization is not maintained, and the solution becomes spatially uniform.While it may still be possible to elicit polarization which is maintained, withinbiologically relevant values this region does not permit a polarization response.The second region reveals Turing patterns (not shown, but see Figure 6.9) inresponse to a small stimulus. The third region is shown in Figure 6.11(b,c), and ap-pears unstable due to a Hopf bifurcation. The bifurcation diagram for ra indicatesa threshold-driven polarization response that must increase in magnitude as ra in-creases. Given a large stimulus, ra = 0.15 with ε = 2 and ra = 0.23 with ε = 6 forparts (b) and (c) respectively, polarization does indeed occur. Interestingly in thisregion, due to the Hopf bifurcation, the pulse itself exhibits transient oscillations.PDE simulations with respect to changes in sa, the half saturation constant forautoactivation, are shown in Figure 6.12. Parts (a) and (b) correspond to a Turingregime. Here it is shown that with sa = 0.05, only one peak of activation forms.Increasing to sa = 0.3 results in the formation of two activation peaks, both ofwhich are stable when simulation times are extended. Part (c) corresponds to thesecond patterning regime predicted by the bifurcation diagram in Figure 6.10(c):a downward threshold regime. While the attempt to initiate polarization in thecorresponding regime for ra failed, here we see that a polarization response results.A stimulus of ε = 0.05 was applied to 5% of the domain, as per the standardstimulus in Equation Pseudopod Centered ModelThe “Pseudopod Centered” RDES were originally developed as an ODE system in[45] to model the extension and retraction of filopods in nerve cell growth cones(see Chapter 2 for biological definitions). It was adapted by [52] as the RDES for achemotaxis model that operates through pseudopod-based feedback. A local acti-vator a and inhibitor c are localized to the pseudopods, whereas a global inhibitorb rapidly diffuses throughout the cell. The PDE version of the system from [52] is90200 400 6000510200 400 6000510x 1226200 400 6000510Timex0.511.5a) c) b) TimeFigure 6.12: Kymographs of u for the Meinhardt RDES in Equation 6.18 test-ing regions in Figure 6.10(c). Parameters are given in Table 6.3, and sais varied as follows: (a) sa = 0.05 (Turing regime producing one peak),(b) sa = 0.3 (Turing regime producing two peaks), and (c) sa = 0.6(Threshold regime, downward perturbation). Notice that the Turingregime may give rise to two stable peaks, as in part (b).as follows:∂a∂ t = Da∂ 2a∂x2 +s(a2/b+ba)(sc + c)(1+ saa2)− raa, (6.20a)∂b∂ t = Db∂ 2b∂x2 +rb|Γ(t)|∫Γ(t)adx− rbb, (6.20b)∂c∂ t = Dc∂ 2c∂x2 +bca− rcc. (6.20c)The first reaction term in Equation 6.20(a) represents autocatalysis of the ac-tivator. Notice that c is present in the denominator; an increase in local inhibitorc effectively leads to a decrease in the local activator a at that particular location.Peaks of activation are not stable in this model, as they continually pop up andthen decay due to local inactivation from c. The global inhibitor b contributes to alowered autocatalysis rate of a.91Both inhibitors experience production rates proportional to the amount of apresent, as evidenced by the first terms of the respective reaction kinetics. That is,with more local activator, both inhibitors increase production. Increase in the localinhibitor c occurs due to local concentrations of a. The global inhibitor b diffusesquickly, and so production rate is proportional to the global accumulation of a onthe membrane, given by the integral. All three components experience decay, givenby raa, rbb, and, rcc.The stimulus is given by the function:s(x) = ra[1+dy cos(2pi(x− j)L)](1+drRnd). (6.21)There is both an asymmetric gradient component representing an external chemoat-tractant as well as noise in that gradient. The strength of the external stimulus isgiven by dy and the direction is represented by j, where j is a location on theboundary of the cell. Since L = 1, we have that j ∈ [0,1]. Rnd is a uniform randomnumber on the interval [0,1] which represents random fluctuations in the externalstimulus; the magnitude of noise is given by dr. See Table 6.4 for parameter values,and note that the stimulus is constantly present in the Pseudopod model.Preliminary examination of the Pseudopod Centered model is shown in Fig-ure 6.13. Parameter values from Table 6.4 and the original form of the stimulus(above) were used to generate the simulations. The kymographs for part (a) show aregion of activation in the center of the domain for the local activator (top) and in-hibitor (bottom). This region features peaks that are not static, but continually riseand fall. The kymograph of the global inhibitor (center) is homogeneous through-out the domain.Figure 6.13(b) features a line plot of the concentration of a, the local activator,at t = 0, 600 and 1000 sec. Here, I featured t = 600 sec in blue, because the dis-tinction between this and the line for t = 1000 sec is more clear. One can see howthe peaks are not stable, but change for different snapshots in time. The kineticsbifurcation diagrams in part (c) were constructed with respect to the average valueof the stimulus s and the decay rate of the local inhibitor rc. The reaction kinet-ics are purely stable; the introduction of diffusion is responsible for the polarizedpattern we see.92Parameter Value DescriptionDa 10−9 Diffusion of local activator (µm2s−1)Db 1 Diffusion of global inhibitor (µm2s−1)Dc 10−8 Diffusion of local inhibitor (µm2s−1)ra 0.02 Decay rate of the activator (s−1)ba 0.1 Production rate of the activator (µM)sa 0.005 Saturation of activator autocatalysis (µM−2)rb 0.03 Production and decay rate of global inhibitor (s−1)bc 0.005 Production rate of the local inhibitor (s−1)rc 0.013 Decay rate of the local inhibitor (s−1)sc 0.2 Michaelis-Menten constant (µM)s 0.0201 Average value of stimulus (µM · s−1)dy 0.02 Stimulus, external asymmetry (unitless)dr 0.01 Stimulus, magnitude of random noise (unitless)j 0.5 Stimulus, direction (µm2)L 1 Cell length (µm)Table 6.4: Parameter values for the Pseudopod Centered model in Equa-tion 6.20. Values are drawn from [45] and [52]. Using the diffusionrates for the PDE model resulted in a polarization response that froze,not the results expected. Through personal correspondence with authorJohn A. Mackenzie, I learned that the deforming boundary domain (afeature of their model) is an important aspect in unfreezing the pattern,given these diffusion rates. In order to get continually changing activa-tion peaks similar to Meinhardt’s ODE model in [45], he suggested toincrease the disparity in diffusion rates between the locally and globallyacting species.930 0.02 0.04 0.0601234567sConcentration of a0 0.02 0.040123456rcConcentration of aTimex  500 1000 1500 200000.512468Timex  500 1000 1500 200000.512.12.2Timex  500 1000 1500 200000.511230 0.5 10246810DomainConcentration of at = 0t = 600t=1000Stimulus (s)0.02050.0185a) b) c) d) StimulusFigure 6.13: Preliminary analysis for the Pseudopod Centered model inEquation 6.20 using the parameter values in Table 6.4. The stimulus isgiven by Equation 6.21, which represents an external chemoattractantwith noise. (a) Kymographs of solutions for a, b, and c from top tobottom respectively. The slowly diffusing variables a and c show con-tinually changing peaks of activation, while b is homogeneous acrossthe domain. (b) Line plot of solutions for a at t = 0, 600, and 1000sec (see legend) along with the stimulus (red, right axis). (c) Kineticsbifurcation plot for s, the average value of the stimulus. (d) Kineticsbifurcation plot for rc, decay rate of the local inhibitor c. In (a) and(b), the stimulus is given in Equation 6.21. In (c) and (d), the stimulusis a constant whose default value is given in Table 6.4.946.4.1 LPA system and bifurcation diagramsSince we do not consider the entire boundary of the cell for LPA, I will make thefollowing simplification:rb|Γ(t)|∫Γ(t)adx≈ rba.For the LPA approximation, I will remove the spatial aspect of the stimulus s(x)and instead use the average value of this function. Taking the average gives theconstant value s(x) ≡ s = 0.0201. With these modifications, and assuming that aand c are slow diffusing while b is fast diffusing, the LPA system for the PseudopodCentered model is given by the following equations:daldt=s(a2l /bg +ba)(sc + cl)(1+ saa2l )− raal, (6.22a)dagdt=s(a2g/bg +ba)(sc + cg)(1+ saa2g)− raag, (6.22b)dbgdt= rbag− rbbg, (6.22c)dcldt= bcal− rccl, (6.22d)dcgdt= bcag− rccg. (6.22e)For the LPA bifurcation diagrams, we are again working with the parameters sand rc, featured in Figure 6.14(a) and (b) respectively. We see from both diagramsthat the default parameters give rise to a Turing instability.In the first diagram, there are essentially two different regions: Turing andstable. There is a local branch which is present in the Turing regime, but, as wehave seen before, its presence can be neglected. A region of stability can be foundwhen s < 0.0033, but its size is negligible and it can also be neglected. The Turingregime spans essentially from s ≈ 0 to the Hopf bifurcation at s = 0.60. Whens > 0.60, the Pseudopod Centered model is stable to any perturbations.In part (b), there are also two patterning regimes represented: a Turing regime950 0.02 0.04 0.060246sConcentration of al0 0.02 0.040246rcConcentration of alBP BP BP BP H H H H H H a) b) Figure 6.14: LPA bifurcation diagrams for the Pseudopod Centered LPAmodel in Equation 6.22. Parameters are (a) the average value of thestimulus s and (b) the decay rate of the local inhibitor rc. Legend: Wellmixed branch (black); Local branch (red); Solid line (stable); Dash-dot line (unstable); BP = branch point; H = Hopf bifurcation. Thereare three regimes in these diagrams: Turing regime, purely stable, anddownward threshold response.and a threshold regime. Again, there is a very small stable regime to the left of theHopf bifurcation on the well-mixed branch, but as it occurs at essentially s = 0, itcan be neglected. The Turing regime extends from rc ≈ 0 to the rightmost branchpoint at rc = 0.027. In this bifurcation diagram, the local branch extends into theregion where the well-mixed branch becomes stable, creating a possible polariza-tion response when a large enough downward perturbation is applied. Therefore,when rc > 0.027, we have a downward threshold regime.6.4.2 PDE simulations and connection to LPAI provide simulations for changes in both the parameters s and rc in Figure 6.15and Figure 6.16 respectively. The Pseudopod Centered model gives rise to inter-esting dynamic patterns, so several different values within the Turing regime arevisualized for this section. The stimulated 5% of the domain is in the center rather962000 400000.51x2000 400000.5105102000 4000 600000.51Timex2000 4000 600000.51Time66.57a) b) c) d) Figure 6.15: Kymographs of the Pseudopod Centered model in Equa-tion 6.20, testing the regimes mapped in Figure 6.14(a). The parameters is modified as follows: (a) s = 0.01 (b) s = 0.02 (c) s = 0.06 (all Tur-ing regime), and (d) s = 0.07 (stable). While the parameters in (c) stillconstitute a Turing regime, the pattern is very close to being stable. Onprogressing to (d), our pattern succumbs to stability. Notice the longersimulations times for (c) and (d). Other parameter values are given inTable 6.4.than the edge. Changing the location of the stimulus provides a more interestingview of the dynamics for this model, but the results still hold when the stimulus ispositioned at the standard edge of the domain.For parts (a), (b), and (c) in Figure 6.15, polarization patterns result from Tur-ing instabilities. The parameter values are s = 0.01, 0.02, and 0.06. At lower valuesfor s (part (a)), the activation peaks remain very localized to the center, in the di-rection of the small perturbation. As s increases (part (b)), the patterns becomemuch more dynamic, and activation peaks spread throughout the domain. Further972000 400000.51x2000 400000.512000 400000.51Timex2000 400000.51Time020400510150505a) b) c) d) Figure 6.16: Kymographs of the Pseudopod Centered model in Equa-tion 6.20, testing the regimes mapped in Figure 6.14(b). Parametervalues are given in Table 6.4, except for rc which takes the values:(a) rc = 0.005 (b) rc = 0.01 (c) rc = 0.02 (all Turing regime), and (d)rc = 0.035 (threshold, downward perturbation). Notice the differentdynamics within the Turing regime. The peaks start very tall and tran-sient, but they eventually work towards a shorter height and longerduration. In (d), a downward perturbation produced a stable trough ofthe active form a for the threshold regime.increasing s (not shown) leads to patterns that are similar again to part (a). Whens is very close to the Hopf bifurcation (as in part (c)), very small peaks result andform repeating oscillations. In part (d), s has passed into the stable regime. Verysmall oscillations result from the stimulus, but unlike part (c), these oscillationseventually die away. Note the longer simulation times for parts (c) and (d).Figure 6.16 follows a similar progression through the two regimes as Fig-ure 6.15, with the first three kymographs representing different values within theTuring regime. The parameter rc progresses from 0.005 to 0.01 to 0.02 for parts98(a-c). In response, the peaks initially are very tall but are sparse and present fora short time. They decrease in height while increasing in width and duration asrc increases. At rc = 0.02 in part (c), the peaks are actually stable and the pat-tern freezes. In part (d), rc = 0.035 is in the threshold regime. A perturbation ofε = 4.3, which reaches just below the threshold defined in Figure 6.14(b), inducesa stable peak of activation.6.5 Balanced Inactivation ModelThe Balanced Inactivation model by Levine et al. [33] was developed to model notjust an amplification of the external chemical gradient, but a switch-like responseleading to polarization. This response causes distinct front and back regions in thecell. The model was developed with a local excitation / global inhibition (LEGI)design; there are no inherent pattern forming regions. Instead, patterning occursonly when a stimulus is present. There are three abstract components to this model:a membrane-bound activator A, and an inhibitor that has membrane-bound Bm andcytosolic B forms. The model is given by the following RDES:∂A∂ t = kaS(x)− k−aA− kiABm, (6.23a)∂Bm∂ t = kbB− k−bBm− kiABm, (6.23b)∂B∂ t = D∂ 2B∂x2 , (6.23c)with the boundary condition for the outward pointing normal derivative of B:D∂B∂n = kaS(x)− kbB. (6.24)Notice the units for the variables (see Table 6.5 for parameter units). A and Bm,which are located at the boundary, have units of molecule/µm. B, which is foundthroughout the domain, has units of molecule/µm2. Taking note of the units isimportant for modifications I later make to this model.In the presence of a stimulus, the term kaS(x) represents equal production rates99Parameter Value Descriptionka 1 Production of A (s−1)k−a 0.2 Spontaneous degradation of A (s−1)ki 10 Mutual inactivation of A and Bm (µm · s−1molecule−1)kb 3 Production of B (µm · s−1)k−b 0.2 Spontaneous degradation of B (s−1)D 10 Diffusion of B in cytosol (µm2s−1)L 10 Cell length (µm)S¯ 10 Average value of stimulus (molecule ·µm−1)p 0.2 Characterizes gradient steepness (dimensionless)Table 6.5: Parameter values for the Balanced Inactivation model Equa-tion 6.23. Here I modify ki from the original paper by Levine et al. [33]to have a value of 10 µm · s−1molecule−1.of A and B. The function S(x) is the stimulus, and approximates a gradient ofchemoattractant. The free form of the inhibitor, B, becomes membrane bound, Bm,at the rate kb. Once the inhibitor is located at the membrane, it can inactivate Athrough the term kiABm. Degradation of A and Bm are given by the terms k−aAand k−bBm respectively. Since A and Bm are located on the cell membrane, it isassumed that diffusion occurs so slowly in comparison with cytosolic diffusionthat it can be ignored. The authors note that the results are robust when membranediffusion (D≈ 1µm2/s) is included. Cytosolic species B diffuses at a rate given bythe diffusion coefficient D.For the original publication, Levine et al. [33] model their stimulus after agradient of chemoattractant. The form of their stimulus S(x) is a straight, gradedline with the values: S f = S¯(1+ p) at the front of the cell and Sb = S¯(1− p) atthe back. The parameter S¯ is the average of the stimulus, and p is a dimensionlessparameter corresponding to gradient steepness. Spoiler alert: we will see that the100general form for the stimulus (i.e., a narrow width, arbitrarily tall, and transientperturbation) that was applied to all previous models does not generate a sustainedpolarization response in the Balanced Inactivation model. Therefore, I will use thestimulus S(x) for most simulations.6.5.1 Model modificationsIn the Balanced Inactivation model, the interpretation for approximating the cell to1D differs from the previous models. In order to properly compare the models, inthis section I modify the original geometry of the Balanced Inactivation model tomatch the 1D transect geometry of the other models discussed so far. As discussedbelow, both the PDES and the boundaries must be modified. See Figure 6.17(a,b)for a visual depiction of where A, Bm, S(x), and B are originally located. Themembrane-bound species A and Bm are found only on the membrane (boundary)where they can interact with each other and the stimulus S(x) through the func-tions presented in Equation 6.23(a,b). Cytosolic species B is found in the reactionequation for Bm because it is free to diffuse anywhere on the cell, including to themembrane. But B is the only chemical species found on the interior of the cell. Itdoes not react here because it is inactive and there is no stimulus. Therefore, theonly place it can react is at the membrane, or boundary, hence the reason that theauthors placed the reaction term for B in the boundary conditions.From previous models, we are used to “no flux” boundary conditions, whichwould appear as Bx(0, t) = 0 = Bx(L, t). In the Balanced Inactivation model, thereaction for B takes place on the membrane (in the boundary conditions) ratherthan in the actual RDE. In its current form, the model is incompatible with LPA. Tounderstand why, consider what happens when we make the LP approximation. Thespatial aspect of the model is removed, and therefore the reaction kinetics for B arenot considered, since they occur on the boundary of the cell. The LP equation forBg becomes B′g = 0 and any reaction information is lost. We want to reformulatethe Balanced Inactivation model in a way that follows the same simplified view ofa cell that the other models follow (see the schematic diagram in Figure 6.17(c)).These modifications will allow us to perform LPA on the Balanced Inactivationmodel with no lost information.101a) b) c) BA, Bm, S(x)1D approximationBoundaryzonesBFrontBackFrontBackB A, Bm, S(x)FrontBackA, Bm, S(x)A, Bm, S(x)Figure 6.17: Schematic diagram for the Balanced Inactivation model. Leg-end: cytoplasm = coloured, membrane = white. (a) A cross sectionalview of a cell, marking the 1D approximation the authors take forconstructing the Balanced Inactivation model. Dashed lines indicatethe 1D domain, and circled in blue are the regions of interaction be-tween B and A,Bm (indicated as boundary zones). (b) A schematicview of the 1D domain. Non-diffusing A, Bm, and S(x) interact witheach other across the domain (indicated by the dashed line). B dif-fuses everywhere, but only interacts with the non-diffusing species atthe boundaries (indicated by the points). (c) The simplified view of thecell after reformulating the model in the same manner as the Wave Pin-ning model. A thin slice of cytoplasm contains B, and is sandwichedbetween the membrane housing A, Bm, and S(x). Interactions for allspecies occur across the domain.I will assume the following approximation from the previous models. A thinstrip of cytoplasm is squished between the membrane, as in a lamellipod, leadingto reactions between membrane bound and cytosolic species occurring throughoutthe domain. The only difference between the two classes now is the diffusion rates.Instead of having the reaction terms in the boundary conditions, we move them tothe RDE for B to get the following equation:∂B∂ t = D∂ 2B∂x2 +kaS(x)L−kbBL. (6.25)102Now, A, Bm, S(x), and B can interact with each other throughout the entire domain.Note that the reaction rates have been divided by the cell length, L. Previously, Bonly reacted at the membrane. Now, since we consider B to react across the entirelength of the 1D cell, the rates need to be adjusted to be consistent with the originalformulation of the model. To clarify, recall that the units for B are molecule/µm2.Dividing by the length L(µm) ensures the units for each term are consistent.The original stimulus was given as a graded line, so we can define our stimulusto have the form S(x) = s0 + s1x. I will take S¯ = 10, s0 = S¯(1− p), and s1 = 2p.These values will result in a line with the endpoints S(0) = Sb and S(L) = S f ,consistent with the stimulus from [33] that I described above. The parameter p canbe modified to change the steepness of the gradient.An additional consideration is the new boundary conditions. To find these, Ifirst solve for the HSS of B using Equation 6.25 with S(x) = s0 + s1x:0 =kas0L+kas1xL−kbBSSL,BSS =kas0kb+kas1xkb.We can find the boundary conditions by taking the first derivative with respect to xand evaluating at the endpoints:∂BSS∂x∣∣∣∣x=0,L=kas1kb.With these modifications, I arrive at the following modified equations for theBalanced Inactivation model:∂A∂ t = kaS(x)− k−aA− kiABm, (6.27a)∂Bm∂ t = kbB− k−bBm− kiABm, (6.27b)∂B∂ t = D∂ 2B∂x2 +kaS(x)L−kbBL, (6.27c)103with constant flux boundary conditions for B:∂B∂x∣∣∣∣x=0=kas1kb=∂B∂x∣∣∣∣x=L, (6.28)and the stimulus given by the function S(x) = s0 + s1x where:S¯ = 10,s0 = S¯(1− p),s1 = 2p.Note that we had “no flux” or “periodic” boundary conditions on the previousmodels, and initially a constant flux boundary condition may feel unsettling. In thecontext of the original version of the model however, recall the original boundarycondition for B in Equation 6.24. In this equation, B is created at the boundary.Therefore, it is understandable for B to have constant flux boundary conditions.There are still “no flux” conditions on A and Bm.Initial investigation was done using the modified Balanced Inactivation modelin Equation 6.27 with the constant flux boundary conditions in Equation 6.28.Since the Crank Nicolson method assumes zero flux boundary conditions, I makeuse of the Thomas Algorithm (outlined in Appendix D) for my numerical simu-lations. Here, the front of the cell occurs at position x = 10, and the stimulus ispresent for the duration of the simulation. From the kymographs and line plotsin Figure 6.18(a,c,d), polarization is quickly established and maintained, with Alocalizing to the front and Bm localizing to the back of the cell. The free form Bkeeps an almost uniform concentration by quickly diffusing through the cell. Thebifurcation diagram in Figure 6.18(b) indicates purely stable reaction kinetics withrespect to the inactivation parameter ki.6.5.2 LPA system and bifurcation diagramsWith the modifications given above, and assuming that A and Bm are the slowvariables while B is the quickly diffusing variable, I arrive at the following LP104TimePosition of A  20 40 60 80051028TimePosition of Bm  20 40 60 80051028TimePosition of B  20 40 60 80051033.540 5 1002468DomainConcentration of At = 0t = 5, 10, 25Stimulus0 5 1002468DomainConcentration  AB mB91113a) b) c) d) 0 5 10012345kiConcentration of AtFigure 6.18: Initial investigation for the modified Balanced Inactivationmodel in Equation 6.27 with parameters from Table 6.5. (a) Kymo-graphs of solutions for all three components: A (top), Bm (centre), andB (bottom). (b) Bifurcation diagram of the reaction kinetics using theinactivation parameter ki. The bifurcation diagram is completely sta-ble. (c) Line plot of solutions for A at t = 0, 5, 10, 25 sec (timesprogress in direction of arrow) with the stimulus (red, right axis). (d)Line plot of solutions for A (black, solid), Bm (red, solid), and B (red,dot-dash) at t = 80 sec. Profiles for A show a switch-like response:the slope is shallow at the back of the cell, while the steepness greatlyincreases for the front half.1050 5 10012345kiConcentration of Al0 5 1001234SConcentration of Ala) b) Figure 6.19: LPA bifurcation diagrams for the Balanced Inactivation systemin Equation 6.29 using bifurcation parameters: (a) ki and (b) S. Noticethat the entire region is stable for both parameters, and there is only aglobal branch.system:dAldt= kaS− k−aAl− kiAlBml, (6.29a)dAgdt= kaS− k−aAg− kiAgBmg, (6.29b)dBmldt= kbBg− k−bBml− kiAlBml, (6.29c)dBmgdt= kbBg− k−bBmg− kiAgBmg, (6.29d)dBgdt=kaS¯L−kbBgL. (6.29e)The Balanced Inactivation model is not conservative, and so this system is used forthe bifurcation analysis in Figure 6.19.As expected from the LEGI design of this model, the bifurcation diagrams forthe LP system are globally stable; there are no Turing or threshold regimes presentthat would lead to patterning. Instead, a polarization response results purely whena stimulus is overlaid on the model, and disappears if the stimulus is removed. In10620 40 60 800510Timex20 40 60 800510Time20 40 60 800510Timea) b) c) 20 40 60 800510Timexd) 28Figure 6.20: Kymographs of the modified Balanced Inactivation model inEquation 6.27. Simulations were conducted with different types ofstimuli: (a) standard stimulus in Equation 6.3 with ε = 10 (b) originalstimulus, becomes uniform (S¯ = 10) at t = 40 (c) original stimulus,direction reversed at t = 40 and (d) original stimulus, occurring at eachdomain edge.this case, LPA confirms that the cell has no regimes in which spontaneous polar-ization or maintenance of the response in the absence of a stimulus would occur.It provides no information about the parameters that cannot be learned from bifur-cation diagrams of the reaction kinetics alone. The steady state value predicted byLPA is for the situation where a uniform stimulus is present. A, Bm and B would behomogeneous across the domain, and decay to zero if the stimulus is removed.6.5.3 PDE simulations and connection to LPAIn this section, I will examine PDE simulations of the Balanced Inactivation model.However as there are no patterning regions in the LP bifurcation diagrams, few con-nections can be made to the information gained from LPA. Instead, I will direct myfocus to the effects on the polarization response under different parameter valuesor types of stimuli.In Figure 6.20, we see the response of the model to different types of stimuli.Moving from (a) through (d), first the general stimulus is applied. This stimulusis of the form given in Equation 6.3 with a value of ε = 10 (corresponding to thevalue of S¯). As this stimulus only exists as a tall, narrow perturbation to the HSSwithin the initial conditions, and effectively disappears immediately, no patterning1070510x0 5 100246810xConcentration of A  p = 0.05p = 0.15p = 0.250510xa) b) 20 40 60 800510Timex28Figure 6.21: Balanced Inactivation simulations for different parameters, us-ing the original stimulus and parameter values from Table 6.5. (a) Lineplots for the gradient magnitude p = 0.05, 0.15, and 0.25 at t = 80 sec.We notice that as an increase in p leads to an increase in the steepnessof A, primarily at the front of the cell. (b) Kymographs for ki = 1, 5,and 15 (top to bottom respectively). Increasing ki appears to decreasethe intensity of the switch-like response.results. In view of this result we will use the persistent, gradient-type stimulus forthe remainder of the simulations.Simulations in parts (b) and (c) begin with the original stimulus, where param-eter values are given in Table 6.5. At t = 40 sec, the type of stimulus switches toeither a uniform value of S¯ or reverses the direction of the gradient. The modelresponds immediately by losing polarization or reversing the direction of polariza-tion. The final kymograph, part (d), is the result of a bi-directional stimulus. Aunique axis of polarization is not maintained and the cell is able to develop twopolarization fronts. From Figure 6.20, it is evident that the model mimics the over-laid stimulus through A, with the modification that A becomes more amplified andswitch-like than S(x).Figure 6.21 shows the result of modifying either (a) the gradient steepness108parameter p or (b) the inactivation parameter ki. The line plots in (a) show solutionsof A at t = 80 sec for the values p = 0.05, 0.15, and 0.25. Increasing p correspondsto an increase in the external stimulus gradient, which in turn causes a greaterresponse within the cell. The front of the cell responds with a greater concentrationof A, while the back changes very little, leading to a response which is qualitativelymore switch-like.The kymographs in Figure 6.21(b) correspond to the values ki = 1, 5, and 15(default value is ki = 10). The resulting pattern changes very little: the switch-likeresponse becomes more pronounced as ki increases but the values of A at x = 0 andx = 10 are consistent across all diagrams. From the LP (and reaction kinetic) bifur-cation diagram, it is expected that the HSS value Ass will change. This predictionis consistent with simulations, with Ass = Bmss respectively taking on the values3.06, 1.39, and 0.81 for each value of ki (under a uniform stimulus). The HSS of Bremains consistent at 3.33 for all simulations.6.6 SummaryThrough this chapter, I analyzed five models using LPA, Turing analysis (for the Ot-suji and Goryachev models), and PDE simulations. LPA bifurcation diagrams wereshown for representative parameters, and they revealed new information about theinitiation of patterns in many of the models. In particular, many models havethreshold regimes (regions in which the model requires a perturbation of a cer-tain height for polarization to be initiated). Turing analysis proved to be importantin defining regions of stability for the Otsuji and Goryachev models. These regionswere not detected by LPA, because it does not account for specific diffusion rates,only the limiting cases of zero or infinite diffusion. PDE simulations confirmed theresults predicted by both LPA and Turing analysis, and also revealed long-term be-haviour. Simulations are important here, since the predictions by LPA and Turinganalysis are only valid for the pattern’s initiation. With these results, we have newinformation regarding the patterning regimes in each model.109Chapter 7Results and comparisonsIn this chapter, I will compile and interpret the information obtained for the fivemodels in Chapter 6, as well as the Wave Pinning model in Chapter 5. Prior to us-ing local perturbation analysis (LPA), comparing these six cell polarization modelswas difficult. Save for computing many partial differential equation (PDE) simula-tions under different combinations of parameters, we did not know exactly whichregimes would be present in the models or how robust those regimes were. Nowwe have a broader view of the parameter space in these six models. LPA creates aclearer picture of the regions in which pattern formation is possible.7.1 Summary of resultsThis section will summarize the results for each of the models. Table 7.1 lists thepatterning regimes and other features (Hopf bifurcations or limit cycles) that werepresent in the LPA bifurcation diagrams.7.1.1 Wave Pinning ModelThe Wave Pinning model was created to simulate polarization of the active andinactive forms of a single Rho GTPase [49]. In the original publication, the au-thors focused on the formation of a polar distribution through wave fronts whichstall. Referring to Table 7.1, we see that each different patterning regime is repre-sented. Mori et al. [49] use several different transient stimuli to test the response of110Model Patterning Regimes Other FeaturesWave Pinning I, II(a,b), III noneGoryachev I , III noneOtsuji I, II(a), III noneMeinhardt ASDM I, II(a,b), III Hopf, Limit cyclePseudopod Centered I, II(b), III HopfBalanced Inactivation I noneTable 7.1: Summary of the results obtained from LPA, Turing analysis, andPDE simulations in Chapter 6. The patterning regimes are as follows: I =stable, II(a) = upward threshold, II(b) = downward threshold, III = Tur-ing. Recall that upward (downward) thresholds are regions that require aperturbation that is larger (smaller) than the HSS.the model. A localized stimulus, graded stimulus, reversing orientation, and largeamplitude noise all result in polar distributions. LPA predicts that, provided thestimulus reaches a height greater than that predicted by the unstable local branch,polarization results. This fact indicates why large amplitude noise is required togenerate a response: it must be large enough to surpass the threshold. My simula-tions show that if the stimulus is only applied to the initial conditions, the heightof the perturbation needs to be of greater magnitude than that predicted by LPA.7.1.2 Goryachev ModelBudding in yeast cells inspired the development of the model by Goryachev andPokhilko [15]. The authors create a detailed biochemical model which featuresthe GTPase Cdc42 along with its regulatory molecules and downstream effectorBem1. Emerging from this model is a Turing-type instability which serves to createa peak of activation. The LPA bifurcation analysis does not reveal any additionalregimes, but a Turing analysis shows that the model exhibits stability at extremeparameter values. In addition to the instability, the authors hypothesize that theinterconversion between cytosolic and membrane forms of Cdc42 result in a unique111bud. Simulations for different parameter values consistently resulted in a tall butnarrow peak of activation, and the Goryachev model only ever formed one peak.The sole effect different parameters seemed to impart was the length of time beforea stable, polarized distribution was reached.7.1.3 Otsuji ModelOtsuji et al. [54] developed their conceptual model following connections theymade between more detailed models. The authors tested the idea that a mass con-servative system with diffusion driven instability was sufficient to generate a single,stable polarization front. LPA bifurcations and PDE simulations in Section 6.2 in-dicate the presence of both a Turing regime and an upward threshold regime. ATuring instability was expected [54], but the presence of a threshold regime is new.In comparison to the Turing regime which is robust, the width of the thresholdregime is small, and the predicted perturbation height grows exponentially as ourbifurcation parameter decreases. In addition, the threshold perturbation height ismuch greater than predicted by LPA. When considering biological relevance, thisregion may be more representative of stability, and could indicate a region in pa-rameter space where the model becomes increasingly stable to stimulation.We know from a Turing analysis that different parameters may result in mul-tiple peaks forming. This fact is confirmed by simulations, where we have foundregimes with one, two, or even three peaks forming. Multiple peaks always resolveinto one polarization front in the Otsuji model, but the lifetime of these transientpeaks changes under different parameter values. At this point, PDE simulations arethe only method to inform us of the lifetime of these peaks, since any predictionsmade by a Turing analysis break down once the peaks grow. Secondary peaks withlong lifetimes represent multiple polarization fronts in a cell, and may restrict thesize of the Turing regime.7.1.4 Meinhardt Activator-Substrate Depletion ModelMeinhardt’s ASDM is based on a mechanism that admits a Turing instability. Thefull system models the localization of FtsZ to the center of a 1D bacterial cellduring cell division, but I focus on the pair of equations responsible for polariza-112tion. The activator refers to membrane-bound FtsZ and it depletes the substrate,cytosolic FtsZ, as it localizes to the membrane. The previous models are all massconservative, but due to unequal production and decay rates, Meinhardt’s ASDM isnot.The bifurcation diagrams in Figure 6.10 reveal new threshold regimes in ad-dition to the default Turing regime. An interesting feature was the presence of aHopf bifurcation, causing an unstable limit cycle and oscillations in the pulse it-self. Exploring the connection between the Hopf bifurcations in LPA and possibleoscillations in the full PDES was beyond the scope of this thesis, but remains aninteresting open problem for further investigation. This model is the only examplein my thesis where a regime predicted by the LP system could not be actualizedthrough PDE simulations. In Figure 6.11, despite a very large downward perturba-tion (both in width and magnitude), the model did not maintain polarize.Simulations reveal that multiple peaks may form in the Turing regime. Un-like the Otsuji model, when these peaks form they are stable. In the context ofMeinhardt’s full reaction-diffusion equation (RDE) system for bacterial cell divi-sion, the additional equations act to resolve any multiple peaks into one centralpeak. However, the two reaction equations that I consider in this thesis are concep-tual, and could potentially be applied to other systems. In fact, equations for theWave Pinning and ASDM are very similar (see Equation 5.1 for Wave Pinning andEquation 6.18 for the ASDM), prompting an anonymous reviewer to suggest the twomodels were exactly the same. Lack of conservation in Meinhardt’s ASDM is likelythe cause for multiple stable peaks, and could be problematic if these equations areapplied to other cell systems that require a single axis of polarity.7.1.5 Pseudopod Centered ModelThe Pseudopod Centered model represents an actual case in which a set of equa-tions were modified for use in different cell systems. An ordinary differential equa-tion (ODE) system was developed by Meinhardt [45] to represent polarization foreither chemotactic cells or growth cones of neurons. It was subsequently modi-fied into RDES representing pseudopod generation and coupled to a 2D deformingboundary in [51, 52].113The original parameter set corresponds to a Turing regime, but LPA bifurcationdiagrams reveal a new downward threshold regime and regions of stability. PDEsimulations for the threshold regime, visualized in Figure 6.15, show one stablepeak whose width and location are purely determined by the perturbation appliedto the model. Here the pattern is not dynamic, but in the Turing regime peaks arecontinually growing and shrinking. When a stimulus is constantly present, forma-tion of peaks is localized in the direction of the perturbation and a stronger stimuluscorresponds to increased localization of the peaks. Diffusion rates were not clearlydefined for this model, so several values were tested. When I increased the dis-parity between diffusion rates in simulations, the accuracy for predictions (e.g.,regime boundaries and threshold values) made by the LPA bifurcations increased.7.1.6 Balanced Inactivation ModelThe final model I analyzed was the Balanced Inactivation model by Levine et al.[33]. This model is based on the local excitation / global inhibition (LEGI) scheme,and does not feature inherent pattern forming mechanisms. Polarization is only ob-served for a graded stimulus; a uniform stimulus results in a rapid uniform increasein components, but the response soon decays [33].In a model such as this one, what do the LPA bifurcations look like? Whereaseach of the previous models exhibited inherent pattern formation as a result ofnonlinear reaction terms and diffusion, the LP bifurcation diagrams for BalancedInactivation in Figure 6.19 were purely stable. Pattern formation in the BalancedInactivation model will not result from any sort of transient perturbation. But whena graded stimulus is present, the components are organized in a polarized patternbased on the gradient, a result not predicted by LPA or linear stability analysis.Since there were no distinct regimes to test, I used PDE simulations to investigatethe patterns that form in response to different stimuli. The Balanced Inactivationmodel reacted with a modified readout of the stimulus, where the response wasswitch-like rather than purely an amplification of the gradient.1147.2 Comparing models according to patterning regimesIn this section, I will compare the models according to the predominant regimes(Turing, threshold, and stable) and the robustness of those regimes. For many ofthe models, the default parameter sets resulted in polarization due to Turing insta-bilities. LPA revealed that several models also contain previously unknown thresh-old regimes. While a region of stability does not result in patterns, it is a regioncommon to most of the models, and is particularly interesting in the Balanced Inac-tivation model. The robustness of a regime is indicated by the amount of variationa parameter may undergo while still remaining in that particular pattern regime. Avery robust regime will constitute a large portion of parameter space.7.2.1 Turing instabilitiesIn a Turing regime, the model is unstable to small amplitude perturbations. Themost robust regimes were found in the Goryachev and Otsuji models. A shortfallof LPA is that it does not take different diffusion rates into account. In the caseof the Goryachev and Otsuji models, a Turing analysis revealed a limit to the Tur-ing regime, not predicted by LPA, when realistic diffusion rates are considered.Meinhardt’s ASDM and Pseudopod Centered models had comparably sized Tur-ing regimes, and the Wave Pinning model had the least robust regime. Among allthese models, Wave Pinning was the only model whose default parameters werenot located in the Turing regime.The Wave Pinning, Goryachev, ASDM, and Otsuji models are all comparable intheir ability to form one stable peak of activation in response to a very small per-turbation. Under different parameter sets, multiple peaks of activation are presentin the Otsuji, ASDM, and Pseudopod Centered models, but each model exhibitsdifferent long-term behaviour with respect to the additional peaks. In the conser-vative Otsuji model, multiple peaks will always resolve into one peak, and theirtransient lifetime is dependent on the parameter values. The ASDM is not a conser-vative model. In the case where two peaks form, both are stable. In certain regionsof the Turing regime, the Pseudopod Centered model is unique in its behaviour.Even with the default stimulus from Equation 6.3, the polarization pattern does notfreeze. Instead, many peaks form and continue to rise and fall akin to pseudopods115extending and retracting.What does it mean when Turing analysis predicts that two peaks can form?In some cases, it is not biological for two stable or long-lived peaks of activationto exist. This fact is particularly true when considering cells that do not movevia pseudopods. But do we dismiss a model biologically just because a Turinganalysis predicts multiple peaks? Recall that Turing analysis is only valid for theinitial formation of the pattern, just like LPA. Once the perturbations start to grow,the assumptions made for the Turing linear stability analysis break down. Turinganalysis cannot then predict the long-term state of the solution, but only the initialrise of patterns.Consider the Otsuji model. Turing analysis reveals that multiple modes areunstable. But, through the combined effects of autocatalysis and conservation,multiple peaks eventually diminish while one peak grows to its maximum size.Biologically, a cell may receive two stimuli and form a transient response to both.Since the cell cannot move in two directions at once, one front becomes dominantwhile the peak of activation at the other front recedes. It may not be necessaryto deem a model as “unbiological” if a Turing analysis reveals multiple unstablemodes.7.2.2 Threshold regimes (upward and downward)In a threshold regime, the model is destabilized by perturbations of a certain size.The perturbation may be larger or smaller than the HSS, and is consequently termedupward or downward respectively. Several models have threshold signatures intheir respective LPA bifurcation diagrams: Wave Pinning (upward and downward),Otsuji (upward), ASDM (upward and downward), and Pseudopod centered (down-ward). The Wave Pinning model is the only system in which this region was pre-dominantly studied by its authors. The presence of a threshold regime in the othermodels was either not described by the authors, or is newly revealed by LPA.I found three different signatures for threshold regimes in the LPA diagrams.The first appears as a stable global branch, followed by an unstable local branch,and then another stable local branch. A perturbation past the unstable local branch,116whether this is above or below the global branch, will settle at the stable localbranch, effectively creating a polarized distribution. This type of threshold regimeis found in the Wave Pinning and Pseudopod Centered models. The other type ofthreshold regime was found in the Otsuji model and Meinhardt’s ASDM. Here, astable global branch had an unstable local branch either above it (Otsuji) or below(ASDM). Due to conservation in the Otsuji model, the pulse settled at a finite value.In Meinhardt’s ASDM, the pulse was pushed to zero. The third type of thresholdregime was unique to the ASDM. Hopf bifurcations were present on both the localand global branches. The global branch was stable, but an unstable limit cycle waspresent on the local branch. A perturbation toward it resulted in oscillations of thepulse itself before reaching a stable configuration.The threshold values predicted by LPA matched some models more closely thanothers. The Otsuji and Wave Pinning models required much larger perturbationsthan predicted by the LPA bifurcations. The Pseudopod centered model was closerto the prediction, and the ASDM required a far smaller perturbation than the un-stable local branch. In terms of the downward threshold regime, I often neededto increase the width of the perturbation from 5% of the domain, because increas-ing the pulse magnitude would result in negative values. In fact, in the case ofone region in the ASDM (see Figure 6.11), I was unsuccessful in obtaining a stablepolarized distribution despite a large width and magnitude in the perturbation.7.2.3 Stable regionsStable regimes indicate a region in parameter space in which the model is impervi-ous to polarization. Biologically, researchers are able to manipulate certain aspectsof the cellular biochemistry which could render the cell unresponsive to stimuli. Itis interesting to consider regions in parameter space in which the models are stable.With the LPA bifurcations, I found large regions of stability in the Wave Pinningand Pseudopod Centered models (Figure 5.2 and Figure 6.14, respectively). Theseindicate that increasing or decreasing a parameter by too much destroys polariza-tion. A Turing analysis revealed regions of stability in the Otsuji and Goryachevmodels (see Section 6.2.2 and Section 6.1.2) that were not predicted by their re-spective LP systems. The transition into a stable regime can be abrupt, as in the117Wave Pinning model, or gradual, as in the Pseudopod Centered model, where thedynamic oscillations become smaller and smaller in magnitude.The Balanced Inactivation model has no inherent polarizing ability; its LPA bi-furcation diagrams are completely stable. However, it still responds when given agraded stimulus. Provided the stimulus is persistent, the cell responds by ampli-fying and modifying the shape of the gradient to reach a polar distribution. Notethat other models like Wave Pinning or Pseudopod Centered do not respond in thesame way in their respective stable regimes. That is, the stimulus is not modified.The models feature different underlying dynamics, and so the stable regimes takeon different meanings.7.3 Comparing models according to qualitative featuresRecall that in Section 3.4, I listed a number of qualitative features of polarizingcells [29], and addressed a few questions concerning cell polarization models.Since we lack biological details, it is helpful to interpret the models accordingto qualitative rather than quantitative details. We do not expect the simulationsto match quantitative biological data. In this section I will compare the modelsusing these qualitative features and address the cell polarization questions. (Notethat the relative robustness of patterns between the models was discussed above inSection 7.2.)Maintenance of polarity, amplification, sensitivity, and adaptation are all fea-tures pertaining to the long-term patterns. The results from Chapter 6 did not revealnew information concerning these features, so we may refer to Table 3.2. With re-spect to maintenance of polarity, peaks of activation in the Pseudopod Centeredmodel will be present in the absence of a stimulus, representing a non-uniform do-main, but there is no specified direction. Whether a model has a unique axis ofpolarity or not can be determined through long-term simulations or Turing analy-sis. Multiple peaks of activation through Turing instabilities were discussed abovein Section 7.2.1. Under multiple sources of stimulation, the Wave Pinning andBalanced Inactivation models will also admit multiple activation peaks.An interesting issue in cell polarization modelling is determining which mod-els feature spontaneous polarization and which models require a stimulus of a cer-118tain size. Models with threshold-type responses can account for resting states incells, in which a cell remains unpolarized until it receives a signal. Turing regimesaccount for spontaneous polarization, or represent cells that are sensitized and po-larize in response to minute perturbations in the environment. Under appropri-ate parameter sets, each model could polarize spontaneously except the BalancedInactivation model. Balanced Inactivation is interesting in that it does not sponta-neously polarize or require a certain threshold, but it is the only model in this thesiswhich accounts for adaption to a uniform increase in chemoattractant. So, while itdoes not necessarily require a threshold-sized value for the stimulus, it does requirea stimulus gradient.Another possibility for achieving a resting state is modification of the reactionrates. Since the networks are very complicated, the reaction rates may be dependenton other signalling proteins in the network. For example, GEFs and GAPs modifythe reaction rates for conversion between the active and inactive forms of GTPases(see the discussion in Section 2.1). Under different circumstances, shifts in theseother proteins can cause shifts in the reaction rates, corresponding to a change inthe pattern regime. So, if a model features a Turing regime with a neighbouringthreshold or stable regime, it could potentially account for resting states. Under thisconsideration, the Otsuji and Pseudopod Centered models would also be included.The Goryachev model features stability, but at such extreme parameter values thatthis regime is not likely relevant.Some motile cells move via pseudopods, while the motion in others is drivenby a lamellipod. Both occur through actin dynamics, but superficially these twostructures appear very different. Pseudopods are constantly being extended andretracted; the lamellipod is a more robust and long-lived structure and the cellmoves in a gliding motion. But is there a distinction in the underlying mechanismsfor these two processes? The unstable nature of pseudopods indicates a long-termpattern that is unstable, and possibly an additional component to aid in destabilizingthe polarized pattern. We can compare the Pseudopod Centered model with modelslike Wave Pinning and Otsuji. Comparing the LPA bifurcation diagrams of thesethree models does not indicate a distinct difference in the underlying dynamics.It appears that the dynamics for pattern initiation are the same, but the minimalmodel must be extended in some way to account for pseudopod dynamics. A119third equation in the Pseudopod Centered model, the “local inhibitor,” serves todestabilize the pattern.7.4 Wave Pinning versus Meinhardt’sActivator-Substrate Depletion ModelAt first sight, it appears that Meinhardt’s ASDM is very similar to the Wave Pin-ning model. The reaction kinetics in both models (Wave Pinning Equation 5.1and Meinhardt’s ASDM Equation 6.18) feature autoactivation and comparable de-cay terms. However, Meinhardt’s model lacks conservation. Simulations underthe default parameter sets are also qualitatively similar (see Figure 5.1 and Fig-ure 6.9). One might conclude that the comparable reaction terms and long-termPDE simulations are evidence enough that the models can be used interchangeably.LP analysis and PDE simulations help to reveal details of the underlying dynamics,and aid in the comparison of these two models. I confirm that under the defaultparameter values, the mechanisms for pattern formation are very different. Thereare more similarities than previously thought when considering a broader range ofparameters, but quantitative differences remain.We know before performing LPA that the default parameters of Meinhardt’sASDM lie in a Turing regime, while polarization in the Wave Pinning model isdriven by the propagation of a wave front. Looking at the corresponding regionsfor the default parameter sets, LPA bifurcation diagrams clearly show that bothmodels require very different types of stimulation to generate a response. Mein-hardt’s ASDM is sensitive to the smallest amount of noise, while the HSS of theWave Pinning model is actually linearly stable. A perturbation to the Wave Pin-ning model must cross a threshold given by the unstable local branch before we areable to achieve polarization.We will now extend our view outside of the default parameter values. Interest-ingly, both models contain regions of Turing instability and threshold behaviour, anew patterning regime for Meinhardt’s ASDM. Both upward and downward thresh-old regimes occur, similarly to the wave pinning model. Conservation is a keydifference between the two models, and Mori et al. [49] note that it is not possiblefor Meinhardt’s ASDM to exhibit wave pinning behaviour due to lack of conserva-120tion. Therefore, the threshold behaviour is qualitatively similar, but not the same asthe Wave Pinning model. Lack of conservation in Meinhardt’s ASDM also allowedfor multiple stable peaks, particularly in the Turing regime. Multiple peaks in theWave Pinning model are not stable, and eventually resolve into one peak.7.5 Connecting patterning regimes to simulation stimuliCells can be stimulated in numerous ways: persistent graded stimulus, transient lo-calized stimulus, random conditions, and multiple stimuli (either competing witheach other or reversing direction). Jilkine and Edelstein-Keshet [29] represent thesedifferent types of stimuli in PDE simulations of the models. With the LPA bifurca-tion diagrams, we are clearly informed of the type of initiating stimulus requiredunder different parameter values. Before diving into PDE simulations, we can ana-lyze the LPA map to understand what types of stimuli will incite polarization. ForSection 7.5.1 and Section 7.5.2 I pose a general question about Turing and thresh-old regimes: how does the type of regime inform us of what types of stimuli themodel will polarize from? With information from the LPA bifurcations, we do nothave to ask if polarization will result from a particular stimulus. The focus can bemore directed to what the polarization response looks like.7.5.1 Turing regimesIn a Turing regime, it doesn’t matter what form the stimulus takes. Our model isalways going to produce patterns. Localized or graded, persistent or transient, andlarge perturbations or small amounts of noise are all sufficient features for the stim-ulus to excite patterning. In terms of PDE simulations, even machine noise initiatespatterns. Biologically, a Turing regime corresponds to cells that are sensitized orto cells that polarize spontaneously. A typical trend that I noticed in simulationsfor this thesis was that a larger perturbation (stimulus) caused patterns to developfaster.As we have seen through the simulations in Chapter 6, the form of the stimulusmay have an effect on the type of long-term pattern that results. In the pseudopodcentered model, the original stimulus localizes the peaks of activation in the center(and changing the direction of the stimulus reorients the peaks). It is continually121present to enforce the location of the activation peaks, removal of the stimulus re-sults in activation peaks located throughout the domain. In contrast, models likeGoryachev and Otsuji which also operate in a Turing regime do not require a con-stant stimulus to maintain the orientation of polarization. Changing the directionof stimulus does not reorient polarization, once formed the pattern is static.7.5.2 Threshold regimesIn a threshold regime, the size of the perturbation does matter. Patterns are initiatedunder many different forms of stimuli whether very narrow or graded, provided thestimulus is large enough in magnitude to pass a certain threshold. Graded stimuliwith a too weak gradient, or noise whose amplitude is too small will not resultin patterns, but increasing the magnitude of these perturbations past the thresholdwill likely lead to patterning. LPA informs us of the height necessary for patterninitiation, and in some cases (such as Wave Pinning) even informs us of the highersteady state the perturbation will settle to.7.6 Comparing simple and detailed modelsBiological signalling networks are a complicated tangle of many interacting com-ponents, as evidenced by the already simplified discussion on an actin polymer-ization pathway in Section 2.1. Oftentimes we lack experimental detail for thesecomplicated networks and cannot parse certain types of feedback, so we createsimple models that are stripped of cumbersome details. These simpler models aremuch easier to analyze, and are helpful for elucidating some of the underlyingmechanisms at play in the networks. However, we want to ensure that these ab-stract models preserve the essential dynamics of the detailed models. With learningabout dynamic behaviour of the actual signalling networks as one goal of compu-tational modelling, we eventually want to return to complex models.Analysis can be much more difficult to approach for complex models, part ofthe reason why simple models are so appealing. However, it is relatively straight-forward to use the LPA method for investigating even complicated models. In [23],Holmes et al. create a model based on experimental work described by [35] and usethe LP method of analysis. Interactions between the GTPases Cdc42, Rac, and Rho122and the three forms of phosphoinositides make up a complex model of nine PDES.The resulting LP system is composed of fifteen ODES. But despite the high dimen-sionality of the system, analysis via bifurcation software is easier than analysis ofmany nonlinear PDE systems. A bifurcation diagram for the basal activation rateof Rac reveals a very similar signature to the Wave Pinning model. In their system,the additional biological detail does not change the core motif for pattern initiation.The authors even use LP bifurcations to understand experimental manipulations byobserving changes to the LPA signature of a single key parameter while allowingother parameters to change.123Chapter 8ConclusionsThere were two major goals of this thesis:1. demonstrating the use of LPA, including its benefits and caveats, and2. analyzing and comparing six cell polarization models using information gainedvia LPA and PDE simulations.In Chapter 3, I provided some background on current methods of analysis for PDESand introduced the method of LPA in Chapter 4. The Wave Pinning and ActinWaves model examples introduced in Chapter 5 demonstrated the use of LPA: howto convert the PDES into the LPA system and the subsequent bifurcation analysison these ODES. This section also illustrated how the LPA bifurcation signatures areaffected when a simple model (Wave Pinning) extends into a model that is morecomplex (Actin Waves). Chapter 6 was the main section for analysis. I analyzedfive polarization models with distinct features through LPA, PDE simulations, andTuring analysis (for the Otsuji and Goryachev models). In Chapter 7, I comparedthe six polarization models to each other and made general observations linkingthe LPA results to cell polarization models.In this final chapter, I will discuss some conclusions concerning the LPA method.What were the advantages of the method and how did it help in revealing under-lying features of the cell polarization models? How does LPA aid us in comparingand classifying the models? What information needs to be gained by other meth-ods of analysis? I will also provide some closing remarks on the cell polarization124models and suggest future directions for this work and for LPA.8.1 Conclusions concerning LPAUsing LPA to reveal more information about pattern initiation in the different mod-els was a major aspect of my thesis. Through this application, I have highlightedsome of the benefits of the LPA method, as well as areas in which other methodsare required to reveal more information. Additionally, situations arose in which Ineeded to modify the models before applying LPA. Here I will address some of thegeneral benefits and caveats of using LPA, listed in Chapter 4, in the context of thecell polarization models.8.1.1 LPA quickly reveals information about pattern dynamicsIf a model satisfies the requirements, namely it is a reaction diffusion system inwhich DuDv, LPA is generally easy to implement. The major step in convertingto the LP ODE system is determining which variables are global and which arelocal. Once the LPA system is built, it is easy to use ODE bifurcation software tocreate bifurcation diagrams for different parameters. A major part of analyzingthe models in Chapter 6 was interpreting these bifurcation diagrams, which revealunderlying information about the initiation of patterns. We can use this informationto make predictions on the model’s behaviour.LPA recovers the number, values, and stability of the steady state(s) in a model,as well as dynamics including threshold responses, Turing instabilities, limit cy-cles, and Hopf bifurcations. Over the six different models I analyzed, all of thesedynamics were represented. Since most models originally operated on parametersets in the Turing regime, much of the information gained was new. However, someof the newly detected regions were not robust to parameter changes or were locatedat extreme parameter values, and so may not be considered biologically relevant.In addition to revealing new dynamics of these models, LPA also informs us ofthe robustness of the regions. For each parameter, we can see the variation allowedin that particular parameter value before the regime changes. Some models, likethose by Otsuji and Goryachev, have very large patterning regimes that are robustto changes in parameters. Others, like the Wave Pinning model, exhibit purely125stable behaviour when a parameter is too large or small, and therefore have a finelytuned region in which polarization can occur. We also immediately know the typeof stimulus that is required for a response. With the information on steady states,pattern-forming dynamics, and parameter robustness, LPA provides a more directroute for understanding and comparing different models.8.1.2 Model modifications and assumptionsModifications were required for some models because LPA does not take spatialvariations into account. I approximated the stimuli in the Otsuji and PseudopodCentered models by a constant, average value, and simplified the boundary in-tegral in the Pseudopod model. The Balanced Inactivation model required moreextensive revision, and was rebuilt in order to retain reaction information locatedin the boundary conditions. This new format for the Balanced Inactivation modelretained all information from the original derivation, but it was reorganized suitLPA and to match the format of the other models. When applying LPA to othermodels, we need to be sure that spatially dependent components are dealt withappropriately.In most cases, determining which variables are slow or fast is straightforward.However, some models containing more than two equations do not always leadto distinct slow and fast categories. For example, in early investigations of thePseudopod Centered model, I had the case where Da < Dc < Db [52]. It wasunclear which class Dc belonged to. I suggested creating two LPA systems, wherec was both slow and fast, and PDE simulations would then determine which versionof the LP system was a closer match to dynamics in the RDE. In fact, in references[23, 41], the phosphoinositides were treated as both slow and fast variables, sincetheir diffusion speeds fall between the slow and fast diffusion speeds of the activeand inactive (respectively) GTPases.8.1.3 Other methods provide additional informationLPA is a useful tool for easily and quickly visualizing when a pulse of a certainheight will grow or decay. It is also helpful in recovering Turing instabilities,which reveal when small perturbations grow into a macroscopic pattern. But once126the pulse becomes wider or the perturbations grow taller, these predictions breakdown. LPA only tells us about a pattern’s initiation, that is, whether a pattern willform and what type of perturbation generates a response. To reveal the long-termdynamics of the system, we require PDE simulations which show the actual polar-ization patterns, how long they take to form, and whether the pattern is frozen ordynamic.PDE simulations are also important for defining the boundaries between regimes,determining the exact height of a threshold perturbation, and visualizing transitionsbetween regimes. LPA is an approximation: we are approximating pulse dynamicsin a PDE system by a system of ODES. Therefore, pulse dynamics at the boundariesor the perturbation size for threshold behaviour may not be correctly predicted.These discrepancies may be minor, where LPA closely approximates the bound-aries, or it can be large. The LPA bifurcations predicted threshold behaviour inthe Otsuji model, but the threshold required to produce polarization in simulationswas much larger than predicted. Another interesting phenomenon is the transitionfrom one region to the next. Transitions could be immediate or gradual, and PDEsimulations can reveal interesting transition patterns when one patterning regimeblends into to the next [43].Another caveat of LPA is that it does not recover an entire Turing analysis. LPApredicts when the HSS will be unstable to small amounts of noise, but we gain noinformation on wave numbers, and therefore do not know how many peaks mayform. Additionally, the conditions for Turing instabilities depend on the diffusionrates, and LPA only provides predictions for extreme diffusion speeds (zero or in-finite diffusion). A full Turing analysis will reveal the unstable wave numbers andrestrictions from diffusion.Knowing the number of peaks that are unstable is important in terms of mod-elling cell polarization, because multiple peaks of activation correspond to multipleprotrusion zones. In some cells, like budding yeast cells or motile eukaryotes thatmove via lamellipods, this property could be detrimental to the cell’s reproductionor survival (however, see Section 7.2.1). The presence of multiple peaks or effectsfrom diffusion speeds not predicted by LPA may reduce the relevant parameter set.1278.2 LPA and cell polarization models8.2.1 Answering polarization modelling questionsIn Section 3.4.2, I presented a few key questions for cell polarization modelling. Ianswered them specifically for each model in Section 7.3, but I will address themmore generally now. LPA is immediately helpful for two of these questions: whichcells polarize spontaneously versus requiring a threshold-sized stimulus and howrobust the polarization response is to changes in parameters. If the model polarizesspontaneously, the LPA diagram will contain a Turing regime. We saw several typesof threshold regimes through the six cell polarization models, with the commonfeature that a stimulus is required to pass a certain threshold before polarization willoccur. Robustness of these regimes is clearly indicated in the bifurcation diagrams.In fact, detection of patterning regimes and their robustness are the two key reasonsLPA is so helpful.The other three questions require additional analysis: does a model featuremaintenance of polarity; is there a difference in the underlying mechanisms forpseudopods versus lamellipods; and is there a unique axis of polarity? Since LPApredictions are typically only valid for the start of a simulation, maintenance of po-larity must be confirmed through PDE simulations. If the model remains polarizedfor a biologically-relevant length of time once the stimulus is removed, we can con-clude that the model features this property. Whether a model will exhibit a uniqueaxis of polarity may be determined through a Turing analysis or simulations. Re-call that additional peaks in the Otsuji model decayed, so while a Turing analysispredicts multiple polarized fronts, PDE simulations are still important to confirmthe results. The question of pseudopod versus lamellipods cannot be answered byLPA, as we may recall that the bifurcations diagrams for the Pseudopod-Centeredmodel had the same appearance as for the other models. It requires simulationsand some insight into the particular cell polarization model.8.2.2 Classifying modelsIn general, we can classify the models according to their LPA signatures. For thisthesis I was predominantly concerned with Turing and threshold regimes, as well128as regions of stability. LPA quickly reveals the presence of these regimes, and wecan further compare the models by each regime’s robustness. Remember that therobustness may differ according to the bifurcation parameter. From here, additionalanalysis provides another level of classification. For example, PDE simulationsreveal the long-term behaviour of the model and its response to different stimuli,and Turing analysis tells more about the unstable wavenumbers or limitations whenactual diffusion rates are taken into account.We may consider the patterning regimes to be the first level of classification forour polarization models. The following list describes some secondary features forclassifying or comparing models.• Turing regime: In this regime, the model is unstable to even the smallestperturbations, so the model indicates a situation in which polarization occursspontaneously, or where the cell is sensitized for very slight stimuli. Multiplewavenumbers can become unstable, and so we may encounter multiple peaksof activation. Parameter values may modify features such as the time forpolarization and the width of the polarization peak.• Threshold regime: The corresponding cell requires a stimulus which sur-passes a certain threshold before polarization will occur, the height of whichis predicted by LPA. Further investigation into the reaction terms may revealthe cause for this response, such as bistability.• Stable regime: A transient stimulus will not yield polarization. It is possiblefor a model with a LEGI design, to exhibit polarization features to a constantstimulus, even if it is inherently stable.• Additional features: LPA may reveal additional features contained in themodel such as Hopf bifurcations or limit cycles. These features can result ininteresting dynamics of the model.• PDE simulations: Simulations can indicate the time until polarization occursand the size of the active zone, whether the model is dynamic or the patternfreezes, and exactly what the polarization pattern looks like. Responses todifferent stimuli, such as alternating direction, are also helpful for classifyingthe models.129From the mathematical analysis (LPA, Turing, PDE simulations and otherwise),we can add a further level of classification according to the qualitative polarizationfeatures the models exhibit: amplification, sensitivity, maintenance, spontaneouspolarization, pseudopod-splitting, single polarity axis, and adaptation (see Sec-tion 3.4.1). Parameter robustness and the time it takes for a model to polarize areimportant biological connections as well. This list is not exhaustive, and dependingon the cell type or research goals, other features become important. Often, mathe-matical analysis will indicate which features a model exhibits, and we can comparethe results to the actual biological system.8.3 Future directionsWith the additional tool of LPA for analyzing cell polarization models, we have abroader view of the underlying features behind these models. In this thesis, I inves-tigated a variety of models with distinct differences to highlight how the underlyingfeatures compare, and the results inspire directions for future work. Instead of an-alyzing distinct models, we may compare models that are more similar, either infunctionality or in the type of the cell they simulate. Like Otsuji et al. [54], inves-tigating the similarities in models can point to common underlying features, andcould even reveal plausible pattern-forming modules for the biological systems.I stress that the LPA bifurcation diagrams examined in this thesis were one di-mensional. That is, they only offer insight into the model’s behaviour over a rangeof values for a single parameter. In order to more fully investigate the parameterspace and pattern robustness, we can perform a multi-dimensional type of analy-sis. We can automatically generate and analyze many LPA bifurcation diagramsfor a single parameter, but each diagram would have a slightly modified parameterset. This process would provide a much clearer picture of how patterning regimeschange as more than one parameter is varied.In the short time since its conception, LPA has been applied primarily to an-alyzing cell polarization models. But, other than satisfying the requirements forits use, there are no restrictions on the types of models LPA can be used for. Thisfact begs the question: how might LPA and its results be applied to other systems?RDES function in modelling microscopic species such as the active and inactive130proteins in this thesis, viruses and whole cells, to the macroscopic organisms inpredator-prey systems or pest control. When any of these systems feature a dispar-ity in diffusion speeds, LPA would be a useful new tool in investigating patterningregimes.131Bibliography[1] N. Andrew and R. H. Insall. Chemotaxis in shallow gradients is mediatedindependently of PtdIns 3-kinase by biased choices between randomprotrusions. Nature cell biology, 9(2):193–200, 2007. → pages 11[2] S. Bodin and M. D. Welch. Plasma membrane organization is essential forbalancing competing pseudopod- and uropod-promoting signals duringneutrophil polarization and migration. Molecular biology of the cell, 16(12):5773–5783, 2005. → pages 9[3] H. R. Bourne, D. A. Sanders, and F. McCormick. The GTPase superfamily:a conserved switch for diverse cell functions. Nature, 348(6297):125–132,1990. → pages 7[4] H. R. Bourne, D. A. Sanders, and F. McCormick. The gtpase superfamily:conserved structure and molecular mechanism. Nature, 349(6305):117–127,1991. → pages 7[5] T. Bretschneider, K. Anderson, M. Ecke, A. Mu¨ller-Taubenberger,B. Schroth-Diez, H. C. Ishikawa-Ankerhold, and G. Gerisch. Thethree-dimensional dynamics of actin waves, a model of cytoskeletalself-organization. Biophysical journal, 96(7):2888–2900, 2009. → pages 20[6] B. A. Camley, Y. Zhao, B. Li, H. Levine, and W.-J. Rappel. Periodicmigration in a physical model of cells on micropatterns. Physical reviewletters, 111(15):158102–158107, 2013. → pages 62, 63, 64[7] P. Devreotes and C. Janetopoulos. Eukaryotic chemotaxis: distinctionsbetween directional sensing and polarization. Journal of biologicalchemistry, 278(23):20445–20448, 2003. → pages 2, 3[8] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov. MATCONT: a MATLABpackage for numerical bifurcation analysis of ODEs. ACM Transactions onMathematical Software, 29(2):141–164, 2003. → pages 41, 66132[9] L. Edelstein-Keshet, W. R. Holmes, M. Zajac, and M. Dutot. From simple todetailed models for cell polarization. Philosophical Transactions of theRoyal Society B: Biological Sciences, 368(1629):doi:10.1098/rstb.2013.0003, 2013. → pages 3, 20, 38, 62[10] S. I. Ellenbroek and J. G. Collard. Rho GTPases: functions and associationwith cancer. Clinical & experimental metastasis, 24(8):657–672, 2007. →pages 8[11] S. I. Fraley, Y. Feng, A. Giri, G. D. Longmore, and D. Wirtz. Dimensionaland temporal controls of three-dimensional cell migration by zyxin andbinding partners. Nature Communications, 3(719):doi:10.1038/ncomms1711, 2012. → pages 62[12] S. Funamoto, R. Meili, S. Lee, L. Parry, and R. A. Firtel. Spatial andtemporal regulation of 3-phosphoinositides by PI 3-kinase and PTENmediates chemotaxis. Cell, 109(5):611–623, 2002. → pages 6[13] G. Gerisch, T. Bretschneider, A. Mu¨ller-Taubenberger, E. Simmeth,M. Ecke, S. Diez, and K. Anderson. Mobile actin clusters and travelingwaves in cells recovering from actin depolymerization. Biophysical journal,87(5):3493–3503, 2004. → pages 20[14] A. Gierer and H. Meinhardt. A theory of biological pattern formation.Kybernetik, 12(1):30–39, 1972. → pages 19[15] A. Goryachev and A. Pokhilko. Dynamics of Cdc42 network embodies aTuring-type mechanism of yeast cell polarity. FEBS letters, 582(10):1437–44, 2008. → pages 24, 67, 68, 111[16] V. Grieneisen. Dynamics of auxin patterning in plant morphogenesis. PhDThesis Utrecht University, 2009. → pages 35[17] A. Hall. Small GTP-binding proteins and the regulation of the actincytoskeleton. Annual review of cell biology, 10(1):31–54, 1994. → pages 5[18] J. M. Haugh, F. Codazzi, M. Teruel, and T. Meyer. Spatial sensing infibroblasts mediated by 3 phosphoinositides. The Journal of cell biology,151(6):1269–1280, 2000. → pages 6[19] I. Hecht, M. L. Skoge, P. G. Charest, E. Ben-Jacob, R. A. Firtel, W. F.Loomis, H. Levine, and W. J. Rappel. Activated membrane patches guidechemotactic cell motility. PLoS computational biology, 7(6):e1002044,2011. → pages 10133[20] W. Holmes, A. Carlsson, and L. Edelstein-Keshet. Regimes of wave typepatterning driven by refractory actin feedback: transition from staticpolarization to dynamic wave behaviour. Physical biology, 9(4):046005,2012. → pages 41, 49, 52, 53[21] W. Holmes, L. Edelstein-Keshet, S. Maree, and V. Gieneisen. The localperturbation analysis: a method for detecting spatio-temporal patternforming capabilities in fast-slow reaction diffusion systems in preparation.2012. → pages 35[22] W. R. Holmes. An efficient, nonlinear stability analysis for detecting patternformation in reaction diffusion systems. Bulletin of mathematical biology,76(1):157–183, 2014. → pages 35[23] W. R. Holmes, B. Lin, A. Levchenko, and L. Edelstein-Keshet. Modellingcell polarization driven by synthetic spatially graded rac activation. PLoScomputational biology, 8(6):e1002366. doi:10.1371/journal.pcbi.1002366,2012. → pages 122, 126[24] P. A. Iglesias and P. N. Devreotes. Navigating through models ofchemotaxis. Current opinion in cell biology, 20(1):35–40, 2008. → pages 3,9, 21[25] P. A. Iglesias and A. Levchenko. Modeling the cell’s guidance system.Science Signaling, 2002(148):DOI: 10.1126/stke.2002.148.re12, 2002. →pages 22[26] R. H. Insall. Understanding eukaryotic chemotaxis: a pseudopod-centredview. Nature Reviews Molecular Cell Biology, 11(6):453–458, 2010. →pages 11[27] J. E. Irazoqui, A. S. Gladfelter, and D. J. Lew. Scaffold-mediated symmetrybreaking by Cdc42p. Nature cell biology, 5(12):1062–1070, 2003. → pages13[28] C. Janetopoulos, L. Ma, P. N. Devreotes, and P. A. Iglesias.Chemoattractant-induced phosphatidylinositol 3, 4, 5-trisphosphateaccumulation is spatially amplified and adapts, independent of the actincytoskeleton. Proceedings of the National Academy of Sciences, 101(24):8951–8956, 2004. → pages 27[29] A. Jilkine and L. Edelstein-Keshet. A comparison of mathematical modelsfor polarization of single eukaryotic cells in response to guided cues. PLoS134computational biology, 7(4):e1001121. doi:10.1371/journal.pcbi.1001121,2011. → pages 3, 9, 18, 22, 28, 66, 68, 77, 118, 121[30] K. Keren, Z. Pincus, G. M. Allen, E. L. Barnhart, G. Marriott, A. Mogilner,and J. A. Theriot. Mechanism of shape determination in motile cells.Nature, 453(7194):475–480, 2008. → pages 10[31] J. R. Kuhn and T. D. Pollard. Real-time measurements of actin filamentpolymerization by total internal reflection fluorescence microscopy.Biophysical journal, 88(2):1387–1402, 2005. → pages 6[32] A. Levchenko and P. A. Iglesias. Models of eukaryotic gradient sensing:application to chemotaxis of amoebae and neutrophils. Biophysical journal,82(1):50–63, 2002. → pages 21, 22[33] H. Levine, D. Kessler, and W. Rappel. Directional sensing in eukaryoticchemotaxis: a balanced inactivation model. Proceedings of the NationalAcademy of Sciences, 103(26):9761–6, 2006. → pages 27, 99, 100, 103, 114[34] L. Li, S. F. Nørrelykke, and E. C. Cox. Persistent cell motion in the absenceof external signals: a search strategy for eukaryotic cells. PLoS One, 3(5):e2093. doi:10.1371/journal.pone.0002093, 2008. → pages 11[35] B. Lin, W. R. Holmes, C. J. Wang, T. Ueno, A. Harwell,L. Edelstein-Keshet, T. Inoue, and A. Levchenko. Synthetic spatially gradedRac activation drives cell polarization and movement. Proceedings of theNational Academy of Sciences, 109(52):E3668–E3677, 2012. → pages 122[36] J. Lutkenhaus and and S. Addinall. Bacterial cell division and the Z ring.Annual review of biochemistry, 66(1):93–116, 1997. → pages 13[37] L. Ma, R. Rohatgi, and M. W. Kirschner. The Arp2/3 complex mediatesactin polymerization induced by the small GTP-binding protein Cdc42.Proceedings of the National Academy of Sciences, 95(26):15362–15367,1998. → pages 6[38] L. Ma, C. Janetopoulos, L. Yang, P. N. Devreotes, and P. A. Iglesias. Twocomplementary, local excitation, global inhibition mechanisms acting inparallel can explain the chemoattractant-induced regulation of PI (3, 4, 5) P3response in Dictyostelium cells. Biophysical journal, 87(6):3764–3774,2004. → pages 22135[39] Y. T. Maeda, J. Inose, M. Y. Matsuo, S. Iwaya, and M. Sano. Orderedpatterns of cell shape and orientational correlation during spontaneous cellmigration. PloS one, 3(11):e3734. doi:10.1371/journal.pone.0003734, 2008.→ pages 11[40] A. F. Mare´e, A. Jilkine, A. Dawes, V. A. Grieneisen, andL. Edelstein-Keshet. Polarization and movement of keratocytes: a multiscalemodelling approach. Bulletin of mathematical biology, 68(5):1169–1211,2006. → pages 6[41] A. F. Mare´e, V. A. Grieneisen, and L. Edelstein-Keshet. How cells integratecomplex stimuli: the effect of feedback from phosphoinositides and cellshape on cell polarization and motility. PLoS computational biology, 8(3):e1002402. doi:10.1371/journal.pcbi.1002402, 2012. → pages 126[42] W. Margolin. Themes and variations in prokaryotic cell division. FEMSmicrobiology reviews, 24(4):531–548, 2000. → pages 13[43] M. A. Mata, M. Dutot, L. Edelstein-Keshet, and W. R. Holmes. A model forintracellular actin waves explored by nonlinear local perturbation analysis.Journal of theoretical biology, 334:149–161, 2013. → pages 49, 51, 53, 54,127[44] M. Matecic, D. L. Smith Jr, X. Pan, N. Maqani, S. Bekiranov, J. D. Boeke,and J. S. Smith. A microarray-based genetic screen for yeast chronologicalaging factors. PLoS genetics, 6(4):e1002044.doi:10.1371/journal.pcbi.1002044, 2010. → pages 10[45] H. Meinhardt. Orientation of chemotactic cells and growth cones: modelsand mechanisms. Journal of Cell Science, 112(1):2867–74, 1999. → pages26, 90, 93, 113[46] H. Meinhardt and P. A. de Boer. Pattern formation in Escherichia coli: amodel for the pole-to-pole oscillations of Min proteins and the localizationof the division site. Proceedings of the National Academy of Sciences, 98(25):14202–14207, 2001. → pages 25, 83, 85[47] A. Millius, S. N. Dandekar, A. R. Houk, and O. D. Weiner. Neutrophilsestablish rapid and robust WAVE complex polarity in an actin-dependentfashion. Current Biology, 19(3):253–259, 2009. → pages 20, 49[48] A. Mogilner, J. Allard, and R. Wollman. Cell polarity: quantitative modelingas a tool in cell biology. Science, 336(6078):175–179, 2012. → pages 16136[49] Y. Mori, A. Jilkine, and L. Edelstein-Keshet. Wave-pinning and cell polarityfrom a bistable reaction-diffusion system. Biophysical journal, 94(9):3684–3697, 2008. → pages 20, 21, 23, 41, 42, 43, 63, 85, 110, 120[50] A. Narang, K. Subramanian, and D. Lauffenburger. A mathematical modelfor chemoattractant gradient sensing based on receptor-regulated membranephospholipid signaling dynamics. Annals of biomedical engineering, 29(8):677–691, 2001. → pages 24[51] M. Neilson, J. Mackenzie, S. Webb, and R. Insall. Modeling cell movementand chemotaxis using pseudopod-based feedback. SIAM Journal onScientific Computing, 33(3):1035–57, 2011. → pages 113[52] M. P. Neilson, D. M. Veltman, P. J. van Haastert, S. D. Webb, J. A.Mackenzie, and R. H. Insall. Chemotaxis: a feedback-based computationalmodel robustly predicts multiple aspects of real cell behaviour. PLoSbiology, 9(5):e1000618. doi:10.1371/journal.pbio.1000618, 2011. → pages26, 90, 93, 113, 126[53] M. D. Onsum and C. V. Rao. Calling heads from tails: the role ofmathematical modeling in understanding cell polarization. Current opinionin cell biology, 21(1):74–81, 2009. → pages 3, 20, 22[54] M. Otsuji, S. Ishihara, K. Kaibuchi, A. Mochizuki, S. Kuroda, et al. A massconserved reaction–diffusion system captures properties of cell polarity.PLoS computational biology, 3(6):e108. doi:10.1371/journal.pcbi.0030108,2007. → pages 24, 76, 82, 112, 130[55] E. M. Ozbudak, A. Becskei, and A. van Oudenaarden. A system ofcounteracting feedback loops regulates Cdc42p activity during spontaneouscell polarization. Developmental cell, 9(4):565–571, 2005. → pages 13[56] C. A. Parent and P. N. Devreotes. A cell’s sense of direction. Science, 284(5415):765–770, 1999. → pages 2[57] H. O. Park and E. Bi. Central roles of small GTPases in the development ofcell polarity in yeast and beyond. Microbiology and Molecular BiologyReviews, 71(1):48–96, 2007. → pages 13[58] T. D. Pollard, L. Blanchoin, and R. D. Mullins. Molecular mechanismscontrolling actin filament dynamics in nonmuscle cells. Annual review ofbiophysics and biomolecular structure, 29(1):545–576, 2000. → pages 6, 7137[59] W. H. Press. Numerical recipes 3rd edition: The art of scientific computing.Cambridge university press, 2007. → pages 33[60] A. J. Ridley. Rho GTPases and actin dynamics in membrane protrusions andvesicle trafficking. Trends in cell biology, 16(10):522–529, 2006. → pages 6[61] L. Rothfield, S. Justice, and J. Garcia-Lara. Bacterial cell division. Annualreview of genetics, 33(1):423–448, 1999. → pages 13[62] Y. Sakumura, Y. Tsukada, N. Yamamoto, and S. Ishii. A molecular modelfor axon guidance based on cross talk between Rho GTPases. Biophysicaljournal, 89(2):812–822, 2005. → pages 24[63] A. T. Sasaki, C. Chun, K. Takeda, and R. A. Firtel. Localized Ras signalingat the leading edge regulates PI3K, cell polarity, and directional cellmovement. The Journal of cell biology, 167(3):505–518, 2004. → pages 27[64] B. D. Slaughter, S. E. Smith, and R. Li. Symmetry breaking in the life cycleof the budding yeast. Cold Spring Harbor perspectives in biology, 1(3):doi:10.1101/cshperspect.a003384, 2009. → pages 13[65] S. Soh, M. Byrska, K. Kandere-Grzybowska, and B. A. Grzybowski.Reaction-diffusion systems in intracellular molecular transport and control.Angewandte Chemie International Edition, 49(25):4170–4198, 2010. →pages 17, 18[66] K. Subramanian and A. Narang. A mechanistic model for eukaryoticgradient sensing: spontaneous and induced phosphoinositide polarization.Journal of theoretical biology, 231(1):49–67, 2004. → pages 24[67] K. F. Swaney, C. H. Huang, and P. N. Devreotes. Eukaryotic chemotaxis: anetwork of signaling pathways controls motility, directional sensing, andpolarity. Annual review of biophysics, 39:265–289, 2010. → pages 1, 2[68] S. Tahirovic and F. Bradke. Neuronal polarity. Cold Spring Harborperspectives in biology, 1(3):doi: 10.1101/cshperspect.a001644, 2009. →pages 12[69] M. The´ry, V. Racine, M. Piel, A. Pe´pin, A. Dimitrov, Y. Chen, J.-B. Sibarita,and M. Bornens. Anisotropy of cell adhesive microenvironment governs cellinternal organization and orientation of polarity. Proceedings of the NationalAcademy of Sciences, 103(52):19771–19776, 2006. → pages 62138[70] A. Turing. The chemical basis of morphogenesis. PhilosophicalTransactions of the Royal Society of London. Series B, Biological Sciences,237(641):37–72, 1952. → pages 19, 34[71] M. Valiente and F. J. Martini. Migration of cortical interneurons relies onbranched leading process dynamics. Cell adhesion & migration, 3(3):278–280, 2009. → pages 12[72] P. Van Damme, E. Bogaert, M. Dewil, N. Hersmus, D. Kiraly,W. Scheveneels, I. Bockx, D. Braeken, N. Verpoorten, K. Verhoeven, et al.Astrocytes regulate GluR2 expression in motor neurons and theirvulnerability to excitotoxicity. Proceedings of the National Academy ofSciences, 104(37):14825–14830, 2007. → pages 10[73] A. B. Verkhovsky, T. M. Svitkina, and G. G. Borisy. Self-polarization anddirectional motility of cytoplasm. Current Biology, 9(1):11–20, 1999. →pages 12[74] M. G. Vicker. Eukaryotic cell locomotion depends on the propagation ofself-organized reaction-diffusion waves and oscillations of actin filamentassembly. Experimental cell research, 275(1):54–66, 2002. → pages 20, 49[75] M. G. Vicker. F-actin assembly in Dictyostelium cell locomotion and shapeoscillations propagates as a self-organized reaction–diffusion wave. FEBSletters, 510(1):5–9, 2002. → pages 20, 49[76] R. Wedlich-Soldner, S. C. Wai, T. Schmidt, and R. Li. Robust cell polarity isa dynamic state established by coupling transport and GTPase signaling.The Journal of cell biology, 166(6):889–900, 2004. → pages 13[77] O. D. Weiner, G. Servant, M. D. Welch, T. J. Mitchison, J. W. Sedat, andH. R. Bourne. Spatial control of actin polymerization during neutrophilchemotaxis. Nature cell biology, 1(2):75–81, 1999. → pages 9[78] O. D. Weiner, W. A. Marganski, L. F. Wu, S. J. Altschuler, and M. W.Kirschner. An actin-based wave generator organizes cell motility. PLoSbiology, 5(9):e221. doi:10.1371/journal.pbio.0050221, 2007. → pages 20,49[79] M. D. Welch, A. Iwamatsu, and T. J. Mitchison. Actin polymerization isinduced by Arp2/3 protein complex at the surface of Listeriamonocytogenes. Nature, 385(6613):265–269, 1997. → pages 6139[80] K. Wong, O. Pertz, K. Hahn, and H. Bourne. Neutrophil polarization:spatiotemporal dynamics of RhoA activity support a self-organizingmechanism. Proceedings of the National Academy of Sciences, 103(10):3639–3644, 2006. → pages 10[81] Y. Xiong, C.-H. Huang, P. A. Iglesias, and P. N. Devreotes. Cells navigatewith a local-excitation, global-inhibition-biased excitable network.Proceedings of the National Academy of Sciences, 107(40):17079–17086,2010. → pages 20[82] J. Xu, F. Wang, A. Van Keymeulen, P. Herzmark, A. Straight, K. Kelly,Y. Takuwa, N. Sugimoto, T. Mitchison, and H. R. Bourne. Divergent signalsand cytoskeletal assemblies regulate self-organizing polarity in neutrophils.Cell, 114(2):201–214, 2003. → pages 10[83] 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 networkreorganization breaks symmetry at the cell rear to spontaneously initiatepolarized cell motility. The Journal of cell biology, 178(7):1207–1221,2007. → pages 11, 12140Appendix AURL referencesFigure 2.3(a): This image is from the classic 16 mm film shot in the 1950s by thelate David Rogers at Vanderbilt University.Andre´s Trevin˜o, 2006. Crawling Neutrophil Chasing a Bacterium. [onlinevideo]. Available at: xh-bkiv c [Accessed 11February 2014].Figure 2.3(f): Micrograph of Escherichia coli by the National Institute ofHealth, United States Department of Health and Human Services. This work is inthe public domain. Available at: library.aspx.141Appendix BCrank-Nicolson methodThe Crank-Nicolson method is a finite difference method for solving PDES such asthe following:∂u∂ t = D∂ 2u∂x2 + f (u). (B.1)Equation B.1 has the same form as the PDES I consider for this thesis, where D isthe diffusion coefficient, f (u) is a reaction function dependent on u, and u≡ u(x, t).The Crank-Nicolson method extends to systems of PDES, but for simplicity I willonly consider one equation for this derivation.Let x ∈ [0,L] represent the spatial domain and t ∈ [0,T ] represent the temporaldomain. We discretize these domains by dividing them into equally spaced inter-vals. Space is divided into I intervals of length ∆x, and time is divided into Nintervals of length ∆t. With this discretization, let xi = i ·∆x and tn = n ·∆t, wherethe indices have the values i = 0,1, . . . I and t = 0,1, . . .N. We represent u(xi, tn) asuni . Note that xI = I ·∆x = L and tN = N ·∆t = T .Under the Crank-Nicolson method, our approximation for the solution of u isbased on two other finite difference methods: the trapezoidal rule for time and142central difference in space. We have:∂u∂ t =un+1i −uni∆t,∂ 2u∂x2 =(uni+1−2uni +uni−1)+(un+1i+1 −2un+1i +un+1i−1 )2(∆x)2.Substituting the above two equations into Equation B.1 results in:un+1i −uni∆t=D2(∆x)2((uni+1−2uni +uni−1)+(un+1i+1 −2un+1i +un+1i−1 ))+ f (uni ).Let r = D∆t/2(∆x)2. We then rearrange the terms with superscripts n and n+1 toopposite sides of the equation:un+1i − r(un+1i+1 −2un+1i +un+1i−1 ) = uni + r(uni+1−2uni +uni−1)+∆t f (uni ).Finally, grouping terms according to the subscripts i−1, i, and i+1 leads to:− run+1i+1 +(1+2r)un+1i − run+1i−1 = runi+1 +(1−2r)uni + runi−1 +∆t f (uni ). (B.3)Equation B.3 provides equations for the interior points of the domain, i= 1 . . . I−1.To find the boundary equations for indices i = 0 and i = I, we need to incorporatethe “no flux” boundary conditions. No flux boundary conditions in the PDE modelhave the following form:∂u∂x∣∣∣∣x=0= 0 =∂u∂x∣∣∣∣x=L. (B.4)To discretize Equation B.4, I will introduce phantom grid points with the indicesi = −1 and i = I + 1. Then, evaluating at i = 0 (corresponding to x = 0), or i = I143(corresponding to x = L). The discretized boundary conditions are given by:un0−un−12∆x= 0,unI+1−unI2∆x= 0,and correspond to un−1 = un0 and unI+1 = unI . Note that the boundary conditions holdfor all times t, so the index n can take any value on its interval (n ∈ [0,N]). If wenow consider the form Equation B.3 has for the boundary points x0 = 0 and xI = L,we have:−run+11 +(1+2r)un+10 − run+1−1 = run1 +(1−2r)un0 + run−1 +∆t f (un1),−run+1I+1 +(1+2r)un+1I − run+1I−1 = runI+1 +(1−2r)unI + runI−1 +∆t f (unI ).Substituting the values un1 = un0 and unI = unI+1 will remove the phantom grid point,and give the boundary equations:(1+2r)un+10 −2run+1−1 = (1−2r)un0 +2run−1 +∆t f (un1), (B.7a)−2run+1I+1 +(1+2r)uI+1i = 2runI+1 +(1−2r)unI +∆t f (unI ). (B.7b)In matrix form, the interior equations (represented by Equation B.3) combined144with the boundary equations (Equation B.7a and Equation B.7b) become:1+2r −2r 0 · · · 0−r 1+2r −r · · · 00 −r 1+2r · · ·...0 0.... . . −r0 0 · · · −2r 1+2run+11un+12...un+1I−1un+1I=1−2r 2r 0 · · · 0r 1−2r r · · · 00 r 1−2r · · ·...0 0.... . . r0 0 · · · 2r 1−2run1un2...unI−1unI+∆tf (un1)f (un2)...f (unI−1)f (unI ).We write the above matrix system as:Aun+1 = Bun + f (un). (B.9)The indices of u correspond to the spatial discretization i = 0,1, . . . I. Since nrepresents the current timestep, we can solve Equation B.9 for un+1 to get thesolution for the next timestep, n+1:un+1 = A−1(Bun + f (un)). (B.10)145Appendix CTuring AnalysisA Turing analysis helps us to detect diffusion-driven instabilities in certain PDEsystems. For a diffusion-driven instability to occur, the steady state must be linearlystable to perturbations when diffusion is absent, but become unstable when diffu-sion is introduced. In this appendix, I will derive conditions leading to diffusion-driven instability. Consider the following two component system:∂u∂ t = Du∂ 2u∂x2 + f (u,v), (C.1a)∂v∂ t = Dv∂ 2v∂x2 +g(u,v). (C.1b)The first set of conditions I will derive are those confirming stability of the HSS inthe absence of diffusion. That is, in vector notation we are looking at the system:ddt(uv)=(f (u,v)g(u,v))(C.2)Let the HSS of the system be given by (u∗,v∗). By definition of the steady state,f (u∗,v∗) = g(u∗,v∗) = 0. Note that it is not necessarily trivial to solve for the HSS.We investigate the stability of (u∗,v∗) by observing the behaviour of small,inhomogeneous perturbations to the steady state. Let these perturbations be given146by:w(x, t) =(u′(x, t)v′(x, t))=(u(x, t)−u∗v(x, t)− v∗). (C.3)We substitute the perturbation Equation C.3 into Equation C.2 to get:w =ddt(u−u∗v− v∗)=(f (u−u∗,v− v∗)g(u−u∗,v− v∗)). (C.4)Now, linearizing Equation C.4 about the HSS (u∗,v∗) leads to the following system:ddt(u−u∗v− v∗)= J(u−u∗v− v∗), (C.5)where J is the Jacobian, and has the formJ =[fu fvgu gv]∣∣∣∣∣ss. (C.6)Since Equation C.5 is a linear system, we can solve to obtain solutions of theformw = α1v1 exp(σ1t)+α2v2 exp(σ2t), (C.7)where α arises from the initial conditions, and v1,2 are the eigenvectors of J corre-sponding to the eigenvalues σ1,2. The eigenvalues are given by:σ1,2 =tr(J)±√tr(J)2−4det(J)2. (C.8)For stable solutions, we require that the perturbation w diminishes to 0 as t → ∞,and so σ1,2 < 0. Thus, we have the following two conditions to guarantee stabilityof the HSS (note that the presence of an ∗ symbol implies “evaluated at the steadystate (u∗,v∗)”):1. Det(J) = f ∗u g∗v− f∗v g∗u > 0,2. Tr(J) = f ∗u +g∗v < 0.147Now that we have regions in which the HSS (u∗,v∗) is stable in the absenceof diffusion, we repeat the linear stability analysis with diffusion included. Thistime, we substitute the small, inhomogeneous perturbations, given by w(x, t) inEquation C.3, into the full PDE system in Equation C.1. Linearizing the resultabout the HSS results in the following system:ddt(u−u∗v− v∗)=[fu fvgu gv]∣∣∣∣∣ss(u−u∗v− v∗)+[Du 00 Dv]∂ 2∂x2(u−u∗v− v∗). (C.9)More simply, we can write wt = Jw+Dwxx. We now have a linear system, andsolutions can be given by:w = αeσtcos(qx). (C.10)As before, this expression for w represents the form for the small perturbations.Therefore, the amplitude α must be small. The eigenvalues are represented byσ , and q = npi/L. Substituting Equation C.10 into Equation C.9 results in thefollowing system:α1σ = f ∗u α1 + f ∗v α2−Duq2α1,α2σ = g∗uα1 +g∗vα2−Dvq2α2.We will rearrange the above system in terms of α1,2:α1(σ − f ∗u +Duq2)+α2(− f ∗v ) = 0,α1(−g∗u)+α2(σ −g∗v +Dvq2)+ = 0,148Since we want a nontrivial perturbation, we require that:det(σ − f ∗u +Duq2 − f ∗v−g∗u σ −g∗v +Dvq2)= 0= (σ − f ∗u +Duq2)(σ −g∗v +Dvq2)−g∗u f ∗v= σ2 +σ( f ∗u −Duq2 +g∗v−Dvq2)+(( f ∗u −Duq2)(g∗v−Dvq2)− f ∗v g∗u).We want to determine if there is some value of q that gives rise to Re(σ)> 0, asituation in which perturbations grow. Notice that this characteristic equation canalso be written as σ2 + tr(JD)σ +det(JD) = 0, where JD has the form:JD = J−Dq2 =[fu fvgu gv]∣∣∣∣∣ss−[Duq2 00 Dvq2]=[f ∗u −Duq2 f ∗vg∗u g∗v−Dvq2].Similarly to our analysis without diffusion:σ = tr(JD)±√tr(JD)2−4det(JD)2,and the following conditions lead to stability:1. det(JD) =(( f ∗u −Duq2)(g∗v−Dvq2)− f ∗v g∗u)> 0,2. Tr(JD) = ( f ∗u −Duq2 +g∗v−Dvq2) < 0.Since the goal is to show that the system in Equation C.1 is unstable when diffusionis present, we must have that either condition (1) or (2) is not satisfied. Noticethat condition (2) is always true, because I have already shown that Tr(J) = f ∗u +g∗v < 0, and Du,Dv,q2 are positive quantities. We now question whether thereare any values for q2 for which condition (1) is is not true. That is, det(J) =(( f ∗u −Duq2)(g∗v−Dvq2)− f ∗v g∗u)< 0.To find conditions for which det(JD) < 0, we will first expand this expressionto see that it represents an upward facing parabola:det(JD) = DuDv(q2)2− ( f ∗u Dv +g∗vDu)q2 +det(J). (C.14)149To find the minimum of the parabola, we will differentiate det(JD) with respectto q2:d(det(JD))dq2= 2DuDvq2− f ∗u Dv +g∗vDu,⇒ q2min =12(f ∗uDu+g∗vDv).For real solutions of qmin, we must have ( f ∗u /Du + g∗v/Dv) > 0. Substituting q2minback into the expression for det(JD) in Equation C.14, and requiring that det(JD)<0, leads to the following inequality:det(J)−Dug∗v +Dv f∗u4DuDv< 0.Rearranging the above inequality leads to our final condition for diffusion-driveninstability:(Dug∗v +Dv f∗u ) >√4DuDvdet(J) > 0.To summarize, we have three conditions for diffusion-driven instability:1. det(J) = f ∗u g∗v− f∗v g∗u > 0,2. tr(J) = f ∗u +g∗v < 0,3. (Dug∗v +Dv f∗u ) >√4DuDvdet(J) > 0.150Appendix DDerivation of the ThomasalgorithmWe are interested in solving the tridiagonal matrix system that results from the dis-cretization of our RDES. The Thomas algorithm is a method for solving tridiagonalmatrix systems in O(n) operations as opposed to the O(n3) operations which arerequired for Gaussian elimination. Another advantage of the Thomas algorithm isthe ease of implementing different boundary conditions. The following is a deriva-tion for the Thomas algorithm, with constant flux boundary conditions.We are considering a tridiagonal matrix system of the form:b1 c1 0 · · · 0a2 b2 c2 · · · 00 a3 b3 · · ·...0 0.... . . cn−10 0 · · · an bnu1u2...un−1un=d1d2...dn−1dn.From the matrix equation, we have a system of equations of the form:anun−1 +bnun + cnun+1 = dn. (D.1)Note that a1 = 0 and cn = 0. Now also assume that:151un−1 = γnun +βn. (D.2)This assumption results from solving the nth matrix equation for un and sub-stituting this expression into the n+ 1th equation. This substitution results in anequation involving only the variables un+1 and un+2. Here we can see un+1 recur-sively defined as a function of un+2 and unknown coefficients γn and βn, which willbe determined next.D.1 Forward sweepTo find the coefficients γn, βn we will substitute our assumption Equation D.2 intothe general matrix equation Equation D.1 as follows:an(γnun +βn)+bnun + cnun+1 = dn(anγn +bn)un + cnun+1 = dn−anβn(anγn +bn)un = dn−anβn− cnun+1un =dn−anβnanγn +bn+−cnanγn +bnun+1.The final expression has the same form as Equation D.2. Thus, we see that thecoefficients are given by the following expression:γn+1 =−cnanγn +bn,βn+1 =dn−anβnanγn +bn.and are recursively defined as functions of γn, βn. We must consider the bound-ary conditions to get values for γ1, β1. Here, we are considering flux boundaryconditions, so we have the case where:D∂u∂x∣∣∣∣0,L= J0. (D.3)152After discretizing, this condition gives:u1−u0 =J0∆xD= F0,u0 = u1−F0.This final equation matches the recursive definition Equation D.2, and we seethat γ1 = 1 and β1 =−F0.D.2 Backward sweepIn order to determine the variables for the vector u we will use Equation D.2 forun. From the forward sweep, all coefficients γn, βn have been determined. Thebackward sweep begins with a value for un, which is determined from the rightboundary condition. As before, we have the flux boundary condition given byEquation D.3, which after discretizing yields:un−un−1 = F0, (D.4)and we have the recursive definition:un−1 = γnun +βn. (D.5)We can substitute the recursive definition Equation D.5 into Equation D.4 toget:un− γnun−βn = F0,un =F0 +βn1− γn.As we already have values for γn and βn, this expression gives the starting valuefor the recursive formula Equation D.2. We can now determine our solution vector,u.153


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items