Science, Faculty of
Mathematics, Department of
DSpace
UBCV
Liu, Yue
Master of Science - MSc
University of British Columbia
Small GTPases are a family of signalling proteins that regulates cell shape through actin assembly. Actin filaments are a major component of the cellular cytoskeleton, a supporting structure of the cell. Patterns of GTPase activity are related to a variety of cell behaviors, including cell polarization, formation of protrusions such as filopodia, and actin waves. These phenomena are in turn related to cell migration, cancer metastasis and viral infection. We seek to investigate pattern formation in the context of GTPase activity through a variant of the reaction-diffusion partial differential equation model, known as the Wave Pinning Model. By determining the possible spatio-temporal patterns, as well as their required conditions, we can better understand the cellular behaviors corresponding to these patterns. We first introduce and motivate the wave pinning model and its extensions involving F-actin feedback and source-sink terms. Next, we explore the behavior of the model with numerical simulations. We have identified patterns of localized spots, travelling waves, pulses, and others. We interpret some of these patterns as cell polarization, filopodia formation, and actin waves, as well as other complex behaviors that do not closely resemble common experimental observations. Next, we delineate distinct parameter regimes with Turing analysis and local perturbation analysis (LPA). Turing analysis is a classical method for determining linear stability of homogeneous steady states. LPA is a recently developed technique that is able to detect instabilities to perturbations of finite size that cannot be detected by Turing analysis. LPA has the caveat that it describes only the limiting behavior of the system as the ratio of diffusion coefficient of the two forms of GTPase goes to zero. We also study the effects stemming from the interaction of GTPase dynamics and the deforming cell boundary with a simplified wave pinning model coupled with cell size, and with Cellular Potts Model simulations coupling the full wave pinning model with the cell boundary. The results of the analysis allow us to comment on the impact of cell geometry, the various terms in our model and their interactions on pattern formation and cell behavior.
Analysis of pattern formation inreaction-diffusion models for cellpolarizationbyYue LiuB.Math. Hons., University of Waterloo, 2017A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2019c© Yue Liu 2019The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled: Analysis of pattern formation in reaction-diffusion models for cell polarization submitted by Yue Liu in partial fulfillment of the requirements for the degree of Master of Science in Mathematics Examining Committee: Leah Edelstein-Keshet, Mathematics Supervisor Eric Cytrynbaum, Mathematics Supervisory Committee Member Wayne Nagata, Mathematics Supervisory Committee Member iiAbstractSmall GTPases are a family of signalling proteins that regulates cell shapethrough actin assembly. Actin filaments are a major component of the cel-lular cytoskeleton, a supporting structure of the cell. Patterns of GTPaseactivity are related to a variety of cell behaviors, including cell polariza-tion, formation of protrusions such as filopodia, and actin waves. Thesephenomena are in turn related to cell migration, cancer metastasis and viralinfection. We seek to investigate pattern formation in the context of GTPaseactivity through a variant of the reaction-diffusion partial differential equa-tion model, known as the Wave Pinning Model. By determining the possiblespatio-temporal patterns, as well as their required conditions, we can betterunderstand the cellular behaviors corresponding to these patterns.We first introduce and motivate the wave pinning model and its exten-sions involving F-actin feedback and source-sink terms. Next, we explorethe behavior of the model with numerical simulations. We have identifiedpatterns of localized spots, travelling waves, pulses, and others. We inter-pret some of these patterns as cell polarization, filopodia formation, andactin waves, as well as other complex behaviors that do not closely resemblecommon experimental observations.Next, we delineate distinct parameter regimes with Turing analysis andlocal perturbation analysis (LPA). Turing analysis is a classical method fordetermining linear stability of homogeneous steady states. LPA is a recentlydeveloped technique that is able to detect instabilities to perturbations offinite size that cannot be detected by Turing analysis. LPA has the caveatthat it describes only the limiting behavior of the system as the ratio ofdiffusion coefficient of the two forms of GTPase goes to zero.We also study the effects stemming from the interaction of GTPase dy-namics and the deforming cell boundary with a simplified wave pinningmodel coupled with cell size, and with Cellular Potts Model simulationscoupling the full wave pinning model with the cell boundary. The resultsof the analysis allow us to comment on the impact of cell geometry, thevarious terms in our model and their interactions on pattern formation andcell behavior.iiiLay SummaryGTPases are a family of signalling proteins that controls cell shape. Bychanging shape, the cell can move and form protrusions such as filopodia.Cell movements are essential for embryonic development and cancer metas-tasis, while filopodia serve a variety of purposes. Therefore, studying thedynamics of GTPases helps us understand these biological phenomena.This project focuses on pattern formation in a reaction-diffusion modeldescribing the changes in levels of activity of GTPase across space and time.I ran numerical simulations to determine the possible patterns in this model.I found a pattern of spikes, which corresponds to filopodia formation, as wellas a pattern of moving waves resembling wave-like cell activities observedin experiments. Next, I analyzed the model with both classical and recentmathematical techniques to find the conditions required by these patterns.Finally, I discussed the interactions between GTPase and cell shape.ivPrefaceThis project is conducted by the author, Yue Liu, under the guidance ofProfessor Leah Edelstein–Keshet.The models examined in this thesis are related to existing works by Moriet al [30], Holmes et al [16] and Verschueren & Champneys [43], and someof the results are reproduced as benchmarks for my analysis and numericalmethods. The numerical simulations in Chapter 3 are produced with myown original code, and the text in this thesis is my own writing.Figure 2.2 was produced by Leah Edelstein–Keshet, used with permis-sion.The simulations shown in Figure 6.3 were run by Elisabeth Rens, whotogether with Zachary Pellegrin developed the program used for the simu-lation. The figures are used here with their permission.The code given in Appendix C are based on skeleton code provided byDr. Timm Trekatis as a part of MATH 521 course material.All other material is my original work, and has not yet been published.Some sections come from parts of my past course projects: parts of Sec. 3.1and 3.8 from MATH 521, and parts of Ch. 5 from MATH 553.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 GTPase, the actin cytoskeleton and cell motility . . . . . . . 11.2 Filopodia and virus invasion . . . . . . . . . . . . . . . . . . 21.3 Actin waves . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.4 Modelling motivation . . . . . . . . . . . . . . . . . . . . . . 31.5 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . 51.6 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . 62 Mathematical models . . . . . . . . . . . . . . . . . . . . . . . 82.1 Representing the cell with two- or one-dimensional domains . 82.2 Wave pinning model . . . . . . . . . . . . . . . . . . . . . . . 92.2.1 Non-dimensionalization . . . . . . . . . . . . . . . . . 102.3 Actin feedback extension . . . . . . . . . . . . . . . . . . . . 112.4 Non-conservative extension . . . . . . . . . . . . . . . . . . . 122.5 Combining the models . . . . . . . . . . . . . . . . . . . . . . 132.6 Mass conservation and well-mixed models . . . . . . . . . . . 13viTable of Contents3 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . 163.1 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . 163.2 Selection of simulation parameters . . . . . . . . . . . . . . . 173.3 Simulations of wave pinning in 1D . . . . . . . . . . . . . . . 183.4 Simple analysis of wave pinning . . . . . . . . . . . . . . . . 193.5 Simulations of actin feedback in 1D . . . . . . . . . . . . . . 223.6 Simulations of non-conservative model in 1D . . . . . . . . . 233.7 Simulations of the combined model in 1D . . . . . . . . . . . 243.8 Simulations of wave pinning and non-conservative extensionin 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.9 Simulations with spatially inhomogeneous parameters . . . . 264 Turing analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1 Overview of Turing analysis . . . . . . . . . . . . . . . . . . 364.2 Turing analysis of wave pinning model and non-conservativeextension in 1D . . . . . . . . . . . . . . . . . . . . . . . . . 374.2.1 Wave pinning . . . . . . . . . . . . . . . . . . . . . . 394.2.2 Non-conservative extension . . . . . . . . . . . . . . . 414.2.3 Predicting the wave length of pattern formed underthe non-conservative model . . . . . . . . . . . . . . . 415 Local perturbation analysis . . . . . . . . . . . . . . . . . . . 465.1 Overview of LPA . . . . . . . . . . . . . . . . . . . . . . . . . 465.2 LPA for Wave pinning . . . . . . . . . . . . . . . . . . . . . . 495.2.1 More on Wave pinning . . . . . . . . . . . . . . . . . 545.3 LPA for non-conservative extension . . . . . . . . . . . . . . 545.4 LPA for actin feedback extension . . . . . . . . . . . . . . . . 585.5 LPA for the combined Model . . . . . . . . . . . . . . . . . . 636 Interaction with changing domain . . . . . . . . . . . . . . . 646.1 Well-mixed wave pinning model in 1D with dynamic domain 646.1.1 Contraction case (Rho) . . . . . . . . . . . . . . . . . 666.1.2 Expansion case (Rac, Cdc42) . . . . . . . . . . . . . . 706.2 PDE simulations with deforming boundary using the CellularPotts Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 717 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.1 Summary of model behaviors . . . . . . . . . . . . . . . . . . 767.2 Biological and mathematical implications . . . . . . . . . . . 787.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 80viiTable of ContentsBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83AppendicesA Code for 1D simulations . . . . . . . . . . . . . . . . . . . . . 88B XPPAUT programs for bifurcation analysis . . . . . . . . . 91C Code for 2D simulations . . . . . . . . . . . . . . . . . . . . . 93viiiList of Tables3.1 Parameter values for numerical simulations . . . . . . . . . . 185.1 Summary of the wave pinning model regimes . . . . . . . . . 525.2 Summary of the non-conservative model regimes . . . . . . . 565.3 Summary of the actin feedback model regimes . . . . . . . . . 596.1 Parameter values for the dynamical cell size model . . . . . . 67ixList of Figures2.1 1D domain illustration . . . . . . . . . . . . . . . . . . . . . . 92.2 Graphical representation of the four models . . . . . . . . . . 143.1 Wave pinning simulation . . . . . . . . . . . . . . . . . . . . . 203.2 Negative wave pinning . . . . . . . . . . . . . . . . . . . . . . 213.3 Steady state of the wave pinning model . . . . . . . . . . . . 223.4 Numerical simulations of the actin feedback model . . . . . . 283.5 Exotic patterns observed under the actin feedback model . . 293.6 Numerical simulations of the non-conservative model . . . . . 303.7 Steady state pattern of the non-conservative model . . . . . . 313.8 Numerical simulations of the combined model . . . . . . . . . 323.9 Numerical simulations of the combined model (2) . . . . . . . 333.10 2D simulations of the wave pinning and non-conservative model 343.11 Simulations for the combined model with spatially-dependentparameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1 Plot of activation vs inactivation terms . . . . . . . . . . . . . 404.2 Turing regime of the non-conservative model . . . . . . . . . 424.3 Dispersion relation for the non-conservative model . . . . . . 434.4 Transition from precursor pattern to the final pattern for thenon-conservative model . . . . . . . . . . . . . . . . . . . . . 445.1 Illustration of the idea behind LPA . . . . . . . . . . . . . . . 485.2 One-parameter bifurcation diagrams of the LPA and well-mixed wave pinning system . . . . . . . . . . . . . . . . . . . 505.3 Two-parameter bifurcation diagrams of the LPA and well-mixed wave pinning system . . . . . . . . . . . . . . . . . . . 515.4 Illustration of LPA regimes for the wave pinning model withplots of f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.5 Bifurcation diagrams of the LPA and well-mixed non-conservativemodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.6 Bifurcation diagrams of the actin feedback extension . . . . . 61xList of Figures5.7 Details of the two-parameter bifurcation diagrams for theactin feedback extension . . . . . . . . . . . . . . . . . . . . . 625.8 Bifurcation diagram of the combined LPA system . . . . . . . 636.1 Bifurcation diagrams for the dynamic cell size model, con-traction case . . . . . . . . . . . . . . . . . . . . . . . . . . . 696.2 Bifurcation diagrams for the the dynamic cell size model, ex-pansion case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 726.3 Simulations of the coupled CPM-GTPase model . . . . . . . . 757.1 Comparison of LPA and Turing regimes for the non-conservativemodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81xiAcknowledgementsI would like to thank my supervisor, Professor Leah Edelstein–Keshet, forher guidance as well as encouragement and support throughout this project.She has taught me a wide range of knowledge, skill and experience thatwill surely benefit me throughout my career, and I’m deeply moved by herpatience and kindness.I would like to thank Elisabeth Rens and Zachary Pellegrin for developingand providing the CPM simulations, and kindly allowing me to use theresults in Chapter 6.I would also like to thank my fellow students, especially Cole Zmur-chok, Alastair Jamieson–Lane and Frederic Paquin–Lefebvre, as well as allmembers of our research group for providing me with ideas during our dis-cussions, as well as advice and support. Moreover, I would like to thankmy course instructors, in particular Professors Michael Ward, Wayne Na-gata, Eric Cytrynbaum, Timm Treskatis and Neil Balmforth, who taughtme much of the knowledge I used in this project and introduced me tofascinating research.Finally, I would like to acknowledge that this work is supported by theNatural Sciences and Engineering Research Council of Canada (NSERC)under the Alexander Graham Bell Canada Graduate Scholarship to me andDiscovery Grant to Leah Edelstein-Keshet.xiiDedicationTo my parentsxiiiChapter 1IntroductionCell polarization is a fascinating biological phenomenon where cells breaktheir internal symmetry and form distinct internal regions in response toexternal signals. We are especially interested in the role of GTPase activity,a family of signalling proteins that regulate cell shape as well as cell polarity.Understanding the possible dynamics of GTPase activity has broad appli-cation areas including cell motility, cancer metastasis and viral infection.There have been many recent developments in quantitative biology, specif-ically mathematical modelling. Most models are ODE, PDE, or agent-based,and some are stochastic. Their strengths and characteristics are discussedin [29], along with examples related to cell polarity.The aim of this paper is to build and analyze a PDE model for GTPasedynamics, identify possible spatio-temporal patterns, as well as conditionsunder which these patterns are possible. This will contribute to understand-ing the essential mechanisms behind cell polarization and related phenom-ena.In this chapter, we will provide a brief overview of the biological back-grounds, provide motivations for our model, discuss potential applicationsof our research, and review existing models of GTPase activity.1.1 GTPase, the actin cytoskeleton and cellmotilitySmall GTPases are a family of signalling proteins that regulate the shape ofthe cell via polymerization of actin filaments [15]. Their main members areRac, Rho and Cdc42. They all have a membrane-bound active form, whichdiffuse much slower than their inactive form.The actin filaments, also known as F-actin, are assembled from monomersknown as G-actin [24, Chapter 19.1] and form a major component in thecellular cytoskeleton. As a whole, the cytoskeleton supports and stabilizesthe shape of the cell, anchors the cell to the extra-cellular matrix, as wellas participates in a variety of cell functions such as movement, cell division,11.2. Filopodia and virus invasionendocytosis and signalling pathways [11].By controlling cell shape, the interaction of GTPases and actin drives aform of cell motility. In the front of the cell, the GTPases Rac and Cdc42up-regulate actin assembly, which extends the cell forward. At the sametime, Rho GTPase activates acto-myosin contraction at the back of the cell,pulling the rear inward. Together this allows the entire cell to crawl forward.Polarization is required for cell migration during embryo development(well known, see for example [13]). It also goes wrong in various stages ofcancer metastasis [33]. Therefore, understanding the conditions under whichcell polarization is possible is of great interest.1.2 Filopodia and virus invasionActin is also responsible for the formation of protruding structures suchas lamellipodia and filopodia [1, Chapter 16]. Filopodia are thin, long,finger-like protrusions formed by tightly-packed bundles of actin filaments.Lamellipodia is a much larger structure containing an actin mesh, insteadof bundles. Although filopodia serve many purposes, they are utilized byseveral types of virus during the infection process [3]. These viruses includedengue [48] and rabies [47], which use filopodia to enter the cell, and HIV-1 [32], which can trigger the formation of filopodia once inside a cell, anduse it to facilitate transmission to another cell. Not all forms of GTPasesare involved in this process. In [48] the authors showed that only Rac andCdc42, but not RhoA, is relevant to virus infection.Consequently, understanding GTPase activity and its implications tofilopodia formation can suggest new ways of combating viral infections.1.3 Actin wavesIn [44], the author observed wave-like F-actin activities in amoeba cells,and hypothesized that the wave pattern, in addition to cell polarization,plays a role in cell motility. A later study [46] experimentally observedwave fronts of the assembly of F-actin, and activities of Rac GTPase and arelated protein (Hem-1) in human immune cells, and found that these wavefronts annihilate upon collision. In addition, the study found strong spatialcorrelation between Rac activity and actin waves, which is in turn stronglycorrelated with the advancement of cell edges, supporting the hypothesisthat actin waves are related to cell motility.21.4. Modelling motivationIn light of these discoveries, it would be interesting to incorporate actinassembly into our model in order to reproduce and explain the formation ofactin waves.1.4 Modelling motivationThe GTPase signalling network is highly complex. There are 40 proteinspresented in [41], and it is still not a complete list. A mechanistically accu-rate model would contain too many variables and parameters for analysis.In order to obtain a mathematical model tractable for analysis, we will worktoward a phenomenological, conceptual model that is simple yet still suffi-cient to produce the behaviors of interest.Our work is mainly based on three previous studies: [30, 43, 16]. Themost basic model, known as the wave pinning model, first proposed by [30],is a highly simplified and abstract model. This model considers only a singleGTPase in isolation, and aims to describe cell polarization. The model usesa reaction-diffusion equation to describe the spatio-temporal dynamics ofthe active and inactive forms of a GTPase. The reaction term describes therate of conversion of one form into the other, with an abstract nonlinearactivation term and a linear deactivation term. As a whole, this modelpreserves the total amount of GTPase. One important biological fact toincorporate is the vastly different diffusion coefficient for these two forms,due to the fact that the active form is membrane-bound whereas the inactiveform can move freely in the fluid cell interior. This turns out to be importantfor polarization and pattern formation. This will be shown in Ch. 4, 5.There are two other papers proposed two separate extensions to the ba-sic model, which we will analyze. The non-conservative extension proposedby [43] includes additional source and sink terms, which represents the pro-duction and decay of GTPase. This results in patterns that we will interpretas filopodia formation. The actin feedback extension proposed by [16] addsF-actin as another variable interacting with the GTPase, which producespatterns similar to actin waves.Our main goal is for each model, to identify any possible spatio-temporalpatterns in GTPase dynamics, and interpret their effect on cell shape. Thenwe determine the mathematical conditions for the existence of these pat-terns, which would correspond to the biological conditions for cell polariza-tion, filopodia formation and actin waves.This analysis will help us answer the following questions of interest:Q1. How do the source and sink effects interact with the feedback effect31.4. Modelling motivationfrom actin assembly? GTPase source/sink terms lead to localized spots[43], whereas feedback from F-actin to GTPase results in travellingwaves [16]. What kind of dynamic patterns arise when the two effectsare combined in a single GTPase model?Q2. Does locally elevated GTPase activation lead to localized pattern for-mation, or will the pattern spread farther away? Is the model consis-tent with viral infection through filopodia?Q3. Do 1D and 2D geometries yield different patterns in GTPase activity?A top-down view of a cell can be approximated by a 2D domain (seeSec. 2.1), but cells can be forced into a 1D geometry using microfluidicchannels or substrates with grooves.Q4. What impact does cell size have on pattern formation? Cells differsignificantly in size, even within the same species.Q5. How do GTPase dynamics and a dynamic domain boundary interact?Given that GTPase regulates cell shape, a system with a dynamicboundary coupled with GTPase activity might yield a more complexbehavior than with a static boundary.Q6. Can the model account for the oscillating spots in GTPase activityobserved in [39]? This was one of the motivations for combining thenon-conservative terms, which results in spots, and actin feedback,which can potentially enable the spots to oscillate.Q7. Does the presence of actin waves produce a qualitatively differentmovement pattern for the cell, compared to simple polarization? Cellmovement can be modelled as the process of cell boundary deformationdue to interaction with GTPase.Q8. Can Turing analysis (introduced in Ch. 4) predict the wavelength ofspatially periodic patterns, if such patterns exist?Q9. How do the results of local perturbation analysis (LPA, introduced inCh. 5) relate to the results of Turing analysis?Q10. How does LPA and Turing analysis perform when there are more thantwo interacting species?41.5. Related works1.5 Related worksBeside our main references [30, 43, 16] which first proposed the wave pinning(WP) model and its extensions, there are other relevant papers. In [31], theauthors used matched asymptotic analysis to explore the PDE bifurcationsof the WP model, and obtained results on the speed and final pinned positionof the wave front. In [27], the authors performed local perturbation analysis(LPA) on the actin feedback model and identified possible boundaries ofpattern-forming parameter regimes.There are many closely related models as well. In [17], the authorscompared the WP model to a related model that involves both Rac andRho GTPases (Rac enhances cell protrusion and expansion, whereas Rhoactivates cell contraction). Both have reaction terms similar to the WPmodel, but slightly modified to describe mutual inhibition. The paper founddistinct regimes of the parameter space that correspond to different cellshapes and sizes, determined by equilibrium levels of Rac and Rho. Theabove model is further extended in [17] by incorporating Cdc42 and a fewother proteins as well, with their own reaction-diffusion equations, resultingin a highly complex model. This model is able to recover wave pinningbehavior. The authors also discussed how the presence of the other proteinsinfluences LPA results.There are models that consider the coupling between GTPase and theenvironment surrounding the cell. In [19], the authors consider a Rac-Rhomodel similar to [17], but with the additional coupling to the extracellularmatrix, and investigated the formation of lamellipodia. In [49], the authorscoupled GTPase dynamics to mechanical tension of the cell, and linked anarray of cells together. This resulted in a variety of interesting collectivecell behaviors, such as waves of oscillatory contraction and relaxation prop-agating through the group of cells. This model is related to our modelin Sec. 6.1.1. In addition, this study also simulated groups of cells withdynamic shapes in 2D.In [45], the authors derived a stochastic version of the WP model bytracking the number of GTPase molecules in discretized spatial compart-ments. This model reduces to the wave pinning PDE in the limit of verylarge number of molecules. It was found that the stochastic model preservesthe wave pinning behavior when there is a large number of molecules andlow randomness. The parameter regime for wave pinning shrinks as thenumber of molecules decreases or randomness increases.All studies mentioned above are ultimately based on the WP model.There are other unrelated models for cell polarization as well. A broad51.6. Thesis overviewcomparison between these models is given in [9, 20].The Otsuji model [34] is another reaction-diffusion model for GTPasedynamics. It obeys mass conservation similar to the WP model, except witha very different reaction function. The authors focused on the behavior ofthe model under a time-dependent inhomogeneous stimulation. The usualbehavior of this model is a single peak that follows the stimulation around,and persists after the stimulation is turned off.Turing’s 1952 paper [42] is the earliest paper that studied diffusion-drivenpattern formation. Rather than focusing on a specific model, the paperdiscussed reaction-diffusion systems in general. Some of its techniques areapplied to our model in Ch. 4.Another non-specific model of pattern formation is the Gierer–Meinhardt(GM) model [12], which is a well-studied model with many related papers.The GM model consists of a short-range activator and a long-range inhibitordiffusing at different rates. This model exhibits polarization, periodic staticpatterns as well as travelling waves, depending on parameters. All of thesebehaviors are present in the wave pinning model (with extensions), shownin Ch. 3. There is an extension to the GM model given in [28].There are other very different approaches to modelling cell polarity.Some focus on proteins other than GTPases, some are more experimen-tal and focus on fitting equations to observed data, and some use agentagent-based models. Many of these are reviewed in [29].Some studies combine experimental observations with modelling. In [36],the authors proposed a model for GTPase dynamics similar to the WP modelwith different reaction and without mass conservation, and compares thesimulations to indirectly measured GTPase activity in plant root hair cells.In [39], the authors measured GTPase and actin assembly in roundwormembryos, and fitted reaction functions to measurements.Finally, there are many other directions used to study cell polarity bothin experimental methods and mathematical modelling. Some of them aresummarized in [37].1.6 Thesis overviewChapter 2 introduces the assumptions and formulations of the models indetail. We will also discuss how, in certain cases, the model obeys massconservation, which simplifies the analysis.Chapter 3 contains the numerical simulations of the models, under avariety of parameters and initial conditions, showing all known patterns. We61.6. Thesis overviewalso discuss the result of numerical experiments on determining conditionsfor existence of the patterns, which will guide the subsequent analysis.Chapter 4 focuses on Turing analysis as a way to determine the linearstability of homogeneous steady state (HSS) solutions for our PDE models.This is of interest since a HSS unstable to spatially heterogeneous perturba-tion is a sufficient condition for pattern formation. We will also discuss howTuring analysis can help determine the wave length of the spatially periodicpatterns, when these patterns exist.Chapter 5 uses local perturbation analysis (LPA), a more recent tech-nique, to detect another form of instability which requires a localized, suffi-ciently large perturbation. LPA is able to determine more distinct pattern-forming parameter regimes than Turing analysis, with the drawback thatit only works in the limit that the ratio of the diffusion coefficients of theactive and inactive forms goes to zero.Chapter 6 discusses and analyzes two ways to couple the wave pinningmodel to a dynamic domain representing the effect of GTPase on the shapeof the cell. The results allow us to more concretely comment on the con-sequence of GTPase dynamics on cell shape and the related biological phe-nomena.7Chapter 2Mathematical modelsIn this chapter, we describe a few increasingly complex reaction-diffusionPDE models for describing the dynamics of GTPase activity in some domainrepresenting a cell. The most basic one, known as wave pinning model, wasproposed by [30], which focuses on explaining cell polarization. There aretwo separate extensions to this model, the first of which is the actin feedbackmodel proposed by [16], which aims to describe a phenomenon known asactin waves. The second is the non-conservative model proposed by [43],which includes additional terms to be more realistic.We will combine these two extensions into a single model, which hasnot been considered before. We will also discuss the choice of domain, andnon-dimensionalization of the models.2.1 Representing the cell with two- orone-dimensional domainsWe will consider both 1D and 2D domains. Many kinds of cells of interesthave a flat shape, which means that the vertical thickness would be too thinto support complex patterns. As such, a 3D model likely will not provideadditional insight beyond a 2D model. For 2D, we consider square and cir-cular domains. Although GTPase activity will affect the actin cytoskeleton,therefore deforming the shape of the domain over time, this effect will notbe considered here. Instead, we assume the domain is static. We will returnto this question in Chapter 6. The boundary conditions are no-flux, sinceGTPases cannot diffuse or be transported across the cell membrane.For the ease of analysis and simulation, we also consider reducing the2D models to 1D by considering some suitable subsets of the cell. The firstway is to take a section across a diameter of the cell. This results in no-fluxboundary conditions for the 1D PDEs at the two ends of the cell. The secondway is to take the perimeter of the cell as our 1D domain. This give riseto periodic boundary conditions. Numerical experiments confirm that thesetwo approaches give mostly analogous results. The subsequent sections will82.2. Wave pinning modelonly use the first approach.Figure 2.1: Diagram illustrating the two ways to reduce the model to 2D.Left: horizontal slice across a diameter, with L being the length of thediameter. Right: the perimeter, with periodic parameterization s.2.2 Wave pinning modelThis model was first proposed by [30] and involves only one GTPase inisolation. It is the basic model that the subsequent ones below built upon.It is intended to capture of polarization behavior necessary for cell motility.Let u(x, t), v(x, t) represent the concentration of active and inactive formof a generic GTPases, respectively. Both quantities diffuse, but the inactiveform diffuses much faster since GTPase can only be active when it is boundto the membrane, whereas the inactive form can move more freely insidethe cytosol. Suppose that there is a basal activation rate, k for v to beconverted to u. Furthermore, u has a positive feedback effect on itself,which can be represented as a Hill-function activation term with height γand Hill coefficient n (which controls its “sharpness”). The choice of theactivation function is not important, as [31] showed. Any function withroughly the same shape will do. Since the activation is already nonlinear,we can assume deactivation is linear and still obtain interesting results. Thisyields the following equations:∂u∂t= Du∇2u+ f(u, v), (2.1a)∂v∂t= Dv∇2v − f(u, v), (2.1b)f(u, v) = (k + γunun0 + un)︸ ︷︷ ︸A(u), activation ratev − η︸︷︷︸inactivation rateu, (2.1c)92.2. Wave pinning modelwhere Du Dv, t > 0, x ∈ Ω, where Ω is a 1D interval [0, L], or a 2Dshape for example a square [0, L]× [0, L] or a circle. The no-flux boundaryconditions are:∂u∂x=∂v∂x= 0 ∀ x ∈ ∂Ω. (2.2)In 1D domain, this can be written as∂u∂x∣∣∣∣x=0=∂v∂x∣∣∣∣x=0=∂u∂x∣∣∣∣x=L=∂v∂x∣∣∣∣x=L= 0 .The variables and the parameters have the following units:[u] = [v] = [u0] =Glm, [Du] = [Dv] =l2t,[k] = [γ] = [η] = t−1, [L] = l, [n] = 1,(2.3)where m = 1 or 2 is the dimension of the domain, G refers to the concen-tration of GTPase, l is length and t is time. These are all strictly positive.Moreover, n ≥ 2 is an integer.2.2.1 Non-dimensionalizationThis model can be scaled in order to reduce the number of parameters. Letx˜ =xL,(y˜ =yLif 2D), u˜ =uu0, v˜ =vu0,t˜ =tL2/Dv=tT, T = L2/Dv.Applying this to the system (2.1) yields:∂u˜∂t˜=DuDv∇˜2u˜+ Tu0f,∂v˜∂t˜= ∇˜2v˜ − Tu0f,Tu0f =Tu0(k + γu˜21 + u˜2)v˜u0 − Tu0ηu˜u0,where ∇˜2 = ∂2∂x˜2+∂2∂y˜2.Now defineδ =DuDv, f˜ =Tu0f, k˜ = Tk, γ˜ = Tγ, η˜ = Tη,102.3. Actin feedback extensionand finally drop all tilde to obtain the scaled form of the system:∂u∂t= δ∇2u+ f(u, v), δ 1 (2.5a)∂v∂t= ∇2v − f(u, v), (2.5b)f(u, v) = (k + γun1 + un)v − ηu . (2.5c)The remaining variables (u, v), quantities (x, t) and parameters (δ, k, γ, η)are all non-dimensional. The domain has been scaled to [0, 1] for 1D and[0, 1]× [0, 1] for 2D.The characteristic behavior of this system when perturbed from ahomogeneous equilibrium is the formation of a travelling wave-front thatslows down and stalls before reaching the boundary, resulting in a polarizedstate, hence the name wave pinning. In the context of cells, this can beinterpreted as cell motility, where one end of the cell expand and the otherend contracts. For simulations of this model, see Sec. 3.3.2.3 Actin feedback extensionThis extension was proposed by [16] in order to capture the well-known actinwave effect observed in cells. Examples of such observations can be foundin [44], for slime mold cells, and in [46], for human white blood cells. Theseexperiments found trains of waves of actin activity propagating throughoutthe cells.It is known that Rho GTPases such as Rac and Cdc42 up-regulate actinpolymerization via a complex network of interacting proteins, including theWAVE complex, WASP and Arp2/3 [41], and in return, F-actin (the poly-mer form of actin) down-regulates Rho GTPase through RGA-3/4, for ex-ample[39]. This creates a negative feedback effect, making it natural toexpect patterns periodic in time.In the spirit of seeking simplicity, the entire network of interactions canbe abstractly represented as follows: F-actin down-regulates GTPase activa-tion while GTPase up-regulates F-actin. This leads to the actin wave model,112.4. Non-conservative extensionwhere F represents F-actin:∂u∂t= δ∇2u+ f(u, v, F ), (2.6a)∂v∂t= ∇2v − f(u, v, F ), (2.6b)∂F∂t= (knu− ksF ), (2.6c)f(u, v, F ) = (k + γun1 + un)v − (η + s F1 + F)u . (2.6d)The above is the scaled version of the system. The scaling procedure issimilar to Sec. 2.2.1. Notice that if s = 0, then F decouples from thedynamics of u, v and the wave pinning system is recovered.  is small, whichrepresents the relatively slow response of actin to GTPase. This model canproduce a variety of behaviors, shown in Sec. 3.5. These behaviors includewave pinning with oscillating front, reflecting pulses and travelling wavetrains, the last of which can be interpreted to correspond to actin wavesobserved in experiments.2.4 Non-conservative extensionIn [43], the authors proposed to add source and sink terms to the GTPasedynamics. The source term represents the production of inactive GTPase.The sink term was originally intended to represent dilution effect due toGTPase causing the cell expanding locally, which is applicable only to Racand Cdc42. Alternatively, it can be interpreted as the natural decay of activeGTPase. The goal was to capture formation of localized active GTPase, andthe idea to add source and sink terms comes from [36]. The scaled systemis:∂u∂t= δ∇2u+ f(u, v)− cθu, (2.7a)∂v∂t= ∇2v − f(u, v) + cα, (2.7b)f(u, v, F ) = (k + γun1 + un)v − ηu . (2.7c)Notice that setting c = 0 recovers the wave pinning model. This modeltends to produce a spatially periodic pattern consisting of a series of evenly-spaced, static spikes. Numerical results can be found in Sec. 3.6. Thesespikes, representing highly localized GTPase activity, can be interpreted asfilopodia formation in the case where the GTPase causes the cell to expand.122.5. Combining the models2.5 Combining the modelsIt is straight-forward to combine these two extensions to form a combinedmodel:∂u∂t= δ∇2u+ f(u, v, F )− cθu, (2.8a)∂v∂t= ∇2v − f(u, v, F ) + cα, (2.8b)∂F∂t= (knu− ksF ), (2.8c)f(u, v, F ) = (k + γun1 + un)v − (η + s F1 + F)u . (2.8d)As before, all parameters are strictly positive (except c, s ≥ 0, and n ≥ 2)and have unit dimension. c and s serve as switches to turn the two extensionson and off.Notice that the Hill function in f is monotone increasing and boundedbetween 0 and 1, so0 < k < ∂vf = k + γun1 + un< k + γ . (2.9)This holds regardless of the other parameters, especially c and s, whichmeans it applies to all models. This observation will be useful later.This model exhibits complex behaviors that can be seen as a mixture ofthe behaviors under the two extensions, unsurprisingly. Similar to the non-conservative extension, usually a series of spikes form. However these spikesdo not remain static. Instead they move around with velocity dependingof s, and they can collide and interact in interesting ways. The simulationresults can be found in Sec. 3.7.2.6 Mass conservation and well-mixed modelsAn important feature to notice is that when c = 0, the total mass of GTPaseis conserved. This can be shown with the standard argument for no-flux132.6. Mass conservation and well-mixed models(a) (b)(c) (d)Figure 2.2: Graphical representation of the four models. Figure producedby Leah Edelstein-Keshet. (a) The wave pinning model. (b) The actin wavemodel. (c) The source & sink model. (d) The combined model. In thesediagrams, u and v represents active/inactive GTPase as in the equations,F-actin corresponds to F , the synthesis and decay arrows represents the cαand cθu terms respectively, and the other arrows indicate feedback effects.boundary conditions and the divergence theorem:∂∂t∫Ω(u+ v)dx =∫Ω(δ∇2u+ f +∇2v − f)dx (substitute in the PDE)= δ∫Ω∇2udx+∫Ω∇2vdx= δ∫∂Ω∇u · nˆ dS︸ ︷︷ ︸=0 due to no-flux BC+∫∂Ω∇v · nˆ dS︸ ︷︷ ︸=0= 0 .(2.10)Consider well-mixed systems, that is (2.8) with homogeneous initial condi-tion. Then the diffusion terms drop out and we are left with a system of142.6. Mass conservation and well-mixed modelsODEs:ut = f(u, v, F )− cθu, (2.11a)vt = −f(u, v, F ) + cα, (2.11b)Ft = (knu− ksF ) . (2.11c)Equilibria of this system are called homogeneous steady states (HSS). Ifc = 0, by mass conservation w = u+ v, the mean concentration of GTPaseacross the domain, is constant. Therefore, we can reduce the system tout = f(u,w − u, F ), (2.12a)Ft = (knu− ksF ) . (2.12b)Furthermore, if s = 0, then we are left with a single ODE:ut = f(u,w − u) = g(u) . (2.13)In this chapter, we introduced four models for GTPase dynamics andshowed simplifications possible for the well-mixed models under mass con-servation. In the next chapter we will demonstrate their behavior usingnumerical simulations.15Chapter 3Numerical simulationsThis chapter contains numerical solutions of the models introduced in Chap-ter 2. The numerical schemes used to produce these simulations are dis-cussed in Sec. 3.1, and considerations for the parameters are given in Sec. 3.2.The simulations are run for long enough such that the system reaches an(apparent) steady state, or until the pattern is obvious.The numerical simulation results of [30, 16, 43] are reproduced as a partof Sec. 3.3, 3.5, 3.6. These results are compared with simulations that usesdifferent initial conditions and parameters that either better illustrate themodel behaviors or lead to distinct behaviors not discussed by the earlierpapers.This chapter is closely related to the next two chapters. Much of theanalysis is guided by the properties of the numerical solutions, while theanalysis identifies possibly distinct parameter regimes and calls for moresimulations.3.1 Numerical methodsFor the 1 spatial dimension case, we use a finite difference scheme in thespatial domain with 400 mesh points based on a modified version of Crank-Nicolson (CN) [4]. CN is an implicit scheme with 2nd order accuracy andreasonably efficient for linear systems (well known, see [22, Ch. 7.2]). How-ever, the nonlinear reaction term f can be computationally expensive forimplicit treatment. The remedy is the IMEX (implicit-explicit) methodproposed by [40]. The idea is to treat the nonlinear reaction term explicitlyto save on complexity, while still treating the linear diffusion term implicitly,resulting in a linear system. It was found empirically that this adaptationmostly retains the stability of CN. Therefore, we use the IMEX-CN schemeas a compromise between stability and efficiency.When actin feedback is present in the model, we use the simple explicitEuler scheme for F since F does not diffuse. This scheme can be summarized163.2. Selection of simulation parametersas follows:(I − 12hδ∆˜)un+1 = un + h[f(un, vn, Fn)− cθun + 12δ∆˜un],(I − 12h∆˜)vn+1 = vn + h[−f(un, vn, Fn) + cα+ 12∆˜vn],Fn+1 = Fn + h(knun − ksFn) .where h = ∆t is the time step and ∆˜ is the discrete Laplacian. This isimplemented in Matlab, and the program is provided in Appendix A.For 2D simulations, we use a finite element method, with method of linesfor spatial-temporal discretization and 1st degree polynomial interpolationfor basis functions. The spatial domain is discretized with a regular mesh of5000 triangles, and we use the same time stepping scheme (IMEX-CN) as inthe 1D case. This is implemented with the FEniCS [2] package in Python,and the code is provided in Appendix C.In both cases, we verified that the chosen step size h is sufficiently smallby running the same simulation with time step h and h/2, producing solu-tions uh(x, t) and uh/2(x, t) (and also for v and F ), and checking that theydo not differ by more than 0.01 at any x ∈ Ω, t > 0.3.2 Selection of simulation parametersFor all simulations, we refer to the combined system of equations (2.8).The default parameters (which produce the most typical behavior for thesystems) are summarized in Table 3.1.Notice in the table that the unscaled domain length L is re-introducedas a parameter, whereas the non-dimensional domain remains 0 ≤ x ≤ 1.Explicitly varying L gives a convenient way to change multiple parameters.The selection of value for most of these parameters is based on [16], withα, θ coming from [43], and some modifications guided by Turing analysis(Ch. 4) and local perturbation analysis (Ch. 5).Although these sets of parameters produce many known behaviors, itmight not yield clean bifurcation diagrams for analysis. Therefore in otherchapters, we use different parameter sets. Additionally, we will deviate fromthe default set in order to discuss the effect of certain parameters.One consideration in choosing between parameter sets that produce thesame qualitative behavior is their robustness, which means we desire to havethe solution reproducible on a wide range of parameters. For example, theauthors of [16] choose n = 3 for the actin wave model because produces a173.3. Simulations of wave pinning in 1DParameter Meaning WP NC AF Combined2δ Diffusion coef. ratio 0.01L Domain length 1 10k Basal activation rate 1.5L2 1L2 − 6L2 1L2γ Nonlinear activation rate 30L2n Hill coefficient 3 2η Inactivation rate 15L2 5L2 15L2 5.2L2c NC terms on/off 0 1 0 1α Source strength 1.5L2θ Sink strength 4.5L2 5.5L2s Actin feedback strength 0 0− 50L2 Actin reaction rate 0.1kn Actin activation rate 24L2ks Actin inactivation rate 7.5L2Table 3.1: The parameters in the combined model, their meaning, and thevalues they take for different simulations. WP: wave pinning (Sec. 3.3);NC: non-conservative extension (Sec. 3.6); AF: actin feed back extension(Sec. 3.5); Combined 2: one of the parameter set used for the com-bined model (Sec. 3.7). All parameters (except L) are scaled to be non-dimensional, as in Sec. 2.2.1.considerably larger pattern-forming regime in the parameter plane comparedto n = 2, which was used by [30, 43]. In Sec. 3.4, we show why a sharpernon-linear function makes wave pinning more robust. For this reason, wechoose n = 3 as the default parameter for most of simulations. However, inthe LPA chapter (Ch. 5), we often use different parameters that result inbifurcation diagrams that are easier to read.3.3 Simulations of wave pinning in 1DFor wave pinning, the parameters are taken from Table 3.1(WP). Theseparameters satisfy the condition that δ is small, and two stable homogeneoussteady states exist when v = 1 as discussed in Sec. 3.4.For initial conditions, the default is to start at v = 1 and u ≈ 0.102,which corresponds to the lower homogeneous equilibrium. The higher equi-librium is at u ≈ 1.812. Various perturbations will be applied to u initiallyto trigger pattern formation. Another possibility is to start at v = 1, u =183.4. Simple analysis of wave pinning0.834·(x), where u = 0.834 is the unstable equilibrium and (x) ∼ Unif [0, 2]is random noise.Fig. 3.1 shows the simulation results in the form of kymographs, wherethe horizontal axis is time, the vertical axis is space, and color indicates thevalue of u, v. Observe that for the first two cases, the initial perturbationeventually smoothens out. However it does not cover the entire domain,resulting in polarization. The random initial condition also results in apolarized steady state. In these simulations, u can vary greatly across thedomain while v is much closer to uniform, as expected given the differencein diffusion rate.In Fig. 3.2, we use similar initial conditions except that the perturbationat one end is negative. The result is similar to the previous case except thatthe role of high/low u is flipped.It is important to mention that if we start at the steady homogeneousequilibrium, and the perturbation is too narrow or too small, then the systemwill eventually be able to reach homogeneity again and polarization will notoccur. With random initial conditions, polarization is also not guaranteed.The proportion of the domain with high/low active GTPase depends on thetotal amount of GTPase, and the polarity (i.e. whether u is high on the leftor right side) is determined by the initial condition.3.4 Simple analysis of wave pinningIn the wave pinning model, the steady state of the typical numerical solution(Fig. 3.3) is a partition of the domain into three regions, a high-u, a low-uand a relatively narrow transition zone, with u nearly constant in the firsttwo. Due to the fast diffusion of v, its level is nearly constant throughout thedomain. This suggests that in both the left and right regions, the solutionis approximately a stable homogeneous steady state (HSS). This turns outto be the case for Fig. 3.3. However, note that the final v value is notnecessarily close to its initial value. Numerical experiments show that thissteady state is stable to small inhemogeneous perturbations. We call thisfinal steady state the “wave pinning solution”.These observations suggests the following necessary condition for wavepinning to occur: for some value of v∗, there are two HSS (u−, v∗), (u+, v∗)with u− < u+. Furthermore, they must be stable in two senses: (1) stableto any uniform perturbation in u (while v is held constant), and (2) stable toperturbation to both u and v, while conserving mass. These can be thoughtof as a common stable equilibrium to the following two ODEs, where u = u±193.4. Simple analysis of wave pinning(a) u (b) v(c) u (d) v(e) u (f) vFigure 3.1: Numerical simulation of the wave pinning model (2.5), withparameters from Table 3.1(WP). Initial condition: v = 1, u = 0.102 on thewhole domain except u = 6 for (a,b) 0 ≤ x ≤ 0.1; (c,d) 0.4 ≤ x ≤ 0.5; (e,f)random noise, u = 0.834 · (x). Note that not all initial conditions resultin wave pinning: a small perturbation from the HSS will simply decay andno pattern forms. The behaviors shown in (a,b,e,f) reproduces the solutionsshown in Fig. 2 of [30]. 203.4. Simple analysis of wave pinning(a) u (b) vFigure 3.2: Simulation of the wave pinning model (2.5) with parameters fromTable 3.1(WP), except γ = 150, and initial condition v = 0.1866, u = 1.3134except u = 0 on 0 ≤ x ≤ 0.1. In this case “negative wave pinning” occurs:the negative perturbation in u leads to formation of a low u region on theleft of the domain.satisfy:∂u∂t= f(u, v∗), where v∗ is a constant parameter, (3.1)∂u∂t= f(u,w − u), where w = u+ v∗ is a constant parameter, (3.2)This can be stated equivalently as a condition on f as follows: there existv∗, u− < u+ such thatf(u−, v∗) = f(u+, v∗) = 0, (3.3)∂f(u±, v∗)∂u< 0, (3.4)df(u,w − u)du=∂f(u±, v∗)∂u− ∂f(u±, v∗)∂v< 0 . (3.5)However, we know from (2.9) that ∂f∂v > 0. Therefore the condition (3.5)is redundant given (3.4). Related work in [31] considers a more generalfamily of functions f of different forms, where both of these conditions arenecessary. The authors found that these two conditions, together with the“Maxwell condition”, is sufficient for wave pinning to occur.The equilibria (u±, v∗) can be visualized by plotting f , which is done213.5. Simulations of actin feedback in 1DFigure 3.3: Final steady state of the wave pinning model showing threedistinct regions: a high-u region on the interval 0 ≤ x ≤ 0.2 where(u, v) ≈ (1.7506, 0.9801), a low-u region on 0.5 ≤ x ≤ 1 with (u, v) ≈(0.1018, 0.9664), and a transition zone in between. In all regions v ≈ 1is nearly constant, and the extremes are close to the high (1.812, 1) andlow (0.102, 1) equilibria. Parameters and initial conditions same as inFig. 3.1(a,b).in Fig. 5.4, and this analysis is further expanded with LPA techniques inSec. 5.2.3.5 Simulations of actin feedback in 1DWe use the parameters from Table 3.1(AF) to simulate the actin feedbackmodel (2.6). In the rest of this section, the L2 factor is implied (e.g. k = 6really means k = 6L2).The default initial conditions are u = 0 except u = 4 for 0 ≤ x ≤ 0.01,v = 2.5, F = 0. We are not initializing near the stable HSS because doing sousually results in no patterning. The patterns observed are quite sensitiveto initial conditions.A variety of behaviors can be observed. Fig. 3.4 displays four qual-itatively different behaviors: (1) wave pinning with oscillating boundary(WPO), where polarization occurs similar to wave pinning, but the frontdividing the polarized regions oscillates back and forth; (2) reflecting pulse(RW), where a single pulse traverses the domain at constant velocity andgets reflected back whenever it reaches the boundary; (3) a single pulse (SP)223.6. Simulations of non-conservative model in 1Dwhich is absorbed at the boundary, the system reaches an HSS afterward; (4)wave train (WT), a train of waves that originates either at a boundary or inthe interior of the domain, and propagates with constant velocity until thewaves are absorbed at boundaries. Additionally, with k = 1.5, s = 5, we canobserve simple wave pinning (WP; not shown, visually similar to Fig. 3.1).Out of these, WT is the closest behavior to the actin waves observed inbiological experiments [44, 46]. WPO appears to be related to simple wavepinning via a Hopf bifurcation for the location of the wave front. This ishard to verify since we do not have a closed-form equation that describesthe movement of the wave front. The WPO behavior is qualitatively similarto the “tango waves” found in a very different reaction-diffusion model by[23]. The relation between these two behaviors is unknown.In general, the spatial profile of F lags behind u as expected since itis a slowly reacting variable depending on u. The pattern in v is usuallyin-phase with u but with opposite sign, i.e, v is high where u is low, andvice versa. Moreover, the gradient of v tends to be much shallower than udue to its fast diffusion.Additionally, there are some more complex patterns possible, shown inFig. 3.5. These two share certain characteristics with the simpler patterns.The patterns shown in (a-c) are quite similar to WPO, except that thedomain is divided into five polarized regions instead of two. The transientbehavior near the beginning is also reminiscent of WT. The patterns in (d-f)can be seen as RW but with a group of four pulses rather than one.Compared to the results in [16], these pictures show a richer range ofpatterns using mostly the same parameter set (with different scaling). Themain difference is that here we used a longer domain L = 10, providing morespace for pattern to develop, whereas [16] used L = 1, which means thepatterns are more confined and boundary effects are much more prominent.3.6 Simulations of non-conservative model in 1DFor the non-conservative extension (2.7), we use the parameters from Ta-ble 3.1(NC). One noticeable difference with the parameter set (WP) is themuch larger domain size L = 10. Since this model tends to produce pe-riodic patterns, L must be long enough to be able to support these pat-terns. The values for α, θ are chosen so that δ = 0.01, γ = 30 puts thesystem in the Turing regime (the pattern-forming regime where the HSS islinearly unstable, see Sec.4.2.2), yet is not too far from the regime’s bound-ary (see Fig. 4.2). Under this set of parameters, the unique equilibrium233.7. Simulations of the combined model in 1D(u∗, v∗) ≈ (0.3333, 1.2315) (see (4.13)).Fig. 3.6 shows some simulation results. In (a,b), we initialize at the HSSand perturb u on the left end of the domain. This perturbation quicklyevolves into a narrow spike. The pattern then propagates throughout thedomain, forming more spikes at regular intervals until running out of space.After that, the spikes shifts slightly to be evenly spaced and eventually reachsteady state shown in Fig. 3.7(a). This is the typical pattern formed underthis model. Notice that the transition between the initial HSS and the finalpattern happened very rapidly.We can also start with small random noise superimposed on the HSS:u = u∗(1 + Unif(−0.01,+0.01)). In this case, after a brief transient phase,periodic patterns similar to Fig. 3.7(a) form. The exact number of peakscan vary depending on initial condition. We will discuss the wave length ofthe pattern in Sec. 4.2.3.We have also experimented with slightly different parameter sets. InFig. 3.6(c,d), we use a set of parameters where the HSS is linearly stable(as found in Sec. 4.2.2), so small perturbations to the HSS simply decays.However, this does not mean that the HSS is globally stable, since a strong,localized perturbation can lead to the formation of a single persistent spike,shown in Fig. 3.6(f). We call this the soliton solution. Its formation will beexplained in Sec. 5.3.Most other variations of the parameter set lead to patterns similar toFig. 3.6(e), except for two interesting cases. If we use a smaller domainlength, for example L = 1, we only get wave pinning, shown in Fig. 3.7(a).This is expected since the domain is likely too small to support a full periodof the pattern. Using η = 15 leads to a periodic series of spikes with the topof the peak bifurcates into two, as shown in Fig. 3.7(b).3.7 Simulations of the combined model in 1DWith both extensions “on”, the system displays much more complex behav-iors. The first set of parameters we use are the same as Table 3.1(NC), butwith s varying. For initial conditions, we start at the unique HSS:u∗ =αθ, F∗ =knksu∗, v∗ =cα+ (η + sF∗1+F∗ )u∗k + γ un∗1+un∗, (3.6)and perturb it with some noise: u = u∗(1 +Unif(−0.1,+0.1)). The resultsare shown in Fig. 3.8. With s low enough, the system behavior is similar to243.8. Simulations of wave pinning and non-conservative extension in 2Dthe s = 0 case, where we get a static, spatially periodic pattern characteris-tic to the non-conservative extension (Fig. 3.8(a), compare to Fig. 3.6(a)).As we increase s, the peaks begin to move, reminiscent of the wave trainbehavior seen under the actin wave extension (Fig. 3.4(j)). These peaksmove with constant velocity by themselves and appears to repel each otherwhen the peaks get close. When s is moderate, this repulsion is strongenough to make the peaks reverse their directions if they are on collisioncourse (Fig. 3.8(b)). At higher s, they collide. With the parameter set fromTable 3.1(NC), both pulses tends to survive when they collide and they sim-ply pass through one another (Fig. 3.8(c)). However, for L = 20, we alsoobserve the formation of a standing pulse that remains stationary for sometime with a persistent oscillating amplitude (Fig. 3.8(d)).We also use another parameter set from Table 3.1(Combined2). Underthese parameters, the peaks sometimes annihilate upon collision (Fig. 3.9(c)),and sometimes only one survives (Fig. 3.9(b)). At even higher s, we observea localized standing wave pattern that oscillates rapidly (Fig. 3.9(d)).3.8 Simulations of wave pinning andnon-conservative extension in 2DWe ran simulations in 2D for the wave pinning and non-conservative modelsusing the same parameters as before. The initial conditions are u = u∗except u = 4u∗ in the corner [0, 0.1]× [0, 0.1], and v = v∗, where (u∗, v∗) =(0.3333, 1.1111) is the HSS for the non-conservative model. The results areshown in Fig. 3.10. The results for wave pinning are not surprising: uinitially spreads out from the corner as a 2D travelling wave, and that iseventually pinned along a front determined by the initial conditions.The non-conservative model is more interesting. Based on 1D simula-tions, we expect evenly-spaced stripes to form around the corner as concen-tric rings, as happens initially. However, these rings quickly break up intospots. The spots shift around to become more evenly spaced, and then theysettle as a steady state is reached. We have not found any parameter setsunder which the ring pattern is stable.These patterns are not sensitive to the shape of the domain. Simulationson circular, rectangular and other domains with simple shapes producedpatterns with the same qualitative characteristics (not shown).253.9. Simulations with spatially inhomogeneous parameters3.9 Simulations with spatially inhomogeneousparametersIn this section we consider cases where the parameters are dependent onspace. The main motivation is that locally up-regulated GTPases are re-quired for certain kinds of viruses to infect cells. Specifically, [48] shows thatthe formation of filopodia, which is regulated by GTPases, is essential fordengue virus to invade human endothelial cells. Other papers have shownthat filopodia play a role in the infection of rabies [47] and HIV [32].The main approach is to start with the combined model (2.8) in 1Dand identify a parameter set that produces behaviors that we interpret asfilopodia formation. We then allow one of the parameters to vary smoothlyin space, such that there are regions where the parameter lies in the pattern-forming regime where the HSS is unstable, and others where it lies in thestable regime. The boundaries for the two regime can be estimated byrunning simulations with constant parameters, as in Sec.3.7. We expect torecover the familiar patterns in the unstable region. It would be interestingto see if the patterns can persist into the stable region, that is, whether wecan get coexistence of the HSS and some stable pattern.For this, we use the parameters from Table 3.1(Combined2), except withL = 20 so that the domain is long enough to support more peaks, whichmakes it easier to see the behavior in the interior of each region. There aretwo parameters we can change to represent local up-regulation of GTPase:the basal activation rate k, and the feedback activation parameter γ. Wepresent the results usingγ(x) =[15 + 15 exp(−(x− 1/2)2/σ2))]L2, (3.7)which is a Gaussian bump centered at the middle of the domain. Makingk spatially dependent instead yields analogous results. From numerical ex-periments, with constant γ, the boundary between the stable and unstableregime is near γ ≈ 18.5. The unstable region where γ(x) > 18.5 is aroundhalf of the domain when σ = 0.2.The results in Fig. 3.11 show that the pattern does not leak into thestable region. Furthermore, if σ is small enough, i.e. the unstable region isnarrow enough, there will be a lone peak confined within it which oscillatesin-place. In general, the system behaves as if there is a boundary nearγ ≈ 18.5 with Dirichlet boundary conditions. This suggests that local up-regulation does not impact the system behavior globally. Relating back tobiology, this means that viral activity can create localized spots of high263.9. Simulations with spatially inhomogeneous parametersGTPase activity. For Rac GTPase, this will likely cause a protrusion thateventually evolves into a filopodia, facilitating viral entrance into the cell.273.9. Simulations with spatially inhomogeneous parameters(a) u, k = 1.5, s = 18 (b) v, k = 1.5, s = 18 (c) F , k = 1.5, s = 18(d) u, k = 1.5, s = 27 (e) v, k = 1.5, s = 27 (f) F , k = 1.5, s = 27(g) u, k = 1.5, s = 36 (h) v, k = 1.5, s = 36 (i) F , k = 1.5, s = 36(j) u, k = 6, s = 30 (k) v, k = 6, s = 30 (l) F , k = 6, s = 30Figure 3.4: Numerical simulations of the actin feedback model (2.6) withparameters from Table 3.1(AF), with s, k as indicated by the labels, and de-fault initial conditions. Each row corresponds to one parameters set, showingu, v, F (left to right). We can observe the four behaviors discussed in thetext as we vary k and s: (a-c) Wave pinning with oscillating front (WPO);(d-f) Reflecting waves (RW); (g-i) Single pulse absorbed at boundary (SP);(j-l) Persistent wave trains (WT). These figures used a longer domain lengthcompared to [16], which gives more room for the patterns to develop andlead to richer behaviors compared to Fig. 3 of [16]. 283.9. Simulations with spatially inhomogeneous parameters(a) u (b) v (c) F(d) u (e) v (f) FFigure 3.5: Some of the more exotic patterns observed under the actinfeedback model (2.6) with parameters from Table 3.1(AF) and varying k, s.(a-c) k = 5, s = 10, default initial conditions. The pattern is similar to WPObut with more than two polarized regions; (d-f) k = 5, s = 30, default initialconditions except that the excitation region is 0 ≤ x ≤ 0.1. The resultingpattern is similar to RW except that we have a group of four pulses traversingthe domain together. These patterns were not discovered by [16], althoughthe authors did find some “exotic” patterns not shown here.293.9. Simulations with spatially inhomogeneous parameters(a) u (b) v(c) u (d) vFigure 3.6: Numerical simulation of the non-conservative model (2.7) with(a,b) default parameters (Table 3.1(NC)); (c,d) γ = 15L2, η = 15L2. Initialcondition: (a,b) u = u∗ except u = 1 on 0 ≤ x ≤ 0.1, v = v∗. (c,d)u = u∗ = 0.33333 except u = 10u∗ on 0.4 ≤ x ≤ 0.41, v = v∗ = 3.19298.In (a,b), the formation of a peak on the left side of the domain triggers theformation of new peaks some regular distance away, until space runs out.Once all peaks are formed, they shift slightly so that they are evenly spaced.In (c,d), the single initial peak persists, but it does not lead to formation ofnew peaks. We call this the soliton solution.303.9. Simulations with spatially inhomogeneous parameters(a) Default parameters (b) γ = 15L2(c) L = 1 (d) η = 15Figure 3.7: Final steady state pattern of the non-conservative model (2.7)with most parameters from Table 3.1(NC), except the parameters indicatedon the labels. (a) and (b) correspond to the steady state of Fig. 3.6(a,b)and (c,d) respectively. In (c) the shortened domain results in wave pinning;(d) Higher deactivation rate η = 15 results in peaks with bifurcated tops.(a,b) reproduces Fig. 5 (a,d) of [43], respectively.313.9. Simulations with spatially inhomogeneous parameters(a) L = 10, s = 3 (b) L = 10, s = 8(c) L = 10, s = 18 (d) L = 20, s = 35Figure 3.8: Numerical simulations of the combined model (2.8) with defaultparameters from Table 3.1(NC) except with varying L, s, and HSS+noiseinitial condition as described in the text. Here we only show the kymographfor u, since both v and F are similar except that F lags a little behind, andv is out-of-phase. Their correspondence can be seen in earlier simulations.Notice the progression from static peaks, to moving but repelling peaks, tocolliding peaks as s increases.323.9. Simulations with spatially inhomogeneous parameters(a) s = 5 (b) s = 8(c) s = 18 (d) s = 35Figure 3.9: Numerical simulations of the combined model (2.8) with pa-rameters from Table 3.1(Combined2) with varying s, and HSS+noise initialcondition as described in text. Notice as we increase the actin feedbackstrength s, the behavior transitions from slowly moving, repelling peaksto colliding peaks, similar to Fig. 3.8. At higher s, we observe a rapidlyoscillating standing wave pattern in some parts of the domain.333.9. Simulations with spatially inhomogeneous parameters(a) Wave pinning, t = 1.25(b) Wave pinning, t = 30(c) Non-conservative extension, t = 0.145(d) Non-conservative extension, t = 0.994Figure 3.10: 2D simulations of the wave pinning (a,b) and non-conservative(c,d) model, using the same parameters as in 1D (Table 3.1(WP) and (NC)),with u shown on the left and v on the right. For each model, two snap shotsare shown: one when the pattern begin to take shape, and another after thesystem reached steady state. 343.9. Simulations with spatially inhomogeneous parameters(a) s = 0, σ = 0.2 (b) s = 0, σ = 0.2(c) s = 10, σ = 0.2 (d) s = 10, σ = 0.2(e) s = 10, σ = 0.03 (f) s = 10, σ = 0.03Figure 3.11: Simulations for the combined model (2.8) with L = 20, andspatially-dependent parameter γ(x) as in (3.7), s, σ as indicated in the labels,and all other parameters from Table 3.1(Combined2). The kymographs ofu, and the snapshots at t = 1 are shown. Notice that the pattern formsin the middle of the domain where γ is high enough to be in the unstableregime, and this pattern does not leak outside to the region where γ is small.This reflects localized filopodia formation in response to viral activity.35Chapter 4Turing analysisIn this chapter, we use Turing analysis to determine the stability of thehomogeneous steady state (HSS) solutions under the wave pinning and non-conservative model. If a HSS is Turing-unstable, then even an infinitesimalperturbation can grow and lead to pattern formation. The results in thissection will help to explain some of the patterns seen in Ch. 3, but notall of them. Other patterns, which require the initial perturbation to besufficiently strong, will be explained by another technique (LPA) in Ch. 5.In [43], the authors performed weakly nonlinear analysis on the non-conservative model, which is a more general method that contains Turinganalysis as its leading order. We show the correspondence between ourresults and [43] in Sec. 4.2.2.4.1 Overview of Turing analysisTuring analysis is a technique for analyzing the stability of HSS solutions forreaction-diffusion systems. This was first done by Turing in [42], which wasthe earliest mathematical work to study diffusion-driven pattern formation.A more modern treatment of the theory can be found in [8, Chapter 11].The main idea of Turing analysis is to linearize the system around a HSS.This HSS is required to be stable in the well-mixed model (WM), so thatdiffusion is necessary for instability. Then, use an ansatz for the solution todetermine the stability of all possible spatially periodic modes. The HSS islinearly stable if and only if all modes decay in time.We will focus on wave pinning and the non-conservative extension inone spatial dimension, which is within the scope of traditional analysis.The results show that Turing regimes exist for both of these cases. Thesimpler techniques no longer apply when actin feedback is added. Althoughthe method can be generalized to those cases, the results will be harder tointerpret.364.2. Turing analysis of wave pinning model and non-conservative extension in 1D4.2 Turing analysis of wave pinning model andnon-conservative extension in 1DFor the wave pinning or non-conservative model in 1D, we start by linearizingthe system (2.7) around a chosen HSS (u∗, v∗). Define(u˜, v˜) = (u, v)− (u∗, v∗),the linearization can be written as∂∂t[u˜v˜]= D∂2∂x2[u˜v˜]+ J[u˜v˜], (4.1)where D is a matrix containing the diffusion coefficients on the diagonal,and J is the Jacobian of the well-mixed (WM) system at (u∗, v∗):D =[δ 00 1], J =[∂f∂u − cθ ∂f∂v−∂f∂u −∂f∂v]∣∣∣∣∣(u,v)=(u∗,v∗). (4.2)We are especially concerned with instabilities that occur only in the presenceof diffusion. Since ∂vf > 0 (2.9), we havedet(J) = −(∂uf − cθ)∂vf + ∂uf∂vf = cθ∂vf ≥ 0 . (4.3)For c > 0, this means the HSS is stable in WM if and only if Tr J < 0.In the case c = 0, the system is degenerate. Using mass conservation, thisdegeneracy can be removed (as in (2.13)) and the Jacobian of the equivalent1D system is exactly Tr J , so the same condition applies.Next, consider small perturbation from the equilibrium, i.e. |u˜|, |v˜| 1,therefore neglecting all non-linear terms. We are interested in whether thisperturbation grows or not. The normal form ansatz is:[u˜v˜]=[αuαv]cos(qx)eσt .Here q > 0 is the wave number, which is the spatial frequency for the mode.Higher q represents a mode with more peaks in the domain. σ ∈ C is thegrowth rate. A mode will grow if and only if its corresponding <(σ) > 0.Only the cosine modes (and not the sine modes) are allowed due to theno-flux boundary conditions, which requires ∂x(u˜, v˜) = 0 at x = 0.374.2. Turing analysis of wave pinning model and non-conservative extension in 1DInsert the ansatz into the PDE. Neglecting the non-linear terms, anddropping the tildes yieldsσ[uv]= −Dq2[uv]+ J[uv]= (J −Dq2)[uv]= M[uv], (4.4)where M = J −Dq2. Since this results from linearization at an equilibrium,(u, v) ≡ (0, 0) is a trivial solution. For any non-trivial solutions to exist,(4.4) implies that (M − σI) must be singular, i.e. σ is an eigenvalue of M .For this 2-dimensional system, σ can be written asσ =12(TrM ±√Tr(M)2 − 4 detM).Since we require that the HSS (u∗, v∗) be stable in the well-mixed model,we know that Tr J < 0. As a result, TrM = Tr J − q2(1 + δ) ≤ Tr J < 0,which means one of the eigenvalues will always have a negative real part.Therefore, the HSS is linearly unstable if and only if detM < 0.In general, detM is a quadratic function in q2 with a non-negative lead-ing coefficient. For our system, it isdetM = q4δ + q2 (−∂uf + cθ + δ∂vf)︸ ︷︷ ︸A+ cθ∂vf︸ ︷︷ ︸B, (4.5)where the partial derivatives are understood to be evaluated at (u∗, v∗). Theminimum with respect to q2 occurs at q2 = −A2δ . However, when A > 0 theminimum really occurs at q2 = 0, since q2 cannot be negative. Therefore,the minimum value isH = minq2detM =−A24δ+B if A ≤ 0B if A > 0. (4.6)The problem is now is to determine whether H < 0, in which case the HSSis unstable, or H > 0, in which case the HSS is stable.Since ∂vf > 0 (2.9), we have B ≥ 0 always, so if A > 0 then there is noinstability. In summary, the conditions for Turing instability are preciselyTrJ < 0, A < 0, B < A2/4δ. The first two can be expanded as:∂uf − cθ − ∂vf < 0,−∂uf + cθ + δ∂vf < 0.384.2. Turing analysis of wave pinning model and non-conservative extension in 1DTherefore, the conditions can be rewritten as:δ∂vf < ∂uf − cθ < ∂vf, (4.7a)4cθδ∂vf < (−∂uf + cθ + δ∂vf)2 . (4.7b)Observe that (4.7a) implies δ < 1, and it is easier to satisfy both inequalitiesif δ 1. This shows that pattern formation requires δ < 1, which is areasonable assumption based on biology (Sec. 2.2.1).We will show that the conditions in (4.7) are in fact equivalent to theclassical Turing analysis results from [8, Ch. 11.5], which states that thenecessary and sufficient conditions for diffusive instability areJ11 + J22 < 0, (4.8a)J11J22 − J12J21 > 0, (4.8b)J11D22 + J22D11 > 2√D11D22(J11J22 − J12J21)1/2 > 0, (4.8c)where Jij , Dij , i, j = 1, 2 are entries of the matrices J,D given in (4.2).The first inequality (4.8a) is simply Tr J < 0, which reduces to the secondhalf of (4.7a). The second inequality (4.8b) translates to detJ > 0, which isalways true for our model due to (4.3). Finally, squaring the last one (4.8c)and expanding the definitions yields exactly (4.7b). The signs are correctafter squaring due to the first inequality of (4.7a).4.2.1 Wave pinningFor the wave pinning model, c = 0, and (4.5) simplifies toδ∂vf < ∂uf < ∂vf . (4.9)There is always at least one equilibrium that satisfies the second inequality.Consider (2.13), and recallg(u) = f(u,w − u) = A(u)(w − u)− ηu, k < A(u) < k + γ . (4.10)It is easy to see that g(0) > 0 and g(w) < 0, and g is continuous. Hence,there must be a point u∗ ∈ (0, w) where g(0) = 0, g′(0) < 0. We can visualizethis better by writing f(u, v) = A(u)v − ηu, then plot the two terms withrespect to u while holding v constant, as Fig. 4.1. Note that this analysisdoes not require a specific form for A(u). Rather it applies to any sigmoidfunction with A(0) > 0. In this case, A′(u) is a bump-shaped function. Forour choice of using a Hill function, as n increases A(u) becomes sharper394.2. Turing analysis of wave pinning model and non-conservative extension in 1DFigure 4.1: Plot of A(u)v (activation, red, with typical parameter valuesand v constant), and ηu (inactivation, blue, for different values of η). Eachintersection corresponds to a HSS. The HSS is stable in WM at the pointswhere ηu crosses from below A(u) to above. Note that generically, we eitherget a single stable HSS, or a triplet of two stable HSS and one unstable HSSin between.and approaches a Heaviside step function H(u− 1) as n → ∞, while A′(u)becomes more concentrated and approaches a Dirac delta function δ(u− 1)as n→∞. The maximum of A′(u) occurs atumax =(n− 1n+ 1) 1n≈0.5773 if n = 20.7937 if n = 31 as n→∞. (4.11)We can expand (4.9) and with a little manipulation, getδA(u) + η < A′(u)v < A(u) + η. (4.12)Suppose (u∗, v∗) is a stable HSS, that is, it satisfies the second inequality.Then in the case where n is large, we must have u∗ close to umax for thefirst inequality to hold (at least for small δ). This gives a rough method ofpredicting Turing instability based on the location of the HSS. Since HSS arenot unique (in fact we have a one-parameter family of HSS given implicitlyby f(u∗, v)∗) = 0), and each HSS have different stability condition, it is hardto identify the stability boundary for all HSS.404.2. Turing analysis of wave pinning model and non-conservative extension in 1D4.2.2 Non-conservative extensionThe non-conservative extension (2.7) has a unique HSS atu∗ =αθ, v∗ =cα+ ηu∗A(u)=cα+ ηu∗k + γ un∗1+un∗. (4.13)This can be found by realizing that any HSS must satisfy0 = f(u∗, v∗)− cθu∗ = f(u∗, v∗)− cα,and u∗ can be uniquely solved by cancelling the f term, and v∗ solved byre-substituting u∗.This time, since the equilibrium is unique, we can simply evaluateH from(4.6) at the equilibrium to determine its stability given a set of parameters.We identify two main bifurcation parameters of interest: δ, the diffusionratio, and γ, the magnitude of the only non-linear term. Keeping the restof the parameters the same as in Table 3.1(NC), we plot the level curve inδ − γ plane where H = 0, which is the bifurcation boundary separating thelinearly stable and unstable regimes, shown in Fig. 4.2. Compare this figureto Fig. 5 of [43], which is a similar diagram for a slightly different parameterset. Notice that the shape of the curves for θ = 4, 4.5 in Fig. 4.2 and theboundary in Fig. 5 of [43] agree with each other.In the unstable regime, as long as the initial condition is not spatiallyuniform, the system will produce the static, spatially periodic pattern asseen in Fig. 3.6(e). However the number of peaks vary depending on ini-tial conditions and parameters. These patterns correspond to the solutionidentified in Fig. 5(a) of [43].In the stable regime, any small perturbation to the HSS simply decays,the system returns to the HSS and no pattern forms. However, in some casesnear the bifurcation boundary, a large, localized perturbation can lead tothe formation of a single persistent spike. In contrast to the unstable regime,this does not lead to the formation of additional spikes. This corresponds tothe soliton solution seen in Fig. 3.6(f), and the solution identified in Fig. 5(d)in [43]. Since the perturbation is not infinitesimal, this behavior cannot bepredicted by the linear Turing analysis. Later in Sec. 5.3 we will explain thesoliton solution using LPA theory.4.2.3 Predicting the wave length of pattern formed underthe non-conservative modelGiven the periodic steady states seen in the non-conservative model (Fig. 3.6(e)),Turing analysis might be able to predict the period of the pattern. Leave γ414.2. Turing analysis of wave pinning model and non-conservative extension in 1DFigure 4.2: Boundary of the Turing-unstable regime in δ−γ plane, with α =1.5, θ varying and all other parameters from Table 3.1(NC). The unstableregime is on the upper-left, while the stable regime is on the lower-right.as a bifurcation parameter, and fix all others as in Table 3.1(NC) (especiallyθ = 4.5, δ = 0.01 with regards to Fig. 4.2).For this section, it is easier to think about the problem if we rescale xso that the domain is 0 ≤ x ≤ L, with L = 10 as given in Table 3.1(NC).With these parameters, the bifurcation boundary is at γ = γ0 ≈ 24.706.Next, we plot q versus the real part of σ (denoted <(σ)), where σ is theeigenvalue of M with largest real part. This plot is known as the dispersionrelation. When the HSS is unstable, there must be a range of q such that thecorresponding <(σ) > 0. Additionally, if we are not far from the bifurcationboundary, this range must be small. The plot is shown in Fig. 4.3.Due to the boundary condition, the underlying Sturm-Liouville problemforces the wave number to be discrete, that is, we must haveq = qn =npiL, n ∈ N.The mode qn repeats its period n/2 times. With L = 10, we then expectthe final pattern to repeatN =n2=qL2pi424.2. Turing analysis of wave pinning model and non-conservative extension in 1DFigure 4.3: The dispersion relation for the non-conservative model withdifferent values of γ, with other parameters from Table 3.1(NC). Noticethat at bifurcation point γ0 = 24.706, the curve is tangent to the x-axis atq ≈ 5.72. For γ = 25, the range of unstable q is 4.52 < q < 7.24periods, with N restricted to (half-)integers. For γ just at γ0, the unstablemodes are near q = 5.72 which corresponds to N = 9.1037. For γ = 25, theprediction is 7.5 ≤ N ≤ 11.5.We run numerical simulations with a random noise initial conditionu = (1 + Unif[−0.01,+0.01])u∗to verify the prediction. Indeed, for γ = 24.7 < γ0 (or smaller), the noisequickly decays and no pattern forms. The same happens for γ just aboveγ0, which is also expected since the unstable eigenvalue is vanishingly small,so patterns can form only after a very long time. With γ = 25, patternsreadily form, but they usually repeats for 3 to 5 periods, contrary to theprediction.However, upon a closer inspection, we found that there is an initialprecursor pattern which forms gradually with the number of periods withinexpectation (usually 9 ≤ N ≤ 11), with amplitude on the order of 0.1. Oncea peak reaches a certain amplitude, it will very rapidly grow to the full size ofthe final pattern (amplitude ≈ 1.3) while suppressing nearby peaks. Otherprecursor peaks farther away from the grown one might survive longer andeventually transition to full size, or be suppressed by another nearby peak434.2. Turing analysis of wave pinning model and non-conservative extension in 1D(a) (b)(c) (d)Figure 4.4: Simulations for the non-conservative model (2.7) with randominitial conditions u = u∗(1 + Unif(−0.01,+0.01)), v = v∗ with default pa-rameters from Table 3.1(NC) except γ = 25 in (a,b,c) and γ = 24.75 in (d).The color range is chosen so that the precursor pattern is more visible. Therapid transition from the shallower, higher frequency precursor pattern tothe final pattern can be clearly seen.which has transitioned sooner. Since the precursor pattern is very short-lived and shallow compared to the final pattern, it is easily overlooked.Fig. 4.4 shows a few simulations focusing on the transition between theinitial and final pattern. Here we see that the initial pattern has the ex-pected number of peaks, while the final pattern has much less. The time oftransition has some variability, and different regions can transition at differ-ent times. With γ closer to the bifurcation boundary, it will take the initialpattern longer to form, but the transition is just as rapid. The effect of agrowing peak suppressing nearby peaks is clearly visible. There is also somemerging between adjacent peaks. Whether a peak survives or not dependson other nearby peaks, making this a non-linear interaction. Therefore the444.2. Turing analysis of wave pinning model and non-conservative extension in 1Dsystem behavior after transition is beyond the scope of the linear Turinganalysis.As a result, Turing analysis cannot directly predict the number of peaksin the final pattern. In order to find out how many peaks are possible, we useperiodic initial condition u = u∗(1 + 0.1 cos(Nx/L)) to “force” a N -peakssolution. The result is that the final pattern can have anything between 3and 11 peaks. Too few peaks leaves a large enough empty space betweenthem for a new peak to form, while too many peaks leads to some peaksbeing suppressed after forming.However, there is another way to use the Turing analysis results withthe “minimum patch size” idea from [35]. The observation that there existan upper bound on the number of peaks implies that there is a minimumpatch size Lm, which is the minimum domain length for L that can supporta single peak. A solution with a single peak corresponds to the mode N = 1,so the Lm is the minimum L, such that q = 2Npi/L is within the unstablerange. For our parameters, the unstable range is 4.52 < γ < 7.24, soLm = 2pi/7.24 ≈ 0.8678. Our domain with length L = 10 can fit 11 sub-domains of length Lm inside, and each can support one peak. Therefore,this calculation predicts that the final pattern can have at most 11 peaks,exactly as observed. This method cannot predict the minimum number ofpeaks, however.45Chapter 5Local perturbation analysisLocal perturbation analysis (LPA) is a method for examining the evolutionof a localized perturbation to a homogeneous steady state (HSS) for a fast-slow diffusion-reaction system. It provides a way to systematically detectcertain forms of nonlinear instabilities that are not detectable by Turinganalysis. LPA was first developed by [14], and has been used in [9, 17, 16,27] to analyze wave pinning and related models. We will compare our resultsto the results these earlier papers in the following sections. A comprehensiveguide to LPA and other examples of application are given in [18].In this section, we will compare the stability results from LPA with anal-ysis of the well-mixed systems and Turing analysis. In particular, LPA pro-vides a rough guide for parameter exploration and hints on possible regimesfor the full PDE model.5.1 Overview of LPAThe main idea behind LPA is to exploit the difference in the rates of diffu-sion, which allows us to go one step further than analyzing the well-mixedsystem (WM, (2.11)). Take the limit of the ratio diffusion coefficients (δin our model) to 0 and consider the problem on an intermediate time-scale,i.e rescale time by t =√δτ . In this case, as a leading order approxima-tion, the fast diffusing quantities (v in our model) diffuse infinitely fast,whereas the slowly diffusing quantities (u, F ) do not diffuse at all. This canbe understood as a sort of “zeroth order” asymptotic approximation in√δ.In this approximation, v will always be spatially homogeneous since anyinhomogeneity will be instantly smoothened out by the diffusion. We pre-scribe an initial condition where u and F are homogeneous across the entiredomain with values (u∗, F∗), except for a localized perturbation at somepoint x0 in the interior of the domain. This perturbation is assumed to have465.1. Overview of LPAinfinitesimal width but finitely positive height. That is,v(x, 0) = v∗,u(x, 0) = u∗ + (u˜− u∗)δ(x− x0),F (x, 0) = F∗ + (F˜ − F∗)δ(x− x0),(5.1)where δ(·) denotes Dirac delta function. Note that this differs from Turinganalysis, which assumes a global perturbation with infinitesimal height. Foreach of the slow variables, we can track the global level (u, F ) and the levelat the local perturbation (uL, FL) with two separate variables that dependsolely on time, with initial conditions:u(t) = u(x, t), x 6= x0, u(0) = u∗,v(t) = v(x, t), v(0) = v∗,F (t) = F (x, t), x 6= x0, F (0) = F∗,uL(t) = u(x0, t), uL(0) = u˜,FL(t) = F (x0, t), FL(0) = F˜ .(5.2)Since the perturbation is interpreted to be very narrow, the local variableswill have vanishing effect on the rest of the system, so only the global levelsappear in the reaction term of v, which in turn affects both local and globalu, F . Hence, in general, the LPA for the combined system (2.8) is:∂u∂t= f(u, v, F )− cθu,∂v∂t= −f(u, v, F ) + cα,∂F∂t= (knu− ksF ),∂uL∂t= f(uL, v, FL)− cθuL,∂FL∂t= (knuL − ksFL).(5.3)The first three equations are exactly the same as the WM system (2.11).The resulting ODE system (5.3) is much easier to analyze than the fullPDE system. For each set of interesting parameters, we produce a series ofbifurcation diagrams with respect to different parameters with XPPAUT [10,7] and Matcont [5]. These diagrams identify parameter regimes such thatthe ODE behavior is qualitatively different.475.1. Overview of LPAIn this chapter, we follow AUTO’s convention on bifurcation diagrams.On one-parameter diagrams, red/black curves indicate positions of stable/unstable equilibria respectively, while green/blue indicate the range of sta-ble/unstable limit cycles. On two-parameter diagrams, red/light blue/darkblue curves trace the position of limit points (fold points)/branch points(transcritical points)/Hopf points, respectively.These results do not translate directly to the full PDE model for δ > 0.In some cases distinct ODE regimes correspond to the same behavior for thePDE. In addition, LPA cannot help to detect any bifurcation with respectto δ, whereas Turing analysis (Ch. 4) can. In general, LPA results need tobe verified with numerical PDE simulation with small δ.Figure 5.1: Illustration of the idea behind LPA, showing a localized per-turbation to a HSS. The fast diffusing variable v is spatially homogeneous,while the slowly diffusing variable u can be decomposed into a global uGand local uL variables which do not depend on x.There are some properties that LPA systems exhibit in general. First,for any solutionu = φ(t), v = ψ(t), F = ϕ(t)of WM,u = φ(t), v = ψ(t), F = ϕ(t), uL = φ(t), FL = ϕ(t)is always a solution of LPA. This means that any features (equlibria branches,limit cycles, bifurcation points) that appear on bifurcation diagrams of WM485.2. LPA for Wave pinningmust also appear on bifurcation diagrams of LPA. We call these “globalbranches”. LPA might (and usually does) have additional branches notpresent in WM, we call these “local branches”.Note that the dynamics of the global variables are the same as the WMsystem. This means the eigenvalues of WM are a subset of the eigenvalues ofLPA. In fact, it is easy to see that the Jacobian of LPA is block upper diag-onal, with the top diagonal block being the Jacobian of WM. Consequently,if an equilibrium/limit cycle is unstable in WM, it must also be unstable inLPA. The converse is not true, however, as stability in WM does not implystability in LPA.The structure of the LPA system allows us to decompose the LPA dy-namics into two steps, at least in the case where a stable equilibrium existsfor WM. First, the global variables reach a stable equilibrium independentof the local variables, with v = v∗. Next, consider the system with only thelocal variables, and set v = v∗ as a constant parameter. If the global vari-ables do not reach any stable equilibrium, they will enter the local equationsas a non-autonomous forcing term. This will be helpful during some of theanalysis.For our models, notice that L, the domain length, has no effect on theLPA dynamics. With the diffusion terms absent, L can be absorbed byrescaling of time, so this chapter assumes L = 1 for simplicity.5.2 LPA for Wave pinningThe LPA system corresponding to the wave pinning model (2.5) is simply(5.3) with c = s = 0, and F, FL equations removed. However, this can besimplified with mass conservation similar to Sec. 2.6. Introduce an additionalconstant parameter w = u + v for mean concentration of GTPases, andeliminate v from the equations. This results in a system of 2 ODEs:∂u∂t= f(u,w − u), (5.4a)∂uL∂t= f(uL, w − u). (5.4b)For this model, we identify γ (the magnitude of the only nonlinear term)and w (total concentration) as primary parameters of interest.In Fig. 5.2 and 5.3, there are seven regimes identified, which are sum-marized in Table. 5.1. The global branch can be easily distinguished, andits stability compared to WM agrees with earlier discussions. Compared to495.2. LPA for Wave pinning(a) WM, w = 2 (b) LPA, w = 2(c) WM, w = 3.5 (d) LPA, w = 3.5(e) WM, w = 4 (f) LPA, w = 4Figure 5.2: Bifurcation diagrams of the wave pinning WM and LPA systems(5.4) using bifurcation parameter γ, and other parameters taking valuesfrom Table 3.1(WP) except w. The purple lines are located at bifurcationpoints, separating the different parameter regimes. Notice that all curvesin the WM diagrams (the global branches) also appear in LPA. Moreover,the global branch in LPA can be unstable in some intervals despite it beingstable in WM, as discussed in Sec. 5.1.505.2. LPA for Wave pinning(a) WM (b) LPA(c) LPAFigure 5.3: Two-parameter bifurcation diagrams of the wave pinning (a)WM and (b,c) LPA systems using bifurcation parameters w, γ, and otherparameters taking values from Table 3.1(WP). Each curve in these diagramstraces the location of a bifurcation point shown in Fig. 5.2, and forms theboundary of a parameter regime. The one-parameter bifurcation diagrams inFig. 5.2 corresponds to vertical cross-sections of the diagrams here. The LPAregimes identified I to VII match with the regimes labelled in Fig. 5.2(b,d,f).These regimes and their meanings are summarized in Table 5.1. (c) A zoominto the region near the cusps in (b). We can compare (b) to Fig.3(a) of[17], which is a LPA diagram for the same model but different parameters.Observe that the two diagrams agree on the fold curves (red), but [17] omitsthe transcritical curve (light blue), which separates several distinct regimes.515.2. LPA for Wave pinningthe earlier study [17], we have additionally traced a branch of transcriticalbifurcation in two-parameter continuation in Fig. 5.3, which allows us toidentify several new regimes.Regime Classification DescriptionI Stable One stable GB, no LBII PolarizableOne stable GB, one stable LBlocated above the GBIII PolarizableOne stable GB, three stable LBslocated on both sides of the GBIV PolarizableTwo stable GBs, three stable LBs: one above bothGBs, one in between, and one below both GBsV PolarizableOne stable GB, one stable LBlocated below the GBVI UnstableThe only GB is unstable, two stable LBslocated on both sides of the GBVII UnstableThree GBs, all unstable, four stable LBslocated on both sides of the GBTable 5.1: Summary of the wave pinning model regimes identified in Fig. 5.2and 5.3. Abbreviations: GB: global branch; LB: local branch. Stable: allstable branches are global branches. Polarizable: there exist both stableglobal and local branches. Unstable: all global branches are unstable, inwhich case there must exist some stable local branches.In Regime I, there is only one global branch, which is stable. In thisregime wave pinning does not occur, and any perturbation to the HSS simplydecays exponentially. We call this the stable regime.In Regimes II/III/V, there is one stable branch of global solution, butstable local branches also exist. Their position in relation to the globalbranch differs by regime: in Regime II the local branch is above the globalbranch, whereas the opposite happens in Regime V. In Regime III there arelocal branches located on both sides of the global branch. The boundariesof the domains of attraction are given by the unstable branches between thestable branches. This means that a sufficiently large local perturbation in theright direction can lead to pattern formation. In Regime II, if the local initialperturbation is positive and above the threshold, then the LPA system willbe attracted to an equilibria with uL > u. The equivalent initial conditionfor the full PDE results in wave pinning, as verified by PDE simulations(Fig. 3.1). Regime V is similar except that the initial perturbation needs to525.2. LPA for Wave pinningbe sufficiently negative. In this regime, simulations for the full PDE modelresult in wave pinning only with a large negative perturbation, as observedin Fig.3.2. In Regime III, wave pinning can occur for both large positive andnegative perturbations. Regime IV is similar to III except that there aretwo stable global branches, which means that wave pinning can occur fortwo different global states. PDE simulations confirms that these predictionsare correct. We call regimes with both stable global and local branches thepolarizable regimes.In Regime VI, the global branch is unstable, which means that anyperturbation to the initial HSS can lead to pattern formation. Regime VIIis the same, except that there are three global branches, all unstable. Thisregime behaves the same as VI for the full PDE. The HSS is also Turing-unstable for small δ, following the analysis in Ch. 4. In fact, we can showthat Regimes VI and VII corresponds exactly to the Turing regime as δ → 0.First, since the global branch (u∗, v∗) is stable in WM, we know:df(u,w − u)du∣∣∣∣(u∗,v∗)=∂f(u∗, v∗)∂u− ∂f(u∗, v∗)∂v< 0 . (5.5)Second, the global branch is unstable in LPA. The Jacobian for the LPAsystem is[∂f(u,w−u)∂u∂f(u,w−u)∂v∂f(uL,w−u)∂u∂f(uL,w−u)∂uL] ∣∣∣∣(u=uL=u∗,v=v∗)=[∂f(u∗,v∗)∂u − ∂f(u∗,v∗)∂v 0−∂f(u∗,v∗)∂v ∂f(u∗,v∗)∂u].(5.6)We have already shown that the first eigenvalue ∂uf −∂vf is negative (5.5).So the global branch is unstable in LPA if and only if the second eigenvalueλ = ∂uf(u∗, v∗) > 0. Hence this regime is characterized by0 < ∂uf(u∗, v∗) < ∂vf(u∗, v∗),which is exactly the condition for Turing (4.9) as δ → 0. This result can besimilarly generalized to other models.Bifurcation diagrams with respect to the other parameters (not shown)tell the same story: that there are Turing regimes, stable regimes, andpolarizable regimes separated by fold and transcritical bifurcations. A fewearlier studies have performed LPA analysis with different parameters: [9]used k as bifurcation parameter, and [17] used w. Both have results andinterpretations agreeing with ours.535.3. LPA for non-conservative extension5.2.1 More on Wave pinningLPA allows us to further understand the wave pinning system. We candecompose of LPA dynamics into “global” and “local” steps as discussedin Sec. 5.1. For wave pinning, this is particular helpful since each step is asingle ODE, allowing us to visualize the stability of the equilibria by plottingthe right hand side of the ODE.For the global dynamics, we plot g(u) = f(u,w − u) against u. Weknow from Sec. 4.2.1 that, generically, g(u) has one or two stable zeros. Thespecial cases where g(u) has a degenerate (i.e. double) root correspond tofold bifurcation points. The zeros corresponds to the global equilibria.Next for the local dynamics, for each root u∗ of g(u), we plot f(u,w−u∗).One of its roots will be u∗, the others (if any) give locations of the localequilibria. The equilibria are stable if the derivative at the correspondingroot is negative. The plots corresponding to each regime can be found inFig. 5.4. This provide a way to enumerate all distinct regimes.With a sharper Hill function (i.e. larger n), the “bumps” in Fig. 5.4(c-g)will be higher (e.g. the bump is higher in (d) than (c)). This means therewill be a larger parameter range such that f crosses the x-axis three times.This explains why larger n makes pattern formation more robust, as claimedin Sec. 3.2.5.3 LPA for non-conservative extensionThe LPA system for the non-conservative extension is (5.3) with c > 0, s = 0and F, FL removed. For this system, we identify γ as the main bifurcationparameter of interest as before, as well as c, which controls the magnitudeof the additional terms.As before, we can decompose LPA dynamics into global and local dy-namics. The global equilibrium (u∗, v∗) (4.13) is unique. Then any localbranches uL∗ must satisfy f(uL∗, v∗) = 0. Expanding this yieldsA(uL∗)cα+ ηu∗A(u∗)− ηuL∗ − cθuL∗ = 0 .After some manipulations, we get(cθ + η)(u∗A(uL∗)A(u∗)− uL∗)= 0, (5.7)which simplifies toA(uL∗)A(u∗)=uL∗u∗. (5.8)545.3. LPA for non-conservative extension(a) (b) (c) (d)(e) (f) (g)Figure 5.4: Plots of the reaction term f(u,w− u) (global) and f(u,w− u∗)(local) (2.5c), with parameters from Table. 3.1(WP) illustrating each LPAregime identified for the wave pinning model. We do not plot the localcurves corresponding to unstable global equilibria since they are irrelevantfor determining dynamics in general. The roots of the global and localcurves correspond to the position of global and local equilibria. Notice thatthe number and location of equilibria agrees with the bifurcation diagramsFig. 5.2, 5.3.Recall thatA(u) = k + γun1 + un, u∗ =αθ.Therefore (5.8) shows that the local branches are independent of both c andη.A few more things can be inferred from (5.8). First, for γ k, the LHS≈ 1 so uL∗ = u∗, which means that there is no local branch for small γ. Wewill show that for γ k, there are always a high and low local branches.The low branch is uL ≈ 0, since with γ → ∞, uL = 0, both sides of (5.8)evaluate to 0. With a bit of further manipulation, we can get (in the limitγ →∞):un−1L∗1 + unL∗=un−1∗1 + un∗. (5.9)The functionh(u) =un−11 + un555.3. LPA for non-conservative extensionsatisfies h(0) = 0 = h(u→∞), and it has a single peak at up ≥ 1 (providedn ≥ 2). Following the discussion in Sec. 4.2.2 we will focus on parameterswith α < θ, that is u∗ < 1. Therefore there exist a point uL∗ > up > 1 suchthat h(uL∗) = h(u∗), which corresponds to the high local branch.In the following bifurcation diagrams, we use parameters from Table 3.1(Combined2) but with η = 5. These parameters yield bifurcation diagramsthat are easy to view (the regimes are not too wide or too narrow). Thesame regimes are present for the parameters from Table 3.1(NC), which wasused for simulations in Sec.3.6, but the resulting diagram is hard to read.Fig. 5.5 identifies four distinct regimes of the LPA system and theycan be similarly interpreted as in the previous section. These regimes aresummarized in Table 5.2.Regime Classification DescriptionI Stable One stable GB, no LBII PolarizableOne stable GB, one stable LBlocated above the GBIII UnstableThe only GB is unstable, two LBslocated on both sides of the GBIV UnstableOne GB, two LBs all unstable,each enclosed by a periodic orbitTable 5.2: Summary of the non-conservative model regimes identified inFig. 5.5. For abbreviations see caption of Table. 5.1.In Regime I, there are no local branches and the only global branch (thehorizontal line in Fig. 5.5a) is stable. For the full PDE, all perturbation tothe HSS decays and no pattern formation is observed. Regime II is a polar-izable regime. Any perturbation with insufficient height decays. However,in the full PDE, if the perturbation is sufficiently tall and narrow, it willpersist, leading to the soliton solution shown in Fig. 3.6(f).Regimes III and IV has the global branch being unstable, so they bothcorrespond to the Turing regime. In full PDE simulations, any perturbationleads to pattern formation, which tends to spread across the entire domain,as seen in Fig. 3.6(a,b,e). Regime IV is bounded by a pair of Hopf bifur-cations. Here the global branch is unstable even in WM. The full PDE’sbehavior in this regime is not qualitatively distinct from Regime III. This isa little surprising, since we might expect to see non-static patterns, such asoscillating spikes given the limit cycles.Fig. 5.5(b) shows that c has no effect on the location of the fold and565.3. LPA for non-conservative extension(a) WM (b) LPA(c) LPAFigure 5.5: Bifurcation diagrams for the non-conservative extension, withother parameters from Table 3.1(Combined2) except η = 5. (a) WM, (b,c)LPA, using bifurcation parameters (a,b) γ, with c = 1, (c) c and γ. Noticethe thin polarizable regime II sandwiched between the stable I and TuringIII regimes. The triplet of Hopf bifurcations, which is not present in wavepinning, turns out to be inconsequential.575.4. LPA for actin feedback extensiontranscritical points in (a). This is expected, since c do not affect the locationor existence of the local branches. The behavior as γ → 0, γ →∞ also agreeswith analysis.5.4 LPA for actin feedback extensionFor the actin feedback extension, we use mass conservation to eliminate vfrom the LPA system as before, and identify s, the strength of actin feedbackand k, the basal rate of activation as main parameters of interest. The resultis a system of four ODEs:∂u∂t= f(u,w − u, F ), (5.10a)∂F∂t= (knu− ksF ), (5.10b)∂uL∂t= f(uL, w − u, FL), (5.10c)∂FL∂t= (knuL − ksFL). (5.10d)It is rather hard to analyze the equilibria of this system by hand. So we willjump right to the bifurcation diagrams, shown in Fig. 5.6 and 5.7. We onlydistinguish between the regimes separated by fold and transcritical curvesand omit the Hopf curves since it is unlikely that they will be significant, asexplained below. We also ignore the very narrow regimes since they are likelyinconsequential overall. The six regimes identified this way are summarizedin Table. 5.3.One interesting characteristic of these diagrams is the unstable periodicorbits that exist for very narrow parameter ranges. They emerge as sub-critical Hopf bifurcations. The unstable cycle enlarges until it collides witha saddle point, turning into a homoclinic orbit to the saddle, then it disap-pears beyond that. This is known as saddle-loop bifurcation, or homoclinicbifurcation (see [21, Ch.6.2]). Since the parameter regimes where the peri-odic solutions exist are very narrow, it is unlikely that a real cell would fallinto one of these regimes, so we believe that the Hopf bifurcations are likelyinconsequential to biology.Compared to Fig. 5 of [16], which is a LPA diagram in k − s planecontaining only one of the Hopf curves (dark blue), our diagram (Fig. 5.7)additionally traces the fold (red) and transcritical (light blue) bifurcationpoints and identifies a larger number of distinct regimes. The authors of[16] identified the large Hofp bifurcation curve that forms a loop in Fig. 5.7,585.4. LPA for actin feedback extensionRegime Classification DescriptionI Stable One stable GB, no LBII PolarizableOne stable GB, one stable LBlocated above the GBIII UnstableThe only GB is unstable, two stable LBslocated on both sides of the GBIV UnstableThe only GB is unstable, one stable LBlocated above the GBV PolarizableTwo stable GBs, three stable LBs: one above bothGBs, one in between, and one below both GBsVI PolarizableOne stable GB, three stable LBslocated on both sides of the GBTable 5.3: Summary of the actin feedback model regimes identified inFig. 5.6. For abbreviations see caption of Table. 5.1.tracing the pair of Hopf points on the global branch in Fig. 5.6(d). However,since we now know that the periodic solutions exist only in a very narrowrange of parameters, the Hopf curve is unlikely to form the only regimeboundary as discussed by [16].Fig. 6(b) of [16], which is a single parameter bifurcation diagram withrespect to k, tells a similar story to our Fig. 5.6(b,d) in the sense that ittoo identified subcritical Hopf bifurcations, as well as stable, polarizable andunstable regimes separated by fold and transcritical bifurcations.The two parameter LPA diagram is investigated further in [27], whichhas computed eigenvalues of both the Turing and LPA systems and identifiedtwo Bogdanov-Takens (BT) and generalized Hopf (GH) bifurcations locatedon the main Hopf curve (the one that traces the pair of Hopf points onthe global branch). The BT points are visible in Fig. 5.7; they are the endpoints of the Hopf curves, where they bifurcate off fold curves. The GHbifurcations are located on the Hopf curves, and not visible in our figures.At a GH point, there is another curve that branches off from the Hopf curve,which corresponds to a fold bifurcation for periodic solutions [21, Ch. 8.3].This is a global bifurcation, which is hard to trace with AUTO. It is likelythat there are codimension-2 bifurcation points lying on the other Hopfcurves not discussed by [27].The interpretations of the LPA diagrams are similar to the wave pinningmodel. Whenever we have a stable local branch, we expect some kind ofpattern formation. Unlike the wave pinning model, there are multiple pos-595.4. LPA for actin feedback extensionsible patterns in this model. In this case, LPA cannot accurately predictthe type of pattern. Especially, the consequence of the subcritical Hopf bi-furcations to the full PDE is not clear, but it might suggest some kind of(quasi-)periodic behavior.We match the regimes shown in Fig. 5.6(e) with the numerical simula-tions in Fig. 3.4 in an attempt to identify distinct patterns with their corre-sponding regimes. Additionally, we use data from [16], which performed aparameter sweep in this plane and automatically classified the behavior intofive categories: the four in Fig. 3.4, as well as no patterning.Unsurprisingly, Regime I results in no patterning. Regime V roughlycorresponds to usual wave pinning with static wave front. In Regime VIwe get wave pinning with oscillating front (WPO, Fig. 3.4(a-c)). The cor-respondence is not exact, as expected since LPA is an approximation forδ → 0, while the full PDE uses a finite positive value for δ.However, Regimes II, III and IV cannot be conclusively identified withany one of the behaviors. The data from [16] suggests that II results in nopatterning, while the regime corresponding to wave trains (WT, Fig. 3.4(j-l)) is contained within III. Regimes III and IV can both produce reflectingwaves (RW, Fig. 3.4(d-f)) and single pulse (SP, Fig. 3.4(g-i)). It is possiblethat one of the Hopf curves omitted in Fig. 5.6 (e) (but shown in Fig. 5.7)forms the boundary of a distinguishable regime, but there are too many Hopfcurves to be tractable and none seems to explain the transition between RWand SP.Additionally, [16] observed exotic patterns near the boundary of theregimes. In Fig. 3.5, we showed some possible non-standard patterns (givensuitable initial conditions) well within the interior of II. Moreover, in RegimeII the system is polarizable (akin to Regime V of the wave pinning model inFig. 5.3,5.2) so we expect for s→ 0 to recover wave pinning behavior, whichis observed in full PDE simulations. This means that there must be anotherbifurcation boundary inside Regime II near s = 0, likely of non-linear nature,that indicates a transition from WP to more complex behaviors.In summary, LPA worked well in identifying no-pattern and WP regimes.It is not very useful for predicting the more complex patterns. Many of thosepatterns involve interacting waves, which suggests that they are non-linear,non-local phenomena, explaining why LPA cannot account for them.605.4. LPA for actin feedback extension(a) WM, k = 1.5 (b) LPA, k = 1.5(c) WM, k = 6 (d) LPA, k = 6(e) LPAFigure 5.6: Bifurcation diagrams of the actin feedback extension. (a-d)Using bifurcation parameter s; (e) two-parameter continuation using k, s.In (e), the Hopf curves are omitted since they are too close to the othercurves and clutter up the diagram. The regimes that are too narrow arenot labelled since they are unlikely to be consequential. Notice the nearlyvertical blue curves indicating unstable periodic orbits.615.4. LPA for actin feedback extensionFigure 5.7: Same as Fig. 5.6(e) but with the Hopf curves included. A fewHopf curves lie very close to one of the other curves for most of their length,creating some very narrow regimes.625.5. LPA for the combined Model5.5 LPA for the combined ModelThe LPA diagrams for the combined model are very complex, and mostlybeyond the scope of interpretation. This is unsurprising given the complexbehavior exhibited by the PDE. The bifurcation diagram shown in Fig. 5.8,which uses parameter values from Table 3.1(Combined2), contains manylimit cycle bifurcations, such as torus and period-doubling. One thing thediagram can provide is the minimum value of s required for any non-staticpatterns (corresponding to the first triplet of Hopf bifurcation in Fig. 5.8).With s below this value, the system behaviour is the same as that of thes = 0 case, which reduces back to the non-conservative extension.Figure 5.8: Bifurcation diagrams for the combined model using bifurcationparameter s, with other parameter values taken from Table 3.1(Combined2).Notice the many branches of periodic solutions and the bifurcations on them.In a parameter range around s = 20, there are no stable equilibria nor stableperiodic solutions even though the system remains bounded, which suggeststhe presence of chaos.63Chapter 6Interaction with changingdomainSo far, the simulations and analysis assume a static domain for the PDE.However, given that the GTPases regulates cell shape and size, more re-alistically the domain should be dynamic as well. This chapter presentstwo approaches for adding the coupling between the changing domain andGTPase concentration to existing models.6.1 Well-mixed wave pinning model in 1D withdynamic domainThe first approach is to consider the well-mixed wave pinning model, thatis (2.13), repeated here:∂u∂t= f(u, v) = g(u), v = w − u . (6.1)The goal is to couple this equation with the size of the domain, L, which isno longer held constant. Here for convenience we interpret the domain tobe one dimensional (so L is for length). The interpretation and analysis areidentical for 2D domains (simply change L to A for area).This approach is related to [49]. The authors considered the couplingbetween Rho GTPase and mechanical tension (which depends on cell sizethrough L− L0, where L0 is the basal cell size) by adding a feedback termto the activation function. We will instead consider the dilution effects ofchanging cell size to GTPase dynamics. We will compare the results at theend of Sec.6.1.1. Another related study is [6], in which the authors alsoconsidered the dilution effects, but instead of GTPase dynamics they useda very different dynamics for force-producing molecules (i.e. actin, myosin).The system (6.1) has been reduced to a single ODE using mass conser-vation. However, u and v are concentrations, but mass conservation is astatement about total mass, which is no longer equivalent should the size ofthe domain be dynamic.646.1. Well-mixed wave pinning model in 1D with dynamic domainWe revise the equations as follows to obtain the correct conservationstatement in this context. Let T be the total mass of GTPase, which isconstant. Define U(t), V (t) to be the total amount of active/inactive GT-Pase, related to old variables via w = T/L, u = U/L, v = V/L. The originalODE should be restated in terms of U and V . Remember f(u, v) has unitof concentration/time, so the revised equation should be:∂U∂t= Lf(U/L, V/L), V = T − U . (6.2)Rewrite in terms of the old variables, and expand the left hand side usingthe product rule of differentiation gives∂u∂t= f(u, v)−uL∂L∂t︸ ︷︷ ︸(∗), v =TL− u . (6.3)Notice that compared to (6.1), the only difference is the additional (*) term,which represents dilution effects. If the domain is expanding, i.e. ∂tL > 0,then the (*) term is negative since u is diluted on an expanding domain,and vice versa: a contracting domain concentrates u, which corresponds toa positive (*) term.We need an equation for L to close the system. The simplest choiceis to assume that the cell behaves like a Hookean spring, with rest lengthdetermined by u:∂L∂t= −κ(L− L0(u)) . (6.4)The rest length L0(u) should either be monotone decreasing if the GTPasecauses contraction (Rho), or monotone increasing if the GTPase causes ex-pansion (Rac/Cdc42). To be consistent with the non-linearity in the reactionterm, which is a Hill function centered at u0, we will also use a Hill functionfor L0, centered at some uc. For simplicity, we can assume all nonlinear656.1. Well-mixed wave pinning model in 1D with dynamic domainterms have the same Hill coefficient. Therefore, the full system is:∂u∂t= g(u)− uL∂L∂t, (6.5a)∂L∂t= −κ(L− L0(u)), (6.5b)g(u) =(k + γunun0 + un)(TL− u)− η = f(u,TL− u), (6.5c)L0(u) =Lb + (Lh − Lb)(1− ununc +un)for Rho (contraction) ,Lb + (Lh − Lb) ununc +un for Rac/Cdc42 (expansion) .(6.5d)Here Lb means basal (minimum) rest length and Lh means heightened(maximum) rest length. We will refer to this system (6.5) as the dynamicalcell size model, and its equilibria must satisfyL = L0(u), u =k + γ unun0+unη + k + γ unun0+unTL0=TL0(1− ηη + k + γ unun0+un).(6.6)Although we have scaled u0 away during non-dimensionalization in Sec. 2.2.1,here we put it back explicitly to remind us of the true dimensions of theterms. Nonetheless we will still take its value to be 1.In the following analysis, we will make use of the sharp switch approxi-mation. The idea is to let n→∞ so that the Hill function terms approachHeaviside step function H, allowing us to make the following approxima-tions:unun0 + un⇒ H(u−u0), ununc + un⇒ H(u−uc), 1− ununc + un⇒ H(uc−u) .(6.7)These approximations make the system piecewise linear, so easier to analyze.It is still not fully linear, so we can still hope that the approximation willretain some of the interesting characteristics.The parameter values we use are summarized in Table. 6.1. We willexplain how to guide parameter selection with analysis based on the sharpswitch approximations.6.1.1 Contraction case (Rho)First, consider the effect of Rho GTPase, which makes the cell contract.As before, we use γ as the bifurcation parameter, and T as a secondary666.1. Well-mixed wave pinning model in 1D with dynamic domainParameter Dimension Meaning CC ECk t−1 Basal activation rate 1γ t−1 Nonlinear activation rate 0− 25η t−1 Inactivation rate 5Lb l Basal rest length 1 5.5Lh l Heightened rest length 3 16.5u0 G Threshold for activation 1uc G Threshold for length 4 0.85T G Total mass of GTPase 0− 25 0− 80κ t−1 Length reaction rate 5n 1 Hill coefficient 0−∞Table 6.1: Parameter values for the dynamical cell size model. CC andEC stand for contraction case (Sec. 6.1.1) and expansion case (Sec. 6.1.2),respectively. Abbreviation for dimensions: t : time; l : length; G : GTPase.parameter, so we can compare the result to Fig. 5.3(a), the model with Lconstant. We will focus on the case uc > u0.Under the sharp switch approximation, there are three possible branchesof equilibria. Their locations and conditions of existence can be found withsimple algebraic manipulations:(a) u < u0. Then L = Lh, andu1 =TLh(1− ηη + k).The u1 branch is independent of γ, and exists only ifT < Lhu0η + kk= Lhu0(1 +ηk) .This means that the low GTPase activity equilibrium can only existwhen there is not too much total GTPase.(b) u0 < u < uc. Then L = Lh, andu2 =TLh(1− ηη + k + γ).Notice that when γ = 0, we have u1 = u2, yet we should have u1 < u2.Therefore the u1 and u2 branches cannot both exist at low γ. The676.1. Well-mixed wave pinning model in 1D with dynamic domainexistence condition for the u2 branch simplifies toLhu0T − Lhu0 η − k < γ <LhucT − Lhuc η − k .Since u0 < uc, this range is non-empty only if T > Lhuc. (Anotherpossibility is T < Lhu0, in which case the entire range lies below 0,so it is not biologically relevant). Hence the medium GTPase activ-ity equilibrium requires that there is enough total GTPase, and theactivation strength is at a moderate level.(c) u > uc. Then L = Lb, andu3 =TLb(1− ηη + k + γ).The existence condition isγ >LbucT − Lbuc η − k .If the right hand side is non-positive, i.e. T < Lbuc(1 + η/k), thenthe lower bound of γ is greater than zero. Otherwise the u3 branchexists for all positive γ. This means that the high GTPase activityequilibrium requires either high total GTPase or high activation, whichmakes sense.We ask whether there are parameter settings that yield interesting behaviors,i.e. as many branches existing as possible, which means we wantLhuc < T < min(Lhu0(1 +ηk), Lbuc(1 +ηk))=(1 +ηk)min(Lhu0, Lbuc) .(6.8)So we needmax(ucu0,LhLb)<(1 +ηk),and κ does not influence equilibrium behavior. The values in Table. 6.1(CC)are reasonable choices that satisfy our condition (6.8).Fix T = 15, which is within the interesting range where all three branchesare possible. Then u1 exists for all γ, u2 exists if 0.25 ≤ γ ≤ 19, andu3 exists if γ ≥ 0.8181. It turns out that all three branches are stable.Bifurcation diagrams are shown in Fig. 6.1, where we compare the sharpswitch approximation (n → ∞) for n = 20, where the nonlinearity is verysharp but still continuous.686.1. Well-mixed wave pinning model in 1D with dynamic domain(a) n→∞, T = 15 (b) n = 20, T = 15(c) n→∞ (d) n = 20Figure 6.1: Bifurcation diagrams for the dynamic cell size model (6.5) in thecontraction case, (a,b) using bifurcation parameter γ; (c,d) two-parametercontinuation using T and γ, with other parameters from Table. 6.1(CC).(a,c) Using sharp switch approximation for the nonlinear terms (limitingcase where n → ∞) (b,d) Hill function nonlinearity with n = 20. The fourfold points labelled A to D in (b) correspond to curves with the same labelsin the two parameter continuation (d). In (d), the numbers indicate thestable branches that exist within each regime, e.g. 1, 3 means u1 and u3exist in this regime, but not u2. The cusp is at T = 16.16, γ = 0.2797. TheC curve crosses the x-axis just outside the frame.696.1. Well-mixed wave pinning model in 1D with dynamic domainAs the diagrams show, the sharp switch approximation retains mostof the behaviors of the smooth nonlinearity. One major difference is thattaking n→∞ causes the fold points to become degenerate. The branches ofunstable solutions are not shown, but conceptually, they should be horizontallines at u = u0, u = uc. The other main difference is that the curve labelledB in Fig. 6.1(d) becomes perfectly vertical in Fig. 6.1(c). As a consequence,whereas the u1 branch exists for all γ for n → ∞, it only exists for a shortrange for n = 20.Compare Fig. 6.1(d) to the earlier results in Fig. 5.3(a) for the samemodel (the well-mixed wave pinning model, (2.13)) but with L treated asa constant (albeit with slightly different parameters, though the generalshape still holds). Observe that whereas Fig. 5.3(a) has one pair of foldpoints joined by a cusp point, Fig. 6.1(d) has two.Notice that we did not find any periodic solution, as expected. Consideran increase in Rho activity, which makes the cell contract. This increasewill further concentrate Rho, creating a positive feedback effect. This isin contrast to the results of [49], which considered the coupling betweenmechanical tension and Rho without the dilution effect. In this approach,the feedback from cell size to Rho is negative, which results in an oscillatoryregime. The middle branch does not exist in this model, and there are twomono-stable regimes where only one stable branch exists. These two regimesare also present in our model.6.1.2 Expansion case (Rac, Cdc42)Next we consider the case where the GTPase makes the cell expand (Rac,Cdc42). In this case, a contracting cell concentrates the GTPase, driving itslevel up, and making it expand. This means that there is a negative feedbackeffect, which potentially leads to periodic solutions. We will consider the caseuc < u0. As before, we have three branches of equilibria for n→∞:(a) u < uc. Then L = Lb, andu1 =TLb(1− ηη + k0),independent of γ, exists only ifT < Lbuc(1 +ηk0).706.2. PDE simulations with deforming boundary using the Cellular Potts Model(b) uc < u < u0. Then L = Lb, andu2 =TLh(1− ηη + k0),again independent of γ, exists only ifLhuc(1 +ηk0) < T < Lhu0(1 +ηk0) .(c) u > u0. Then L = Lh, andu3 =TLh(1− ηη + k0 + γ).The existence condition isγ >Lhu0T − Lhu0 η − k0, T > Lhu0 .As before, the existence conditions are within expectation: the low GTPaseactivity equilibrium requires low total GTPase, and so on. We choose theset of parameters from Table. 6.1(EC), which are in line with the aboverequirement. For our choice of parameters, u1 exists for T < 51.425, u2exists for 154.275 < T < 181.5, and with T = 40, u3 exists for γ > 6.02128.We fix T = 40 for the bifurcation diagrams. For the sharp switch ap-proximation, all three branches are stable , and the bifurcation diagram isuninteresting: it looks the same as Fig. 6.1(a) except that the middle branchis missing. However, for n = 8, the system has a pair of Hopf bifurcation,see Fig. 6.2, and periodic solutions exist for an intermediate range of γ, asexpected.6.2 PDE simulations with deforming boundaryusing the Cellular Potts ModelIn this section, we discuss another approach where we couple the PDE forGTPase dynamics with the Cellular Potts Model (CPM) for the boundaryof the cell in 2D. The results in this section are a collaboration with twoother members of our group. The program we used for the simulation wasdeveloped by Zachary Pellegrin, and Elisabeth Rens who also adapted thecode for our models and ran the simulation. I provided the equations andthe parameter values.716.2. PDE simulations with deforming boundary using the Cellular Potts Model(a) T = 40 (b) T = 65(c)Figure 6.2: Bifurcation diagrams for the the dynamic cell size model (6.5) inthe expansion case, (a,b) using bifurcation parameter γ; (c,d) two-parametercontinuation using T and γ, with other parameters from Table. 6.1(EC). (a)T = 40, which is left of both cusps in (c). Notice that AUTO had troublekeeping track of the periodic branch on the right side. The two blue curveslikely form a closed loop. (b) T = 65, which is in between the two cusps.(c) shows that the pairs of fold and Hopf points are joined by a cusp in theT − γ plane.726.2. PDE simulations with deforming boundary using the Cellular Potts ModelAn in-depth introduction to CPM can be found in [25]. The basic ideaof CPM is to describe the cell and its surrounding environment on a regularsquare lattice, where each square can be in one of two states: occupied bythe cell, or not. We then define an energy function of the following form for2D simulations:E(p,A) = Jp+ αp(p− p0)2 + αa(A−A0)2, (6.9)where p and A are the perimeter and area of the cell, respectively, p0 andA0 are the “rest” perimeter and area (akin to the rest length of a spring),and J, αp, αa are constant parameters.At each time step, we consider a local shape change at every point onthe boundary of the cell, where the cell can either expand or contract byone square (a free square becomes occupied or vice versa). Whether thischange takes place is decided randomly via a Monte Carlo algorithm, withits probability depending on the change in energy ∆E if the change were totake place:p(∆E) = min(1, exp(−∆E + YT)), (6.10)where T is the “temperature” parameter describing how much randomnesswe allow. If the change decreases the energy by at least Y , then it is alwaysaccepted. Otherwise the change is accepted with probability determined bya Boltzmann function, which decreases monotonically as ∆E increases. Ahigher temperature means energetically disfavoured changes are more likely.Next, we couple CPM with the combined PDE model (2.8) for GTPasedynamics. The details of the method can be found in [38]. We use thesame lattice for both CPM and finite difference computation for the PDE.The energy function (6.9) is modified with an extra term depending on u,such that GTPase activity makes either contractions (Rho) or expansions(Rac, Cdc42) more likely. At each iteration, we solve the PDE for 1000 timesteps, with ∆t = 10−5, then update the cell boundary according to CPM.The no-flux boundary conditions for the PDE are enforced by modifying thediscrete Laplacian for each possible local boundary shapes. When the cellexpands into a free square, the GTPase activity in the new square is set tobe the same as the average of the adjacent squares. In order to enforce massconservation, the GTPase activity in the local neighborhood is re-normalizedto have the same total GTPase mass as before. In the case of contraction,the GTPase contained in the square to be freed is redistributed evenly tothe neighboring squares.In the simulations, we use parameters similar to the ones we used fornumerical simulations in Ch. 3. The results are shown in Fig. 6.3, which736.2. PDE simulations with deforming boundary using the Cellular Potts Modelfocuses on the Rho (contraction) case. Notably, we obtained the three mainpatterns: wave pinning, spots and actin waves patterns in accordance withour expectations given the static boundary simulations from Ch. 3.In Fig. 6.3(a,b), we obtain wave pinning. By comparing this figure to1D simulations in Fig. 3.1 and 2D simulations with static boundary inFig. 3.10(a,b), we notice the familiar defining characteristics of wave pin-ning: the cell is polarized into regions of high and low GTPase activity withthe wave front stationary with respect to the cell. The cell as a whole slowlymoves to the right as expected.In Fig. 6.3(a,b), we obtain the spots pattern observed in 2D simulationscorresponding to the non-conservative model in Fig. 3.10(c,d). Notice thateven though we initialized the system at the unique HSS without any per-turbation, the pattern still forms. The reason is that the HSS is unstable,which means the perturbation due to the random change in the boundaryis sufficient to trigger pattern formation. The spots cause the boundary tocave in, which in turn pushes the spots away. In some cases we observedspots being split into two, or two spots being pushed together and mergingdue to boundary changes. It would be interesting to repeat the simulationfor the expansion (Rac, Cdc42) case, where we expect the spots to causeformation of thin protrusions corresponding to filopodia.In Fig. 6.3(e,f), we observe travelling waves that we interpret to representthe “actin waves” as described in Sec. 1.3. The dynamics of these waves arevisually similar to videos of actual actin waves observed in experiments, forexample the video supplements of [26]. We used the same parameters as inFig. 3.4(j-l), where we obtained the persistent wave train behavior.746.2. PDE simulations with deforming boundary using the Cellular Potts Model(a) Wave pinning, t = 0.5 (b) Wave pinning, t = 5(c) Spots, t = 0.06 (d) Spots, t = 2(e) Actin wave, t = 0.1 (f) Actin wave, t = 1Figure 6.3: Simulations of the coupled CPM-GTPase model. For each setof parameters, we show two snapshots of GTPase activity u in the cell, onejust after the pattern has taken shape, and another one sometimes laterto show how the pattern evolves and its effects on the shape and locationof the cell. (a,b) Wave pinning, with c = s = 0, δ = 0.04 and all otherparameters from Table. 3.1(Combined2). Initial conditions: u = 0, exceptu = 3 on a region comprising 10% of the cell on the left edge, and v = 2.(c,d) Spots pattern, with c = 1, s = 0, δ = 0.04 and all other parametersfrom Table. 3.1(Combined2). Initial conditions: u = 0.3333, v = 2.5278.(e,f) Actin waves, with s = 30L2, k = 6L2 and all other parameters fromTable. 3.1(AF). Initial conditions: u = 0, except u = 5 on 10% of the cellon the left edge, and v = 2, F = 0.75Chapter 7DiscussionWe have shown a variety of possible behaviors in the wave pinning modeland its extensions in Ch. 3. In Ch. 4 and 5, we analyzed the bifurcationproperties of the model to determine the number and locations of distinctparameter regimes, and their correspondence with observed patterns. InCh. 6, we investigated the interplay between GTPase dynamics and cellshape. In this chapter, we summarize the mathematical results and interprettheir significance in answering the biological questions we posed in Sec. 1.4.7.1 Summary of model behaviorsIn one spatial dimension, we observed several patterns in the wave pinningmodel combined with the two extensions, some of which were reproduced in2D. Here is a summary of their characteristics, along with the conditions fortheir existence and interpretations. For all patterns, a small diffusion ratioδ between active and inactive GTPase is required.B1. Static wave pinning (Fig. 3.1, 3.2): the domain is partitioned intotwo or more alternating regions of high/low GTPase activity, withrelatively narrow transition zones between them. We interpret thisbehavior as the most basic form of cell polarization. For this to happen,we need the feedback activation γ to be strong enough, and totalamount of GTPase T at a moderate level to allow the existence ofboth low and high GTPase activity equilibria.B2. Wave pinning with oscillating front (Fig. 3.4(a-c)): like (1) but thetransition zone oscillates back and forth. This behavior has similarexistence conditions as (1), but requires non-zero, weak actin feedbacks. This might be a more realistic picture for cell polarization.B3. Reflecting waves (Fig. 3.4(d-f)): a travelling pulse reflected back atboundaries. I am not aware of any corresponding experimental ob-servations. Given a large enough γ and moderate T , this behaviorrequires an intermediate range of s.767.1. Summary of model behaviorsB4. Persistent wave trains (Fig. 3.4(j-k) in 1D, reproduced in 2D with mov-ing boundary conditions in Fig. 6.3(e,f)): a train of travelling wavesabsorbed at the boundaries. We interpret this as representing actinwaves. Beside the high γ and moderate T , this behavior is observedonly when both the basal activation k and actin negative feedback sare in an intermediate range, achieving a balance between them.B5. Single pulse (Fig. 3.4(g-i)): perturbation to the HSS results in a singletravelling pulse, which is eventually absorbed at the boundary and thesystem returns to HSS. The regime on which this occurs is next to theregimes for (3) and (4), and we could not identify a clear boundaryseparating these regimes.B6. Static, spatially periodic spikes (Fig. 3.6 (a,b), Fig. 3.7(a)): A steadystate consisting of a series of equally spaced thin spikes. In 2D, thecorresponding pattern is a 2D lattice of spots. We interpret this to rep-resent filopodia formation, since such a spike in Rac activity will likelycause a thin cellular protrusion. This behavior requires the sourceand sink terms, as well as a high enough γ to put the system into theTuring regime.B7. Lone static spike (soliton, Fig. 3.6 (c,d), Fig. 3.7(b)): a strong andhighly localized perturbation leads to a single spike while the rest ofthe domain remains at the HSS. This too can correspond to filopodiaformation. The corresponding parameter regime is a narrow range ofγ sandwiched between the stable regime and the Turing regime for (6).B8. Repelling pulses (Fig. 3.8(a)): a series of spikes similar to (6), but thespikes are mobile and repel each other when they get too close. Thisbehavior does not resemble any known observations. This happenswhen both the source and sink terms are present and the actin feedbacks is positive but weak.B9. Colliding wave trains (Fig. 3.8(b)): travelling pulses that emerges ei-ther at the boundary or in the interior of the domain, and collide andinteract in complex manner. The overall behavior is similar to (4),but much more unpredictable. The conditions are similar to (8) butrequire a high s.B10. Localized rapidly oscillating standing waves (Fig. 3.9(d)): this is ob-served in a region of the domain occupied by a standing wave with a777.2. Biological and mathematical implicationsshort wave length and high oscillating frequency, and is usually tran-sient. We are not aware of any corresponding experimental observa-tions. The conditions are similar to (9) but require a very high s.B11. Localized oscillating spikes (Fig. 3.11(e,f)): spikes that are confinedlocally but oscillate in amplitude. This occurs when the parameters arespatially inhomogeneous. Further, in a narrow region of the domainthe parameters lie in a pattern-forming regime, and in the rest of thedomain the parameters lie in the stable regime. This behavior is acaricature of filopodia formation in response to viral infection.B12. Other “exotic” behaviors observed near the boundary of parameterregimes in the actin feedback model (Fig. 3.5). We have mostly ignoredthese as they are rather hard to characterize.Note that there are a few behaviors that we could not identify with anyknown experimental observations. There are a few possible explanations.First, the parameter regime might be outside of the biologically relevantrange. For example, (10) requires a very high s, which might not be possibleto achieve. Second, the parameter regime might be too narrow, in whichcase a real cell would be unlikely to fall in this regime. Finally, the realsystem has greater complexity, heterogeneity and stochasticity, which is notcaptured by our simplified model. Some of the regimes, especially the narrowones, might no longer exist after these factors are taken into consideration.Altogether, we have identified behaviors that correspond to polarization,filopodia formation, and actin waves. These results will help us to explainsuch biological phenomena.7.2 Biological and mathematical implicationsWe now return to the questions we posed at the end of Sec. 1.4.Q1. The non-conservative terms result in multiple spikes of high GTPaseactivity (B6). Compared to simple polarization resulting from wavepinning (B1), the spots are much thinner than the polarized regions,but both patterns are static. Actin feedback turns static patternsto travelling pulses and waves. Together, this results in the spikescolliding with each other with various outcomes (B9), as well as thestanding waves pattern (B10), both are quite distinct from static spikesand actin waves. The biological implication of these results are unclear.787.2. Biological and mathematical implicationsQ2. Locally elevated GTPase activation leads to the localized formation ofhot-spots of GTPase activity, which do not spread to other parts ofthe domain (B11). These hot-spots could be associated with filopo-dia formation. This suggests that if a virus can manipulate GTPaseactivation in its vicinity, either directly or through some upstream sig-nalling proteins, then it can trigger the formation of a single filopodiumwhich the virus then exploits to enter the cell.Q3. The periodic spikes pattern in 1D (B6) becomes a lattice of spotsin 2D, whereas we initially expected a stripe pattern given the 1Dbehavior. Waves trains (B4) in 2D can form curves and spirals, whichcannot happen if the system is constrained to a 1D geometry. On theother hand, wave pinning (B1) behaves similarly in both 1D and 2Ddomains.Q4. Formation of most patterns require the cell to be large enough. Forexample, the pattern of spikes (B6) requires the domain to be largeenough to support at least one period of the pattern. The exact sameparameters in a smaller cell can only result in wave pinning. If thecell is too small, then even wave pinning is not possible. This suggeststhat GTPase patterns are simpler in smaller cells.Q5. In the dynamic cell size model, we showed that the coupling betweenGTPase and cell size can result in additional equilibria, which corre-spond to additional possible resting states for the cell. The Hofp bifur-cation observed in the GTPase promoting cell expansion case showsthat oscillations in both GTPase activity and cell size are possible.This behavior cannot be captured with a static domain. The CPMsimulations in 2D showed that the wave pinning, spots, and wave trainpatterns retain their characteristics in a coupled model. However, thefeedback from the deforming boundary can turn an otherwise staticpattern into a dynamic one, such as the spots splitting and mergingdue to boundary movements. Altogether, the behavior is much morecomplex.Q6. The models do not account for the experimentally observed oscillatingspots without heterogeneity in the parameters. By allowing the pa-rameters to vary in space in a very specific way, we can obtain spotsthat oscillate in-place (B11), but the oscillations are much smaller thanobserved in [39]. This suggests that our models lack certain ingredientsthat are essential to this behavior.797.3. Future workQ7. The CPM simulations (Fig. 6.3) indicate that both simple polarizationand actin waves are consistent with cell motility. However, whereassimple polarization results in a regular movement at constant speedand direction, actin waves allow the cell to protrude in random direc-tions, possibly allowing for local “exploration” and random motion.Q8. Unlike LPA, Turing analysis can predict the wavelength of the initialprecursor pattern. However, as the pattern evolves and the spikesinteract, the system leaves the linear regime and the Turing resultsare no longer valid. We can use Turing analysis to indirectly predictthe minimum wavelength of the final pattern.Q9. The unstable regimes identified by LPA correspond to Turing-unstableregimes in the limit δ → 0. In the unstable regimes, any infinitesimalperturbation to the HSS leads to pattern formation. LPA additionallyidentifies the polarizable regimes, which are pattern-forming regimesthat are Turing-stable. Pattern formation in these regimes require asufficiently large perturbation to the HSS. On the other hand, Turinganalysis can describe system behavior for δ > 0, whereas LPA can-not. In Fig. 7.1, we illustrate the relation between LPA and Turingbifurcation diagrams.Q10. The addition of feedback from actin assembly to the system results inmany more bifurcations in LPA diagrams. We were forced to ignoresome features, such as the Hopf curves, to keep the analysis tractable.In the presence of both actin and non-conservative GTPase terms inthe model, the LPA diagrams become hopelessly complex to analyze.This limits the utility of LPA when there are many interacting species.LPA was also not very helpful for predicting the final pattern whenmultiple patterns are possible in the model.7.3 Future workThere are a few directions for future work in this project.We are interested in the properties of the soliton solution, which repre-sents localized GTPase activity. The plan is to use asymptotic analysis toobtain a series expansion of the soliton solution in terms of the small param-eter δ. This would allow us to determine its stability, as well as width andheight. LPA helped to determine the zone of attraction of a solitary peakin the limit δ → 0. In Fig. 5.5(b), Regime II, the black curve determines, in807.3. Future work(a) LPA (b) TuringFigure 7.1: Comparison of LPA and Turing bifurcation diagrams for thenon-conservative model. (a) is a zoom of the LPA diagram from Fig. 5.5(b).(b) is the Turing bifurcation diagram similar to Fig. 4.2, except using thesame parameters as (a). Observe that both the LPA-stable (I) and the LPA-polarizable (II) regimes in (a) located to the left of γc = 16.765 correspondto the Turing-stable regime below the blue curve in (b). The LPA-unstableregimes (III, IV) correspond to the Turing-unstable regime above the curve.The curve passes through δ = 0, γ = γc. The bifurcation boundary betweenRegimes I and II, and between III and IV cannot be detected by Turinganalysis. Given that numerical simulations have shown that the PDE pro-duces the same behavior (Fig. 3.7(a)) in both Regimes III and IV, it ispossible that these are not distinct regimes for the PDE. Overall, the LPAdiagram (a) can be seen as a vertical slice of the Turing diagram (b) atδ = 0, with additional bifurcation boundaries that separates the LPA-stableand LPA-polarizable regimes.this limit, the size of perturbation to the HSS needed to produce the soliton.Asymptotic analysis would provide correction terms to this estimation for afinitely small δ.There are also experimental results such as [39] that fitted the activationfunction of GTPase to observational data. We can use this function toreplace the Hill function term in our model. This would arguably be morebiologically relevant. It would be interesting to see how this replacementchanges the bifurcation behaviors.Another direction is to combine our work in Sec. 6.1 with the work of [49],by considering the dilution effect from cell size and the feedback from me-chanical tension together. Since these two have opposite effects on GTPase817.3. Future workactivity, their interaction would potentially be complex and interesting.In Sec. 6.2, we only ran simulations for Rho (contraction case) so far. Inthe future, we plan to run simulations for the expansion case, as well as theparameter regimes that correspond to other behaviors.82Bibliography[1] Alberts, B. et al. Molecular Biology of the Cell. 5th. New York: Gar-land Science, 2008.[2] Alnaes, M. S. et al. “The FEniCS project version 1.5”. In: Archive ofNumerical Software 3.100 (2015), pp. 9–23.[3] Chang, K. et al. “Filopodia and viruses: an analysis of membraneprocesses in entry mechanisms”. In: Frontiers in Microbiology 7 (2016),p. 300.[4] Crank, J. and Nicolson, P. “A practical method for numerical evalua-tion of solutions of partial differential equations of the heat-conductiontype”. In: Mathematical Proceedings of the Cambridge PhilosophicalSociety. Vol. 43. 1. Cambridge University Press. 1947, pp. 50–67.[5] Dhooge, A., Govaerts, W., and Kuznetsov, Y. A. “MATCONT: aMATLAB package for numerical bifurcation analysis of ODEs”. In:ACM Transactions on Mathematical Software (TOMS) 29.2 (2003),pp. 141–164.[6] Dierkes, K., Sumi, A., Solon, J., and Salbreux, G. “Spontaneous os-cillations of elastic contractile materials with turnover”. In: PhysicalReview Letters 113.14 (2014), p. 148102.[7] Doedel, E. J. “AUTO: A program for the automatic bifurcation anal-ysis of autonomous systems”. In: Congressus Numerantium 30 (1981),pp. 265–284.[8] Edelstein-Keshet, L. Mathematical Models in Biology. Vol. 46. Philadel-phia: SIAM, 1988.[9] Edelstein-Keshet, L., Holmes, W. R., Zajac, M., and Dutot, M. “Fromsimple to detailed models for cell polarization”. In: Philosophical Trans-actions of the Royal Society of London B: Biological Sciences 368.1629(2013), p. 20130003.[10] Ermentrout, B. Simulating, Analyzing, and Animating Dynamical Sys-tems: A Guide to XPPAUT for Researchers and Students. Philadel-phia: SIAM, 2002.83BIBLIOGRAPHY[11] Fletcher, D. A. and Mullins, R. D. “Cell mechanics and the cytoskele-ton”. In: Nature 463.7280 (2010), p. 485.[12] Gierer, A. and Meinhardt, H. “A theory of biological pattern forma-tion”. In: Kybernetik 12.1 (1972), pp. 30–39.[13] Gilbert, S. Developmental Biology. 9th. Sunderland: Sinauer Asso-ciates, 2010.[14] Grieneisen, V. “Dynamics of auxin patterning in plant morphogenesis-A multilevel model study”. PhD thesis. University Utrecht, 2009.[15] Hall, A. “Rho GTPases and the actin cytoskeleton”. In: Science 279.5350(1998), pp. 509–514.[16] Holmes, W. R., Carlsson, A. E., and Edelstein-Keshet, L. “Regimes ofwave type patterning driven by refractory actin feedback: transitionfrom static polarization to dynamic wave behaviour”. In: Physical Bi-ology 9.4 (2012), p. 046005.[17] Holmes, W. R. and Edelstein-Keshet, L. “Analysis of a minimal Rho-GTPase circuit regulating cell shape”. In: Physical Biology 13.4 (2016),p. 046001.[18] Holmes, W. R., Mata, M. A., and Edelstein-Keshet, L. “Local pertur-bation analysis: a computational tool for biophysical reaction-diffusionmodels”. In: Biophysical Journal 108.2 (2015), pp. 230–236.[19] Holmes, W. R., Park, J., Levchenko, A., and Edelstein-Keshet, L. “Amathematical model coupling polarity signaling to cell adhesion ex-plains diverse cell migration patterns”. In: PLoS Computational Biol-ogy 13.5 (2017), e1005524.[20] Jilkine, A. and Edelstein-Keshet, L. “A comparison of mathematicalmodels for polarization of single eukaryotic cells in response to guidedcues”. In: PLoS Computational Biology 7.4 (2011), pp. 1–15.[21] Kuznetsov, Y. Elements of Applied Bifurcation Theory. 3rd. New York:Springer Science & Business Media, 2004.[22] Larsson, S. and Thomee, V. Partial Differential Equations with Nu-merical Methods. Vol. 45. Heidelberg: Springer Science & BusinessMedia, 2008.[23] Li, Y.-X. “Tango waves in a bidomain model of fertilization calciumwaves”. In: Physica D: Nonlinear Phenomena 186.1-2 (2003), pp. 27–49.84BIBLIOGRAPHY[24] Lodish, H. et al. Molecular Cell Biology. 5th. New York: WH FreemanNew York, 2004.[25] Magno, R., Grieneisen, V. A., and Mare´e, A. F. “The biophysicalnature of cells: potential cell behaviours revealed by analytical andcomputational studies of cell surface mechanics”. In: BMC Biophysics8.1 (2015), p. 8.[26] Marchesin, V., Montagnac, G., and Chavrier, P. “ARF6 promotes theformation of Rac1 and WAVE-dependent ventral F-actin rosettes inbreast cancer cells in response to epidermal growth factor”. In: PLoSOne 10.3 (2015), e0121747.[27] Mata, M. A., Dutot, M., Edelstein-Keshet, L., and Holmes, W. R.“A model for intracellular actin waves explored by nonlinear localperturbation analysis”. In: Journal of Theoretical Biology 334 (2013),pp. 149–161.[28] Meinhardt, H. “Orientation of chemotactic cells and growth cones:models and mechanisms”. In: Journal of Cell Science 112.17 (1999),pp. 2867–2874.[29] Mogilner, A., Allard, J., and Wollman, R. “Cell polarity: quantita-tive modeling as a tool in cell biology”. In: Science 336.6078 (2012),pp. 175–179.[30] Mori, Y., Jilkine, A., and Edelstein-Keshet, L. “Wave-pinning andcell polarity from a bistable reaction-diffusion system”. In: BiophysicalJournal 94.9 (2008), pp. 3684–3697.[31] Mori, Y., Jilkine, A., and Edelstein-Keshet, L. “Asymptotic and bifur-cation analysis of wave-pinning in a reaction-diffusion model for cellpolarization”. In: SIAM Journal on Applied Mathematics 71.4 (2011),pp. 1401–1427.[32] Nikolic, D. S. et al. “HIV-1 activates Cdc42 and induces membraneextensions in immature dendritic cells to facilitate cell-to-cell viruspropagation”. In: Blood (2011), blood–2010.[33] Olson, M. F. and Sahai, E. “The actin cytoskeleton in cancer cellmotility”. In: Clinical & Experimental Metastasis 26.4 (May 2008),p. 273.[34] Otsuji, M. et al. “A mass conserved reaction–diffusion system cap-tures properties of cell polarity”. In: PLoS Computational Biology 3.6(2007), e108.85BIBLIOGRAPHY[35] Painter, K. J. and Hillen, T. “Spatio-temporal chaos in a chemotaxismodel”. In: Physica D: Nonlinear Phenomena 240.4-5 (2011), pp. 363–375.[36] Payne, R. J. and Grierson, C. S. “A theoretical model for ROP local-isation by auxin in Arabidopsis root hair cells”. In: PLoS One 4.12(2009), e8337.[37] Rappel, W.-J. and Edelstein-Keshet, L. “Mechanisms of cell polariza-tion”. In: Current Opinion in Systems Biology 3 (2017), pp. 43–53.[38] Rens, E. and Edelstein-Keshet, L. “From Energy to Force in the Cel-lular Potts Model”. In preparation.[39] Robin, F. B., Michaux, J. B., McFadden, W. M., and Munro, E. M.“Excitable RhoA dynamics drive pulsed contractions in the early C.elegans embryo.” In: BioRxiv (2016), p. 076356.[40] Ruuth, S. J. “Implicit-explicit methods for reaction-diffusion prob-lems in pattern formation”. In: Journal of Mathematical Biology 34.2(1995), pp. 148–176.[41] Schwartz, M. “Rho signalling at a glance”. In: Journal of Cell Science117.23 (2004), pp. 5457–5458.[42] Turing, A. M. “The chemical basis of morphogenesis”. In: Philosoph-ical Transactions of the Royal Society of London. Series B, BiologicalSciences 237.641 (1952), pp. 37–72.[43] Verschueren, N. and Champneys, A. “A model for cell polarizationwithout mass conservation”. In: SIAM Journal on Applied DynamicalSystems 16.4 (2017), pp. 1797–1830.[44] Vicker, M. G. “F-actin assembly in Dictyostelium cell locomotionand shape oscillations propagates as a self-organized reaction–diffusionwave”. In: FEBS Letters 510.1-2 (2002), pp. 5–9.[45] Walther, G. R., Mare´e, A. F., Edelstein-Keshet, L., and Grieneisen,V. A. “Deterministic versus stochastic cell polarisation through wave-pinning”. In: Bulletin of Mathematical Biology 74.11 (2012), pp. 2570–2599.[46] Weiner, O. D. et al. “An actin-based wave generator organizes cellmotility”. In: PLoS Biology 5.9 (2007), e221.[47] Xu, H. et al. “Real-time imaging of rabies virus entry into living verocells”. In: Scientific Reports 5 (2015), p. 11753.86[48] Zamudio-Meza, H., Castillo-Alvarez, A., Gonzalez-Bonilla, C., andMeza, I. “Cross-talk between Rac1 and Cdc42 GTPases regulates for-mation of filopodia required for dengue virus type-2 entry into HMEC-1 cells”. In: Journal of General Virology 90.12 (2009), pp. 2902–2911.[49] Zmurchok, C., Bhaskar, D., and Edelstein-Keshet, L. “Coupling me-chanical tension and GTPase signaling to generate cell and tissue dy-namics”. In: Physical Biology 15.4 (2018), p. 046004.87Appendix ACode for 1D simulationsThe following program written in Matlab is used in Ch. 3 for running 1Dsimulations. The mathematical background is discussed in Sec. 3.1. Partsof the code for setting preferences and saving output has been removed forsimplicity. A more complete version of the code will be published on Github:https://github.com/liuyue002/Wave-pinning-model.Program A.1 Solver for the combined model(2.8) in 1D%% Numerics parametersdrawperframe =50; % update the animation once per 50 i t e r a t i o n snx=400; % number o f s p a t i a l d i s c r e t i z a t i o n po in t sdx=1/nx ;T=1; % end time f o r the s imu la t i ondt =0.0002; % time d i s c r e t i z a t i o nnt=T/dt +1;nFrame=c e i l ( (T/dt )/ drawperframe ) ;%% Model parametersLL=20ˆ2; d e l t a =0.01; k0=1.0∗LL ; gamma=18.5∗LL ; n=2; eta =5.2∗LL ;c =1; theta =5.5∗LL ; alpha =1.5∗LL ;s=40∗LL ; e p s i l o n =0.1 ; kn=24∗LL ; ks =7.5∗LL ;f = @(u , v ,F) ( k0+gamma.∗u . ˆ n . / ( 1 + u . ˆ n ) ) . ∗ v − ( eta+s ∗F. / ( 1 + F ) ) . ∗ u ;% the unique equ i l i b r i umu0=alpha / theta ;f 0=u0∗kn/ ks ;v0=(c∗ alpha+(eta+s ∗ f 0 /(1+ f0 ) )∗ u0 )/ ( k0+gamma∗u0ˆn/(1+u0ˆn ) ) ;x=l i n s p a c e (0 , 1 , nx ) ’ ;u=ze ro s (nx , 1 ) ; v=ze ro s (nx , 1 ) ; F=ze ro s (nx , 1 ) ;88Appendix A. Code for 1D simulationso=ones (nx , 1 ) ;A=spd iags ( [ o −2∗o o ] , [−1 0 1 ] , nx , nx ) ;A(1 ,1)=−1; % f o r no−f l u x BCA(nx , nx)=−1;A=A/( dx ˆ 2 ) ; % A i s the d i s c r e t e Laplac ian%% i n i t i a l c o n d i t i o n su=u0 ∗ ( ( rand ( s i z e (u ) ) ) ∗ 0 . 2 + 0 . 9 ) ;v ( : )= v0 ;F( : )= f0 ;%% Set up f i g u r ef i g=f i g u r e ( ) ;x l a b e l ( ’ x ’ ) ; y l a b e l ( ’u , v ,F ’ ) ;a x i s ( [ 0 1 0 6 ] ) ;uu=ze ro s (nFrame , nx ) ; %f o r p l o t t i n g kymographsvv=ze ro s (nFrame , nx ) ;f f=ze ro s (nFrame , nx ) ;%% Main s o l v e rth =0.5 ; % 0 : forward Euler , 0 . 5 : Crank−Nicolson , 1 : backward EulerTu=speye ( nx)−th∗dt∗ d e l t a ∗A;Tv=speye ( nx)−th∗dt∗A;f o r t i =1:1 : nti f (mod( t i , drawperframe ) == 1)c l a ;hold on ;p l o t (x , u ) ; p l o t (x , v ) ; p l o t (x ,F ) ;l egend ( ’u ’ , ’ v ’ , ’F ’ ) ;hold o f f ;drawnowuu( iFrame , : )= u ;vv ( iFrame , : )= v ;f f ( iFrame , : )=F ;endurhs = u + dt ∗( f (u , v ,F) + (1−th )∗ d e l t a ∗A∗u − c∗ theta ∗u ) ;unew = Tu\urhs ;89Appendix A. Code for 1D simulationsvrhs = v + dt∗(− f (u , v ,F) + (1−th )∗A∗v + c∗ alpha ) ;vnew = Tv\vrhs ;Fnew = F +(dt∗ e p s i l o n )∗ ( kn∗u − ks∗F ) ;F=Fnew ; u=unew ; v=vnew ;endf i g u=plot kymograph (uu , nFrame , nx ,T) ;f i g v=plot kymograph ( vv , nFrame , nx ,T) ;f i gF=plot kymograph ( f f , nFrame , nx ,T) ;The next program is for producing the kymographs.Program A.2 plot kymograph.mfunc t i on [ f i g ] = plot kymograph (uu , nFrame , nx ,T)tTick = ( 0 : 0 . 2 : 1 )∗ nFrame ;tTickLabel=c e l l f u n ( @num2str , num2cel l ( ( 0 : 0 . 2 : 1 ) ∗T) , ’ un ’ , 0 ) ;xTick = ( 0 : 0 . 2 : 1 )∗ nx ;xTickLabel=c e l l f u n ( @num2str , num2cel l ( ( 0 : 0 . 2 : 1 ) ) , ’ un ’ , 0 ) ;f i g=f i g u r e ( ) ;a x i s ( [ 0 nFrame 0 nx ] ) ;% Ignore the i n i t i a l t r a n s i e n t c o n d i t i o n s f o r c a l c u l a t i n g c o l o r rangeumax=max(max(uu ( 1 0 : end , : ) ) ) ; umax=c e i l (umax∗10)/10 ;umin=min ( min (uu ( 1 0 : end , : ) ) ) ; umin=f l o o r ( umin ∗10)/10 ;%umin =0.32;umax=0.35; %tmpu c o l o r t i c k =[umin , umax ] ;imagesc (uu ’ , u c o l o r t i c k ) ;s e t ( gca , ’ YDir ’ , ’ normal ’ ) ;c o l o rba r ( ’ FontSize ’ ,40 , ’ TickLabels ’ , u c o l o r t i c k , ’ Ticks ’ , u c o l o r t i c k ) ;s e t ( gca , ’ XTick ’ , tTick ) ; s e t ( gca , ’ XTickLabel ’ , tTickLabel ) ;s e t ( gca , ’ YTick ’ , xTick ) ; s e t ( gca , ’ YTickLabel ’ , xTickLabel ) ;x l a b e l ( ’ t ’ ) ; y l a b e l ( ’ x ’ ) ;end90Appendix BXPPAUT programs forbifurcation analysisThe following XPPAUT [10] programs were used to produce some of thebifurcation diagrams in Ch. 5 and Sec. 6.1. Although the combined modelcovers all extensions of the wave pinning model, we should always use massconservation to eliminate degeneracy whenever possible. Otherwise, we willobtain incorrect stability results due to the system’s degeneracy.Matcont [5] was used for some of the other bifurcation diagrams. TheMatcont scripts are fairly long so they are not included here, but they willbe available on Github, along with our utility functions for plotting.Program B.1 LPA for the wave pinning model (5.4), see Sec. 5.2du/dt=F(u,w-u)duL/dt=F(uL,w-u)F(G,GI)= (k+gamma*G^n/(1+G^n))*GI-eta*Gpar w=2,k=1.5,gamma=0.01,n=3,eta=15init u=1,uL=1.5Program B.2 LPA for the actin feedback model (5.10), see Sec. 5.4duL/dt=F(uL,w-u,FactinL)du/dt=F(u,w-u,Factin)dFactin/dt=eps*(kn*u-ks*Factin)dFactinL/dt=eps*(kn*uL-ks*FactinL)F(a,b,c)=(k0+gamma*a^n/(1+a^n))*b-(eta+s*c/(1+c))*apar w=2.5,eps=0.1,kn=24,ks=7.5,k0=6,gamma=30,eta=15,s=0.01,n=3init u=1,uL=1,Factin=5,FactinL=591Appendix B. XPPAUT programs for bifurcation analysisProgram B.3 LPA for the combined model, see Sec. 5.5. Set c = 0 for thenon-conservative model.duL/dt=f(uL,v,FactinL)-c*theta*uLdu/dt=f(u,v,Factin)-c*theta*udv/dt=-f(u,v,Factin)+ c*alphadFactin/dt=eps*(kn*u-ks*Factin)dFactinL/dt=eps*(kn*uL-ks*FactinL)f(G,GI,FACT)=(k0+gamma*G^n/(1+G^n))*GI - (eta+s*FACT/(1+FACT))*Gpar k0=1,gamma=30,eta=5,n=3,theta=5.5,alpha=1.5,c=1par s=0.01,kn=24,ks=7.5,eps=0.1init u=1,uL=1,v=5,Factin=1,FactinL=1Program B.4 The dynamic cell size model (6.5). Expansion case shown,contraction case is similar.du/dt=(b+gamma*(u^n/(1+u^n) ))*(w/L-u)-eta*u+u*k*(L-(Lb+Ldiff*(u^n /(uc^n+u^n))))/LdL/dt=-k*(L-(Lb+Ldiff*(u^n/(uc^n+u^n))))par n=8,k=5,w=40,gamma=9,eta=10,uc=0.85,Ldiff=11,b=1,Lb=5.5init u=1,L=192Appendix CCode for 2D simulationsThe following program written in Python was used to run 2D simulationsin Sec. 3.8. We make use of the FEniCS package [2], which is a well-knownlibrary for finite element methods. The output can be viewed with othersoftwares, such as Paraview.This program is built upon the code provided by Dr. Timm Treskatisas a part of the teaching material for MATH 521: Numerical Analysis ofPartial Differential Equations.Program C.1 2D finite element code for the non-conservative model (2.7)from f u t u r e import p r i n t f u n c t i o nfrom f e n i c s import ∗from mshr import ∗p = 1 #order o f i n t e r p o l a t i o nmesh = RectangleMesh ( Point ( 0 . 0 , 0 . 0 ) , Point ( 1 . , 1 . ) , 80 , 80 , ” c ro s s ed ” )V = FunctionSpace (mesh , ”P” ,p)# No−f l u x BC i s automat i ca l l y implemented by d e f a u l t## Numerical methods parameterstheta = 0 .5 # 1 : Backward Euler , 0 : Forward Euler , 0 . 5 : C−RT = 3 0 . ; t = 0 . ; dt = 1e−2; nt = i n t (T/dt )drawPerIter = 5## Model parametersL = 1 . ; Du = 0 . 0 1 ; Dv = 1 . ; k=1.5∗(L∗∗2 ) ; gamma=30.∗(L∗∗2 ) ;eta =15∗(L∗∗2 ) ; c =0. ; t h e t a l o s s =4.5∗(L∗∗2 ) ; alpha =1.5∗(L∗∗2 ) ;## Unique equ i l i b r iumueq=alpha / t h e t a l o s s ;veq=(c∗ alpha+eta ∗ueq )/ ( k+gamma∗( ueq∗∗2 / (1.+ ueq ∗ ∗ 2 ) ) ) ;## I n i t i a l c o n d i t i o n sv o ld = i n t e r p o l a t e ( Constant ( veq ) , V)u o ld = Express ion ( ’ x [ 0 ] < 0 .1 ? 4 .∗ ueq : ueq ’ , degree=p , ueq=ueq )93Appendix C. Code for 2D simulationsu o ld = i n t e r p o l a t e ( u old , V)## Assemble the l i n e a r systemu = Tria lFunct ion (V)v = Tria lFunct ion (V)z = TestFunction (V)B1 = (u∗z + theta ∗dt∗Du∗dot ( grad (u ) , grad ( z ) ) )∗ dxB2 = ( v∗z + theta ∗dt∗Dv∗dot ( grad ( v ) , grad ( z ) ) )∗ dxu = Function (V, name=’u ’ ) ; u . a s s i g n ( u o ld )v = Function (V, name=’ v ’ ) ; v . a s s i g n ( v o ld )fo ldername = ’ OutputFolder ’r e s u l t s u = F i l e ( fo ldername+’ /CN u . pvd ’ )r e s u l t s u << (u , t )r e s u l t s v = F i l e ( fo ldername+’ /CN v . pvd ’ )r e s u l t s v << (v , t )sourceterm = i n t e r p o l a t e ( Constant ( c∗ alpha ) , V)f o r k in range ( nt ) :t += dtf = ( k+gamma∗( u o ld ∗∗2 ) / (1 . + u o ld ∗∗2))∗ v o ld − eta ∗ u o ldRHS = ( u o ld ∗z − (1.− theta )∗Du∗dt∗dot ( grad ( u o ld ) , grad ( z ) ) + \dt∗ f ∗z − dt∗c∗ t h e t a l o s s ∗ u o ld ∗z )∗dxs o l v e (B1 == RHS, u)RHS = ( v o ld ∗z − (1.− theta )∗Dv∗dt∗dot ( grad ( v o ld ) , grad ( z ) ) − \dt∗ f ∗z + dt∗ sourceterm ∗z )∗dxs o l v e (B2 == RHS, v )i f ( k+1) % drawPerIter == 0 :r e s u l t s u << (u , t )r e s u l t s v << (v , t )u o ld . a s s i g n (u)v o ld . a s s i g n ( v )94
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Analysis of pattern formation in reaction-diffusion models for cell polarization
Text