Interplay Between Mechanosensitive Adhesions and Membrane Tension Regulates Cell Motility

The initiation of directional cell motion requires symmetry breaking that can happen with or without external stimuli. During cell crawling, forces generated by the cytoskeleton and their transmission through mechanosen-sitive adhesions to the extracellular substrate play a crucial role. In a recently proposed one-dimensional model [P. Sens, Proc. Natl. Acad. Sci. USA 117 , 24670 (2020)], a mechanical feedback loop between force-sensitive adhesions and cell tension was shown to be sufﬁcient to explain spontaneous symmetry breaking and multiple motility patterns through stick-slip dynamics, without the need to account for signaling networks or active polar gels. We extend this model to two dimensions to study the interplay between cell shape and mechanics during crawling. Through a local force balance along a deformable boundary, we show that the membrane tension coupled with shape change can regulate the spatiotemporal evolution of the stochastic binding of mechanosensitive adhesions. Linear stability analysis identiﬁes the unstable parameter regimes where spontaneous symmetry breaking can take place. Using simulations to solve the fully coupled nonlinear system of equations, we show that, starting from a randomly perturbed circular shape, this instability can lead to keratocyte-like shapes. Simulations predict that different adhesion kinetics and membrane tension can result in different cell motility modes including gliding, zigzag, rotating, and sometimes chaotic movements. Thus, using a minimal model of cell motility, we identify that the interplay between adhesions and tension can select emergent motility modes. DOI: 10.1103/PRXLife.1.023007


