From Signalling to Cell Behaviour: Modelling Multi-ScaleOrganization in Single and Collective Cellular SystemsbyCole ZmurchokB.Sc., University of Alberta, 2012M.Sc., University of Alberta, 2014a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoral studies(Mathematics)The University of British Columbia(Vancouver)July 2018c© Cole Zmurchok, 2018The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the dissertation entitled:From Signalling to Cell Behaviour: Modelling Multi-Scale Organization in Single andCollective Cellular Systemssubmitted by Cole Zmurchok in partial fulfillment of the requirements of the degree ofDoctor of Philosophy in MathematicsExamining Committee:Leah Edelstein-Keshet, Mathematics (Supervisor)Eric Cytrynbaum, Mathematics (Supervisory Committee Member)Andrew Rechnitzer, Mathematics (University Examiner)Geoffrey Wasteneys, Botany (University Examiner)Additional Supervisory Committee Members:James J. Feng, Mathematics (Supervisory Committee Member)iiAbstractIndividually and collectively, cells are organized systems with many interacting parts. Math-ematical models allow us to infer behaviour at one level of organization from information atanother level. In this thesis, I explore two biological questions that are answered throughthe development of new mathematical approaches and novel models.(1) Molecular motors are responsible for transporting material along molecular tracks(microtubules) in cells. Typically, transport is described by a system of reaction-advection-diffusion partial differential equations (PDEs). Recently, quasi-steady-state (QSS) methodshave been applied to models with linear reactions to approximate the behaviour of the PDEsystem. To understand how nonlinear reactions affect the overall transport process at thecellular level, I extend the QSS approach to certain nonlinear reaction models, reducingthe full PDE system to a single nonlinear PDE. I find that the approximating PDE is aconservation law for the total density of motors within the cell, with effective diffusion andvelocity that depend nonlinearly on the motor densities and model parameters. Cell-scalepredictions about the organization and distribution of motors can be drawn from theseeffective parameters.(2) Rho GTPases are a family of protein regulators that modulate cell shape and forcesexerted by cells. Meanwhile, cells sense forces such as tension. The implications of thistwo-way feedback on cell behaviour is of interest to biologists. I explore this question bydeveloping a simple mathematical model for GTPase signalling and cell mechanics. Themodel explains a spectrum of behaviours, including relaxed or contracted cells and cellsthat oscillate between these extremes. Through bifurcation analysis, I find that changes insingle cell behaviour can be explained by the strength of feedback from tension to signalling.When such model cells are connected to one another in a row or in a 2D sheet, waves ofcontraction/relaxation propagate through the tissue. Model predictions are qualitativelyconsistent with developmental-biology observations such as the volume fluctuations in acellular monolayer. The model suggests a mechanism for the organization of tissue-scalebehaviours from signalling and mechanics, which could be extended to specific experimentalsystems.iiiLay SummaryCells are complex systems with many interacting parts. Multi-scale mathematical modelsare used to explore how organization emerges from constituent parts. I answer two biologicalquestions by developing a new approximation method and by formulating a novel model.In order to distribute cellular cargo, cells employ proteins called molecular motors fortransport. The interactions and movement of motors within a cell is described using asystem of partial differential equations. To more easily understand and determine the effectof nonlinear interactions on the overall transport of cargo, I develop a new approximationmethod to reduce the system of equations to a single equation.Protein signalling is responsible for controlling cell size, and can also be affected bymechanical tension experienced by a cell. I develop a new model incorporating the two-way feedback between cell signalling to cell mechanics to explore and explain single andcollective cell behaviours observed in experiments.ivPrefaceChapter 1 and Chapter 4 were written by the author, Cole Zmurchok, for the purposes ofthis thesis.Chapter 2 and Appendix A are based on the published work by Zmurchok, Small, Ward,and Edelstein-Keshet [106], of which I am the first author. I supervised Tim Small, whowas an undergraduate research assistant involved in the computational and mathematicalanalysis in the early stages of this project. Michael Ward and Leah Edelstein-Keshet arejoint supervisors for this project, and assisted in writing parts of the paper. I performedthe research and wrote the paper under their supervision. Michael Ward suggested thenumerical shooting scheme and used this method to produce figures, and performed theboundary layer analysis.Chapter 3 and Appendix B is based on published work by Zmurchok, Bhaskar, andEdelstein-Keshet [107], of which I am the first author. Dhananjay Bhaskar was involvedin the 2D implementation of the model using the cellular Potts framework (CPM), andis responsible for producing the figures using the CPM. I collaborated with DhananjayBhaskar throughout this process, and we worked closely to adapt, implement, and codethe 2D implementation (including the Kuramoto order parameter study). I supervisedJim Shaw, an undergraduate research volunteer, who was responsible for performing anadditional numerical bifurcation analysis of the Rac-Rho-tension model. Leah Edelstein-Keshet supervised the research and assisted in writing the paper.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Multi-scale Modelling in Cellular Systems . . . . . . . . . . . . . . . . . . . 11.2 Intracellular Transport by Molecular Motors . . . . . . . . . . . . . . . . . . 11.2.1 A Primer on Quasi-steady-state Methods . . . . . . . . . . . . . . . 41.3 The Interplay Between Cell Signalling and Cell Mechanics . . . . . . . . . . 81.3.1 Dynamical Systems for Cell Signalling and Cell Mechanics . . . . . . 101.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Application of Quasi-Steady-State Methods to Nonlinear Models of In-tracellular Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.1 Intracellular Transport by Molecular Motors . . . . . . . . . . . . . . . . . . 172.2 Models of Intracellular Transport . . . . . . . . . . . . . . . . . . . . . . . . 212.2.1 Kinesin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.2.2 Kinesin-Dynein Model . . . . . . . . . . . . . . . . . . . . . . . . . . 242.2.3 Myosin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26vi2.3 Quasi-steady-state Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4 Examples of the QSS Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 312.4.1 QSS Reduction: Kinesin Model . . . . . . . . . . . . . . . . . . . . . 322.4.2 QSS Reduction: Kinesin-Dynein Model . . . . . . . . . . . . . . . . 422.4.3 QSS Reduction: Myosin Model . . . . . . . . . . . . . . . . . . . . . 452.5 Boundary Layer Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 552.5.1 The Kinesin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 592.5.2 The Kinesin-Dynein Model . . . . . . . . . . . . . . . . . . . . . . . 592.5.3 The Myosin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 622.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633 Coupling Mechanical Tension and GTPase Signalling to Generate Celland Tissue Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 683.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 683.2 Minimal Model for a Single Mechanochemical Cell . . . . . . . . . . . . . . 713.2.1 Model Equations and Definitions . . . . . . . . . . . . . . . . . . . . 713.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.3 Mechanical Coupling in a 1D Array of Cells . . . . . . . . . . . . . . . . . . 763.3.1 Tissue Dynamics in 1D Depend on Mechanical Feedback Strength . 763.4 Cell shape and cell-cell interactions in 2D epithelial sheets . . . . . . . . . . 793.4.1 Adapting the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 793.4.2 Single Cell Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 803.4.3 Coupling CPM Cells in 2D . . . . . . . . . . . . . . . . . . . . . . . 823.4.4 Waves of Contraction and GTPase Activities in 2D Model Tissue . . 833.5 Rac and Rho GTPase Model . . . . . . . . . . . . . . . . . . . . . . . . . . 843.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 894 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97A Supporting Materials for the Application of Quasi-Steady-State Meth-ods to Nonlinear Models of Intracellular Transport . . . . . . . . . . . . 107A.1 Convergence of the Kinesin Model to the QSS for ε→ 0 . . . . . . . . . . . 107A.2 Numerical Methods for the QSS . . . . . . . . . . . . . . . . . . . . . . . . . 108A.3 Microtubule Density and Binding by Motor Complexes . . . . . . . . . . . . 110A.3.1 Kinesin Model with Nonuniform MT Density . . . . . . . . . . . . . 110A.3.2 Kinesin-Dynein Model and the Function Q(x) . . . . . . . . . . . . . 110viiA.4 Scaling the QSS Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111A.4.1 The Kinesin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111A.4.2 Kinesin-Dynein Model Scaling . . . . . . . . . . . . . . . . . . . . . 114A.4.3 Myosin Model Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . 115A.4.4 Non-spatial Myosin Model . . . . . . . . . . . . . . . . . . . . . . . . 116B Supporting Materials for Coupling Mechanical Tension and GTPaseSignalling to Generate Cell and Tissue Dynamics . . . . . . . . . . . . . 118B.1 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118B.2 Scaling the GTPase Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118B.3 1D Methods: Multicellular Simulations . . . . . . . . . . . . . . . . . . . . . 120B.4 2D Methods: Cellular Potts Model . . . . . . . . . . . . . . . . . . . . . . . 120B.5 2D Methods: Patch Size and Synchronization . . . . . . . . . . . . . . . . . 122B.6 2D Methods: Large-Tissue Simulations . . . . . . . . . . . . . . . . . . . . . 122B.7 Additional 2D Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122viiiList of FiguresFigure 1.1 Motor based intracellular transport in fungi. . . . . . . . . . . . . . . . 3Figure 1.2 Effective diffusion in a system with membrane and cytosolic proteins. . 6Figure 1.3 GTPases act as switches. . . . . . . . . . . . . . . . . . . . . . . . . . . 9Figure 1.4 Mechanochemical coupling between GTPase signalling and cell mechanics. 10Figure 1.5 Mechanical model from Odell et al. [60]. . . . . . . . . . . . . . . . . . . 12Figure 1.6 Phase-plane for the Odell-Oster Model. . . . . . . . . . . . . . . . . . . 14Figure 1.7 Phase-plane for the GTPase-tension Model. . . . . . . . . . . . . . . . . 14Figure 1.8 A schematic of a typical CPM. . . . . . . . . . . . . . . . . . . . . . . . 16Figure 2.1 A schematic diagram of kinesin-based intracellular transport in a 1D cellof length L0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 2.2 A schematic diagram of kinesin-dynein motor-based intracellular trans-port in a 1D cell of length L0. . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.3 Geometry of the QSS approximation. . . . . . . . . . . . . . . . . . . . 31Figure 2.4 Effect of nonlinear binding and microtubule polarity. . . . . . . . . . . . 37Figure 2.5 Effect of the relative magnitudes of binding and unbinding rates kb, ku. 38Figure 2.6 Effects of two spatially dependent Microtubule (MT) bias functions, P (x). 40Figure 2.7 Effect of the Hill function parameters K and n. . . . . . . . . . . . . . . 41Figure 2.8 Comparison of the full solution with the QSS solution. . . . . . . . . . . 43Figure 2.9 The effect of parameters Q, v, ka, k on the total density. . . . . . . . . . 44Figure 2.10 Full numerical vs. asymptotic solutions to the myosin model. . . . . . . 48Figure 2.11 Region of solution existence (unshaded). . . . . . . . . . . . . . . . . . . 50Figure 2.12 Effect of the (scaled) binding rate, kb. . . . . . . . . . . . . . . . . . . . 51Figure 2.13 Effect of the (scaled) stalling rate, kbw. . . . . . . . . . . . . . . . . . . 51Figure 2.14 Effect of treadmilling speed, v. . . . . . . . . . . . . . . . . . . . . . . . 52Figure 2.15 pB(x, t) converges to Type I QSS. . . . . . . . . . . . . . . . . . . . . . . 53Figure 2.16 Steady-state behaviour of the regularized myosin model. . . . . . . . . . 54Figure 2.17 Myosin model initial condition dependence. . . . . . . . . . . . . . . . . 55ixFigure 2.18 Qualitative analysis of boundary layer behaviour of the kinesin-dyneinmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 2.19 Qualitative analysis of boundary layer behaviour of the myosin model. . 63Figure 3.1 The minimal model for coupled GTPase activity and cellular-tension. . 72Figure 3.2 Bifurcation diagrams for the minimal model . . . . . . . . . . . . . . . . 74Figure 3.3 Dynamics of the minimal model for a single cell with one GTPase andfeedback from tension. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 3.4 Cell interactions in a 1D array of “model cells”. . . . . . . . . . . . . . . 76Figure 3.5 1D tissue dynamics result from mechanochemical interactions. . . . . . 78Figure 3.6 Single cell oscillations in 2D cells simulated with the Cellular Potts Model(CPM). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Figure 3.7 Simulation of a 2D tissue with intermediate adhesion. . . . . . . . . . . 85Figure 3.8 Simulation of a 2D tissue with strong adhesion. . . . . . . . . . . . . . . 86Figure 3.9 Simulation of a 2D tissue with randomly chosen β in each cell. . . . . . 87Figure 3.10 Bifurcation diagrams for the minimal Rac-Rho-tension model. . . . . . . 88Figure 3.11 Dynamics of the single-cell Rac-Rho-tension model. . . . . . . . . . . . 89Figure A.1 Convergence of the QSS approximation to the full model for ε→ 0. . . 108Figure A.2 Phase-plane analysis of the non-spatial myosin model. . . . . . . . . . . 117Figure B.1 Additional 1D tissue dynamics. . . . . . . . . . . . . . . . . . . . . . . . 121Figure B.2 Initial conditions for 9 independent oscillatory cells. . . . . . . . . . . . 123Figure B.3 Synchronization does not appear in the uncoupled oscillating cells. . . . 124Figure B.4 Initial conditions for 9 mechanically coupled oscillatory cells. . . . . . . 125Figure B.5 Synchronization in the low adhesion regime. . . . . . . . . . . . . . . . . 126Figure B.6 Synchronization in the intermediate adhesion regime. . . . . . . . . . . . 127Figure B.7 Synchronization in the high adhesion regime. . . . . . . . . . . . . . . . 128Figure B.8 A single relaxed cell in 2D. . . . . . . . . . . . . . . . . . . . . . . . . . 128Figure B.9 A single cell in 2D can exhibit damped oscillations. . . . . . . . . . . . . 129Figure B.10 A single cell in 2D can enter a small amplitude limit cycle. . . . . . . . 129Figure B.11 Stochastic switching between the low amplitude limit cycle and the highamplitude limit cycle for a single cell in 2D. . . . . . . . . . . . . . . . . 129Figure B.12 Simulation of a 2D tissue with weak adhesion. . . . . . . . . . . . . . . 130xGlossaryBC boundary conditionCPM Cellular Potts modelFP Fokker-PlanckIVP initial value problemIBVP initial-boundary value problemMCS Monte-Carlo stepMT MicrotubuleODE ordinary differential equationPDE partial differential equationQSS quasi-steady-statexiAcknowledgmentsI would like to thank and acknowledge Leah Edelstein-Keshet for her continued mentorshipand constant encouragement. Thank you, Leah. I am certain that your lessons and guidancewill have an even greater and lasting impact in the coming years.I would like to thank the many members of the UBC math biology, applied math,and math education community. I would especially like to thank Michael Ward for hismentorship and guidance, and Eric Cytrynbaum and Jimmy Feng for their mentorshipand help as members of my supervisory committee. I would also like to acknowledge themembers of the Keshet-Feng research group for their questions, suggestions, and for weeklygroup meetings.xiifor AnniexiiiChapter 1Introduction1.1 Multi-scale Modelling in Cellular SystemsThe broad question that this thesis aims to address is how biological mechanisms combineto organize cell behaviours. Biological experiments and observations, allow for us to under-stand the mechanisms, facts, and theories about how various proteins, molecular motors,cell signalling molecules, and biophysical structures interact with each other. Multi-scalemodelling and analysis from applied mathematics can serve as a tool for understanding howcell-scale or tissue-scale organization emerges from the interactions of the many differentconstituent parts.This thesis is divided into two self-contained parts, Chapter 2 and Chapter 3. Each partdescribes a different example of multi-scale modelling in cell biology and utilizes differentmathematical approaches. In this introduction, I provide context and background for eachpart with respect to the broader field, and preview the main questions, methods, and results.I end the introduction with a brief description of the thesis contents.1.2 Intracellular Transport by Molecular MotorsIn many cellular scenarios, diffusion alone is insufficient to transport cargo over the distancesrequired to maintain cellular function. This is especially important when cargo needs to betransported over long distances, such as in neurons or fungal cells (Figure 1.1(a)). For aneuron to function, it is necessary to transport cargo over a long distance from the cell bodydown the axon to the synapses. In fungi, it is necessary to transport cargo to maintain andexpand the cell wall during growth. Figure 1.1(b)–(c) depicts the distribution of cargo inboth yeast-like and hyphal fungi cells in green against the cell membrane in red. In bothneurons and fungi, the transport of cargo is mediated by proteins called molecular motors.Molecular motors are proteins that utilize energy in order to transport cargo by “walking”1along cellular “tracks”. Molecular motors typically transport cargo in vesicles and movealong protein filaments of the cytoskeleton (as in Figure 1.1(d)). Cytoskeletal filamentssuch as microtubules (MTs) or actin filaments give a cell structure and provide roads forthe molecular motors to use [10]. Microtubules and actin filaments are asymmetric filamentswith distinct “plus” and “minus” ends.Different types of molecular motors such as kinesin, dynein, and myosin motors havedifferent biophysical properties. For example, kinesin motors, which largely walk towardsthe plus ends of polarized MTs can also exist as unbound, cytosolic forms [5]. Dyneinmotors walk towards the minus ends of MTs. Motors may also be present in complexes,where more than one motor of any type can be bound to cargo. Myosin motors walk onactin filaments instead of MTs, and are mostly associated with muscle contraction. Somemyosin motors are responsible for transporting cellular cargo [86, 101]. The interactionsof many molecular motors, the protein filaments that they move along, and other cellularfactors result in cell-scale distribution of molecular motors (and/or cargo).Nonetheless, it is difficult to understand how a variation in a particular biophysical pa-rameter will affect the cell-scale motor distribution in such a complicated system, more soif nonlinear effects are important. Many mathematical approaches have had success un-derstanding the transport phenomena at different scales, e.g., several motors attached toa single cargo [4, 31, 40, 48, 55] or at the cell-scale using stochastic or partial differentialequation approaches [11, 66, 72, 85]. A system of reaction-advection-diffusion partial dif-ferential equation (PDE) is often used to quantitatively describe transport by molecularmotors within cells [16, 26, 58]. To understand the effect of parameters in such a system,I develop an approximation method to reduce the dimensionality of the system and revealhow the effective transport parameters depend on the biophysical parameters of the molecu-lar motors. The approximation method relies on existence of two separate time-scales in thesystem: fast binding, unbinding and other reaction interactions, but slow spatial processesand conservation of molecular motors. Using methods from asymptotic analysis, I finda quasi-steady-state PDE, which acts as an approximation to the full reaction-advection-diffusion system.The quasi-steady-state (QSS) approximation developed here extends the methods ofNewby and Bressloff [58] and Bressloff and Newby [7], who primarily developed and usedthe quasi-steady-state approximation to determine the mean first passage time for molecularmotor intermittent search, i.e., the average time for cargo to be transported to a specificcellular site by motors. From a PDE transport model for early endosome organizationin fungal hyphae [26] with linear state transitions, Dauvergne and Edelstein-Keshet [16]utilized quasi-steady-state approximation methods to determine the effective velocity andeffective diffusion parameters which describe the overall transport process. A limitation of2Figure 1.1: An overview of how molecular motor based intracellular transport supportsgrowth in fungi. In (a), the two main growth modes of fungi are shown: yeast-like cellsand filamentous fungi. In (b)–(c), the cell membrane is shown in red, while myosin motorstransporting material necessary for the construction of the cell wall is shown in green. In(d), the roles of molecular motors and cytoskeleton are illustrated. Molecular motors areresponsible for transporting vesicles containing cargo. Reprinted from Current Opinionin Microbiology, 14, Gero Steinberg, Motors in fungal morphogenesis: cooperation versuscompetition, 660–667, Copyright 2011, with permission from Elsevier.this approach is that the methodology only applies to linear state transitions. Nonlinearinteractions may be better motivated in some biological situations. For example, rates ofmotor binding to microtubules could be limited by competition for binding sites, or couldbe cooperative. Mass-action kinetics or other nonlinear interactions, such as “traffic-jam”cubic nonlinearities could play an important role in the organization of the transport process[102]. Here, I extend the quasi-steady-state methods from Newby and Bressloff [58] andDauvergne and Edelstein-Keshet [16] to a class of models with nonlinear reaction kinetics.I apply the QSS methodology to three different transport systems. These models consist3of (1) a model for kinesin motors with a saturating binding rate, (2) a model for kinesin-dynein-cargo complexes whose interactions on a microtubule can change the direction oftransport, and (3) a model for unconventional myosin motors that stop moving (stall) uponencountering other stalled motors. Although the motors move at known speeds, I showthat their interactions with MTs and each other lead to slower overall effective “speed”(transport rate) as well as spread (similar to diffusion) in their spatial distribution. Thebiological contribution is that I am able to relate the effective rate of transport and diffusionto the details of binding, MT polarity and nonlinear motor interactions. Mathematically,this manifests through a scalar conservation law for the density of motors within the cell,parametrized in terms of the density of motors in one of the states, e.g., the freely diffusingstate. In this so-called “QSS PDE”, the flux determines an effective velocity and effectivediffusion coefficient that depend nonlinearly on the biophysical parameters and the motordensity. From the effective velocity and effective diffusion, the overall cell-scale transportcan be understood as a function of the motor-scale model parameters.In the next section, I provide a short mathematical primer on the ideas behind thequasi-steady-state approximation.1.2.1 A Primer on Quasi-steady-state MethodsThe main mathematical idea behind the quasi-steady-state approximation is the notion ofseparation of time-scales. In a dynamical system, the time-scale associated with one variableor process may be slower than that of others. Consider, for example, the following ordinarydifferential equation (ODE) system for two species u(t) and v(t):dudt= f(u, v), (1.1)dvdt=1ε(v − a), (1.2)where f is some function describing the rate of change of u, a is the steady-state of speciesv, and 0 < ε 1 is a small parameter so that 1ε 1 is large. Provided that the function fis not small, i.e., f = O(1), then u varies on a time-scale much faster than v. The detailsof the contribution of various terms in the dynamical system should be obtained throughcareful non-dimensionalization. As such, the quasi-steady-state approximation isv ≈ a (1.3)4since v rapidly relaxes to its steady-state a. From this approximation, it is possible to studyonly the reduced system, which has fewer dynamical variables:dudt= f(u, a). (1.4)Theoretically, the parameter dependence of solutions u(t) on a could be more easily deduced.This is the main idea behind the derivation of classical Michaelis-Menten kinetics (see [83]for a careful case study or Chapter 7 of [22] for an overview).To illustrate the utility of the QSS approximation in a cellular context and to explain theutility of the approximation, I will discuss a variant of the QSS approximation from Mare´eet al. [49]. In this work, the authors use a multi-scale modelling approach to understandcell polarization and movement. In part of the model, regulatory proteins are found tobind and unbind from the cytosol onto the cell membrane, and the diffusion coefficients forproteins in either state are different—membrane bound proteins diffuse more slowly thanin the cytosol. Owing to the fact that switching between the states is rapid, the authorsuse a QSS approximation to describe the effective diffusion of the protein. I review part oftheir approximation below.Let GC(x, t) and GM(x, t) be the amounts of regulatory proteins in the cytosolic andmembrane-bound states, respectively, and G(x, t) = GC(x, t) + GM(x, t) the amount atspatial location x. Suppose that the cytosolic proteins bind to the membrane with ratecoefficient kon and unbind with rate coefficient koff, and that the two classes of proteinshave different diffusion coefficients, DC and DM, respectively. This scenario is illustratedin Figure 1.2.For simplicity, consider a 1D geometry, i.e., x ∈ R, with no-flux boundary conditionsat either end of the cell. In this case, the dynamics of the proteins can be captured by thefollowing reaction-diffusion system of partial differential equations:∂GM∂t= DM∂2GM∂x2+ konGC − koffGM, (1.5)∂GC∂t= DC∂2GC∂x2− konGC + koffGM. (1.6)Note that the total amount of protein within the entire cell is conserved by these equations,given no-flux boundary conditions, and that the proteins can only transition between themembrane bound and cytosolic states. This is a similarity which the models for molecularmotor transport in Chapter 2 share. As this section is a “primer”, I will not proceed formallyand introduce a small parameter ε into this problem; instead, only focus on the idea. Owingto the separation of time-scales (slow spatial processes and rapid binding and unbinding),it is possible to think of the terms in the PDE system (1.5) contributing on different time-5cell membraneGCcytosolkon koffGMDCDMFigure 1.2: Some proteins convert between active, membrane-bound and inactive, freelydiffusing cytosolic states (GM and GC, respectively). Proteins bind to the membrane andunbind from the membrane with rate coefficients koff and kon respectively. The diffusioncoefficient DM in the cell membrane is typically smaller than in the cytosol DC. The QSSreduction method can be applied to a system of this type to derive an effective diffusion,which describes the diffusion of the total amount of protein.scales. First, at any given spatial location, the amount of cytosolic and membrane-boundproteins will equilibrate on a short time-scale with essentially no diffusion. Second, on aslower time-scale, the diffusion will become significant. As such, it is helpful to think of ashort time-scale denoted by τ (versus t), where the protein dynamics are dominated by thereaction-terms:dGMdτ= konGC − koffGM, (1.7)dGCdτ= −konGC + koffGM. (1.8)The fact that these equations operate on a fast time-scale means that the variables rapidlyequilibrate. On the short time-scale τ , this motivates the search for the steady-state of thissystem of ODEs. The steady-state of this ODE system is given by the linear relationshipGM =konkoffGC, hence the amount of membrane-bound proteins can be calculated from theamount of cytosolic proteins and vice versa. From this, I can make a quasi-steady-stateapproximation in the full system:GM(x, t) ≈ konkoffGC(x, t). (1.9)Using this approximation, I find that G(x, t) = kon+koffkoff GC(x, t) =kon+koffkonGM(x, t) since6G(x, t) = GM(x, t)+GC(x, t). For the same reason, note that adding the PDE in (1.5) gives∂G∂t= DM∂2GM∂x2+DC∂2GC∂x2. (1.10)Using the quasi-steady-state approximation G = kon+koffkoff GC =kon+koffkonGM in the PDE forG, I find that G satisfies a diffusion equation:∂G∂t= DMkonkon + koff∂2G∂x2+DCkoffkon + koff∂2G∂x2(1.11)=(DMkonkon + koff+DCkoffkon + koff)∂2G∂x2(1.12)= Deff∂2G∂x2, (1.13)where the effective diffusion coefficient, Deff is a combination of the reaction rate coefficients,kon and koff, and the diffusion coefficients of proteins in either states, DM and DC. Increasingkoff, for example, would increase the contribution of the diffusion in the cytosolic state, tothe overall distribution of proteins. In this case, the effective diffusion coefficient can alsobe interpreted as a weighted average of the original diffusion coefficients with the weightsas the mean fraction residence time in either state. To see this, identify τM =1koffas themean residence time on the membrane (before unbinding with rate koff) and τC =1konas themean residence time in the cytosol (before binding to the membrane with rate kon). Withthese definitions, note thatDeff =konkoffkon + koff(DMkoff+DCkon)(1.14)=1τM + τC(DMτM +DCτC) (1.15)= DMτMτM + τC+DCτCτM + τC. (1.16)This calculation shows that the effective diffusion is a combination of diffusion on themembrane and diffusion in the cytosol, with the relative importance of each weighted bythe fraction of time spent in the given state. The QSS not only reduces the complexity ofthe model but also provides insight into the biophysical processes.In Chapter 2, the QSS approximation is applied to a model similar to (1.5); however,instead of a reaction-diffusion model, I will extend the QSS methodology to a class ofthree-component reaction-advection-diffusion system with nonlinear state transitions. Inthis case, I will show that the quasi-steady-state approximation is a PDE written in termsof one of the molecular motor states, denoted by α(x, t). The QSS PDE will be a balance7equation for the total amount of motors in the cell, with effective diffusion and effectivetransport terms that depend on the model parameters and on the amount of motors in thereference state α:∂∂ty(α) = − ∂∂x(V(α)−D(α)∂α∂x). (1.17)In contrast with previous work, the dependence on α can be nonlinear. Here, y(α) is thetotal density of the motors in the cell, V(α) is the effective transport term, and D(α) isthe effective diffusion. The QSS PDE is supplemented with no-flux boundary conditions ateither end of the cell.1.3 The Interplay Between Cell Signalling and CellMechanicsIn the second part of the thesis, Chapter 3, instead of studying how sub-cellular interactionslead to cell-scale organization, I “zoom out” and study how cell-cell interactions can leadto emergent organization and behaviour at the tissue level.The Rho-family GTPase proteins are central regulators within signalling networks of eu-karyotic cells that are largely responsible for coordinating downstream signalling leading tochanges in cell shape mediated by the actin cytoskeleton [73]. GTPases act as switches andexist in either an active (membrane-bound) or inactive (freely-diffusing cytosolic) state, asillustrated in Figure 1.3. When activated, GTPases transmit signals to other proteins thateventually lead to changes in cell behaviour. The activation and deactivation of GTPasesis mediated by other proteins known as GTPase activating proteins (GAPs) and guaninenucleotide exchange factors (GEFs), respectively. There is a large body of literature ad-dressing how GTPases spontaneously segregate to the front or back in a cell and how thiscan lead to cell polarization and movement [25, 54, 59, 64, 95, 99]. Here, I focus on adifferent aspect of GTPase activity and cell behaviour, namely how GTPase activity cancause a cell to contract or spread [3, 14, 78], and the resulting feedback with tension andmechanical forces.Rho GTPases cycle between membrane-bound active forms, and freely-diffusing inactiveforms. Their activation and deactivation is controlled by a large signalling network respon-sible for coordinating cellular responses to stimuli. The Rho-family GTPases Rac1, RhoA,and Cdc42 are central regulators of this network. Rac1 and RhoA (henceforth Rac andRho) are downstream of cell-surface receptors that are sensitive to many different stimuli,including mechanical tension [18, 74, 97]. Moreover, there is increasing evidence for thespecific molecular mechanisms and pathways by which Rac and Rho respond to mechanicalsignals [37, 62]. Of note is the effect that mechanical tension can have on activating or deac-tivating GTPase activity [15, 35]. Additionally, the feedback between mechanical signalling8cell membraneGAPs GEFsinactive GTPaseactive GTPasecytosoleffectorssignalsFigure 1.3: GTPases act as switches. They exist in a membrane-bound active or a freelydiffusing inactive state. Activation and deactivation of GTPases are regulated by otherproteins (GAPs and GEFs), which are activated by cellular signals. Once active, GTPasessignal to downstream effectors.and GTPase behaviour is two-way. While forces such as mechanical tension can influenceGTPase activity, GTPase activity modifies cell shape through downstream signalling, whichwill subsequently affect the mechanical forces acting on a cell. This change in forces willthen subsequently affect GTPase activity, and so on. It is this interplay—between signallingand mechanics—that I seek to understand.Many recent experimental and modelling studies have provided evidence for the ideathat cell behaviours (such as change in cell shape or polarity) can be explained as emergentproperties of small subsets of large signalling networks [3, 9, 14, 78]. Mechanochemicalinteractions have been included in mathematical models of cell behaviour [60, 61, 63], andmodels consisting of GTPase modules can explain cell polarization [54], cell shape [33], andcell migration patterns [34, 65].Using a minimal model for GTPase activity within a cell, I will explore the implicationsof the idea that mechanical tension on the cell is responsible for modifying the GTPaseactivation rate on cell behaviour. Specifically, I assume that GTPase can cycle betweeninactive and active forms, with some positive feedback upon activation. This is illustratedin Figure 1.4 with the active GTPase denoted by G and inactive by Gi. I also assume thatactive GTPase is responsible for contraction in the cell, and that mechanical tension on thecell is responsible for increasing the rate of activation of GTPase (purple boxes in Figure 1.4).Active GTPase leads to cell contraction, while tension is assumed to increase the activationrate of GTPase. If the cell is initially stretched, high tension will lead to the activation ofGTPase, leading to contraction in turn. Once contracted, the initial stretch will no longercontribute to increasing the activation rate of GTPase, and the cell will return to rest. Toimplement this mechanochemical coupling, I use a mechanical model consisting of springsand dashpots to model cell length, and I propose and analyze a minimal GTPase-tensionmodel. The simple two dimensional ODE model describes the dynamics of the GTPase9contractiontensionGiGFigure 1.4: GTPase cycles between active and inactive forms, G and Gi, respectively, withpositive feedback upon the active state. Active GTPase leads to cell contraction, whiletension is assumed to increase the activation rate of GTPase. If the cell is initially stretched,high tension will lead to the activation of GTPase, leading to contraction in turn. Oncecontracted, the initial stretch will no longer contribute to increasing the activation rate ofGTPase, and the cell will return to rest.activity and the length of the cell. Thanks to the model simplicity it is possible to studythe solution behaviour easily using numerical bifurcation analysis (as in Chapter 3), or bystudying the phase plane (as illustrated in the next section). In short, the model exhibitsthree distinct behaviours dependent on the strength of coupling between tension and theGTPase activation rate: long, relaxed cells; short, contracted cells; and cells that oscillatebetween these two extremes. Building on this understanding, I next consider the dynamicsof the minimal model when many cells are coupled together in a 1D array (representing anepithelial sheet), and finally using the Cellular Potts model (CPM) (reviewed in the nextsection) to explore the dynamics in a 2D epithelial tissue. In this way, the behaviour of thecell collective or tissue can be understood as emergent behaviour from cell-cell mechanicalcoupling and mechanochemical signalling in each cell.In the next section, I provide a short mathematical primer for the modelling approachesused in Chapter 3, and draw comparisons to the new model I have developed with themechanochemical model from the 1980s by Odell et al. [60].1.3.1 Dynamical Systems for Cell Signalling and Cell MechanicsIn this section, I provide some preliminary results and explain some of the mathematicalmodels used in Chapter 3. In particular, I will outline (1) a minimal model for GTPaseactivity, (2) compare my results with those of Odell et al. [60], and (3) briefly describe theCellular Potts model (CPM).10A Minimal Model for GTPase ActivityI adapt the minimal GTPase model in Chapter 3 from the modelling work in [33, 34, 54],but ignore spatial effects and only consider the well-mixed model. The basic form of theequation isdGdt= (rate of activation)Gi − (rate of inactivation)G, (1.18)whereG is the amount of active GTPase andGi is the amount of inactive GTPase. Followingprevious work, I assume that the total amount of GTPase in the cell is roughly constantover the timescale of interest, i.e., that GT = G+Gi is constant, and that there is positivefeedback from active GTPase to itself. This leads to an equation of the formdGdt=(b+ γGn1 +Gn)(GT −G)−G, (1.19)where b is the basal rate of activation, and γ gives the amplitude of the positive feedback.The Hill function Gn1+Gn is a saturating function of G. When G is much greater than 1,Gn1+Gn ≈ 1, but when G is much less than 1, Gn1+Gn ≈ 0. Note that time has been scaled sothat the rate coefficient of deactivation of GTPase appears to be 1. This model has therequisite features for GTPase activity: a high-activity and a low-activity steady-state, andis bistable for a range of parameters. Bistability means that the GTPase activity G couldtend toward the high-activity steady-state or low-activity steady-state depending on initialconditions. Moreover, the presence of bistability indicates the possibility of hysteresis, i.e.,history-dependent transitions. Suppose, for example, that a stimuli can increase the GTPaseactivity. If the cell is at the low-activity steady-state, then a small, single stimuli wouldnot automatically cause the cell to jump to the high-activity steady-state. Instead, thanksto the presence of bistability, the small stimuli would disappear as the cell returns to thelow-activity steady-state. Only if the stimuli is sufficient to push the GTPase activity intothe regime of high-activity will the cell transition away from the low-activity steady-state.The presence of bistability and hysteresis can also give rise to a relaxation oscillation, wherethe behaviour of the system slowly transitions through the low or high-activity steady-stateregion before “jumping” to the other (provided the system is coupled to another variable,such as cell length). Indeed, in Chapter 3, I will modify the activation rate to include atension-dependent activation rate (later called f(T )), and the feedback from tension willdrive the system into a low-activity, high-activity, or oscillatory state.The Mechanical Basis of MorphogenesisIn the 1980s, Odell et al. [60] developed a mechanical model for cells in an epitheliumduring morphogenesis (development). The tissue undergoes shape changes as the cells bend11ODELL ET AL. Mechanical Basis of Morphogenesis 447 apical surface basal surface FIG. 1. Schematic of an epithelial cell layer in cross section showing apical junctions. insist that the cells adhere to one another laterally; one can also allow lateral slippage between cells. The me- chanical properties of the cells are the crucial feature of the model, which we now state in the form of hy- potheses. (a) Beneath the apical surface, each cell contains a network of contractile fibers (microfilaments) which are anchored to the plasma membrane at the lateral pe- riphery of the cell (e.g., the zonula adherens, cf. Spooner, (1975)). This is shown schematically in Fig. 2a. For the purposes of the model, we need not specify the precise configuration of these elements; we specify only that their contraction shortens the apical circumference of the cell and reduces the apical surface area (cf. Fig. 2b). Thus, a cell can deform to a trapezoidal cross section as a result of the “purse-string” contraction occurring in the apical cortex. Whether or not apical constriction actually occurs depends on whether the contractile fil- aments generate forces large enough to overcome in- ternal viscous forces and tractions applied externally by neighboring cells. (b) We shall assume that, during the contraction process, the volume of each cell remains essentially con- stant. Thus, the constant volume constraint will induce the cell to elongate basally in response to an apical contraction. We postpone specifying how the contractile fibers and structural reinforcing microtubules are dis- tributed elsewhere in the cell. We suppose only that their distribution leads to cell boundaries that act like passive viscoelastic structures, and bulk interior cyto- plasm that acts like a passive viscoelastic solid (cf. Marsland, 1956; Taylor and Condelis, 1979). (c) The apical contractile filaments constitute an ac- tive (i.e., excitable) system as follows. (i) If an apical fiber is stretched a small amount, by the drawing apart of the apical surface, it acts as an elastic material and when released; contracts back to its original length. (ii) If, however, the apical bundle is stretched beyond a certain point, the contractile system “fires” and an active contraction is triggered. This “rapid” contrac- tions works to shrink the apical bundle as shown in Fig. 2b. The system does not return to its original config- uration; instead, it remains “frozen” in a new, con- tracted state with an apical surface area smaller than before. The viscoelastic properties of the apical filament bun- dle implied by hypothesis (c) are summarized in Fig. 3. The mathematical model underlying Fig. 3 is given in Appendix 1. 2.2. Description of Mechanical Properties The description of the mechanical properties of the microfilament bundle described in Fig. 3 can be sche- matized as the mechanical model shown in Fig. 4. The internal viscoelastic properties of the apical bundle are represented by a dashpot with viscosity cc, and a spring of elasticity k. The elastic restoring force of the spring is a function of the difference between the actual and rest lengths, F elas = -ML - LJ, (1) while the viscous drag force is proportional to the ve- locity of contraction, F. =- dL wse ’ dt . (2) The motion of this mechanical system must obey Newton’s laws: mass X acceleration = sum of forces. That is, md2L=F dt2 elas + Fvisc + Floadr (3) where Fload is the force exerted on the bundle, parallel to the bundle by neighboring systems, and m represents the net mass moved by a length change of the band. A simple argument, sketched in Appendix 2, dem- onstrates the following crucial fact. For virtually all embryological processes we can completely neglect the effects of inertial forces. That is, if we compare the magnitude of the acceleration term with the viscous term by using the dimensionless ratio apical filiment bundle (a) lb) FIG. 2. (a) Network of contractile filaments in the apical region of an epithelial cell. (b) “Purse-string” contraction of the apical circum- ference by the apical bundle. 448 DEVELOPMENTAL BIOLOGY VOLUME 85, I381 we find that R-C 10m5. Therefore, the motion of the fil- load ament bundle model is always such that viscous forces k exactly balance elastic and external loading forces. Thus the equation describing the system is (3) with its -L- left-hand side set to zero and Eqs. (1) and (2) substi- FIG. 4. Schematic of viscoelastic unit used to model a filament tuted. This can be written as bundle. (5) That inertial forces are not a factor in morphogenetic processes is a crucial and little appreciated fact, and leads to some surprising and counterintuitive phenom- ena which will be discussed below. To complete the description of the model we need to describe how the rest length, L+, varies. This can be defined by a differential equation specifying how L,, varies in time when the spring is stretched: L=L, L /- T2 TLa Lo2 LOI FIG. 3. Our postulates concerning the viscoelastic properties of the apical circumferential bundle of contractile filaments. The vertical axis (L) represents the actual instantaneous length of the cell’s apical bundle. The horizontal axis (Lo) represents the equilibrium (rest) length of the bundle; i.e., the length it would assume in the absence of stretching forces. The apical bundle is stress-free only when the actual length equals the equilibrium length: L = &. In the simplest case, we shall assume that the apical bundle has only two stable equilibrium lengths: long (J&) and short (&). By stable, we mean that following a small displacement in the actual length (L), the sys- tem returns to the same equilibrium point. For example, a dilation of the fiber bundle length from 4, to a point a, returns to 4, along a trajectory, Z’i. Separating the two stable equilibria is a “firing threshold.” A dilation from 4, to a point b which exceeds the firing threshold will not return to Lvl, but will contract along a trajectory T2 to the shorter equilibrium length L,,n. It is important to note that the trajectories traced out in (&, L) space describe the visoelastic response of an isolated bundle of contractile filaments, i.e., one ex- periencing no resisting forces other than those generated internally by their own deformation. $$ = G(L, L,,). The exact form of the function G(L,L,J used in our cal- culations is given in Appendix 1; however, the results depend only on its qualitative features. Equations (5) and (6) describe the dynamic behavior of a single apical filament bundle. Figure 3 describes the qualitative features of the fiber behavior for the case where the neighbor forces, (Fload), vanish. 2.3. Construction of Mechanical Model The next step is to construct a mechanical model of a cell using the viscoelastic filament model. We do this as follows. Figure 5 shows a cross section of a typical cuboidal epithelial cell. Each face of the cell is repre- sented by one of the viscoelastic elements described above. For most of the applications we shall address, only the apical element need be active. The diagonal elements are required to model the in- ternal viscoelastic properties of the cell’s microtubular cytoskeleton. Finally, we shall assume that however the cell changes its shape, the internal volume remains con- stant (cf. Appendix 3). If the apical surface of the cell is stretched beyond the firing threshold, then the cell will undergo a cycle of shape change, as shown in Fig. 6. The shape history of each cell is computed by solving a collection of dif- Apical Surface Basal Surface FIG. 5. Cross section of an epithelial cell showing the arrangement of viseelastic units modeling the cytoskeleton. In our model, only the apical region is mechanically excitable. Figure 1.5: The apical filament bundle of an epithelial cell can contract (a)–(b). A viscoelas-tic element consists of a spring and a dashpot arranged in parallel and is used to model theapical filament bundle mechanics (right panel). Reprinted from Developmental Biology, 85,G.M. Odell, G. Oster, P. Alberch, and B. Burnside, The Mechanical Basis of Morp ogen-esis I. Epithelial Folding and Invagination, 446–462, Copyright 1981, with permission fromEls vi r.and deform them elves in order to undergo morphological processes such as invaginationand neurulation. The main idea was to represent the deformable actin-based cortex of anepithelial cell as a spring-dashpot system with a rest-length that would depend on a proteinor biochemical signalling. The actin cortex and the mechanical model used by Odell t al.[60] is illustrated in Figure 1.5. In (a)–(b), the apical filament bundle is illustrated atthe top of an epithelial cell. This apical filament bundle can contract. Odell et al. [60]hypothesized that when the actin-based cortex is stretched b yond critical length, theprotein or biochemical signalling would become active and change the rest-length of thecell, making it smaller. Thus, when one c ll is stretched, it will subsequently contrac . Ina tissue, when many such cells are coupled together, the contraction of one cell will imposea stretch on the neighbours and consequently induce size changes in neighbouring cells.The spring-dashpot system used by Odell et l. [60] is what I adopt a a mechanicalmodel in Chapter 3 (illustrated in the right panel of Figure 1.5). In such a system, thelength f the cell is governed by the ollowing equation, derived from Hooke’s Law and fr mthe fact that inertia is negligible in the deformation of cells:dLdt= −kλ(L− L0), (1.20)where L is the cell length, k is the spri g constant, λ is the viscous damping coefficient, andL0 is the rest-length of the spring. Note that µ is used for λ in Figure 1.5. This equationdescribes how a stretched or compressed spring will relax back towards its rest-length. Odell12et al. [60] assume that the baseline rest-length ε is decreased by the presence of a proteinsignal, C:L0 = ε+11 + σCn(1.21)where σ and n are model parameters. At the same time, the protein signal C is assumed tobe produced at a rate proportional to the cell length, decays linearly, and is auto-catalytic(positive feedback):dCdt=αC21 + βC2− νC + γL. (1.22)In this equation, the Hill function αC21+βC2describes the auto-catalytic production of theprotein signal C, and −νC and γL describes the decay and production of the signal, re-spectively.To understand the dynamics of the two-compartment ODE system above, it is possibleto study the dynamics in the phase-plane (Figure 1.6). In this phase-plane the steady-statesare identified as the intersections of the C and L nullclines (curves for which dCdt = 0 anddLdt = 0, respectively). Using parameters from [60], the system has three steady-states:(1) a high L low C stable steady-state, (2) a saddle point, and (3) a low L high C stablesteady-state. The key feature of the system is that the stable manifold of the saddle point(purple curve in Figure 1.6) separates the phase-plane into two regions. Also shown aretrajectories (grey curves) that start with different lengths, but no protein activity C. Onlywhen the initial conditions are such that the length, L(0) is above the stable manifold, willthe cell “fire” and end in the contracted state (green “firing trajectory”).In Chapter 3, I use the same spring-dashpot system to model cell mechanics and asimilar rest-length dependence. Instead of an unknown protein signal, C, I consider theeffect of GTPase signalling G on cell length. I also assume that tension, proportional to thedifference of the current length and the rest-length T ∝ L−L0, will increase the activationrate of the GTPase. As such, I suppose thatdGdt=(b+ f(T ) + γGn1 +Gn)(GT −G)−G, (1.23)where f(T ) describes how tension increase the activation rate. Depending on the strength ofthe feedback, f(T ), the coupled GTPase-length system can exist in a low-activity, relaxed-length steady-state; a high-activity, contracted steady-state; or continuously cycle betweenthese extremes (limit cycle). I illustrate this limit cycle behaviour in Figure 1.7 which showsthe phase-plane for the GTPase-tension model. The L nullcline has essentially the sameshape as in Figure 1.6. For this set of parameters, the GTPase-length system is oscillatorywith all trajectories (grey) converging to the limit cycle (shown in green). In Chapter 3, Istudy the effect of feedback strength on the behaviour of the GTPase-tension model.1310-1 100 101 10200.511.522.53Figure 1.6: Phase-plane for the mechanochemical model from Odell et al. [60]. Note thelogarithmic scale for the C-axis. Steady-states are found at the intersections of the C andL nullclines. Several trajectories are shown starting with various cell lengths L and with noprotein signal C. Only if the cell is sufficiently stretched (above the stable manifold) willthe cell end in the contracted state (high C, small L).0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.811.2Figure 1.7: Phase-plane for the GTPase-tension model. For this set of parameters, there isone unstable steady-state (at the intersection of the G and L nullclines). There is a stablelimit cycle (green trajectory) which corresponds to an oscillatory cell.14The Cellular Potts ModelThe Cellular Potts model (CPM) is a lattice-based model for modelling cell behaviour.Cells are represented as a collection of lattice sites, which can grow and shrink by addingor removing sites. I use a freely available implementation of the CPM called CompuCell3D[90]. In the CPM, the movement of cells is controlled by a total system energy, called aHamiltonian, H. To simulate cells of a certain size, for example, a volume-dependent energyterm is added to the Hamiltonian. When cells deviate from the target volume specified, thevolume energy is high. The CPM then uses a Monte-Carlo method to accept changes thatdecrease the Hamiltonian, until the energy is minimized. Such changes include invasion orretreat at the edge of some region consisting of lattice sites that we identify as a cell. Tomimic random fluctuations, even some changes that increase the Hamiltonian are acceptedwith some probability that is set by a parameter analogous to thermal energy. In this way,the cells attain their target volume.A schematic of a typical CPM is shown in Figure 1.8. Each cell is a collection of latticesites. At each step in the Monte-Carlo method (called a Monte-Carlo step (MCS)), oneor more lattice sites are selected to change. If the proposed change decreases the overallenergy of the system, i.e., ∆H < 0, then the change is accepted. If the proposed changeincreases the overall energy of the system, i.e., ∆H ≥ 0, then the change is accepted asa small noise-induced fluctuation with probability exp(−∆H/T ), where T is a numericalparameter known as the temperature. These changes are accepted to capture the noisy,stochastic nature of biophysical systems, and to avoid getting “trapped” in local energyminima of the Hamiltonian H.Other energies can be added to the Hamiltonian H. For example, in Chapter 3, I willalso include an adhesion energy that describes how cells “stick” together. The idea is thatthe adhesion energy depends on cell-cell contacts. For those cells which have large interfaceswith their neighbours have stronger adhesion and therefore lower adhesion energy. In thisway, the cells in the CPM tend to group together and remain contiguous. Additional detailsregarding the CPM simulations used in Chapter 3 can be found in Appendix B.1.4 Thesis OutlineUsing these mathematical and computational tools, I discuss two examples of multi-scalemodelling in cell biology in the next chapters. In Chapter 2, I apply quasi-steady-statemethods to models of intracellular transport by molecular motors. This chapter is self-contained (it has its own introduction and discussion) and is supplemented by AppendixA, which contains some additional details referenced in the main chapter text. In Chapter3, I use a dynamical systems approach to explore the interplay between cell signalling and151 11112222222222222222221 11111 1111333333333333333333333 3 3 3 32 2Initial cell configuration Proposed change Change acceptedFigure 1.8: A schematic of a typical CPM. Each cell, labelled 1, 2, and 3, is a collection oflattice sites. A proposed change is accepted if the change reduces the overall energy in thesystem or with some small probability if it increases the overall energy.cell mechanics, and the resulting implications on cell behaviour. Likewise, this chapter isself-contained and is supplemented by Appendix B, which contains some additional details.In Chapter 4, I conclude the thesis with a summary of the results, a discussion of thesignificance of the work, and present some questions and ideas for future work.16Chapter 2Application of Quasi-Steady-StateMethods to Nonlinear Models ofIntracellular Transport2.1 Intracellular Transport by Molecular MotorsDiffusion is a fast transport mechanism on the length scale of a typical cell, a few tensof micrometers. However, some specialized cells, including neurons, are up to 1 metrein length. This length scale imposes dramatic constraints on the transport of structural,metabolic, and signalling components from the neuronal cell body (the soma) to the ends ofdendrites or axons. Molecular diffusion is extremely inefficient for transport at such lengthscales. Fortunately, cells have evolved active transport mechanisms consisting of molecularmotors that bind to microtubule tracks and convey cargo packaged in vesicles across thecell [10].Microtubules (MTs) are asymmetric, having distinct “plus” and “minus” ends. Thetwo major types of molecular motors, kinesin and dynein, walk on microtubules in oppositedirections: kinesin walks towards the plus ends, while dynein walks towards the minus endsof MTs. Although kinesin motors can also work towards the minus ends and have otherroles within cells, I consider only those kinesin motors which walk towards the plus ends ofMTs. Both motors exist in several states, including unbound, cytoplasmic forms [5], andMT-bound as well as bound singly or in groups to cargo. The overall transport of motorsacross the cell depends on the polarity and configuration of MTs, the rates of binding toand unbinding from MTs, and the motor speeds while bound. Transport also depends onmolecular diffusion in the cytosol.One convenient experimental system is Ustilago maydis, a fungus whose long filamentous17hyphae contain MTs of mixed polarity [23, 80, 81, 86, 87]. Microtubules of mixed polar-ity also occur in the proximal regions of neuronal dendrites [2, 8, 88]. In these systems,particularly in the fungal hyphae, motors have been observed to move bidirectionally: firsttowards one cell end, and then towards the opposite end. This observation can be explainedin one of two ways. Either multiple motors (dynein and kinesin) bound to the same cargocan “take turns” pulling the load, or else a single motor, by detaching and binding to a MTof opposite polarity, would then change its direction of motion.Modelling Motor-based TransportAn intriguing question is how to approach the multi-scale problem of bridging betweenthe rates and events at the molecular level (binding, unbinding, and motor speeds) andthe overall cargo distribution and effective transport speed at the cellular level [84]. Thishas motivated the development of a number of mathematical models at various levels ofdetail. A number of efforts have dealt with the tug-of-war or teamwork of several motorsattached to a single cargo [4, 31, 40, 48, 55]. In many cases, such models mandate stochasticand computational approaches, that consider multiple states (n,m motors of distinct typesattached to a cargo, etc.). Other approaches simplify the problem to consider only a fewstates, and formulate transport equations [85] or derive such PDEs from a master-equationapproach to the stochastic motor behaviour. Examples of such approaches include (1) ananalysis and mean-field approximation of the dynamics of the totally asymmetric simpleexclusion process with Langmuir kinetics [66], (2) a study of the spontaneous formationof traffic “jams” resulting from transport on two parallel lanes (two parallel microtubuletracks) [72], (3) the incorporation of a kinetic model for motor stepping dynamics, anda study of the resulting effects on collective transport [11]. The approach here followsthe novel and elegant linear theory developed by Bressloff and Newby [7, 57] for importantinsights into motor function by deriving a quasi-steady-state (QSS) Fokker-Planck equation.The Fokker-Planck equation describes the overall transport in the system through effectivediffusion and effective velocity, which depend on the model parameters. Although thislinear theory is based on simplifications and assumptions (e.g. that the binding/unbindingkinetics are fast on the timescale of transport across the cell), it provides a useful way togain insight into the role of various parameters in determining the overall functionality ofthe transport system.In recent work, [26] used the PDE approach to model the transport of early endo-somes (cargo transported by kinesin and dynein) inside Ustilago maydis, arriving at goodagreement with experimental observations, and posing several hypotheses for further ex-perimental studies. A followup paper [16] applied the methods of [7, 57] to the examplesmotivated by [26]. In both these recent works, the models included microtubules of mixed18polarity, with and without a bias towards one end of the (1D) cell, and linear rates ofbinding and unbinding from the MT. Results in [16], for example, demonstrate that theeffective velocity of transport is the average of motor velocities, weighted by the fraction oftime spent in a given state, whereas the effective diffusivity is similarly such an average, butincludes an additional term that represents the variance in velocities of motors in differentstates.Linearity of the binding rates presumes that there is no interaction between groups ofmotors, and that binding sites are ample and unlimited. But in many biological situations,such assumptions are unwarranted. Complicated, possibly nonlinear, features have beenobserved in molecular motor traffic jams [44], and exclusion of one motor by others [79].Another case is the effect of microtubule associated proteins (MAPS) such as tau that mod-ulate the ability of motors to bind to MTs or to stay bound [19, 53]. MTs can also havevarious post-translational modifications that affect the availability or affinity of bindingsites to motors. For example, kinesin-1 binds with higher affinity to MT that have beenmodified by acetylation [70]. Considering such effects leads to models in which the bind-ing or unbinding is nonlinear and saturating, or to models that include mass-action-typereaction terms. The effect of spatially varying parameters resulting from non-homogeneousMT polarity, ATP gradients, and MAPS have been investigated in the context of intracel-lular transport in neuronal cells using quasi-steady-state methodology [57], yet the effectof nonlinear kinetic terms has been largely unexplored analytically. The need to generalizeprevious analysis to include models with such nonlinearities motivates the approach in thischapter.The Quasi-steady-state Reduction MethodThe main mathematical focus, discussed in detail in §2.3, is to extend the quasi-steady-state(QSS) reduction method introduced in [57] for reaction–advection–diffusion systems withlinear reaction kinetics to a class of problems where the kinetics are nonlinear, but wherea conservation condition is satisfied. The latter represents the fact that motors transitbetween states, but are conserved overall. The QSS method relies on the assumption thatthe nonlinear kinetics occur on a faster time-scale than the diffusion and advection processes.Owing to the conservation condition, in this limit of fast reaction kinetics, a one-parameterfamily of quasi-steady-state solutions is obtained from the equilibrium state of the kinetics.Simply put, this means that it is possible to approximate the solution of the full systemwith a single variable that will vary in space and time. When there are no eigenvalues ofthe linearization of the kinetics along this one parameter family that lie in the unstableright half-plane, this quasi-steady-state solution is referred to as a slow solution manifoldfor the full reaction–advection–diffusion system. When this condition on the Jacobian of19the nonlinear kinetics is satisfied, it is possible to use an asymptotic expansion together witha Fredholm alternative condition to derive a single scalar quasi-steady-state PDE, whicheffectively parameterizes the slow solution manifold, and approximates the solutions to thefull system.In §2.4, the asymptotic formalism of §2.3 is applied to analyze three specific nonlinearsystems for the binding and unbinding of molecular motors. These models are formulatedin §2.2 and consist of (1) a model for a single motor (“kinesin”) transiting between motionalong right-pointing MTs, diffusion in the cytosol, and motion along left-pointing MTs(with transitions only through the cytoplasmic pool), (2) a model for kinesin-dynein-cargocomplexes moving left or right along MTs or diffusing in the cytosol (interactions on a MTare assumed to lead to motor swaps that also change the direction of motion), and (3) amodel for motors (“unconventional myosin”) whose encounters with each other on an actinfilament lead to stalling. In all three cases, motors exchange between cytosolic diffusiblestates and states bound to a track (MT or actin). Nonlinearity stems from saturated bindingkinetics in (1), mass-action motor interactions leading to swaps in (2), and from stalling in(3).Overall, the QSS PDE is used to analyze the behaviour of steady-state solutions ofthe full reaction–advection–diffusion system as parameters are varied, and the results arethen interpreted biologically. The main conclusion is that in all three cases studied, theresulting QSS PDE is a conservation law for the total motor density within the cell, witheffective velocity and effective diffusion that depend nonlinearly on the model parametersand motor density. Predictions about the full model behaviour are made using the analyticalinsight gained through the QSS reduction, and the effective velocity and effective diffusionfunctions. To verify the QSS method and analysis, the steady-state and time-dependentbehaviour of both the full models and the QSS PDE are studied numerically.Summary of Results for the Molecular Motor ModelsIn the kinesin motor model, the nonlinear interactions depend on the density of cytosolicmotors. The QSS PDE describes the bulk motor distribution through effective velocity anddiffusion coefficients. These effective coefficients are related to the original velocity anddiffusion coefficients weighted by the time spent in the directed-movement and random-movement states. Moreover, the polarity distribution of MTs affects the bulk motor distri-bution by changing the sign of the velocity, which can bias the distribution of motors to theright- or left-end of the cell.Unlike the kinesin motor model, the nonlinear kinetics in the kinesin-dynein motorcomplex model arise due to a mass-action law which describes the rate at which motorcomplexes turn in response to motor complexes heading in the other direction. In this20case, the resulting QSS PDE is again a conservation law for the total density of motorcomplexes, with the advection speed dependent on the motor complex speed, the turningrate, and which motors in the complex are active. In addition to these parameters, theresulting diffusion coefficient is dependent on the binding affinity of the motor complex toMTs. A sufficiently high turning rate can reverse the distribution of motor complexes fromone end of the cell to the other, even if the probability of moving to one end of the cell ishigh.In the myosin motor model, two different QSS PDE arise from the nonlinear reactionkinetics. In the first case, the motors equilibrate between freely diffusing and walking onMT, without any motors in the stalled state. In the second case, there are some motors inthe stalled state. In the first case, the QSS PDE is linear, with effective diffusion coefficientand effective velocity mediated by the binding rate of myosin motors. We find that theasymptotic solution compares favourably with full numerical simulations of the myosinmodel. In the second case, the resulting QSS PDE is nonlinear, but is a conservation lawfor the total density of myosin motor. The effective transport rate depends on the densityof stalled motors, the velocity of stalled motors due to actin treadmilling, and the stallingrate. The effective rate of diffusion depends on all model parameters except for the velocityof stalled motors due to treadmilling. The second QSS is only valid for a range of parameterspace and stalled motor density. Outside of this range, the QSS PDE is ill-posed. A furthernovel feature of the myosin model is that the full system always converges to the first QSS,where there are no stalled motors. Through a boundary layer analysis (§2.5), I determinethat this results from the boundary conditions. I develop an alternate myosin model whichhas the same QSS approximations but, depending on initial conditions, can realize eithercase.2.2 Models of Intracellular TransportI model the cell as a 1D tube of length L0, with its left-end at x = 0. The densities ofmotors are described as number per unit cell length, with the cross-sectional area of thecell assumed to be constant. Molecular motors exist in any number, n, of possible stateswithin the cell, with pi(x, t) denoting the density of motors in state i.A reaction–advection–diffusion system describes the evolution of the vector density p ≡(p1, . . . , pn)T of motors as∂p∂t= M(p) + f(p), (2.1)where f ≡ (f1, . . . , fn)T describes the state transition rates, and M is a diagonal matrixof linear differential operators characterizing the advection and diffusion of motors in eachstate. For example, a term on the diagonal in row i, Mii = −vi ∂∂x +Di ∂2∂x2, would describe21the advection (vi) and the diffusion (Di) of motors in state fi. We assume that the endsof the cells are closed, and hence impose an overall zero-flux condition at the cell ends. Inaddition, we assume that the motors are exchanged between states in such a way that thereis no net loss or gain of motors, i.e., thatn∑i=1fi = 0. (2.2)These two assumptions result in conservation of the total density of molecular motor in thecell.The goal is to develop a theoretical framework to analyze models of the form (2.1) wherethe reaction term f is nonlinear and the reactions occur on a time-scale that is fast relativeto the time-scale of the advection and diffusion processes. This theory is then applied tothree specific nonlinear binding mechanisms. In §2.2.1 and §2.4.1 a nonlinear kinesin modelis considered, in §2.2.2 and §2.4.2 a nonlinear kinesin-dynein model is considered, whilea nonlinear myosin model is considered in §2.2.3 and §2.4.3. This analysis extends theprevious analysis for linear reaction models developed in [7, 16, 57] to allow for nonlinearreaction mechanisms.2.2.1 Kinesin ModelIn hyphae of the fungus, Ustilago maydis, for example, kinesin motors walk along micro-tubules within the cell or diffuse freely in the cytosol [16, 26, 81, 81, 86]. The densitypR(x, t) (respectively pL(x, t)) represents the population of kinesin bound to right-polarized(respectively left-polarized) MTs walking toward the end of the cell at x = L0 (respectivelyx = 0). The population of freely diffusing cytosolic kinesin is modelled by the densitypU(x, t) (U for unbound). Inside this 1D domain, 0 ≤ x ≤ L0, the density of MTs is con-stant and nonzero, with the MT distribution described by 0 ≤ P (x) ≤ 1, representing thefraction of MTs pointing to the right at a point x. Since kinesin always walks towards aMT plus end, it can reverse its direction of motion only by unbinding from a given MT andrebinding to a MT of opposite polarity. For this reason, we can assume that, in this model,motor transitions occur only through the cytosolic state. We describe the spatiotemporalevolution of the kinesin densities by the transport equations (see the schematic diagram in22x = 0 x = L0polarized MTright-moving kinesin motorleft-moving kinesin motordiffusing kinesin motorstate transitionsFigure 2.1: A schematic diagram of kinesin-based intracellular transport in a 1D cell oflength L0. Kinesin motors can bind to polarized microtubules (MTs, blue arrows), andmove to the right (purple circles with right-pointing arrows) or to the left (green circleswith left-pointing arrows). While unbound, kinesin motors are free to diffuse in the cell’scytoplasm (red circles with right and left-pointing arrows). State transitions (orange dashedarrows) occur through the freely diffusing cytosolic state.Figure 2.1):∂pR∂t= −v∂pR∂x+ Pkbg(pU)− kupR, (2.3a)∂pL∂t= v∂pL∂x+ (1− P )kbg(pU)− kupL, (2.3b)∂pU∂t= D0∂2pU∂x2− kbg(pU) + kupR + kupL. (2.3c)In Eqs. (2.3), bound kinesin moves to the right or to the left with velocity v, and D0is the diffusion coefficient for cytosolic kinesin. The unbinding rate constant is ku, whilethe binding rate constants for kinesin binding to right-polarized and left-polarized MTsare kbPg(pU) and kb(1 − P )g(pU), respectively. Here, P = P (x) is the fraction of MTspolarized towards the right in the cell. Here we have assumed a constant density of MTsacross the cell (absorbed into the constant kb). We discuss a generalization to nonuniformMT density m(x) in Appendix A.3.1. The function g(pU), possibly nonlinear, describes howother processes such as competition for binding sites or binding co-operativity are modelled.For instance, saturated binding due to a limited number of binding sites could be depictedby a term of the formg(pU) = gmpUK + pU, (2.4)for some parameters K > 0 and gm > 0. Forms such as (2.4) are obtained by assumptionstypical of Michaelis-Menten kinetics. Conservation of the kinesin motors within the cellimplies that zero-flux boundary conditions are required to model the impermeable cell23ends: (vpR − vpL −D0∂pU∂x)∣∣∣∣x=0,L0= 0. (2.5)The two additional boundary conditions are that there is no right-moving kinesin at theleft endpoint of the cell and no left-moving kinesin at the right endpoint. These boundaryconditions result from the fact that to create a flux of right-moving kinesin at a given point,there had to be a kinesin bound to a MT to the left of that point—which is impossible atx = 0, the leftmost point in the cell. A similar argument at the rightmost point in the cellestablishes the right endpoint. Thus, the following two Dirichlet conditions must hold:vpR(0) = 0 and vpL(L0) = 0. (2.6)2.2.2 Kinesin-Dynein ModelThe three-state kinesin model, formulated in §2.2.1, is a simplification of intracellular cargotransport. Cargo in fungal hyphae is typically bound to one dynein and four or five kinesinmotors at a time [81]. In this case, the entire kinesin-dynein-cargo complex may be trans-ported toward or away from the cell tip, depending on which motors are actively involvedin the transport process and the polarity of the MTs to which they are bound. In thissection, I describe a simple model for the organization and transport of cargo bound to akinesin-dynein motor complex.The populations of kinesin-dynein-cargo complexes are divided into right-moving, left-moving, and freely diffusing sub-classes, regardless of the molecular motors active in thetransport process. In the three-state kinesin model, the nonlinearities are restricted tobinding and unbinding interactions. To explore the effect of nonlinear interactions betweenmotors in distinct sub-classes, consider linear binding and unbinding interactions, but allowfor nonlinear interaction terms between the right- and left-moving species when they arein proximity on a MT. Yochelis et al. [101, 102] have recently used a model with a similarnonlinear interaction to describe the spatial organization and dynamics of unconventionalmyosin motors in actin-based cellular protrusions. I ask whether the QSS theory can beapplied to a model of this type.The population of right-moving (respectively left-) motor complexes walking toward theend of the cell at x = L0 (respectively x = 0) is described by density pR(x, t) (respectivelypL(x, t)). The population of freely diffusing cytosolic motor complexes is described bydensity pU(x, t). Here, a “binding bias” function, Q, represents the probability that whena free motor complex binds to a MT, it becomes a right-moving motor complex. Then(assuming no stalled states on the MT) the probability of becoming a left-moving motorcomplex, upon binding to a MT, is (1−Q). Since kinesin walks towards the plus ends, while24dynein walks towards the minus ends of MTs, the function Q(x) actually comprises severalbiological quantities, including local MT polarity, ratio of kinesin and dynein moleculesin a complex, as well as respective affinities to MT of these two motors. In AppendixA.3.2, I discuss how this simplification by a single function can be related to such biologicalfactors. An important distinction between this and the previous model is that now direction-switching can take place on a microtubule, and does not require unbinding into the cytosol.The above simplification allows for the detailed study of a nonlinear interaction betweenright- and left-moving populations. One possible interaction between these two populationsis that a direction change results upon an encounter with a motor-complex travelling in theopposite direction. Assume that when a right-moving complex meets a left-moving complex,the right-moving complex changes direction with rate coefficient krl. Similarly, when a left-moving complex meets a right-moving complex, the right-moving complex changes directionwith rate coefficient klr. These direction changes are due to a swap between a motor thatis actively walking, e.g., dynein, and its passive partner motor kinesin, or vice versa, in thegiven complex.Freely diffusing motor complex binds to MTs at rate kb, and diffuses in the cytosolwith diffusion coefficient D0. Bound motor complexes can move to the right (or left) withvelocity vr (or vl) or they can unbind from MTs with rate ku. These assumptions lead tothe following reaction-diffusion-advection system on 0 ≤ x ≤ L0 (see the schematic diagramin Figure 2.2):∂pR∂t= −vr∂pR∂x+ kbQpU − kupR − krlpRpL + klrpLpR, (2.7a)∂pL∂t= vl∂pL∂x+ kb(1−Q)pU − kupL + krlpRpL − klrpLpR, (2.7b)∂pU∂t= D0∂2pU∂x2− kbpU + ku(pR + pL). (2.7c)With kc ≡ krl − klr, this model can be written as∂pR∂t= −vr∂pR∂x+ kbQpU − kupR − kcpRpL, (2.8a)∂pL∂t= vl∂pL∂x+ kb(1−Q)pU − kupL + kcpRpL, (2.8b)∂pU∂t= D0∂2pU∂x2− kbpU + ku(pR + pL). (2.8c)As before, conservation of the motor complexes within the cell implies that zero-flux bound-25x = 0 x = L0MTright-moving motor complexleft-moving motor complexdiffusing motor complexstate transitionsFigure 2.2: As in Figure 2.1 but for the kinesin-dynein motor complexes. Color code asbefore for MT, and for left-moving, right-moving, or diffusing complexes. A new feature isthat state transitions can also occur through the collision of a left- and right-moving motorcomplex (orange dashed arrows, right).ary conditions are required to model the impermeable cell ends:(vrpR − vlpL −D0∂pU∂x)∣∣∣∣x=0,L0= 0. (2.9)The remaining two boundary conditions are thatvrpR(0) = 0 and vlpL(L0) = 0. (2.10)2.2.3 Myosin ModelLike kinesin and dynein motors, unconventional myosin motors are also responsible forintracellular transport in actin-based cellular protrusions, such as filopodia and stereocilia[56]. Filopodia are long, thin cellular protrusions with actin filaments at their core. Thesestructures are involved in cell motility, adhesion, and communication [51]. Stereocilia arehighly organized protrusions on hair-cells of the inner ear, responsible for hearing [82].The actin-based filamentous scaffold that supports these protrusions is known to undergoturnover. The actin-based scaffold is maintained by the delivery of new actin monomersubunits to the distal ends of the protrusions, and the disassembly of the actin bundle at itsbase [76]. The apparent motion of the actin filament bundle due to continual assembly anddisassembly at opposite ends is called treadmilling. The transport of those monomers andother material is facilitated by unconventional myosin motors [56]. In [101, 102], a reaction–advection–diffusion model was employed to describe the self-organization of waves and pulsetrains in myosin motor distribution along cell protrusions. Inspired by this model, I considera simplified system with the same nonlinear cross-species interaction term to demonstratethat the QSS method can be applied.26Consider three populations of myosin motors: bound (pB), walking (pW), and unbound(freely diffusing) (pU) in a 1D geometry. Suppose that the base of the protrusion of lengthL0 is at x = 0, but assume that the protrusion is self-contained and impose zero total-flux boundary conditions at both ends. Adapted from [101, 102], the myosin dynamics aredescribed by the following set of reaction–diffusion–advection equations on 0 ≤ x ≤ L0:∂pW∂t= −vw∂pW∂x− kˆbw(pB)2pW + kˆbpU − kupW, (2.11a)∂pB∂t= vb∂pB∂x+ kˆbw(pB)2pW − kupB, (2.11b)∂pU∂t= Df∂2pU∂x2− kˆbpU + ku(pB + pW). (2.11c)Due to actin treadmilling, bound (stalled) motors are effectively transported toward thebase of the actin bundle with the treadmiling velocity vb. Bound motors unbind with rateku and walking motors can become bound if they encounter a sufficiently high density ofbound motors (kˆbw(pB)2pW). Walking motors, on the other hand, move to the distal endof the cell protrusion with velocity vw. Walking motors may also unbind to become freelydiffusing motors. The freely diffusing motors have diffusion coefficient Df, and can reattachto an actin filament and transition to a walking motor with rate coefficient kˆb.Assume that the total flux of myosin is zero at either end of the protrusion, which givesthe boundary condition (vwpW − vbpB −Df∂pU∂x)∣∣∣∣x=0,L0= 0. (2.12)As before, there are two additional boundary conditions:vwpW(0) = 0 and vbpB(L0) = 0, (2.13)which ensures that there is no right-moving and left-moving myosin at the left and rightendpoints, respectively.2.3 Quasi-steady-state ReductionThe quasi-steady-state (QSS) reduction method, developed in [7] for the case where thevector f of state transitions is linear, will be extended to allow for nonlinear f . Here, thenonlinearities can encode nonlinear biological phenomena, such as the affect of saturatedbinding due to competition for binding sites or traffic jam-style interactions as mentioned in§2.2. In this asymptotic approach, the key assumption is that the timescale associated withtransitions between states, represented by binding and unbinding mechanisms, is short27relative to the time it takes for motors to move across the cell. This is warranted, forexample, in long cells such as fungal hypha. The separation of time scales introduces asmall dimensionless parameter ε ≈ v/(L0k), where v is the motor velocity, L0 is the celllength, and k is a typical transition rate.Using the QSS approximation, the aim is to reduce the system of transport equationsto a scalar nonlinear PDE describing the dynamics of the system for small ε. To this end,consider rescaling space and time so that the non-dimensional length of the cell is 1 and sothat one of the motor subpopulations moves with non-dimensional speed 1, scaling distanceby the cell length, and scaling time by the time it takes for a walking motor to move acrossthe cell. That is, introduce the new dimensionless variablesx? =xL0, t? =tviL0. (2.14)Under this scaling, and with the assumption that the timescale associated with transitionsbetween states is short, the system (2.1) can be written, upon dropping the stars, as∂p∂t= M(p) +1εf(p), (2.15)where f(p) represents the O(1) nonlinear motor state transition kinetics. Here M is thelinear n× n matrix differential operator in the re-scaled coordinates, with zero off-diagonalentries, so that Mij = 0 for i 6= j, and diagonal entries Mii = −vi∂/∂x + Di∂2/∂x2 fori = 1, . . . , n, with vi possibly not all unity if the right- and left-moving motors have differentspeeds. Details of the scaling leading to (2.15) for the three specific models are given inAppendix A.4.The QSS reduction method exploits the assumed small parameter ε in (2.15). On ashort time scale, where t = O(ε) so that τ = t/ε, (2.15) yields∂p∂τ= f(p) +O(ε). (2.16)Ignoring O(ε) terms, this nonlinear ODE system describes the spatially-decoupled dynamicsto leading-order on a short time-scale. Define the quasi-steady-state, p0, of (2.15) to bethe steady-state of this system, i.e., f(p0) = 0. For a general nonlinear function f , asolution to f(p0) = 0 is not guaranteed, and here consider only f such that (2.16) hasa steady-state solution. Due to the conservation (2.2) of motors within the cell, to solvef1(p0) = ... = fn(p0) = 0, it suffices to solve the under-determined algebraic systemf1(p0) = ... = fn−1(p0) = 0, and automatically find fn(p0) = 0. As such, generically, whena steady-state exists it can be written parametrically as p0 = p0(α) in terms of some scalarquantity α = α(x, t). As there may be more than one solution to f(p0) = 0, in order to28ensure that the system converges to a steady-state, I introduce the following concept of aslow manifold:Definition 2.3.1 Let p0(α) be a solution to f1 = . . . = fn−1 = 0. Then p0(α) is a slowmanifold of (2.15) provided that the Jacobian matrixJ = J(α) ≡f1p1 ... f1pn.... . ....fnp1 ... fnpn∣∣∣∣∣∣∣∣p=p0(α), (2.17)has all eigenvalues satisfying <(λ) ≤ 0 for all α on the range of definition. Moreover, λ = 0is always an eigenvalue of J for any α, i.e. Jφ = 0 for some φ 6= 0.To motivate the need for such a criterion, consider the new time-scale τ = t/ε, so that(2.15) reduces to leading-order to∂p∂τ= f(p). (2.18)In order for the ODE dynamics (2.18) to have the limiting behaviorlimτ→∞p(τ) = p0(α0), (2.19)at least for initial conditions near the slow manifold p0, where α0 is determined by theinitial condition, the eigenvalues of the Jacobian J(α) must satisfy <(λ) ≤ 0 for all valuesof α. By differentiatingf(p0(α)) = 0,with respect to α, it follows that J must always have a zero eigenvalue, i.e., thatJφ = 0, where φ =dp0dα(α). (2.20)The remaining eigevnalues of J(α) must satisfy Re(λ) < 0, which leads to the key assump-tion on the nonlinearity f .Assumption 2.3.2 Assume that the vector f of state transitions is such that there is exactlyone solution branch p0(α) to f = 0 for which the condition on the Jacobian J in Definition2.3.1 holds. Further, assume that the zero eigenvalue of J has multiplicity one for any αon its range of definition.With this assumption, it is now possible to derive a nonlinear PDE for the evolution ofα(x, t) in the quasi-steady-state p0(α). To do so, I propose an asymptotic expansion for29the solution to the full model consisting of correction terms of smaller order about the QSS.That is, I expand p as a series in ε about the quasi-steady-state asp = p0(α) + εp1 + · · · . (2.21)Substituting this expansion into (2.15) givesp0t + εp1t + · · · =1εf(p0 + εp1) + Mp0 + εMp1 + · · · . (2.22)With a Taylor expansion for the nonlinear term, together with the fact that f(p0) = 0, theO(1) terms result in:Jp1 = p0t −Mp0. (2.23)By Assumption 2.3.2, there exists a unique (up to scalar multiple) φ such that Jφ = 0.Since the eigenvalues of J and JT are identical, λ = 0 is also an eigenvalue of JT ofmultiplicity one. This guarantees the existence of a unique (up to scalar multiple) ψ suchthat ψTJ = 0T . In fact, the eigenvalue ψ = (1, . . . , 1)T is easily identified, as a result of thefact that (2.2) holds. From the Fredholm alternative, a solution to (2.23) exists if and onlyif ψT (p0t −Mp0) = 0. This solvability condition yieldsψTp0t = ψTMp0, (2.24a)which is a scalar nonlinear PDE for α(x, t). This PDE (2.24a) for α(x, t) is called theQSS PDE and the boundary conditions for α can be readily obtained from a conservationcondition (see the examples in §2.4.1, §2.4.2, and §2.4.3 below). In terms of α(x, t), theleading-order asymptoticsp ∼ p0(α(x, t)) +O(ε), (2.24b)then provides an approximate solution to the full system (2.15) when t = O(1) and awayfrom any boundary layers near the endpoints x = 0, 1. The system (2.24b) is supplementedby appropriate boundary conditions (BCs). For the three-component molecular motorssystems of §2.4, appropriate BCs are presented below. A boundary layer analysis for thesemodels is presented in §2.5.For the case where f is linear, as studied in [57] and [16], the O(ε) term in (2.24b) canbe calculated explicitly. However, in the extension of the theory to allow for a nonlinearf , it is in general analytically intractable to calculate this correction term. This differenceresults from the nonlinear interactions in the class of models studied here. To highlight thedifference, the geometry of the QSS approximation with both linear and nonlinear reactionsis illustrated in Figure 2.3. Here, p1, p2, and p3 generically reference the three different30p1p2p3f(p0) = Ap0 = 0yp0p w ∼ "w1(a) Geometry with linear reactionsp1p2p3f(p0) = 0(p01(α); p02(α); p03(α))(b) Geometry with nonlinear reactionsFigure 2.3: Geometry of the QSS approximation in the linear and nonlinear cases. In(a), with linear interactions, the solution to the full model p can be decomposed into acomponent satisfying f(p0) = Ap0 = 0, for reaction matrix A into a small correction,w ∼ εw1. In (b), with nonlinear reactions, the full model converges to the QSS f(p0) = 0quickly. On the slow manifold, the solution is parametrized by α(x, t).motor states considered in the models in §2.2. In the linear case (Figure 2.3(a)), the quasi-steady-state satisfies a linear system of equations f(p0) = Ap0 = 0. Due to conservation,A has a one-dimensional kernel. As such, the solution p can be decomposed into twocomponents: one component in the kernel, yp0, and a correction term orthogonal to thekernel, w. In the linear case, the QSS approximation describes the time evolution of y andexplains how the correction term w = O(ε). In the nonlinear case, this projection methodno longer applies, and the solution to f(p0) = 0 is more complicated. Nonetheless, the QSSapproximation in the nonlinear case suggests that the solution to the full model will quicklyconverge to the slow manifold described by f(p0) = 0, and that the slow manifold can beparametrized (due to conservation) by some scalar quantity α. The QSS approximation,in the nonlinear case, describes the spatiotemporal evolution of the parameter α(x, t) fromwhich the distribution of motors among the three states p1, p2, and p3 can be ascertained.2.4 Examples of the QSS TheoryIn this section, the QSS reduction method is appplied to the molecular motor models thatwere described in §2.2.312.4.1 QSS Reduction: Kinesin ModelAs shown in Appendix A.4.1, the kinesin model (2.3) of §2.2.1 can be scaled to a system ofthe form (2.15) wherep =pRpLpU , f(p) = kaP (x)g(pU)− pRka(1− P (x))g(pU)− pL−kag(pU) + pR + pL , M =−∂∂x 0 00 ∂∂x 00 0 D ∂2∂x2 ,(2.25)where D = D0vL0 , ε =vL0ku, and ka =kbkuif g is linear and ka =kbgmkuρor ka =kbgmkuKif g iseither a Hill or Michaelis-Menten nonlinearity, respectively. For the case of unbiased MTdistribution, with P (x) = 0.5, and linear binding function g, ka = kb/ku. In this case, karepresents the ratio of time spent in the unbound (diffusive) state to the time spent in thebound state (directed motor motion on MTs). As shown in (A.16), if g is nonlinear, thenthat ratio gets modified by other parameters reflecting the nonlinear interactions.Following the method described in §2.3, the quasi-steady-state p0(α) is found from thecondition that f = 0. Set f1 = f2 = 0 in (2.25) to getpR = P (x)kag(pU), pL = (1− P (x)) kag(pU), (2.26)which are two nonlinear equations in three unknowns. It is convenient to parameterize thefree variable by a scalar, and I set pU = α. This gives the quasi-steady-state solution branchasp0(α) = P (x)kag(α)(1− P (x))kag(α)α , (2.27)where the parameter α = α(x, t) is the unknown cytosolic motor density. A calculation ofthe Jacobian J in Definition 2.3.1 shows that J has the eigenvaluesλ = 0, λ = −1, λ = −1− kag′(α). (2.28)Therefore, a sufficient condition for p0 to be a slow manifold in the sense of Definition 2.3.1is that g is a monotonically increasing function. This condition is biologically sensible, andimplies that the rate of motors binding to MT increases with the cytosolic motor concen-tration: the more motors are in the cytosol, the more binding can take place (increasing,possibly up to some saturation level).32To derive the QSS PDE for α(x, t), I use the solvability condition (2.24) to find(1, 1, 1)∂∂t P (x)kag(α)(1− P (x))kag(α)α = (1, 1, 1)M P (x)kag(α)(1− P (x))kag(α)α .By using (2.25) for the matrix differential operator M, this expression reduces to∂∂t(kag(α) + α) = − ∂∂x(P (x)kag(α)) +∂∂x((1− P (x))kag(α)) +D∂2α∂x2,which yields the QSS PDE∂∂t(kag(α) + α) =∂∂x(D∂α∂x− (2P (x)− 1)kag(α)). (2.29)As shown in (2.80) of Appendix 2.5, to determine the boundary conditions for (2.29),substitute (2.21) into the original boundary conditions (2.5) and retain terms up to O(ε).This leads to (D∂α∂x− (2P (x)− 1)kag(α))∣∣∣∣x=0,1= 0, (2.30)which are zero-flux boundary conditions for the QSS PDE (2.29). Moreover, by integratingthe PDE (2.29) across the domain, and using the boundary conditions, the QSS PDE canbe recognized as a conservation law for the total density of kinesin motors:∂∂t∫ 10y(x, t) dx =∂∂t∫ 10(kag(α) + α) dx = 0, (2.31)where, with e ≡ (1, . . . , 1)T , I have definedy(x, t) ≡ eTp0(α(x, t)) = kag(α(x, t)) + α(x, t), (2.32)as the total density of kinesin motor in any state at (x, t). Therefore, from (2.31), the totalmotor mass satisfies∫ 10 y(x, t) dx =∫ 10 y(x, 0) dx.The QSS PDE (2.29) describes the bulk behaviour of cytosolic motors, pU = α, through-out the cell, but away from any boundary layers near the domain endpoints, when ε 1.In terms of α, using (2.27) in (2.24b) determines the behaviour of the densities of right- andleft-moving kinesin motors in the bulk region away from any boundary layers near eitherx = 0 or x = 1. The boundary-layer analysis, given in Appendix 2.5, and summarized in(2.85) for the kinesin model, shows that the right-moving and left-moving motors have aclassic boundary layer structure near x = 0 and x = 1, respectively, with a boundary layer33width of O(ε).In the case where P (x) = P is constant, the QSS PDE (2.29) reduces to∂α∂t= V(α)∂α∂x+D(α)∂2α∂x2, (2.33a)where the effective velocity V(α) and effective diffusion coefficients D(α) are defined byV(α) ≡ (1− 2P )kag′(α)kag′(α) + 1, and D(α) ≡ Dkag′(α) + 1. (2.33b)If P (x) is a smooth spatially varying function, then an additional nonlinear source/sinkterm in α, proportional to P ′(x), would appear in (2.33a).For a general g(α), the QSS PDE (2.33) provides an opportunity to make predictionsregarding the bulk behaviour of the molecular motors within the cell. The effective veloc-ity and effective diffusion coefficients V(α) and D(α) are velocity and diffusion coefficientsweighted by the fraction of time spent in directed (motor) and random (diffusive) mo-tion, respectively. These effective velocity and diffusion coefficients depend on the modelparameters.A bias in the MT polarity proportion, P , results in a corresponding bias in the effectivevelocity V(α), in such a way that V(α) is positive when P > 12 and is negative when P < 12 .Although α represents the density of cytosolic motors, it influences the behavior in the otherstates due to the assumption of rapid transitions between states. This bias agrees with theintuition that in areas where more MTs are biased to the right, more motors will be directedtowards the right end of the cell. When the MT polarity is unbiased, i.e., P = 12 , then theQSS PDE (2.33) reduces, as expected, to a nonlinear diffusion equation with no advection.In addition, when g(α) is monotone increasing, V(α) is a saturating function of kaand D(α) is a saturating function of 1/ka. Increasing kb, corresponding to increasing ka,increases the effective velocity V(α), while decreasing the effective diffusion coefficient D(α).Similarly, increasing ku, which decreases ka, causes an increase in the effective diffusion,but decreases the effective velocity. In the molecular motor system, when kb ku, so thatka 1, the expectation is that the advective processes to dominate over diffusion as motorsspend more time being transported on MTs than diffusing in the cytosol. Conversely, whenku kb, so that ka 1, we expect diffusion to dominate over advective processes, asthe motors spend less time walking on MTs than diffusing in the cytosol. The parameterdependence of V(α) and D(α) on ka in the QSS PDE (2.33) reflects this tradeoff.In the following subsections, I will explore how specific choices of the interaction functiong(α) and the MT polarity P (x) affects the QSS PDE, and further explore the parameter-dependencies discussed briefly above.34Saturated binding modelConsider the kinesin model with a saturated binding rate:g(α) ≡ α1 + cα. (2.34)This choice models the basic Michaelis-Menten biochemical kinetics with 1/c representingthe motor density at which the binding rate is 1/2 of its maximal magnitude. This choiceof g represents the idea that binding sites on MTs are limited. As cytosolic motor densityα increases, those MT sites become saturated so that g → 1. When c = 0, the binding rateis linear and the model reduces to that studied in [16].From (2.27), the quasi-steady-state p0(α) for this saturated binding kinesin model withconstant polarity P isp0(α) =P kaα(1+cα)(1− P ) kaα(1+cα)α . (2.35)Since g(α) is monotone increasing, the condition in Definition 2.3.1 holds, and p0(α) is aslow manifold. Therefore, from (2.29), the QSS PDE for α(x, t) reduces to∂∂t(kaα(1 + cα)+ α)=∂∂x(D∂α∂x− (2P − 1) kaα(1 + cα)). (2.36)Using (2.30), and as shown in (2.80) of §2.5, this QSS PDE inherits its zero-flux boundaryconditions from the full system asD∂α∂x− (2P − 1) kaα(1 + cα)= 0, at x = 0, 1. (2.37)To compare the QSS approximation with numerical approximations of the full kinesinmodel (2.15) with (2.25) and (2.34), the initial condition α(x, 0) = α0 needs to be chosensuch that the total density y is the same for the full system and the QSS PDE. Conservationof mass with the initial condition pR = 0, pL = 0 and pU = 1 at t = 0 for the full systemimplies that ∫ 10(pR(x, t) + pL(x, t) + pU(x, t))dx = 1, (2.38)for all t. Recall that the QSS PDE is a conservation law for y(x, t) = kag(α) + α, which isthe total amount of kinesin in the cell. Therefore, one correct initial condition is to chooseα0 to be the unique solution ofy(x, 0) = kag(α0) + α0 = 1. (2.39)35The steady-state solution α(x) of the QSS PDE (2.36) is the solution to the nonlocalproblem∂α∂x=kaD(2P − 1)g(α),∫ 10(kag(α) + α) dx = 1, (2.40)where g(α) is defined in (2.34). There are a few special cases for which explicit solutions to(2.40) can be found. Explicit solutions can be found if 1g(α) , P (x), and g(α) are integrable(using separation of variables to solve the differential equation). In particular, when P = 0.5,so that α = αc, where αc is a constant, (2.40) reveals thatαc =1ka + 1, (c = 0); αc =12c(c− (ka + 1) +√(c− (ka + 1))2 + 4c), (c > 0).(2.41)For the linear binding case c = 0, where ka = kb/ku, observe that the expression for αc isαc =1kb1ku+ 1kb.which represents the fraction of time spent in the unbound state (kb gives the rate at whicha freely diffusing motor binds to MTs, so 1kb gives the mean residence time in the unboundstate). In addition, for linear binding where c = 0 so that g(α) = α, then α(x) = αceβxwhere β ≡ (2P − 1)ka/D. Substituting this form into the nonlocal condition of (2.40) givesα(x) = αceβx, where αc =β(ka + 1)1(eβ − 1) , β =(2P − 1)kaD. (2.42)When the MT polarity P 6= 0.5 is also constant across the cell, this case reduces to simpleexponential distributions of all kinesin states; that distribution is biased towards the left(P < 0.5) or towards the right (P > 0.5), as previously described in [16]. However, ingeneral, the solution to the nonlocal problem (2.40) must be obtained numerically. Asshown in Appendix A.2, by recasting this nonlocal problem into an initial value problem,its solution can be computed using a simple numerical shooting procedure.In Figure 2.4(a-d), numerical approximations of the steady-state solution to the fulltransport model (dashed) and the QSS PDE (solid) for both linear binding (c = 0) andsaturated nonlinear binding (c = 1), for two constant values of the MT polarity are shown.For P = 0.5, and for c = 0 and c = 1, the advection term in (2.33) vanishes, leaving purelydiffusive motion. For P = 0.6, the MT polarity is biased to the right. Consequently, thedistributions of bound and cytosolic motors are also biased towards the right end of thecell at x = 1. From Figure 2.4(c,d), observe that the saturated binding term with c = 1slows the rate at which kinesin leaves the cytosolic compartment, causing more kinesin tobe sequestered in the middle of the cell. Further observe that the QSS approximation is not360 0.5 1x0.360.380.40.420.44αpU(a) P = 0.5, c = 00 0.5 1x00.511.5αpU(b) P = 0.6, c = 00 0.5 1x0.460.480.50.520.54αpU(c) P = 0.5, c = 10 0.5 1x00.511.5αpU(d) P = 0.6, c = 1Figure 2.4: Effect of nonlinear binding and microtubule polarity. A comparison for P = 0.5(unbiased MT polarity, left panels) and for P = 0.6 (MT biased to the right, right panels)of the steady-state cytosolic density pU(x) (dashed curves) of the full model (2.15), (2.25),and (2.34), with the steady-state α(x) (solid curves) from the QSS PDE (2.36). (a,b) linearbinding (c = 0). (Results in agreement with [16]). (c,d) Saturated nonlinear binding withc = 1. The parameters are ka = 5/3, ε = 0.02, and D = 0.1. The total mass was initiallyfixed at∫ 10 y(x, t) dx = 1, and is preserved in time. Notice the different vertical scalesbetween (a) and (b), and between (c) and (d). The QSS approximation describes the bulkbehaviour of the system well, but does not capture the boundary behaviour.valid in thin boundary layers near the two edges of the cell. These boundary layers resultfrom the reduction of the full three-equation model with four boundary conditions, to asingle PDE with two boundary conditions. The results from the boundary layer analysisgiven in (2.85) of §2.5 show that the unbound kinesin motor density pU near the twoboundaries differs from its outer approximation pU ∼ α by an error O(ε/D).In Figure 2.5, I compare the steady-state solution to the QSS approximation in thelinear binding (c = 0) and saturated binding case (c = 1), for the parameter range where370 0.5 10.90.9511.051.11.15xy(totaldensity) c = 0c = 1(a) ka = 0.10 0.2 0.4 0.6 0.8 101234xy(totaldensity) c = 0c = 1(b) ka = 5/3Figure 2.5: Effect of the relative magnitudes of binding and unbinding rates kb, ku. Steady-state solutions y(x) = eTp0(α(x)), obtained from the steady-state α(x) of the QSS PDE(2.36) with linear binding (c = 0, solid) and saturated binding (c = 1, dashed) when ka < 1(a) and ka > 1 (b). The other parameters are P = 0.6, and D = 0.1, and the total masswas∫ 10 y(x) dx = 1. In general, saturated binding results in a shallower gradient of motorsacross the cell. The steady-state behavior illustrates the effects of kb and ku. For example,for large ku (relative to kb) as in (a) where ka = 0.1, the effective velocity, V(α), is muchsmaller than the effective diffusion coefficient, D(α). This leads to a comparatively moreuniform density of motors than in (b), where kb is larger than ku, and the advection termdominates.ka < 1 (a) and ka > 1 (b) with P = 0.6. In general, saturated binding results in a shallowergradient of cytosolic motors across the cell. This result agrees with the intuition thatsaturated binding restricts the rate of binding for large motor density. This consequentlyrestricts the total number of motors walking to the right-end of the cell (P = 0.6), and inturn, saturated binding restricts the total number of motors that accumulate at the cellend.Saturated binding with a spatially variable MT polarityNext, consider a spatially varying MT polarity throughout the cell, P = P (x), in thecorresponding system of transport equations for the case of saturated binding. In this case,the quasi-steady-state p0(α) isp0(α) = P (x)kag(α)(1− P (x))kag(α)α , g(α) = α1 + cα. (2.43)38The QSS PDE, from (2.29), is∂∂t(kaα(1 + cα)+ α)=∂∂x(D∂α∂x− (2P (x)− 1) kaα(1 + cα)). (2.44)Observe that the sign of the advection term depends only on the sign of (2P (x) − 1). IfP (x) < 0.5, then advection is to the left, while if P (x) > 0.5, then advection is to the right.Biologically, if the MT polarity changes across the cell, the bulk molecular motor behaviourwill change correspondingly. If P (x) > 0.5 on some subinterval, the MT bias is to the right.This leads to a collection of motors walking to the right in this subinterval. Moreover, ifP (x) < 0.5 on some subinterval, then the bulk movement of motors in this subinterval is tothe left.To explore the effect of non-constant P (x) on the QSS PDE (2.44), consider two hypo-thetical MT polarity functions. First, considerP (x) =12[1− tanh(x− 12)], (2.45)for which P (0) ≈ 1, P (1) ≈ 0, P (12) = 12 , and P ′(x) = −12 sech2(x − 12). For x ∈ [0, 12), wehave P (x) > 12 , which indicates that the MT polarity is biased to the right in the left part ofthe cell. Similarly, for x ∈ (12 , 1], P (x) < 12 , which indicates that the MT polarity is biasedto the left in the right part of the cell. As a result of this MT polarity bias, the effectivevelocity coefficient in the QSS PDE changes signs at x = 12 . From this, it is expected thatkinesin will walk toward the centre of the cell and become “trapped” there. In Figure 2.6(a) depicts an aggregation of kinesin motors in the centre of the cell at steady-state aspredicted by the QSS PDE for both linear (c = 0) and saturated binding (c = 1). Thesteady-state problem was solved numerically by the shooting method outlined in AppendixA.2.Following [16] and [26], where molecular motor movement in the hyphae of the fungusUstilago maydis was studied, the second choice is to consider a MT polarity bias near x = 0and x = 1 that is polarized towards these cell ends, while the MTs near the cell centre pointto the right and to the left with (roughly) equal probability. As a model of such a polarityconsiderP (x) =12(1 + tanh[2(x− 12)]). (2.46)From the numerical computations of the steady-state of the QSS PDE, as shown in Figure2.6(b), observe that with such a P (x) most of the kinesin motors are pushed towards theboundaries of the cell for both linear (c = 0) and saturated binding (c = 1). This resultsfrom the highly left-biased region at the left end of the cell and the highly right-biased39region at the right end of the cell. Moreover, saturated binding sequesters more kinesinin the cytosolic compartment in the middle of the cell with a non-zero density persistingthroughout the cell at steady-state.0 0.5 100.511.52xy(totaldensity) c = 0c = 1(a) P (x) = 12(1− tanh (x− 12))0 0.5 10.511.522.5xy(totaldensity) c = 0c = 1(b) P (x) = 12(1 + tanh(2(x− 12)))Figure 2.6: Effects of two spatially dependent MT bias functions, P (x). Steady-statesy(x) = eTp0(α(x)), obtained from the steady-state α(x) of the QSS PDE (2.44) withspatially varying MT polarity where (a) MT “point towards” the cell center (describedby P (x) in (2.45)) and (b) “point towards” the cell ends (P (x) as given in (2.46)). Bothpanels depict linear (c = 0, solid) and saturated binding (c = 1, dashed). In (a), observe anaccumulation of kinesin at the center of the cell whereas in (b) the accumulation is at thecell ends. Saturated binding sequesters more kinesin motors in the cytosolic compartment,which results in the shallower, diffusion-dominated, motor distributions in the case c = 1 inboth (a) and (b). Other parameters are ka = 5/3, and D = 0.1. The total mass was fixedat∫ 10 y(x) dx = 1.Hill function bindingNext, consider a general Hill function for the binding rate, g(α), given byg(α) =αnKn + αn, (2.47)where n ≥ 1 and K > 0. Hill functions with n ≥ 2 are typically used to model positivefeedback or cooperative binding in biological systems. In this case, suppose that kinesinmotors binding cooperatively to the MTs in such a way that for low densities of motorsthe binding rate is slow, at intermediate densities (α ≈ K) binding is rapid, while forhigh densities of motors the binding rate saturates to some maximal level. The parameterK describes the value of α at which g(α) reaches half of its maximum value, while theparameter n describes the “sharpness” of the switch.With this choice (2.47) of monotonically increasing g(α), the quasi-steady-state slow400 0.5 100.511.522.5xy(totaldensity) n = 2n = 3n = 5n = 7(a) K = 10 0.5 100.511.522.53xy(totaldensity) n = 2n = 3n = 5n = 7(b) K = 0.50 0.2 0.4 0.6 0.8 101234xy(totaldensity) n = 2n = 3n = 5n = 7(c) K = 0.1Figure 2.7: Effect of the Hill function parameters K and n. Steady-states y(x) =eTp0(α(x)), obtained from the steady-state α(x) of the QSS PDE (2.44) with a Hill func-tion binding rate (2.47) for P = 0.6. The parameter K represents the density of motorspU that leads to g(pU) = 1/2 whereas the Hill coefficient n governs the “sharpness” of theHill function. Other parameters are ka = 5/3 and D = 0.1. The total mass was fixed at∫ 10 y(x) dx = 1.manifold is given in terms of g(α) by (2.43). In addition, the QSS PDE is given by (2.29) withboundary conditions (2.30). Below, I numerically examine the role of the Hill parameters nand K, and discuss the effects that these parameters have on the bulk-behaviour of kinesinwithin the cell.Figure 2.7 depicts numerical approximations to the steady-state solution of the QSSPDE for different values of n and K when P is fixed at P = 0.6. In particular, in Figure2.7(a), steady-state solutions are shown for a fixed K = 1 and for increasing n. At motordensity α = K the binding rate is half-maximal, so that g(α) < 12 for α < K. This impliesthat the advection term, ka(2P − 1)g(α), remains relatively small for α < K. This makessense, since motors hardly bind to MT at that low density. As K decreases from panel (a)to (c) of Figure 2.7, the switch to rapid binding is made possible wherever α exceeds K.For α > K, the advection term is near maximal resulting in an aggregation of kinesin motorat the right-end of the cell. Hence, decreasing K from the value 1 shifts the system fromslow-advection to fast-advection, as seen by a comparison of the bulk distribution of motorsacross the cell in panels (a), (b) and (c). The parameter n controls the “sharpness” of thetransition zone near α ≈ K in the Hill function. As n increases, the approximation g(α) ≈ 0for α < K and g(α) ≈ 1 for α > K improves. In Figure 2.7(b), K = 0.5. As n increases, theswitch from slow-advection to fast-advection becomes shaper. Hence, for large n, in regionswhere the cytosolic motor density α is larger than K, advection dominates over diffusion.Increasing n results in a sharper distribution of motors across the cell in the steady statesolution.412.4.2 QSS Reduction: Kinesin-Dynein ModelAs shown in Appendix A.4.2, the kinesin-dynein model (2.8) of §2.2.2 can be scaled to asystem of the form (2.15), wherep =pRpLpU , f(p) = kaQpU − pR − kpRpLka(1−Q)pU − pL + kpRpLpR + pL − kapU , M =−∂∂x 0 00 v ∂∂x 00 0 D ∂2∂x2 .(2.48)Here the positive dimensionless parameters v, ka, k, and D, are defined in terms of theoriginal parameters of (2.8) byv ≡ vlvr, D ≡ D0vrL0, ε ≡ vrkuL0, ka ≡ kbku, k ≡ kcρku=(krl − klr)ρku. (2.49)Without loss of generality, assume that krl > klr, so that k > 0, since the cell ends areinterchangeable. It is convenient to parameterize the quasi-steady-state solution in terms ofpL = α. In Appendix A.4.2, I determine that there is a unique quasi-steady-state solutionsatisfying f = 0 given byp0(α) =pRpLpU =Qαkα+1−Qα1ka(α+ Qαkα+1−Q) . (2.50)To determine whether this quasi-steady-state solution is a slow manifold in the sense ofDefinition 2.3.1, it is necessary to calculate the eigenvalues λ of the Jacobian of f at p = p0.The first eigenvalue is λ = 0 and the other two eigenvalues λ± satisfy the quadratic equationλ2−σ1λ+σ2 = 0; σ1 ≡ −2−ka+k(pR − pL) , σ2 ≡ 1+ka+k(1+ka)(pL−pR). (2.51)By using (2.50) for pL and pR, σ1 and σ2 are given explicitly:σ1 = −2− ka − kαH(Q), σ2 = 1 + ka + kα(1 + ka)H(Q); (2.52)where H(Q) ≡ 1 − Q1+kα−Q . A necessary and sufficient condition for Re(λ±) < 0 is thatσ1 < 0 and σ2 > 0 in (2.52). In Appendix A.4.2, I show that these inequalities hold for anyQ on 0 ≤ Q ≤ 1. Therefore, p0 is a slow manifold in the sense of Definition 2.3.1.Next, to determine the QSS PDE for α(x, t) governing the dynamics on the slow man-ifold, it is necessary to calculate the terms in the solvability condition (2.24a). This leads420 0.2 0.4 0.6 0.8 100.20.40.60.81xpR(a) pR versus x0 0.2 0.4 0.6 0.8 100.20.40.60.81xpL(b) pL versus x0 0.2 0.4 0.6 0.8 100.511.52xy(totaldensity)(c) y versus xFigure 2.8: Comparison of the full solution with the QSS solution. Shown are the steady-state of the full model (2.15) and (2.48) (dashed curve) for ε = 0.02 and the solution ofthe QSS PDE (2.53) (solid curve) for pR (a), pL (b), and the total density y at position x.The parameters are D = 0.1, ka = 2, k = 2, Q = 0.9, and v = 0.5. The total mass in thecell was fixed at∫ 10 y(x) dx = 1. The QSS approximation agrees well with the full solutionexcept near the boundary layers at the ends of the cell.to the QSS PDE for α(x, t), given by∂∂t((1 +1ka)(kα+ 1kα+ 1−Q)α)=∂∂x(V(α)α+D(α)∂α∂x), (2.53a)where the “effective transport rate” and the “effective rate of diffusion” are given byV(α) =(v − Qkα+ 1−Q), D(α) = Dka(1 +(1−Q)Q(kα+ 1−Q)2), (2.53b)together with the zero-flux boundary conditions (see (2.80) of §2.5)Vα+D∂α∂x= 0, at x = 0, 1. (2.53c)In Figure 2.8, numerical results for the motor densities pR, pL, and the total density y,in the steady-state solution of the full transport model [(2.15) and (2.48) with ε = 0.02]and in the corresponding steady-state of the QSS PDE (2.53) are compared. As shown,the full solution and the QSS solution agree well in the middle of the cell, but, as before,the QSS does not capture the boundary layer behavior near the cell ends. §2.5 provides aqualitative phase-plane analysis of the boundary-layer solutions and, in particular, predictsthat pR ≈ 0.82 at x = 1, which agrees well with the result in Figure 2.8(a).How do solutions of the QSS PDE (2.53) behave? What is the role of parameters in theoriginal model in the overall transport process? First observe from (2.50) that the density430 0.5 1x0510y(totaldensity) Q = 0.9Q = 0.7Q = 0.5Q = 0.3Q = 0.1(a) Q0 0.5 1x0510y(totaldensity)v = 2v = 1v = 0.75v = 0.5v = 0.1(b) v0 0.5 1x012345y(totaldensity) ka = 16ka = 8ka = 4ka = 2ka = 1(c) ka0 0.5 1x02468y(totaldensity) k = 8k = 4k = 2k = 1k = 0.5(d) kFigure 2.9: The effect of parameters Q, v, ka, k on the total density. The total density isy(x) = eTp0(α(x)), obtained from the steady-state α(x) of the QSS PDE (2.53). Baselineparameters are k = 2, ka = 2, D = 0.1, v = 0.5, Q = 0.9. This value of Q biases the bulkmotor distribution to the right (top labelled curve in (a)). In (a), decreasing the bindingbias Q, (probability of binding to the right) results in a shift in right-biased movement toleft-biased movement. In (b), an increase in v (the ratio of the velocites of left-movingto right-moving complexes) biases net movement towards the left. In (c), an increase inka (which represents the ratio of binding to unbinding rates kb/ku) sharpens the interfacebetween the regions of high- and low-density of motors. In (d), increasing the turning rateconstant, k, also biases the net movement to the left end of the cell. The total mass wasset to∫ 10 y(x) dx = 1.of freely diffusing motors is a weighted average of the left-moving and right-moving motorswith weight 1/ka (ratio of mean time spent bound to mean time spent freely diffusing).The density of right-moving motors at QSS, given by Qαkα+1−Q , saturates up to Q/k, as thedensity of left-moving motors, α, increases. From (2.53a) the sign of the effective transportvelocity V in (2.53b) determines the direction of motion, with the motion being to the left44if this quantity is positive. The net movement is to the left when the density of left-movingmotors, α, exceeds a threshold, i.e., when α > v(Q−1)+Qvk . For example, with fixed v andk, changing Q (which is the probability that a freely diffusing motor complex binds intothe right-moving state) will change this condition. Lowering Q increases the probabilitythat a freely diffusing motor binds into the left-moving state, which should bias the netadvection to the left. The “effective diffusivity” D of the system in (2.53b) is influenced bythe parameters D, ka, k and Q. Increasing ka decreases the effective diffusion coefficient in(2.53a), which should lead to steeper solution profiles across the cell (as usual, increasing Dhas the opposite effect). Increasing the turning parameter k also decreases the diffusivity ofthe motors. The binding bias parameter Q appears in the diffusion coefficient in two ways.First, as Q → 0 or Q → 1, the diffusion coefficient approaches the limiting value D/ka.Second, there exists a critical Q-value that maximizes the effective rate of diffusion, givena fixed motor density α and fixed k (this critical Q-value is kα+12αk+1).In Figure 2.9, I plot steady-state solutions to the QSS PDE (2.53) for a range of valuesof several parameters. These steady-states are readily calculated numerically by using anumerical shooting method (see Appendix A.2). The top labelled curve in panel (a) isproduced with a baseline parameter set (k = 2, ka = 2, D = 0.1, v = 0.5, Q = 0.9) towhich parameter variations can be compared. The total mass of kinesin-dynein complex isfixed as∫ 10 y(x) dx = 1, where y(x) = eTp0(α(x)) and p0 is defined in (2.50). Decreasingthe probability, Q, of binding to the right-moving state (panel (a)) allows for more freelydiffusing motors to bind to the left-moving state, and a shift in right-biased movement toleft-biased movement. Increasing the velocity ratio of left-moving to right-moving motorcomplexes v (panel (b)), biases net movement towards the left end of the cell, as expected.In (c), an increase in ka, which decreases the “effective diffusivity” D, sharpens the interfacebetween the regions of high- and low-density of stalled motors. In (d), increasing the turningrate constant, k, also biases the net movement to the left end of the cell. Note that highvalues of k are required to shift the behaviour from right-biased to left-biased due to thehigh baseline Q-value (Q = 0.9).2.4.3 QSS Reduction: Myosin ModelNext, I study the QSS reduction of the myosin model given in (A.24). The analysis of thismodel will differ from that of the previous two models in that there are two possible quasi-steady-state solutions. In addition, the boundary-layer behaviour will play a nontrivial rolein the dynamics.As shown in Appendix A.4.3, the myosin model (2.11) of §2.2.3 can be scaled to a system45of the form (2.15) byp =pWpBpU , f(p) =−kbw(pB)2pW + kbpU − pWkbw(pB)2pW − pBpB + pW − kbpU , M =−∂∂x 0 00 v ∂∂x 00 0 D ∂2∂x2 ,(2.54)where the dimensionless parameters v, D, ε, kbw, and kb are defined byv ≡ vbvw, D ≡ DfvwL0, ε ≡ vwkuL0, kbw ≡ kˆbwρ2ku, kb ≡ kˆbku. (2.55)Upon setting the nonlinear kinetics in the scaled myosin model (A.24a) and (A.24c) tozero, one obtains the two equations:kbw(pB)2pW − pB = 0, −kbpU + pB + pW = 0. (2.56)The two possible solutions to the first equation in (2.56) are pB = 1/[kbwpW]and pB = 0.In the latter case, the motors equilibrate between freely diffusing and walking on MT, withno motors in the bound, stalled state. In the former case, there is some proportion of motorsthat are stalled. I analyze each of these cases in turn.Type I quasi-steady-states: pB ≡ 0In the case, with pB ≡ 0, let pU be the free parameter and set pU = β(x, t). This yields thequasi-steady-statep0(β) =pWpBpU =kbβ0β . (2.57)For p0, I readily calculate that the eigenvalues λ of the Jacobian of the kinetics f(p) atp = p0 are λ = 0, λ = −1, and λ = −1 − kb. Therefore, (2.57) is a slow manifold in thesense of Definition 2.3.1. The QSS PDE for β(x, t) is calculated by expanding the solvabilitycondition (2.24a). This yields the linear PDE(kb + 1)∂β∂t=∂∂x[D∂β∂x− kbβ], 0 < x < 1; D∂β∂x= kbβ, on x = 0, 1. (2.58)The steady-state solution βs(x) of (2.58) having a unit mass, so that∫ 10 (kb + 1)β dx = 1, issimplyβs(x) =(kb(kb + 1)D)ekb(x−1)/D1− e−kb/D , (2.59)46which determines the steady-state p0[βs(x)] from (2.57).Moreover, since the time-dependent QSS PDE (2.58) is linear, it is readily solved byseparation of variables asβ(x, t) = βs(x) + ekbx/D∞∑n=1cne−λnDt/(kb+1)Φn(x), (2.60)where cn for n ≥ 1 are coefficients defined in terms of the initial data β(x, 0). Here λ =λn > 0 and Φ = Φn(x) are the positive eigenvalues and eigenfunctions of the Sturm-Liouvilleproblem(w(x)Φ′)′+ λw(x)Φ = 0, 0 < x < 1; Φ′(0) = Φ′(1) = 0, w(x) ≡ ekbx/D. (2.61)Since the myosin model (A.24) is linear when pB ≡ 0, the boundary-layer analysis nearx = 0 and x = 1 is routine for this quasi-steady-state. At steady-state, and with pB = 0in (A.24), it follows from (2.82) that there is no boundary-layer near x = 1. By solvingthe boundary layer equations (2.83) near x = 1, the leading-order uniform steady-stateapproximation ispW = kbA(ekbx/D − e−x/ε), pU = Aekbx/D, where A ≡ kbD(kb + 1)e−kb/D1− e−kb/D .(2.62)By (2.62), pU is an exponentially increasing function. By comparison, pW has a rapidlydecaying correction factor (since 1/ε is large in the second exponential), which produces asmall “knee” in its graph, Figure 2.10 (a), close to the origin.Numerical results reveal that the steady-state (2.62) with pB = 0 is realizable from thelong-time dynamics of the full transport model (A.24) with different initial states for pW,pB, and pU at t = 0. Figure 2.10 depicts the numerical solution pW and pU to (A.24) att = 130 for the parameter values ε = 0.02, kb = 0.3, kbw = 0.5, and D = 0.1, when theinitial densities are spatially uniform and equally-partitioned as pW = pB = pU = 1/3 att = 0. The full dynamics quickly drives pB to zero as t increases. From Figure 2.10, notethat at t = 130 the computed motor densities pW and pU from the full model agree wellwith the steady-state asymptotic result (2.62).470 0.2 0.4 0.6 0.8 100.20.40.60.8xpW (Full numerics)(Asymptotics)(a) pW versus x0 0.2 0.4 0.6 0.8 100.511.522.5xpU (Full numerics)(Asymptotics)(b) pU versus xFigure 2.10: Full numerical vs. asymptotic solutions to the myosin model. Shown aresteady-state motor densities (solid curves) pW (a) and pU (b) (shown at t = 130) computedfrom the full time-dependent myosin transport model (A.24) for ε = 0.02 and with thespatially uniform initial condition pW = pB = pU = 1/3 at t = 0, so that the total mass isunity. The parameters are kb = 0.3, kbw = 0.5, and D = 0.1. Although pB > 0 at t = 0, thedynamics quickly drives pB to zero as t increases. The dashed curves in (a) and (b) are theasymptotic results (2.62) for the steady-state, which compare favorably with the numericalresults.Type II quasi-steady-states: pB > 0It is also possible to let pB 6= 0 be the free parameter, and define pB = α(x, t). Upon solving(2.56) for pW and pU, the quasi-steady-state solution for (A.24) is given byp0(α) =pWpBpU =1kbwαα1kb(α+ 1kbwα) . (2.63)The eigenvalues λ of the Jacobian of the kinetics f(p) at p = p0 reveal whether p0 is aslow manifold in the sense of Definition 2.3.1. One eigenvalue is λ = 0, while the remainingtwo eigenvalues λ± satisfy the quadratic equation λ2 − σ1λ+ σ2 = 0, where σ1 and σ2 aregiven byσ1 = −2− kb + 2kbwpBpW − kbw(pB)2, (2.64a)σ2 =(1− 2kbwpBpW)(1 + kb + kbw(pB)2) + 2k2bw(pB)3pW − kb + kb(1 + kbw(pB)2),(2.64b)with pB and pW as given by the entries in (2.63). Upon using (2.63) for p0, σ1 and σ2 areσ1 ≡ −kb − α2kbw, σ2 ≡ (kb + 1)(α2kbw − 1). (2.65)48Since σ1 < 0, a necessary and sufficient condition for Re(λ±) < 0 is that σ2 > 0 in(2.65). From the expression for σ2 in (2.65), it follows that p0 is a slow manifold wheneverkbw > 1/α2.For kbw > 1/α2, the QSS PDE results from the solvability condition (2.24a). This yieldsthat(1, 1, 1)∂∂t1kbwαααkb+ 1kbkbwα = (1, 1, 1)M1kbwαααkb+ 1kbkbwα . (2.66)By calculating the various terms in this expression, the following nonlinear QSS PDE isobtained for α(x, t):∂∂t((kb + 1)(kbwα2 + 1)kbkbwα)=∂∂x(V(α)α+D(α)∂α∂x), (2.67a)where the “effective transport rate” and the “effective rate of diffusion” are given byV(α) = vα− 1kbwα, D(α) = D (kbwα2 − 1)kbkbwα2. (2.67b)From (2.80) of §2.5, the zero-flux boundary conditions for this conservation law areVα+D∂α∂x= 0, at x = 0, 1, (2.67c)which are exactly zero-flux boundary conditions for the QSS PDE (2.67). From (2.67b)we observe that the advection direction depends on the sign of V. In particular, if α <1/(√vkbw), the net movement is to the right. By integrating the QSS PDE over the domain,and by using (2.67c), we obtain a conservation law for y(x, t) = eTp0[α(x, t)], where p0(α)is defined in (2.63). For all t > 0, we obtain in terms of α(x, t) that∫ 10y(x, t)dx =∫ 10y(x, 0) dx , y(x, t) ≡ (kb + 1)kbkbw(kbwα2 + 1)α. (2.68)We remark that on the range kbwα2 − 1 > 0 for which p0 is a slow manifold for thedynamics, the QSS PDE (2.67a) is well-posed in that the diffusion coefficient in (2.67a) ispositive. In fact by expanding (2.67a), we obtain that (2.67a) is equivalent to the followingPDE with a constant diffusivity D/(kb + 1),∂α∂t=Dkb + 1∂2α∂x2+kbwkb(kb + 1)(kbwα2 − 1)((vα2 +1kbw)∂α∂x+2Dαkbkbw(∂α∂x)2). (2.69)Alongside the transport term involving ∂α∂x , the source term2Dαkbkbw(∂α∂x)2describes how49gradients in α can lead to an increase in motor density, especially for low densities (so that1/α is large).Steady-state solutions to the QSS PDE (2.67) are solutions to the nonlocal problemdαdx= −kbD(vkbwα2 − 1)kbwα2 − 1 α ,(kb + 1)kbkbw∫ 10(kbwα2 + 1)αdx = 1 , (2.70)provided that kbwα2 − 1 > 0 on 0 ≤ x ≤ 1. Here, the total mass has been fixed as∫ 10 y(x, 0) dx = 1. It is possible to use the numerical shooting method described in AppendixA.2 to solve (2.70) and, further, to numerically identify the region in the kbw versus kbparameter space where kbwα2 − 1 > 0 on 0 < x < 1. For D = 0.1, this region is shown inFigure 2.11(a) and in Figure 2.11(b) for v = 0.1 and v = 0.5, respectively.2 4 6 8 10020406080kbkbw No Existence(a) v = 0.12 4 6 8 1001020304050kbkbw No Existence(b) v = 0.5Figure 2.11: Region of solution existence (unshaded). Shown are the regions in the kbwversus kb parameter space where a steady-state to the myosin model Type II QSS PDE(2.67) for D = 0.1 exists when (a) v = 0.1 and (b) v = 0.5. In the shaded regions, there isno steady-state to the Type II QSS PDE. On the boundary of these regions α = 1/√kbwat x = 0. The total mass was fixed at∫ 10 y(x, 0) dx = 1. The points marked in the leftand right panel are parameter values corresponding to solutions shown in Figure 2.13 andFigure 2.12, respectively.Figure 2.12(a)-(c) depicts the QSS motor-densities pB(x), pW(x), and pU(x), for threevalues of kb corresponding to taking a horizontal slice at fixed kbw = 12 through theparameter plane of Figure 2.11(b) with v = 0.5. In terms of α(x), these densities are givenby (2.63). From Figure 2.12(a)-(b), observe that as kb increases there is an accumulationof bound myosin motors, with a corresponding decrease in walking myosin motors near theleft end of the cell. From Figure 2.12(c), observe that as kb increases, there is a decrease inunbound freely diffusing motors in the cytosolic compartment in the middle of the cell.Figure 2.13(a-c) depicts the QSS motor-densities pB(x), pW(x), and pU(x), for threevalues of kbw corresponding to taking a vertical slice at fixed kb = 3.0 through the phase-500 0.2 0.4 0.6 0.8 1012345xpB=α kb = 1.6kb = 2.6kb = 3.6(a) pB versus x0 0.2 0.4 0.6 0.8 100.050.10.150.20.25xpW kb = 1.6kb = 2.6kb = 3.6(b) pW versus x0 0.2 0.4 0.6 0.8 100.511.5xpU kb = 1.6kb = 2.6kb = 3.6(c) pU versus xFigure 2.12: Effect of the (scaled) binding rate, kb. The QSS densities pB (a), pW (b), andpU (c), computed from (2.70) and (2.63), are plotted for three values of kb corresponding totaking a horizontal slice through the parameter space of Figure 2.11(b) with fixed kbw = 12.Other parameters are D = 0.1, and v = 0.5. The total mass was fixed at∫ 10 y(x) dx = 1.0 0.2 0.4 0.6 0.8 10.20.40.60.81xpB=α kbw = 19.43kbw = 22kbw = 25(a) pB versus x0 0.2 0.4 0.6 0.8 100.050.10.150.20.25xpW kbw = 19.43kbw = 22kbw = 25(b) pW versus x0 0.2 0.4 0.6 0.8 10.20.250.30.35xpU kbw = 19.43kbw = 22kbw = 25(c) pU versus xFigure 2.13: Effect of the (scaled) stalling rate, kbw. As in Figure 2.12 but for three valuesof kbw corresponding to taking a vertical slice through the parameter space of Figure 2.11(a)with fixed kb = 3.0.diagram of Figure 2.11(a) with v = 0.1. Observe from Figure 2.13(a-b) that as the transitionrate kbw between walking to bound motors increases, there is a decrease in walking motors,with a corresponding increase in bound motors near the left end of the cell.Finally, Figure 2.14(a,b) depicts the QSS motor densities for v = 0.1 and v = 0.5,respectively, for the parameters kb = 3, kbw = 20, and D = 0.1. As the treadmillingspeed, v, increases from v = 0.1 to v = 0.5, note that the system switches from right-biasedadvection to left-biased advection. This matches the observation that net movement is tothe right if pB ≡ α < 1/√vkbw. For small treadmilling velocity v, this condition is moreeasily satisfied since the quantity 1/√vkbw is large.Two notable features distinguish the myosin model from previous models discussedherein. The first is existence of two possible QSS approximations, as shown. A secondfeature pertains to the boundary-layer behavior near x = 0 and x = 1. This is analyzed in510 0.2 0.4 0.6 0.8 100.20.40.60.81xmotordensities PBPWPU(a) v = 0.10 0.2 0.4 0.6 0.8 10123456xmotordensities PBPWPU(b) v = 0.5Figure 2.14: Effect of treadmilling speed, v. QSS densities pB, pW, and pU, computed from(2.70) and (2.63), for (a) v = 0.1 and (b) v = 0.5. As v increases, the system switches fromright-biased advection to left-biased advection. Other parameters are kbw = 20, kb = 3 andD = 0.1. The total mass is∫ 10 y(x, 0) dx = 1.detail in §2.5.3 based on the full myosin transport model (A.24) near x = 0. There, usingphase-plane analysis, I explain that it is always possible to insert a boundary layer nearx = 0 to satisfy pW = 0 at x = 0. However, §2.5.3 also shows that there is no steady-stateboundary-layer solution near x = 1 that allows the extra boundary condition pB = 0 atx = 1 to be satisfied. This difficulty results from the fact that pB = 0 is the slow manifoldfor the Type I solutions of §2.4.3. Since no steady-state boundary layer solution exists inthe full model, any non-zero density of stalled motors pB will tend to 0 via a backwardspropagating wave that leaves pB = 0 in its wake. The full myosin model converges to a TypeI QSS (2.57) regardless of the initial condition. An example of this behaviour is shown inFigure 2.15, where kbw = 25, kb = 3, D = 0.1, v = 0.5, and ε = 0.02. As a result, drawingconclusions about the behavior of the full system from the QSS PDE becomes difficult.This leads to the question of which QSS PDE, Type I or Type II, better describes the bulksystem dynamics.One possible regularization to overcome this problem with the boundary-layer near x = 1is to add an asymptotically small diffusion term ε1pBxx to (A.24), where ε1 = O(ε). Thisregularization term does not affect the quasi-steady-states at leading order. The additionof such a small “regularizing” diffusion term also appears in the traveling-wave analysis of[102]. The fully scaled model is as in (A.24c), but with the additional small diffusion term520 0.5 1x00.10.20.30.4t012Figure 2.15: pB(x, t) converges to Type I QSS. The density of bound motors, pB(x, t), tendsto zero behind a wave propagating backwards from x = 1. The full myosin model convergesto a Type I QSS as no non-trivial steady-state solution satisfies the boundary conditionpB(1) = 0. Parameters are kbw = 25, kb = 3, D = 0.1, v = 0.5, and ε = 0.02.in the pB equation:∂pW∂t= −∂pW∂x+1ε(−kbw(pB)2pW + kbpU − pW), (2.71a)∂pB∂t= ε1∂2pB∂x2+ v∂pB∂x+1ε(kbw(pB)2pW − pB), (2.71b)∂pU∂t= D∂2pU∂x2+1ε(pB + pW − kbpU). (2.71c)The boundary conditions are as before, (2.12) and (2.13), but instead of pB(1, t) = 0, it isnecessary to impose that∂pB∂x(0, t) = 0 and∂pB∂x(1, t) = 0, (2.72)for conservation of mass.In this case, both Type I and Type II QSS PDE are valid approximations of the fullsystem, and it is possible to add steady-state boundary layers near x = 0 and x = 1 forthe regularized model (2.71). However, it is intractable anaytically to analyze the globalbehavior of time-dependent solutions for (2.71), so as to predict which of the two types ofQSS PDEs will result from an an arbitrary initial state. Figure 2.16(a) depicts that thefull model (2.71) with asymptotically small diffusion term ε1pBxx has a steady-state withnon-zero pB, and that solutions can converge to the Type II QSS (as compared with Figure2.15). In this case, as shown in Figure 2.16(b), solutions to the full myosin model and theType II QSS PDE agree as expected.Due to the existence of two QSS solutions, it is expected that the initial condition for53(2.71) determines whether the full myosin model converges to the Type I or Type II QSS. Toelucidate this hypothesis, I fix the model parameters kbw = 25, kb = 3, D = 0.1, v = 0.5,ε = 0.02, and ε1 = 0.005 and numerically determine to which QSS the regularized fullsystem of PDE’s (2.71) converges for a range of spatially homogenous initial conditions:pW(x, 0) = c1, pB(x, 0) = c2, with 0 ≤ c2, c2 ≤ 1, c1 + c2 ≤ 1 and pU(x, 0) = 1 − c1 − c2(to ensure conservation of total mass). In Figure 2.17, the results of this exploration areshown in a phase-diagram. For a given pair of spatially homogenous initial conditions,(pB(x, 0), pW(x, 0)), a circle indicates that the model (2.71) converges to a Type I QSS, whilea cross indicates that the model (2.71) converges to a Type II QSS. The line on the phase-diagram indicates the unstable manifold which emanates from a saddle-point steady-statein the myosin-model reaction kinetics (the non-spatial myosin-model). For a phase-planeanalysis of the non-spatial model, see §A.4.4. Below this unstable manifold, the solutionsconverge to a steady-state with pB = 0, similar to a Type I QSS. Above this unstablemanifold, solutions converge to a steady-state with pB > 0, similar to a Type II QSS. Thediscrepancy between the unstable manifold computed from the non-spatial model and thephase-diagram from the fully-spatial model indicates that the spatial processes enlarge theregion of attraction for Type II QSS with non-zero pB.0 0.5 1x00.511.52pWpBpU(a) Steady-state myosin densities0 0.5 1x02468y(totaldensity)QSSFull system(b) QSS and myosin model comparisonFigure 2.16: Steady-state behaviour of the regularized myosin model. (a) The steady-statebehaviour of the full myosin model with small pB diffusion term. Note that pB > 0 for all x.(b) A comparison of the total density of myosin in the Type II QSS approximation and in thefull model. Note that this Type II QSS behaviour approximates the full system dynamicswell. Parameters are kbw = 25, kb = 3, D = 0.1, v = 0.5, ε = 0.02, and ε1 = 0.005. Thetotal mass was fixed at∫ 10 y(x) dx = 1.540 0.5 100.51(a) Phase-diagram0 0.1 0.200.050.10.150.2(b) Enlarged version of (a)Figure 2.17: Myosin model initial condition dependence. The steady-state behaviour ofthe full myosin model with pB diffusion depends on the initial conditions. For a given pairof spatially homogenous initial conditions,(pW(x, 0), pB(x, 0))= (c1, c2), with pU(x, 0) =1 − c1 − c2, the solution will converge to a Type I steady-state, with∫ 10 pB(x) dx = 0(indicated by a circle), or to a Type II steady-state (cross), with∫ 10 pB(x) dx > 0. The line inthe phase-diagram represents the unstable manifold computed from the non-spatial myosinmodel (§A.4.4). In the non-spatial model, the solution converges to a Type I steady-state(pB = 0) with initial conditions below this unstable manifold, while the solution convergesto a Type II steady-state (pB > 0) with initial conditions above this unstable manifold.Parameters are kbw = 25, kb = 3, D = 0.1, v = 0.5, ε = 0.02, and ε1 = 0.005. The totalmass was fixed at∫ 10 y(x) dx = 1.2.5 Boundary Layer AnalysisIn this section, the appropriate boundary conditions for the QSS PDEs are determined, andthe boundary layers in the full models near x = 0, 1 are analyzed. In particular, I explainhow the QSS PDE inherits the boundary condition from the full reaction-advection-diffusionPDE system. Using the method of matched asymptotics, I also analyze the boundary layerbehaviour of solutions to the full PDE system to explain estimate the error in the QSSapproximation.The discussion begins with general three-component systems on 0 ≤ x ≤ 1 of the formp1t = −v1p1x + f1ε, p2t = v2p2x +f2ε, p3t = Dp3xx +f3ε, (2.73a)where v1, v2, D are positive O(1) constants, ε 1, and the kinetics fj = fj(p1, p2, p3) for55j = 1, 2, 3, satisfy the conservation conditionf1 + f2 + f3 = 0. (2.73b)The three models, kinesin, kinesin-dynein, and myosin are systems of PDE of the form(2.73a). Imposing the mass constraint ∂t∫ 10 (p1 + p2 + p3) dx = 0, and setting p1(0, t) =p2(1, t) = 0, reveals the following boundary conditions for (2.73a):Dp3x + v2p2 − v1p1 = 0, at x = 0, 1; p1(0, t) = 0, p2(1, t) = 0. (2.73c)Matched asymptotic analysisAs in the QSS reduction in the previous sections, assume that there is a unique one-parameter family p0(α) ≡ (p01(α), p02(α), p03(α))T of solutions to the leading-order problemf = (f1, f2, f3)T = 0, and that p0 is a slow manifold for (2.73) in the sense of Definition2.3.1. Then, as was shown in §2.3, α = α(x, t) satisfies the QSS PDE (2.24a), which can bewritten as∂t(p01 + p02 + p03)= ∂x(−v1p01 + v2p02 +D∂xp03) . (2.74)The QSS solution is known as the outer solution, which is valid away from the boundariesx = 0, 1.To determine an appropriate boundary condition for (2.74) as x → 0+, I analyze theboundary layer structure for (2.73) near the left endpoint x = 0. As x → 0+, then it isexpected that the solution to the full system will agree with the outer solution (from theQSS PDE) with errors that depend on x:p1 = p010 +O(x), p2 = p020 +O(x), p3 → p030 + xdp03dx∣∣x=0+ · · · , (2.75)where p0j0 ≡ p0j (α(0, t)) denotes the QSS solution evaluated at x = 0. for j = 1, 2, 3. Sincethe cell-ends are interchangeable, only an analysis of the boundary layer near x = 0 ispresented—a similar analysis can be done near x = 1.To determine the width of the boundary layer, consider the dominant balances in thefull system. For t = O(1) the two possible balances for the spatial derivatives in (2.73a)near x = 0 are x = O(√ε) and x = O(ε) (considering the presence of a first-order and asecond-order spatial dertivatives and the 1ε multiplying the nonlinear function f). On the56wider scale, use the change of variables ξ = x/√ε to obtain from (2.73a) thatp1t = − v1√εp1ξ +f1ε, p2t =v2√εp2ξ +f2ε, p3t =Dεp3ξξ +f3ε. (2.76)In this case, for ε 1, the leading-order contributions are from those terms with factor 1ε .At leading order, f1 = f2 = 0 (from the first two equations in (2.76)). This implies thatf3 = 0 thanks to conservation (2.73b). Since the QSS solution satisfies f1 = f2 = f3 = 0,it follows that on this wide-scale that p1 ∼ p010, p2 ∼ p020, and p3 ∼ p030 for x = O(√ε). Inother words, the QSS approximation is still valid on this wide-scale when x = O(√ε).The other dominant balance for spatial derivatives in (2.73a) is x = O(ε). To study thisregion, introduce η ≡ x/ε, and obtain from (2.73a) thatεp1t = −v1p1η + f1, εp2t = v2p2η + f2, εp3t = Dεp3ηη + f3. (2.77a)From (2.73c), the boundary conditions for this system areDεp3η + v2p2 − v1p1 = 0, at η = 0; p1(0, t) = 0, (2.77b)while the asymptotic matching conditions, as obtained from (2.75), are thatp1 ∼ p010, p2 ∼ p020, p3 ∼ p030 + εηdp03dx∣∣x=0, as η →∞. (2.77c)For t = O(1), neglect the asymptotically negligible left-hand sides of (2.77a) to obtain− v1p1η = −f1, v2p2η = −f2, Dεp3ηη = −f3. (2.78)By adding the equations in (2.78), using the conservation condition (2.73b), and afterintegrating in η, for all η > 0,Dεp3η − v1p1 + v2p2 = A, (2.79)where A is independent of η. Evaluating this expression at η = 0, (2.77b) yields that A = 0.With A = 0, evaluating (2.79) as η →∞ by using the matching condition (2.77c) yieldsDdp03dx− v1p01 + v2p02 = 0, at x = 0. (2.80a)This key result shows that to obtain the boundary condition at x = 0 for the QSS PDEfor α(x, t), it is possible to substitute the outer approximation p1 = p01(α), p2 = p02(α), and57p3 = p03(α), into the first condition of (2.73c). In this sense, the QSS PDE inherits theno-flux boundary condition (2.73c) at x = 0. Note that a similar analysis can be done nearx = 1, with the analogous result thatDdp03dx− v1p01 + v2p02 = 0, at x = 1. (2.80b)Next, to fully characterize the beahviour of the boundary layer and to estimate the errormade with the QSS approximation, I complete the boundary layer analysis near x = 0 byasymptotically expanding:p3 = p030 +εDP3 + · · · , (2.81)and obtain from the first two equations in (2.78), together with (2.79) with A = 0, thefollowing boundary-layer problem on 0 < η <∞:v1p1η = f1(p1, p2, p030); p1(0) = 0, p1 → p010 as η →∞, (2.82a)v2p1η = −f2(p1, p2, p030); p2 → p020 as η →∞, (2.82b)P3η = v1p1 − v2p2; P3η ∼ Ddp03dx∣∣x=0as η →∞. (2.82c)Here, the first two equations result from the dynamics of the full model in the region wherex = O(ε), i.e., with spatial coordinate η. Boundary conditions for these equations resultfrom the fact that p1(0) = 0, and that both p1 and p2 should match with the outer solution(from the QSS) as η →∞. Finally, the boundary conditions at x = 0 (2.79) imply that thecorrection term P3 must satisfy the P3η = v1p1 − v2p2, with the derivative matching thederivative of the outer solution as η →∞.Although the first two equations for p1 and p2 are uncoupled from P3, in general it isnot possible to calculate p1 and p2 analytically when f1 and f2 are nonlinear in p1 and p2.However, the system for p1 and p2 can be readily studied qualitatively in the phase-plane.A similar boundary layer analysis can be done near x = 1. To study this boundarylayer, define η = (1− x)/ε to find, in place of (2.82a) and (2.82b), thatv1p1η = −f1(p1, p2, p031); p1 → p011 as η →∞, (2.83a)v2p1η = f2(p1, p2, p031); p2(0) = 0, p2 → p021 as η →∞. (2.83b)Here p0j1 ≡ p0j (α(1, t)), for j = 1, 2, 3.Next, I will study the phase-plane behaviour of the boundary layer solution for thekinesin, kinesin-dynein, and the myosin models.582.5.1 The Kinesin ModelFor the kinesin model (2.25) of §2.4.1, the boundary layer system (2.82) can be solvedexplicitly. With the QSS approximation p0, as given in (2.27), note that v1 = v2 = 1,p1 = pR, p2 = pL, and p3 = pU. From the QSS for the kinesin model, (2.27), the outersolution gives p010 = kaP (0)g(α0), p020 = ka [1− P (0)] g(α0), and p030 = α0, where α0 ≡α(0, t). Therefore, using the reaction kinetics in (2.25), the boundary layer problem (2.82)becomesp1η = p010 − p1, p2η = −p020 + p2, P3η = p1 − p2. (2.84)The solution with p1(0) = 0, p1 → p010 and p2 → p020 as η →∞, is simply p1 = p010(1− e−η),and p2 = p020. Then, P3 is obtained up to a constant by integrating the last equation in(2.84). In this way, the boundary layer solution for x = O(ε) is thatpR ∼ p010(1− e−x/ε), pL ∼ p020, pU = p030 +εD(ηDdαdx∣∣x=0+ p010e−η +B),(2.85)where the constant B can only be determined from a two-term outer QSS solution, thatis intractable analytically. This analysis shows two key features. First, the right-movingmotors have a classic boundary-layer behaviour when x = O(ε). Second, for x = O(ε) theunbound kinesin motor density pU differs from its outer approximation only by an errorO(ε/D). A similar calculation can be done for the boundary layer near x = 1 using (2.83).2.5.2 The Kinesin-Dynein ModelFor the kinesin-dynein model (2.48), I study the boundary-layers equations (2.82) for thelayer near x = 0 qualitatively in the phase-plane. Using f in (2.48), and setting v1 = 1 andv2 = v, (2.82a) and (2.82b) on 0 < η <∞ becomep1η = −p1 − kp1p2 + kaQp030, p1(0) = 0, p1 → p010 ≡Qα0kα0 + 1−Q as η →∞,(2.86a)p2η = −1v[ka(1−Q)p030 − p2 + kp1p2], p2 → α0 as η →∞, (2.86b)where p030 = (kα0 + 1)α0/[ka(kα0 + 1−Q)]. To analyze (2.86) in the phase-plane, it isconvenient to introduce new variables q1(η) and q2(η) defined byp1 =r2kq1, p2 =r1kq2, where r1 = kα0, r2 ≡ Qr1r1 + 1−Q. (2.87)59In terms of q1 and q2, (2.86) transforms to the two-component dynamical systemq1η = g1(q1, q2) ≡ (1− q1) + r1(1− q1q2), q1(0) = 0, q1 → 1 as η →∞,(2.88a)q2η = g2(q1, q2) ≡ −1v[1− q2 + r2(q1q2 − 1)] , q2 → 1 as η →∞. (2.88b)This system for q1 and q2 is more easily studied qualitatively, as it has a equilibrium solutionat q1 = q2 = 1, and the asymptotic matching conditions require that q1 and q2 tend 1 asη → ∞, to match with the outer solution. Moreover, the initial condition requires thatq1(0) = 0. As such, in the phase-plane, I seek to show the existence of a trajectory fromthe q2 axis for η = 0 that converges to q1 = q2 = 1. Below, I will argue that the equilbrium(1, 1) is a saddle point for the dynamics and demonstrate that such a trajectory does exist.This trajectory will correspond to the boundary layer solution after changing coordinatesback to p1 and p2.Note that r2 depends on r1. With r2 considered as a function of r1: r2 = 0 whenr1 = 0; r2 → Q < 1 as r1 → ∞; and r2 is monotone increasing in r1 since dr2/dr1 =[Q(1−Q)]/(r1 + 1−Q)2 > 0 holds for 0 < Q < 1. It follows that 0 < r2 < 1 for anyr1 > 0.The determininant of the Jacobian Jg of g1 and g2 at the equilibrium state q1 = q2 = 1isdet(Jg) = − 1v (kα0 + 1−Q)[(1−Q)(1 + 2kα0) + kα20]< 0,revealing that q1 = q2 = 1 is a saddle point for the dynamics. In Figure 2.18(a), I plotthe phase portrait q2 versus q1 and nullclines for (2.88) for representative values r1 = 2,r2 = 0.5, and v = 0.5. The q2 nullcline intersects the q2 axis at q2 = 1 − r2 ∈ (0, 1) since0 < r2 < 1. This plot indicates the existence of a unique value q2(0) = q02 > 1 − r2 forwhich (2.88) has a solution with (q1, q2) → (1, 1) as η → ∞. This solution corresponds tothe stable manifold of the saddle point, which intesects the q2-axis so that q1(0) = 0. Thesolution is the boundary-layer solution after returning to the original variables p1 and p2.This qualitative analysis confirms the existence of a boundary-layer solution near x = 0 forthe kinesin-dynein model for a range of parameters.A similar phase-plane analysis can be done to analyze the boundary-layer system (2.83)600 0.5 1 1.5 200.511.52(a) r1 = 2.0, r2 = 0.5, v = 0.50 0.5 1 1.5 200.511.52(b) r1 = 1.69, r2 = 0.85, v = 0.5Figure 2.18: Qualitative analysis of boundary layer behaviour of the kinesin-dynein model.Phase portraits of q2 versus q1 for boundary layer solutions of the kinesin-dynein modelnear x = 0 (a) and near x = 1 (b) from (2.88) and (2.89), respectively. In (a) there is aunique value q2 = q02 at q1 = 0 for which (2.88) has a solution with (q1, q2) → (1, 1) asη → +∞. In (b) there is a unique value q1 = q01 at q2 = 0 for which (2.89) has a solutionwith (q1, q2) → (1, 1) as η → ∞. The parameter values r1, r2, and v for (b) are thoseconsistent with Figure 2.8.near x = 1. In place of (2.88), the boundary-layer system is nowq1η = −g1(q1, q2) ≡ − [(1− q1) + r1(1− q1q2)] , q1 → 1 as η →∞, (2.89a)q2η = −g2(q1, q2) ≡ 1v[1− q2 + r2(q1q2 − 1)] , q2(0) = 0 q2 → 1 as η →∞,(2.89b)where, in place of (2.87), r1 and r2 are now defined by r1 = kα1 and r2 ≡ Qr1/(r1 + 1−Q),and α1 = α at x = 1. Also note that the asymptotic matching conditions for η → ∞ aredifferent. Here, the solution must satisfy q2(0) = 0. As such, I seek a solution that intersectsthe q1 axis and converges to (1, 1). Figure 2.18(b) depicts the phase portrait and nullclinesfor (2.89) for r1 = 1.69, r2 = 0.85, and v = 0.5. This corresponds to the parameter valuesused in Figure 2.8. The phase portrait shows the existence of a unique value q1(0) = q01for which (2.89) has a solution with (q1, q2) → (1, 1) as η → ∞. As before, the solutioncorresponds to the stable manifold of the saddle point, and yields q01 ≈ 1.95. In terms ofthe original variables this yields p1 ≈ 0.83 at x = 1 (from (2.87)), which agrees with thenumerical approximations in Figure 2.8.612.5.3 The Myosin ModelFor the full myosin transport model (A.24), the boundary-layer equations (2.82a)–(2.82b)can be studied qualitatively. In this section, I will explain that a boundary layer solutionexists near x = 0, however, there is no boundary-layer solution for the myosin model nearx = 1. This result will explain why the Type II QSS is not realized by numerical solutionsto the full myosin model as in §2.4.The boundary-layer near x = 0 can also be studied qualitatively in the phase-plane.Upon setting v1 = 1 and v2 = v, (2.82a) and (2.82b) on 0 < η <∞ becomep1η = −kbwp1p22 − p1 + kbp030, p1(0) = 0, p1 → p010 ≡1kbwα0, as η →∞,(2.90a)p2η = −1v(kbwp1p22 − p2), p2 → α0 as η →∞, (2.90b)where p030 = (α0 + 1/[kbwα0]) /kb and α0 = α(0, t). As before, introduce new variables q1and q2 defined byp1 =1kbwα0q1, p2 = α0q2, (2.91)so that in terms of r ≡ kbwα20, (2.90) becomesq1η = g1(q1, q2) ≡ −r(q1q22 − 1)+ 1− q1, q1(0) = 0, q1 → 1 as η →∞,(2.92a)q2η = g2(q1, q2) ≡ −1v(q1q22 − q2), q2 → 1 as η →∞. (2.92b)As in the last section, the equilibrium is a saddle point and the boundary layer solutioncorresponds to the stable manifold of the saddle point. At the equilibrium state q1 = q2 =1, the determinant of the Jacobian Jg of g1 and g2 is det(Jg) = (1− r)/v. Therefore,det(Jg) < 0 and q1 = q2 = 1 is a saddle-point if r ≡ kbwα20 > 1. Figure 2.19(a) depictsphase portrait of q2 versus q1 and nullclines for (2.92) for the representative values r = 5and v = 0.5. Observe that there is a unique value q2(0) = q02 for which (2.92) has a solutionwith (q1, q2) → (1, 1) as η → ∞. As such, there is always a boundary-layer solution nearx = 0 for the myosin model.A similar boundary-layer system near x = 1 can be obtained from (2.83) for the myosin620 0.5 1 1.5 200.511.52(a) r = 5, v = 0.50 1 2 3 4 5 6 700.511.52(b) r = 5, v = 0.5Figure 2.19: Qualitative analysis of boundary layer behaviour of the myosin model. Phaseportraits of q2 versus q1 for boundary layer solutions of the myosin model near x = 0 (a) andnear x = 1 (b) from (2.92) and (2.93), respectively. In (a) there is a unique value q2 = q02at q1 = 0 for which (2.92) has a solution with (q1, q2) → (1, 1) as η → +∞. However, forthe right boundary-layer, the phase-plane in (b) there is no value q1 = q01 > 0 at q2 = 0 forwhich (q1, q2)→ (1, 1) as η →∞.model. In place of (2.92), the system isq1η = −g1(q1, q2) ≡ r(q1q22 − 1)− 1 + q1, q1 → 1 as η →∞, (2.93a)q2η = −g2(q1, q2) ≡ 1v(q1q22 − q2), q2(0) = 0 q2 → 1 as η →∞, (2.93b)where r is now defined by r = kbwα21 with α1 = α(1, t). Although the equilibrium pointq1 = q2 = 1 is a saddle point of (2.93) whenever r > 1, the phase portrait in the q2 versusq1 plane in Figure 2.19(b) shows that there is no value q1(0) = q01 > 0 on q2 = 0 for which(q1, q2)→ (1, 1) as η →∞.As such, for the Type II QSS approximation (2.63) in the myosin model there is nosteady-state boundary-layer solution near x = 1 that allows the extra boundary conditionpB = 0 at x = 1 to be satisfied.2.6 DiscussionThe quasi-steady-state reduction method for molecular motor transport was introduced in[57] for reaction–advection–diffusion systems with linear reaction kinetics. Here, I havegeneralized this method to a class of problems where the kinetics are nonlinear, but wherea conservation condition is satisfied. The QSS method relies on the assumption that the63nonlinear kinetics occur on a faster time-scale than the diffusion and advection processes.In this limit of fast reaction kinetics, and under a condition on the eigenvalues of theJacobian of the kinetics, the full system dynamics were shown to be well-approximated bythe dynamics on a slow solution manifold, which consists of a single scalar quasi-steady-state PDE. This asymptotic formalism was used to analyze three specific nonlinear modelsfor the binding and unbinding of molecular motors.Three models I used as case-studies contained two distinct types of nonlinear reaction.(1) The kinesin model has a nonlinearity in the binding rate of motors to MT (due tosaturation, with and without binding cooperativity). This model reduces (with parameterc = 0) to the linear binding case considered in a previous study [16], and is used here asa basic “control” to validate our method. Typical nonlinear responses such as Michaelis-Menten or Hill function kinetics were used to describe the dependence of binding rate onthe free motor density (represented by the increasing and saturating function g). Here thenonlinearity was a function of a single state-variable. (2) In the second class of models,nonlinearity stemmed from interaction between motors in different states, such as collisionsthat lead to direction changes or stalling while bound to a MT. Both the kinesin-dyneincomplex model and the myosin motor model shared such aspects.Each model satisfied a conservation law, namely the total density of motors was fixed inthe cell (the total density was fixed at 1 for numerics throughout, as discussed in AppendixA.4). This constraint served an important purpose, as it was used to reduce the systemfrom n to n − 1 states (where n = 3 for all our models). In each case, the populationof motors in various states was defined in terms of one reference state (denoted by α(x)).The choice for that reference state was merely a matter of convenience of calculations, andspecific to each case.Many elements of the linear QSS theory carry over to the nonlinear analysis here. How-ever, the geometry of projections in the linear case (as developed in [7, 57]) no longer holds,suggesting that obtaining higher order terms in asymptotic solutions is no longer tractable.Obtaining expressions for such correction terms remains an open problem. Moreover, inmany cases, the diffusion coefficient in the unbound state is taken to be O(ε). If this is thecase, in those particular cases where the drift term vanishes, our QSS PDE would simplyreduce to a conservation law for the total density of motors, and fail to describe the dy-namics of the system. To avoid this, it is necessary to assume that the diffusion coefficientin the unbound state is O(1).For all such models, the QSS reduction of (2.15) leads to new scalar nonlinear PDEs,are not easily amenable to analytical solution techniques. Although it was still necessary tosolve these QSS PDEs numerically, the QSS reduction does effectively eliminate the smallparameter ε from the full model and avoids the more challenging numerical task of having64to compute solutions to the full nonlinear vector system (2.15) of PDEs at each small ε.The QSS analysis permits the formulation of conclusions about the overall rate of trans-port (advection velocity) of the system that results from the combination of motors walkingon MT, diffusing while unbound, and kinetics of binding, unbinding, switching directions,and/or stalling. Additionally, the QSS PDE was shown to provide insight into the be-haviour of the steady-state solutions as parameters are varied. This insight was used tointerpret cell-level behaviours resulting from various specific molecular-motor-level assump-tions. I now summarize some of the major conclusions and their implications for each ofthe case-studies.Kinesin modelHere the cytosolic motor state was used as the reference state α, and a Fokker-Planck (FP)equation (2.29) was derived for the total motor density. In the special case of spatiallyconstant microtubule bias, this reduced further to the FP equation (2.33a) for the cytosolicstate from which we can draw several conclusions. (a) The overall transport directiondepends on the sign of (1−2P ). (b) When (1−2P ) 6= 0 (which means that more MTs pointto one end of the cell than to the other), I predict an exponential spatial motor distributioncorresponding to MTs bias. (c) Both the effective diffusion and the effective transportrates are (essentially) averages of the diffusion and transport rates in the underlying states,weighted by the fraction of time spent in each of those states. These conclusions areconsistent with results of the linear models in [16]. (d) When MT polarity bias P (x) isspatially nonuniform, there arises the possibility for motors to pile up either at cell endsor in the middle of the cell, as shown in Figure 2.6. This reflects the earlier results forthe QSS reduction of a model with spatially varying parameters. In this case, the resultingQSS PDE had spatially dependent effective diffusion and velocity [58]. (e) The overall effectof nonlinear binding in this case is that more kinesin motors are sequestered in the freelydiffusing class, which results in a shallower motor density across the cell. The shallowersolution profile results from the fact that the binding rate is limited in both the saturatedbinding and Hill function binding cases. (f) Hill function binding (which could representcooperative motor binding interactions) creates ‘kinks’ and inflection points in the spatialmotor distribution, since the Hill function turns binding on or off more sharply than doesMichaelis-Menten kinetics.Kinesin-Dynein modelHere the nonlinearity involves a product of two state variables (left and right moving com-plexes), a composite left-right bias function Q(x), and possibly distinct velocities whenmoving right or left (see Appendix A.3.2 for the relationship of the function Q to the un-65derlying biological details). Here the left-moving motor variable was used as reference stateα. Both the effective transport rate and effective diffusion rate are “density dependent”functions of α. The effective transport rate depends intuitively on the model parameters.Increasing the velocity of left-moving complexes, decreasing the probability of binding tothe right-moving state, or increasing the right-to-left turning rate all result in biasing trans-port towards the left-end of the cell. The effective diffusion rate is scaled by 1/ka (ka is theassociation constant), which intuitively modulates how many molecular motor complexesremain in the cytosolic vs. bound states. The effective diffusion rate is further increasedfrom baseline through the “tug-of-war” that the motor complex exerts on its cargo. Thisincrease results from the product (1−Q)Q, which gives the probability of binding into theleft-moving and right-moving state. Although a motor cannot simultaneously bind into theleft-moving and right-moving state, I find that the competition between right-moving andleft-moving states increases the effective diffusion of the system—this makes sense, as anyrapid switching between right- and left-moving states is similar to a diffusive mechanism.Myosin modelThe motor interference was assumed to cause stalling with a higher-degree nonlinearity((pB)2pW) than in the kinesin-dynein motor complex model, which was inspired by thenonlinear interactions in a model for myosin aggregations [101, 102]. Moreover, the stalledand walking myosin motors have different velocities, with the stalled motors being trans-ported due to actin treadmilling. Interestingly, this higher-degree nonlinearity gave rise totwo distinct QSS solutions, one of which was characterized by the absence of stalled motors(pB = 0, “Type I QSS”). In this case, the QSS PDE is linear and the steady-state solutionwas found explicitly. For the second QSS solution with pB 6= 0, I identified a nonlinearFP equation with diffusivity D/(1 + kb), a density-dependent effective transport term, andan additional term proportional to (∂α/∂x)2. The latter (“Type II QSS”) exists only fora subset of parameters (Figure 2.11). Moreover, solutions to the full system converge tothe Type I solution, unless the full model is corrected by an asymptotically small diffusionterm for the stalled motors. Interestingly, such a term had been included in the model in[101, 102]. There, the inclusion of this small diffusion term was justified physically as a smallrandom motion of stalled motors, yet the analysis here reveals a mathematical justification.This peculiar effect stems from an issue with the boundary layer at the cell end x = 1.The small diffusive correction term changes the pB equation from hyperbolic to parabolic,allowing the model to be consistent with boundary conditions that the uncorrected modelcannot satisfy. In this case, the existence of two QSS solutions suggested further investiga-tion of the behaviour of the full myosin model. Through extensive numerical simulations, itwas possible to determine which QSS PDE would better describe the dynamics of the full66system (Figure 2.17). In the end, a phase-plane analysis of the non-spatial kinetic modellargely suggested which QSS PDE would be valid for a given set of spatially homogenousinitial conditions.All in all, the QSS analysis is generalizable to nonlinear models for molecular motors.That said, the examples discussed herein are simplified prototypes and caricatures of actualmolecular motor behaviour. For example, a caveat of the kinesin model (2.3) is that thenonlinear binding function, g(pU), may not accurately describe biological effects such ascompetition for binding sites on a single MT. As formulated with a saturating functionfor g(pU), the model implies that crowding in a region of the cell is responsible for limitingthe binding rate of motors to MTs, rather than explicit competition for binding sites. Inthe kinesin model, it is necessary to interpret the saturated binding rate as a result ofcompetition or crowding for binding sites on a single MT.In reality, many more states and interactions between states could occur, making thebiological system more realistic and interesting, but also much more complicated to analyzemathematically. In this analysis, I have not considered the cases of heterogenous multi-motor complexes composed of a distribution of motor types, nor the additional interactionswith cargo such as vesicles or early endosomes. It remains unclear at present whethersimilar methods would lead to insights in more realistic and complex models. The QSSmethodology has also been extended to 2-dimensional models in the context of a searcheralternating between ballistic and diffusive movement phases [6] with linear kinetics. Themethod presented here for 1D nonlinear models, should extend to two dimensional nonlinearmodels, provided the conditions on the kinetic terms are met, although it remains an openproblem for which classes of nonlinear kinetics and in which spatial dimensions it is possibleto analytically find an approximating QSS PDE.67Chapter 3Coupling Mechanical Tension andGTPase Signalling to GenerateCell and Tissue Dynamics3.1 IntroductionThe Rho-family GTPase proteins are central regulators within signalling networks of eukary-otic cells. While their effects extend to nearly all cell functions, a primary well establishedrole is to control the actin cytoskeleton, actomyosin assembly and myosin contraction [73].This fact makes Rho GTPases important in regulating cell shape in single cells and inepithelia. Rac1 promotes cell spreading by activating downstream signalling that leads toactin polymerization and cellular protrusions such as lamellipodia. RhoA activates differentdownstream signalling that in turn activates myosin-induced cell contraction. Hence, whileRac1 promotes cell spreading, RhoA counteracts this by stimulating cell contraction. Whileprevious studies have addressed how GTPases spontaneously segregate to front or back ina cell [54, 64, 95, 99], and how this leads to cell polarization and motility [25, 59], here Ifocus primarily on the effect of GTPase activity on cell contraction or spreading, and ontheir interplay with tension and mechanical forces experienced by cells.Rho GTPases cycle between active and inactive forms: they are activated by guaninenucleotide-exchange factors (GEFs), and inactivated by GTPase-activating proteins (GAPs)[73]. GTPases signalling proteins are interconnected, with crosstalk via a host of proteins.The proteins Rac1, RhoA and Cdc42 are central regulators, downstream of cell-surface re-ceptors that sense a host of stimuli, including small ligand gradients [74], adhesion molecules,extracellular matrix (ECM), substrate stiffness [18], as well as forces and mechanical tension[97].68It has been known for many years that mechanical tension can stimulate cells andlead to signal transduction, but details of the connections were poorly understood. Morerecently, techniques for measuring forces felt by cells [20, 77] have been used in coordinationwith methods for observing activity of GTPases [68]. This kind of experimental work hasrevealed a direct connection between mechanical tension and GTPase activity in cells. Forexample, Weiner and coworkers [35] showed that the aspiration (which increases tension)of a neutrophil membrane by a micropipette directly inhibits Rac1 activity. When tensionis released, Rac1 activity resumes in the cell. Compressing cells was shown to activateRhoA [30] in a rapid and reversible way. Isotropic stretching of vascular smooth musclecells on an elastic substrate was shown to inhibit Rac (timescale of 5 minutes, recovery over45 minutes) in [37]. The authors also quantified Rac activity versus % stetch, showing adecrease by about 50% in response to a 15% stretch (Figure 2B in [37]) How cells sensemechanical forces is reviewed in [20, 97], and the identity of multiple Rho GEFs and twoGAPs involved in mechanotransduction is summarized in [62].Cells have diverse mechanosensory mechanisms, and the molecular details of the linkbetween mechanical tension and GTPase activities are still emerging. Specific examples ofmechanosensory mechanisms include Rap1 as a tension-sensor and its effect on Rac1 [24] ortension-sensitive calcium ion channels that produce signals to the Rho GTPases [30]. Cel-lular adhesions and related structures (integrins, vinculin, and talin) act as mechanosensorsthat funnel signals to central regulators [41, 69]. Membrane tension is known to affect actinassembly directly by limiting polymerization and through a signalling pathway that inhibitsactin nucleation via a protein called WAVE2 [21]. Cell-substrate and cell-cell adhesion, cy-toskeleton, and their effects on Rho proteins is reviewed in [62]. Other proteins, such asmerlin, can act as a mechanochemical transducer by localizing to cortical cell-cell junctionswhen pulling forces are transmitted from cell to cell in epithelial tissue [15]. Focal adhesionkinase (FAK), for example, inhibits RhoA and activates it in response to tension-dependentintegrin reorganization, facilitating cyclic activation of RhoA and Rac1 [91]. Finally, specificproteins for sensing cell membrane curvature can also regulate the cycles of cell protrusionand retraction by controlling Rac1 through GAPs [92].The connection between mechanical forces and intracellular signalling is a two-waystreet. On one hand, mechanical tension can influence GTPase activity. On the other hand,GTPases lead to cell deformation (spreading or contraction) that exerts pulling, stretching,or contractile forces on the cell, the local extracellular matrix, and/or neighboring cells.This two-way feedback between chemical and mechanical signalling merits investigation,which is the main focus of this chapter. While mechanochemical interactions have beenconsidered in previous mathematical models for cell behaviours [34, 60, 61, 63, 65], to myknowledge, this is the first instance that mechanochemical interactions are applied to a69GTPase signalling cell and tissue dynamical system.Rho GTPase are embedded in complex signalling networks, with many effectors, in-terconnections, and inputs, but proteins such as Rac1, RhoA and Cdc42 play a centralregulatory role. Moreover, the details of such networks vary from one cell type to another,and adapt to cell state and environment. Several recent experimental and modelling stud-ies have provided evidence for the hypothesis that certain cell behaviours can be explainedas emergent properties of relatively small subsets of these networks, consisting of GTPasemodules. Examples of these simple modules include the bistability and hysteresis of cellshape [3, 14, 78] and cell motility behaviour [9], as well as diverse motility phenotypes inmelanoma cells on patterned adhesion surfaces [34, 65]. Using simple underlying models forGTPase “circuits” guides the approach here. Instead of attempting to describe the com-plexity of a large signalling network, which may vary from cell to cell, I restrict attention toa minimal GTPase signalling model and a simplified physical model for cell tension. Fromthis simple, conceptual model, a range of emergent behaviour can be explained. Furtherfine-tuning and adaptations to specific cells and experimental systems is left to the nextmodelling step, when mechanistic mathematical models can be studied in conjunction withexperiments.The first step is to consider a single GTPase, such as RhoA, associated with actomyosincontraction. I present a minimal model for RhoA activity, capable of bistable dynamicsand link it to feedback from mechanical tension. High RhoA activity leads the cell tocontract, which generically results in the reduction of tension from any applied stretch (Ido not assume a specific biophysical tension-sensing mechanism). I study this conceptualmechanochemical model in a single cell, without spatial effects. I characterize high or lowGTPase activities and transitions between these, and the coupled dynamics of cell tension.Owing to the simplicity of this two ordinary differential equation model, it is possible fullycharacterize parameter-dependence and delineate regimes of behaviour through numericalbifurcation analysis. In a the “single-cell GTPase-tension model”, there exist regimes of(1) high and (2) low RhoA (corresponding to contracted or relaxed cells) separated by (3)regimes of spontaneous, persistent cycling between these states which correspond to cyclesof contraction and relaxation in the cell.I then consider the dynamics of the minimal model when many cells are coupled togethermechanically in a 1D role or in a 2D epithelial sheet. In a collective, when one cell changesshape, forces are transmitted to neighbouring cells which are then transduced into GTPasesignalling. Despite the simplicity of the conceptual model, multicellular systems exhibit avariety of interesting behaviour. Tissues can contract or oscillate as a whole, or waves andspatially correlated dynamical patterns of activity and size can emerge in large 1D or 2Dtissues.70In the final step, I also consider a related GTPase circuit consisting of both Rac1 andRhoA (henceforth Rac and Rho) and mechanical tension. Mutual inhibition between thesehas been found in a number of cell types [46, 75] and highlighted in recent biological lit-erature for both normal and malignant cells [9, 28, 65, 67, 78]. The effect of such GTPaseinteractions on cell shape has been explored theoretically [34, 49], but the two-way feedbackbetween cell mechanics and cell signalling is the main theme that motivates the work herein.3.2 Minimal Model for a Single Mechanochemical CellFirst consider the simplest case, where the mechanosensitivity of a single cell affects itsGTPase activity (Figure 3.1(A)), which, in turn, affects a contractile actomyosin meshworkin the cell. The minimal model tracks the activity of a GTPase such as RhoA over time ina single cell. (While RhoA is known to redistribute intracellularly, I ignore spatial varia-tions within a cell, so as to build a first working multicellular model.) RhoA acts throughRho-associated protein kinase (ROCK) to phosphorylate myosin light chain, leading to ac-tomyosin contraction. Consequently, to capture the mechanical contraction, associate amechanical Kelvin-Voigt element (spring-dashpot system) with the cell size. In one dimen-sion, cell size is represented by a length, L (Figure 3.1(C)). A cell at mechanical equilibriumhas some constant “rest-length”, L = L0 (Figure 3.1(B)). To couple the signalling with themechanical tension, assume that cell tension, T , proportional to (L− L0), enhances RhoAactivation. Thus, if a resting cell is temporarily stretched, RhoA activity increases. Inturn, active RhoA results in contraction of the cell (I assume that active RhoA decreasesthe rest-length of the cell). As the cell contracts towards its rest-length, the effect of thetemporary increase in length is removed (resolving the tension). This reduces the GTPaseactivity to a lower level. The overall paradigm of the model is shown as a cycle throughstates following the purple arrows in Figure 3.1(B).3.2.1 Model Equations and DefinitionsFor the activity of the GTPase, I adopt the generic equationdGdt= (Tension-dependent rate of activation)Gi − (Rate of inactivation)G, (3.1)where Gi is the level of inactive GTPase. Ignoring spatial variation, and assuming thatthe total GTPase GT is roughly constant over the timescale of interest (GT = G + Gi =constant), leads to a single equation, (3.1), with Gi = GT − G. This equation is a directadaptation of the 1-GTPase spatially-uniform version of the model in [33], with the additionof the mechanical coupling.In the case of a linear equation (3.1), that is, if terms in braces are constant, no interest-71(B)(A)contractiontensionGiGL0resting cellexternal forceLstretched cellRho activated by tensioncontractioncontracted cellrelaxationRho inactivated(C)λ λkLFigure 3.1: The minimal model for coupled GTPase activity and cellular-tension. (A)Schematic of our minimal model for a GTPase “mechanochemical cell”. Typical GTPasecycling between active (G, orange) and inactive (Gi, blue) forms. Black arrows denote inter-conversion (solid), and positive feedback (dashed) from the active GTPase and from tensionto GTPase activation. Purple elements (in (A) and (B)) represent mechanical effects. Weassume that Gi = GT −G by conservation. Active RhoA results in cell contraction, whichreduces tension. Tension is assumed to increase the activation rate of RhoA. (B) A rest-ing cell (rest length L0, top left) is stretched by an external force to length L (bottomleft); the “spring” schematic represents contractile actomyosin). Tension T ∝ (L − L0)in the stretched cell activates RhoA (inset, lower right, color scheme as in (A)), leadingto a coupled mechanochemical system. RhoA activity results in actomyosin-powered cellcontraction, which eventually reduces cell tension as the cell approaches its new contractedrest-length. As RhoA is inactivated by the loss of tension (upper right), the cell relaxes.(C) Mechanical representation of the actomyosin cell cortex as a Kelvin-Voigt element.ing behaviour is found. Some feedback is needed to obtain the nonlinearities that generatebistability and allow for non-trivial dynamics. It is typical to assume positive feedback fromactive GTPase to its own activation [33, 54] (see [36] for the equivalence of other assump-tions). Furthermore, based on the prevalence of GEF-associated mechanotransduction [62],I include the tension-dependent feedback f(T ) in the activation rate. This leads to a modelequation of the formdGdt=(b+ f(T ) + γGn1 +Gn)(GT −G)−G. (3.2a)where b is basal activation rate (in the absence of feedback from mechanics) scaled by theconstant inactivation rate, and γ is a similarly scaled rate of feedback activation. (Detailsof the scaling are provided in the Appendix.) In this model, cell tension depends on the“size” L of the cell relative to its concurrent rest-length L0. I considered several forms of72f(T ), as described in the Appendix; however, I concentrate on the case thatf(T ) = β11 + exp[−αT ] , where T = L− L0. (3.2b)The parameter β governs the strength of feedback from tension to GTPase activation. Theso-called squashing function in (3.2b) means that the mechanical input has no effect if LL0, but builds up to a maximal level of β for L L0. Consequently, the model GTPase, G,is sensitive to a pulling force, but not to a squeezing or contractile force. It is straightforwardto generalize this minimal assumption to other mechanosensory mechanisms, but I onlyconsider the effect of tension here. The parameter α governs the sharpness of the GTPaseactivation response to cell stretching. It is worth remarking that this form of mechanicalmodel (3.2a) with (3.2b) bears a close resemblance to the equation proposed in [30] forthe dynamics of active RhoA in human fibrosarcoma cells that are exposed to mechanicaltension. The squashing function has the same basic property of switch-like activation asthe Hill-function dependence on T in [30].For the mechanical coupling, assume that GTPase activity (e.g. RhoA activating ROCK,which activates myosin light chain) effectively shortens the rest length of the “corticalactomyosin spring” promoting contraction. This is described through the following equationfor the cell size L:dLdt= −ε(L− L0), where L0 = `0 − φ GpGph +Gp, (3.2c)and ε = 2k/λ is the rate of contraction. This model assumes that the cell acts as anover-damped elastic spring with Hookian spring constant k, and viscous coupling to a fixedsubstrate (viscosity λ). The rest length, L0, is assumed to decrease from a fixed rest length,`0, depending on the amount of active GTPase within the cell, G. The dependence on G isrepresented by a Hill function with amplitude φ, half-maximum GTPase activity Gh, andpower p. For large GTPase activity G, the rest length approaches L0 ≈ `0 − φ, while forlow GTPase activity, the rest length remains near `0. A switch occurs close to the activitylevel G = Gh. The larger p, the sharper the transition between small and large L0 values.Equation (3.2c) presumes the over-damped regime, where inertial forces are negligible, asappropriate for modelling cell-scale behaviour. Equation (3.2c) follows from a force balanceat the two cell ends:λdx1dt= k(L− L0) and λdx2dt= −k(L− L0), where L = x2 − x1, (3.3)for x1, x2 positions of the left and right cell boundaries (in 1D), and fromdLdt =dx2dt − dx1dt .730.0 0.2 0.4 0.6b0.000.250.500.751.001.25Gcontractionrelaxation(a) Bistability in GTPase module.0.0 0.1 0.2 0.3 0.40.20.40.60.8L(b) GTPase-tension coupling results in new bi-furcations.Figure 3.2: Bifurcation diagrams for the minimal model (3.2). (a) On its own the GTPasemodel (3.2a) (with f(T ) = 0) exhibits bistability with respect to the activation rate b.(Other parameters: γ = 1.75, n = 4, GT = 2). Mechanical tension affects the GTPaseactivation rate, leading to the possibility of a relaxation oscillator (hysteresis loop) shownin this diagram. (b) Bifurcation diagram for the coupled GTPase-tension minimal modelEquation (3.2), showing how cell length L varies with the strength of coupling (β) of tensionto GTPase activation. L can be long (small β, solid blue line), oscillatory (middle values ofβ, magenta line), or short (large β, solid yellow line). In both, red points are saddle nodebifurcations, and the black point corresponds to a Hopf bifurcation (after which stableoscillations emerge). Other parameters are b = 0.1, γ = 1.5, GT = 2, φ = 0.75, Gh = 0.3,ε = 0.1, α = 10, `0 = 1, and n = p = 4.3.2.2 ResultsThe single GTPase model on its own (with β = 0), is bistable for a range of parameters.As the basal activation rate b increases, the system transitions from a monostable statewith low GTPase activity, through a bistable regime, and finally to a monostable state withhigh GTPase activity (Figure 3.2(a) and [33, 34]). With mechanical feedback (β 6= 0) asdescribed in Section 3.2.1, there are three regimes of behaviour: (1) for small β, the cellremains relaxed with low GTPase activity, (2) for large β, the cell becomes contracted withhigh GTPase activity, and (3) for intermediate β, the cell dynamics tends to a stable limitcycle with GTPase activity cycling between low and high levels. The bifurcation diagramof Figure 3.2(b) shows these three regimes of behavior, displaying steady state cell length,L, as a function of the coupling feedback strength, β. For this choice of parameters, thethree regimes of behavior occur for different intervals of β; however, for different parametervalues, it is possible that the limit cycle and contracted steady-state can both be stable forthe same value of β.740 50 100 150 200t0.000.250.500.751.00(a) β = 0.050 50 100 150 200t0.000.250.500.751.00(b) β = 0.20 50 100 150 200t0.000.250.500.751.00GLL0(c) β = 0.3Figure 3.3: Dynamics of the minimal model for a single cell with one GTPase (“RhoA”) andfeedback from tension to GTPase activation, Equation (3.2). In (a), the feedback strengthfrom tension (β = 0.05) is weak, and the cell remains relaxed. In (b), the feedback (β = 0.1)is of intermediate strength, and limit cycle oscillations arise. In (c), the coupling is so strongthat GTPase activity is always high, and the cell stays in a contracted state. Parametersare b = 0.1, α = 10, γ = 1.5, n = p = 4, GT = 2, φ = 0.75, Gh = 0.3, ε = 0.1 and`0 = 1. When the GTPase activity level is close to G = Gh, the cell rest length changessharply from L0 ≈ `0 to L0 ≈ `0−φ (green dash-dotted), resulting in the dramatic changesin cell length seen in (B). Some lag stem from the slower dynamics of L, due to the slowmechanical response (small parameter ε).The dynamics of the cell size (L, solid), rest-length (L0, dash-dotted curve) and GTPaseactivity (G, dashed curve) is shown in Figure 3.3 for each of these regimes. When thefeedback from tension upon GTPase activation is small, (a) β = 0.05, the cell remainsrelaxed. As the feedback parameter increases, (b) β = 0.2, the cell oscillates, or (c) forβ = 0.3, the cell contracts and maintains a small length.The results can be understood based on known dynamical systems behaviour of abistable system (the GTPase activity) with slow negative feedback (the mechanical con-traction). The coupling can constrain the bistable system to either its low or its highsteady state levels, or, for intermediate coupling, lead to a trajectory around a hysteresisloop. In the latter case, the system behaves as a relaxation oscillator due the separationof time scales between fast G and slow L (note the small parameter ε in the L ODE). Asshown by the hysteresis loop in Figure 3.2(a), the activation rate is increased when the cellis stretched, eventually leading to a transition from low to high GTPase activity. At thispoint, the GTPase activity leads to cell contraction, effectively decreasing the rest lengthL0. As the cell contracts, L approaches L0, and tension decreases, reducing GTPase acti-vation rate and transiting to the low GTPase state. This resets the rest-length to a largervalue. With the appropriate relative timescales of mechanics and chemical signalling, thiscycle repeats, setting up the limit cycle oscillations.75(A)(B)kλxj+2xj+1xj−1 xjFigure 3.4: Cell interactions in a 1D array of “model cells”. (A) The contraction-relaxationof each cell affects the force of pulling on its neighbours. Each cell has its own internalGTPase signalling. (B) The array behaves much like a system of over-damped springs inseries. The GTPase signalling affects the rest-lengths of the springs Lj , and the dynamicsthen moves the nodes xj that represent cell borders.3.3 Mechanical Coupling in a 1D Array of CellsHaving characterized the minimal “model cell”, next consider the behaviour of a coupledarray of such cells. As a first step, consider coupled cells mechanically in one spatialdimension (1D), as shown in Figure 3.4. Here the lengths of the cortical actomyosin Kelvinelements are simply the distances Lj = xj+1 − xj , j = 1, . . . , N − 1 between “nodes”(edges of cells along a 1D axis). Each cell has its own internal GTPase signalling, followingEquation (3.2a), and only responds to neighbouring cells through mechanical force. Hence,the motion of the cell ends, xj , is prescribed by the following system of ODE:λdx1dt= k(L1 − L1,0), (3.4a)λdxjdt= −k(Lj−1 − Lj−1,0) + k(Lj − Lj,0), (3.4b)λdxNdt= −k(LN−1 − LN−1,0), (3.4c)with j = 2, . . . , N − 1 giving the index j of the N − 1 cells. The rest-length in each cell,Lj,0, is coupled to GTPase signalling according to Equation (3.2c).3.3.1 Tissue Dynamics in 1D Depend on Mechanical Feedback StrengthWhen many cells are coupled together, new tissue-level behaviours emerge. For example, asone cell is displaced or contracts, its neighbours are stretched or squeezed. This change in76length then affects tension, T , and ultimately the GTPase activity, G of the neighbour(s),that can similarly affect their neighbours, and so on. The emergent behaviour depends onthe signalling parameters of the individual cells. For example, if the strength of feedbackfrom mechanics to GTPase activity, β, is sufficiently small or sufficiently large everywhere,the entire tissue will be relaxed (and long) or highly contracted (and short), respectively.Examples of these behaviours are shown in Figure 3.5(a) and (c).For β in the single-cell oscillatory regime (β = 0.2, as in Figure 3.3(b)) a small arrayof cells (N = 10) can exhibit synchronous oscillations, as shown in Figure 3.5(b). In thiscase, as each cell expands or contracts, the force exerted on its nearest neighbours inducesa change in the chemical signalling, which results in the coordination of the entire group(possibly excluding the cells at either end). This shows up as coherent bands of colour inFigure 3.5(b) while the total length of the array (vertical dimension in Figure 3.5) oscillates.I next asked whether a propagating wave of contraction, similar to that obtained in thework of Odell et al. [61] (as discussed in Chapter 1), could be obtained with this model.To explore this, I set up a simulation in which each cell would sit in the relaxed-lengthsteady-state, but provide an initial GTPase perturbation to one or more cells at the end ofa tissue. I was unable to reproduce a wave of contraction in this way. However, using amodification of the model with linear feedback (instead of the squashing and hill functions)between GTPase signalling and tension, I was able to produce wave-like behaviour. Thespecific parameters and feedback functions for this linear-feedback model are collected inSection B.3. The wave of contraction behaviour is shown in Figure 3.5(d) where two cellsat the right end of the row are initially “stimulated” with high GTPase activity, whilethe rest of the cells are at their relaxed steady-state. Contraction of the stimulated cellsstretches their immediate neighbours to the left, which activates new GTPase signalling inthose neighbours and subsequent contraction. In this way, a unidirectional wave of GTPaseactivity and contraction sweeps across the entire row of cells. As discussed later, this waveof contraction resembles a wave associated with zippering in the neural tube closure of anascidian embryo [29].As the number of cells increases, the spatial extent of the mechanical force transductioncan no longer span the entire “tissue”, and appears to become localized to some neighbour-hoods. Then, patches of contraction and relaxation emerge; these can propagate throughoutthe tissue as waves of contraction-relaxation. Typical examples for 50 and 100 cells in such1D arrays are shown in Figure 3.5(e) and (f). In such large arrays, GTPase activity is alsoseen to form wave patterns that sweep back and forth across the 1D domain. This leadsto the slanted bands of colour in Figure 3.5(e) and (f). Even though the GTPase activityis not directly coupled between cells, the mechanical coupling effectively leads to GTPasecoordination on some spatial and temporal scale.77(a) β = 0.05, 10 cells (b) β = 0.2, 10 cells(c) β = 0.3, 10 cells (d) Wave of contraction, 14 cells(e) β = 0.2, 50 cells (f) β = 0.2, 100 cellsFigure 3.5: 1D tissue dynamics result from mechanochemical interactions. Kymographsshow the 1D position of cell edges (vertical axis, black curves; suppressed for clarity in (e),(f)) with colour indicating the GTPase activity within each cell. Parameters as in Figure3.3(b), except in (d). For a small number of cells (N = 10), the tissue can be (a) relaxed,(b) oscillatory, or (c) contracted. (d) An initial perturbation of GTPase activity at one endof the row can propagate a wave of GTPase activity and contraction throughout the wholerow of cells. Larger number of cells: (e) N = 50, (f) N = 100: waves of contraction andrelaxation propagate across the tissue and model details in Section B.3.78I did not explore all possible behaviours of this 1D cell collective, and an analytical ornumerical analysis of the behaviour of the cell collective is warranted to delineate possibleregimes of behaviour. Nonetheless, the examples provided in this section illustrate aninteresting range of behaviour. A limitation of the 1D approach is that each cell onlyinteracts with 2 neighbours and cell shape does not play a role. As such, in the nextsection, I study a realization of the GTPase-tension model in a 2D setting.3.4 Cell shape and cell-cell interactions in 2D epithelialsheets3.4.1 Adapting the ModelIn order to describe cell expansion and/or contraction in 2D, I modified the model torepresent changes in projected cell area, A, rather than cell length. Generalizing from the1D model, assign a “resting cell area” A0 to the cell, and assume that positive (A − A0)corresponds to an average cell-stretching tension that has an effect in 2D similar to (L−L0)in the 1D model cell. This assumption could be modified, scaled according to A ≈ cL2,or adapted to experimental data. In the context of the simple conceptual model, the maineffect, preserved by these assumptions, is that GTPase activity and mechanical tensionswitch one another on or off.To simulate cell shape and intracellular chemistry, I worked with Dhananjay Bhaskar touse a publicly available software package, CompuCell3D, that represents cell shapes usingthe Cellular Potts model (CPM) formalism [90] as introduced in Chapter 1. Briefly, thepixel-based motion of a cell edge outwards (expansion) or inwards (contraction) is governedby a Hamiltonian, H, describing the total energy in the system. The Hamiltonian includesadhesion energies, and volume constraints (area constraint in 2D). At each time step (calleda Monte-Carlo step (MCS)) in the simulation, several small changes are introduced, calledpixel-copy or spin-copy attempts. The CPM algorithm accepts such changes if this decreasesthe Hamiltonian (overall energy of the system), or accepts it randomly otherwise as a smallnoise-induced fluctuation. While CPM does not explicitly track forces, it has recentlybeen shown to be consistent with other simulations where forces are made explicit [47], forexample, vertex-based cell models.Several aspects of the simulations were adapted to the technical requirements of theCellular Potts model (CPM). The timescale τ and notion of a “target area”, AT , wasintroduced. The actual cell area A and the GTPase-governed target cell area AT are trackedin each CPM cell. The target area for a cell is determined by a system of ODEs that couplesub-cellular biochemistry (assuming that the cell is well-mixed) to cell mechanics. This79leads the following model equations for the GTPase activity G and target area AT :1τdGdt=(b+ f(T ) + γGn1 +Gn)(GT −G)−G, (3.5a)1τdATdt= −ε(AT −A0(G)), where A0(G) = a0(1− φ GpGph +Gp). (3.5b)Here, a0 is the constant baseline cell area. The target area AT approaches A0 on thetimescale τ that can be controlled to increase or decrease the speed of the feedback. Notethat cell area A does not approach target area AT instantaneously, but through the additionor removal of lattice sites over several MCS. That is, A, is updated stochastically to approachthe target area AT by the CPM. Here, I also assumedf(T ) = βAmAmT +Am, where T = A−AT . (3.5c)Tension is defined as T = A−AT , which is a “delayed” form of A−A0. In turn, the functionf(T ), describes the feedback on GTPase activation from tension. This Hill function has theproperty that as m increases, its shape is fundamentally similar to that of the squashingfunction used in the 1D GTPase model (Equation (3.2b)). For this exploration, I againassume that GTPase activity is uniform inside a given cell through variable across theentire collection of cells. Stochasticity in the CPM leads to interesting behaviour (e.g.,stochastic switching) which is not observed in deterministic numerical solutions.As before, assume that increasing tension (represented as the difference between targetarea and actual cell area), can increase GTPase activity via Equation (3.5c). To appropri-ately calibrate the CPM to observe the same oscillatory dynamics as in the one-cell singleGTPase model, it is necessary to choose the timescale and the time step for numericallyintegrating the ODEs, τ and ∆t, respectively, so that each MCS is τ∆t = 2000 · 0.001 = 2units of time t. In the next section, I consider the GTPase-tension model in a single 2-dimensional cell, and later consider the coupling between interconnected cells.3.4.2 Single Cell DynamicsWith appropriate calibration, the 2D Cellular Potts model (CPM) implementation recapitu-lates the behaviour found in the single cell model. As shown in Figure 3.6(a) and AppendixFigure B.10, a parameter set corresponding to 1D cell oscillations also led to single-celloscillations in the 2D CPM cell. The CPM also produces relaxed cells and contracted cellsfor corresponding parameter sets in the 1D model, see Appendix Figure B.8 and B.9 forrelaxed cells.Cellular Potts model simulations have inherently stochastic behaviour due to the al-80(a) Single cell oscillations, β = 0.2.(b) Single cell, stochastic oscillations, β = 0.5Figure 3.6: Single cell oscillations in 2D cells simulated with the Cellular Potts model (CPM)CompuCell 3D software [90]. In (a) and (b), cell color represents a (spatially uniform)GTPase activity level from low (blue) to high (yellow and orange), as shown in the colorbar. Cell shape changes over time as indicated by the progression of snapshots numbered bythe Monte-Carlo step (MCS) of the CPM (MCS increase left to right and top to bottom).Cell target area (green) and actual area (blue) as well as GTPase activity is plotted over250 MCS of the CPM. In (a), β = 0.2 results in a single oscillatory cell. In (b), the cellstochastically switches between high GTPase steady state (corresponding to β = 0.5) anda large amplitude limit cycle. Other parameters were: τ = 2000, b = 0.1, m = 10, γ = 1.5,n = p = 4, GT = 2, ε = 0.1, a0 = 400, φ = 0.75, and Gh = 0.3.81lowable random fluctuations mentioned above. As a result, I found new behaviour thatwas not found in the deterministic 1D cell simulations, namely that spontaneous cyclesof high to low GTPase level (and low to high cell areas) could occur, even in parametersets consistent with monostable states. An example of this type is shown in Figure 3.6(b).Here, parameter values were set to the stable small-size single-cell regime in the 1D model(β = 0.5). The cell was in a contracted state for some time, but displayed two cycles ofcontraction-relaxation, at MCS 100 and 150, before returning to its quasi-quiescent state.In Appendix Figure B.11, there is an example of cells switching between small and largelimit cycle oscillations for β = 0.175.3.4.3 Coupling CPM Cells in 2DWhat happens when there are multiple cells in the 2D simulation? To answer this question,consider two types of situations in which N cells are present, where each cell is governedby its own set of 2D equations (see Apppendix) with the same set of parameters but withrandom initial conditions. As shown in Appendix Figure B.2 for N = 9 cells, the firstimplementation was of cells that have no direct mechanical coupling. As expected, in thiscase, cells behave independently with distinct and uncorrelated copies of the dynamics. TheRho GTPase levels inside such cells (top, Appendix Figure B.3) remain unsynchronized,as detected by the Kuramoto order parameter and the variance in the phase (details inAppendix Section B.5).Next, consider cells that were contiguous and can interact via adhesion terms in the CPMHamiltonian (details in Appendix Section B.4). Essentially, cells that have larger interfaceswith their neighbours have stronger adhesion (and lower adhesion energy). An example ofsuch simulations are shown in Appendix Figure B.4. As a cell changes size, neighbouringcells are affected through cell-cell adhesion. As one cell area contracts, its neighbours arestretched, causing their tension, proportional to (A − AT ), to increase. This promotes aneighbour’s GTPase activity, and leads it to contract. In this way, mechanical forces arepropagated throughout the tissue and affect GTPase signaling in each cell. As seen inAppendix Figure B.5, GTPase activities rapidly synchronize in the entire group of 9 cells,with a few small fluctuations in phase seen occasionally. In collaboration with DhananjayBhaskar, we used the Kuramoto order parameter to quantify the degree of synchronizationbetween the 9 cells. The Kuramoto order parameter is a complex number whose magnitudemeasures the phase-coherence of oscillators and can vary between 0 and 1. When oscillatorsare close in phase, the Kuramoto order parameter is closer to 1 ([42], [89] for a review).The larger Kuramoto order parameter and lower variance in the phase (note scaling on thevertical axis) as compared with simulation for independent cells (Appendix Figure B.3) alsoconfirms this synchronization.82The strength of cell-cell adhesion can affect the results. As shown in the sequenceof Appendix Figure B.5-B.7, as cells become more adhesive to one another (“low cell-celladhesion energy”) than to the surrounding “medium”, the mechanical coupling is stronger,and the synchronization of cell oscillations is more regular.3.4.4 Waves of Contraction and GTPase Activities in 2D Model TissueNext, again with Dhananjay Bhaskar, we asked how larger numbers of cells, also in 2D CPM,would behave when coupled mechanically through their adhesion. To probe this question,consider a simulation with a circular tissue composed of 373 contiguous cells with initialareas randomly chosen. As before, parameters of each cell are in the oscillatory single-cellregime. Results are shown in Figure 3.7 for the case of intermediate adhesion, in Figure 3.8for the case of strong adhesion, and in Appendix Figure B.12 for the case of weak adhesion.Here the 2D tissue is much larger than a few cell diameters. Figures 3.7 and 3.8 show twoviews of the same “tissue”, one (a) indicating cell area on a colour scheme of blue (small)to yellow (large), and a second (b) representing the concurrent GTPase activity level fromlow (blue) to high (orange).In contrast to the case of few cells, where synchronized oscillations were observed overthe entire population, synchronization is limited to patches in larger tissues. Moreover,waves of contraction/relaxation and GTPase activity propagate throughout the tissue. Thisbehaviour can be seen in the successive snapshots in Figure 3.7(a) and (b) for the case ofintermediate cell-cell adhesion. Bands of highly contracted cells (dark blue) are noteworthyin several panels in Figure 3.7 (a), and coincide with interfaces between zones of high andlow GTPase activity in Figure 3.7 (b).The strength of cell-cell adhesion affects the strength of coupling and extent of synchro-nization. In the case of weak cell-cell adhesion (Appendix Figure B.12(a) and (b)) relativelysmall patches are seen, and cells tend to detach from the periphery of the tissue. Waves ofcontraction and expansion are observed. As cell-cell adhesion is increased from the baselinesimulation in Figure 3.7, cells are more likely to favour adhering to each other. They thenexperience larger changes in area as one of their neighbours shrinks or grows. This resultsin nearly the entire tissue of cells expanding and contracting, though we still tend to see awave of synchronization spreading from the centre to the edge of the tissue, as in Figure 3.8.This leads to the conjecture that the patch size (number of cells in a group with coordinatedbehaviour) increases with the strength of cell-cell adhesion.In the final CPM experiment, consider the case that cells are heterogeneous, with arange of values of the feedback parameter β coupling mechanics to the GTPase activation.Consider a large simulated tissue with values of β assigned randomly to each cell. Resultsare shown in Figure 3.9. With the range of values of β, individual cells could be either83contracted, oscillatory, or relaxed. Due to the presence of some oscillatory cells in the tissue,those cells which would normally be quiescent at either relaxed or contracted steady-states,undergo oscillations due to pulling by oscillatory neighbours. Patches of activity in thetissue persist, though with a heterogenous GTPase activity. These forced oscillations aresuggestive of a mechanism by which tissue dynamics can be driven by a few pace-maker cells,whose phenotype is oscillatory. I discuss how these model predictions relate to epithelialdynamics in biological systems in §3.6.3.5 Rac and Rho GTPase ModelIn this section, I determine whether some of the lessons learned from the single-GTPasemodel would carry over to similar conclusions in a slightly expanded Rac-Rho GTPasecircuit. It is well-known that Rac1 and RhoA are mutually inhibitory under many situations[9, 67, 78]. Here, the analysis starts from the well-mixed variants of the Rac-Rho modeldescribed in [33], and in the melanoma-based modelling of [34]. In the latter case, couplingof front and rear compartments of a cell (through extracellular matrix signalling) was foundto lead to the possibility of distinct behavioural regimes, including stable high Rho or Rac,or cycling between those levels. Here the mechanical coupling has an effect similar to theECM coupling in that paper.In the mutually inhibitory Rac-Rho model, the total level of Rac (Rho) GTPase, RT(ρT ) is assumed to be roughly constant over the timescale of interest. Hence only the activeforms of the GTPases need to be tracked. Assuming that each of Rho and Rac inhibits theactivation of the other, consider the set of equationsdRdt=bR1 + ρn(RT −R)− δR, (3.6a)dρdt= (bρ + f(T ))11 +Rn(ρT − ρ)− ρ, (3.6b)where f(T ) represents the activation of Rho GTPase by tension T as in (3.2b). As before,Rho GTPase activity decreases the rest-length of the cell:dLdt= −ε (L− L0) , where L0 = `0 − φ ρpρph + ρp. (3.6c)Note that while Rac is a candidate for mechanosensory inputs and also has an affect on cellsize, I initially assume that only Rho is affected by tension and has an affect on cell length.On its own, without feedback from mechanics, the minimal Rac-Rho mutual inhibitionmodel has a region of bistability, as shown in Figure 3.10(a). As either the basal activationrates, bρ or bR, increase, the system transitions from a monostable state with either species84(a) Cell area in a 2D tissue over time.(b) GTPase activity in the same 2D tissue over time.Figure 3.7: Simulation of a 2D “tissue” (N = 373 cells) in the intermediate adhesion scenariousing CompuCell3D [90]. Individual cells satisfy the minimal GTPase-tension model, withT ∝ (A − AT ), where A is cell area, and AT is the target area. Cell-medium adhesionenergy (80) is equal to cell-cell adhesion energy (80) in the Hamiltonian, H. In (a), cellsare coloured based on their current cell area, while in (b), cells are coloured based on theuniform level of GTPase activity within each cell. Cells with smallest area (dark blue in(a)) are correlated with an interface between high (orange) and low (blue) GTPase activityin (b). Waves of contracting cells and relaxing cells are observed throughout patches in thetissue. Parameters listed in Appendix Section B.4 and B.6.85(a) Cell area in a 2D tissue over time.(b) GTPase activity in the same 2D tissue over time.Figure 3.8: As in Figure 3.7 but for the strong adhesion scenario. Cell-medium adhesionenergy (80) is greater than cell-cell adhesion energy (30) in the Hamiltonian, H. The entiretissue is synchronized. In (b), cells are coloured by area, while in (b), cells are coloured byGTPase activity. Notice that cells at the outer edge are first to expand/contract as they areless constrained by neighbours, so that expansion/contraction is ‘outside-in’. Parametersare the same as in Figure 3.7, and are listed in Appendix Section B.4 and B.6.at a high steady-state (and the other at a low steady-state) or from a coexistence state, intothe bistable regime. Assume, as before, that stretching a cell would increase the activationrate of RhoA. With that assumption, the same regimes of behaviour as in the single GTPasemodel (Section 3.2.1) occur in the Rac-Rho model. These three regimes of behaviour dependon the strength of feedback from tension to Rho activation (a parameter denoted γρ). Thedependence is shown in the bifurcation diagram in Figure 3.10(b). For small γρ, the cellremains long and relaxed, with high levels of Rac activity and low levels of Rho activity86Figure 3.9: As in Figure 3.7 (b), but with the parameter β (feedback strength from tensionto GTPase activation) initially randomly chosen for each cell. Cells are coloured by GTPaseactivity. Cells in steady state are forced to oscillate due to mechanical coupling with cellsthat are in the limit cycle regime. In this case, the baseline area parameter is increaseda0 = 600 (resulting in larger variation in cell area), and temperature parameter of Pottsmodel T = 15 is decreased from baseline. In the Hamiltonian, H, cell-cell adhesion energyis 60, and cell-medium adhesion energy is 80. Other parameters are as in Appendix B.4and B.6.(Figure 3.11(a)). For large γρ, the cell is contracted, (small L) with low levels of Rac activityand high levels of Rho activity (Figure 3.11(c). For intermediate γρ, limit cycle oscillationsarise (Figure 3.11(b)). There is a regime of parameter space where a stable limit cycle andstable steady-state coexist (approximately 12.29 ≤ γρ ≤ 15.61). In this parameter regime,depending on initial conditions, the cell may either end up in the oscillatory regime, or atthe contracted cell state.In the above Rac-Rho model, I considered only coupling between mechanical tensionand Rho activity, ignoring possible direct effects of mechanosensing on Rac activity. Rac isknown to cause cell spreading via actin assembly, an effect that I had similarly omitted. Tocheck the possible outcomes of such additional factors, I briefly explored several variantsof the above default Rac-Rho-tension model. Specifically, I experimented with inclusionof (1) the effects of compression (as opposed to tension) sensing with feedback to Racactivation and (2) the effect of Rac activity on cell size, modelled as an increase in therest-length L0. Feedback from mechanics to Rac activation can be interpreted as a changein bR in Figure 3.10(a). This can push the underlying Rac-Rho signalling model into a870 10 20 30bρ0102030bR(a) Bistability (inside curve) in the Rac-Rho model without tension.10 15 20 25 30ρ0.20.40.60.8L(b) GTPase-mechanics coupling resultsin a range of possible dynamics.Figure 3.10: Bifurcation diagrams for the minimal Rac-Rho model of Equations (3.6). (a)On its own, the Rac-Rho model (with β = 0) exhibits bistability (inside red-borderedregion) with respect to the activation rates bR and bρ. Mechanical tension affects the RhoGTPase activation rate, leading to the possibility of a relaxation oscillator by traversingthe bistable region (grey arrows). (b) Bifurcation diagram for the coupled Rac-Rho-tensionminimal model (3.6), showing how cell length L varies with the strength of coupling (γρ)of tension to Rho activation. L can be long (small γρ, solid black line), oscillatory (middlevalues of γρ, magenta curve), or short (large γρ, solid black line). As opposed to the singleGTPase-tension model, it is possible for a stable limit cycle to coexist alongside a stablesteady-state (for 12.29 ≤ γρ ≤ 15.61). Here, the red points correspond to saddle node (fold)bifurcations, and the black point to a Hopf bifurcation. Other parameters are bR = 15,bρ = 5, RT = ρT = 4, δ = 1, n = p = 3, γL = 0.75, ε = 0.1, α = 10, ρh = 1 and `0 = 1.regime of different behaviour—high Rac, bistability, coexistence, or low Rac—and alter theresulting cell behaviour accordingly (in a mechanical feedback-dependent manner). In theseadditional numerical experiments, behaviour similar to Figure 3.11 arises; albeit with Racactivity increasing the cell length and the relaxation oscillation arising from compressioninstead of tension.Aside from the above complementary Rac-feedback-only model, I guided an undergrad-uate advisee, Jim Shaw, to experiment with mixed Rac and Rho feedbacks and antagonisticeffects on cell size. With Jim Shaw, I also considered feedback from tension and/or com-pression to GAPs as well as GEFs (inactivation versus activation terms in the GTPaseequations). Overall, similar regimes of behaviour can be found in many such examples,within smaller or larger regions of parameter space. Cases of specific interest should hence-forth be linked to specific biological examples where the mechanical coupling to knownGEFs or GAPs is of interest.880 50 100 150 200t01234(a) γρ = 100 50 100 150 200t01234(b) γρ = 140 50 100 150 200t01234ρRL(c) γρ = 18Figure 3.11: Dynamics of the model (3.6) for a single cell with two GTPases (“Rac1” and“RhoA”) and feedback from tension to Rho activation. In (a), the feedback from tension(γρ = 10) is weak, and the cell remains large and relaxed (L ≈ 1) with high Rac and lowRho activities. In (b), the feedback (γρ = 14) is of intermediate strength, and limit cycleoscillations arise, provided that the initial conditions send the system to the stable limitcycle, instead of the contracted-cell steady-state (see Figure 3.10). In (c), the coupling isso strong (γρ = 18) that RhoA activity is always high, Rac is low, and the cell stays in acontracted state (L 1). Parameters are as in Figure 3.10.3.6 DiscussionFeedback between biochemical signalling and mechanical forces plays a vital role in devel-opmental biology and morphogenesis. Given the increasing biological evidence for the rolethat mechanical forces play in signalling networks, such as GTPases, mathematical andcomputational approaches are relevant and important to elucidate behaviours and suggesthypotheses. For example, a recent review [27] highlights how diffusion-driven patterns,differential adhesion, buckling instabilities in growing layers, and flows in active materials(cytoskeleton and motor proteins) lead to patterning. Following the experimental work of[35], a single-cell model was also developed to describe the inhibition of cell polarization bymembrane tension [98]. The authors used a more sophisticated spatio-temporal model ofthe GTPases in a 2D cell (based on the idea of wave-pining [54]), its downstream effect onactin, and a cell boundary represented using the phase-field method. The model was ableto account for observations on how build-up of tension in a neutrophil (by aspiration intoa micropipette) and sudden release of tension (by severing a long cellular protrusion) affectthe level and distribution of GTPase activity.While the effect of tension on GTPase activity was studied previously [30, 35], this is thefirst model that links GTPase-induced cell contraction to tension-induced GTPase activityin single cells and in a 1D and 2D tissue. Interestingly, a model based on mechanochem-ical coupling of some (indeterminate) signalling chemical and cell length that was studiedmathematically and computationally decades ago [60] bears resemblance to the model pre-sented here. That previous model was aimed at understanding folding and invagination89of epithelia, for example in the process of gastrulation [61, 63]. It was shown then thata localized stimulus in one cell could result in active localized contraction in some neigh-bourhood, creating the first fold in an early embryo. As GTPases and their effects wereyet to be characterized, this early modelling work was theoretical and speculative. Now,such work has additional relevance within the context of mechanical tension and GTPasesignalling and illustrated herein. More recently, single and collective cellular oscillationswere accounted for by a generic oscillator model for turnover of force-producing material(such as myosin motors) contracting against an elastic element [17]. Similar to the resultshere, varying the mechanical and kinetic properties of the system can transition the cellbehaviour from relaxed cell, to oscillatory, or to contracted cell length and collective cellbehaviour from unsynchronized to synchronized oscillations [17].Through coupling a simple GTPase bistable model without spatially distributed activitywithin a cell to a simple elastic (Kelvin-Voigt element) cell, I found three distinct regimesof behaviour, including high and low GTPase activity, with coordinated cell tension, andpersistent periodic cycling between those states. Here the dynamic pattern of contractileactivity stems exclusively from cell size fluctuations, amplified by tension-dependent GT-Pase activity. I did not include chemical diffusion (each cell is assumed to hold a uniformGTPase level), nor explicit cytoskeletal flows.Since actual signalling networks are incredibly complex and intricate, it is a significantchallenge to understand how cell behaviour can emerge from underlying components andproperties. Nonetheless, large networks have been studied theoretically, e.g., by [32, 39, 50].For example, Boolean models of cell signalling including tens of interacting species have beenused to show oscillatory activity of Rac and Rho GTPases [32]. Given this level of detail, itis easy to understand specific biological mechanisms yet difficult to explain the connectionsbetween network features and parameters on one hand, and emergent cell behaviour on theother. For this reason, stripped-down conceptual models that concentrate on key topologiesand regulators have a role to play in the theoretical understanding of cell behaviour [93].This principle motivates the analysis of small models for the GTPases.Assumptions made for the purpose of simplification can be modified substantially with-out changing the overall conclusions. For example, while our equation for Rho activationresembles that of [30], the authors’ assumption about the Hill function-dependence on ten-sion can be modified to another switch-like function that turns on at some critical forcemagnitude. Furthermore, while there is so far evidence for the multiplicity of GEFs involvedin mechanotransduction (relative to GAPs; see [62] Table 1), the model works equally wellwith GAP-sensitive responses as with the GEF-based GTPase response assumed here.One of the key findings is that simple coupling between GTPase activity and tension isconsistent with a range of biologically-relevant cell behaviour. The simplest model already90produces contracted or relaxed cells as well as cyclic fluctuations between these states.In the Rac-Rho circuit, such dynamic oscillatory regimes can coexist with static steadystates in the same parameter range, highlighting the dependence on stimuli and/or initialconditions. Oscillations in cell size are observed under laboratory conditions in epithelialmonolayers [103, 104]. While the link between tension and GTPases may be just one factoroperating in such systems, the model suggests experiments that could be used to test theconnection. In particular, inhibitors of ROCK that would abrogate the connection betweenRhoA and actomyosin contraction, or of Rac that would inhibit the antagonism of Rac toRho could be used to test the effect on the presence, frequency, and synchrony of cell volumeoscillations.Hashimoto el al. investigated the “zippering” in the neural/epidermal boundary of thesea squirt (Ciona intestinalis) embryo, part of the process that sets up neurulation overa time frame of about 2 hours [29]. Zippering involves successive shortening of cellularjunctions, one after the other, in a unidirectional wave of contractions up the zippering axis.The contraction was shown to be powered by the localization of active myosin, along theboundary, and to be dependent on Rho GTPase activity [29]. Furthermore, their paper wasaccompanied by kinematic simulations that reproduced the sequence of contractions, basedon assigned tensions (for pre- and post-contraction cells) and assigned time intervals. Here,I propose a simple model that aims at closing the gap between kinematics and dynamics.Briefly, instead of the manual assignment of forces and time intervals, a closed-loop chemical-mechanical system can give rise to the wave-like pattern of sequential cell contractions. Asshown in Figure 3.5(d), under suitable conditions, a unidirectional wave of Rho activationcontraction is supported by the minimal model. That said, while the conceptual model herecan account for the formation of wave, it is clear that other factors such as communicationbetween cells reaching across the “zipper” sides play vital roles not considered here (EdwinMunro, personal communication).In other organisms, such as Drosophila, Rac GTPase is known to have multiple rolesin early morphogenesis [100]. During Drosophila dorsal closure, over- and under-expressingRac results in the excess assembly of lamellipodia or disrupts the assembly of an actin cable(and subsequent zippering) and cell protrusions. While GTPases such as Rac regulate cellbehaviour during these morphogenetic processes, it is likely that cell and tissue mechanicsalso play an important role. Upstream mechanical signalling to Rho GTPases may occuras cells move and forces are transmitted, or as cell-cell junctions are rearranged. In thecase of zippering in sea squirt embryos, or in Drosophila dorsal closure, further validationof the mechanisms and/or completion of other essential elements remains as a future step.Nonetheless, with these mechanosensing assumptions, it is possible that feedback betweensignalling and mechanics can account for diverse single and collective cell behaviour in91these morphogenetic processes. Extending the conceptual model here to specific organismsby connecting to GTPase signalling and cell mechanics to data from experiments remainsa direction for future work.I focused on cell size (expanded or contracted), but it is also of interest to considerhow cell polarization is affected by mechanical cues in isolated cells and in cell collectives.See [43] for some theoretical background and review. Importantly, the results point toparameter regimes in which cells oscillate between compression and relaxation (in the caseof a single Rho-like GTPase), or compression and stretching (for Rac-Rho). But it is knownthat such cyclic stretching can itself change the properties of cells, reorganizing stress-fibers,for example in a Rho-dependent manner in endothelial cells [38]. It would be of interestto explore such polarity and directionality in future 2D models of this type, as well as toconsider how the feedback between GTPases and tension operate in collective cell migration[71]. There is evidence that GTPases also affect the cell-cell adhesion [52, 96] and tight-junctions [105], which would affect the coupling of mechanical transduction between cellsin a tissue. This could be of interest in future models.92Chapter 4ConclusionsBroadly, this thesis has addressed how biological mechanisms combine to organize cell be-haviours using models and techniques from applied mathematics. It is possible to gain anunderstanding of mechanisms, facts and theories about how a host of biological playersinteract with each other from biological experiments and observations. Building on theseobservations, multi-scale modelling and analysis can serve as a platform for understandinghow cell-scale and tissue-scale organization emerges from the interactions of the many bi-ological players. Along these lines, I have discussed two examples in this thesis. As thefirst example, I extended asymptotic quasi-steady-state (QSS) reduction methods to a classof nonlinear reaction-advection-diffusion PDE systems satisfying a conservation conditionwhich describe molecular motor transport within cells. As the second example, I developeda model to explore the interplay between GTPase signalling and cell mechanics, which couldgenerate a wide range of single and collective cell behaviour.In this section, instead of summarizing the results from each part (results are summarizedand discussed in sections 2.6 and 3.6, respectively), I will comment on the significance,broader contributions, and future applications of each part.Intracellular Transport by Molecular MotorsThe contribution of the work in Chapter 2 was the extension of existing QSS methodsfor molecular motor transport to include nonlinear reaction kinetics. These methods werepreviously limited to models with linear reaction kinetics, but were successfully appliedto understand molecular motor based transport in neurons [57] and in fungal hyphae [16].Here, I have illustrated that even with the incorporation of nonlinear reactions, which maybetter represent biological interactions, the QSS approximation methodology can serve as abridge between molecular events, interactions, and motor speeds and the overall transportat the cell scale. Biological insight into the effective transport and effective diffusion rates93is obtained through the QSS methodology. For example, in the kinesin-dynein model, theQSS PDE revealed that changes in any number of parameters could switch the distributionof motors within the cell from one cell end to the other (see Figure 2.9). The biologicalsignificance of this finding is that a mutation or other change in one of the biophysicalproperties of the molecular motors can drastically alter the overall, cell-scale, distributionof motors. Such a mutation would have a detrimental impact on the normal function of thetransport process.The three-case studies illustrated the generality of the QSS approximation method toa variety of nonlinear interactions; however, it remains an open question whether the QSSapproximation is possible for more realistic models with a greater number of motor states;or in models with higher degree nonlinearities. Moreover, the QSS approximation relies onthe assumption that there is a separation of time-scales between binding/unbinding andmotion across the cell.Recent work by Ciocanel et al. [12] has estimated biophysical motor parameters fromactive transport data and a PDE approximation method for a system of transport equationswith linear reactions. Future work should use data from specific experimental systemsto validate the QSS approximation. Although the biological systems considered here aretruly three-dimensional (in space), the 1D models considered here are a sufficiently goodapproximation for fungal hyphae (or neurons) which are long thin cells. Nonetheless, itwould be interesting to extend the QSS approximation to higher spatial dimensions withnonlinear reactions (the QSS approximation method has been studied in the context ofmodels in two spatial dimensions with linear reactions [6]). A final interesting directionis the open question of what happens in systems with multiple quasi-steady-states, as inthe myosin model in Chapter 2. Of particular interest is a system with three quasi-steady-states in a bistable arrangement (two stable, with an unstable state in between). In such amodel, the first questions are to determine which QSS better approximates solutions to thefull PDE system, especially if the system is attracted to different QSS in different spatialregions within the cell.The Interplay Between Cell Signalling and Cell MechanicsThe contribution of the work in Chapter 3 was the first exploration of the interplay betweenGTPase-induced cell contraction and tension-induced GTPase activity in single cells and in a1D and 2D tissue. The main biological result is that a wide variety of cell behaviour emergesfrom this signalling feedback: contracted cells, relaxed cells, oscillatory cells, synchronizedcontractile tissue, and waves of contraction in a large tissue. This behaviour is consistentwith a range of biologically relevant cell behaviours. Although many models have lookedat the implications of GTPase signalling for cell behaviour [33, 49] and many others have94sought to understand the implications of mechanochemical interactions for cell behaviours[17, 34, 60], the work here stands as one of the first steps to understanding the interplaybetween GTPase signalling and cell mechanics within a mathematical model.The successes of the GTPase-tension modelling work are (1) the emergence of a varietyof cell behaviour from the interplay of signalling and mechanics (2) the theoretical under-standing possible from the use of a stripped-down, conceptual model and (3) the possible,immediate, extensions to a variety of specific experimental systems such as Dropsophiladorsal closure or “zippering” in the neural/epidermal boundary in the Ciona intestinalisembryo. The simplicity of the two ODE model allowed for numerical bifurcation analysis tocharacterize the solutions to the single cell model, but this was limited to the single cell sys-tem and could be extended to the multi-cell system. Using the Cellular Potts model (CPM)revealed the implications of 2D neighbour interactions in the GTPase-tension model, butleaves the question of how to translate between a Monte-Carlo step and a unit of time t.As future work, I can suggest the following projects.• Extending the specific GTPase-tension system to specific, data-driven, systems (asmentioned above) to explore the interplay between signalling and mechanics in regu-lating morphogenetic processes.• Obtaining a continuum limit partial differential equation of the 1D multi-cellularsystem to understand the transition between synchronization and waves of contractionas a bifurcation in mechanical or signalling parameters (as in [17]).• Utilizing coarse-graining methods to derive partial differential equation models thatapproximate the behaviour of the 2D CPM. This could provide a link to the sub-cellular signalling and the population-level outcomes. For example, Alber et al. [1],Lushnikov et al. [45], Turner et al. [94] have used coarse-graining methods to derivePDE descriptions of cell behaviour from cellular Potts models of cell migration, celladhesion and chemotaxis. In these examples, numerical simulations of the PDE modelmatch the CPM simulations, and the coefficients of the PDE are derived from the CPMparameters.• Determining how cell polarity and directionality in 2D cell migration models are af-fected by mechanical signalling and GTPase signalling.Multi-scale Modelling in Cell BiologyUpon reflecting on the conclusions, significance, and contribution of the results in thesis,it is evident that multi-scale mathematical modelling is useful for explaining qualitativebehaviour and suggesting plausible hypotheses for experimentally observed phenomena. To95illustrate, consider the following examples. First, from the QSS methods in Chapter 2,the effective transport and effective diffusion rates were found to depend on the biophys-ical properties of the molecular motors. In an experiment where the overall, cell-scale,distribution of motors or cargo is observed the quantitative measures would describe theeffective transport and diffusion rates instead of specific molecular motor behaviours. Fromsuch data (and further experiments to alter or inhibit different players in the transportprocess) it would be possible to obtain a realistic biophysical understanding of the motorsand interactions that control transport. Second, a result from Chapter 3 is the emergenceof oscillations of cell size in a large epithelial tissue. Such oscillations are also observed inepithelial monolayers [103, 104]. Experimental manipulations to increase or decrease theadhesions of cells (through drugs that target integrins, focal adhesions, or other cell-celljunctions, for example) and quantitative analysis of the resulting behaviour could test thehypothesis that increasing adhesion leads to increased synchronization within the mono-layer.Nonetheless, conceptual multi-scale models such as those in this thesis can have limita-tions. A few limitations of the work here are the (1) lack of connections to specific experi-mental data, (2) the usage of deterministic models, (3) over-simplification of cell signalling,(4) assuming that cells are homogenous internally, and (5) the myriad of other biologicalfactors that I have ignored. When coupled to data, model parameters can be estimated,models can test hypotheses and suggest new experiments (as suggested above). Noise, whichis inherent to may biophysical systems, has largely been ignored in this thesis. Althoughthe reaction-advection-diffusion PDE systems studied in Chapter 2 can be understood asaverages or approximations of the noisy agent-based behaviour of many molecular motors,and the Cellular Potts model (CPM) revealed additional single cell behaviour (stochas-tic switching between steady-states) in Chapter 3, a more careful treatment of the noisybiological processes may be helpful in a more mechanistic modelling approach.Even with the simplifications, assumptions, and limitations of the specific models, it isapparent that a wide variety of organization can be understood from interactions at differentlevels. In Chapter 2, I illustrated how, from sub-cellular interactions, quasi-steady-statemethods can serve as a tool understand cell-scale distribution of molecular motors and theeffects that particular biophysical parameters have on the resulting motor distributions. InChapter 3, I connected signalling and mechanics at the cell scale to generate and understandcellular and multi-cellular behaviours.96Bibliography[1] M. Alber, N. Chen, P. M. Lushnikov, and S. A. Newman. Continuous macroscopiclimit of a discrete stochastic model for interaction of living cells. Physical ReviewLetters, 99(16), Oct 2007. doi:10.1103/physrevlett.99.168102. → page 95[2] P. W. Baas, J. S. Deitch, M. M. Black, and G. A. Banker. Polarity orientation ofmicrotubules in hippocampal neurons: uniformity in the axon and nonuniformity inthe dendrite. Proceedings of the National Academy of Sciences, 85(21):8335–8339,Nov 1988. doi:10.1073/pnas.85.21.8335. → page 18[3] C. Bakal, J. Aach, G. Church, and N. Perrimon. Quantitative morphologicalsignatures define local signaling networks regulating cell morphology. Science, 316(5832):1753–1756, Jun 2007. doi:10.1126/science.1140324. → pages 8, 9, 70[4] D. Bhat and M. Gopalakrishnan. Effectiveness of a dynein team in a tug of warhelped by reduced load sensitivity of detachment: evidence from the study ofbidirectional endosome transport ind. discoideum. Physical Biology, 9(4):046003,Jun 2012. doi:10.1088/1478-3975/9/4/046003. → pages 2, 18[5] T. L. Blasius, N. Reed, B. M. Slepchenko, and K. J. Verhey. Recycling of kinesin-1motors by diffusion after transport. PLoS ONE, 8(9):e76081, Sep 2013.doi:10.1371/journal.pone.0076081. → pages 2, 17[6] P. C. Bressloff and J. M. Newby. Quasi-steady-state analysis of two-dimensionalrandom intermittent search processes. Physical Review E, 83(6), Jun 2011.doi:10.1103/physreve.83.061139. → pages 67, 94[7] P. C. Bressloff and J. M. Newby. Stochastic models of intracellular transport.Reviews of Modern Physics, 85(1):135–196, Jan 2013.doi:10.1103/revmodphys.85.135. → pages 2, 18, 22, 27, 64[8] P. R. Burton. Dendrites of mitral cell neurons contain microtubules of oppositepolarity. Brain Research, 473(1):107–115, Nov 1988.doi:10.1016/0006-8993(88)90321-6. → page 18[9] K. M. Byrne, N. Monsefi, J. C. Dawson, A. Degasperi, J.-C. Bukowski-Wills,N. Volinsky, M. Dobrzyn´ski, M. R. Birtwistle, M. A. Tsyganov, A. Kiyatkin, andet al. Bistability in the Rac1, PAK, and RhoA signaling network drives actincytoskeleton dynamics and cell motility switches. Cell Systems, 2(1):38–48, Jan2016. doi:10.1016/j.cels.2016.01.003. → pages 9, 70, 71, 8497[10] D. Chowdhury, A. Schadschneider, and K. Nishinari. Traffic Phenomena in Biology:From Molecular Motors to Organisms. Traffic Granul. Flow, 2(4):223–238, 2005. →pages 2, 17[11] L. Ciandrini, M. C. Romano, and A. Parmeggiani. Stepping and crowding ofmolecular motors: Statistical kinetics from an exclusion process perspective.Biophysical Journal, 107(5):1176–1184, Sep 2014. doi:10.1016/j.bpj.2014.07.012. →pages 2, 18[12] M.-V. Ciocanel, J. A. Kreiling, J. A. Gagnon, K. L. Mowry, and B. Sandstede.Analysis of active transport by fluorescence recovery after photobleaching.Biophysical Journal, 112(8):1714–1725, Apr 2017. doi:10.1016/j.bpj.2017.02.042. →page 94[13] R. Clewley, W. Sherwood, M. LaMar, and J. Guckenheimer. PyDSTool, a softwareenvironment for dynamical systems modeling, 2007. → page 118[14] S. Cooper, A. Sadok, V. Bousgouni, and C. Bakal. Apolar and polar transitionsdrive the conversion between amoeboid and mesenchymal shapes in melanoma cells.Molecular Biology of the Cell, 26(22):4163–4170, Aug 2015.doi:10.1091/mbc.e15-06-0382. → pages 8, 9, 70[15] T. Das, K. Safferling, S. Rausch, N. Grabe, H. Boehm, and J. P. Spatz. A molecularmechanotransduction pathway regulates collective migration of epithelial cells.Nature Cell Biology, 17(3):276–287, Feb 2015. doi:10.1038/ncb3115. → pages 8, 69[16] D. Dauvergne and L. Edelstein-Keshet. Application of quasi-steady state methodsto molecular motor transport on microtubules in fungal hyphae. Journal ofTheoretical Biology, 379:47–58, Aug 2015. doi:10.1016/j.jtbi.2015.04.033. → pages2, 3, 18, 19, 22, 30, 35, 36, 37, 39, 64, 65, 93[17] K. Dierkes, A. Sumi, J. Solon, and G. Salbreux. Spontaneous oscillations of elasticcontractile materials with turnover. Physical Review Letters, 113(14):10, Oct 2014.doi:10.1103/physrevlett.113.148102. → pages 90, 95[18] D. E. Discher. Tissue cells feel and respond to the stiffness of their substrate.Science, 310(5751):1139–1143, Nov 2005. doi:10.1126/science.1116995. → pages 8, 68[19] R. Dixit, J. L. Ross, Y. E. Goldman, and E. L. F. Holzbaur. Differential regulationof dynein and kinesin motor proteins by tau. Science, 319(5866):1086–1089, Feb2008. doi:10.1126/science.1152993. → page 19[20] A. Diz-Mun˜oz, D. A. Fletcher, and O. D. Weiner. Use the force: membrane tensionas an organizer of cell shape and motility. Trends in Cell Biology, 23(2):47–53, Feb2013. doi:10.1016/j.tcb.2012.09.006. → page 69[21] A. Diz-Mun˜oz, K. Thurley, S. Chintamen, S. J. Altschuler, L. F. Wu, D. A. Fletcher,and O. D. Weiner. Membrane tension acts through PLD2 and mTORC2 to limitactin network assembly during neutrophil migration. PLOS Biology, 14(6):e1002474,Jun 2016. doi:10.1371/journal.pbio.1002474. → page 6998[22] L. Edelstein-Keshet. Mathematical Models in Biology. Society for Industrial andApplied Mathematics, Jan 2005. doi:10.1137/1.9780898719147. → page 5[23] G. Fink. Dynein-dependent motility of microtubules and nucleation sites supportspolarization of the tubulin array in the fungus Ustilago maydis. Molecular Biology ofthe Cell, 17(7):3242–3253, Apr 2006. doi:10.1091/mbc.e05-12-1118. → page 18[24] S. A. Freeman, S. Christian, P. Austin, I. Iu, M. L. Graves, L. Huang, S. Tang,D. Coombs, M. R. Gold, and C. D. Roskelley. Applied stretch initiates directionalinvasion through the action of Rap1 GTPase as a tension sensor. Journal of CellScience, 130(1):152–163, May 2017. doi:10.1242/jcs.180612. → page 69[25] M. Fukata, M. Nakagawa, and K. Kaibuchi. Roles of Rho-family GTPases in cellpolarisation and directional migration. Current Opinion in Cell Biology, 15(5):590–597, 2003. doi:10.1016/S0955-0674(03)00097-8. → pages 8, 68[26] J. Gou, L. Edelstein-Keshet, and J. Allard. Mathematical model with spatiallyuniform regulation explains long-range bidirectional transport of early endosomes infungal hyphae. Molecular Biology of the Cell, 25(16):2408–2415, Jun 2014.doi:10.1091/mbc.e14-03-0826. → pages 2, 18, 22, 39[27] P. Gross, K. V. Kumar, and S. W. Grill. How active mechanics and regulatorybiochemistry combine to form patterns in development. Annual Review ofBiophysics, 46(1):337–356, May 2017. doi:10.1146/annurev-biophys-070816-033602.→ page 89[28] C. Guilluy, R. Garcia-Mata, and K. Burridge. Rho protein crosstalk: another socialnetwork? Trends in Cell Biology, 21(12):718–726, Dec 2011.doi:10.1016/j.tcb.2011.08.002. → page 71[29] H. Hashimoto, F. B. Robin, K. M. Sherrard, and E. M. Munro. Sequentialcontraction and exchange of apical junctions drives zippering and neural tubeclosure in a simple chordate. Developmental Cell, 32(2):241–255, Jan 2015.doi:10.1016/j.devcel.2014.12.017. → pages 77, 91[30] L. He, J. Tao, D. Maity, F. Si, Y. Wu, T. Wu, V. Prasath, D. Wirtz, and S. X. Sun.Role of membrane-tension gated Ca2+ flux in cell mechanosensation. Journal of CellScience, 131(4):jcs208470, Jan 2018. doi:10.1242/jcs.208470. → pages 69, 73, 89, 90[31] A. G. Hendricks, E. Perlson, J. L. Ross, H. W. Schroeder, M. Tokito, and E. L.Holzbaur. Motor coordination via a tug-of-war mechanism drives bidirectionalvesicle transport. Current Biology, 20(8):697–702, Apr 2010.doi:10.1016/j.cub.2010.02.058. → pages 2, 18[32] J. H. R. Hetmanski, E. Zindy, J.-M. Schwartz, and P. T. Caswell. A MAPK-drivenfeedback loop suppresses Rac activity to promote RhoA-driven cancer cell invasion.PLOS Computational Biology, 12(5):e1004909, May 2016.doi:10.1371/journal.pcbi.1004909. → page 9099[33] W. R. Holmes and L. Edelstein-Keshet. Analysis of a minimal Rho-GTPase circuitregulating cell shape. Physical Biology, 13(4):046001, Jul 2016.doi:10.1088/1478-3975/13/4/046001. → pages 9, 11, 71, 72, 74, 84, 94, 120[34] W. R. Holmes, J. Park, A. Levchenko, and L. Edelstein-Keshet. A mathematicalmodel coupling polarity signaling to cell adhesion explains diverse cell migrationpatterns. PLOS Computational Biology, 13(5):e1005524, May 2017.doi:10.1371/journal.pcbi.1005524. → pages 9, 11, 69, 70, 71, 74, 84, 95, 120[35] A. R. Houk, A. Jilkine, C. O. Mejean, R. Boltyanskiy, E. R. Dufresne, S. B.Angenent, S. J. Altschuler, L. F. Wu, and O. D. Weiner. Membrane tensionmaintains cell polarity by confining signals to the leading edge during neutrophilmigration. Cell, 148(1-2):175–188, Jan 2012. doi:10.1016/j.cell.2011.10.050. → pages8, 69, 89[36] A. Jilkine, A. F. M. Mare´e, and L. Edelstein-Keshet. Mathematical model forspatial segregation of the Rho-family GTPases based on inhibitory crosstalk.Bulletin of Mathematical Biology, 69(6):1943–1978, 2007.doi:10.1007/s11538-007-9200-6. → page 72[37] A. Katsumi, J. Milanini, W. B. Kiosses, M. A. del Pozo, R. Kaunas, S. Chien, K. M.Hahn, and M. A. Schwartz. Effects of cell tension on the small GTPase Rac. TheJournal of Cell Biology, 158(1):153–164, Jul 2002. doi:10.1083/jcb.200201105. →pages 8, 69[38] R. Kaunas, P. Nguyen, S. Usami, and S. Chien. From the cover: Cooperative effectsof Rho and mechanical stretch on stress fiber organization. Proceedings of theNational Academy of Sciences, 102(44):15895–15900, Oct 2005.doi:10.1073/pnas.0506041102. → page 92[39] T.-H. Kim, N. Monsefi, J.-H. Song, A. von Kriegsheim, D. Vandamme, O. Pertz,B. N. Kholodenko, W. Kolch, and K.-H. Cho. Network-based identification offeedback modules that control RhoA activity and cell migration. Journal ofMolecular Cell Biology, 7(3):242–252, Mar 2015. doi:10.1093/jmcb/mjv017. → page90[40] S. Klumpp and R. Lipowsky. Cooperative cargo transport by several molecularmotors. Proceedings of the National Academy of Sciences, 102(48):17284–17289, Nov2005. doi:10.1073/pnas.0507363102. → pages 2, 18[41] A. Kumar, M. Ouyang, K. Van den Dries, E. J. McGhee, K. Tanaka, M. D.Anderson, A. Groisman, B. T. Goult, K. I. Anderson, and M. A. Schwartz. Talintension sensor reveals novel features of focal adhesion force transmission andmechanosensitivity. The Journal of Cell Biology, 213(3):371–383, May 2016.doi:10.1083/jcb.201510012. → page 69[42] Y. Kuramoto. Chemical Oscillations, Waves, and Turbulence, volume 19. SpringerBerlin Heidelberg, 1984. doi:10.1007/978-3-642-69689-3. → pages 82, 122100[43] B. Ladoux, R.-M. Me`ge, and X. Trepat. Front–rear polarization by mechanical cues:From single cells to tissues. Trends in Cell Biology, 26(6):420–433, Jun 2016.doi:10.1016/j.tcb.2016.02.002. → page 92[44] C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez, and J. Howard.Molecular crowding creates traffic jams of kinesin motors on microtubules.Proceedings of the National Academy of Sciences, 109(16):6100–6105, Mar 2012.doi:10.1073/pnas.1107281109. → page 19[45] P. M. Lushnikov, N. Chen, and M. Alber. Macroscopic dynamics of biological cellsinteracting via chemotaxis and direct contact. Physical Review E, 78(6), Dec 2008.doi:10.1103/physreve.78.061904. → page 95[46] M. Machacek, L. Hodgson, C. Welch, H. Elliott, O. Pertz, P. Nalbant, A. Abell,G. L. Johnson, K. M. Hahn, and G. Danuser. Coordination of Rho GTPaseactivities during cell protrusion. Nature, 461(7260):99–103, Aug 2009.doi:10.1038/nature08242. → page 71[47] R. Magno, V. A. Grieneisen, and A. F. Mare´e. The biophysical nature of cells:potential cell behaviours revealed by analytical and computational studies of cellsurface mechanics. BMC Biophysics, 8(1):8, May 2015.doi:10.1186/s13628-015-0022-x. → pages 79, 121[48] R. Mallik, A. K. Rai, P. Barak, A. Rai, and A. Kunwar. Teamwork in microtubulemotors. Trends in Cell Biology, 23(11):575–582, Nov 2013.doi:10.1016/j.tcb.2013.06.003. → pages 2, 18[49] A. F. M. Mare´e, A. Jilkine, A. Dawes, V. A. Grieneisen, and L. Edelstein-Keshet.Polarization and movement of keratocytes: A multiscale modelling approach.Bulletin of Mathematical Biology, 68(5):1169–1211, Jun 2006.doi:10.1007/s11538-006-9131-7. → pages 5, 71, 94[50] A. F. M. Mare´e, V. A. Grieneisen, and L. Edelstein-Keshet. How cells integratecomplex stimuli: The effect of feedback from phosphoinositides and cell shape oncell polarization and motility. PLoS Computational Biology, 8(3):e1002402, Mar2012. doi:10.1371/journal.pcbi.1002402. → page 90[51] P. K. Mattila and P. Lappalainen. Filopodia: molecular architecture and cellularfunctions. Nature Reviews Molecular Cell Biology, 9(6):446–454, May 2008.doi:10.1038/nrm2406. → page 26[52] J. McCormack, N. J. Welsh, and V. M. M. Braga. Cycling around cell-cell adhesionwith Rho GTPase regulators. Journal of Cell Science, 126(2):379–391, Jan 2013.doi:10.1242/jcs.097923. → page 92[53] D. P. McVicker, L. R. Chrin, and C. L. Berger. The nucleotide-binding state ofmicrotubules modulates kinesin processivity and the ability of tau to inhibitkinesin-mediated transport. Journal of Biological Chemistry, 286(50):42873–42880,Oct 2011. doi:10.1074/jbc.m111.292987. → page 19101[54] Y. Mori, A. Jilkine, and L. Edelstein-Keshet. Wave-pinning and cell polarity from abistable reaction-diffusion system. Biophysical Journal, 94(9):3684–3697, May 2008.doi:10.1529/biophysj.107.120824. → pages 8, 9, 11, 68, 72, 89[55] M. J. I. Mu¨ller, S. Klumpp, and R. Lipowsky. Motility states of molecular motorsengaged in a stochastic tug-of-war. Journal of Statistical Physics, 133(6):1059–1081,Dec 2008. doi:10.1007/s10955-008-9651-7. → pages 2, 18[56] R. Nambiar, R. E. McConnell, and M. J. Tyska. Myosin motor function: the ins andouts of actin-based membrane protrusions. Cellular and Molecular Life Sciences, 67(8):1239–1254, Jan 2010. doi:10.1007/s00018-009-0254-5. → page 26[57] J. Newby and P. C. Bressloff. Local synaptic signaling enhances the stochastictransport of motor-driven cargo in neurons. Physical Biology, 7(3):036004, Aug2010. doi:10.1088/1478-3975/7/3/036004. → pages 18, 19, 22, 30, 63, 64, 93[58] J. M. Newby and P. C. Bressloff. Quasi-steady state reduction of molecularmotor-based models of directed intermittent search. Bulletin of MathematicalBiology, 72(7):1840–1866, Feb 2010. doi:10.1007/s11538-010-9513-8. → pages 2, 3, 65[59] C. D. Nobes and A. Hall. Rho GTPases control polarity, protrusion, and adhesionduring cell movement. The Journal of Cell Biology, 144(6):1235–1244, Mar 1999.doi:10.1083/jcb.144.6.1235. → pages 8, 68[60] G. Odell, G. Oster, B. Burnside, and P. Alberch. A mechanical model for epithelialmorphogenesis. Journal of Mathematical Biology, 9(3):291–295, May 1980.doi:10.1007/bf00276030. → pages ix, 9, 10, 11, 12, 13, 14, 69, 89, 95[61] G. Odell, G. Oster, P. Alberch, and B. Burnside. The mechanical basis ofmorphogenesis. Developmental Biology, 85(2):446–462, Jul 1981.doi:10.1016/0012-1606(81)90276-1. → pages 9, 69, 77, 90[62] K. Ohashi, S. Fujiwara, and K. Mizuno. Roles of the cytoskeleton, cell adhesion andRho signalling in mechanosensing and mechanotransduction. Journal ofBiochemistry, 161(3):mvw082, Jan 2017. doi:10.1093/jb/mvw082. → pages8, 69, 72, 90[63] G. F. Oster, J. D. Murray, and A. K. Harris. Mechanical aspects of mesenchymalmorphogenesis. Development, 78:83–125, 1983. → pages 9, 69, 90[64] M. Otsuji, S. Ishihara, C. Co, K. Kaibuchi, A. Mochizuki, and S. Kuroda. A massconserved reaction–diffusion system captures properties of cell polarity. PLoSComputational Biology, 3(6):e108, 2007. doi:10.1371/journal.pcbi.0030108. → pages8, 68[65] J. Park, W. R. Holmes, S. H. Lee, H.-N. Kim, D.-H. Kim, M. K. Kwak, C. J. Wang,L. Edelstein-Keshet, and A. Levchenko. Mechanochemical feedback underliescoexistence of qualitatively distinct cell polarity patterns within diverse cellpopulations. Proceedings of the National Academy of Sciences, 114(28):E5750–E5759, Jun 2017. doi:10.1073/pnas.1700054114. → pages 9, 69, 70, 71102[66] A. Parmeggiani, T. Franosch, and E. Frey. Totally asymmetric simple exclusionprocess with langmuir kinetics. Physical Review E, 70(4), Oct 2004.doi:10.1103/physreve.70.046101. → pages 2, 18[67] M. Parri and P. Chiarugi. Rac and Rho GTPases in cancer cell motility control.Cell Communication and Signaling, 8(1):23, 2010. doi:10.1186/1478-811x-8-23. →pages 71, 84[68] O. Pertz. Designing biosensors for Rho family proteins – deciphering the dynamicsof rho family gtpase activation in living cells. Journal of Cell Science, 117(8):1313–1318, Mar 2004. doi:10.1242/jcs.01117. → page 69[69] M. Pines, R. Das, S. J. Ellis, A. Morin, S. Czerniecki, L. Yuan, M. Klose,D. Coombs, and G. Tanentzapf. Mechanical force regulates integrin turnover indrosophila in vivo. Nature Cell Biology, 14(9):935–943, Aug 2012.doi:10.1038/ncb2555. → page 69[70] N. A. Reed, D. Cai, T. L. Blasius, G. T. Jih, E. Meyhofer, J. Gaertig, and K. J.Verhey. Microtubule acetylation promotes kinesin-1 binding and transport. CurrentBiology, 16(21):2166–2172, Nov 2006. doi:10.1016/j.cub.2006.09.014. → page 19[71] M. Reffay, M. C. Parrini, O. Cochet-Escartin, B. Ladoux, A. Buguin, S. Coscoy,F. Amblard, J. Camonis, and P. Silberzan. Interplay of RhoA and mechanical forcesin collective cell migration driven by leader cells. Nature Cell Biology, 16(3):217–223, Feb 2014. doi:10.1038/ncb2917. → page 92[72] T. Reichenbach, E. Frey, and T. Franosch. Traffic jams induced by rare switchingevents in two-lane transport. New Journal of Physics, 9(6):159–159, Jun 2007.doi:10.1088/1367-2630/9/6/159. → pages 2, 18[73] A. J. Ridley. Rho family proteins: coordinating cell responses. Trends in CellBiology, 11(12):471–477, 2001. doi:10.1016/S0962-8924(01)02153-5. → pages 8, 68[74] A. J. Ridley. Rho GTPase signalling in cell migration. Current Opinion in CellBiology, 36:103–112, Oct 2015. doi:10.1016/j.ceb.2015.08.005. → pages 8, 68[75] K. Rottner, A. Hall, and J. Small. Interplay between Rac and Rho in the control ofsubstrate contact dynamics in the control of substrate contact dynamics. CurrentBiology, 9(12):640–S1, Jun 1999. doi:10.1016/s0960-9822(99)80286-3. → page 71[76] A. K. Rzadzinska, M. E. Schneider, C. Davies, G. P. Riordan, and B. Kachar. Anactin molecular treadmill and myosins maintain stereocilia functional architectureand self-renewal. The Journal of Cell Biology, 164(6):887–897, Mar 2004.doi:10.1083/jcb.200310055. → page 26[77] B. Sabass, M. L. Gardel, C. M. Waterman, and U. S. Schwarz. High resolutiontraction force microscopy based on experimental and computational advances.Biophysical Journal, 94(1):207–220, Jan 2008. doi:10.1529/biophysj.107.113670. →page 69103[78] H. Sailem, V. Bousgouni, S. Cooper, and C. Bakal. Cross-talk between Rho and RacGTPases drives deterministic exploration of cellular shape space and morphologicalheterogeneity. Open Biology, 4(1):130132–130132, Jan 2014.doi:10.1098/rsob.130132. → pages 8, 9, 70, 71, 84[79] M. E. Schneider, A. C. Dose, F. T. Salles, W. Chang, F. L. Erickson, B. Burnside,and B. Kachar. A new compartment at stereocilia tips defined by spatial andtemporal patterns of myosin IIIa expression. Journal of Neuroscience, 26(40):10243–10252, Oct 2006. doi:10.1523/jneurosci.2812-06.2006. → page 19[80] M. Schuster, S. Kilaru, G. Fink, J. Collemare, Y. Roger, and G. Steinberg.Kinesin-3 and dynein cooperate in long-range retrograde endosome motility along anonuniform microtubule array. Molecular Biology of the Cell, 22(19):3645–3657, Aug2011. doi:10.1091/mbc.e11-03-0217. → page 18[81] M. Schuster, R. Lipowsky, M.-A. Assmann, P. Lenz, and G. Steinberg. Transientbinding of dynein controls bidirectional long-range motility of early endosomes.Proceedings of the National Academy of Sciences, 108(9):3618–3623, Feb 2011.doi:10.1073/pnas.1015839108. → pages 18, 22, 24[82] M. Schwander, B. Kachar, and U. Mu¨ller. The cell biology of hearing. The Journalof Cell Biology, 190(1):9–20, Jul 2010. doi:10.1083/jcb.201001138. → page 26[83] L. A. Segel and M. Slemrod. The quasi-steady-state assumption: A case study inperturbation. SIAM Review, 31(3):446–477, Sep 1989. doi:10.1137/1031091. → page5[84] G. Shubeita and S. Gross. Intracellular transport: Relating single-moleculeproperties to in vivo function. Comprehensive Biophysics, 4:287–297, 2012.doi:10.1016/b978-0-12-374920-8.00418-5. → page 18[85] D. Smith and R. Simmons. Models of motor-assisted transport of intracellularparticles. Biophysical Journal, 80(1):45–68, Jan 2001.doi:10.1016/s0006-3495(01)75994-2. → pages 2, 18[86] G. Steinberg. Motors in fungal morphogenesis: cooperation versus competition.Current Opinion in Microbiology, 14(6):660–667, Dec 2011.doi:10.1016/j.mib.2011.09.013. → pages 2, 18, 22[87] G. Steinberg, R. Wedlich-Soldner, M. Brill, and I. Schulz. Microtubules in the fungalpathogen Ustilago maydis are highly dynamic and determine cell polarity. Journalof Cell Science, 114(3):609–622, 2001. → page 18[88] M. C. Stone, F. Roegiers, and M. M. Rolls. Microtubules have opposite orientationin axons and dendrites of drosophila neurons. Molecular Biology of the Cell, 19(10):4122–4129, Jul 2008. doi:10.1091/mbc.e07-10-1079. → page 18[89] S. H. Strogatz. From Kuramoto to Crawford: exploring the onset of synchronizationin populations of coupled oscillators. Physica D: Nonlinear Phenomena, 143(1-4):1–20, Sep 2000. doi:10.1016/s0167-2789(00)00094-4. → pages 82, 122104[90] M. H. Swat, G. L. Thomas, J. M. Belmonte, A. Shirinifard, D. Hmeljak, and J. A.Glazier. Multi-Scale Modeling of Tissues Using CompuCell3D, chapter Multi-ScaleModeling of Tissues Using CompuCell3D, pages 325–366. Elsevier, 2012.doi:10.1016/b978-0-12-388403-9.00013-8. → pages 15, 79, 81, 85, 118, 120[91] A. Tomar and D. D. Schlaepfer. Focal adhesion kinase: switching between GAPsand GEFs in the regulation of cell motility. Current Opinion in Cell Biology, 21(5):676–683, Oct 2009. doi:10.1016/j.ceb.2009.05.006. → page 69[92] K. Tsujita, T. Takenawa, and T. Itoh. Feedback regulation between plasmamembrane tension and membrane-bending proteins organizes cell polarity duringleading edge formation. Nature Cell Biology, 17(6):749–758, May 2015.doi:10.1038/ncb3162. → page 69[93] M. A. Tsyganov, W. Kolch, and B. N. Kholodenko. The topology design principlesthat determine the spatiotemporal dynamics of G-protein cascades. MolecularBioSystems, 8(3):730, 2012. doi:10.1039/c2mb05375f. → page 90[94] S. Turner, J. A. Sherratt, K. J. Painter, and N. J. Savill. From a discrete to acontinuous model of biological cell movement. Physical Review E, 69(2), Feb 2004.doi:10.1103/physreve.69.021910. → page 95[95] A. Van Keymeulen, K. Wong, Z. A. Knight, C. Govaerts, K. M. Hahn, K. M.Shokat, and H. R. Bourne. To stabilize neutrophil polarity, PIP3 and Cdc42augment RhoA activity at the back as well as signals at the front. The Journal ofCell Biology, 174(3):437–445, Jul 2006. doi:10.1083/jcb.200604113. → pages 8, 68[96] F. M. Vega, M. Thomas, N. Reymond, and A. J. Ridley. The Rho GTPase RhoBregulates cadherin expression and epithelial cell-cell interaction. Cell Communicationand Signaling, 13(1):6, 2015. doi:10.1186/s12964-015-0085-y. → page 92[97] V. Vogel and M. Sheetz. Local force and geometry sensing regulate cell functions.Nature Reviews Molecular Cell Biology, 7(4):265–275, Feb 2006.doi:10.1038/nrm1890. → pages 8, 68, 69[98] W. Wang, K. Tao, J. Wang, G. Yang, Q. Ouyang, Y. Wang, L. Zhang, and F. Liu.Exploring the inhibitory effect of membrane tension on cell polarization. PLOSComputational Biology, 13(1):e1005354, Jan 2017. doi:10.1371/journal.pcbi.1005354.→ page 89[99] K. Wong, O. Pertz, K. Hahn, and H. Bourne. Neutrophil polarization:Spatiotemporal dynamics of RhoA activity support a self-organizing mechanism.Proceedings of the National Academy of Sciences, 103(10):3639–3644, Feb 2006.doi:10.1073/pnas.0600092103. → pages 8, 68[100] S. Woolner, A. Jacinto, and P. Martin. The small GTPase Rac plays multiple rolesin epithelial sheet fusion—dynamic studies of Drosophila dorsal closure.Developmental Biology, 282(1):163–173, Jun 2005. doi:10.1016/j.ydbio.2005.03.005.→ page 91105[101] A. Yochelis, S. Ebrahim, B. Millis, R. Cui, B. Kachar, M. Naoz, and N. S. Gov.Self-organization of waves and pulse trains by molecular motors in cellularprotrusions. Scientific Reports, 5(1), Sep 2015. doi:10.1038/srep13521. → pages2, 24, 26, 27, 66[102] A. Yochelis, T. Bar-On, and N. S. Gov. Reaction–diffusion–advection approach tospatially localized treadmilling aggregates of molecular motors. Physica D:Nonlinear Phenomena, 318-319:84–90, Apr 2016. doi:10.1016/j.physd.2015.10.023. →pages 3, 24, 26, 27, 52, 66[103] S. M. Zehnder, M. Suaris, M. M. Bellaire, and T. E. Angelini. Cell volumefluctuations in MDCK monolayers. Biophysical Journal, 108(2):247–250, Jan 2015.doi:10.1016/j.bpj.2014.11.1856. → pages 91, 96[104] S. M. Zehnder, M. K. Wiatt, J. M. Uruena, A. C. Dunn, W. G. Sawyer, and T. E.Angelini. Multicellular density fluctuations in epithelial monolayers. PhysicalReview E, 92(3):1–8, Sep 2015. doi:10.1103/physreve.92.032729. → pages 91, 96[105] C. Zihni and S. J. Terry. RhoGTPase signalling at epithelial tight junctions:Bridging the GAP between polarity and cancer. The International Journal ofBiochemistry & Cell Biology, 64:120–125, Jul 2015. doi:10.1016/j.biocel.2015.02.020.→ page 92[106] C. Zmurchok, T. Small, M. J. Ward, and L. Edelstein-Keshet. Application ofquasi-steady-state methods to nonlinear models of intracellular transport bymolecular motors. Bulletin of Mathematical Biology, 79(9):1923–1978, Jul 2017.doi:10.1007/s11538-017-0314-1. → page v[107] C. Zmurchok, D. Bhaskar, and L. Edelstein-Keshet. Coupling mechanical tensionand GTPase signaling to generate cell and tissue dynamics. Physical Biology, Feb2018. doi:10.1088/1478-3975/aab1c0. → page v106Appendix ASupporting Materials for theApplication of Quasi-Steady-StateMethods to Nonlinear Models ofIntracellular TransportA.1 Convergence of the Kinesin Model to the QSS for ε→ 0As discussed in §2.6, one advantage of the QSS methodology is that it is not necessary to per-form the possibly numerically expensive task of computing time-dependent or steady-statesolutions to the full PDE system for each small value of the parameter ε. Computationalsavings would be amplified if the methodology was extended to 2 or 3 dimensions or themethodology was required in an experimental context or being fit to data. Since numericalcomputations for the full models and the QSS PDEs herein require approximately the sameamount of computational time, there is hardly a computational advantage to using theQSS PDE as a proxy for the full model. Nonetheless, repeated time-dependent numericalsimulations of the regularized myosin model for the creation of Figure 2.17 did require asignificant amount of computational time. The need for repeated numerical simulationsfor different initial conditions or for different parameters is suggestive of the computationaladvantages of the QSS approximation.Using the QSS as a proxy for the full PDE system does incur errors that depend on thesize of ε. To illustrate this, I compare the steady-state solution of the QSS PDE (2.29) and(2.30) with corresponding steady-state results computed from the full model (2.15) with(2.25) for a few values of ε in Figure A.1. These results show that as ε decreases the QSSPDE accurately predicts the steady-state of the full model.1070 0.5 1x0246QSSε = 0.1ε = 0.01ε = 0.001Figure A.1: The steady-state solution to the full nonlinear model (2.15) with (2.25)converges to the steady-state of the QSS PDE (2.29) with (2.30) as ε → 0. Hereg(α) = α/(1 + α), P = 0.6, kb = 0.5, ku = 0.3, and D = 0.01.A.2 Numerical Methods for the QSSIn this appendix, I show how to numerically compute the steady-state solution of the QSSPDEs by recasting the nonlocal problem into an initial-boundary value problem (IBVP),which is amenable to a numerical shooting method.For the QSS PDE associated with the kinesin model (2.29) of §2.4.1, the steady-stateproblem isdαdx=kaD[2P (x)− 1] g(α),∫ 10(kag(α) + α) dx = 1, (A.1)where g(α) is either the saturated binding model (2.34) or the Hill function (2.47). Toreformulate (A.1), define N(x) byN(x) ≡∫ x0(kag[α(η)] + α(η)) dη − 1. (A.2)Then, (A.1) is equivalent to the ODE systemdαdx=kaD[2P (x)− 1] g(α), dNdx= kag(α) + α, (A.3)with N(0) = −1. In order to find a solution α(x) that satisfies (A.1), it is necessary tofind α(0) such that when the initial value problem (IVP) for α(x) is solved, N(1) = 0.Specify α(0) = β, where β is a value to be determined. Next, solve the IBVPs (A.3) forvarious values of β and output the quantity N(1;β). In this numerical shooting procedure,Newton’s method on β is then used to satisfy the required terminal constraint N(1;β) = 0.Once this β is found, the steady-state solution α(x) can be calculated by solving the IVPwith the initial condition α(0) = β.108A similar approach can be used to compute steady-state solutions of the QSS PDE (2.53)for the kinesin-dynein model of §2.4.2 subject to the total mass constraint ∫ 10 y(x) dx = 1.In place of (A.3), obtaindαdx= −kaD[v(kα+ 1−Q)−Q](kα+ 1−Q)2 +Q(1−Q) (kα+ 1−Q)α,dNdx=(1 +1ka)(kα+ 1)αkα+ 1−Q,(A.4)with N(0) = −1 and α(0) = β. As before, β is a shooting parameter determined numericallyby satisfying the terminal constraint N(1;β) = 0.Finally, we consider steady-state solutions of the QSS PDE (2.67a) for the myosinmodel of §2.4.3 subject to the total mass constraint ∫ 10 y(x) dx = 1. In place of (A.3),obtaindαdx= −kbD(vkbwα2 − 1)kbwα2 − 1 α,dNdx=(kb + 1)kbkbw(kbwα2 + 1)α, (A.5)with N(0) = −1 and α(0) = β, again where β is computed numerically to satisfy theconstraint N(1;β) = 0. A steady-state solution exists only when kbwα2 > 1 on 0 ≤ x ≤ 1.To numerically determine the boundary in parameter space where kbwα2 > 1 holds on0 ≤ x ≤ 1 for the steady-state when 0 < v < 1, it is convenient to reformulate (A.5). DefineA(x) ≡ √kbwα(x) to transform (A.5) todAdx= −c1(vA2 − 1)A2 − 1 A ,dNdx= c2(A2 + 1)A, where c1 ≡ kbD, c2 ≡ kb + 1kb√kbw.(A.6)A steady-state solution to the QSS PDE exists only when A(x) > 1 on 0 ≤ x ≤ 1. Since(A.6) implies that A(x) is monotonic in x whenever A > 1, then it is possible that A→ 1+only for x→ 0+ or x→ 1−. However, since A→ 1/√v > 1 on the infinite line as x→∞, itfollows that A→ 1+ as x→ 0+. To determine the local behaviour as A→ 1+ and x→ 0+,note that (A.6) implies dA/dx ∼ c1(1− v)/[2(A− 1)] and dN/dx ∼ 2c2. This yields thelocal behaviourA ∼ 1 +√c1(1− v)x, N ∼ −1 + 2c2x, as x→ 0+. (A.7)For a fixed v and D > 0, with 0 < v < 1, the region in the parameter space kbw versuskb where A(x) > 1 on 0 ≤ x ≤ 1, is determined as follows. Fix c1 in (A.6), numericallyintegrate the IBVPs (A.6) with the local behavior (A.7) imposed at some x = δ, with0 < δ 1, and numerically shoot on the value of c2 for which N(1; c2) = 0. From (A.6),this determines kb and kbw as kb = c1D and kbw = [(kb + 1)/(kbc2)]2.109A.3 Microtubule Density and Binding by Motor ComplexesA.3.1 Kinesin Model with Nonuniform MT DensityTo explicitly incorporate the possibility that MT density, m(x) (as well as fraction of MTpointing to the right, P (x)) varies across the cell, it is possible to write the kinesin-modelequations as∂pR∂t= −v∂pR∂x+ P (x)kbmm(x)g(pU)− kupR, (A.8a)∂pL∂t= v∂pL∂x+ (1− P (x))kbmm(x)g(pU)− kupL, (A.8b)∂pU∂t= D0∂2pU∂x2− kbmm(x)g(pU) + kupR + kupL. (A.8c)This modification of the model introduces another factor into coefficients that are alreadyspatially-dependent, but otherwise leaves the model structure unchanged. Hence, the tech-niques in the paper apply as before with kbmm(x) replacing the parameter kb.For the purposes of the proof-of-concept QSS reduction, now restrict attention to uni-form MT density so that m(x) ≡ m0 is a constant. Then the model for kinesin is given by(A.9) as below, with the assignmentkb = kbmm0.That is, the binding constant kb is understood to represent the net rate of binding, whichincludes both the per-MT-binding rate and the MT density.A.3.2 Kinesin-Dynein Model and the Function Q(x)The kinesin-dynein model simplifies the binding of free motor complexes into states thatmove right with probability Q(x), and left with probability 1−Q(x). I consider the case ofmotor complexes that all have nk kinesin and nd dynein components (the case of complexeswith a variety of motor numbers can be handled by considering the mean composition ofa complex or the mean ratio between the two motor types). Also define the parameterskbd and kbk as the binding rates for a (single) dynein and for a (single) kinesin to a MT,and consider m(x) as the local MT density. Then the quantity kbQ in the model can bedecomposed as follows:kbQ(x) = m(x) [P (x)nkkbk + (1− P (x))ndkbd] .110This relates the aggregate binding rate to the probability that a kinesin binds to right-pointing MT and that dynein binds to left-pointing MT. Similarly,kb(1−Q(x)) = m(x) [(1− P (x))nkkbk + P (x)ndkbd] .Since such details merely substitute one spatially-dependent function for another, the QSSanalysis previously described carries over as before.A.4 Scaling the QSS ModelsIn this section, details of the scaling of the molecular motor transport models is presented.A.4.1 The Kinesin ModelConsider the kinesin model with uniform MT density. This system is∂pR∂t= −v∂pR∂x+ Pkbg(pU)− kupR, (A.9a)∂pL∂t= v∂pL∂x+ (1− P )kbg(pU)− kupL, (A.9b)∂pU∂t= D0∂2pU∂x2− kbg(pU) + kupR + kupL. (A.9c)Define T to the total amount of motors inside the cell:T ≡∫ L00(pR(x) + pL(x) + pU(x))dx ≡∫ L00y(x) dx,and ρ = T/L0 to be the average density of motors in the cell.Scale space, time, and densities as follows:x? =xL0, t? =tvL0, pJ? =pJρ, y? =yρwhere y? = pR?+pL?+pU?is the total scaled density. Here, distance has been scaled by thecell length and time by the time that a motor takes to walk across the cell. The densitiesof motors in each state is scaled by the average motor density across the cell.With this scaling, the total amount of motors isT =∫ 10(ρpR?(x?) + ρpL?(x?) + ρpU?(x?))d(L0x?).111Taking out the constant factor of ρL0 ≡ T from the integral results inT = ρL0∫ 10(pR?(x?) + pL?(x?) + pU?(x?))dx?,which leads to ∫ 10y?dx? =∫ 10(pR?(x?) + pL?(x?) + pU?(x?))dx? = 1.With this scaling, the integral of the total scaled density is unity, which is assumed through-out the numerical computations above.Substituting the scaled variables into the PDE system (A.9) leads tovL0∂(ρpR?)∂t?=−vL0∂(ρpR?)∂x?+ ku(P (x)kbkug(ρpU?)− (ρpR?)), (A.10a)vL0∂(ρpL?)∂t?=vL0∂(ρpL?)∂x?+ ku((1− P (x))kbkug(ρpU?)− (ρpL?)), (A.10b)vL0∂(ρpU?)∂t?=D0L20∂2(ρpU?)∂x?2+ ku(ρpR?+ ρpL? − kbkug(ρpU?)). (A.10c)This leads to two cases, depending on whether the function g is linear or not.Case I: g is linear. In this case, it is possible to eliminate the factor ρ from everyterm. Dividing each term in the equations by vρ/L0 and dropping the stars leads to∂pR∂t= −∂pR∂x+1ε(P (x)kapU − pR) , (A.11a)∂pL∂t=∂pL∂x+1ε((1− P (x))kapU − pL), (A.11b)∂pU∂t= D∂2pU∂x2+1ε(pR + pL − kapU), (A.11c)where D, ε, and ka are defined byD ≡ D0vL0, ε ≡ vL0ku, ka ≡ kbku. (A.12)In this case, these dimensionless parameters represent, respectively, the ratio of (time to betransported : time to diffuse) across the cell (D), the ratio of (time spent unbound : timeto walk) across the cell (ε), and the ratio of (time spent unbound : time spent bound) (ka).Case II: g is Michaelian or Hill. In this case,g(p) = gmpnKn + pn, n = 1, 2, . . . .112Then, (A.10) becomesvL0∂(ρpR?)∂t?=−vL0∂(ρpR?)∂x?+ ku(P (x)kbkugm(ρpU?)n[Kn + (ρpU?)n]− (ρpR?)), (A.13a)vL0∂(ρpL?)∂t?=vL0∂(ρpL?)∂x?+ ku((1− P (x))kbkugm(ρpU?)n[Kn + (ρpU?)n]− (ρpL?)), (A.13b)vL0∂(ρpU?)∂t?=D0L20∂2(ρpU?)∂x?2+ ku(ρpR?+ ρpL? − kbkugm(ρpU?)n[Kn + (ρpU?)n]). (A.13c)Define a new constant A ≡ K/ρ. This constant is the ratio of the motor concentration atwhich the binding rate is half-maximal to the average motor density in the cell. Dividenumerator and denominator of the Hill function by ρn. Further, divide every term in theequations by vρ/L0 as before. After rearranging and dropping the starred notation, obtain∂pR∂t= −∂pR∂x+1ε(P (x)ka(pU)n[An + (pU)n]− pR), (A.14a)∂pL∂t=∂pL∂x+1ε((1− P (x))ka (pU)n[An + (pU)n]− pL), (A.14b)∂pU∂t= D∂2pU∂x2+1ε(pR + pL − ka (pU)n[An + (pU)n]), (A.14c)where D and ε are as before, but ka now depends on whether g is a Michaelis-Menten ora Hill function. This holds for any Hill coefficient n. Note that, in particular, for the casen = 1, which is the Michaelian case considered,∂pR∂t= −∂pR∂x+1ε(P (x)kapU[1 + cpU]− pR), (A.15a)∂pL∂t=∂pL∂x+1ε((1− P (x))ka pU[1 + cpU]− pL), (A.15b)∂pU∂t= D∂2pU∂x2+1ε(pR + pL − ka pU[1 + cpU]), (A.15c)where c ≡ 1/A = ρ/K. In (A.14) and (A.15) ka is defined byka ≡ kbgmkuρ(Hill), ka ≡ kbgmkuK, (Michaelis-Menten). (A.16)In either case, the parameter ka describes the ratio of time spent bound to the time spentunbound, mediated by the nonlinear binding kinetics.113Finally, scale the boundary conditions in (2.5) to get(pR − pL −D∂pU∂x)∣∣∣∣x=0,1= 0, (A.17)together withpR(0, t) = 0 and pL(1, t) = 0. (A.18)A.4.2 Kinesin-Dynein Model ScalingDefine kc ≡ krl − klr. Then the model can be written as∂pR∂t= −vr∂pR∂x+ kbQpU − kupR − kcpRpL, (A.19a)∂pL∂t= vl∂pL∂x+ kb(1−Q)pU − kupL + kcpRpL, (A.19b)∂pU∂t= D0∂2pU∂x2− kbpU + ku(pR + pL). (A.19c)Scale all variables as before. Then terms of the form (kc/ku)pRpL will lead to the form(kc/ku)ρpR?ρpL?, so that what remains, after canceling out a factor of vrρ/L0 from everyterm in each equation, and dropping the starred quantities, is∂pR∂t= −∂pR∂x+1ε(kaQpU − pR − kpRpL) , (A.20a)∂pL∂t= v∂pL∂x+1ε(ka(1−Q)pU − pL + kpRpL), (A.20b)∂pU∂t= D∂2pU∂x2+1ε(pR + pL − kapU), (A.20c)where the parameters arev ≡ vlvr, D ≡ D0vrL0, ε ≡ vrkuL0, ka ≡ kbku, k ≡ kcρku=(krl − klr)ρku. (A.21)Here ρ is the average density of motors inside the cell. These dimensionless parametersrepresent, respectively, the (left:right) walking speed ratio (v), the ratio of (time to betransported : time to diffuse) across the cell (D), the ratio of (time spent unbound : timeto walk) across the cell (ε), the ratio of (time spent unbound : time spent bound) (ka),and the turning parameter k, which represents the ratio of (net right-left direction switches: unbinding rate). Note that the average density of motors ρ enters into the turning rateparameter due to the nonlinearity of the model with respect to the turning of motors when114they collide on a MT.Details of QSS Reduction of Kinesin-Dynein ModelThis section provides some details of the QSS reduction of the kinesin-dynein model. Uponsetting f2 = f3 = 0 in (2.48), obtain the two equationskpRpL = pL − ka(1−Q)pU, −kapU + pR + pL = 0. (A.22)It is convenient to let pL be the free variable and parameterize the quasi-steady-state interms of pL = α. By solving (A.22) for pR and pU, the quasi-steady-state solution p0 asgiven in (2.50). The non-zero eigenvalues λ± of the Jacobian of the kinetics satisfy thequadratic equation given in (2.51) and (2.52). A necessary and sufficient condition forRe(λ±) < 0 is that σ1 < 0 and σ2 > 0 in (2.52). To establish this result, consider someproperties of H(Q) defined in (2.52). First observe that H(0) = 1, so that trivially σ1 < 0and σ2 > 0 when Q = 0. Then, since H′(Q) = −(1 + kα)/(1 + kα−Q)2 < 0, it followsthat σ1 < 0 and σ2 > 0 on 0 ≤ Q ≤ 1 provided that σ1 < 0 and σ2 > 0 when Q = 1.These inequalities do hold at Q = 1, since by using H(1) = (kα− 1)/(kα), one finds thatσ1 = −1 − ka − kα and σ2 = kα(1 + ka) > 0 when Q = 1. This proves that Re(λ±) < 0for any Q in 0 ≤ Q ≤ 1. As a result, p0 defined in (2.50) is a slow manifold in the senseof Definition (2.3.1) for any Q in 0 ≤ Q ≤ 1. Finally, by using p0 and the operator M , asdefined in (2.48), in the solvability condition (2.24), the QSS PDE can be derived (2.53).A.4.3 Myosin Model ScalingThe myosin model is∂pW∂t= −vw∂pW∂x− kˆbw(pB)2pW + kˆbpU − kupW, (A.23a)∂pB∂t= vb∂pB∂x+ kˆbw(pB)2pW − kupB, (A.23b)∂pU∂t= Df∂2pU∂x2− kˆbpU + ku(pB + pW). (A.23c)Using the scaling as before, the terms(pB)2pW will lead to the forms(ρpB?)2(ρpW?). Thiswill result in a constant factor ρ2 that remains after canceling out ρ from all terms in the115equation. As a result, upon dropping the starred quantities,∂pW∂t= −∂pW∂x+1ε(−kbw(pB)2pW + kbpU − pW), (A.24a)∂pB∂t= v∂pB∂x+1ε(kbw(pB)2pW − pB), (A.24b)∂pU∂t= D∂2pU∂x2+1ε(pB + pW − kbpU), (A.24c)where the dimensionless parameters v, D, ε, kbw, and kb are defined byv ≡ vbvw, D ≡ DfvwL0, ε ≡ vwkuL0, kbw ≡ kˆbwρ2ku, kb ≡ kˆbku. (A.25)Recall that ρ is the average density of motors inside the cell. These dimensionless parametersrepresent, respectively, the bound:walking motor speed ratio (v), the ratio of (time to betransported : time to diffuse) across the cell (D), the ratio of (time spent unbound : time towalk) across the cell (ε), the interaction parameter kbw, which represents the ratio of (netrate of collisions that result in direction change : unbinding rate), and the ratio of (timespent unbound : time spent bound) (kb). Note that the average density of motors ρ entersinto the interaction rate parameter due to the nonlinearity of the model with motor-motorinteraction.A.4.4 Non-spatial Myosin ModelIn §2.4.3, I seek to determine whether the Type I or Type II QSS PDE better approximatesthe behaviour of the full myosin system. To understand the behaviour, I study the non-spatial myosin model kinetics through a phase-plane analysis, where the advection anddiffusive processes in (A.24) are neglected.The non-spatial myosin model kinetics are described by the following system of ODEs:dpWdt= −kbw(pB)2pW+kbpU−pW , dpBdt= kbw(pB)2pW−pB , dpUdt= pB+pW−kbpU ,(A.26)where time has been scaled to remove the ε-dependence. Due to conservation of mass, it ispossible to write pU = 1 − pW − pB. This facilitates the reduction of this system of threeequations to a system of two equations:dpWdt= −kbw(pB)2pW + kb(1− pW − pB)− pW, (A.27a)dpBdt= kbw(pB)2pW − pB. (A.27b)116With kbw = 25 and kb = 3, a phase-plane analysis (see Figure A.2) reveals the existence of anunstable manifold which divides the (pW, pB) plane into two regions. For initial conditionsbelow this unstable manifold, the system converges to a steady-state with pB = 0, butpW > 0, as in Type I QSS. For initial conditions above this unstable manifold, the systemconverges to a steady-state with pB > 0, as in Type II QSS.Figure A.2: A phase-plane analysis of the non-spatial myosin model (A.27) reveals theexistence of an unstable manifold that divides (pW, pB) space into two regions. For initialconditions below the unstable manifold, the system tends to a steady-state with pB = 0,but for initial conditions above the unstable manifold, the system tends to a steady-statewith pB > 0.117Appendix BSupporting Materials for CouplingMechanical Tension and GTPaseSignalling to Generate Cell andTissue DynamicsB.1 Numerical MethodsNumerical integration of the single cell models and bifurcation analysis was preformedusing PyDSTool [13]. Numerical integration of the single-cell and multicellular models waspreformed using MATLAB 2017a (The MathWorks, Inc. Natick, Massachusetts, UnitedStates). Cellular Potts model simulations were produced with CompuCell3D [90].B.2 Scaling the GTPase ModelThe dynamics of active GTPase are governed by the following differential equation:dGdt=(bˆ+ γˆGnGn0 +Gn)(GT −G)− δG. (B.1)Here, bˆ is a basal activation rate, γˆ gives the magnitude of the positive feedback upon theactivation rate, and G0 describes the concentration of GTPase at which positive feedbackreaches its half-maximal effect. GT −G gives the total concentration of inactive GTPase.To reduce the size of parameter space, scale GTPase concentration by the half-maxquantity G0, and scale time by the active GTPase residence time 1/δ, respectively. The118equations becomedGdt=(b+ γGn1 +Gn)(GT −G)−G. (B.2)The mechanical stimulus term f(T ) was added to the activation rate, i.e., I assumed that itoperates via a GEF. I considered several forms of mechanical feedback from cell deformationto GTPase activity. Based on the idea that the difference between the current cell “length”L and the current cell “rest-length” L0 creates the tension that stimulates mechanosensitivepathways, express the feedback in terms of L and L0. Consequently, I experimented witheach of the following forms:f0(L) = β(L− L0), Linear case;f1(L) = βLmLm0 + Lm, Hill function;f2(L) = β11 + exp[α(L− L0)] , Squashing function;f3(L) = β11 + exp[α (L−L0)L0] , Strain-dependent squashing function.All four cases share the property that GTPase activation is amplified if L L0. Thelinear function f0 has the property that both stretching (L > L0) and compression (L < L0)affect GTPase activation, albeit in opposite ways (stretching increases while compressiondecreases the GTPase activation rate.) The squashing function f2 is predominately unidi-rectional, i.e., only L > L0 has a significant effect, so stretching, but not compressing a cellaffects its signalling. This function was used in the minimal model and has the advantageof specifically tracking tension. At the same time, f1 has a similar effect as f2, and wasto a large extent indistinguishable in the dynamical results obtained (see Appendix FigureB.1(b) for a bifurcation diagram of the single-cell model with the Hill function responsef1). The noticeable difference occurs in the synchronization of large tissue simulations.Compare, for example, the simulation with the squashing function f2 in Figure 3.5(e) andthe simulation with the strain-dependent squashing function f3 in Appendix Figure B.1(c).The Hill function and strain-dependent squashing function, f1 and f3 have a similarshape for all L and L0 and are approximately equal for the parameters used herein. Thechange that f1 and f3 can affect in the GTPase activation rate is relative to the currentrest-length of the cell L0. This is different from the squashing function f2, which assumes amechanosensing mechanism by which tension can activate GTPase signalling regardless ofthe current rest length, L0.The Rac-Rho mutual inhibition model in Section 3.5 Equation (3.6), is the scaled model.Details of this scaling are similar to scaling for the single GTPase model. The reader is119referred to [33, 34].B.3 1D Methods: Multicellular SimulationsEquations (3.4) and (3.2a) were implemented in a collection of cells in 1D. Each cell hasits own GTPase activity Gj , which is described by Equation (3.2a), with lengths given byLj = xj+1 − xj . Numerical integration was done using MATLAB 2017a (The MathWorks,Inc. Natick, Massachusetts, United States) for all 1D multicellular simulations.For Figure 3.5(a)-(c), (e), and (f), and B.1(a) and (c), GTPase activity in each cell, Gj ,affects the rest length throughLj,0 = `0 − φGpjGph +Gpj. (B.3)Tension is assumed to affect the GTPase activation rate through the squashing functionresponse to tension (f2 above, also Equation 3.2b). Parameter values for these simulationsare b = 0.1, γ = 1.5, n = p = 4, GT = 2, α = 10, `0 = 1, φ = 0.75, Gh = 0.3, k = 1, λ = 10,and β varies. Initial conditions are Lj(0) = 0.7 and Gj(0) = 1 for the N = 10 simulations,and random initial lengths with Gj(0) = 1 for all the simulations with N > 10 and forthe N = 10 case with one oscillatory cell, Figure B.1. Instead of the squashing functionresponse to tension (f2), the strain-dependent squashing function (f3) was used in FigureB.1(c).For Figure 3.5(d), we simulated N = 14 cells, and assumed linear responses for bothGTPase-activation from tension (f0) and for rest-length from GTPase activity:Lj,0 = `0 − φGj . (B.4)Initial conditions for this simulation were Lj(0) = 0.68, Gj(0) = 0.45 for all j with theexception of Gj(0) = 1.2 for j = 13, 14. Other parameters were b = 0.3, β = 0.35, φ = 0.7,n = 4, GT = 1.75, `0 = 1, k = 1, and λ = 10.B.4 2D Methods: Cellular Potts ModelCompuCell3D, an open-source implementation of the cellular Potts model, is used for 2Dsimulations of the GTPase-tension model [90]. The cellular Potts model is an individualcell-based model where each cell occupies one or more discrete lattice sites. Cells can expandoutwards or contract inwards by adding or removing lattice sites at the cell perimeter.The dynamics of each “cell” is governed by a Hamiltonian energy function, H. TheHamiltonian for the 2D simulation consists of an area constraint term (also called volume120(a) All 10 cells have β = 0.3except for one randomly chosenoscillatory cell with β = 0.2.0.0 0.1 0.2 0.3 0.4β0.20.40.60.8L(b) Single cell bifurcation di-agram with Hill function re-sponse from tension f1(T )(c) β = 0.2, 50 cells, withGTPase activation rate f(T ) =f3(T ).Figure B.1: Additional 1D tissue dynamics result from mechanochemical interactions. Ky-mographs show the 1D position of each cell (vertical axis) with color indicating the GTPaseactivity within each cell. In (a), a single oscillatory cell with β = 0.2 can induce tissue-leveloscillations among a population of contracted cells with β = 0.3. In (b), single cell dynamicswith the Hill function response from tension, f1(T ), qualitatively resemble the dynamicswith the squashing function f2(T ). In (c), waves of contraction propagate through the tissueof 50 oscillatory cells with the strain-dependent feedback, f3. See also SI Movies 10 and 11,for (a) and (c) respectively.deformation term in a general 3D context) and an adhesion energy term. The area constraintis implemented in terms of a (time-varying) target area. Target area represents the area(number of lattice sites) that each cell would occupy in an optimal lattice configuration.The adhesion energies specify the interactions between different cells and the surroundingmedium (extra-cellular space, or “medium”). Additionally, a connectivity constraint isimposed that penalizes the Hamiltonian if lattice sites for each cell do not form a connecteddomain. This avoids fragmentation of the “cells”.Lattice sites are added or removed from cells in “spin-copy attempts”. A spin-copyattempt is accepted if it decreases the overall energy of the system, as defined by theHamiltonian. A spin-copy attempt is also accepted with a non-zero probability if it resultsin a small increase in the Hamiltonian. The temperature parameter in the Boltzmann dis-tribution of accepted unfavourable spin-copies controls the degree of exploration of energet-ically unfavourable lattice configurations. Given N lattice sites, a collection of N spin-copyattempts constitutes one Monte-Carlo step (MCS) of the simulation. The Metropolis algo-rithm is used to determine the quasi-deterministic kinetics of lattice configurations evolvingunder the Hamiltonian. While CPM does not explicitly track forces, it has recently beenshown to correspond to other vertex-based simulations where forces are made explicit [47].In the case of single cells, the model parameters are τ = 2000, b = 0.1, m = 10, γ = 1.5,n = p = 4, GT = 2, ε = 0.1, a0 = 400, φ = 0.75, and Gh = 0.3. Single cell simulationsran for 250 MCS, with temperature parameter 30. The cell-cell and cell-medium adhesion121energies are set to 0.1 and 80 in the Hamiltonian H, respectively, and I did not impose aperimeter (surface) constraint. The area constraint parameter in the Hamiltonian H wasset to: λA = 1 and initial conditions were set as G(0) = 1, AT (0) = 320, and A(0) = 320.B.5 2D Methods: Patch Size and SynchronizationTo explore the idea that adhesion strength could affect the extent of synchrony amongthe cells in the tissue, I varied the adhesion energy in the Hamiltonian H among a smalltissue of 9 cells, with each cell in the oscillatory regime. Treating the tissue as a system ofcoupled oscillators and it is possible to numerically quantify the level of synchrony usingthe Kuramoto order parameter and the variance in the distribution of phase angles ofthe oscillators. The Kuramoto order parameter describes the degree of synchronizationin a collection of coupled oscillators (see [42], or [89] for a review). As adhesion-strengthincreases, the oscillators are more synchronized, with an apparent increase in the Kuramotoorder parameter, and a decrease in variance in the distribution of phase angles. See FigureB.3 and B.5-B.7.To determine the Kuramoto order parameter for the small tissue, the dominant fre-quency of each oscillating cell (i.e., frequency with highest magnitude) is determined overtime using a sliding window with the real-valued Fourier transform (RFFT). The fixedsize window contains time series data roughly equivalent to 3 periods of oscillation. TheKuramoto order parameter and variance in phase is calculated by determining the phasecorresponding to the dominant frequency for all nine oscillators.Model parameters are as in Section B.4, except the initial conditions for cell areas arerandomly chosen, initial conditions for GTPase are randomly chosen between 0 and 1 andinitial target area is also randomly chosen between 350 and 450.B.6 2D Methods: Large-Tissue SimulationsA circular tissue consisting of 373 cells with randomly chosen area (and β = 0.2, corre-sponding to the oscillatory regime) is used as the initial lattice configuration. Initial targetarea was set to the initial cell area for each cell. The simulation was carried out for 2000MCS. Initial target area is set equal to the initial area. Initial GTPase concentration israndomly chosen between 0 and 1. The remaining model parameters are as in Section B.4.B.7 Additional 2D ResultsIn this section, some additional 2D CPM simulation results are presented. Model parametersare as before, outlined in Section B.4. Each figure shows 8 snapshots of the cell behaviour,with the colour indicating the GTPase activity. Also shown are the cell area, target area,122Figure B.2: CPM initial lattice configuration for 9 oscillating cells with no mechanicalcoupling or adhesion and randomly chosen initial cell area and GTPase activity.and the GTPase activity over time. These results include:1. Figure B.8: a single relaxed cell with large constant area and low GTPase activitywith β = 0.05.2. Figure B.9: damped oscillations occur for β = 0.1.3. Figure B.10: a small amplitude limit cycle with β = 0.15.4. Figure B.11: stochastic switching between a low amplitude limit cycle and high am-plitude limit cycle with β = 0.175.Also, a large tissue simulation similar to Figure 3.7 and 3.8 but with weak adhesion isshown in Figure B.12.123Figure B.3: Time series of GTPase activity, Kuramoto order parameter and variance inphase for 9 independent oscillators shown Figure B.2. Note that initial cell area is equal toinitial cell target area, hence initial conditions are not fully randomized. Variance increaseswith time due to the stochastic nature of spin copy attempts that offset the initial conditions.124Figure B.4: CPM initial lattice configuration for 9 oscillating mechanically coupled cellswith randomly chosen initial cell area and GTPase activity. Each lattice site can onlybe occupied by one cell and overlap is not allowed. Cell-cell and cell-medium adhesionparameters govern the mechanical interactions of the cells. Strength of adhesion betweencells is higher if cell-cell adhesion energy is lower compared to cell-medium adhesion energy,and vice-versa.125Figure B.5: Synchronized oscillations in the low adhesion regime in a simulation as in FigureB.4 over 1000 Monte Carlo “time steps”. Cell-medium adhesion energy (40) is less thancell-cell adhesion energy (80) in the Hamlitonian H. Adhesion strength is low, which impliesless entrainment/synchrony.126Figure B.6: Synchronization in the intermediate adhesion regime. Cell-medium adhesionenergy (80) is equal to cell-cell adhesion energy (80) in the Hamlitonian H.127Figure B.7: Synchronization in the high adhesion regime. Cell-cell adhesion energy (60) isless than cell-medium adhesion energy (80) in the Hamlitonian H. This implies high degreeof adhesion strength between cells, leading to entrainment.Figure B.8: Relaxed cell with low Rho GTPase activity, β = 0.05. Cells are coloured byGTPase activity. Cell area, target area, and GTPase activity are plotted over time.128Figure B.9: Damped oscillations for β = 0.1. Cells are coloured by GTPase activity. Cellarea, target area, and GTPase activity are plotted over time.Figure B.10: Small amplitude limit cycle, β = 0.15. Cells are coloured by GTPase activity.Cell area, target area, and GTPase activity are plotted over time.Figure B.11: Stochastic switching between low amplitude limit cycle and high amplitudelimit cycle, β = 0.175. Cells are coloured by GTPase activity. Cell area, target area, andGTPase activity are plotted over time.129(a) Cell area in a 2D tissue over time.(b) GTPase activity in the same 2D tissue over time.Figure B.12: As in Figure 3.7, but in the weak adhesion scenario. In (A), cells are colouredbased on their current cell area, while in (B), cells are coloured based on the uniform levelof GTPase activity within each cell. In the Hamiltonian, H, cell-medium adhesion energy(60) is less than cell-cell adhesion energy (80). Notice that some cells detach from the tissuedue to low adhesion strength. Patches of synchronized cell oscillations are still observed.130
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- From signalling to cell behaviour : modelling multi-scale...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
From signalling to cell behaviour : modelling multi-scale organization in single and collective cellular… Zmurchok, Cole 2018
pdf
Page Metadata
Item Metadata
Title | From signalling to cell behaviour : modelling multi-scale organization in single and collective cellular systems |
Creator |
Zmurchok, Cole |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Individually and collectively, cells are organized systems with many interacting parts. Mathematical models allow us to infer behaviour at one level of organization from information at another level. In this thesis, I explore two biological questions that are answered through the development of new mathematical approaches and novel models. (1) Molecular motors are responsible for transporting material along molecular tracks (microtubules) in cells. Typically, transport is described by a system of reaction-advection-diffusion partial differential equations (PDEs). Recently, quasi-steady-state (QSS) methods have been applied to models with linear reactions to approximate the behaviour of the PDE system. To understand how nonlinear reactions affect the overall transport process at the cellular level, I extend the QSS approach to certain nonlinear reaction models, reducing the full PDE system to a single nonlinear PDE. I find that the approximating PDE is a conservation law for the total density of motors within the cell, with effective diffusion and velocity that depend nonlinearly on the motor densities and model parameters. Cell-scale predictions about the organization and distribution of motors can be drawn from these effective parameters. (2) Rho GTPases are a family of protein regulators that modulate cell shape and forces exerted by cells. Meanwhile, cells sense forces such as tension. The implications of this two-way feedback on cell behaviour is of interest to biologists. I explore this question by developing a simple mathematical model for GTPase signalling and cell mechanics. The model explains a spectrum of behaviours, including relaxed or contracted cells and cells that oscillate between these extremes. Through bifurcation analysis, I find that changes in single cell behaviour can be explained by the strength of feedback from tension to signalling. When such model cells are connected to one another in a row or in a 2D sheet, waves of contraction/relaxation propagate through the tissue. Model predictions are qualitatively consistent with developmental-biology observations such as the volume fluctuations in a cellular monolayer. The model suggests a mechanism for the organization of tissue-scale behaviours from signalling and mechanics, which could be extended to specific experimental systems. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-07-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0368799 |
URI | http://hdl.handle.net/2429/66413 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_september_zmurchok_cole.pdf [ 8.38MB ]
- Metadata
- JSON: 24-1.0368799.json
- JSON-LD: 24-1.0368799-ld.json
- RDF/XML (Pretty): 24-1.0368799-rdf.xml
- RDF/JSON: 24-1.0368799-rdf.json
- Turtle: 24-1.0368799-turtle.txt
- N-Triples: 24-1.0368799-rdf-ntriples.txt
- Original Record: 24-1.0368799-source.json
- Full Text
- 24-1.0368799-fulltext.txt
- Citation
- 24-1.0368799.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0368799/manifest