Cell size, mechanical tension, and GTPase signaling in the Single Cell
mmanuscript No. (will be inserted by the editor)
Cell size, mechanical tension, and GTPase signaling in the Single Cell
Andreas Buttenschön · Yue Liu · Leah Edelstein-Keshet
Received: date / Accepted: date
Abstract
Cell polarization requires redistribution of specific proteins to the nascent front and back of aeukarytotic cell. Among these proteins are Rac and Rho, members of the small GTPase family that regulatethe actin cytoskeleton. Rac promotes actin assembly and protrusion of the front edge, whereas Rho activatesmyosin-driven contraction at the back. Mathematical models of cell polarization at many levels of detail haveappeared. One of the simplest based on “wave-pinning”, consists of a pair of reaction-diffusion equations fora single GTPase. Mathematical analysis of wave-pinning so far is largely restricted to static domains in onespatial dimension. Here we extend the analysis to cells that change in size, showing that both shrinkingand growing cells can lose polarity. We further consider the feedback between mechanical tension, GTPaseactivation, and cell deformation in both static, growing, shrinking, and moving cells. Special cases (spatiallyuniform cell chemistry, absence or presence of mechanical feedback) are analyzed, and the full model isexplored by simulations in 1D. We find a variety of novel behaviors, including “dilution-induced” oscillationsof Rac activity and cell size, as well as gain or loss of polarization and motility in the model cell.
Keywords
Cell size, motility, GTPases, Rac and Rho, cell polarization, reaction-diffusion equations,growing 1D domain
Acknowledgements
LEK and coauthored are (partially) supported by NSERC (Natural Sciences and Engineering ResearchCouncil) Discovery grant 41870 to LEK. YL was partially supported by an NSERC postgraduate fellowship. AB was partiallysupported by an NSERC PDF Fellowship. We are grateful to the Pacific Institute for Mathematical Sciences for providing spaceand resources for AB’s postdoctoral research. The authors thank C. Zmurchok and W. R. Holmes and E. Rens for valuablediscussions.A. ButtenschönDepartment of Mathematics, University of British Columbia, Vancouver V6T 1Z2, BC, CanadaE-mail: [email protected]. LiuDepartment of Mathematics, University of British Columbia, Vancouver V6T 1Z2, BC, CanadaE-mail: [email protected]. Edelstein-KeshetDepartment of Mathematics, University of British Columbia, Vancouver V6T 1Z2, BC, CanadaE-mail: [email protected] a r X i v : . [ q - b i o . CB ] A ug Andreas Buttenschön et al.
Cell polarization and motility play a vital role in living organisms across all size scales. In normal devel-opment, wound healing, and immune surveillance, cells undergo directed migration tightly controlled byintracellular signaling pathways. Hence, these cellular phenomena have attracted attention of modelers, andgained representation at a whole spectrum of levels of detail, in 1D (Otsuji et al. 2007; Mori et al. 2008;Dawes and Edelstein-Keshet 2007; Holmes et al. 2012), 1.5D (the edge of a 2D domain) (Meinhardt 1999;Neilson et al. 2011), as well as 2D (Marée et al. 2006; Wolgemuth and Zajac 2010; Marée et al. 2012; Camleyet al. 2013), and 3D (Cusseddu et al. 2018) geometries.An elementary model for cell polarization by “wave-pinning” was introduced in (Mori et al. 2008) in1D, based on a single signaling protein (GTPase) inside a stationary cell of fixed length. The wave-pinningmodel was a simplification of the more biologically detailed models of GTPase role in polarity captured in(Marée et al. 2012; Jilkine and Edelstein-Keshet 2011). Generalization of this model to a pair of mutuallyantagonistic (Rac-Rho) GTPases and use of approximation methods demonstrated coexistance of cell statessuch as polarized, low uniform and high uniform protein activity levels (Holmes and Edelstein-Keshet 2016).As interesting aspect of the GTPases Rac and Rho is that they tend to concentrate at cell edges, wherethey drive actin assembly and edge protrusion (Rac) or actomyosin contraction (Rho), which moves the celledge outwards or inwards. The cell tends to resist deformation, which means that when one edge moves,the opposite edge is pulled or pushed accordingly. At the same time, protrusion or retraction may cause cellvolume to change, particularly in combination with aquaporin channels that allow water to flow in or outof the cell. Hence, cells can undergo fluctuations in size (Zehnder et al. 2015) that potentially affect theirinternal concentration.In (Mori et al. 2008, 2011; Jilkine and Edelstein-Keshet 2011; Holmes and Edelstein-Keshet 2016), cellsize was assumed to be static, whereas here we lift this restriction. The effect of cell size was briefly consideredin modeling the polarization response of HeLa cells to Rac stimulation in microfluidic channels (Holmes et al.2012). It is also incorporated in a companion paper by (Zmurchok and Holmes 2019), where the focus is oncell shape diversity arising from autoactivation as well as mutual inhibition in Rac-Rho interactions. Cellswere deformable in (Zmurchok et al. 2018), but possibility dilution was neglected, whereas here we considersuch effect. Finally, in (Zmurchok et al. 2018; Bui et al. 2019), feedback from tension due to cell deformationwas included, but only in model cells with spatially uniform GTPase activity. Here we generalize to fullspatially distributed GTPase activity, which introduces new interesting behavior.Mathematical analysis is possible for simple geometries and models with few components, and becomesincreasingly forbidding as the size of the signaling network or the complexity of the geometry increases. Forthis reason, we concentrate on a toy model in 1D geometry to describe the interplay between cell size, cellpolarization/motility and mechanical feedback. We show that this interplay can lead to cycles of expansionand contraction in certain cases and/or initiation of motility after some delay time.The model is simple enough that some analytical results can be written to characterize regimes of be-haviour. We can apply several methods, including linear stability analysis, approximations of kinetic termsby piecewise linear functions (where steady states can be found in closed form), as well as bifurcation analy-sis and numerical PDE solutions. Using these methods, we here generalize and extend some of the previousmathematical findings about wave-pinning, and help to increase insight about how such mechanisms for cellpolarization could operate.
We briefly outline previous work (Mori et al. 2008; Holmes and Edelstein-Keshet 2016) as a starting pointfor our paper. Concentrations of active ( G ) and inactive ( G i ) GTPases are modeled by reaction-diffusion ell size, mechanical tension, and GTPase signaling in the Single Cell 3 Mechanicalfeedback CellprotrusionCellcontraction
Polarization model
Fig. 1
The model considered here has two key modules: the polarization model, and the one-dimensional visco-elastic cellmodel. Here a Kelvin-Voigt element is placed between the cell endpoints ( x ( t ) and x ( t )). While other visco-elastic models arepossible, here we focus on the interplay between mechanics and signalling. The polarization model captures the spatio-temporaldynamics of one or more small GTPases affecting the cell’s cytoskeleton. Here their effects i.e. assembly of cell protrusions (Rac)or cell contractions (Rho) are captured by applying forces on the cell’s endpoints. Here we consider three different signallingmodules. For a spatially uniform single GTPase module (either Rac or Rho) we present detailed phase-plane and bifurcationanalysis in Section 7. In particular, we demonstrate the existence of “dilution-induced” oscillations of Rac concentration and cellsize. In Section 8 we show that a spatially heterogenous GTPase (Rac or Rho) allows the cell to polarize and move. However,if a cell is to soft i.e. it undergoes large size changes, it may undergo “dilution-induced” depolarize. Finally, in Section 8.2 weconsider a polarization module consiting of antagonistic Rac and Rho. We capture all features of the single GTPase models,but also gain more complicated spatio-temporal patterns. (RD) equations, G t = D∆G + A ( G ) G i − IG, ( G i ) t = D i ∆G i − A ( G ) G i + IG. (1)Active GTPase is membrane-bound, so the rates of diffusion satisfy D (cid:28) D i . A ( G ) is a rate of activationand I a rate of inactivation. Feedback from the active GTPase to its kinetics is assumed to be directedthrough the rate of activation of GTPase by GEFs (guanine nucleotide exchange factors), whereas the rateof inactivation by GAPs (GTPase activating proteins) is taken as constant. Since the GTPase cannot enteror leave through the cell edge, the PDEs are supplemented by Neumann (no-flux) boundary conditions.Equations (1) conserve the total amount of GTPase ( G T ) in a cell, that is G T := Z ( G + G i ) d x. As in previous papers, we consider both the single GTPase (“wave-pinning”) model of (Mori et al. 2008), aswell as the Rac-Rho model of (Holmes and Edelstein-Keshet 2016).
Single GTPase model:
The standard functional form for the rate of activation, A ( G ), is A ( G ) = (cid:18) β + γ G n G n (cid:19) . The Rac-Rho (two GTPase) model:
Let G R , G ρ represent the activities of Rac and Rho GTPases, re-spectively. We assume mutual inhibition. That is, the rate of activation of Rac is reduced by Rho andvice versa. The standard function form for the activation rate of Rac is A R ( G ρ ) = (cid:18) β + γ
11 + G nρ (cid:19) . Andreas Buttenschön et al.
Fig. 2
Left: Top-down view of a eukaryotic cell. Middle: the nucleus has been removed, so that the cell fragment is a thin sheetwith approximately uniform thickness ( < µ m). A 1D transect across the cell is shown. Right: the front and rear positions ofthe cell edge are labeled. We consider cell motion such that x , x , L are all possibly time-dependent. The region x ≤ x ≤ x is the domain inside which GTPase reaction-diffusion takes place. Mechanically, the 1D cell is represented as a spring dashpot(Kelvin-Voigt) element. Fig. 3
A summary of the main results of this paper. We discuss how changes in cell size can result in dilution of internalchemical concentration or activity (a). We show that size fluctuations can lead to loss (or gain) of chemical polarization bywave-pinning, as in (Mori et al. 2008) (b). We then introduce mechanical force of protrusion (Rac) or contraction (Rho). Inthe uniform Rac case (c), cell size fluctuations can occur by virtue of dilution (and loss of Rac activity) when the cell expands.Finally, we show how the mechanochemical model linking GTPase activity to cell protrusion can lead to cell motility. and similarly for the dependence A ρ ( G R ) of Rho activation rate on the activity of Rac (Holmes andEdelstein-Keshet 2016). The geometry for the single-cell model in Fig. 2 represents a transect across the diameter of the cell.We keep track of the coordinates for the front and rear cell edges x ( t ) , x ( t ), where (WLOG) x ( t )
There are two new terms that appear, namely an effective “advection term”: v · ∇ u which corresponds toelemental volumes moving with the flow due to local growth and a dilution term, u ∇ · v due to the localvolume change.For example, if only the centroid moves e.g. x c ( t ) = βt , and R ( t ) = R . Then v = β , so no dilutiontakes place. Alternately, if only the radius of the cell increases with x c ( t ) = 0, and R ( t ) = R exp ( kt ), then v ( x, t ) = kx . In this case we have both advection and diffusion. In particular, v · ∇ u = kx ∇ u, and u ∇ · v = uk. Refer to Fig. 4 for an example of this case.3.2 Mapping to a constant domain for numerical simulationsTo create numerical simulations of our 1D model (while avoiding the complexities of a moving boundaryproblem), it proves helpful to map the moving boundary problem (2) defined on Γ ( t ) to a stationary domain¯ x ∈ [ − , x, t ) (¯ x, ¯ t ) = (cid:18) x − x c ( t ) R ( t ) , t (cid:19) . (3)In Appendix B, we show that the resulting equation for u is u ¯ t = DR ( t ) u ¯ x ¯ x + f ( u ) − u (cid:18) ˙ R ( t ) R ( t ) (cid:19) , in [ − , . (4)Interestingly, the change of coordinates eliminates the advection term, while keeping the dilution term. Atthe same time, the diffusion coefficient becomes time dependent. The advantage of equation (4) is thatit is formulated on a fixed domain. This makes its numerical solutions using a method of lines approachstraightforward compared to equation (2) (see Appendix B).The single GTPase model, so transformed to the domain [ − , G ¯ t = DR ( t ) G ¯ x ¯ x + A ( G ) G i − IG − D ( t ) G, in [ − , , (5a)Inactive: ( G i ) ¯ t = D i R ( t ) ( G i ) ¯ x ¯ x − A ( G ) G i + IG − D ( t ) G i , in [ − , , (5b)where D ( t ) = (cid:18) ˙ R ( t ) R ( t ) (cid:19) = (cid:18) ˙ L ( t ) L ( t ) (cid:19) is an effective “dilution term”, and where A ( G ) is the activation rate, as before. In the case that cell volumeis simply redistributed, while only the membrane component expands or contracts, the term containing D ( t )would be omitted from the second equation (see also Appendix D). The dilution term can be expressedalternately in the above radius or length-dependent forms, used interchangeably in what follows. ell size, mechanical tension, and GTPase signaling in the Single Cell 7 G d t = (cid:18) β + γ G n G n (cid:19) G i − IG − G (cid:18) ˙ L ( t ) L ( t ) (cid:19) , d G i d t = IG − (cid:18) β + γ G n G n (cid:19) G i − G i (cid:18) ˙ L ( t ) L ( t ) (cid:19) . Using mass conservation we reduce this set of two equations to one equation by eliminating the inactiveGTPase variable: G i = G T L ( t ) − G, (6)where G T is the total GTPase. This leads tod G d t = (cid:18) β + γ G n G n (cid:19) (cid:18) G T L ( t ) − G (cid:19) − IG − G (cid:18) ˙ L ( t ) L ( t ) (cid:19) . (7)We later show that this case has interesting consequences when coupled with mechanics. A model of thistype, but without the dilution effect was also studied in (Holmes and Edelstein-Keshet 2016). We first consider results when domain size affects the signaling, but without the effect of signaling on themotion of the cell edges. We show how the classic wave-pining is modified in such cases.4.1 Parameter valuesSince some of our results require numerical simulations, we need biological parameters for the GTPasesignaling equations. Typical parameter for the GTPase dynamics are shown in Table 1.
Parameter Symbol Value Range Reference
Basal activation rate β [s − ] 0 . − . γ [s − ] 1 −
25 Jilkine and Edelstein-Keshet; Holmes and Edelstein-KeshetDecay-rate I [s − ] 1 Jilkine and Edelstein-KeshetEC parameter G [µ m ] 1 Michaelson et al.; Mori et al.; Jilkine and Edelstein-KeshetActive GTPase diffusion D [µm · s − ] 0 . D i [µm · s − ] 100 Jilkine and Edelstein-KeshetTotal GTPase G T [µ m ] 1 − n −
25 Holmes and Edelstein-Keshet
Table 1
Parameters for the GTPase intracellular signaling models. Note that G is the GTPase level for half-max GTPaseautoactivation in the Hill function of equation (7), prior to non-dimensionalization. A range of up to 1000-fold has been assumedfor the diffusion of GTPase protein in the membrane (active) and cytosol (inactive). The upper bound for the Hill coefficient ischosen such that there is a sufficiently large regime in which oscillations are possible. Andreas Buttenschön et al. Fig. 4
Effect of domain size on polarization. The GTPase model with changing length (simulated using the mapping to [ − , L = 55µm and parameters values: n = 4 , D i = 1 , D = 0 . , γ = 2 , b = 0 . , I = 1 , G T = 2 . R ( t ) = R exp ( − kt ), as the cell shrinks G T increases. Middle: A static domain with R ( t ) = R .Right: A growing domain with R ( t ) = R exp ( kt ), (with k = 10 − ): as the cell grows, G T decreases. In both cases in whichthe domain size changes, G T eventually leaves the parameter regime in which wave pinning is possible. The transition is veryrapid. T = 800. Note the rapidityof the transition from polarized to uniform activity. Loss of polarization if the cell becomes too small isconsistent with bifurcation results of (Mori et al. 2011). Loss of polarization when the cell grows too muchcorresponds to the case where mean total GTPase concentration in the cell drops, also consistent with (Moriet al. 2011). In short, either increase or decrease of cell size eventually moves the system out of the regimeof parameters in which wave-pinning can be sustained.Up to this point, we considered only the GTPase distributions inside the cell, and the effect of changesin cell size on that internal level. To describe the feedback between chemistry and mechanics, we next endowthe cell with additional properties that link its size to mechanical forces. The basic ideas follow those of(Zmurchok et al. 2018), but with the adjusted mass balance in place. We now specify assumptions and equations for the endpoints of the cell, which then determine its size. Asshown in Fig. 2, the cell is modeled as a 1D contractile element of length L and rest length L , analogous toa Kelvin-Voigt element with constitutive relation for the stress, σ given by σ ( t ) = E(cid:15) ( t ) + η d(cid:15)dt , ell size, mechanical tension, and GTPase signaling in the Single Cell 9 where E is the cell’s Young modulus, η the viscous component and (cid:15) the strain. Cell motion is in anoverdamped regime where inertial terms in the Newtonian force balance are negligibly small. In this regime,the equation of motion becomes F = µv , where µ is the cell’s frictional coefficient. We denote the forcesapplied at nodes x , x by F , F respectively. The cell’s strain is given by ε = L − L = x − x − L . Parameter Symbol Value Range Reference
End point friction µ [pN · s · µm − ] 10 − Drasdo and HöhmeCell viscosity η [pN · s · µm − ] 6 × PalssonCell Young’s modulus E [pN · µm − ] 5 × Drasdo and Höhme; PalssonCell rest length L [µm] 10 − F max [pN] 10 − Cell force half-rate constant G c [µ m ] 1 . G Cell force Hill coefficient m − Table 2
Mechanical parameters for the model cell. G is defined in Table 1. Parameter values from 3D models are used toinform our 1D model. x , x , satisfy the following equationsof motion µ d x d t = F + E ( x − x − L ) + η (cid:18) d x d t − d x d t (cid:19) , (8a) µ d x d t = F − E ( x − x − L ) − η (cid:18) d x d t − d x d t (cid:19) . (8b)Solving the above equation for the time derivatives, we obtain two differential equations for the cell’sendpoints ˙ x ˙ x = 1 µ ( µ + 2 η ) µ ( F + Eε ) + η ( F + F ) µ ( F − Eε ) + η ( F + F ) . (9)5.2 Equivalent system for centroid and radiusWe can also write this system in terms of the cell’s centroid, x c ( t ) = ( x + x ) /
2, and cell’s radius, R ( t ) = L/ x − x ) / ˙ x c ˙ R = 12 µ ( µ + 2 η ) ( µ + 2 η )( F + F ) µ ( F − F − Eε ) , where ε = 2 R − L . In simplest form, the latter is written˙ x c = F + F µ , ˙ R = F − F − E (2 R − L )2( µ + 2 η ) . From the second equation we recognize the familiar spring-length equation (also directly obtained bysubtracting (8)(b-a)): ( µ + 2 η ) d L d t = ( F − F ) − E ( L − L ) . (10)In the Appendix, we also show that imposing external forces is equivalent to changing the rest-length ofthe spring, as in (Zmurchok et al. 2018; Liu 2019).5.3 Special casesFrom the equations for x c ( t ) , R ( t ), it is easy to see that in the case of equal forces on the two cell ends, weget F = F = F ⇒ ˙ x c = F/µ, ˙ R = 0 , so the cell moves at speed F/µ without deforming. In the case of equal and opposite forces on the two cellends F = − F = | F | ⇒ ˙ x c = 0 , ˙ R = ( F − Eε ) / ( µ + 2 η ) . Hence, the cell does not move, but its size changes.
We now link signaling inside the cell to cell mechanical properties. In (Zmurchok et al. 2018), this couplingwas done for the spatially uniform GTPase case without the dilution terms. We extend that model byconsidering a fully nonuniform GTPase distribution inside a cell correct mass balance for cell size changes.6.1 GTPase activity affects cell contraction and expansionThe Rac-induced cell protrusion and Rho-induced cell contraction can be modeled in various ways. Here wewill assume that Rac and Rho (by assembling actin or activating actomyosin) create forces at cell endpoints.These forces either push the endpoint outwards (Rac) or pull it inwards (Rho). Forces are taken to be afunction of the local GTPase level. In the spatially dependent GTPase cases, each cell endpoint can experiencea distinct force.Denote by G ∈ { G R , G ρ } the active form of either GTPase. Then the forces F , on the cell’s endpointsare taken to be F k ( t ) = ± F max G m ( x k , t ) G mc + G m ( x k , t ) , (11)where k = 1 ,
2. The sign of the force would correspond to an outwards force for Rac and an inwards forcefor Rho.According to (11), force builds up continuously as G ( x k ) approaches a characteristic level, G c . Theforce saturates to a maximal level F max and the sharpness of the response increases with m . (Note that ourframework here appears to be different from (Zmurchok et al. 2018; Liu 2019), where GTPase activity affectedthe cell rest-length; we show the mathematical equivalence of these two model variants in the Appendix.) ell size, mechanical tension, and GTPase signaling in the Single Cell 11 G ¯ t = DR ( t ) G ¯ x ¯ x + f ( G, G i ) − G (cid:18) ˙ R ( t ) R ( t ) (cid:19) , in [ − , , (12a)Inactive GTPase: ( G i ) ¯ t = D i R ( t ) G i ¯ x ¯ x − f ( G, G i ) − G i (cid:18) ˙ R ( t ) R ( t ) (cid:19) , in [ − , , (12b)Kinetics, BCs: f ( G, G i ) = (cid:18) β + γ G n G n (cid:19) G i − IG, G ¯ x = G i ¯ x = 0 for ¯ x = ± , (12c)Cell centroid, radius: ˙ x c = F + F µ , ˙ R = F − F − E (2 R − L )2( µ + 2 η ) , (12d)Forces on cell ends: F = ± F max G m ( − G mc + G m ( − , F = ± F max G m (1) G mc + G m (1) . (12e)In the following sections, we study variants of this model with and without the diffusion terms, and withthe protrusive (Rac, G R ) and contractile (Rho, G ρ ) GTPases acting on their own and together. In this case, we drop the diffusion terms in the first two PDEs in Eqs. (12), so that active GTPase G is timedependent only, uniform throughout the cell. The conservation statement for GTPase reduces to Z L ( t )0 ( G ( t ) + G i ( t )) dx = L ( t )[ G ( t ) + G i ( t )] = G T > G T is the constant total amount of GTPase in the cell. Since G, G i > G ( t ) < G T L ( t ) . This implies that the region of interest in the L − G plane is bounded above by the curve G = G T /L . Indeedwe show below (Result 1) that this forms an invariant region, i.e. that trajectories do not cross this curve.Furthermore, since G is spatially uniform, and hence the same at both cell ends, we have that forces on theends are equal in magnitude but opposite in sign i.e. F = − F , and the cell tension, F T ( G ) = F ( G ) − F ( G ),satisfies | F T | = 2 | F | . This also means that the cell’s centroid remains fixed i.e. ˙ x c = 0.Coupling the length of the cell to the GTPase dynamics, Eqn. (7), leads to the reduced model variant tostudy: d L d t = − E ( L − L ) + F T ( G ) µ + 2 η , (13a) F T ( G ) = ± F max G m G mc + G m , (13b)d G d t = (cid:18) b + γ G n G n (cid:19) (cid:18) G T L − G (cid:19) − IG − G (cid:18) ˙ LL (cid:19) . (13c) The sign on the force terms in (13b) is + for a protrusive GTPase ( G = G R ) and − in the contractileGTPase case ( G = Gρ ). We have reduced the original system (12) to the two coupled ODES (13a) and (13b)for L ( t ) , G ( t ). This system can be studied in the L − G phase plane. Below we explore the detailed behaviorof this system.To ensure that the cell length remains positive in the contractile (Rho) case we require that EL /F max >
1. In addition, we assume that in this section the force switches “ON” at G c >
1; the case in which G c < Result 1 (Invariant Region)
The region defined by ∆ := { ( L, G ) :
L > G ≥ G ≤ G T /L } is invariant with respect to the flow given by equations (13) .Proof Since EL /F max >
1, we have that L > L >
0, and similarly G > G = 0. So flowcannot leave the region through the axes. Now consider the direction of the flow along the boundary curve f ( L ) = G T /L . The normal vector along the curve is given by (1 , L /G T ), then computing the component offlow along this normal direction, we find( L ( t ) , G ( t )) · (1 , L /G T ) < . so the flow is into the region, and cannot leave through this boundary. Hence flow is trapped in the region ∆ , which is thus invariant.In Fig. 5, we show typical phase plane behavior of the system (13) for the contractile Rho ( G = G ρ , left)and protrusive Rac ( G = G R , right) cases. In both cases we find the region ∆ bounded by the (thick black)curve G = G T /L . The level of total GTPase, G T controls the position of this curve, so that the region islarger for larger values of G T . The nullclines, obtained by setting d L/ d t = 0 in (13a) and d G/ d t = 0 , ˙ L = 0in (13c) have shapes with the following features: the L nullcline (black) is sigmoidal, but its orientation flipsbetween the two cases. The G nullcline (grey) is “Z-shaped”. It remains to determine how many possiblesteady states can occur at intersections of such curves, and which parameters govern possible bifurcations.7.1 Sharp-switch analysisTo help explore this question, we use a common shortcut, the sharp-switch approximation (Holmes andEdelstein-Keshet 2016) in the L equation. Next we estimate the number of possible steady states. Result 2 (Number of steady states)
Consider ODE (13c) with L ( t ) set to either of the two steady statesof the ODE (13a) . Then ODE (13c) has either one or three positive steady states.Proof In the contractile (Rho) case, the L nullcline is given by: L eq = L − F max E G m G mc + G m . In the protrusive (Rac) case, the minus sign is replaced by +.We assume that the switch is sharp and approximate L by L ∞ ∈ { L , L − F max /E } . The nullcline for G can be written as: L [ A ( G ) + I ] + ˙ L = A ( G ) G T G , L, G > , where A = ( b + γh ( G )) and h ( G ) is the Hill-function in (13a). Substituting the expression for ˙ L from theODE (13c), leads to a more complicated nullcline equation which we do not display, but which leads to the ell size, mechanical tension, and GTPase signaling in the Single Cell 130 . . . . . . . G ]3 . . . . . . . . C e ll L e n g t h [ L ] Contractile-Case (Rho) . . . . . . . G ]4 . . . . . . . . . C e ll L e n g t h [ L ] Protrusive-Case (Rac)
Fig. 5
Phase-plane plots for the contractile (left: Rho) and protrusive (right: Rac) cases of the dynamical system (13). In solidblack L nullcline, and in grey the G nullcline. The steady states are shown by red dots. Note the invariant region, bounded bythe curve L = G T /G . Here G T = 15 , G c = 1 . , γ = 2 . , n = 25 , m = 16 , L = 6, and other parameter values are as in Table 1,and Table 2. In both cases, the most left ( G ) and most right ( G ) steady states can be approximated using a sharp switchapproximation, the middle steady states are saddles. sigmoidal curve plotted in Fig. 5 (black) and in Fig. 6 (blue). Observe that the direction of this switch flipsbetween the Rac and the Rho cases, as shown in the top and bottom panels of Fig. 6, respectively.As we are seeking steady states, we may set ˙ L = 0, simplifying the resulting equation to L = G T ( b + γh ( G )) G ( I + b + γh ( G )) ≡ n ( G ) , where we have defined n ( G ) as the rational function that appears. We note that asymptotically n ( G ) behavesas: lim G → n ( G ) = ∞ , lim G →∞ n ( G ) = 0 . To find intersections of n ( G ) with L ∞ we must solve: G n +1 L ∞ ( I + b + γ ) − G n G T ( b + γ ) + GL ∞ ( I + b ) − bG T = 0 . Let p be the number of positive solutions of this polynomial. Since the sequence of the polynomial’s coefficientsfeatures three sign changes, we have by Descartes’ rule that 3 − p must be non-negative and even (for astatement of Descartes’ rule we refer the reader to Appendix 2 of (Murray 1989)). We conclude that p = 1or 3.Result 2 implies that in most cases we have either one or three solutions. Note that in the most generalcase when L is not set equal to either of its two steady-states, Descartes’ implies that there may be 1 , , L = L ∞ these additional . . . . Rho concentration [ G ] C e ll L e n g t h [ L ] Low GTPase [ G T = 4 ] L nullcline G nullcline Rho concentration [ G ] C e ll L e n g t h [ L ] Medium GTPase [ G T = 12 ] L nullcline G nullcline Rho concentration [ G ] C e ll L e n g t h [ L ] High GTPase [ G T = 225 ] L nullcline G nullcline . . . . Rac concentration [ G ] C e ll L e n g t h [ L ] Low GTPase [ G T = 4 ] L nullcline G nullcline Rac concentration [ G ] C e ll L e n g t h [ L ] Medium GTPase [ G T = 50 ] L nullcline G nullcline . . . Rac concentration [ G ] C e ll L e n g t h [ L ] High GTPase [ G T = 150 ] L nullcline G nullcline Fig. 6
Possible nullcline configurations. Top: contractile GTPase (Rho) case; Bottom: protrusive GTPase (Rac) case. Top: ForRho, the left panel shows the existence of G . The middle panel shows the bistable case, i.e. two steady states in which the cellis large and GTPase is small ( G , and G ). The right panel shows the existence of a single steady state G . Other parameters γ = 2 . , b = 0 . , I = 1 , L = 6. The dotted lines are the upper and lower sharp switch approximations of the sigmoidal L nullcline. Bottom: For Rac, the middle panel demonstrates the case where a stable limit cycle oscillation is expected, a situationthat cannot occur in the Rho case. Other parameters γ = 0 . , b = 0 . , I = 1 , L = 6. For the values of the approximate steadystates G , G , and G see Appendix C.1. steady states must occur in the thin region G ∼ G c in which L transitions between the two possible valuefor L ∞ . Thus in the vast majority of cases we only have 1 or 3 steady states. See also Fig. 6.The sharp-switch approximation leads to explicit values for some of the steady states of the system,their eigenvalues, and hence full characterization of their stability and conditions for existence. The detailedcalculations are provided in Appendix C. A brief summary is given in Table 3. This approximation does notallow easy computation of steady states that occur when nullclines intersect along the upswing portion ofthe switch. However, the geometry of such intersections allows us to determine stability nonetheless, usinggeometric methods described in (Murray 1989) (see also Appendix C for linear stability analysis).In the contractile GTPase (Rho) case, we find three steady states, all stable nodes, with possible co-existence of all three (These are separated by steady states that occur on the portion of the switch notcaptured by the approximation). In the protrusive (Rac) case there are either two coexisting stable steady ell size, mechanical tension, and GTPase signaling in the Single Cell 15 GTPase Level L eq G eq Stability
G < min(1 , G c ) L bG T L ( b + I ) stable nodetransition ? ? Rho: saddle. Rac: unstable, saddle1 < G < G c L ( b + γ ) G T L ( b + γ + I ) stable nodetransition ? ? Rho: saddle. Rac: unstable, saddle G > max( G c , L ( b + γ ) G T L ( b + γ + I ) stable node Table 3
Sharp switch approximated steady state summary. Here L = L − ( F max /E ) for Rho and L = L + ( F max /E ) forRac. The transitions in the table can contain one or more steady states. In case of a single steady state at the transition, thatstate is unstable. Note that the number of coexisting steady states depends on the parameters. See text and Appendix C fordetails. states (separated by a steady state not captured by the approximate) or no coexisting stable steady statesat all. When no stable steady states exist we must have a stable limit cycle, since ∆ is invariant.7.2 SimulationsWe simulated the model for the spatially uniform Rho, and Rac case (see Fig. 7), starting with a cell at itsrest-length. Results are shown in the form of a kymograph, with time on the horizontal axis, and a profileof the cell and its internal active GTPase concentration on the vertical axis.In both the Rho and Rac case of the single-GTPase model, when the total amount of the GTPase, G T inthe domain is too low (left panels), there is no response. For high G T (right panels), the system is triggeredto switch to the high GTPase steady state (shaded yellow). For Rho (top panels in Fig. 7), the cell contractsas a result, and the high activity state becomes amplified even more due to the smaller volume in which it isconcentrated. In contrast, in the Rac case (Bottom panels of Fig. 7), expansion caused by high Rac resultsin a larger cell, but also dilutes Rac activity, depressing it to a lower level.In the Rac, but not the Rho case, the conflicting effects of expansion and dilution can set up an oscillatorybehavior (bottom middle panel, Fig. 7). The cell alternates between long and short phases, with fluctuationsin its Rac activity.7.3 Bifurcation diagramsFinally, we characterized the regimes of behavior in the spatially uniform cases by creating bifurcation plotsof model (13). Results are shown in panels of Fig. 8 for Rho (top) and Rac (bottom). It is evident that aregion of tristablity can exist in the Rho case (labeled C), whereas the Rac case has a regime in which limitcycle oscillations can exist. Panels on the left should be compared with the 2-parameter bifurcation plot forthe single-GTPase case (Fig. 3a) in (Holmes and Edelstein-Keshet 2016). The distinction here is that wehave taken cell size into account, which increases the number of possible steady states and introduces thepossibility of Hopf bifurcations. Fig. 7
Simulation of a cell with uniform internal GTPase. Top: Rho, the GTPase that leads to contraction. Bottom: Rac, theGTPase that leads to cell protrusions. Left panels: Low total amount of GTPase, G T = 2 .
5: no change takes place in eithercase. Right panels: High total GTPase, G T = 50: the cell becomes uniformly high in Rho/Rac and hence shrinks/expands,respectively. Bottom Middle: Cell with intermediate total Rac and parameters such that oscillation is possible. Same parametersas the bottom middle of Fig. 6 ( E = 250 [pN · µm − ]). Initially, the cell is at its “rest-length” L (0) = L , and GTPase is uniformlyset to G T /L . In the Rac case the cell’s stiffness is E = 500 [pN · µm − ], while in the Rho case it is E = 1500 [pN · µm − ]. Thischoice ensures that the cell does not vanish i.e. EL /F max > ell size, mechanical tension, and GTPase signaling in the Single Cell 17 Fig. 8
Bifurcation diagrams for the spatially uniform GTPase cases. Top: the contractile (Rho) case, showing the possibility oftristability. Bottom: The protrusive (Rac) case. Left panels: one-parameter bifurcation with respect to the feedback activationrate, γ . Right panels: two-parameter bifurcation plot with respect to total GTPase, G T and γ . The fold bifurcations in the Rhocase (top left) labeled A, B, C, D, are tracked by curves with the same labels in the two-parameter plot (top right). Numbers1,2,3 indicate existence of the corresponding steady states in regions of the same plot. (For example, “1,2,3" indicates existenceof all three equilibria G , G , G , as discussed in Appendix C. The parameters are b = 1 , I = 5 , m = n = 20 , L = 3 , Eµ +2 η =5 , F max µ +2 η = 10 , G c = 4 , G T = 15. Bottom: Left: a supercritical Hopf bifurcation is evident at G T = 40. Right: the inset shows alimit cycle that exists in the range of G T = 65 and 3 . < γ < .
6. Limit cycles are possible only in the Rac case. The parametersare b = 1 , I = 10 , m = n = 8 , L = 5 . , Eµ +2 η = 5 , F max µ +2 η = 55 , G c = 0 . Here we numerically solve the full PDE model (12) with dynamic cell size, and forces created at the cellends by local GTPase. We take into account dilution, and examine both contractile and protrusive GTPase.Details of the numerical implementation is given in Appendix B.8.1 Single GTPaseResults are shown as kymographs in Fig. 9. Cells were initialized at rest-length, L (0) = L , with elevatedGTPase (step function) at one edge. We tested both stiff (left) and soft (right) cells, by varying the parameter E . Fig. 9
Full spatiotemporal simulations of a single GTPase in the 1D cell. (Shown is the concentration of active GTPase.) Top:the contractile (Rho) case. Bottom: the protrusive (Rac) case. Initially, cell length is L (0) = L . Active GTPase initialized withstep function at the top edge of the cell. Left: A stiff cell with G T = 25; length hardly changes in either case and polarizationis maintained. Right: A soft cell with G T = 10, showing significant length change. Once the cell shrinks or expands too much,it loses polarization. For the Rho case (top row of Fig. 9), motility is by “pushing” or “squeezing” from the rear. That is, asthe edge of the cell with high Rho activity contracts inwards, the cell’s elastic “spring” pushes the other edgeoutwards. If the cell is stiff, both edges then translocate in a nearly linear fashion. A soft cell never managesto compensates for the contraction of the rear edge. Eventually, the cell becomes too small to supporta polarized GTPase pattern. Then polarization is lost and the whole cell contracts and keeps shrinking,amplifying the Rho concentration.The Rac GTPase case is shown in the bottom row of Fig. 9. The cell moves in the direction of the Racgradient (edge protrudes outwards). As before, for stiff cells (left), the rear keeps up with the front due tothe restoring spring-force; a tight GTPase gradient is formed and persists in the moving cell. In the soft cell(right), the rear edge is only weakly pulled along by the elastic restoring force. Once the cell is “too long”,polarization is lost, Rac activity collapses abruptly to a low uniform level, and the cell shrinks.From the simulation results (and many replicates, not shown) we draw the following overall conclusions.(1) The direction of motion depends on initial conditions and GTPase type. Cells move towards the directionof Rac and away from the direction of Rho gradient in the single GTPase simulations. (2) Results demon-strate polarization (gain or loss) as well as size changes and the interplay between these. (3) Polarized cells ell size, mechanical tension, and GTPase signaling in the Single Cell 19 can become non-polar by growing or shrinking too much - both extremes drive the system away from thepolarization parameter regime. (4) Softer cells are particularly prone to depolarization as restoring cell sizetakes longer than the timescale of loss of polarization: see the bottom row of Fig. 9. (Equivalently, cells withlarge protrusive or contractile forces F max are also prone to depolarization.) (5) Cells with two peaks of Racat opposite ends tend to stretch and lose the ability to polarize. Two peaks of Rho cause cells to shrinktoo much and lose polarization. (6) Changing the cell’s internal viscosity does not strongly affect the aboveresult, but it changes the time scale of cell deformation.8.2 Mutual antagonism: Rac and Rho in the same cellFinally, we simulated the Rac-Rho mutually antagonistic full model as described in the “Two GTPase model”of Section 2. The forces of protrusion or retraction on the cell end x k are assumed to depend on both Racand Rho as follows: F k ( t ) = ( − k F max (cid:18) R m ( x, t ) G mc + R m ( x, t ) − ρ m ( x, t ) G mc + ρ m ( x, t ) (cid:19) (cid:12)(cid:12)(cid:12) x = x k , for k = 1 ,
2. In the full model equations (12), this replaces equations (12e). Note that since Rac and Rho aremutually antagonistic, one will dominate at each of the cell endpoints.Cells were initialized with the same parameters, spatially uniform GTPase, with superimposed Gaussiannoise in inactive Rho. Sample results are shown in Fig. 10.These results (including others, not shown) can be summarized as follows: (1) The evolving GTPasepattern is determined by the initial conditions. (2) Rac and Rho segregate to mutually exclusive regionsin the cell. (3) Polarized motile cells spontaneously gain active Rac at the “front” and Rho at the “rear”.(4) Cells with more than a single transition layer display rich dynamics, in particular as plateaus of activeGTPase transition from the interior to a cell edge. This causes some cells to stall, and others to start moving,as shown in Fig. 10. (5) Cells with high Rac at both ends get larger, whereas those with high Rho at bothends shrink. Such metastable states can persist for a long time.In short, the full model with Rac-Rho mutual inhibition displays a diverse set of behaviors not seen inmodel variants that neglect the full spatio-temporal GTPase distribution, or that assume static cell size.
Interplay between the activity of GTPases, the contraction or expansion of a cell, and resultant feedback toGTPase activation was previously explored in a related framework by (Zmurchok et al. 2018). Our approachhere differs in two main respects: (1) GTPase is no longer assumed to be spatially uniform, and (2) cell sizedynamics is considered in the mass balance of the GTPase.Accounting for the spatial distribution of the GTPases inside a cell leads to new features that were absentin (Zmurchok et al. 2018). These include (a) delay while the internal GTPase activity reorganizes spatially,(b) the possibility of polarization and net motion, rather than just cell expansion or contraction (c) caseswhere the cell stalls when plateaus of GTPase form at both cell edges, and (d) loss of polarization if the cellgets too large or too small.Despite the simplicity of the model, initially comprised of a single GTPase (with active/inactive forms)in 1 spatial dimension, we find a variety of possible behaviours even in the spatially homogeneous case,depending on the action of the GTPase. The simplicity of the model then permits analysis of steady states,stability, and Hopf bifurcations. Analysis is considerably more challenging in the more detailed models.There are two ways that cell size affects the results of bifurcation analysis. 1) through the definitionof length and time scales; and 2) by diluting or concentrating the internal signaling activity, as shown bydetailed mass balance, and known classically. This in turn, also affects whether the cell can have a polarized
Fig. 10
Various possible behaviours in the full Rac-Rho model with mutual inhibition. Shown is the profile of active Rac.(Active Rho profile is complementary.) Parameters are L = 10µm , m, n = 4 , γ = 5 , G T = 4. Plots differ in initial inactive Rho(small Gaussian noise about the uniform steady state); other GTPase initially uniform across the cell. (A) Polarized cell withRac front moves at constant rate. (B) Stationary cell with high Rac at both ends expands. (C) Cell with high Rho at both endscontracts without moving. (D) Cell with two peaks of Rac becomes uni-polar (at around 900 s) and continues to move in thesame direction. (E) Similar to D, but with reorganized Rac eventually stalling the cell. (F) Initially static cell with an internalregion of high Rac, which moves toward a boundary, thus mobilizing the cell. Insets: expanded views of regions of interest shownin grey rectangular portions of the trajectories. state. The single GTPase model of Eqs. (1) can sustain “polarized” wave-pinning solutions, but if the domaingrows (or shrinks) too much these are lost. This is consistent with (Mori et al. 2011) where the bifurcationparameter (cid:15) depended on cell size.We showed that the cell could undergo cycles of expansion and contraction in the case of a Rac-likeGTPase that pushes out the cell edges. Recent work on cell oscillations (in the context or confined epithelialsheets) is given in (Peyret et al. 2019), and a review of cycles of protrusion and retraction in experimentaland theoretical models of the cell edge is found in (Ryan et al. 2012). Cycles also occur in paper by (Dierkeset al. 2014), who consider a similar treatment of dilution (of myosin concentration) when the cell stretches,and a mechanical contractile unit driven by myosin. Unlike our model, their chemical equation is simplylinear return to a baseline level c , with a similar dilution term ( c/‘ )( d‘/dt ). Their oscillations stem from anassumption that the “spring” governing cell length has a nonlinear (cubic) term ( K ( ‘ )).One can argue that some cells deform and redistribute their protrusions while maintaining a fixed volume.For this reason, we also considered a model variant with fixed volume. Qualitatively, the two variants aresimilar, through some features differ. See Appendix for a comparison. It is easy to generalize the model (8)to the version in the recent paper by (Bui et al. 2019) who used distinct values of the frictional coefficient µ at the front and the rear of each cell. ell size, mechanical tension, and GTPase signaling in the Single Cell 21 While this study demonstrated the behaviour of a simple GTPase system in 1D, it prepares the groundfor 2D simulations where cell shape can be explored more fully. It also opens the way to considering how cellsinteract with one another and how mechanochemical feedback between those cells affects emergent behaviorat the scale of a tissue.While it is straightforward to carry out analysis in the uniform case, such analysis is more challengingin the spatially distributed cases. Since in our model the forces on the cell ends solely depend on the activeGTPase concentrations at the cells ends, a worth while question for future research is whether it is possibleto replace (in some limit) the full PDE by two ODEs tracking the amount of active GTPase at the cell ends.A distinct challenge is the emergence of more complicated transient dynamics (see Fig. 10), which persiston time-scales significant for cellular behaviour. This highlights the need to not only better understand theasymptotic behaviour of such equations but also their transient dynamics. This becomes even more urgentin view of future models containing more details and multiple cells.
A Mass balance on growing or shrinking domain
Consider a prototypical reaction diffusion equation for some concentration u ( x, t ) in conservative form, u t = −∇ · j ( u ) + f ( u ) , in Γ ( t ) , where j ( u ) is the particle flux. Further, assume no-flux boundary conditions on the boundary of the domain, ∂Γ ( t ). In integralform, we write a conservation law for a small volume element V ( t ) = [ x ( t ) , x ( t ) + ∆x ( t )]dd t Z V ( t ) u ( x, t ) d x = Z V ( t ) [ −∇ · j + f ( u ) ] d x. Using the Reynolds transport theorem the left hand side is equivalent todd t Z V ( t ) u ( x, t ) d x = Z V ( t ) [ u t + ∇ · ( v u ) ] d x, where v is the velocity of the material point x ( t ). We then obtain the following reaction-transport-diffusion equation u t + ∇ · ( v u ) = D∆u + f ( u ) , in [ x ( t ) , x ( t )] . (14)Similarly, the boundary conditions on Γ ( t ) become − D ∇ u + v u = 0 , on ∂Γ ( t ) . We have two new terms1. An advection term, v · ∇ u , which corresponds to elemental volumes moving with the flow due to local growth.2. A dilution term, u ∇ · v , due to local volume changes. B Numerical Implementation
B.1 Mapping to a constant domain
It is most convenient to discretize and numerically solve the PDEs on a constant domain. Hence a mapping is used to transformthe growing/shrinking cell to the interval [ − , x, t ) (¯ x, ¯ t ) = (cid:16) x − x c ( t ) R ( t ) , t (cid:17) . Note that ¯ x is the location in the new domain [ − , x is located in the cell domain [ x ( t ) , x ( t )]. Carefully applyingthe chain-rule leads to u ¯ x = ∂u∂x ∂x∂ ¯ x = u x R ( t ) u ¯ t = ∂u∂t ∂t∂ ¯ t + ∂u∂x ∂x∂t ∂t∂ ¯ t = u t + u x (cid:0) ¯ xR ( t ) + x c ( t ) (cid:1) , where we have used x = ¯ xR ( t ) + x c ( t ) and t = ¯ t in the transformation. Substituting these results into Eqn. (2), we obtain u ¯ t − u ¯ x ¯ xR ( t ) + x c ( t ) R ( t ) + x c ( t ) + ¯ xR ( t ) R ( t ) u ¯ x + u R ( t ) R ( t ) = DR ( t ) + f ( u ) , in [ − , . Cancelling two terms on the LHS. and rearranging leads to u ¯ t = DR ( t ) u ¯ x ¯ x + f ( u ) − u (cid:18) ˙ R ( t ) R ( t ) (cid:19) , in [ − , . The total amount of GTPase on the interval [ − ,
1] becomes time dependent through G T = Z R ( t ) − R ( t ) ( G + G i ) d x = R ( t ) Z − ( G + G i ) d y. ell size, mechanical tension, and GTPase signaling in the Single Cell 23 B.2 Changing rest length is equivalent to imposing additional force
It is easy to see, from Eqn. (10) that adjusting the cell’s rest length, L is equivalent to assuming that additional forces act onthe cell: ( µ + 2 η ) d L d t = ( F − F ) − E ( L − L ) = − E ( L − ( L + ∆L )) , (15)where ∆L = ( F − F ) / E . In other words a change of ∆L is equivalent to an additional force of − E∆L .In (Zmurchok et al. 2018) and (Liu 2019) only the spatially uniform G model variants were considered. There, cell lengthwas modeled by (15), with ∆L ( G ) = ± b G m G mc + G m . Rac (+) is assumed to increase L and Rho (-) to decrease it, so the sign of the Hill function depends on the given GTPase.As shown here, the assumptions that GTPases work through modification of the cell’s effective rest length is mathematicallyequivalent to assumptions we make about GTPase-dependent load forces on the cell ends. B.3 Numerical methods
Each cell introduces a set of mechanical equations for its endpoints i.e. equation (9). Next, each cell contains K GTPases, thatis we include K sets of reaction-diffusion equations (1). The cell interior is divided into a cell-centered grid with uniform length h = 1 /N , where N is the number of grid cells per unit length. Here we use N = 2 . The second order derivative is discretizedusing a second order cell centered finite difference. For more detailed information on the numerical method, see (Gerisch 2001).This results in M total grid points per GTPase species, and K GTPase species. The resulting state vector is given by y = x x G , G , ... G K,M where G j,k denotes the concentration of the k -th GTPase at the j -th spatial grid point. The temporal evolution of the state ofthe cell can be written as: d y d t = f ( t, y ) , where f ( t, y ) is a function whose first two entries contain equation (9) and the final KM entries are the result of the spatialdiscretization of the GTPase PDEs.The resulting system of ordinary differential equations are integrated using the ROWMAP integrator (Weiner et al. 1997),an specialized integrator for large stiff ODEs (such systems are commonly the result of spatially discretizing PDEs). Here weuse the ROWMAP implementation provided by the authors of (Weiner et al. 1997) . The ROWMAP-integrator (written inFortran) is wrapped for easy use in a Python environment using f2py (f2py provides a connection between the Python andFortran languages)m and integrated into scipy ’s integrator class (a package providing several different techniques to integrateODEs). The spatial discretization (right hand side of ODE) is implemented using NumPy . The error tolerance of the integratoris set to v tol = 10 − . C Linear stability of spatially uniform GTPase model
In this section we consider the steady states and their linear stability of equation (13). \https://docs.scipy.org/doc/numpy-dev/f2py/ L eq , G eq ) is J ( L eq , G eq ) = − Eµ +2 η ± F max µ +2 η mG mc G m − eq ( G mc + G meq ) − G T L eq (cid:16) b + γ G neq G neq (cid:17) γnG n − eq ( G neq +1 ) (cid:16) G T L eq − G eq (cid:17) − (cid:16) b + I + γ G neq G neq (cid:17) . Model parameters are all positive, and L eq , G eq > J Rac = − + − sgn( d ) ! , J Rho = − −− sgn( d ) ! . where d := γnG n − eq (cid:0) G neq + 1 (cid:1) (cid:18) G T L eq − G eq (cid:19) − (cid:18) b + I + γ G neq G neq (cid:19) . Note that since (
L, G ) ∈ ∆ , where ∆ is the invariant region in Result 1, we have that d is the difference of two positive numbers.It is straightforward to conclude the following (see for instance (Murray 1989)). In the Rho case, the steady states are eitherstable nodes or saddles and no oscillations are possible. In the Rac case the steady states are stable when sgn( d ) <
0, andoscillations are possible when sgn( d ) > C.1 Sharp-switch limit
Using a sharp-switch approximation we can easily compute two steady states:1. When both nonlinear switches are off, i.e. G (cid:28)
1, the steady states are: L eq = L , G = bG T L ( b + I ) . Since the linearization of (13) is lower-triangular, the eigenvalues are readily found λ = − Eµ + 2 η , λ = − ( b + I ) . Since b, I, E, µ, η >
0, we conclude that this steady state is a stable node. Finally, since G eq (cid:28)
1, we obtain an existencecondition G T < L (cid:16) Ib (cid:17) . Therefore, the steady state of a large cell, with low total active GTPase can only exist if the total amount of GTPase issmall.2. In the case G c < G <
1, i.e. one of the nonlinear switches is turned on, the steady states are: L eq = L ± F max E , G = bG T L eq ( b + I ) . Since the linearization of (13) is lower-triangular as before, the eigenvalues are readily found λ = − Eµ + 2 η , λ = − ( b + γ + I ) . Since b, I, γ, E, µ, η >
0, we conclude that this steady state is a stable node. A condition for the existence of this steadystate is L ( b + γ + I ) b + γ < G T < G c L ( b + γ + I ) b + γ . Thus we must have some intermediate amount of GTPase.3. In the case 1 < G < G c , i.e. one of the nonlinear switches is turned on, the steady states are: L eq = L , G = ( b + γ ) G T L ( b + γ + I )ell size, mechanical tension, and GTPase signaling in the Single Cell 25Eigenvalues are readily found, as before, λ = − Eµ + 2 η , λ = − ( b + γ + I ) . From b, I, γ, E, µ, η > L ( b + γ + I ) b + γ < G T < G c L ( b + γ + I ) b + γ , L G T − L I − b < γ < L G c G T − L G c I − b. For these ranges to be non-empty we require that G T > L G c . Thus we must have some intermediate amount of GTPase.4. In the case, when both nonlinear switches are turned on, i.e. G (cid:29) G c , the steady states are L eq = L ± F M E , G = ( b + γ ) G T L eq ( b + γ + I ) . Then the eigenvalues are λ = − Eµ + 2 η , λ = − ( b + γ + I ) . Since b, I, γ, E, µ, η >
0, this steady state is a stable node. The condition for existence of this steady state is given by G T > G c L eq (cid:16) Ib + γ (cid:17) , or γ > L b G c G T − L b G c I − b. Thus only cells with a sufficiently amount of GTPase can be large or small.In the bifurcation diagram of Fig. 8, the regions of existence of G , G or G are identified. D The case of a deforming cell with conserved volume
Here we consider the case that the volume of a cell is roughly constant, while the cell deforms and its membrane size changed.This idea suggests a slight modification to our model. The modified mass conservation reads G T = G i V + βGL, where V is the constant volume of the cell, and β is a constant of proportionality with units of [length] . We can solve for G i to get a relation analogous to (6): G i = G T V − βGLV . Notice that after substituting this into the well-mixed system, we can set V = 1 by rescaling b and γ . Therefore, we getd G d t = (cid:16) b + γ G n G n (cid:17) ( G T − βGL ) − IG − G (cid:18) ˙ LL (cid:19) . (16)This, together with (13a) and (13b), forms our new system.Much of the analysis of this system is similar to the earlier case. The changes to the nullclines in the G − L phase planecorrespond to a vertical shift in the sharp switch limit, which can be compensated by adjusting the parameters. This suggeststhat the behaviors of the new system is similar to the original (13). The equilibrium branches and the conditions for theirexistence in the sharp switch limit are:1. For G (cid:28) L eq = L , G eq = bG T bβL + I exists provided G T < βL + Ib .
2. For G c < G < L eq = L ± F max E , G eq = bG T bβL eq + I exists provided G c (cid:16) βL eq + Ib (cid:17) < G T < βL eq + Ib .
3. For 1 < G < G c : L eq = L , G eq = ( b + γ ) T ( b + γ ) βL + I exists provided βL + Ib + γ < G T < G c (cid:16) βL + Ib + γ (cid:17) . Fig. 11
Bifurcation diagrams as in Fig. 8, but for the constant volume spatially uniform GTPase model, (16). Top: thecontractile (Rho) case, with parameters k = 1 , I = 5 , L = 3 . , Eµ +2 η = 3 , F max µ +2 η = 10 , G c = 4 , β = 1 , G T = 5 , m = n = 10.Bottom: the protrusive (Rac) case, with parameters k = 1 , I = 10 , L = 2 . , Eµ +2 η = 8 , F max µ +2 η = 55 , G c = 0 . , β = 1 , G T =10 , n = 20 , m = 10. While the Rac case is still tristable, a fold point that previously connected the middle and high branch (inthe orignal model) is now missing. These branches persists as γ → ∞ whereas the lower branch exist only for low γ . In the Raccase, the Hopf curve in the two-parameter bifurcation digram (lower right) differs significantly from the corresponding panel inFig. 8, which yield parameter regimes where the limit cycle exists for γ → G (cid:29) G c : L eq = L ± F max E , G eq = ( b + γ ) G T ( b + γ ) βL eq + I exists provided G T > G c (cid:16) βL eq + Ib + γ (cid:17) . In the Rho (contraction) case, tri-stability still occurs. In this model, we can more readily find parameter sets where themiddle (1 < G < G c ) and high branch persist as γ → ∞ , whereas the low branch terminates at a fold point (Fig. 11). (This isalso possible in the original model.) This suggests that even with a very large positive feedback in GTPase activity, and whenGTPase activity is high enough to trigger that feedback, the cell can still stay at a large size. In contrast, in the original model,the regime where only the middle and high branch coexist is usually bounded above by some maximum γ , as seen in the regimelabelled “2,3" in Fig. 8 (top right).In the Rac (expansion) case, periodic orbits are still possible. The regime in which limit cycles exist is much wider that inthe original model. We can identify parameter sets where a limit cycle exists even at γ = 0, which hints that in this model, apositive feedback in GTPase activity is not neccessary for an oscillating cell.Overall, the model variant with constant cell volume has the same regimes as the original model. However, the location andsize of these regimes has shifted.ell size, mechanical tension, and GTPase signaling in the Single Cell 27The Jacobian of the ODE system at a steady state ( L eq , G eq ), is J ( L eq , G eq ) = − Eµ +2 η ± F max µ +2 η mG mc G m − eq ( G mc + G meq ) − Eµ +2 η L eq γnG T G n − (1+ G n ) − βL (cid:0) b + γ G n G n + G (cid:1) − I . Model parameters are all positive, and L eq , G eq > J Rac = − + − sgn( d ) ! , J Rho = − −− sgn( d ) ! . where d := γnG T G n − (1 + G n ) − βL (cid:16) b + γ G n G n + G (cid:17) − I. The sign patterns are the same as in the original model, and thus the same conclusions apply.