I. INTRODUCTION
The mechanism of cell crawling on substrates is important for understanding numerous biological processes such as morphogenesis and wound healing [1].Experiments have revealed that several main subcellular processes, including actin polymerization [2], adhesion [3], and myosin contraction [4], are spatially and temporally orchestrated to generate coherent cellular motion.These experiments have also lent themselves to systematic theoretical and computational modeling [5].Among these different subprocesses, the initiation of motion is of particular interest because it can happen both due to external cues and spontaneously due to intrinsic biochemical or mechanical instabilities, in a process known as cell self-polarization [6,7].A polarized cell undergoes distinct molecular processes at the front and rear, such as the distribution of Rho family GTPases which regulate the actin protrusion and adhesion formation [8]; this distribution specifies a direction for motility.Interestingly, despite the vast number of molecular players involved [8,9], cell migration is essentially a mechanical process [10,11].The integration of cell signaling into mechanical processes has led to different scales of biophysical models [11][12][13].
Depending on the cell type, migrating cells can assume different shapes.For example, fibroblasts have multiple Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
protrusions [14][15][16], while fast-moving keratocytes are characterized by their flat, smooth, fan-shaped leading edges [17].Despite these differences, a common mechanism for cells to undergo directional motion can be summarized as follows: the leading edge of the cell protrudes as a result of the speed difference between actin polymerization and retrograde flow, the rear retracts to keep up with the front [4,18], and membrane tension plays a crucial role in this process because it can coordinate the protrusion and retraction as a global regulator for cell shape change and motility [19][20][21].With this picture, many phenomenological models have been proposed to explain the underlying mechanisms for motility and shape determination.These include models constructed based on the graded radial extension hypothesis [22], viscoelastic actin network and myosin transport [23], force balance between treadmilling actin filaments and membrane tension [17,24], two-phase fluids with actin polymerization [25], and many other redundant mechanisms [26].Later, more comprehensive free-boundary models incorporated other features such as the discrete stick-slip adhesions [27,28], the orientational order of the actin filament network [29][30][31][32], and the feedback loop between actin flow, myosin, and adhesion [33] and are reviewed in Ref. [34].Most of these simulations either start with a crescent shape to match experimental observations [22,23], with perturbations or polarity fields along a specified direction [27][28][29][30][31][32][33], or with a prescribed front [25].Therefore, how the cell shape transitions from random fluctuations to persistent motile shapes remains unclear from these models.
Apart from steady-moving states, keratocytes can also undergo more complex motility modes such as bipedal motion and spontaneous turning, or a combination of both [35][36][37][38].Experimental and modeling studies reveal that the stick-slip adhesive sites at the rear and their coupling with the cytoskeleton dynamics are crucial for such unsteady motions [37][38][39].However, the models used to investigate this coupling need prior knowledge of the positions of adhesion sites and the broken symmetries.As a result, they cannot predict how the distributions of adhesions are formed in the first place.Other simulations that obtain bipedal motion without prescribed adhesion sites take the substrate deformations into account [40,41], but treat adhesion dynamics phenomenologically.
Mogilner et al. suggested that cell polarization and turning may share similar mechanisms that involve the feedback between actomyosin flows and stick-slip dynamics of adhesions [5].Recently, Sens proposed a one-dimensional (1D) mechanical feedback loop between the binding and unbinding of cell-substrate adhesions and linear cell tension that can lead to the spontaneous symmetry breaking [10].Based on a mean-field approximation of the molecular clutch model for adhesions, this theoretical model also predicted various one-dimensional cell locomotion behaviors such as steady crawling, bistability, and bipedal motion as observed in experiments [37,38], bridging the gap between microscopic factors and whole cell locomotion.Combining the mechanosensitive unbinding of adhesions [10] that are based on first principles with cell shape change can provide a natural description for the formation and the stick-slip dynamics of adhesions and how they can lead to unsteady motions.However, this model does not explore the link between the mechanical feedback loop and cell shape and motility modes.Such an exploration requires a two-dimensional (2D) formulation with a deformable boundary.Here, we asked whether the mechanical feedback loop between adhesion and cell shape change is sufficient to explain the dynamics of cell shape change and motility modes.Specifically, we extend the 1D model in Ref. [10] to two dimensions with deformable boundaries to investigate how the coupling between membrane tension and the stick-slip dynamics of adhesions determines cell shape change and the spatial distribution of traction forces resulting in sustained motion.Linear stability analysis of the model shows that uniform steady states are unstable in certain parameter regimes due to the stick-slip nature of the adhesions.Numerical simulations in two dimensions predict that the force balance between the membrane tension, mechanosensitive adhesions, and the contractility near the cell edge is sufficient for a circular-shaped cell to spontaneously initiate and sustain motion [Fig.1(a)], and the various motility modes mentioned above can be captured [Fig.1(b)] even without the consideration of the complex reorganization of the cytoskeleton.Thus, our 2D model is able to capture the initiation of cell polarization and the different motility modes (directional, bipedal, and turning) by tuning the physical parameters, and demonstrates the crucial role for interplay between adhesions, contractility, and membrane tension in cell motility.
At cellular length scales, inertia is negligible [12], so that forces are balanced everywhere at any instant of time.When the actin filaments treadmill at the edge of the cell, the membrane imposes an opposing force while actomyosin contraction generates a contractile force [Fig.1(c)].These two forces lead to the retrograde flow of actin filaments away from the cell edge and are locally balanced by a friction force pointing outwards and created by transmembrane adhesions that can stochastically bind or unbind to the flat substrate [42]: Here, we apply a coarse-grained approach and assume the focal adhesions are concentrated in a narrow region near the cell periphery with width l 1 , as proposed in Ref. [10].
The cell boundary is thus treated as a one-dimensional curve (t ) evolving in a two-dimensional plane (see Fig. 1 of the Supplemental Material [43]), where the adhesion clusters are material points containing a collection of ρl 1 adhesive linkers on average per unit length along the curve.The position vector of a material point on the boundary at time t is represented by x(α, t ), where α ∈ [0, 2π ] is a Lagrangian parameter (a detailed description of the geometry and parametrization can be found in the Supplemental Material [43]).At each material point, the total number of adhesive linkers for each adhesion cluster is assumed to be conserved [42].Each of the adhesive linkers can unbind and rebind between the substrate and the sliding actin filaments with an off rate k off and an on rate k on [Fig.1(d)], and the effect of thermal fluctuations is modeled with an effective diffusivity D. Denoting the current arclength as s(α, t ) ∈ [0, L(t )], the fraction of bound linkers n(α, t ) at each material point evolves by the kinetic equation According to the Bell-Evans formula [44,45], the mechanical force f b felt by a given linker will lower the energy barrier for it to unbind from the substrate, such that the off rate increases exponentially with the force as where k 0 off is the off rate under zero force and f 0 is a molecular force scale for the linker to rupture with a typical order of several piconewtons.For simplicity, the on rate k on is assumed to be a force-independent constant.
The mechanosensitivity of the adhesions allows us to capture the biphasic relation between the friction force and the retrograde velocity v r [46][47][48][49].Here, we adopt a minimal mean-field approximation [10] where the average extension of a linker is approximated by the retrograde velocity v r times its average lifetime 1/k off .Each linker is viewed as an elastic spring with spring constant k b .By Hooke's law, the force experienced by a single linker is given by f b = k b v r /k off .Hence, the retrograde velocity and the dimensionless off rate r = k off /k 0 off can be related by v r = v β r log r, where v β = k 0 off f 0 /k b is a characteristic velocity scale related to the mechanosensitive unbinding process.In addition to the friction generated by the mechanosensitive adhesions, viscous dissipation between the actin flow and the substrate contributes to a linear friction ζ 0 v r .Assuming that protrusions as well as the retrograde flow are locally normal to the cell boundary [50,51], the total friction force is given by where n is a unit outward normal vector on the cell edge (t ) and ζ 1 = ρl 1 k b /k 0 off .To highlight the interplay between the mechanosensitive adhesions and the membrane tension, we simply treat the contractile force per unit length as a constant force pointing inwards along the normal direction, F contraction = −σ c n.The timescale for the rebinding and unbinding of adhesions (seconds) [52] is much longer than the timescale for the membrane force to equilibrate (milliseconds) [17], so the membrane tension, σ m , is assumed to be spatially uniform along the cell boundary.The force per unit length integrated over the lamellipodium height is then given by F membrane = −2hσ m Hn, where H = 1/h + κ is the total curvature given by the in-plane curvature κ and the lamellipodium radius h along the vertical direction to the substrate.
We further assume that the membrane tension depends linearly on the total area A(t ) of the cell as it evolves in time, where k σ is an effective stiffness that accounts for the extensibility of the membrane.Finally, the shape of the cell evolves according to the kinematic boundary condition where the normal velocity is determined by the difference between the actin polymerization velocity v p and the retrograde velocity v r .Here, we treat the polymerization velocity as a constant along the cell edge following Refs.[53,54], although it can be a function of the myosin density as hypothesized in some more detailed modeling approaches [51].Note that our model neglects the dynamics of the cell cytoskeleton and the cytoplasmic flows occurring inside the bulk, as the relevant forces are assumed to be balanced inside and the retrograde flows tend to concentrate near the cell periphery [33].

