Spots, strips, and spiral waves in models for static and motile cells
JJournal of Mathematical Biology manuscript No. (will be inserted by the editor)
Spots, strips, and spiral waves in models for static and motile cells
GTPase patterns in cells
Yue Liu · Elisabeth G. Rens · Leah Edelstein-Keshet
Received: date / Accepted: date
Abstract
The polarization and motility of eukaryotic cells depends on assembly and contraction of theactin cytoskeleton and its regulation by proteins called GTPases. The activity of GTPases causes assemblyof filamentous actin (by GTPases Cdc42, Rac), resulting in protrusion of the cell edge. Mathematical modelsfor GTPase dynamics address the spontaneous formation of patterns and nonuniform spatial distributionsof such proteins in the cell. Here we revisit the wave-pinning model for GTPase-induced cell polarization,together with a number of extensions proposed in the literature. These include introduction of sources andsinks of active and inactive GTPase (by the group of A. Champneys), and negative feedback from F-actinto GTPase activity. We discuss these extensions singly and in combination, in 1D, and 2D static domains.We then show how the patterns that form (spots, waves, and spirals) interact with cell boundaries to createa variety of interesting and dynamic cell shapes and motion.
Keywords
Pattern formation, intracellular signaling, GTPase, wave-pinning, local perturbation analysis,static and moving boundary computation
The dynamics of the actin cytoskeleton determines internal cell structure, cell shape, and cell motility. Byaccumulating at a cell edge, filamentous actin (F-actin) produces outwards protrusion. Actin assembly isregulated by signaling networks. Central in those networks are the small GTPases, Rac, Cdc42, and Rho.Rac promotes assembly of F-actin, whereas Rho activates myosin motors. The interactions of Rac, Rho,Cdc42, and other molecular players has been modeled in previous work (Mori et al 2008; Verschueren andChampneys 2017; Holmes et al 2012a; Holmes and Edelstein-Keshet 2016; Zmurchok et al 2018; Walther
Y. LiuDepartment of Mathematics, University of British Columbia, Vancouver V6T 1Z2, BC, CanadaE-mail: [email protected]
Present address:
Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UKE. G. RensDepartment 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 ] S e p Yue Liu et al. et al 2012; Edelstein-Keshet et al 2013; Jilkine and Edelstein-Keshet 2011; Otsuji et al 2007) both in 1D and2D. These studies made different modelling decisions and ranged from simple (Mori et al 2008) to detailed(Mar´ee et al 2008). It is challenging to determine parameter sensitivity and map out regimes of behaviorof the more detailed models. This motivates studying minimal models that showcase the possible realms ofpredicted behavior.It was shown previously that the biology of GTPases permits a single member of this family to sponta-neously polarize (i.e. form spatial regions of high vs low activity). This idea was the basis of the wave-pinningmodel (Mori et al 2008, 2011), and depends on the large difference in diffusion of the active (slow) and in-active(fast) forms of a GTPase.Several models have been examined mathematically to describe how a single GTPase coupled to othereffectors or influences could results in spatio-temporal patterns. These include a GTPase with sources andsinks (Verschueren and Champneys 2017), with feedback from F-actin (Holmes et al 2012a; Mata et al 2013),with mechanical tension (Zmurchok et al 2018) and with effects of changing cell size (Buttensch¨on et al 2019).Many of these were explored in reaction-diffusion (RD) equations within a 1D static single cell domain orwith spatially uniform distribution in each of many cells (Zmurchok et al 2018). Some of the behaviors foundin such models include, traveling waves, pulses, or oscillating fronts (Holmes et al 2012a; Mata et al 2013),or localized peaks and “solitons” (Verschueren and Champneys 2017).Here we have two main purposes: (1) to explore what happens when two distinct minimal models arecoupled, and whether this leads to new behavior, (2) to study these systems in 2D domains to determinewhether they produce spots or stripes, and (3) to simulate the same models on a deforming 2D domaindepicting the shape and motility of a cell.Biological motivation for this work comes from several sources. (A) Waves of actin are observed in anumber of experimental systems (Inagaki and Katsuno 2017). In some of these, such waves are seen to causecell edge to cyclically protrude outwards (as the waves impinge on the cell edges). We wondered whether amodel for Rac interacting with F-actin could mimic this kind of behavior. (B) The GTPase model generalizedby the group of Alan Champneys in (Verschueren and Champneys 2017) converts the polarizing cell behaviourinto multiple coexisting peaks. We wondered how such peaks would interact with cell boundaries, and, inparticular, whether they would be associated with smaller protrusions such as filopodia. (C) In some cells,notably the embryos of C. elegans, localized Rho-associated actin clusters are seen to “blink” (oscillatetemporally while maintaining a fixed location) (Robin et al 2016). We asked whether the combined F-actin-Rho model with localized sources could account for such behavior.We first briefly review the three classes of minimal models, show results for the combined model, andthen demonstrate the novel 2D behaviors that are observed once these models are simulated in the deforming2D cell. pots, strips, and spiral waves in models for static and motile cells 3
Fig. 1: Schematic diagram of the models. The original wave-pinning model consists of GTPase (circles) inthe active, (membrane-bound) form, u and inactive form v , with positive feedback (curved grey arrow) from u to its own activation (upwards white arrow). The F-actin extension model (Holmes et al 2012a) includesGTPase activation of F-actin assembly and GTPase inactivation by F-actin (dashed arrow). The source-sink(nonconservative) extension by (Verschueren and Champneys 2017) includes removal of active GTPase andsynthesis of inactive GTPase so that the total amount is no longer conserved. Our model is a system of reaction-diffusion partial differential equations (PDEs) based on the wave pinningmodel first proposed by (Mori et al 2008). The model is extended with a source and sink terms following(Verschueren and Champneys 2017), and feedback from actin, proposed by (Holmes et al 2012a).2.1 Model equationsThe dimensionless form of the model combining both extensions can be written as: ∂u∂t = δ ∇ u + f ( u, v, F ) − cθu, (1a) ∂v∂t = ∇ v − f ( u, v, F ) + cα, (1b) ∂F∂t = (cid:15) ( k n u − k s F ) , (1c) f ( u, v, F ) = A ( u ) v − (cid:18) η + s F F (cid:19) u, A ( u ) = k + γ u n u n , (1d) ∂u∂ n (cid:12)(cid:12)(cid:12)(cid:12) ∂Ω = 0 , ∂v∂ n (cid:12)(cid:12)(cid:12)(cid:12) ∂Ω = 0 , x ∈ Ω, t ≥ . Here u ( x, t ) and v ( x, t ) represent the active and inactive GTPase, respectively. F ( x, t ) represents filamentousactin (F-actin). δ (cid:28) f ( u, v, F ) describes the net rate of GTPase activation, with A ( u )representing activation rate. Parameters k, γ, η, s are the basal activation rate, self-feedback activation, basalinactivation and actin-feedback inactivation rates, respectively. Neumann boundary conditions are used torepresent the fact that GTPases and F-actin do not leak out of the cell edges. Yue Liu et al.
Setting c = s = 0 reduces the system to the original wave pinning (WP) model (Mori et al 2008), whichconserves the total u + v inside the domain. The model has been analyzed in detail elsewhere (Mori et al 2011),but we briefly mention its key property: under specific parameter settings, the WP model sustains wavesthat decelerate and stall in the domain, leading to a stable spatially heterogeneous steady state distributionof u (the “pinned wave”).When c = 1 , s = 0, the system corresponds to the non-conservative (NC) model of (Verschueren andChampneys 2017). When c = 0 , s >
0, we have the actin feedback (AF) model of (Holmes et al 2012a).While each of the above models has been studied previously, here we will also be concerned with theirunion, i.e. the so-called “combined model” (CM) with c = 1 , s >
0. The four models of interest are then(I) WP, (II) NC, (III) AF, and (IV) CM. These four models all have very distinct characteristic behaviors.We will consider these models in several settings (A) a 1D spatial domain, as previously described in theliterature, (B) a static 2D spatial domain where we can distinguish between spots and stripes, and finally(C) a deforming domain whose boundary dynamics is coupled to the evolving solution u (or F ) of the PDE.2.2 GeometryIn many previous papers, simulations were restricted to 1D (Holmes et al 2012a; Mata et al 2013), but avariety of actin wave models exist in more detailed geometries, including 2D (Doubrovinski and Kruse 2011)and 3D (Bretschneider et al 2009). Here, for for simplicity in the 2D static domain case we consider a unitsquare. For the deforming domain, we use the Cellular Potts Model (CPM) to simulate a dynamic 2D cell.The methods and results are introduced in Sec. 5.3.For the ease of analysis and identification of distinct patterns, we first discuss and examine results ina 1D spatial version of the models. We can interpret this 1D geometry in one of two classic ways: (1) asa cross-section along the diameter of a cell. This cross-section neglects any variation in the cell thicknessand includes both the intracellular volume (cytosol) and the top and bottom membranes at every point.Neumann (no flux) boundary conditions are used for the endpoints of the interval. (2) Alternatively, anothercommon assumption is a 1D cell perimeter. In this case, the region considered is close to the cell membrane,with periodic boundary conditions. Here we adhere to the first approach. The case of 1D dynamic cell sizeis considered in (Buttensch¨on et al 2019). We briefly describe methods used to analyse the models. We use local perturbation analysis (LPA) to studythe bifurcation behavior of each model, and compare with results from Turing Linear stability analysis. Afull description of these methods is found in the MSc thesis by one of us (YL) (Liu 2019).3.1 Local perturbation analysisLocal perturbation analysis is a method for examining the evolution of a localized perturbation to a homoge-neous steady state (HSS) for a fast-slow diffusion-reaction system. It provides a way to systematically detectcertain forms of nonlinear instabilities that are not detectable by the more traditional Turing analysis. LPAwas first developed by AFM Mar´ee and V Grieneisen (Cardiff University) (Grieneisen 2009), and has beenused in (Edelstein-Keshet et al 2013; Holmes and Edelstein-Keshet 2016; Holmes et al 2012a; Mata et al2013) and elsewhere to analyze wave pinning and related models.The basic idea of LPA is to take the limit where the slow diffusion coefficients goes to 0 and the fastdiffusion coefficients goes to infinity. We then consider an initial condition where the system is at HSS witha localized perturbation in the form of a spike of infinitesimal width but finite height. The behavior of the pots, strips, and spiral waves in models for static and motile cells 5
PDE can then be captured with an ODE system with “global variables” representing the levels of the PDEvariables away from the spike, and “local variables” for the slow PDE variables at the spike. For example,using subscript L to denote local variables, the LPA system for our combined model (1) is: ∂u∂t = f ( u, v, F ) − cθu, (2a) ∂v∂t = − f ( u, v, F ) + cα, (2b) ∂F∂t = (cid:15) ( k n u − k s F ) , (2c) ∂u L ∂t = f ( u L , v, F L ) − cθu L , (2d) ∂F L ∂t = (cid:15) ( k n u L − k s F L ) . (2e)In the cases where c = 0 or s = 0, we will use mass conservation to remove irrelevant equations and eliminatedegeneracy. This allow us to easily produce bifurcation diagrams using AUTO (Doedel 1981) and delineateparameter regimes. Notice that the LPA system (2) contains the well-mixed system (i.e. the system withoutlocal variables), so any features (branches and bifurcations) of the well-mixed system will also be presentin the LPA system. Hence we can obtain any information that can be gained by analyzing the well-mixedsystem through LPA.3.2 Bifurcation analysisWe refer to branches of equilibria and periodic solutions in the LPA system that are also present in the well-mixed model as “global” branches, as they correspond to solutions in which the local variables are equal tothe global variables and the spike disappears, i.e. a homogeneous solution. The others branches are referredto as “local” branches; they correspond to some kind of pattern.We classify the parameter regimes into three categories: (a) stable, where only global branches are stable.In this regime, no pattern can arise from localized perturbation; (b) polarizable, where stable global and localbranches coexist. In this regime, patterns can form only if the perturbation is sufficiently strong. Finally, (c)unstable, where all global branches are unstable. In this regime even infinitesimal perturbations can lead topattern formation. In Appendix 8, we show that this is equivalent to the classical Turing regime.The sets of parameters for each model are listed in Table 4. For the cases where the total amount ofGTPase is conserved, we define the total mass of GTPase in the cell, w = (cid:90) Ω ( u + v ) dx, as an additional constant parameter. This allow us to eliminate v from the equations by writing it in termsof u and w .All bifurcation diagrams follow AUTO’s conventions. On one-parameter diagrams, red/black curves indi-cate positions of stable/unstable equilibria respectively, while green/blue indicate the range of stable/unstablelimit cycles. On two-parameter diagrams, red/light blue/dark blue curves trace the position of limit points(fold points)/branch points (transcritical points)/Hopf points, respectively. We next apply the methods to compare the behaviors of the four models of interest.
Yue Liu et al. γ (the magnitude of the only nonlinear term) and w (total concentration)as primary parameters of interest. Extending the earlier study (Holmes and Edelstein-Keshet 2016), we alsotrace a branch of transcritical bifurcation in two-parameter continuation in Fig. 3. This allows us to identifyseveral new regimes. We verify that LPA predictions in each regime are indeed correct with simulations ofthe full PDEs. Regime Classification DescriptionI Stable One stable GB, no LBII Polarizable One stable GB, one stable LBlocated above the GBIII Polarizable One stable GB, three stable LBslocated on both sides of the GBIV Polarizable Two stable GBs, three stable LBs: one above bothGBs, one in between, and one below both GBsV Polarizable One stable GB, one stable LBlocated below the GBVI Unstable The only GB is unstable, two stable LBslocated on both sides of the GBVII Unstable Three GBs, all unstable, four stable LBslocated on both sides of the GB
Table 1: Summary of the wave pinning (WP) regimes identified in Fig. 2 and 3. GB: global branch; LB: localbranch. Stable: all stable branches are global branches. Polarizable: there exist both stable global and localbranches. Unstable: all global branches are unstable, so some local branches have to be stable.4.2 Non-conservative (NC) modelIn addition the main bifurcation parameter γ , we also take c , the parameter that controls the magnitude ofthe source/sink terms. This model possess a unique global equilibrium: u ∗ = αθ , v ∗ = cα + ηu ∗ A ( u ∗ ) = cα + ηu ∗ k + γ u n ∗ u n ∗ . Any local branches u L ∗ must satisfy f ( u L ∗ , v ∗ ) = 0. After expanding and some manipulations, we obtain A ( u L ∗ ) A ( u ∗ ) = u L ∗ u ∗ . (3)Since neither A ( u ) nor u ∗ involve c and η , we conclude that the local branches are independent of theseparameters. Furthermore, for γ (cid:28) k , the LHS of (3) ≈ u L ∗ = u ∗ , which means that there is no localbranch for small γ .We will show that for γ (cid:29) k , there are always a high and low local branches. The low branch is u L ≈ γ → ∞ and u L = 0, both sides of (3) evaluate to 0. With a bit of further manipulation, we get(in the limit γ → ∞ ): h ( u L ∗ ) = h ( u ∗ ) , where h ( u ) = u n − u n . pots, strips, and spiral waves in models for static and motile cells 7(a) WM, w = 2 (b) LPA, w = 2(c) WM, w = 3 . w = 3 . w = 4 (f) LPA, w = 4 Fig. 2: Bifurcation diagrams of the well-mixed (WM) and LPA wave pinning system with respect to therate of activation parameter γ . Other parameters as in Table 4(WP) except w . The purple lines are locatedat bifurcation points separating the distinct regimes. Note that the “global branches” (curves in the WMdiagrams) also appear in LPA, though their stability can be different in LPA over certain intervals. Yue Liu et al.(a) WM (b) LPA(c) LPA
Fig. 3: Two-parameter bifurcation plots of the wave pinning (WP) model with respect to parameters w, γ . (a)Well-Mixed (WM) and (b,c) LPA system. Other parameters as in Table 4(WP). Each curve in these diagramstraces the location of a bifurcation point shown in Fig. 2, and forms the boundary of a parameter regime.The one-parameter bifurcation diagrams in Fig. 2 correspond to vertical cross-sections of the diagrams here.The LPA regimes I - VII match with the regimes in Fig. 2(b,d,f). See summary in Table 1. (c) A zoom intothe cusps in (b). (Compare (b) to LPA Fig. 3(a) of (Holmes and Edelstein-Keshet 2016) for the same modelwith different parameter values: our figures agree on the (red) fold curves but ours includes an additionaltranscritical curve (light blue) separating several distinct regimes.) pots, strips, and spiral waves in models for static and motile cells 9
The function h ( u ) satisfies h (0) = 0 = h ( u → ∞ ), and it has a single peak at u p ≥ n ≥ α < θ , that is u ∗ <
1, there exists a point u L ∗ > u p > h ( u L ∗ ) = h ( u ∗ ), which corresponds to the high local branch.In the bifurcation diagrams in Fig. 4, we use parameters from Table 4 (CM2) but with η = 5. (Theseparameters yield visually optimized bifurcation diagrams whose regimes are neither too wide nor too narrow;the same regimes are present for parameters from Table 4(NC) used for PDE simulations, but the resultingbifurcation diagram is harder to read.) Fig. 4 identifies four distinct regimes in the LPA system whoseinterpretation is as in the previous section. The location of the branches agrees with earlier analysis. Theregimes are summarized in Table 2.In Regime I, no pattern forms, as expected. In Regime II a small perturbation decays, but a sufficientlylarge perturbation will persist. In the full PDEs, such perturbation leads to the soliton solution shown inFig. 10(c,d). Both Regime III and IV are unstable, and any perturbation leads to a Turing-type patternconsisting of a series of evenly spaced, static spikes in the full PDE, as in Fig. 10(a,b). The limit cycles inRegime IV suggests that the spikes might oscillate, but this, in fact, does not occur: we found that the PDEbehavior is qualitatively indistinguishable in Regime III and IV. This suggests that the Hopf bifurcationsin the LPA diagram may not necessarily correspond to actual bifurcations for the full PDE, pointing to alimitation of LPA. Regime Classification DescriptionI Stable One stable GB, no LBII Polarizable One stable GB, one stable LBlocated above the GBIII Unstable The only GB is unstable, two LBslocated on both sides of the GBIV Unstable One GB, two LBs all unstable,each enclosed by a periodic orbit
Table 2: Summary of the non-conservative (NC) model regimes identified in Fig. 4. Abbreviations as inTable. 1.4.3 Actin feedback (AF) modelWe use mass conservation to eliminate v from the LPA system as before. The strength of actin feedback s and the basal rate of activation k were our bifurcation parameters. LPA for this model was previouslydiscussed in (Holmes et al 2012a; Mata et al 2013), but here we traced more bifurcations in greater detail.The results are shown in Fig. 5 and 6. We only distinguish between the regimes separated by fold andtranscritical curves and omit the Hopf curves, as explained below. We also ignore some very narrow regimes,to concentrate on six major regimes as summarized in Table. 3.One interesting characteristic of these diagrams is the presence of unstable periodic orbits that emerge assubcritical Hopf bifurcations and exist for very narrow parameter ranges. The unstable cycle enlarges untilit collides with a saddle point, turning into a homoclinic orbit to the saddle, and then disappearing. Thisis known as saddle-loop bifurcation, or homoclinic bifurcation (see (Kuznetsov 2004, Ch.6.2)). Parameterregimes where the periodic solutions exist are very narrow. Hence, while Hopf bifurcations occur, they areunlikely to be playing a major role in the biological application of this model.We can compare our results to those of (Holmes et al 2012a) (Fig. 5, a LPA diagram in k − s planecontaining only one of the Hopf curves). The Hopf curve in (Holmes et al 2012a) corresponds to the dark bluecurve on our diagram, which traces the pair of Hopf points on the global branch in Fig. 5(d)). Furthermore, Fig. 4: Bifurcation diagrams for the non-conservative (NC) model, with parameter values from Table 4(CM2)except η = 5. (a) WM, (b,c) LPA, using bifurcation parameters (a,b) γ , with c = 1, (c) c and γ . Athin polarizable regime II is sandwiched between the stable I and Turing III regimes. The triplet of Hopfbifurcations (not present in WP) does not show up as new behavior in the full PDE simulations. pots, strips, and spiral waves in models for static and motile cells 11 our diagram (Fig. 6) traces the fold (red) and transcritical (light blue) bifurcation points and hence identifiesa larger number of distinct regimes.Interpreting the LPA diagrams (as in the WP model), we can conclude that a stable local branch in LPAcorresponds to a regime of pattern formation in the PDE. Unlike WP, there are multiple possible patternsin this AF model. LPA cannot accurately predict the type of pattern. In particula, the consequence of thesubcritical Hopf bifurcations to the full PDE is unclear, possibly suggesting some kind of (quasi-)periodicbehavior that we did not fully characterize. In (Holmes et al 2012a; Mata et al 2013), a parameter scan ofthe PDE system was included with the LPA diagrams. As previously noted, PDE regimes are not exactlyaligned with LPA regimes since δ (cid:54) = 0 in the full PDEs.In summary, in our hands, LPA worked well in identifying no-pattern and WP regimes, but was lessuseful for predicting the emergence of more complex patterns. Many of those patterns involve interactingwaves, which suggests that they are non-linear, non-local phenomena, explaining why LPA cannot accountfor them. Regime Classification DescriptionI Stable One stable GB, no LBII Polarizable One stable GB, one stable LBlocated above the GBIII Unstable The only GB is unstable, two stable LBslocated on both sides of the GBIV Unstable The only GB is unstable, one stable LBlocated above the GBV Polarizable Two stable GBs, three stable LBs: one above bothGBs, one in between, and one below both GBsVI Polarizable One stable GB, three stable LBslocated on both sides of the GB
Table 3: Summary of the actin feedback (AF) model regimes identified in Fig. 5. For abbreviations seecaption of Table. 1.4.4 Combined model (CM)The LPA diagrams for the combined model are very complex, and mostly beyond the scope of interpretation(see Appendix 7.) This is unsurprising given the complex behavior exhibited by the PDE. The bifurcationdiagram shown in Fig. 20, which uses parameter values from Table 4(CM2), contains many limit cyclebifurcations, such as torus and period-doubling. One thing the diagram can provide is the minimum valueof s required for any non-static patterns (corresponding to the first triplet of Hopf bifurcation in Fig. 20).With s below this value, the system behaviour is the same as that of the s = 0 case, which reduces back tothe non-conservative (NC) model.4.5 Comparison with linear (Turing) stability analysisLinear stability analysis (LSA) was previously applied by (Mori et al 2008) and (Verschueren and Champneys2017) for the wave pinning and non-conservative models respectively. The relative merits of LSA and LPAhave been described in (Mata et al 2013; Holmes 2014) and we briefly summarize some of these in theAppendix.LPA is only valid in the limit of δ →
0. In this limit, LPA contains the Turing stability properties: abranch that is LPA-unstable is also Turing-unstable. See Fig. 21, where we show how the LPA regimes from k = 1 . k = 1 . k = 6 (d) LPA, k = 6(e) LPA (f) LPA Fig. 5: Bifurcation diagrams of the actin feedback (AF) model with respect to parameter s (a-d) and withrespect to k, s in (e,f). (In (e), the Hopf curves are omitted for clarity of the diagram. They are then includedin (f).) The narrow regimes are not labelled. The nearly vertical blue curves indicate unstable periodic orbits. pots, strips, and spiral waves in models for static and motile cells 13 Fig. 6: Same as Fig. 5(f) with the Hopf curves included, and with an indication of patterns in several regimes.A few Hopf curves lie very close to one of the other curves for most of their length, creating some very narrowregimes. The simulation results from Fig. 8 are identified with their corresponding regions on the parameterplane.Fig. 2(b) line up with Turing regimes. LPA can detect instabilities that require a perturbation of sufficientmagnitude (the polarizable regimes), which cannot be detected by Turing analysis. This means LPA canpotentially find more types of pattern.LPA does not predict details of the pattern. We saw this most evidently in the actin feedback (AF)model, where many possible patterns and a large number of parameter regimes exist. Turing analysis predictspattern initiation, but often fails to specify the final pattern that depends on nonlinear interactions. We givean example of this type for the NC model in Fig. 22. We also indicate how the “minimal patch size” ideafrom (Painter and Hillen 2011) can be used to help predict the final pattern using LSA.
We simulated the model for a static cell in 1D (0 ≤ x ≤
1) and 2D (0 ≤ x, y ≤ for numerical simulations are summarized in Table 4. The selection of values for most of these parametersis based on (Holmes et al 2012a), with α, θ coming from (Verschueren and Champneys 2017), and somemodifications guided by LPA and Turing analysis. In contrast to (Holmes et al 2012a) we use a much largerdomain size L , corresponding to a larger cell and allowing for more complex patterns to develop. Parameter Meaning WP NC AF CM2 δ Diffusion coefficient ratio 0 . L Domain length 1 10 k Basal activation rate 1 . L L − L L γ Nonlinear activation rate 30 L n Hill coefficient 3 2 η Inactivation rate 15 L L L . L c NC terms on/off 0 1 0 1 α Source strength 1 . L θ Sink strength 4 . L . L s Actin feedback strength 0 0 − L (cid:15) Actin reaction rate 0 . k n Actin activation rate 24 L k s Actin inactivation rate 7 . L Table 4: The parameters in the combined model, their meanings and values for various simulations. WP:wave pinning; NC: non-conservative extension; AF: actin feedback extension; CM2: one of the parametersets used for the combined model (CM). All parameters (except L ) are scaled to be non-dimensional.5.1 Simulations in a fixed 1D domainWhile 1D simulations for the WP, NC and AF appear in previous works (Mori et al 2008; Verschueren andChampneys 2017; Holmes et al 2012a; Mata et al 2013), we present them here as comparison to the combinedmodel and the 2D case. Results are shown as kymographs, with time on the horizontal axis and is space onthe vertical axis. Color indicates the levels of u and v and/or F (if s > u can varygreatly across the domain while v becomes nearly uniform, as expected given its much faster rate of diffusion.Fig. 8 shows the results for the actin feedback (AF) model. The default initial conditions are u = 0 except u = 4 for 0 ≤ x ≤ . v = 2 . F = 0. We are not initializing near a stable HSS because doing so usuallydoes not result in patterning. The patterns observed are quite sensitive to initial conditions. In addition tosimple wave pinning observed at low s (not shown), the system displays four qualitatively different behaviors:(1) wave pinning with oscillating boundary (WPO), where polarization occurs as in wave pinning, but withan oscillating front position; (2) reflecting pulse (RW), where a single pulse traverses the domain at constantvelocity and gets reflected back at the boundary; (3) a single pulse (SP) that is absorbed at a boundary,before the system returns to HSS; (4) a wave train (WT), that originates either at a boundary or in theinterior of the domain, propagates with constant velocity and gets absorbed at a boundary.In general, the spatial profile of F lags behind u , as expected, since it is a slow variable depending on u .The pattern in v is usually opposite that of u , i.e, v is high where u is low, and vice versa. Moreover, thegradient of v tends to be much shallower than u due to the faster diffusion of v . pots, strips, and spiral waves in models for static and motile cells 15 Some other more complex patterns are shown in Fig. 9. These share some characteristics with the simplerpatterns. The patterns shown in (a-c) are similar to (WPO), but the domain is divided into five regionsinstead of two, with an initial transient reminiscent of (WT). The patterns in (d-f) can be seen as a group offour reflecting pulses (similar to RW) rather than one. Compared to (Holmes et al 2012a), we find a richerrange of patterns using a similar parameter set (with different scaling). The main difference is that the largerdomain used here, L = 10, allows more space for pattern to develop. (In (Holmes et al 2012a), L = 1, sopatterns are more confined and boundary effects are prominent.)Fig. 10 shows two typical patterns in the NC model: a static, Turing-type pattern consisting of a seriesof evenly spaced spikes, and a single spike “soliton” pattern. The final profiles of these two patterns areshown in Fig. 11(a,b). The domain length, L , must be large enough to support such patterns. If L is toosmall to support a full period of the pattern, the result would be simple polarization similar to wave pinning(Fig. 11(c)). Using a higher rate of inactivation η , or a smaller diffusion ratio δ can result in spikes that splitinto two, as shown in Fig. 11(d).For the combined model, we use parameters from Table 4(CM2), mostly similar to the NC case. In Fig. 12,we show the effect of increasing s (strength of actin feedback) on system behavior. With s low enough, thesystem behavior resembles the s = 0 case of a static, spatially periodic pattern, as in the NC case. Forincreasing s , the peaks begin to move with constant velocity by themselves, repelling one another when tooclose. For moderate values of s , the peak repulsion is strong enough that peaks reverse their direction ofmotion if on a collision course (Fig. 12(a)). At higher s , they collide (Fig. 12(b)). At even higher s , we observea localized standing wave pattern that oscillates rapidly in Fig. 12(c), and even more prominently in (d).5.2 Simulations in a fixed 2D domainIn two spatial dimensions, we use the same parameters as in 1D. For the WP model, we start at HSS andperturb one corner of the domain. The pattern we observe (Fig. 13(a,b)) is a direct analogue to the 1Dcase (compare to Fig. 7(a,b)): u initially spreads out from the corner as a 2D travelling wave, and that iseventually pinned along a front determined by the initial conditions.We use a similar initial condition for the NC model. Based on 1D simulations, we expect evenly-spacedstripes to form around the corner as concentric rings, as happens initially (Fig. 13(c)). However, these ringsquickly break up into spots (Fig. 13(d)). The spots spread out, and then settle into a steady state. Wehave not found any parameter sets for stable ring patterns. The patterns are insensitive to the shape of thedomain. Simulations on circular, rectangular and other domains with simple shapes produced patterns withthe same qualitative characteristics (not shown).For CM, we initialize the system at HSS and perturb with noise. Fig. 15 shows the simulation results.With a low s , the pattern is indistinguishable from the static spots under the non-conservative model. As s increases, the spots become mobile and repel each other as in the 1D case. In 2D, as s is increased further,the spots transitions to spiral waves.We also arrived at the AF model by initializing the CM model at HSS plus global noise and c = 1. After apattern starts to form, we gradually decreased c to 0 to arrive at the AF model. (In our hands, this producedmore robust results, with patterns that persisted.) Fig. 14 shows a few snapshot of the simulations. For low s , the pattern resembles slowly drifting and deforming blobs. As s increases, the pattern transitions intospiral waves with decreasing width.5.3 Simulations in a 2D deforming domainAs a final set of numerical experiments, we simulate the models in an evolving 2D domain. The boundariesdeform in response to the chemical levels close to the boundary. We use the Cellular Potts Model (CPM) forthese examples. Fig. 7: Simulation of the wave pinning model (WP), with parameters from Table 4(WP). Initial condition: v = 1, u = 0 .
102 with perturbation u = 6 for (a,b) 0 ≤ x ≤ .
1; (c,d) 0 . ≤ x ≤ .
5; (e,f) random noise, u = 0 . · (cid:15) ( x ). Note that not all initial conditions result in wave pinning: a small perturbation from theHSS will simply decay and no pattern forms. The behaviors shown in (a,b,e,f) correspond to solutions shownin Fig. 2 of (Mori et al 2008). pots, strips, and spiral waves in models for static and motile cells 17(a) u , k = 1 . , s = 18 (b) v , k = 1 . , s = 18 (c) F , k = 1 . , s = 18(d) u , k = 1 . , s = 27 (e) v , k = 1 . , s = 27 (f) F , k = 1 . , s = 27(g) u , k = 1 . , s = 36 (h) v , k = 1 . , s = 36 (i) F , k = 1 . , s = 36(j) u , k = 6 , s = 30 (k) v , k = 6 , s = 30 (l) F , k = 6 , s = 30 Fig. 8: Simulations of the actin feedback model (AF) with parameters from Table 4(AF) ( s, k as indicated onlabels), and default initial conditions. Each row corresponds to one parameters set, showing u, v, F (left toright). We observe four behaviors by varying k and s : (a-c) Wave pinning with oscillating front (WPO); (d-f)Reflecting waves (RW); (g-i) Single pulse absorbed at boundary (SP); (j-l) Persistent wave trains (WT). Weused a larger domain length than (Holmes et al 2012a), leading to a richer set of patterns. Fig. 9: Exotic patterns observed in the actin feedback (AF) model. Parameters as in Table 4(AF) but varying k, s . (a-c) k = 5 , s = 10, default initial conditions. The pattern resembles WPO but with several subregions;(d-f) k = 5 , s = 30, default initial conditions with excitation region 0 ≤ x ≤ .
1. The resulting pattern issimilar to RW but with a group of four pulses traversing the domain.Full details of the CPM can be found elsewhere (Mar´ee et al 2007) and are briefly summarized in theAppendix. The essential feature of the CPM is its ability to track an evolving shape such as morphologyof a motile biological cell (Scianna et al 2013). (In 2D, the cell is “viewed from above” as it migrates ona flat 2D surface.) The neighbourhood of each point inside the shape represents a 2D projection of somesmall cylinder in 3D, containing both membrane and cytosol. Hence, active and inactive GTPases ( u and v )coexist at every point inside the given shape, as they do in our fixed domain 2D simulations.)Commonly, for the CPM, a scalar Hamiltonian, analogous to a potential is assumed to depend on thearea and perimeter of the cell, as well as the interface contact with other cells or empty space. Changesto the boundary of the cell are accepted or rejected stochastically, according to the net changes in theHamiltonian, as described in the Appendix. Our simulations include the following additional features: (1)solving the reaction-diffusion PDEs inside the evolving domain with Neumann boundary conditions at thecell boundaries and (2) modifying the Hamiltonian to depend on the local RD variables.In real cells, actin polymerizes into F-actin, and promotes protrusion of a cell edge. Hence, we link theF-actin variable F in the model to forces on the cell boundary, (by superimposing a chemically-dependentpotential H = ± βF for retractions(+) vs extensions(-) on the basic Hamiltonian, see Appendix). In variantsof the model that do not explicitly track F-actin, we assume that the GTPase u plays a similar role (i.e.,that u , like the GTPase Rac, locally promotes cytoskeleton assembly, creating a protrusive force at the celledge).Simulations are initiated with a circular cell and internal variables close to HSS but with a randomlyplaced peak of active GTPase, u somewhere inside the cell. Figure 16 shows a time series of a CPM simulationwith parameters that produced the absorbing wave simulations in the static domain. We observe three newtypes of dynamics (indicated by arrows 1, 2 and 3) resulting from cell movement. An initial burst in the pots, strips, and spiral waves in models for static and motile cells 19(a) u (b) v(c) u (d) v Fig. 10: Simulation of the non-conservative model (NC) with (a,b) default parameters (Table 4(NC)); (c,d) γ = 15 L , η = 15 L . Initial condition: (a,b) u = u ∗ except u = 1 on 0 ≤ x ≤ . v = v ∗ . (c,d) u = u ∗ =0 . u = 10 u ∗ on 0 . ≤ x ≤ . v = v ∗ = 3 . γ = 15 L (c) L = 1 (d) η = 15 Fig. 11: Final steady state pattern of the non-conservative model (NC) with most parameters from Ta-ble 4(NC), except the parameters indicated on the labels. (a) and (b) correspond to the steady state ofFig. 10(a,b) and (c,d) respectively. In (c) the shortened domain results in wave pinning; (d) Higher inac-tivation rate η = 15 results in bifurcating peaks. (a,b) corresponds to Fig. 5 (a,d) of (Verschueren andChampneys 2017), respectively.that broaden. We find a protrusion that is much broader than in Figure 16 (Arrow 2). Wave absorption islower than in Figure 16, so the cell edge is pushed further out.Figure 18 shows results for parameters corresponding to reflecting waves in a static domain. Here, becausethe cell boundary moves outwards, the waves are usually absorbed, rather than reflected. Occasionally, ifthe wave hits the cell edge tangentially, it is reflected (e.g. at 19 sec in the movie, upper left corner). We pots, strips, and spiral waves in models for static and motile cells 21(a) s = 5 (b) s = 8(c) s = 18 (d) s = 35 Fig. 12: Simulations of the combined model (CM). Parameters as in Table 4(CM2) but varying s , andHSS+noise initial condition as described in text. As we increase the actin feedback strength s , the behaviortransitions from slowly moving, repelling peaks to colliding peaks. At higher s , there is a rapidly oscillatingstanding wave pattern in some parts of the domain.furthermore observe three new wave dynamics in a moving cell with random bursts. A wave can break apartwhen it hits a burst (A), waves can merge (B), or avoid each other (C).As a last experiment, we simulate the formation of spots in the NC model (Figure 19). Since this modelhas no F-actin variable, we base the edge protrusion on the active GTPase u (assumed to act like Rac inpromoting local cytoskeleton assembly). The spots are highly dynamic and, as expected, lead to the formationof small protrusions (“filopodia”) (C). Furthermore, edge deformation also causes the spot pattern to change.When a protrusion forms (stochastically or by locally elevated u ), the spot in a region close to the protrusioncan split into two, one of which moves into the protrusion (A). We also see formation of new spots insideprotrusions (B). In summary we have explored extensions of the wave-pinning model (WP) (Mori et al 2008, 2011) thatcoupled the non-conservative variant proposed by (Verschueren and Champneys 2017) (NC) and the actinfeedback (AF) model of (Holmes et al 2012b). We found that the combined model (CM) borrows featuresfrom both, with moving peaks and wave trains, as well as more complex hybrid dynamics. At the same time,we were unable to find blinking localized spots as observed experimentally in (Robin et al 2016). Despitethe fact that the work of (Robin et al 2016) points to interactions of F-actin with the GTPase Rho, otherunknown factors, missing in our model, should be considered to explain such behavior.We used the local perturbation analysis (LPA) on each model variant. As noted before (Holmes et al2015), LPA recovers Turing analysis. States that are LPA stable are also Turing stable. LPA helps to identifypotentially interesting parameter ranges, including polarizable regimes that cannot be detected by Turinganalysis. At the same time, LPA does not predict details of patterns that emerge, nor accurate bifurcationpoints in the full PDE system. We also encountered examples where LPA identified apparent bifurcationsthat did not materialize as true regimes of behavior in the PDEs.As a second innovation, we simulated all model variants in 2D on both a static and a deforming domain.Previous work (Verschueren and Champneys 2017; Holmes et al 2012a) was concerned with fixed 1D domainsfor the PDEs. We hence showed that the patterns for the nonconservative (NC) model were primarily spots,not bands, whereas the AF model, while appearing to be less robust, produced a variety of moving peaks,bands, and waves, including spiral waves.Finally, our simulations of the models in a deforming domain mimicking a motile cell allowed us toconsider the connection to experimentally observed waves of actin (Inagaki and Katsuno 2017). We showedthat the internal dynamics of the models (and in particular the actin feedback model) have an interestingconsequence on the motility of a “model cell”. Indeed, the waves of high and low signaling levels led toformation of cellular protrusions, and gave rise to nontrivial motion in the deforming cell.
Acknowledgements:
LEK is funded by a Natural Sciences and Engineering Research Council of Canada(NSERC) Discovery Grant. YL was also supported by an NSERC postgraduate Fellowship. We are gratefulto Zachary Pellegrin for his development of a CPM reaction-diffusion solver. We thank members of theFeng-Keshet groups for feedback.
The LPA diagram for the nonconservative model (NC) is given in Figure 20. Due to the many intertwinedbifurcations, it is difficult to interpret this diagram, and we present it here only to demonstrate the limitationsof LPA.
Linear stability analysis (i.e. Turing analysis) (Turing 1952) is a more traditional method of analyzing thestability of reaction-diffusion systems which focuses on the stability of a HSS in response to perturbationsin the form of a global noise with infinitesimal height.Here we examine the relation between Turing and LPA. This was previously done by (Mata et al 2013)for the two-variable case in terms of eigenvalues, and for the general case by Theorem 4.1 of (Holmes 2014).(The proof for this theorem is quite involved.) Here we present a much more elementary argument for thetwo-variable case. pots, strips, and spiral waves in models for static and motile cells 23
Suppose in a general reaction-diffusion PDE with a slow-diffusing quantity u and a fast-diffusing quantity v . Assume that time has been rescaled so that the diffusion coefficient of the fast quantity is 1. ∂u∂t = δ ∇ u + f ( u, v ) ,∂v∂t = ∇ v + g ( u, v ) . The corresponding well-mixed (WM) system is ∂u∂t = f ( u, v ) ,∂v∂t = g ( u, v ) . The LPA system is the above plus an additional equation for the local variable: ∂u L ∂t = f ( u L , v ) . Suppose the PDE system has a homogeneous steady state (HSS) ( u ∗ , v ∗ ). The Jacobian of the well-mixedsystem and LPA system at the corresponding equilibrium is given by J W M = (cid:20) ∂ u f ∂ v f∂ u g ∂ v g (cid:21) , J LP A = ∂ u f ∂ v f ∂ u g ∂ v g ∂ v g ∂ u f = (cid:20) J W M ∗ ∂ u f (cid:21) , where the partial derivatives are understood to be evaluated at the HSS, and ∗ denote entries that areunimportant for later analysis. Notice that eigenvalues of J LP A are the two eigenvalues of J W M , plus ∂ u f .This is due to J LP A being block-lower-triangular. We will show that saying the HSS is Turing-unstable inthe limit δ → M ( q ) = J W M − Dq , D = (cid:20) δ
00 1 (cid:21) δ → −−−→ (cid:20) (cid:21) Suppose the HSS satisfy the condition for Turing instability, which has three conditions (see, for example,(Edelstein-Keshet 1988, Ch. 11.4)): Tr( J W M ) < , (4a)det( J W M ) > , (4b)det( M ( q )) < q > . (4c)Conditions (4a), (4b) are equivalent to saying that the HSS is a stable equilibrium for WM. Next, in thelimit of δ →
0, we computedet( M ) = det( J W M − Dq ) = det (cid:18) J W M − (cid:20) q (cid:21)(cid:19) = det( J W M ) − q ∂ u f Notice that by setting δ = 0, this equation is linear in q instead of quadratic. This means (4c) is equivalentto ∂ u f > ∂ u f is an eigenvalue for J LP A , this means the HSS is unstable in the LPA system.Conversely, suppose that the HSS is stable in WM and unstable in LPA. This means the eigenvalues of J W M all have negative real part, and ∂ u f >
0, which as shown above is equivalent to Turing-unstable.
Compared to Turing analysis, LPA has several advantages and disadvantages. LPA is essentially a “zeroth-order” expansion in δ , so it is only valid in the limit of δ → δ >
0, as opposed to Turing analysis. However, in this limit, LPA contains the Turing stability of thesystem. In particular, LPA-unstable is the same as Turing-unstable, whereas LPA-polarizable and LPA-stable regimes are Turing-stable. This follows the analysis above, and also (Mata et al 2013; Holmes 2014).This correspondence is illustrated in Fig. 21, where we show how the LPA regimes from Fig. 2(b) lines upwith Turing regimes.In pattern-forming regimes, LPA does not help predict the exact form of the pattern. This is mostrelevant to the actin feedback model, when there are many different possible patterns and a large number ofparameter regimes. Turing analysis also cannot predict the final pattern, but it does allow us to predict theinitial precursor pattern that forms and exists only for small t . In Fig. 22, we simulated the non-conservativemodel with parameters chosen such that only a single wave number is unstable. This results in a periodic,shallow precursor pattern which has the exact frequency as the unstable wave number.As the system continue to evolve, once a peak of the precursor pattern reaches a certain amplitude, itvery rapidly grow to the full size of the final pattern while suppressing nearby peaks. Other precursor peaksfarther away from the grown one might survive longer and eventually transition to full size, or be suppressedby another nearby peak which has transitioned sooner. These non-linear interactions cannot be captured byTuring analysis.In the special case that a static, periodic pattern forms, such as in the non-conservative model, the“minimal patch size” idea from (Painter and Hillen 2011), which is based on Turing analysis, can give anupper bound on the number of periods the pattern can have. Bifurcation plots were produced with XPPAUT (Ermentrout 2002; Doedel 1981) and Matcont (Dhooge et al2003). PDEs in 1D were solved using finite difference methods with CrankNicolson time stepping in Matlabwith ∆x = 0 . , ∆t = 0 . https://github.com/liuyue002/Wave-pinning-model .
10 Appendix: Cellular Potts Model simulations
In the CPM, a biological cell is represented by a set of contiguous lattice site in 2D (or 3D) all assigned anindex σ . (For a single cell, the index is 1 and the surrounding medium is given an index of 0.) Here we focusedon a 2D CPM model cell, representing a top-down view of a “biological cell” attached to a flat surface. Weuse the classic Hamiltonian H = λ a ( A − a ) + λ p ( P − p ) + JP. (5)Here the three terms represent the energetic cost for change of area A (cell contraction/expansion) awayfrom the preferred “rest area” a , a cost for elongation or shortening of the cell perimeter P away from apreferred “rest perimeter” p , and a term that describes an interfacial energy associated with the cell-mediuminterface. The weighting factors, λ a , λ p , J are adjusted to set the relative importance of the various energyterms.The perimeter P is approximated as in (Magno et al 2015). For each cell site, we calculate the numberof lattice sites within a certain radius r (here r = 3) in contact with the medium Then we take the sum overall boundary sites and rescale by ξ ( r ) (here ξ ( r ) = 18) to obtain a perimeter approximation: P = 1 ξ ( r ) (cid:88) x : σ ( x )=1 (cid:88) y : {| x − y | 005 to each pixelimplies that cell area is 0.3. We rescaled the parameters η, k, γ, α, θ, k n , k s with L , where L = 10. CPMparameters a, p are in terms of pixels (see figure captions).The initial conditions are v = 2 . , u = 0 and F = 0 everywhere. We introduce a randomly placed circularspot (radius 3 pixels) of u = 5 to represent an initial random burst of active GTPase. Then every 100 MCSwe add another randomly placed spot of elevated GTPase u + = 15 (higher than the initial burst to preventdecay) to depict stochastic bursts of GTPase activation in the cell.For the model without F-actin, we have slightly different initial GTPase fields: v = 1 . , u = 0 . u = 0 . dt = 1 e − 6, so 1000 iterations per MCS.Within a MCS, after every accepted membrane extension or retraction we update the GTPase fields asfollows. After an extension we set u (target)= u (source) and subsequently rescale the level of u throughoutthe cell to conserve mass. After retraction we set u(target)=0 and then u(target) is equally distributedthroughout the whole cell. The same operations are carried out for v . For F-actin, we do the same but donot redistribute, i.e. F (target)= F (source) for extensions and F (target)=0 for retractions. References Alnaes MS, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richardson C, Ring J, Rognes ME, Wells GN (2015) The fenicsproject version 1.5. Archive of Numerical Software 3(100):9–23Bretschneider T, Anderson K, Ecke M, M¨uller-Taubenberger A, Schroth-Diez B, Ishikawa-Ankerhold HC, Gerisch G (2009) Thethree-dimensional dynamics of actin waves, a model of cytoskeletal self-organization. Biophysical journal 96(7):2888–2900Buttensch¨on A, Liu Y, Edelstein-Keshet L (2019) Cell size, mechanical tension, and gtpase signaling in the single cell. arXivpreprint arXiv:190810840Buttensch¨on A, Liu Y, Edelstein-Keshet L (2019) Cell size, mechanical tension, and gtpase signaling in the single cell. arXivpreprint: 190810840Dhooge A, Govaerts W, Kuznetsov YA (2003) Matcont: a matlab package for numerical bifurcation analysis of odes. ACMTransactions on Mathematical Software (TOMS) 29(2):141–164Doedel EJ (1981) Auto: A program for the automatic bifurcation analysis of autonomous systems. Congressus Numerantium30:265–284Doubrovinski K, Kruse K (2011) Cell motility resulting from spontaneous polymerization waves. Physical review letters107(25):258,103Edelstein-Keshet L (1988) Mathematical Models in Biology, vol 46. SIAMEdelstein-Keshet L, Holmes WR, Zajac M, Dutot M (2013) From simple to detailed models for cell polarization. PhilosophicalTransactions of the Royal Society of London B: Biological Sciences 368(1629):20130,003Ermentrout B (2002) Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers andStudents. SIAMGrieneisen V (2009) Dynamics of auxin patterning in plant morphogenesis-a multilevel model study. PhD thesis, UniversityUtrechtHolmes WR (2014) An efficient, nonlinear stability analysis for detecting pattern formation in reaction diffusion systems.Bulletin of Mathematical Biology 76(1):157–183Holmes WR, Edelstein-Keshet L (2016) Analysis of a minimal rho-gtpase circuit regulating cell shape. Physical Biology13(4):046,001Holmes WR, Carlsson AE, Edelstein-Keshet L (2012a) Regimes of wave type patterning driven by refractory actin feedback:transition from static polarization to dynamic wave behaviour. Physical Biology 9(4):046,0056 Yue Liu et al.Holmes WR, Lin B, Levchenko A, Edelstein-Keshet L (2012b) Modelling cell polarization driven by synthetic spatially gradedrac activation. PLoS Computational Biology 8(6):e1002,366Holmes WR, Mata MA, Edelstein-Keshet L (2015) Local perturbation analysis: a computational tool for biophysical reaction-diffusion models. Biophysical Journal 108(2):230–236Inagaki N, Katsuno H (2017) Actin waves: Origin of cell polarization and migration? Trends in cell biology 27(7):515–526Jilkine A, Edelstein-Keshet L (2011) A comparison of mathematical models for polarization of single eukaryotic cells in responseto guided cues. PLoS Computational Biology 7(4):1–15Kuznetsov Y (2004) Elements of Applied Bifurcation Theory, 3rd edn. Springer Science & Business MediaLiu Y (2019) Analysis of pattern formation in reaction-diffusion models for cell polarization. MSc thesis, University of BritishColumbiaMagno R, Grieneisen VA, Mar´ee AF (2015) The biophysical nature of cells: potential cell behaviours revealed by analytical andcomputational studies of cell surface mechanics. BMC Biophysics 8(1):8Mar´ee AF, Grieneisen VA, Hogeweg P (2007) The cellular potts model and biophysical properties of cells, tissues and morpho-genesis. In: Single-cell-based models in biology and medicine, Springer, pp 107–136Mar´ee AF, Komba M, Finegood DT, Edelstein-Keshet L (2008) A quantitative comparison of rates of phagocytosis and digestionof apoptotic cells by macrophages from normal (balb/c) and diabetes-prone (nod) mice. Journal of applied physiology104(1):157–169Mata MA, Dutot M, Edelstein-Keshet L, Holmes WR (2013) A model for intracellular actin waves explored by nonlinear localperturbation analysis. Journal of Theoretical Biology 334:149–161Mori Y, Jilkine A, Edelstein-Keshet L (2008) Wave-pinning and cell polarity from a bistable reaction-diffusion system. Bio-physical Journal 94(9):3684–3697Mori Y, Jilkine A, Edelstein-Keshet L (2011) Asymptotic and bifurcation analysis of wave-pinning in a reaction-diffusion modelfor cell polarization. SIAM Journal on Applied Mathematics 71(4):1401–1427Otsuji M, Ishihara S, Kaibuchi K, Mochizuki A, Kuroda S, et al (2007) A mass conserved reaction–diffusion system capturesproperties of cell polarity. PLoS Computational Biology 3(6):e108Painter KJ, Hillen T (2011) Spatio-temporal chaos in a chemotaxis model. Physica D: Nonlinear Phenomena 240(4-5):363–375Robin FB, Michaux JB, McFadden WM, Munro EM (2016) Excitable rhoa dynamics drive pulsed contractions in the early c.elegans embryo. BioRxiv p 076356Scianna M, Preziosi L, Wolf K (2013) A cellular potts model simulating cell migration on and in matrix environments. Mathe-matical Biosciences & Engineering 10(1):235–261Turing AM (1952) The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London Series B,Biological Sciences 237(641):37–72Verschueren N, Champneys A (2017) A model for cell polarization without mass conservation. SIAM Journal on AppliedDynamical Systems 16(4):1797–1830Walther GR, Mar´ee AF, Edelstein-Keshet L, Grieneisen VA (2012) Deterministic versus stochastic cell polarisation throughwave-pinning. Bulletin of Mathematical Biology 74(11):2570–2599Zmurchok C, Bhaskar D, Edelstein-Keshet L (2018) Coupling mechanical tension and gtpase signaling to generate cell andtissue dynamics. Physical Biology 15(4):046,004pots, strips, and spiral waves in models for static and motile cells 27(a) Wave pinning, t = 1 . t = 30(c) Non-conservative extension, t = 0 . t = 0 . Fig. 13: 2D simulations of the wave pinning (a,b) and non-conservative (c,d) models, using the same param-eters as in 1D (Table 4(WP) and (NC)). Left: u , Right: v . For each model, two snap shots are shown: onewhen the pattern begin to take shape, and another after the system reached steady state. s = 12(b) s = 18(c) s = 27(d) s = 36 Fig. 14: 2D simulations of the actin feedback (AF) model, with parameters from Table 4(AF) and initialconditions described in the text. These snapshots are taken after the patterns have fully developed. As s increases, blobs transitions into thinner and thinner spiral waves. See movies at https://imgur.com/a/61GwiA9 . pots, strips, and spiral waves in models for static and motile cells 29(a) s = 8(b) s = 12(c) s = 18(d) s = 4 Fig. 15: Simulations of the combined model (CM) in 2D, with parameters from Table 4(CM2) and HSS+ noise initial condition. There is a transition from spots to spiral waves near s = 12. See movies at https://imgur.com/a/a0u57GQ . Fig. 16: Snapshots of 2D CPM simulation with parameters from the absorbing waves in a static domain.Visualized is F-actin ( F ) that promotes protrusions ( H = ± βF ). Arrows indicate examples of interestingdynamics. Snapshots are 20 MCS apart. Parameters are: D a = 0 . η = 15, k = 6, n = 3, γ = 30, (cid:15) = 0 . k s = 7 . k n = 24, s = 30. CPM Parameters are: a = 12000, λ a = 2, p = 500, λ p = 20, J = 50, r = 3 , ξ ( r ) = 18, β = 150, T = 100. Movie link https://imgur.com/a/7OmgctR . pots, strips, and spiral waves in models for static and motile cells 31 Fig. 17: Snapshots of 2D CPM simulation with parameters from the oscillating waves in a static domain.Visualized is F-actin ( F ) that promotes protrusions ( H = ± βF ) Arrows indicate examples of interestingdynamics. Snapshots are 20 MCS apart. Parameters are as in Fig. 16, but with k = 1 . s = 18. CPMparameters are as in Fig. 16, but with β = 50. Movie link https://imgur.com/a/eIAjr59 Fig. 18: Reflecting wave parameter set. Visualized is F-actin ( F ) that promotes protrusions ( H = ± βF ).Snapshots in A,B,C are 10,20,20 MCS apart respectively. Parameters are as in Fig. 16, but with k = 1 . s = 27. CPM parameters are as in Fig. 16, but with β = 50. Movie link https://imgur.com/a/FDCn3NY pots, strips, and spiral waves in models for static and motile cells 33 Fig. 19: Spots. Visualized is Rac ( u ) that promotes protrusions ( H = ± βu ). Snapshots in A,B,C are 10,5,25MCS apart respectively. Parameters are: D a = 0 . η = 60, k = 6, n = 3, γ = 30, c = 1, θ = 18, α = 6.CPM Parameters are: a = 10000, λ a = 0 . p = 1000, λ p = 0 . J = 40, r = 3 , ξ ( r ) = 18, β = 200, T = 20.Movie link https://imgur.com/dADHOnS . Fig. 20: LPA bifurcation diagram for the combined model (CM) with respect to the parameter s . Otherparameters as in Table 4(CM2). There are many apparent branches of periodic solutions. In a parameterrange around s = 20, there are no stable equilibria nor stable periodic solutions even though the systemremains bounded, which suggests the presence of chaos. pots, strips, and spiral waves in models for static and motile cells 35(a) LPA (b) Turing Fig. 21: Comparison of LPA and Turing bifurcation diagrams for the non-conservative model. (a) is a zoom ofthe LPA diagram from Fig. 4(b). (b) is the Turing bifurcation diagram reproduced from Fig. 5 of (Verschuerenand Champneys 2017), using the same parameters. Observe that both the LPA-stable (I) and the LPA-polarizable (II) regimes in (a) located to the left of γ c = 16 . 765 correspond to the Turing-stable regimebelow the blue curve in (b). The LPA-unstable regimes (III, IV) correspond to the Turing-unstable regimeabove the curve. The curve passes through δ = 0 , γ = γ c . The bifurcation boundary between Regimes I andII, and between III and IV cannot be detected by Turing analysis. Given that numerical simulations haveshown that the PDE produces the same behavior (Fig. 11(a)) in both Regimes III and IV, it is possible thatthese are not distinct regimes for the PDE. Overall, the LPA diagram (a) can be seen as a vertical slice ofthe Turing diagram (b) at δ = 0, with additional bifurcation boundaries that separates the LPA-stable andLPA-polarizable regimes. Fig. 22: Simulations for the non-conservative model with random initial conditions u = u ∗ (1 +Unif( − . , +0 . , v = v ∗ with default parameters from Table 4(NC) except γγ