# Cell motility modes are selected by the interplay of mechanosensitive
adhesion and membrane tension

Yuzhu Chen, David Saintillan, Padmini Rangamani

###### Abstract

The initiation of directional cell motion requires symmetry breaking that can
happen both with or without external stimuli. During cell crawling, forces
generated by the cytoskeleton and their transmission through mechanosensitive
adhesions to the extracellular substrate play a crucial role. In a recently
proposed 1D model (Sens, PNAS 2020), a mechanical feedback loop between force-
sensitive adhesions and cell tension was shown to be sufficient 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 extended this model to 2D 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 identified 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.

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 [12, 11, 13].

Depending on the cell type, migrating cells can assume different shapes. For
example, fibroblasts have multiple 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 [24, 17], 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 [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.

Figure 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_{\\_}{\text{on}}$ and unbind with a force-dependent rate
$k_{\\_}{\text{off}}$ as the actin filaments move.

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, 40].
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 [41, 42], 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 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 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 [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
during the initiation of motility and resulting 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 2D predict that the interplay between membrane
tension and mechanosensitive adhesions is sufficient for a circular-shaped
cell to spontaneously initiate and sustain motion (Fig. 1A), and the various
motility modes mentioned above can be captured (Fig. 1B) 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
and membrane tension in cell motility.

## Model Development

### Governing equations

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. 1C). 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 [43]:

$\displaystyle\boldsymbol{F}_{\\_}{\text{contraction}}+\boldsymbol{F}_{\\_}{\text{membrane}}+\boldsymbol{F}_{\\_}{\text{friction}}=\boldsymbol{0}.$
(1)

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 by [10]. The cell boundary is thus treated as a one-dimensional
curve $\Gamma(t)$ evolving in a two-dimensional plane (see Fig. S1 of the SI
Appendix), where the adhesion clusters are material points containing a
collection of $\rho 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 $\boldsymbol{x}(\alpha,t)$, where
$\alpha\in[0,2\pi]$ is a Lagrangian parameter (a detailed description of the
geometry and parameterization can be found in the SI). At each material point,
the total number of adhesive linkers for each adhesion cluster is assumed to
be conserved [43]. Each of the adhesive linkers can unbind and rebind between
the substrate and the sliding actin filaments with an off-rate
$k_{\\_}{\text{off}}$ and an on-rate $k_{\\_}{\text{on}}$ (Fig. 1D), and the
effect of thermal fluctuations is modeled with an effective diffusivity $D$.
Denoting the current arclength as $s(\alpha,t)\in[0,L(t)]$, the fraction of
bound linkers $n(\alpha,t)$ at each material point evolves by the kinetic
equation