B. Nondimensionalization
We nondimensionalize the governing equations with the following scales.The characteristic timescale is given by the off rate under zero force load 1/k 0 off .The mechanosensitive unbinding process provides a characteristic velocity scale v β as mentioned in the previous section.The cell size is characterized by its average radius R 0 .The dimensionless variables are listed as follows: The dimensionless parameters governing the system can be classified into two groups, namely, parameters that characterize the stick-slip dynamics of the adhesions and parameters that characterize the general cell properties.We estimate their TABLE I. Range of dimensionless model parameters.See the Supplemental Material [43] for dimensional parameter values [10,23,33,42,47,52,[55][56][57][58][59][60][61].

Dimensionless parameter
Approximate value orders of magnitude from experiments and previous modeling approaches (see the Supplemental Material [43]), as summarized in Table I.The first group contains the dimensionless relative adhesion strength ζ * 1 = ζ 1 /ζ 0 that compares the contribution of the adhesive friction (sticking) and the viscous friction (slipping), the dimensionless on rate r on = k on /k 0 off , and the dimensionless diffusion coefficient The second group contains the dimensionless lamellipodium height h * = h/R 0 , the dimensionless length = v β /k 0 off R 0 that compares the length scale provided by the off rate and the retrograde velocity with the cell size, the dimensionless effective stiffness k * σ = 2k σ R 0 /ζ 0 k 0 off that compares the cell stiffness to the adhesive friction, the dimensionless actomyosin contraction σ * c = σ c /ζ 0 v β , and the dimensionless polymerization velocity v * p = v p /v β .The dimensionless governing equations can be summarized as follows, where the stars have been dropped for simplicity:

C. Numerical implementation
To describe the geometry, the tangent angle, θ , and the arclength derivative, s α , are introduced as independent variables instead of x [62].By specifying the arclength derivative as s α = L/2π , the mesh points are kept equally spaced in arclength at every time step.The derivatives with respect to s and α can be exchanged through ∂ α = s α ∂ s : this enables us to apply a finite difference scheme in the fixed α-parametric domain, to which the curve is mapped as a circle and uniformly discretized in α as α i = 2π (i − 1)/N, i = 1, . . ., N + 1.The method was validated by comparison to the linear stability results at short times (see Fig. 2 of the Supplemental Material [43]), and by comparison with a different scheme based on spline interpolation; excellent agreement was found between the two schemes, with the θ -L method providing enhanced computational speed.
All simulations start from a circular configuration with a uniform fraction of bound linkers, which is perturbed initially, along with other corresponding physical quantities, using the first 100 Fourier modes with amplitudes varying randomly from −10 −4 to 10 −4 .Given the geometric and physical variables at time step t n , we update their values at t n+1 by the following steps: (1) Update the shape of the curve x with the θ -L formulation using an explicit Euler scheme.Compute other geometric quantities such as the normal vector and the curvature.
(2) Update the membrane tension σ m with the explicit Euler scheme.
(3) Update the fraction of bound linkers n with a Crank-Nicolson scheme.
(4) Update the off rate by solving the force balance (10) iteratively using Newton's method and obtain the normal velocity.
The spatial derivatives are discretized by central finite difference and the integrals are computed by the trapezoidal rule with end correction.The number of grid points is taken as N = 2000 and the time step is set to t = 10 −3 .

A. Initiation of cell motion through a stick-slip instability
We first analyze the linear stability of the model to uncover a mechanism for cell spontaneous symmetry breaking and motility initiation through an instability arising from the stick-slip dynamics of the mechanosensitive adhesions.A similar instability was discussed in Ref. [10] for the 1D case as a "stick-slip instability," which gave rise to persistent oscillations between protrusion and retraction phases by switching between sticking and slipping states.Here, we show that, when cell shape change is considered, this instability still exists but distinct modes of deformation are subject to distinct instability criteria due to the spatial effects.
We take the base state for the analysis to be a stationary circle with radius R = 1 where the retrograde velocity balances with the protrusion velocity v p = r log r everywhere, and the fraction of bound linkers is uniformly distributed along the edge with the steady-state value n = r on /(r on + r).Henceforth, overbars are used to denote base-state variables.From the force balance (10), the base-state membrane tension then satisfies We consider small perturbations around the base state with the ansatz φ = φ Here φ represents the base-state value of variable φ, and φ k and λ k denote the initial magnitude of perturbation and corresponding dimensionless growth rate of the kth normal mode, respectively.Note that, in our model, the membrane tension is assumed to be uniform along the cell edge, so perturbations of the membrane tension with modes k = 0 are set to zero.Inserting the ansatz into the governing equations and linearizing the system for small perturbations yields eigenvalue problems for the growth rates λ k .For k = 0, and for k = 0, The governing equations are nonlinearly coupled and, consequently, they provide constraints on the initial perturbations (see the Supplemental Material [43] for details).These constraints satisfy a quadratic equation for the perturbation amplitudes with two possible conjugate imaginary solutions, corresponding to two independent modes of perturbation, analogous to the one-dimensional case [10] where the two independent modes are the symmetric and the antisymmetric modes.
The various Fourier modes represent different modes of deformation.Typical dispersion relations for parameter choices relevant to physiological condition values are plotted in Fig. 2. We find that the precise choice of parameters does not affect the qualitative behavior of the dispersion relation, which can be summarized as follows.The mode k = 0 describes a spatially homogeneous perturbation and corresponds to the global dilation or contraction of the cell.This is the only mode to be affected by the effective membrane stiffness k σ , and we find it to be always linearly stable at physiological values of the stiffness (k σ ∼ 10 3 ).Thus, membrane elasticity always acts to maintain the cell area at its base value in the linear regime.The mode k = 1 is the only mode that is not center symmetric and captures translational motion of the center of mass; its growth rate is found to be purely real regardless of the choice of model parameters.Interestingly, that growth rate is independent of the membrane tension, a consequence of the fact that the k = 1 mode, which describes translation of the center of mass, does not affect the local curvature at linear order (see Eq. ( 3) of the Supplemental Material [43]).When this mode becomes unstable, the corresponding perturbation is amplified and leads to a symmetry breaking in space that singles out a certain direction to initiate locomotion.Note that, different from the translational invariance mentioned in Ref. [54], the cell shape in our model is always coupled to the distribution of the fraction of bound linkers.As a result, a perturbation of the k = 1 mode not only denotes a rigid translocation of the geometry, but also induces a front-rear gradient of the forcesensitive adhesions.This mechanism is similar to the global polarization-translation coupling mode proposed by Lavi et al. [63].Subsequent modes with k > 1 describe various shape deformations with increasingly shorter wavelengths.At relatively low wave numbers, the dispersion relation typically displays two real positive growth rates.As k increases, these give way to two complex conjugate growth rates, suggesting that there can exist oscillations and traveling waves propagating along the edge, known as "stick-slip" waves, which are similar to the lateral waves predicted along a flat edge [10].Finally, all modes beyond a certain wave number become stable, indicating that high-frequency perturbations will decay.
The dependence of the growth rates on the relevant system parameters can also be gleaned from Fig. 2. We first note that a Hopf bifurcation can occur for modes with high wave numbers as the parameters vary, where the real growth rates become growth rates with conjugate imaginary parts, similar to the 1D case [10].We focus on the role of the three parameters that characterize adhesion kinetics.The relative adhesion strength ζ 1 and the on rate r on are found to affect both the magnitude of the unstable growth rates as well as the number of unstable modes, while the effect of the diffusion coefficient D is to damp high-order modes.As ζ 1 increases, both the magnitude of the positive real growth rates and the number of unstable modes increase.The effect of the on rate r on , on the other hand, is nonmonotonic and exhibits a biphasic behavior.When the on rate increases, the magnitudes of unstable growth rates first increase and then decrease.This suggests the possibility of a stick-slip instability similar to the 1D case [10], and provides a mechanism for the cell to oscillate between protruding (sticking) and retracting (slipping) states at the two distinct ends.Here, in two dimensions, certain modes of perturbation are amplified by the mechanosensitive nature of the adhesion clusters and their coupling with the retrograde flow, leading to the local adhesion sites switching between sticking and slipping motions in distinct regions.To make the biphasic relationship between friction and retrograde flow possible, the relative adhesion strength needs to be strong enough for the adhesive friction to play a role compared to the linear friction.Moreover, the rebinding rate cannot be either too small or too large, since in both cases the adhesive linkers either dissociate or rebind to the substrate too quickly.As a result, there is no possible stick-slip transition and the system remains linearly stable.That optimal locomotion efficiency occurs for intermediate adhesions was demonstrated in experiments by varying the concentration of integrins and blocking integrins [64].
Our model assumed a uniform contractile force along the cell edge as a way to single out the role of mechanosensitive adhesions in driving symmetry breaking.The effect of actomyosin contraction on stability is shown in Fig. 2(d), where it is found to have a weak but destabilizing effect on the system.This finding is qualitatively consistent with previous contraction-driven cell motility models and with various experimental observations [33,65,66].

B. Instability of mode k = 1 leads to the formation of the front and rear in the early stage of locomotion
As already mentioned above, the instability of the k = 1 mode provides a mechanism for the self-polarization of the cells and initiation of motility.This is further demonstrated in our numerical simulations.A typical temporal evolution of a cell shape at short times is illustrated in Fig. 3(a), where the cell edge is colored by the fraction of bound linkers.Starting from a circular shape and a randomly perturbed initial distribution of bound linkers at t = 0, high-frequency fluctuations are found to decay very rapidly (t = 0.5) as predicted by the stability analysis.The k = 1 mode, which has the largest growth rate, grows simultaneously, resulting in the formation of a potential cell front where more linkers are bound to the substrate and of a potential rear where fewer linkers are bound (t = 1).The polymerization velocity then exceeds the retrograde velocity in the front while the opposite occurs at the rear, leading to expansion in the front, shrinkage at the rear, and translocation of the center of mass (t = 1.5).The front further keeps protruding while the rear keeps retracting, and the cell ultimately evolves to a fanlike shape with a smooth leading edge (t = 2) very similar to the characteristic shape of coherent keratocytes [37].
The spatiotemporal evolutions of the fraction of bound linkers, off rate, friction force, and curvature during these early stages are plotted as kymographs in Figs.3(b)-3(e).As a result of the growth of the k = 1 mode, the retrograde velocity is low in the front of the cell and high at the rear.Consequently, strong adhesion develops in the front where the off rate is low and hence most linkers are bound to the substrate.Meanwhile the off rate is high with most linkers unbound along regions at the back and sides, similar to the traction stress distribution revealed by experiments in migrating cells [67].Comparing the spatiotemporal distribution of the friction and the curvature [Figs.3(d) and 3(e)] shows that regions with high friction coincide with regions with high curvature as required by the force balance equation.This is also consistent with experiments [67] where traction stress concentrates at the sides.Note that, unlike many previous modeling approaches [37] where the breaking of front-rear symmetry and adhesion distribution were treated as inputs to trigger the initiation of motion and turning, our model provides a route for this distribution of adhesions to emerge spontaneously.

C. The stability diagram predicts various motility modes
Next, we turn our focus to long-time dynamics and the effect of adhesion kinetic parameters.While the initiation of locomotion through the instability of the k = 1 mode is generic, various cell motility modes are observed at long times depending on the relative adhesion strength ζ 1 and on rate r on .We categorize them in a phase diagram in Fig. 4(a), in which the blue shaded regions highlight the linearly unstable regimes for the Fourier modes with wave numbers k = 1, 2, 5, 10, 15, while the dots denote the long-time motility modes.Each dot PRX LIFE 1, 023007 (2023) represents five simulations with different random initial perturbations.A systematic exploration of the parameter space allowed us to classify the motility modes by the following criteria: (i) if the cell moves in a straight line with no turns over the entire simulation time range (t = 100), we call it a "gliding" mode [Figs.4(b) and 4(c)]; (ii) if the cell turns only towards one direction, we call it a "rotational" mode [Fig.4(d)]; (iii) if the cell turns alternatively between left and right resulting in a lateral oscillation, we call it a "zigzag" mode [Figs.4(e) and 4(f)]; (iv) for some combinations of parameters and for different initial conditions, we observed "mixed" modes where the cell can exhibit different types of motions including the three cases mentioned above; and (v) finally, there exist more complicated types of motion that do not fit into any of the cases above and that we label as "others" [Figs.4(g)-4(j)].Most of the motility modes described here have been observed in experiments as well as in previous models: the gliding motion is well known in fast-moving fish keratocytes [68], the zigzag motion is analogous to the "bipedal" motion of oscillating keratocytes [37,38], and spontaneous turning has also been observed in keratocytes even in the absence of external cues [35,39].
Although the dynamics at long times are essentially nonlinear, Fig. 4(a) shows a clear correlation between the stability diagram and the motility modes.In regions of the parameter space with fewer unstable modes, where the on rate is either relatively low or high, the motion tends to be uniaxial with the cell remaining polarized with left-right symmetry.More complex dynamics arise for intermediate on rates, which we attribute to the stick-slip instability.Rotational and zigzag modes arise in this regime, and are characterized by broken left-right symmetry in addition to front-rear asymmetry, resulting in quasiperiodic turns in the cell trajectory.We take a closer look at these turns in Figs.5(a) and 5(b), where we plot the temporal evolution of the membrane tension together with the corresponding shapes and center-of-mass velocities within one oscillation period for the rotational and zigzag modes, respectively.At the initiation of a turn, the distribution of bound linkers is first observed to become asymmetric, with an increase in the fraction of bound linkers on one side of the cell.This results in enhanced sticking on that side along with slipping on the opposite side, allowing the cell to turn in that direction, and this mechanism is reminiscent of the "sticking wave" predicted in the 1D flat case [10].As the cell shape becomes asymmetric and the cell rotates, the membrane tension is found to increase and reaches a peak value before relaxing again in the later stage of the turn.In the zigzag case, the high adhesion region alternatively switches between the left and right sides of the cell, while for the rotational cases the direction in which the high adhesion regions form remains fixed, with equal probabilities for a cell to turn clockwise or counterclockwise.The Hopf bifurcation of higher-order modes may be responsible for this transition from uniaxial translation to oscillatory motion.Moreover, we find that there is a strong correlation between cell shapes and motility modes.For the gliding cases, the cell adopts either a triangular-like shape [for small on rates, Fig. 4(b)] or a nearly circular shape [for large on rates, Fig. 4(c servation is that, regardless of the randomness in the initial perturbations, the emerging shapes corresponding to a given set of parameters are robust; we further elaborate on this point below.

D. Robustness of motility modes
The observations above suggest that the feedback loop between mechanosensitive adhesions and membrane tension uniquely determines the emergent motility mode for a given set of system parameters, irrespective of the initial condition or history of the system.To further demonstrate the robustness of these modes, we vary the parameters during a simulation to analyze transitions between modes.A typical transition is shown in Fig. 6(a): starting from a cell performing a gliding motion, we abruptly vary the value of r on from 16 to 30 at t 0 = 25, causing it to switch to a zigzag mode.As shown in the snapshots, the cell glides smoothly at first, and, after the change in r on , its shape quickly adjusts and starts undergoing zigzags.The oscillation frequency and cell morphology after the transition are similar to the corresponding zigzag motions with r on = 30 starting from a randomly perturbed initial condition.The time instant t 0 at which we start altering the parameters, the period of time t within which we gradually alter the parameters, and the intermediate states do not significantly affect the zigzag dynamics [Fig.6(b)].To quantitatively compare the resulting zigzag modes emerging through different routes, we plot the temporal evolution of the membrane tension in Fig. 6(b) for different choices of t 0 and t.After increasing the on rate and the adhesion strength, the cell expands, leading to an increase in the membrane tension towards the tension value of the new state, regardless of the history of states.This is in agreement with the one-dimensional case [10] as well as experiments [19] where the membrane tension in motile cells is determined by the adhesion strength and cytoskeletal forces, and increases as cell-substrate adhesion strengthens-in our model, as either the on rate or relative adhesion strength increases.
Similar transitions are also observed when we change parameters between other modes (see the Supplemental Material [43], Movies S10-S12).In all cases, after changing the parameters, the cell first goes through a transient state and then soon evolves into the motility mode selected by the new set of parameters in the phase diagram of Fig. 4.This further suggests that these motility modes are attractors of the dynamical system, with the model solutions approaching either fixed points or stable limit cycles in phase space.

E. Effect of adhesion parameters
Finally, we investigate the effect of varying adhesion kinetic parameters on cell geometry, mechanics, and locomotion.For various combinations of r on and ζ 1 , we calculate the time-averaged circularity 4π A/L 2 (with a maximum value of 1 corresponding to a perfect circle), the average membrane tension, and the average distance of the center of mass from the origin x c = M i=1 |x (i) c − x (0) c |/M, and plot them as phase diagrams in Fig. 7 where the motility modes are also labeled.Note that the precise boundary between each mode changes slightly as we vary the diffusion coefficient D, but the qualitative behavior remains unchanged.For the gliding modes with low on rates, the cell adopts a nearly triangular shape with low circularity, with only a weak dependence on ζ 1 .As the on rate increases, the shapes all become more circular regardless of motility mode as seen in Fig. 7(a).A possible explanation for this behavior is that, as the on rate increases, the difference in the fraction of bound linkers between the front and the rear decreases, and therefore the difference in the sticking and slipping velocity also decreases, resulting in less-deformed shapes.The average membrane tension shows little correlation with the various motility modes in Fig. 7(b), and is mainly determined by the on rate and the relative adhesion strength, consistent with the discussion above [10,19].
The ability of a cell to explore space depends strongly on its motility mode as illustrated in Fig. 7(c), where we show the phase diagram for the average distance traveled by the cell.The gliding and the zigzag modes are more unidirectional and thus allow the cell to travel for a longer distance than in the rotational mode.In the case of zigzag trajectories, the left-right oscillations can be accompanied by a circular motion when the relative adhesion strength and the on rate increase, leading to a decrease in the distance traveled.Ultimately, the zigzag mode gives way to trajectories labeled as "others" in Fig. 7(c), which are characterized by irregular turns and may result from the superposition of multiple oscillatory periods.

IV. DISCUSSION
In this work, we have extended the one-dimensional model coupling actin polymerization, adhesion, membrane tension, and shape change previously proposed by Sens [10] to two dimensions to elucidate the role of friction originating from adhesion kinetics in the initiation of migration, the determination of cell shape, and the selection of motility modes through its nonlinear coupling to the membrane tension.With this minimal model, we showed that motile cells can display rich dynamical behaviors by relying on a relatively simple set of physical mechanisms and couplings.We first performed a linear stability analysis and demonstrated that the k = 1 mode, describing translocation of the center of mass, is always the most unstable and is responsible for the cell to obtain its polarity and identify a direction in space for migration by spontaneous symmetry breaking.Our linear analysis of this process [Eq.( 13)] suggests that the initiation of polarization is driven by the interplay of the adhesion kinetics and retrograde flow.We then conducted nonlinear numerical simulations and showed that, driven by a constant actin polymerization rate, the model was able to capture various motility modes commonly seen in experiments, such as unidirectional gliding, bipedal motion, and turning.We identified the nonlinear coupling between stochastic adhesion and membrane tension as the key mechanism involved in the selection of motility modes.Specifically, the relative adhesion strength and on rate were shown to govern membrane tension, which in turn connects spatially distributed adhesions and thus couples the retrograde flow with shape change and adhesion kinetics.
The basic physical mechanisms involved in this process are summarized in Fig. 8.A local increase in the off rate r can have two distinct effects: a decrease in the fraction of bound linkers n and a local retraction of the cell edge.The latter will lead to a global decrease in cell area and, consequently, a decrease in membrane tension due to the membrane elasticity.The force balance (10) dictates that the decrease in n will cause an increase in the off rate r, whereas the decrease in membrane tension will cause r to decrease.Therefore, there is a competition between the membrane tension and the fraction of bound linkers in the regulation of the off rate, tuned by the relative adhesion strength ζ 1 .When ζ 1 is large, the effect of n dominates, resulting in a positive feedback on r, rendering the system unstable.On the other hand, when ζ 1 is small, the situation is the opposite and the system is stable.These predictions are borne out by the results of our stability analysis and numerical simulations.Though membrane tension does not affect the stability of the k = 1 mode, it does affect the stability of higher-order modes and was shown to be involved in cell turns and associated shape changes in the nonlinear regime (Fig. 5).
Our model predictions qualitatively recapitulate experimental observations of symmetry breaking and cell motility modes.Using experimental measurements and mechanical models, Barnhart et al. showed that two feedback loops-one between actin flow and adhesions and the other between actin flow and myosins-are required to initiate motility in fish keratocytes [33].Similarly, using a phase-field approach, Shao et al. [28] showed that coupling between adhesions, actin flow, and myosin contraction is required to recapitulate different experimental observations of keratocyte shape change.In the present work, we showed that the "stick-slip" instability initially proposed by Sens in one dimension [10] can lead to spontaneous initiation of motility by symmetry breaking in two dimensions, and that the feedback between adhesion near the cell edge and membrane tension is sufficient to capture the emergence of complex motility modes.Our study demonstrates that cell polarization, as well as nonlinear motility characteristics and cell shape selection, can potentially be explained by sole consideration of the balance of forces and fluxes applied at the cell edge, independently of bulk cytoskeletal dynamics.An interesting and as yet unexplained experimental observation that our model may help shed light on is that, at low temperatures, the trajectories of fast-moving keratocytes tend to be more unidirectional, while their motion is more circular and less persistent at high temperatures [5].Previous studies have proposed a "steering wheel" mechanism [39], whereby the turning of a cell is caused by asymmetrically distributed adhesion sites at the rear, but the relationship between the turning and temperature is unclear.This effect can be explained by the molecular clutch approach for adhesions in our model.According to Bell's theory [44], the reaction rate increases with temperature, so binding and unbinding processes should be more active at high temperatures.By this effect, an increase in temperature drives an increase in the on rate and relative adhesion strength, and thus drives the transition from gliding to zigzag, rotational, and other modes.
While our model captures many experimental observations and makes testable predictions on the interaction between membrane tension and adhesion, it has certain limitations.In particular, it does not account for the complex molecular machinery underlying actomyosin contraction or associated signaling pathways and only contains a single mechanical feedback loop.We also note that even though we obtain fan-shaped cells in the early stages of motion, the shapes ultimately become triangular for gliding cells, which could be a result of the constant actin polymerization velocity assumed in our model and of ignoring actin remodeling events [5,11,13,51].The contraction generated by the distinct distribution of myosin motors in crawling cells is also neglected in our model.In some studies, contraction alone can be shown to generate spontaneous symmetry breaking and motions [65].
Finally, like other cell motility models [27,29,33], we have assumed that membrane tension is spatially uniform, yet spatial variations in tension have been observed in cells in experiments [69].An extension of our approach to incorporate some of these additional details as well as internal active stresses has the potential to further enrich the model.

PRX LIFE 1 ,
FIG. 1.(a) Schematic describing self-polarization in cells.Fluctuations of certain wavelengths are amplified by the intrinsic instability of the feedback loop between adhesions and membrane tension, leading to self-polarization.(b) Schematic describing the various motility modes.Directional motion is the uniaxial motion of the cell in one preferred direction.Bipedal motion involves antiphase retraction of the left-right trailing edge and lateral oscillation of the cell.In turning motion, the cell rotates in a certain direction with left-right asymmetry.(c) Sketch of the two-dimensional model (top view).The cell shape is determined by the difference between the polymerization velocity v p and the retrograde velocity v r , where v r is related to the off rate of adhesive bonds distributed within the lamellipodium with width l 1 along the cell boundary.The force balance between friction, membrane tension, and contraction is maintained at the cell edge.(d) Sketch of the two-dimensional model (side view).The adhesive bonds bind to the substrate with a constant rate k on and unbind with a force-dependent rate k off as the actin filaments move.

FIG. 2 .
FIG. 2. Dispersion relations showing the real and imaginary parts of the growth rate λ k as functions of wave number k for various choices of the parameters related to adhesion kinetics: (a) growth rates under different relative adhesion strengths ζ 1 with r on = 20, D = 0.01, σ c = 100; (b) growth rates under different on rates r on with ζ 1 = 600, D = 0.01, σ c = 100; (c) growth rates under different diffusion coefficients D with ζ 1 = 600, r on = 20, σ c = 100; and (d) growth rates under different contraction strengths σ c with ζ 1 = 600, r on = 20, D = 0.01.In all cases, the parameters related to the general cell properties are fixed as = 0.001, h = 0.01, k σ = 1000, r p = 50.

023007- 5 FIG. 3 .
FIG. 3. Spatiotemporal evolution of a fan-shaped cell during the initiation of motion (t = 0 to t = 2) from a nonlinear numerical simulation.(a) Time evolution of the cell shape, where the edge is colored by the fraction of bound linkers.The two black dots show the initial and current center-of-mass positions.[(b)-(e)] Kymographs of the fraction of bound linkers, n, off rate r, total friction r log r + ζ 1 n log r, and local curvature κ, as functions of normalized arclength s/L and time t.Parameters values: = 0.001, h = 0.01, k σ = 1000, σ c = 100, r p = 50, r on = 25, ζ = 600, D = 0.05.
FIG. 5. Temporal evolution of the membrane tension and corresponding cell shapes and center-of-mass velocities as functions of time during one period of oscillation for the (a) rotational mode with (r on , α 1 ) = (24, 540) and (b) zigzag mode with (r on , α 1 ) = (26, 600).The cell shapes are colored by the local fraction of bound linkers n, and the orange arrows show the instantaneous center-ofmass velocity.In both cases, other parameters are fixed as = 0.001, h = 0.01, k σ = 1000, σ c = 100, r p = 50, D = 0.05.

PRX LIFE 1 ,
FIG. 6. Transition between different motility modes.(a) Snapshots of a gliding cell switching into a zigzag mode as the adhesion parameters (r on , α 1 ) abruptly change from (16,580) to (30,580) at t = 25.Other parameters are fixed as = 0.001, h = 0.01, k σ = 1000, r p = 50, D = 0.05.(b) Temporal evolution of the membrane tension as the parameters (r on , α 1 ) are varied to cause changes in motility modes.Here, t denotes the duration of the time interval over which the parameters are gradually varied.See movies S10-S12 in the Supplemental Material [43] for videos of the dynamics.

FIG. 7 .
FIG. 7. Phase diagrams in the (r on , ζ 1 ) parameter space showing color maps of the time-averaged (a) cell circularity 4π A/L 2 , (b) membrane tension σ m , and (c) average distance of the center of mass from the origin, x c .The diagrams show averages over five simulations of duration T = 95.Symbols indicate the corresponding motility modes previous identified in Fig. 4. Other parameters are fixed as = 0.001, h = 0.01, k σ = 1000, σ c = 100, r p = 50, D = 0.05.