Dynamics of a cell motility model near the sharp interface limit
DDynamics of a cell motility model near the sharpinterface limit (cid:73)
Nicolas Bolle Department of Mathematics and Statistics, The College of New JerseyEwing Township, NJ
Matthew S. Mizuhara ∗ Department of Mathematics and Statistics, The College of New JerseyEwing Township, NJ
Abstract
Phase-field models have recently had great success in describing the dynamicmorphologies and motility of eukaryotic cells. In this work we investigate theminimal phase-field model introduced in [7]. Rigorous analysis of its sharp in-terface limit dynamics was completed in [17, 18], where it was observed thatpersistent cell motion was not stable. In this work we numerically study thepre-limiting phase-field model near the sharp interface limit, to better under-stand this lack of persistent motion. We find that immobile, persistent, androtating states are all exhibited in this minimal model, and investigate the lossof persistent motion in the sharp interface limit.
Keywords: phase-field, keratocyte motion, traveling wave, rotating cell
1. Introduction
Eukaryotic cell motility underlies numerous biological processes includingthe immune response and cancer metastasis. Cell motion is initiated and main-tained by an evolving cytoskeleton comprised of actin and myosin proteins ca-pable of driving a wide range of motility modes. In recent years a varietyof modeling techniques have been highly successful in replicating, explaining,and predicting cell morphologies observed in experimental settings, see, e.g.,[1, 3, 8, 9, 19, 20, 21, 23, 25]. In this work we focus on the bridge between two (cid:73) https://doi.org/10.1016/j.jtbi.2020.110420. This work is licensed under CC BY-NC-ND4.0. To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-nd/4.0 ∗ Corresponding author
Email addresses: [email protected] (Nicolas Bolle), [email protected] (Matthew S.Mizuhara) Present address: Department of Mathematics, The Ohio State University, Columbus, OH
Preprint submitted to Elsevier August 3, 2020 a r X i v : . [ q - b i o . CB ] J u l pecific modeling paradigms: free boundary problems and phase-field models.Free boundary problems track the cell boundary via a curve (2D) or surface(3D) whose evolution is governed a geometric evolution equation, often dictatedby boundary data of a differential equation solved on the interior. Phase-fieldmodels, on the other hand, use an evolving order parameter whose finite widthtransition layer between phases tracks the cell boundary. Phase-field modelsavoid difficulties of explicitly discretizing and tracking the moving interface,making them ideal for numerical simulation.We are motivated by the 2D phase-field model for keratocyte fragments (e.g.,lacking a nucleus) studied in [7]. It is a minimal version of a more general modelintroduced by Ziebert, et al. in [31]. The original model in [31] has been ex-tended to include spatial adhesion dynamics [30], non-homogeneous substrateeffects [13, 16, 24], and interacting dynamics of multiple cells [14]. These ex-tended models exhibit a wide variety of dynamical modes and can be used tounderstand the complex morphologies of dynamics cells. On the other hand, theminimal model of [7] allows for rigorous mathematical analyses of the model.For example, in 1D, necessary conditions for the existence and stability of per-sistent motion were proven [7]. In 2D, sufficient conditions for existence andnon-existence of persistent motion were studied in [17, 18].However, a numerical exploration of the simplified phase-field model remainsunexplored in 2D. Thus, our primary goal of this work is to numerically studythe minimal 2D phase-field model introduced in [7] over a range of parameters.We find that this simplified version is surprisingly capable of exhibiting a rangeof motions in qualitative agreement with more sophisticated models, includingstationary, persistently traveling, and rotating modes.The minimal model admits a non-trivial sharp interface limit: an asymp-totic reduction of the phase-field model in the limit that the width of the diffuseinterface (the location of the cell boundary) tends to zero, transforming themodel into a free boundary problem. Rigorous analysis of the sharp interfacelimit was completed in both 1D [7] and in 2D [17, 18], where sufficient condi-tions for existence of traveling wave solutions were proved. These analyses werethus able to provide insight into the minimal biophysical mechanisms that arenecessary to drive these motility modes, and which modes require more complexmechanisms.However, previous numerical simulations of the 2D sharp interface limitshowed that persistent motion was unstable, and more crucially, does not ex-ist when certain symmetry is present (see Section 4 for details). Therefore, asecondary goal of our work is to explore persistent motion in the phase-fieldequation, and understand its existence as we approach the sharp interface limit.
2. Model
Keratocytes are prototypical for the experimental and mathematical studyof cell motion. Their characteristic cell length/width is two orders of magnitude2arger than the height while motile, hence they are well described by 2D mod-els (for recent advances in 3D models see, e.g., [10, 26, 29]). Keratocytes areadditionally able to exhibit persistent motion over many times the cell lengthwith approximately constant cell shape, making them ideal starting points forthe study of motion [20, 27].We recall the following key factors contributing to cell motion, and referthe reader to [19] for a more detailed review. A crawling cell maintains self-propagating motion via internal forces generated by actin polymerization. Actinmonomers bind together to form filaments which create a dense network at theleading edge of the cell, known as the lamellipod. The cell’s leading edge pro-trudes via growth of actin laments at the cell membrane and degradation of thefilaments towards the interior of the cell, a process known as actin treadmilling.Intercellular adhesion complexes form ligand bonds to the substrate in order totransform this propulsion force into traction forces. Myosin motors interact withactin filaments to generate contractile forces. Acto-myosin interactions contractof the rear part of the cell, pulling the rear of the cell, and allowing for persistentmotion. In idealized mathematical settings, such persistent motion is describedby traveling wave solutions. Such motility remains a rich area of study: per-sistent motion been observed in both myosin-inhibited [11], and actin-inhibited[12, 22] cells. It has additionally been observed that random fluctuations aresufficient to spontaneously switch a cell from a symmetric, non-motile state tomotile states [2, 10].Moreover, by varying biophysical parameters (such as actin polymerizationstrength, substrate adhesion/elasticity properties, or myosin motor strength),cells additionally exhibit a wide range of motility modes beyond persistent mo-tion, such as stick-slip (oscillations in translational velocities) and bipedal (leftand right sides alternating forward motion) motions [4, 5]. As such, there is adeep need to understand the interactions of various biophysical pathways leadingto such a variety of behaviors.Finally, of particular interest are “rotating” cells, experimentally observed in[15], where cells remained essentially stationary but experienced lateral waves ofprotrusions of the membrane. This was caused by the expression of a particularkinase (MLCK) leading to an increase of myosin activity in the cell’s lamellipod.The increased myosin activity was hypothesized to enhance actin depolymeriza-tion (or alternatively, to sweep the actin network toward the center of the cell),thus resulting in shorter edge protrusion lifetimes. Shorter edge lifetimes cannotsufficiently polarize the cell to generate motion in a single direction, resultingin lateral waves of actin protrusion at the cell boundary.
The following 2D phase-field model for cell motility was originally introducedin [31]; we study a slightly simplified form (e.g., omitting some non-linearities)which is amenable to rigorous analysis (see Section 2.4): ∂ t ρ = D ρ ∆ ρ − τ − F ( ρ ) − α ( ∇ ρ ) · P (2.1) ∂ t P = D P ∆ P − τ − P − ζ ∇ ρ. (2.2)3arameter Description Typical values D ρ interface stiffness 1 µm /sτ time scale of curvature motion 1 sα adhesion site formation rate .5 - 3 µm/sD P actin diffusion . µm /sτ actin depolymerization 10 sζ actin polymerization rate 1 - 2 µm/sR cell size (radius) 5 - 18 µm Table 1: Parameter descriptions of original phase-field model (2.1) - (2.2) together with typicalvalues, from [31]. Here F ( ρ ) = (1 − ρ )( δ ( ρ ) − ρ ) ρ, (2.3) δ ( ρ ) = 12 + µ (cid:20)(cid:90) Ω ( ρ ( x, y, t ) − ρ ( x, y, dxdy (cid:21) . (2.4)The phase-field variable ρ = ρ ( x, y, t ) takes value ρ ≈ ρ ≈ A ( t ) = (cid:90) Ω ρ ( x, y, t ) dxdy, (2.5)is approximately preserved due to penalization from δ , see [31] for more details.Descriptions of parameters and typical values for keratocytes are presented inTable 1.The interface stiffness D ρ represents the ratio of surface tension of the cellmembrane to friction with the substrate [25, 31]. Together D ρ and τ dictatethe rate of passive curvature driven motion and the size of the diffuse interface.In particular the width of the diffuse interface scales with (cid:112) D ρ τ [32]. Activetransport of ρ occurs along the vector field P = P ( x, y, t ) which represents theactin network from a macroscopic point of view. At each ( x, y ) the vector fieldpoints in the average local orientation of actin filaments, and the magnitudeof the vector represents the degree of orientation (e.g., many well-aligned actinfilaments result in a larger magnitude vector) [31].Physically, advection requires formation of substrate adhesions to gener-ate traction forces. To that end, a transmembrane complex of proteins formligand bonds to the substrate and connect to the actin network. In [30] ad-hesion dynamics are explicitly modeled via a PDE for an auxiliary scalar field A = A ( x, y, t ) whose dynamics additionally encode substrate deformations (asa visco-elastic medium), and transport requires adhesion formation so that ∂ t ρ ∼ A P · ∇ ρ . For simplicity, we assume that adhesion is formed instanta-neously and uniformly with the vector field P .Dynamics of actin filaments P are regularized by diffusion and experienceglobal decay due to depolymerization. As actin polymerization is localized to4he boundary of the cell, the source term for P is given by − ζ ∇ ρ , where ζ isthe rate of actin polymerization.Competition between advection by P and curvature motion flow from theAllen-Cahn contribution constitute the main dynamics of interest: one expectsthat if | P | is sufficiently small then the cell remains immobile and if | P | issufficiently large then the cell has sufficiently many active internal forces togenerate motion.In [24], numerical simulations of the more general phase-field model, includ-ing additional non-linear effects from heterogeneous myosin contraction, non-linear dynamics of adhesion complex formation, and substrate viscoelasticity.In that more complex setting, they observe several motility modes includingseveral types of rotating lamellipod solutions. To study the sharp interface limit, one first non-dimensionalizes and choosesan appropriate scaling to give rise to a non-trivial limit. To that end introduceparameters ε := (cid:112) D ρ τ R , ˜ α := αRD ρ , ˜ ζ := ζτ R , β := ˜ α ˜ ζε , (2.6)where R is the characteristic size of the cell (say the radius). Then, undercertain parameter assumptions, (2.1)-(2.2) can be written in non-dimensionalform [7]: ∂ t ρ = ∆ ρ − ε W (cid:48) ( ρ ) − P · ∇ ρ + λ (2.7) ∂ t P = ε ∆ P − ε P − β ∇ ρ, . (2.8)See Appendix A for a detailed derivation. Using typical values of parametersfrom Table 1 we see ε ∼ . β ∼
100 for keratocytes.The non-dimensional parameter ε represents the relative size of the diffuseinterface width (to the cell size). As actin polymerization occurs only where |∇ ρ | (cid:54) = 0, the parameter ε also then dictates the width of the region wherepolymerization occurs. While one would ideally like to decouple surface tensionfrom actin polymerization, this would require complications of the model beyondthe scope of this work. The non-dimensionalized polymerization β is a ratio ofnon-dimensional adhesion site formation and actin polymerization rates to ε .In original parameters β = αζR D ρ and thus it can be thought of as analogous toa P´eclet number, with advection driven by actin polymerization.In (2.7)-(2.8) the phase-field parameter ρ now has an O ( ε ) thick transitionlayer and the sharp interface limit can be analyzed as ε →
0. Since analysisin the sharp interface limit assumes constant β , in the sharp interface limitwe require that the product ˜ α ˜ ζ →
0. For example, in the limit of vanishingcell surface tension we have D ρ → ε → | ζ | ∼ | D ρ | to ensure β is kept constant. Naively, one may expect that if ζ → α → β is sufficiently large (see Section 2.4).We remark that for simplicity of analysis, area preservation in (2.7)-(2.8) isenforced via the Lagrange multiplier λ ( t ) = 1 | Ω | (cid:90) Ω ε W (cid:48) ( ρ ) + P · ∇ ρ dx, (2.9)so that W is fixed to be the standard Allen-Cahn double-well potential W ( z ) = z (1 − z ) . We additionally note that λ = λ ( t ) passively encodes myosinmotor effects, by virtue of enhancing contraction of the cell membrane (e.g.,when λ <
00 00 0055 55 55 55 5555 t=4.8t=5.0
Figure 1: Sample snapshots of the three long time behaviors arising from (2.7) - (2.8) at twonearby times: (left) stationary ( ε = . β = 80), (center) rotating ( ε = . β = 110), and(right) persistent motion ( ε = . β = 130). Color indicates value of ρ and arrows representthe vector field P . One reason we study the system (2.7)-(2.8), rather than the more generalmodel, e.g. in [24], is that the system (2.7)-(2.8) is amenable to rigorous math-ematical analysis: it is possible to derive dynamics in the limit ε → sharp interface limit . In [7], both the pre-limiting and sharp in-terface limit dynamics in 1D are rigorously studied. In particular necessaryconditions for the existence and stability of traveling wave solutions are estab-lished, corresponding to persistently traveling cells. Additionally, Berlyand etal. establish the 2D sharp interface limiting equation, recovering a geometricevolution equation for planar curves, Γ( s, t ) representing the boundary of the6ell. They show that the normal velocity V of the curves at each moment evolveby V = κ + β Φ( V ) − | Γ( t ) | (cid:90) Γ κ + β Φ( V ) ds, (2.10)where κ = κ ( s, t ) represents the curvature at location s and time t , and Φ : R → R is a fixed non-linear function whose form is explicit and depends on the double-well potential W . The integral term, as before, enforces volume preservation.Analysis and numerical simulation of (2.10) was completed in [17, 18]; we brieflyreview the relevant results.Due to the importance of persistently moving cells, the authors consideredtraveling wave solutions of (2.10). These are solutions of the formΓ( s, t ) = Γ ( s ) + v t, where v is a fixed vector representing the velocity of the cell and Γ ( s ) is theunknown cell shape. As expected by physical considerations, we require suffi-ciently large β in order for traveling wave solutions to exist, as β captures therelative strength of actin polymerization. Surprisingly, it was also shown thatif Φ was an even function then traveling waves could not exist, regardless of thevalue of β . In particular, the standard Allen-Cahn potential, as considered inthe current work, results in even Φ.It is thus natural to suspect that the simplified phase-field model (2.7)-(2.8) cannot support persistent motion, perhaps because we have eliminatedsymmetry breaking effects of myosin motors. However, as we will find below, theminimal phase-field model is capable of exhibiting a wide range of motions forvarious parameter regimes, thanks to the finite width transition layer allowingfor non-trivial dynamics of the vector field P .
3. Numerical simulation of the phase-field model near the sharp in-terface
To understand dynamics near the sharp interface limit, we investigate thelong-time dynamics of the phase-field model for various values of β and forsmall values of ε . Simulations of (2.7)-(2.8) are done using an explicit finitedifference method with centered space steps ( h = . , with periodic boundary conditions. Time steps ( dt = 8e-5)are taken sufficiently small to ensure convergence of the simulations: takingsmaller time steps did not qualitatively affect any results. Additionally thecell area is small compared to the domain size to ensure that there were noboundary effects caused by the periodic boundary conditions. Again, takinglarger domain size does not qualitatively change dynamics. The non-local termis approximated by Riemann sums (left or right hand choice was irrelevant dueto simulations being on a torus). We numerically tracked total enclosed area, A ( t ) = (cid:82)(cid:82) ρ ( x, y, t ) dxdy , to ensure that it was conserved over time. As ε varies,while the diffuse interface width varies, we emphasize that the enclosed area7oes not change, as A ( t ) ≈ A (0) is determined by initial conditions alone. Thetotal time of integration is T = 5.Initial conditions for the phase-field are a circular cell. We consider bothpolarized and non-polarized initial conditions for the actin field: for polarizedinitial conditions we assume that the actin field on the interior of the cell isconstant and pointing to the right, with a small random perturbation to ensurerobustness of the results. For non-polarized initial conditions we take randominitial conditions for the actin field with | P | (cid:28)
1. Moreover, we assume suffi-ciently long time integration to ignore transient effects. We find that, regardlessof initial conditions, the long time behavior was qualitatively the same. Thus,we report only on the results of the polarized initial conditions.After integrating past a transient time for the cell to reach stable behavior,we record the resultant dynamics. The resultant data is summarized in Figure2. We note that taking smaller values of ε became computationally expensive,due to the singular nature of the evolution of (2.7)-(2.8).For all solutions which had non-trivial motion, we track the center of mass ofthe cell in order to study its trajectory. We then calculate the average radius ofcurvature of this trajectory in order to distinguish between persistently movingcells (straight line trajectory) and rotating cells (circular trajectory). Thereis an arbitrary distinction between rotating and persistently moving cell, as apersistently moving cell may have a slight turn due to numerical artifact ornon-symmetric initial conditions. Thus, to distinguish the two cases we set athreshold radius of curvature to be 10 (cid:112) | Ω | : any trajectories with larger radiiof curvature are defined to be persistently traveling. ImmobileRotatingPersistent motion .02 .0570150 ε β .02 .05 ε β Figure 2: Long time behavior of center of mass for a range of parameters. We observepersistent motion (cells moving with constant shape in a straight line), stationary states (cellsrelaxing to a circular shape without any motion), and rotating motion (asymmetric cells whosecenter of mass traces a circular path). (left) Trajectories of the center of mass after a transientperiod and (right) classification of type of motion.
Over the parameter range considered, we observe three modes of motion: immobile, persistent, and rotating solutions . Figure 1 presents a snapshotand parameter values giving rise to each type of motion.8 mmobile.
We find that for any fixed ε , there is a critical β cr ( ε ) >
0, so thatfor any β < β cr no motion is possible. This agrees qualitatively with the theorydeveloped in the sharp interface limit, as discussed in Section 2.4. Heuristically,stationary states are expected for sufficiently small β , since taking β = 0, themodel simplifies to the volume preserving Allen-Cahn equation (as P = 0)which asymptotically relaxes to circular steady states. Then, small values of β constitute a regular perturbation and so stationary states are expected topersist for sufficiently small β . Persistent motion.
Persistent motion represents a cell traveling with con-stant shape. Since we use a symmetric double-well potential W , it is perhapssurprising, given the theory developed in the sharp interface limit, that anypersistent motion is possible. We do observe that for sufficiently small values of ε these persistent motions no longer exist in our model, agreeing qualitativelywith the results of [18]. We conclude that stabilization of persistent motion inthe sharp interface limit requires additional myosin motor effects [28], or perhapsalternative scalings to be considered in (2.7)-(2.8). Rotating states.
Rotating states indicate that the effective actin polymer-ization strength is not sufficiently strong to overcome surface tension, creatinga lateral wave of actin propagation along the cell boundary and resulting ina rotating wave of protrusion in the cell, which results in a rotating solution.Simulations show that rotating states emerge when ε is relatively small. Since ε scales with surface tension and ˜ α ˜ ζ = βε , we see that when ε is relatively smallwe have, e.g., ˜ ζ = O ( ε ) so that effective actin polymerization is a further or-der of magnitude smaller. Traveling actin waves leading to such rotating stateshave been investigated in the more complex version of the phase-field model in[24]. Indeed, in [24], the formation of such waves was explained via shockwavesfrom a Burger-like equation, whose non-linear shocks are driven by a quadraticcontribution ∼ |∇ ρ | .Moreover, rotating cells in experiments were explained via a shortening ofedge protrusion lifetime [15]. In numerical simulations, since protrusions aregenerated exclusively from the vector field P , we assume that the edge protru-sion lifetime is related to the time-scale of non-trivial dynamics of P . To thatend, note that on the interface |∇ ρ | = O (1 /ε ), so ε |∇ ρ | = O (1). Thus freezing ∇ ρ for simplicity and writing ∂ t P = ε ( ε ∆ P − P − βε ∇ ρ ), we see that thetime-scale of P relaxing to equilibrium is also O ( ε ). Thus small values of ε correspond to short edge protrusion lifetimes. Again, we emphasize that ε dic-tates both the time scale of dynamics of P as well as the width of the interfacelayer. As aforementioned, the particular scalings of ε in the model are requiredfor comparison to the sharp interface limit. To better understand and explainthe onset of rotational motion in the phase-field model, one must relax thesecoefficient restrictions to decouple these two effects. This is beyond the scope ofthe present work (whose focus is the relationship to the sharp interface limit),though this is an important question for future work.We remark that we did not observe more complicated morphologies suchas ameobid or two wave solutions with both actin waves traveling in the samedirection, as observed in [24]. This suggests that such dynamics are driven9y more complicated biophysical mechanisms, including hetereogeneities in themyosin motor density or in the adhesion complex formation. However, it issurprising that rotating states exist in our simplified model. Our work thussuggests that even with severely weakened myosin contraction, such motilitymay already be possible in cells. While we have already established that persistent motion does not seempossible for any fixed β as ε →
0, we further analyze the trajectory data ofsimulation results as ε → ε or β results in a decrease in the radius of the curvature of the trajectory, seeFigure 3. This suggests that as ε →
0, the cell has a stationary center of mass.While certainly this does not preclude the existence of other motile cells withstationary center of mass, we did not observe any such modes in our simulations.These data in particular provide evidence for why the sharp interface limit(2.10) may not exhibit persistent traveling waves: for small ε we see that theonly stable states which seem to survive are stationary states and rotating states.As ε →
0, even rotating states have center of mass trajectories with smaller andsmaller radii. So, in the limit ε → =100=110=120=130=140=150 ε r = 80= 90= 100= 110= 120= 130= 140= 150 ε Figure 3: (left) Average radius of circular trajectory for rotating solutions. The data showthat, for fixed β , decreasing ε results in a smaller radius of curvature. That is, rotating cellstrace smaller circles. This suggests that for fixed β , as ε →
0, the cell becomes stationary.(right) Distance traveled for both rotating and traveling states show non-monotone behavioras ε is varied: at ε = .
035 we observe (for sufficiently large β ) that distance traveled ismaximized. On the other hand, for fixed ε dependence on β is monotone. .2. Non-monotone dependence of distance traveled on parameters We investigate the dependence of cell speed on model parameters. To thatend, we integrate (2.7)-(2.8) to the end time T = 5 and omit the first 75% ofdata (to disregard transient effects). Then calculate the distance traveled bythe center of mass’s trajectory, s ( t ): d = (cid:88) | s ( t i +1 ) − s ( t i )) | , see Figure 3. We observe that the relationship between distance traveled andparameters is dependent on the motility mode: over parameter ranges wherethe cell is traveling persistently, if one fixes β and increases ε then the distancetraveled increases monotonically. Indeed, at the conclusion of each simulation,we also compute (cid:107) P (cid:107) := (cid:90) (cid:90) | P ( x, y ) | dxdy. This quantity can be thought of as a measure of the total forces from actinpolymerization in the cell. We indeed found that for fixed β , an increase in ε resulted in an increase in (cid:107) P (cid:107) , see Figure 4.However, over parameter ranges where the cell is rotating there is non-monotone dependence on the distance traveled. In that case, the distance trav-eled is maximized when ε = . β is large enough. As ε is the ratioof diffuse interface width (itself a ratio of membrane surface tension to substratefriction) to cell size, this suggests that a ratio of the cell size to the cell surfacetension may be an optimizable quantity for maximizing cell speed. Since ournumerics suggest that the optimal ε is constant for large enough β , we con-jecture that this optimal ratio may be independent from actin polymerizationand adhesion site formation rates (provided they are sufficiently large). To ourknowledge, this has not been explored experimentally, and thus merits study. =70=80=90=100=110=120=130=140=150 ε || P || Figure 4: The L P at steady state (after transient effects) shows that as ε increases,the total amount of actin increases. We interpret this as larger lamellipod size resulting in alarger region of actin polymerization and thus stronger protrusion forces. . Conclusion In this work we numerically investigate the phase-field system (2.7)-(2.8),with particular interest in the dynamics near the sharp interface limit ε → δ > ε , but we change notation herefor clarity), whose limit as δ → δ > β < β cr then all dynamics relaxed to circular states. For β > β cr if traveling wave solutions existed they were unstable. Instead, all dynamicsresulted in one of two long term dynamics. For larger values of δ , simulationsgave rise to cells moving in a bipedal fashion [5]. For smaller values of δ , theintermediate system gave rise to rotating cells. Moreover the rotating cells be-came more circular as δ →
0. Importantly, all non-trivial motions arose onlywhen considering asymmetric double-well potentials.Our numerics agree qualitatively with previous results of the intermediatesystem in many ways. We find that for sufficiently small β all cells converge tocircular stationary states. We also observe rotating cells which become closer tocircular near the sharp interface limit, providing further evidence that the sharpinterface limit cannot support non-trivial motion. On the other hand there aresome significant differences in the observed dynamics: the current phase-fieldmodel does not present bipedal cells, but we do observe persistently travelingcells which are not possible in the intermediate system. Most importantly,traveling waves are possible in the phase-field model with symmetric double-wellpotential W . Analysis of (2.10)shows that symmetric potential W (as in (2.7))cannot give rise to traveling wave solutions in the sharp interface limit, andthey did not appear in the intermediate system at all. Moreover, as this systemis a minimal version of the work in [31] there was no evidence that it couldsupport persistently traveling solutions. Thus their existence in the currentmodel is particularly surprising. A more careful quantitative comparison of theintermediate system with the phase-field model would be of particular interest.However, current numerical methods to solve the intermediate system requireweeks to run for a single parameter set, and as such direct comparisons remainunfeasible.Rotating cells were also previously observed in the more sophisticated modelof [24], but it is perhaps surprising to be able to capture them in the currentmodel without heterogeneity of myosin motors or adhesion dynamics. We note12owever that more complex rotating modes, such as two-wave rotating states,were exhibited in [24]. By carefully taking symmetric initial conditions, wewere able to observe such two-wave rotating states, see Figure 5. However,they are unstable, as a perturbation of initial conditions resulted in single-waverotating states. This suggests that heterogeneity of myosin/adhesion is requiredto stabilize such complex states. t=4.6t=4.4t=4.2 t=4.0t=3.8t=3.6 Figure 5: By taking non-random initial conditions ( ρ circular and P = (1 ,
0) on the interiorof the cell), the model exhibits a two-wave state: at t = 3 . t = 4 . t = 4 . ε = .
03 and β = 130. Color indicates value of ρ and arrows represent thevector field P . Finally, our analysis showed that distance traveled (and thus cell velocity) ismaximized for a fixed value of ε , over a range of values of β . Recalling originalparameters, this suggests that a ratio of cell size to surface tension plays a moredominant role in controlling rotating cell speed than actin polymerization andadhesion site formation alone. Previous experiments on Caenorhabditis elegans sperm cells have found that membrane tension can optimize (traveling) cellspeed [6]: heuristically, if membrane tension is too small, polymerization occursin many directions, impeding motion, whereas if tension is larger, polymerizationcan be focused to one direction.We conjecture that by reincorporating heterogeneous myosin motor effectswe can stabilize persistent motion over a wider range of parameter ranges so asto also obtain an optimal velocity occuring during persistent motion.In order to relate simulations to the sharp interface limit equation, we wererestricted to specific scalings related to ε . In future work it would be interestingto consider more general physical parameters in this simple model. Relaxing13hese restrictions would no longer relate the model to the sharp interface limit,though it would certainly allow for a wider range of quantitative biologicalcomparisons.Finally, more rigorous bifurcation analysis between motility modes would bedesirable: we anticipate existence of a Hopf bifurcation occuring from immobileto rotating states, and saddle-node bifurcations occuring between immobile andtraveling wave states.All simulations were completed in MATLAB and run on The College of NewJersey’s high performance cluster, ELSA (Electronic Laboratory for Science andAnalysis). Codes used to generate data can be provided upon reasonable requestto the authors.
5. Acknowledgments
Portions of this research were completed using the high performance comput-ing cluster (ELSA) at the School of Science, The College of New Jersey. Fund-ing of ELSA is provided in part by National Science Foundation OAC-1828163.MSM was additionally supported by a Support of Scholarly Activities Grantat The College of New Jersey. The authors would like to thank Dr. NicholasBattista for helpful discussions during the writing of this manuscript and Dr.Mykhailo Potomkin for insightful discussions on the non-dimensionalization ofthe model. The authors also thank the referees for their helpful feedback whichgreatly improved the exposition. NB additionally would like to thank The Col-lege of New Jersey Mathematics & Statistics Department for its continuingsupport of undergraduate research, helping to propel his career.
Declarations of interest: none
Appendix A. Non-dimensionalization of the phase-field equations
Starting from ∂ t ρ = D ρ ∆ ρ − τ − W (cid:48) ( ρ ) − α ( ∇ ρ ) · P (A.1) ∂ t P = D P ∆ P − τ − P − ζ ∇ ρ, (A.2)introduce non-dimensional variables T and X via the diffusive scaling t = T R D ρ ,and x = XR , where R is a characteristic length scale of the cell, say the cellradius. Then ∂ T ρ = ∆ X ρ − R D ρ τ W (cid:48) ( ρ ) − αRD ρ ( ∇ X ρ ) · P (A.3) ∂ T P = D P D ρ ∆ X P − R D ρ τ P − ζRD ρ ∇ X ρ, (A.4)where ∇ X and ∆ X are derivatives with respect to the scaled variable X . Wedrop X notation for clarity. 14ntroduce the dimensionless parameters ε := (cid:112) D ρ τ R , ˜ α := αRD ρ , ˜ ζ := ζτ R .
By rescaling P → ˜ α P , we then have ∂ T ρ = ∆ ρ − ε W (cid:48) ( ρ ) − ∇ ρ · P (A.5) ∂ T P = D P D ρ ∆ P − τ τ ε P − ˜ α ˜ ζ ε ∇ ρ (A.6)To obtain the model of [7] requires the assumptions D P ∼ εD ρ , τ ∼ ετ , (A.7)and to define β := ˜ α ˜ ζε = αζR D ρ , (A.8)resulting in ∂ T ρ = ∆ ρ − ε W (cid:48) ( ρ ) − ∇ ρ · P (A.9) ∂ T P = ε ∆ P − ε P − β ∇ ρ. (A.10)The scalings D P (cid:28) D ρ and τ (cid:28) τ are consistent with experimental values forkeratocytes reported in [31], see also Table 1.15 eferences [1] I. S. Aranson, Physical models of cell motility , Springer, 2016.[2] E. Barnhart, K.-C. Lee, G. M. Allen, J. A. Theriot, and A. Mogilner,
Balance between cell-substrate adhesion and myosin contraction determinesthe frequency of motility initiation in fish keratocytes , PNAS (2015),5045 – 5050.[3] E. L. Barnhart, J. Allard, S. S. Lou, J. A. Theriot, and A. Mogilner,
Adhesion-dependent wave generation in crawling cells , Current Biology (2017), no. 1, 27–38.[4] E. L. Barnhart, G. M. Allen, F. J¨ulicher, and J. A. Theriot, Bipedal loco-motion in crawling cells , Biophys J. (2010), 933–942.[5] E. L. Barnhart, K. Lee, K. Keren, A. Mogilner, and J. A. Theriot, Anadhesion-dependent switch between mechanisms that determine motile cellshape , PLoS Biol (2011), e1001059.[6] Ellen L Batchelder, Gunther Hollopeter, Cl´ement Campillo, XavierMezanges, Erik M Jorgensen, Pierre Nassoy, Pierre Sens, and Julie Plastino, Membrane tension regulates motility by controlling lamellipodium organiza-tion , Proceedings of the National Academy of Sciences (2011), no. 28,11429–11434.[7] L. Berlyand, M. Potomkin, and V. Rybalko,
Sharp interface limit in a phasefield model of cell motility , Networks & Heterogeneous Media (2017),no. 4, 551.[8] B. A. Camley, Y. Zhao, B. Li, H. Levine, and W.-J. Rappel, Crawling andturning in a minimal reaction-diffusion cell motility model: Coupling cellshape and biochemistry , Phys. Rev. E (2016).[9] Y. Cao, R. Karmakar, E. Ghabache, E. Gutierrez, Y. Zhao, A. Groisman,H. Levine, B. A. Camley, and W.-J. Rappel,
Cell motility dependence onadhesive wetting , Soft matter (2019), no. 9, 2043–2050.[10] R. J. Hawkins, R. Poincloux, O. B´enichou, M. Piel, P. Chavrier, and R. Voi-turiez, Spontaneous contractility-mediated cortical flow generates cell mi-gration in three-dimensional environments , Biophysical journal (2011),1041–1045.[11] Marc Herant and Micah Dembo,
Form and function in cell motility: fromfibroblasts to keratocytes , Biophysical journal (2010), no. 8, 1408–1417.[12] S. Hirsch, A. Manhart, and C. Schmeiser, Mathematical modeling of myosininduced bistability of lamellipodial fragments , Journal of mathematical bi-ology (2017), no. 1-2, 1–22. 1613] J. L¨ober, F. Ziebert, and I. S. Aranson, Modeling crawling cell movementon soft engineered substrates , Soft Matter (2014), 1365–1373.[14] , Collisions of deformable cells lead to collective migration , ScientificReports (2015).[15] S. S. Lou, A. Diz-Mu noz, O. D. Weiner, D. A. Fletcher, and J. A. The-riot, Myosin light chain kinase regulates cell polarization independently ofmembrane tension or rho kinase , J. Cell Biol. (2015), 275–288.[16] M. S. Mizuhara, L. Berlyand, and I. S. Aranson,
Minimal model of directedcell motility on patterned substrates , Physical Review E (2017), no. 5,052408.[17] M. S. Mizuhara, L. Berlyand, V. Rybalko, and L. Zhang, On an evolutionequation in a cell motility model , Physica D (2016), 12–25.[18] M. S. Mizuhara and P. Zhang,
Uniqueness and traveling waves in a cellmotility model , Discrete & Continuous Dynamical Systems-B (2019),no. 6, 2811.[19] A. Mogilner, Mathematics of cell motility: have we got its number? , J.Math. Biol. (2009), 105–134.[20] A. Mogilner and K. Keren, The shape of motile cells , Curr. Biol. (2009),R762–R771.[21] A. Mogilner and A. Manhart, Intracellular fluid mechanics: Coupling cyto-plasmic flow with active cytoskeletal gel , Annual Review of Fluid Mechanics (2018).[22] R. Pierre, P. Thibaut, and L. Truskinovsky, Mechanics of motility initationand motility arrest in crawling cells , Journal of the Mechanics and Physicsof Solids (2015), 469 – 505.[23] F. Raynaud, M. E. Amb¨uhl, C. Gabella, A. Bornert, I. F. Sbalzarini, J. J.Meister, and A. B. Verkhovsky, Minimal model for spontaneous cell polar-ization and edge activity in oscillating, rotating and migrating cells , NaturePhysics (2016), no. 4, 367–373.[24] C. Reeves, B. Winkler, F. Ziebert, and I. S. Aranson, Rotating lamel-lipodium waves in polarizing cells , Communications Physics (2018), no. 1,1–11.[25] D. Shao, W.-J. Rappel, and H. Levine, Computational model for cell mor-phodynamics , Physical review letters (2010), no. 10, 108104.[26] E. Tjhung, A. Tiribocchi, D. Marenduzzo, and M. E. Cates,
A minimalphysical model captures the shapes of crawling cells , Nature communications (2015), no. 1, 1–9. 1727] A. B. Verkhovsky, T. M. Svitkina, and G. G. Borisy, Self-polarization anddirectional motility of cytoplasm , Current Biology (1999), no. 1, 11–S1.[28] A. K. Wilson, G. Gorgas, W. D. Claypool, and P. De Lanerolle, An increaseor a decrease in myosin ii phosphorylation inhibits macrophage motility. ,The Journal of Cell Biology (1991), no. 2, 277–283.[29] B. Winkler, I. S. Aranson, and F. Ziebert,
Confinement and substrate topog-raphy control cell migration in a 3d computational model , CommunicationsPhysics (2019), no. 1, 1–11.[30] F. Ziebert and I. S. Aranson, Effects of adhesion dynamics and substratecompliance on the shape and motility of crawling cells , PLoS ONE (2013),e64511.[31] F. Ziebert, S. Swaminathan, and I. S. Aranson, Model for self-polarizationand motility of keratocyte fragments , J. R. Soc. Interface (2012), 1084–1092.[32] Falko Ziebert and Igor S Aranson, Computational approaches to substrate-based cell motility , npj Computational Materials2