$\displaystyle\dfrac{\partial n}{\partial
t}(\alpha,t)=k_{\\_}{\text{on}}(1-n)-k_{\\_}{\text{off}}n+D\dfrac{\partial^{2}n}{\partial
s^{2}}.$ (2)

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
$k_{\\_}{\text{off}}=k_{\\_}{\text{off}}^{0}\text{exp}(f_{\\_}b/f_{\\_}0)$,
where $k_{\\_}{\text{off}}^{0}$ 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 pN. For simplicity, the on-rate $k_{\\_}{\text{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_{\\_}{\text{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_{\\_}bv_{\\_}r/k_{\\_}{\text{off}}$. Hence, the retrograde
velocity and the dimensionless off-rate
$r=k_{\\_}{\text{off}}/k_{\\_}{\text{off}}^{0}$ can be related by
$v_{\\_}r=v_{\\_}\beta r\log r$, where
$v_{\\_}\beta=k_{\\_}{\text{off}}^{0}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 $\zeta_{\\_}0v_{\\_}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

$\displaystyle\boldsymbol{F}_{\\_}{\text{friction}}=\left(\zeta_{\\_}0v_{\\_}r+\zeta_{\\_}1\frac{nv_{\\_}r}{r}\right)\boldsymbol{n},$
(3)

where $\boldsymbol{n}$ is a unit outward normal vector on the cell edge
$\Gamma(t)$.

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,
$\boldsymbol{F}_{\\_}{\text{contraction}}=-\sigma_{\\_}c\boldsymbol{n}$. The
time scale for the rebinding and unbinding of adhesions (seconds) [52] is much
longer than the time scale for the membrane force to equilibrate
(milliseconds) [17], so the membrane tension, $\sigma_{\\_}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
$\boldsymbol{F}_{\\_}{\text{membrane}}=-2h\sigma_{\\_}mH\boldsymbol{n}$, where
$H=1/h+\kappa$ is the total curvature given by the in-plane curvature $\kappa$
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,

$\displaystyle\frac{\mathrm{d}\sigma_{\\_}m}{\mathrm{d}t}(t)=k_{\\_}\sigma\frac{\mathrm{d}A}{\mathrm{d}t},$
(4)

where $k_{\\_}\sigma$ 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

$\displaystyle\dfrac{\partial\boldsymbol{x}}{\partial
t}(\alpha,t)=(v_{\\_}p-v_{\\_}r)\boldsymbol{n},$ (5)

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 [53, 54], although it can be a function of the myosin density as
hypothesized in some more detailed modeling approaches [51].

### Non-dimensionalization

We nondimensionalize the governing equations with the following scales. The
characteristic time scale is given by the off-rate under zero force load
$1/k_{\\_}{\text{off}}^{0}$. The mechanosensitive unbinding process provides a
characteristic velocity scale $v_{\\_}\beta$ as mentioned in the previous
subsection. The cell size is characterized by its average radius $R_{\\_}0$.
The dimensionless variables are listed as follows:

$\displaystyle t^{\ast}=tk_{\\_}{\text{off}}^{0},\quad
v_{\\_}r^{\ast}=\frac{v_{\\_}r}{v_{\\_}\beta},\quad\boldsymbol{x}^{\ast}=\frac{\boldsymbol{x}}{R_{\\_}0},\quad\sigma_{\\_}m^{\ast}=\frac{2\sigma_{\\_}m}{\zeta_{\\_}0v_{\\_}\beta}.$
(6)

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 orders of magnitude from experiments and previous modeling
approaches (see SI), as summarized in Table 1. The first group contains the
dimensionless relative adhesion strength
$\zeta_{\\_}1^{\ast}={\zeta_{\\_}1}/{\zeta_{\\_}0}$ that compares the
contribution of the adhesive friction (sticking) and the viscous friction
(slipping), the dimensionless on-rate
$r_{\\_}{\text{on}}=k_{\\_}{\text{on}}/k_{\\_}{\text{off}}^{0}$, and the
dimensionless diffusion coefficient
$D^{\ast}=D/k_{\\_}{\text{off}}^{0}R_{\\_}0^{2}$. The second group contains
the dimensionless lamellipodium height $h^{\ast}=h/R_{\\_}0$, the
dimensionless length $\epsilon=v_{\\_}\beta/k_{\\_}{\text{off}}^{0}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_{\\_}\sigma^{\ast}=2k_{\\_}\sigma
R_{\\_}0/\zeta_{\\_}0k_{\\_}{\text{off}}^{0}$ that compares the cell stiffness
to the adhesive friction, the dimensionless actomyosin contraction
$\sigma_{\\_}c^{\ast}=\sigma_{\\_}c/\zeta_{\\_}0v_{\\_}\beta$, and the
dimensionless polymerization velocity $v_{\\_}p^{\ast}=v_{\\_}p/v_{\\_}\beta$.

Table 1: Range of dimensionless model parameters. See SI Appendix for dimensional parameter values. Dimensionless parameter | Approximate value
---|---
$\zeta_{\\_}1^{\ast}=\zeta_{\\_}1/\zeta_{\\_}0$ | $400-600$
$r_{\\_}{\text{on}}=k_{\\_}{\text{on}}/k_{\\_}{\text{off}}^{0}$ | $10-60$
$D^{\ast}=D/k_{\\_}{\text{off}}^{0}R_{\\_}0^{2}$ | $0.01-0.1$
$\epsilon=v_{\\_}\beta/k_{\\_}{\text{off}}^{0}R_{\\_}0$ | $0.001$
$h^{\ast}=h/R_{\\_}0$ | $0.01$
$\sigma_{\\_}c^{\ast}=\sigma_{\\_}c/\zeta_{\\_}0v_{\\_}\beta$ | $100$
$k_{\\_}\sigma^{\ast}=2k_{\\_}\sigma R_{\\_}0/\zeta_{\\_}0k_{\\_}{\text{off}}^{0}$ | $1000$
$v_{\\_}p^{\ast}=v_{\\_}p/v_{\\_}\beta$ | $200$

The dimensionless governing equations can be summarized as follows, where the
stars have been dropped for simplicity:

$\displaystyle\dfrac{\partial\boldsymbol{x}}{\partial
t}(\alpha,t)=\epsilon(v_{\\_}p-r\log r)\boldsymbol{n},$ (7)
$\displaystyle\dfrac{\partial n}{\partial
t}(\alpha,t)=r_{\\_}{\text{on}}(1-n)-rn+D\dfrac{\partial^{2}n}{\partial
s^{2}},$ (8)
$\displaystyle\frac{\mathrm{d}\sigma_{\\_}m}{\mathrm{d}t}(t)=k_{\\_}\sigma\frac{\mathrm{d}A}{\mathrm{d}t}(t),$
(9) $\displaystyle\sigma_{\\_}c+\sigma_{\\_}m(1+\kappa h)=r\log
r+\zeta_{\\_}1n\log r.$ (10)

### Numerical implementation

To describe the geometry, the tangent angle, $\theta$, and the arclength
derivative, $s_{\\_}\alpha$, are introduced as independent variables instead
of $\boldsymbol{x}$. By specifying the arclength derivative as
$s_{\\_}\alpha=L/2\pi$, the mesh points are kept equally spaced in arclength
at every time step. The derivatives with respect to $s$ and $\alpha$ can be
exchanged through $\partial_{\\_}\alpha=s_{\\_}\alpha\partial_{\\_}s$: this
enables us to apply a finite difference scheme in the fixed
$\alpha$-parametric domain, to which the curve is mapped as a circle and
uniformly discretized in $\alpha$ as
$\alpha_{\\_}i=2\pi(i-1)/N,\;i=1,\cdots,N+1$. The method was validated by
comparison to the linear stability results at short times (see Fig. S2 of the
SI Appendix), and by comparison with a different scheme based on spline
interpolation; excellent agreement was found between the two schemes, with the
$\theta$–$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 $\boldsymbol{x}$
with the $\theta$–$L$ formulation using an explicit Euler scheme. Compute
other geometric quantities such as the normal vector and the curvature. (2)
Update the membrane tension $\sigma_{\\_}m$ with the explicit Euler scheme.
(3) Update the fraction of bound linkers $n$ with a Crank-Nicholson 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 $\Delta t=10^{-3}$.

## Results

### 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 [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
$\bar{R}=1$ where the retrograde velocity balances with the protrusion
velocity $v_{\\_}p=\bar{r}\log\bar{r}$ everywhere, and the fraction of bound
linkers is uniformly distributed along the edge with the steady-state value
$\bar{n}=r_{\\_}{\text{on}}/(r_{\\_}{\text{on}}+\bar{r})$. Henceforth,
overbars are used to denote base-state variables. From the force balance (10),
the base state surface tension then satisfies

$\displaystyle\bar{\sigma}_{\\_}{m}=\frac{1}{1+h/\bar{R}}(\bar{r}\log\bar{r}+\zeta_{\\_}1\bar{n}\log\bar{r}-\sigma_{\\_}c).$
(11)

We consider small perturbations around the base state with the ansatz
$\phi=\bar{\phi}+\delta\phi=\bar{\phi}+\sum_{\\_}k\phi_{\\_}k\text{exp}(\mathrm{i}k\theta+\lambda_{\\_}kt)$.
Here $\bar{\phi}$ represents the base-state value of variable $\phi$, and
$\phi_{\\_}k$ and $\lambda_{\\_}k$ denote the initial magnitude of
perturbation and corresponding dimensionless growth rate of the
$k^{\mathrm{th}}$ normal mode, respectively. Note that in our model, the
surface tension is assumed to be uniform along the cell edge, so perturbations
of the surface tension with modes $k\neq 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
$\lambda_{\\_}k$. For $k=0$,

$\displaystyle\Big{(}1+\log\bar{r}+\zeta_{\\_}1\frac{\bar{n}}{\bar{r}}\Big{)}\lambda_{\\_}0^{2}+\Big{\\{}\Big{(}1+\log\bar{r}+\zeta_{\\_}1\frac{\bar{n}}{\bar{r}}\Big{)}(r_{\\_}{\text{on}}+\bar{r})$
$\displaystyle+\\!\Big{[}\Big{(}1+\frac{h}{\bar{R}}\Big{)}2\pi\bar{R}k_{\\_}\sigma\\!-\\!\frac{\epsilon
h\bar{\sigma}_{\\_}{m}}{\bar{R}^{2}}\Big{]}(1+\log\bar{r})\\!-\\!\zeta_{\\_}1\bar{n}\log\bar{r}\Big{\\}}\lambda_{\\_}0$
$\displaystyle+\\!\Big{[}\Big{(}1+\frac{h}{\bar{R}}\Big{)}2\pi\bar{R}k_{\\_}\sigma\\!-\\!\frac{\epsilon
h\bar{\sigma}_{\\_}{m}}{\bar{R}^{2}}\Big{]}(1\\!+\\!\log\bar{r})(r_{\\_}{\text{on}}\\!+\\!\bar{r})=0,$
(12)

and for $k\neq 0$,

$\displaystyle\Big{(}1\\!+\\!\log\bar{r}\\!+\\!\zeta_{\\_}1\frac{\bar{n}}{\bar{r}}\Big{)}\lambda_{\\_}k^{2}\\!+\\!\Big{[}\Big{(}1\\!+\\!\log\bar{r}\\!+\\!\zeta_{\\_}1\frac{\bar{n}}{\bar{r}}\Big{)}\Big{(}r_{\\_}{\text{on}}\\!+\\!\bar{r}\\!+\\!\frac{Dk^{2}}{\bar{R}^{2}}\Big{)}$
$\displaystyle+\frac{\epsilon
h\bar{\sigma}_{\\_}{m}}{\bar{R}^{2}}(k^{2}-1)(1+\log\bar{r})-\zeta_{\\_}1\bar{n}\log\bar{r}\Big{]}\lambda_{\\_}k$
$\displaystyle+\frac{\epsilon
h\bar{\sigma}_{\\_}{m}}{\bar{R}^{2}}(k^{2}-1)(1+\log\bar{r})\Big{(}r_{\\_}{\text{on}}+\bar{r}+\frac{Dk^{2}}{\bar{R}^{2}}\Big{)}=0.$
(13)

The governing equations are nonlinearly coupled and, consequently, they
provide constraints on the initial perturbations (see the SI 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 anti-symmetric
modes.

The various Fourier modes represent different modes of deformation. Typical
dispersion relations for parameter choices relevant to physiological
conditions 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_{\\_}\sigma$, and we find it to be always linearly
stable at physiological values of the stiffness ($k_{\\_}\sigma\sim 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. Therefore, 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. Subsequent modes with
$k>1$ describe various shape deformations with increasingly shorter
wavelengths. At relatively low wavenumbers, 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 wavenumber become stable,
indicating that high-frequency perturbations will decay.

Figure 2: Dispersion relations showing the real and imaginary parts of the
growth rate $\lambda_{\\_}k$ as functions of wavenumber $k$ for various
choices of the parameters related to adhesion kinetics: (A) growth rates under
different relative adhesion strengths $\zeta_{\\_}1$ with
$r_{\\_}\text{on}=20,\;D=0.01,\;\sigma_{\\_}c=100$; (B) growth rates under
different on-rates $r_{\\_}\text{on}$ with
$\zeta_{\\_}1=600,\;D=0.01,\;\sigma_{\\_}c=100$; (C) growth rates under
different diffusion coefficients $D$ with
$\zeta_{\\_}1=600,\;r_{\\_}\text{on}=20,\;\sigma_{\\_}c=100$; (D) growth rates
under different contraction strengths $\sigma_{\\_}c$ with
$\zeta_{\\_}1=600,\;r_{\\_}\text{on}=20,\;D=0.01$. In all cases, the
parameters related to the general cell properties are fixed as
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;r_{\\_}p=50$.

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 wavenumbers 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 $\zeta_{\\_}1$ and the on-
rate $r_{\\_}\text{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 $\zeta_{\\_}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_{\\_}\text{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 [55].

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 [56, 33, 57].

### Instability of mode $k=1$ is responsible for 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 fan-like shape with a smooth leading edge ($t=2$) very similar to the
characteristic shape of coherent keratocytes [37].

Figure 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+\zeta_{\\_}1n\log r$, and local curvature
$\kappa$, as functions of normalized arclength $s/L$ and time $t$. Parameters
values:
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;\sigma_{\\_}c=100,\;r_{\\_}p=50,\;r_{\\_}\text{on}=25,\;\zeta=600,\;D=0.05$.

The spatiotemporal evolution of the fraction of bound linkers, off-rate,
friction force and curvature during these early stages are plotted as
kymographs in Fig. 3$B$–$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, in
agreement with the traction stress distribution revealed by experiments in
migrating cells [58]. Comparing the spatiotemporal distribution of the
friction and the curvature (Fig. 3$D$ and $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 [58] where traction
stress concentrates at the sides. Note that unlike many previous modeling
approaches [37, 40] 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 through the instability of the $k=1$ mode and to become
further amplified by the feedback loop between mechanosensitive adhesions,
membrane tension, and geometry.

Figure 4: (A) Stability diagram and motility modes in the
$(r_{\\_}{\text{on}},\zeta_{\\_}1)$ parameter space: blue shadings show
regions of instability for modes $k=1,\;2,\;5,\;10,\;15$ (from light to dark
blue), and symbols indicate the motility modes observed in nonlinear
simulations, whose dynamics are illustrated in the remaining panels. Five
simulations were performed for each data point. Points labeled as ‘mixed’
display multiple types of dynamics depending on the initial condition. Points
labeled as ‘others’ display irregular dynamics that are not easily
categorized. (B) Gliding mode, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(16,600)$,
$\Delta t=10$. (C) Gliding mode,
$(r_{\\_}{\text{on}},\alpha_{\\_}1)=(56,560)$, $\Delta t=10$. (D) Rotational
mode, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(24,540)$, $\Delta t=10$. (E) Zigzag
mode, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(26,600)$, $\Delta t=10$. (F) Zigzag
mode, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(36,560)$, $\Delta t=15$. (G)
Others, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(38,580)$, $\Delta t=20$. (H)
Others, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(40,540)$, $\Delta t=20$. (I)
Others, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(40,600)$, $\Delta t=20$. (J)
Others, $(r_{\\_}{\text{on}},\alpha_{\\_}1)=(44,560)$, $\Delta t=20$. In all
cases, the other parameters were set to
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;\sigma_{\\_}c=100,\;r_{\\_}p=50,\;D=0.05$.
The cells contours are colored by the fraction of bound linkers $n$, and the
purples curves show the center-of-mass trajectories in each case. See Movies
S1–S9 in the SI Appendix for videos of the dynamics.

### 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 $\zeta_{\\_}1$ and on-rate
$r_{\\_}{\text{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 wavenumbers $k=1,\;2,\;5,\;10,\;15$, while the dots denote
the long-time motility modes. Each dot 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 (Fig. 4$B$
and $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
(Fig. 4$E$ and $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;
(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” (Fig. 4$G$–$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 [59]; 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
quasi-periodic turns in the cell trajectory. We take a closer look at these
turns in Fig. 5A and 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. 4B) or a nearly
circular shape (for large on-rates, Fig. 4C), while for the zigzag and
rotational motions, the cells are more fan-shaped with curved fronts (Fig.
4D–F). Another key observation is that, regardless of the randomness in the
initial perturbations, the emerging shapes corresponding to a give set of
parameters are robust; we further elaborate on this point below.

Figure 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_{\\_}{\text{on}},\alpha_{\\_}1)=(24,540)$, and (B) zigzag mode with
$(r_{\\_}{\text{on}},\alpha_{\\_}1)=(26,600)$. The cell shapes are colored by
the local fraction of bound linkers $n$, and the orange arrows show the
instantaneous center-of-mass velocity. In both cases, other parameters are
fixed as
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;\sigma_{\\_}c=100,\;r_{\\_}p=50,\;D=0.05$.

### 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. 6A: starting from a cell
performing a gliding motion, we abruptly vary the value of
$r_{\\_}{\text{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_{\\_}{\text{on}}$, its shape quickly adjust
and start undergoing zigzags. The oscillation frequency and cell morphology
after the transition are similar to the corresponding zigzag motions with
$r_{\\_}{\text{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 $\Delta t$ within which we gradually alter the parameters, and
the intermediate states do not significantly affect the zigzag dynamics (Fig.
6B). To quantitatively compare the resulting zigzag modes emerging through
different routes, we plot the temporal evolution of the membrane tension in
Fig. 6B for different choices of $t_{\\_}0$ and $\Delta 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 and relative adhesion strength increases.

Similar transitions are also observed when we change parameters between other
modes (see SI Appendix, 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 either
approaching fixed points or stable limit cycles in phase space.

Figure 6: Transition between different motility modes. (A) Snapshots of a
gliding cell switching into a zigzag mode as the adhesion parameters
$(r_{\\_}\text{on},\alpha_{\\_}1)$ abruptly change from $(16,580)$ to
$(30,580)$ at $t=25$. Other parameters are fixed as
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;r_{\\_}p=50,\;D=0.05$. (B)
Temporal evolution of the membrane tension as the parameters
$(r_{\\_}\text{on},\alpha_{\\_}1)$ are varied to cause changes in motility
modes. Here, $\Delta t$ denotes the duration of the time interval over which
the parameters are gradually varied.

### 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 $\zeta_{\\_}1$, we calculate the time-average circularity
$4\pi 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=\sum_{\\_}{i=1}^{M}|\boldsymbol{x}_{\\_}c^{(i)}-\boldsymbol{x}_{\\_}c^{(0)}|/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
$\zeta_{\\_}1$. As the on-rate increases, the shapes all become more circular
regardless of motility mode as seen in Fig. 7A. 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. 7B, 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. 7C, 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. 7C, which are characterized by irregular turns and
may result from the superposition of multiple oscillatory periods.

Figure 7: Phase diagrams in the $(r_{\\_}{\text{on}},\zeta_{\\_}1)$ parameter
space showing color maps of the time-averaged (A) cell circularity $4\pi
A/L^{2}$, (B) membrane tension $\sigma_{\\_}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
$\epsilon=0.001,\;h=0.01,\;k_{\\_}\sigma=1000,\;\sigma_{\\_}c=100,\;r_{\\_}p=50,\;D=0.05$.

## Discussion

In this work, we have extended the one-dimensional model coupling actin
polymerization, adhesion, membrane tension, and shape change described in [10]
to two dimensions to elucidate the role of the coupling between tension and
adhesion in the determination of cell shape, the initiation of migration, and
the selection of motility modes. 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. 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 turns connects spatially
distributed adhesions and thus couples the retrograde flow with shape change
and adhesion kinetics. Certain fluctuation modes are amplified by this
feedback loop, resulting in directional motion and oscillatory behaviors.

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, in a decrease in membrane tension due to the membrane
elasticity. The force balance (10) dictates that the decrease in $n$ will
cause an increase 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 $\zeta_{\\_}1$. When
$\zeta_{\\_}1$ is large, the effect of $n$ dominates resulting in a positive
feedback on $r$, rendering the system unstable and driving the onset of
motility. On the other hand, when $\zeta_{\\_}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.

Figure 8: A simple mechanical feedback loop between membrane tension and
adhesions determines cell shape and direction of motility.

Our model predictions are consistent with 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 feedback between adhesion and membrane tension, for constant
actin flow, is sufficient to capture both the spontaneous initiation of
motility by symmetry breaking and the emergence of complex motility modes. 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, 40] 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 [13, 11, 5, 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 [56]. Finally, like other cell
motility models [27, 29, 33], we have assumed that membrane tension is
spatially uniform, yet spatial variations in tension has been observed in
cells in experiments [60]. An extension of our approach to incorporate some of
these additional details as well internal active stresses has the potential to
further enrich the model.

## References

  * [1] Jesse R Holt, Wei-Zheng Zeng, Elizabeth L Evans, Seung-Hyun Woo, Shang Ma, Hamid Abuwarda, Meaghan Loud, Ardem Patapoutian, and Medha M Pathak.  Spatiotemporal dynamics of piezo1 localization controls keratinocyte migration during wound healing.  eLife, 10:e65415, 2021.
  * [2] Thomas D Pollard and Gary G Borisy.  Cellular motility driven by assembly and disassembly of actin filaments.  Cell, 112(4):453–465, 2003.
  * [3] Revathi Ananthakrishnan and Allen Ehrlicher.  The forces behind cell movement.  Int. J. Biol. Sci., 3(5):303, 2007.
  * [4] Patricia T Yam, Cyrus A Wilson, Lin Ji, Benedict Hebert, Erin L Barnhart, Natalie A Dye, Paul W Wiseman, Gaudenz Danuser, and Julie A Theriot.  Actin–myosin network reorganization breaks symmetry at the cell rear to spontaneously initiate polarized cell motility.  J. Cell Biol., 178(7):1207–1221, 2007.
  * [5] Alex Mogilner, Erin L Barnhart, and Kinneret Keren.  Experiment, theory, and the keratocyte: An ode to a simple model for cell motility.  In Seminars in cell & developmental biology, volume 100, pages 143–151. Elsevier, 2020.
  * [6] Louise P Cramer.  Forming the cell rear first: breaking cell symmetry to trigger directed cell migration.  Nat. Cell Biol., 12(7):628–632, 2010.
  * [7] Wouter-Jan Rappel and Leah Edelstein-Keshet.  Mechanisms of cell polarization.  Curr. Opin. Syst., 3:43–53, 2017.
  * [8] Anne J Ridley, Martin A Schwartz, Keith Burridge, Richard A Firtel, Mark H Ginsberg, Gary Borisy, J Thomas Parsons, and Alan Rick Horwitz.  Cell migration: integrating signals from front to back.  Science, 302(5651):1704–1709, 2003.
  * [9] J Thomas Parsons, Alan Rick Horwitz, and Martin A Schwartz.  Cell adhesion: integrating cytoskeletal dynamics and cellular tension.  Nat. Rev. Mol. Cell Biol., 11(9):633–643, 2010.
  * [10] Pierre Sens.  Stick–slip model for actin-driven cell protrusions, cell polarization, and crawling.  Proc. Natl. Acad. Sci. USA, 117(40):24670–24678, 2020.
  * [11] Padmini Rangamani, Marc-Antoine Fardin, Yuguang Xiong, Azi Lipshtat, Olivier Rossier, Michael P Sheetz, and Ravi Iyengar.  Signaling network triggers and membrane physical properties control the actin cytoskeleton-driven isotropic phase of cell spreading.  Biophys. J., 100(4):845–857, 2011.
  * [12] Gaudenz Danuser, Jun Allard, and Alex Mogilner.  Mathematical modeling of eukaryotic cell migration: insights beyond experiments.  Annu. Rev. Cell Dev. Biol., 29:501–528, 2013.
  * [13] Yuguang Xiong, Padmini Rangamani, Marc-Antoine Fardin, Azi Lipshtat, Benjamin Dubin-Thaler, Olivier Rossier, Michael P Sheetz, and Ravi Iyengar.  Mechanisms controlling cell size and shape during isotropic cell spreading.  Biophys. J., 98(10):2136–2146, May 2010.
  * [14] Chun-Min Lo, Hong-Bei Wang, Micah Dembo, and Yu-li Wang.  Cell movement is guided by the rigidity of the substrate.  Biophys. J., 79(1):144–152, 2000.
  * [15] B J Dubin-Thaler, G Giannone, H G Döbereiner, and M P Sheetz.  Nanometer analysis of cell spreading on matrix-coated surfaces reveals two distinct cell states and STEPs.  Biophys. J., 86(3):1794–1806, March 2004.
  * [16] Benjamin J Dubin-Thaler, Jake M Hofman, Yunfei Cai, Harry Xenias, Ingrid Spielman, Anna V Shneidman, Lawrence A David, Hans-Gunther Dobereiner, Chris H Wiggins, and Michael P Sheetz.  Quantification of cell edge velocities and traction forces reveals distinct motility modules during cell spreading.  PLoS One, 3(11):e3735, January 2008.
  * [17] Kinneret Keren, Zachary Pincus, Greg M Allen, Erin L Barnhart, Gerard Marriott, Alex Mogilner, and Julie A Theriot.  Mechanism of shape determination in motile cells.  Nature, 453(7194):475–480, 2008.
  * [18] Cyrus A Wilson, Mark A Tsuchida, Greg M Allen, Erin L Barnhart, Kathryn T Applegate, Patricia T Yam, Lin Ji, Kinneret Keren, Gaudenz Danuser, and Julie A Theriot.  Myosin ii contributes to cell-scale actin network treadmilling through network disassembly.  Nature, 465(7296):373–377, 2010.
  * [19] Arnon D Lieber, Shlomit Yehudai-Resheff, Erin L Barnhart, Julie A Theriot, and Kinneret Keren.  Membrane tension in rapidly moving cells is determined by cytoskeletal forces.  Curr. Biol., 23(15):1409–1417, 2013.
  * [20] Alba Diz-Muñoz, Daniel A Fletcher, and Orion D Weiner.  Use the force: membrane tension as an organizer of cell shape and motility.  Trends Cell Biol., 23(2):47–53, 2013.
  * [21] Pierre Sens and Julie Plastino.  Membrane tension and cytoskeleton organization in cell motility.  J. Condens. Matter Phys., 27(27):273103, 2015.
  * [22] Boris Rubinstein, Ken Jacobson, and Alex Mogilner.  Multiscale two-dimensional modeling of a motile simple-shaped cell.  Multiscale Model Simul., 3(2):413–439, 2005.
  * [23] Boris Rubinstein, Maxime F Fournier, Ken Jacobson, Alexander B Verkhovsky, and Alex Mogilner.  Actin-myosin viscoelastic flow in the keratocyte lamellipod.  Biophys. J., 97(7):1853–1863, 2009.
  * [24] Thomas E Schaus, Edwin W Taylor, and Gary G Borisy.  Self-organization of actin filament orientation in the dendritic-nucleation/array-treadmilling model.  Proc. Natl. Acad. Sci. USA, 104(17):7086–7091, 2007.
  * [25] Marc Herant and Micah Dembo.  Form and function in cell motility: from fibroblasts to keratocytes.  Biophys. J., 98(8):1408–1417, 2010.
  * [26] Charles W Wolgemuth, Jelena Stajic, and Alex Mogilner.  Redundant mechanisms for stable cell locomotion revealed by minimal models.  Biophys. J., 101(3):545–553, 2011.
  * [27] Danying Shao, Wouter-Jan Rappel, and Herbert Levine.  Computational model for cell morphodynamics.  Phys. Rev. Lett., 105(10):108104, 2010.
  * [28] Danying Shao, Herbert Levine, and Wouter-Jan Rappel.  Coupling actin flow, adhesion, and morphology in a computational cell motility model.  Proc. Natl. Acad. Sci. USA, 109(18):6851–6856, 2012.
  * [29] Falko Ziebert, Sumanth Swaminathan, and Igor S Aranson.  Model for self-polarization and motility of keratocyte fragments.  J. R. Soc. Interface., 9(70):1084–1092, 2012.
  * [30] Elsen Tjhung, Davide Marenduzzo, and Michael E Cates.  Spontaneous symmetry breaking in active droplets provides a generic route to motility.  Proc. Natl. Acad. Sci. USA, 109(31):12381–12386, 2012.
  * [31] Elsen Tjhung, Adriano Tiribocchi, Davide Marenduzzo, and Michael E Cates.  A minimal physical model captures the shapes of crawling cells.  Nat. Commun., 6(1):5420, 2015.
  * [32] Wieland Marth, Simon Praetorius, and Axel Voigt.  A mechanism for cell motility by active polar gels.  J. R. Soc. Interface., 12(107):20150161, 2015.
  * [33] Erin Barnhart, Kun-Chun Lee, Greg M Allen, Julie A Theriot, and Alex Mogilner.  Balance between cell- substrate adhesion and myosin contraction determines the frequency of motility initiation in fish keratocytes.  Proc. Natl. Acad. Sci. USA, 112(16):5045–5050, 2015.
  * [34] William R Holmes and Leah Edelstein-Keshet.  A comparison of computational models for eukaryotic cell shape and motility.  PLoS Comput. Biol., 8(12):e1002793, 2012.
  * [35] Tim Oliver, Micah Dembo, and Ken Jacobson.  Separation of propulsive and adhesive traction stresses in locomoting keratocytes.  J. Cell Biol., 145(3):589–604, 1999.
  * [36] Rachael A Ream, Julie A Theriot, and George N Somero.  Influences of thermal acclimation and acute temperature change on the motility of epithelial wound-healing cells (keratocytes) of tropical, temperate and antarctic fish.  J. Exp. Biol., 206(24):4539–4551, 2003.
  * [37] Erin L Barnhart, Greg M Allen, Frank Jülicher, and Julie A Theriot.  Bipedal locomotion in crawling cells.  Biophys. J., 98(6):933–942, 2010.
  * [38] Alex J Loosley and Jay X Tang.  Stick-slip motion and elastic coupling in crawling cells.  Phys. Rev. E, 86(3):031908, 2012.
  * [39] Greg M Allen, Kun Chun Lee, Erin L Barnhart, Mark A Tsuchida, Cyrus A Wilson, Edgar Gutierrez, Alexander Groisman, Julie A Theriot, and Alex Mogilner.  Cell mechanics at the rear act to steer the direction of cell migration.  Cell Syst., 11(3):286–299, 2020.
  * [40] Kun Chun Lee, Greg M Allen, Erin L Barnhart, Mark A Tsuchida, Cyrus A Wilson, Edgar Gutierrez, Alexander Groisman, Julie A Theriot, and Alex Mogilner.  Modeling cell turning by mechanics at the cell rear.  bioRxiv, 2020.
  * [41] Falko Ziebert and Igor S Aranson.  Effects of adhesion dynamics and substrate compliance on the shape and motility of crawling cells.  PloS One, 8(5):e64511, 2013.
  * [42] Jakob Löber, Falko Ziebert, and Igor S Aranson.  Modeling crawling cell movement on soft engineered substrates.  Soft Matter, 10(9):1365–1373, 2014.
  * [43] Ulrich S Schwarz and Samuel A Safran.  Physics of adherent cells.  Rev. Mod. Phys., 85(3):1327, 2013.
  * [44] George I Bell.  Models for the specific adhesion of cells to cells: a theoretical framework for adhesion mediated by reversible bonds between cell surface molecules.  Science, 200(4342):618–627, 1978.
  * [45] Evan A Evans and David A Calderwood.  Forces and bond dynamics in cell adhesion.  Science, 316(5828):1148–1153, 2007.
  * [46] AE Filippov, J Klafter, and M Urbakh.  Friction through dynamical formation and rupture of molecular bonds.  Phys. Rev. Lett., 92(13):135503, 2004.
  * [47] Benedikt Sabass and Ulrich S Schwarz.  Modeling cytoskeletal flow over adhesion sites: competition between stochastic bond dynamics and intracellular relaxation.  J. Condens. Matter Phys., 22(19):194112, 2010.
  * [48] Ying Li, Prabhakar Bhimalapuram, and Aaron R Dinner.  Model for how retrograde actin flow regulates adhesion traction stresses.  J. Condens. Matter Phys., 22(19):194113, 2010.
  * [49] Erin M Craig, Jonathan Stricker, Margaret Gardel, and Alex Mogilner.  Model for adhesion clutch explains biphasic relationship between actin flow and traction at the cell leading edge.  Phys. Biol., 12(3):035002, 2015.
  * [50] Juliet Lee, Akira Ishihara, Julie A Theriot, and Ken Jacobson.  Principles of locomotion for simple-shaped cells.  Nature, 362(6416):167–171, 1993.
  * [51] Masoud Nickaeen, Igor L Novak, Stephanie Pulford, Aaron Rumack, Jamie Brandon, Boris M Slepchenko, and Alex Mogilner.  A free-boundary model of a motile cell explains turning behavior.  PLoS Comput. Biol., 13(11):e1005862, 2017.
  * [52] Clarence E Chan and David J Odde.  Traction dynamics of filopodia on compliant substrates.  Science, 322(5908):1687–1691, 2008.
  * [53] AC Callan-Jones, J-F Joanny, and J Prost.  Viscous-fingering-like instability of cell fragments.  Phys. Rev. Lett., 100(25):258106, 2008.
  * [54] Carles Blanch-Mercader and J Casademunt.  Spontaneous motility of actin lamellar fragments.  Phys. Rev. Lett., 110(7):078102, 2013.
  * [55] Sean P Palecek, Joseph C Loftus, Mark H Ginsberg, Douglas A Lauffenburger, and Alan F Horwitz.  Integrin-ligand binding properties govern cell migration speed through cell-substratum adhesiveness.  Nature, 385(6616):537–540, 1997.
  * [56] P Recho, T Putelat, and L Truskinovsky.  Contraction-driven cell motility.  Phys. Rev. Lett., 111(10):108102, 2013.
  * [57] C Alex Safsten, Volodmyr Rybalko, and Leonid Berlyand.  Asymptotic stability of contraction-driven cell motion.  Phys. Rev. E, 105(2):024403, 2022.
  * [58] Maxime F Fournier, Roger Sauser, Davide Ambrosi, Jean-Jacques Meister, and Alexander B Verkhovsky.  Force transmission in migrating cells.  Journal of Cell Biology, 188(2):287–297, 2010.
  * [59] Alexander B Verkhovsky, Tatyana M Svitkina, and Gary G Borisy.  Self-polarization and directional motility of cytoplasm.  Curr. Biol., 9(1):11–S1, 1999.
  * [60] Zheng Shi, Zachary T Graber, Tobias Baumgart, Howard A Stone, and Adam E Cohen.  Cell membranes resist flow.  Cell, 175(7):1769–1779.e13, December 2018.