## Abstract

To understand shapes and movements of cells undergoing lamellipodial motility, we systematically explore minimal free-boundary models of actin-myosin contractility consisting of the force-balance and myosin transport equations. The models account for isotropic contraction proportional to myosin density, viscous stresses in the actin network, and constant-strength viscous-like adhesion. The contraction generates a spatially graded centripetal actin flow, which in turn reinforces the contraction via myosin redistribution and causes retraction of the lamellipodial boundary. Actin protrusion at the boundary counters the retraction, and the balance of the protrusion and retraction shapes the lamellipodium. The model analysis shows that initiation of motility critically depends on three dimensionless parameter combinations, which represent myosin-dependent contractility, a characteristic viscosity-adhesion length, and a rate of actin protrusion. When the contractility is sufficiently strong, cells break symmetry and move steadily along either straight or circular trajectories, and the motile behavior is sensitive to conditions at the cell boundary. Scanning of a model parameter space shows that the contractile mechanism of motility supports robust cell turning in conditions where short viscosity-adhesion lengths and fast protrusion cause an accumulation of myosin in a small region at the cell rear, destabilizing the axial symmetry of a moving cell.

## Author summary

To understand shapes and movements of simple motile cells, we systematically explore minimal models describing a cell as a two-dimensional actin-myosin gel with a free boundary. The models accounts for actin-myosin contraction balanced by viscous stresses in the actin gel and uniform adhesion. The myosin contraction causes the lamellipodial boundary to retract. Actin protrusion at the boundary counters the retraction, and the balance of protrusion and retraction shapes the cell. The models reproduce a variety of motile shapes observed experimentally. The analysis shows that the mechanical state of a cell depends on a small number of parameters. We find that when the contractility is sufficiently strong, cells break symmetry and move steadily along either straight or circular trajectory. Scanning model parameters shows that the contractile mechanism of motility supports robust cell turning behavior in conditions where deformable actin gel and fast protrusion destabilize the axial symmetry of a moving cell.

**Citation: **Nickaeen M, Novak IL, Pulford S, Rumack A, Brandon J, Slepchenko BM, et al. (2017) A free-boundary model of a motile cell explains turning behavior. PLoS Comput Biol13(11):

e1005862.

https://doi.org/10.1371/journal.pcbi.1005862

**Editor: **Martin Meier-Schellersheim, National Institutes of Health, UNITED STATES

**Received: **April 7, 2017; **Accepted: **October 31, 2017; **Published: ** November 14, 2017

**Copyright: ** © 2017 Nickaeen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability: **Numerical codes used to solve the model equations are in the file RotatingCell.zip at http://cims.nyu.edu/~mogilner/codes.html.

**Funding: **This work was supported by the National Institutes of Health grants: 2 P41 RR013186-15 from the National Center for Research Resources and 9 P41 GM103313-15 and U54 GM64346 from the National Institute of General Medical Sciences. AR and JB were supported by the National Science Foundation grant DMS 1460967. AM was supported by the National Institutes of Health grant GM068952 and by US Army Research Office grant W911NF-17-1-0417. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests: ** The authors have declared that no competing interests exist.

## Introduction

Cell motility is a fundamental biological phenomenon that underlies many physiological processes in health and disease, including wound healing, embryogenesis, immune response, and metastatic spread of cancer cells [1], to name a few. Understanding the full complexity of cell motility, exacerbated by complex biochemical regulation, poses enormous challenges. One of them is multiple, sometimes redundant, sometimes complementary or even competing, mechanisms of motility [2]. Many researchers hold the view, which we share, that the way to face this challenge is to study all these mechanisms thoroughly, and then proceed with a more holistic approach.

One of the best studied types of motility is the lamellipodial motility on flat, hard and adhesive surfaces [3], in which broad and flat motile appendages–lamellipodia–spread around the cell body. Biochemical regulation plays an important role in the lamellipodial dynamics, but minimal mechanisms of the lamellipodial motility, such as growth and spreading of a flat actin network wrapped in plasma membrane and myosin-powered contraction of this network, are mechanical in nature [3]. While many cell types exhibit the lamellipodial motility, one model system, the fish epithelial keratocyte cell, contributed very prominently to the understanding of lamellipodial mechanics, due to its large lamellipodium, streamlined for rapid and steady locomotion [4, 5].

There are at least three distinct mechanical states of this system. The cells can be in a stationary symmetric state, with a ring-like lamellipodium around the cell body [6]. Spontaneously, even if slowly, the cells self-polarize, so that the lamellipodium retracts at the prospective rear and takes on a fan-like shape, upon which the cell starts crawling with a constant speed and steady shape [6, 7]. Often, cell’s trajectory changes from straight to circular–the cells start turning [8].

Mechanics of keratocyte movements has been studied extensively [4, 5, 7, 9]. Two principal mechanisms enable the keratocyte motility. First, polymerization of the polarized actin network at the front pushes forward the membrane at the leading edge, stretching the membrane and creating membrane tension at the sides; the membrane then snaps at the rear and pulls forward the depolymerizing actin network [10]. Second, contractile forces generated by myosin, lagging behind in a moving cell, hold the cell sides and retract the rear, allowing the front to protrude [5]. This and stick-slip dynamics of adhesions were recently shown to generate the cell self-polarization [7].

One of the fundamental questions of cell motility concerns dynamics of the cell shape: how do the actin-myosin mechanics in the cell bulk interact with actin growth and membrane mechanics at the boundary to shape, stabilize and propel the cell? This question requires mathematical insight, and in the last two decades, keratocyte mechanics were extensively modeled mathematically. The mathematical problem arising in these models is generally challenging, given that the motile cell is a free-boundary object, in which deformations of the cell shape depend on, and in turn affect, the actin-myosin movements and forces inside the cell. The history of the free-boundary cell modeling was recently reviewed in [11].

To reduce the mathematical complexity of the problem, one can ignore the mechanics in the lateral cross section and consider a simplified one-dimensional (1D) model, essentially representing the cell as a 1D strip of an actin-myosin gel. Mathematical models of this kind [12–15] provided valuable insight into conditions for symmetry breaking, motility initiation, and stabilization of the anterior-posterior length of the moving cell. Modeling the front-to-rear cell mechanics is not the only 1D approach: one can also disregard the bulk of the actin-myosin network and hypothesize that the essential dynamics is concentrated at the very edge of the cell; this allows one to approximate the cell shape by a 1D dynamic contour, which protrudes or retracts locally according to some set of rules. A number of such models [16–19] revealed that a small set of the boundary deformation rules can generate an unexpected variety of dynamic cell shapes mimicking a number of observed motile cell types. The first such model was a Graded Radial Extension mathematical model [20], which integrated experimental data and posited that actin polymerization at the lamellipodial boundary results in protrusion of the cell front and sides in the direction locally normal to the boundary, with spatially graded rate maximal at the center of the leading edge and decreasing towards the sides.

A more accurate mathematical rendering of the lamellipodia is achieved via a full 2D free-boundary model. Its general concept, first introduced in [21], is as follows. Actin-myosin contraction in the bulk of the 2D lamellipodium generates a centripetal actin flow that redistributes myosin powering the contraction; this feedback results in a spatially graded flow that tends to retract the lamellipodial boundary. Actin growth at the boundary results in protrusion countering the retraction, and the balance of protrusion and retraction shapes the lamellipodium, feeding back to the actin-myosin contraction in the bulk of the lamellipodium. The question is: what kind of cell shapes and movements does this model predict?

To address this question, 2D models of actin-myosin mechanics were employed to reproduce steady-state shapes and speed, as well as self-polarization, of a motile cell [5, 7, 22], but a number of important issues have not been adequately explored, including turning behavior and dependence of the motile behavior on the model parameters and boundary conditions. In this paper, we resolve these issues by simulating numerically a minimal free-boundary model described in the next section. We find that 1) cells are stationary when contractility is weak, 2) when contractility is strong, cells break symmetry and move steadily along either straight or circular trajectory, 3) cells exhibit turning behavior when fast protrusion destabilizes the axial symmetry of a planar myosin distribution or cell shape, and 4) motile behavior of a cell is sensitive to conditions of force balance at the cell boundary.

## Models and methods

Mechanics plays a dominant role in keratocyte motility, while the role of biochemical regulation is less clear, and is probably of less importance [3]; thus we focus on mechanical formulation. Moreover, at least for the self-polarization phenomenon, the contractile mechanism of motility is dominant compared to the graded actin polymerization [7], and so we concentrate on the myosin generated forces and movements, and for simplicity assume that the actin growth is uniform around the cell boundary.

The model consists of force-balance and myosin transport equations,

(1)

for the velocity of actin flow **U**(**X**,*t*) and myosin concentration *M*(**X**,*t*) defined locally for **X** ∈ *Ω*(*t*), where *Ω*(*t*) is a moving 2D domain representing cell geometry. Note, that all model equations are formulated in the lab frame of references. We discuss the validity of approximating the flat lamellipodium as a 2D domain in the thin-cell limit in the Supplemental Material. The model is similar to an actomyosin dynamic model suggested in [5, 23], and to active gel models of the soft matter physics [24]. In the force balance equation, the first term describes the force due to passive viscous stresses in the deforming actin network, where *η* is the effective actin viscosity. The form of this term corresponds to the viscous shear stress in the Stokes equation of hydrodynamics; we emphasize that the actin polymer mesh is compressible (fluid cytoplasm can be squeezed easily into the dorsal direction in the cell [23]), and so there is no incompressibility condition. In the Supplemental Material, we discuss the conditions under which the effects of hydrostatic pressure and Darcy flow are negligible in the lamellipodium. Because the movement of the cell takes place on the slow time scales, we do not consider viscoelastic effects [23]. Also, as was done in many modeling studies [25], we ignore for simplicity subtle and complex interplay between bulk and shear viscous stresses.

The second term in the force-balance equation describes the divergence of the myosin contractile stress. As in [5, 23, 24], we assume the stress to be isotropic and proportional to the myosin density, with *σ* denoting the force per unit of myosin density. The third term describes the effective viscous drag arising from creeping movement of F-actin relative to a substrate, mediated by dynamic adhesions and characterized by the viscous drag coefficient *ξ*. The linear dependence of this drag on actin velocity is a standard assumption made in cell mechanics models; for some cases, this assumption was verified experimentally.

The second of Eq (1) describes myosin transport. Kinetics of myosin can be interpreted in terms of transitions between two states, a state of free myosin diffusing in the cytoplasm and a state in which myosin is bound to the actin network [23]; clusters of the bound myosin both contract the actin network and move with it. For the transitions occurring on a fast time scale, the overall transport is well approximated by a diffusion-advection equation, and we assume this limit in our model here. Additional discussion of the myosin transport is in the Supplemental Material. For low viscosity and slow diffusion, however, using **U** as an advection velocity and constant diffusion coefficients *D* results in singular solutions, in which *M* and **U** develop Dirac-delta singularities. The effect is reminiscent of the collapsing phenomenon in a 2D version of the Keller-Segel chemotaxis model [26], which is mathematically similar to our mechanical model. The singular solutions are obviously unrealistic, given that myosin molecules have a finite size. The excluded volume effect can be taken into account by introducing effective velocities [27], **U**_{eff} = **U**(1−*M*/*M*_{max,u}), which approach zero when *M* → *M*_{max,u} But in a free-boundary problem, the actual maximum of myosin concentrations may significantly exceed *M*_{max,u}. This is because in motile solutions, myosin accumulates at the rear of the cell, where it is also swept forward by a moving boundary; mathematically, this effect originates from the Rankine-Hugoniot boundary condition described in the next paragraph. Thus, the effect of molecular crowding on myosin velocity should generally be written as: **U**_{eff} = **U**(1−*M*/*M*_{max,u}), if *M* < *M*_{max,u}, and zero otherwise. Because diffusion is also affected by the crowding, the effective diffusivity in Eq (1) is *D*_{eff} = *D*(1−*M*/*M*_{max,d}), with a value of the cut-off *M*_{max,d} that exceeds the actual maximal densities of myosin [28]. Overall, using a tighter myosin cut-off for advection, *M*_{max,u} < *M*_{max,d}, helps avoiding the singularities and numerical instabilities associated with them in a wide range of model parameters. We also explored the possibility that the anti-crowding effects come from an attenuation of the myosin stress when the myosin density is too high, described in detail in the Supplemental Material.

One of our goals is to investigate how boundary conditions, specified in a free boundary problem at a moving cell boundary *∂Ω*(*t*), affect the model behavior. We explore two types of conditions for the force-balance equation. One of them is the zero actin velocity at the boundary, **U**|_{∂Ω(t)} = 0, which assumes a very narrow band of very strong adhesions near the cell edge that rapidly adjust their positions to the instantaneous location of the boundary. Many experimental observations indicate the presence of such a band. We will term this version of the model a zero-velocity (ZV) model. Alternatively, in the absence of the sticky band, the force balance is reflected by a zero stress condition, introduced earlier in [5, 23]: , where is the unit tensor and **n** is the outward normal. This assumes that the membrane tension is small relative to the contractile and viscous stresses. We will call this variant of the model a zero-stress (ZS) model. Both versions of the model share a no-flux Rankine-Hugoniot boundary condition for myosin, **n** ⋅(−*D*_{eff} ∇_{X}*M* + (**U**_{eff} – **V**_{f}(*M*) |_{∂Ω(t)} = 0, where **V**_{f} is the local boundary velocity.

Kinematics of the boundary is modeled by a superposition of locally normal protrusion powered by actin growth and retraction stemming from the centripetal actin flow: **V**_{f} = *V*_{p}**n** + **U**|_{∂Ω(t)}. The approximation of the speed of normal protrusion *V*_{p} is somewhat different in the two versions of the model, as described below.

In the ZS model, *V*_{p} is defined uniformly along the boundary but depends on the cell size: *V*_{p} = *V*_{0}(*A*_{0}/*A*)−*K*(*A*−*A*_{0}(*A*_{0}/*A*)^{n}), where *A* = |*Ω*(*t*)|, *A*_{0} = |*Ω*(0)|, and *n* = 2. The first term represents the rate of membrane displacement due to actin growth with a rate constant *V*_{0}. The cell-size dependence of this term reflects an effective drop in actin concentration in an expanding cell, but this term alone would still produce an infinite cell expansion for large *V*_{0}. Realistically, cell stretching is limited by membrane tension, which is represented in *V*_{p} by −*KA* term; this is consistent with previous experiments and modeling [4, 5, 29]. The term ∝(*A*_{0}/*A*)^{n} reflects cytoplasmic resistance to contractile forces and thus excludes collapsing of the cell in the model with small *V*_{0}. Mathematically, the quadratic nonlinearity in the density dependence appears to be the lowest nonlinearity preventing the cell collapse in the ZS model. Overall, the second term in *V*_{p}, combining the effects of membrane tension and cytoplasm resistance to contraction, plays an area-preserving role (with the parameter *K* describing sensitivity of *V*_{p} to changes in *A*). Indeed, if *A*<*A*_{0}, the membrane tension decreases, whereas actin polymerization accelerates and the resistance to further contraction rapidly grows. On the other hand, if *A*>*A*_{0}, the membrane tension increases rapidly stopping the actin growth.

For the ZV model, where **U**|_{∂Ω(t)} = 0 and **V**_{f} = *V*_{p}**n**, there must be a nontrivial variation of *V*_{p} along the boundary, since for a uniform *V*_{p}, the cell centroid is always stationary. Based on experimental observations and models showing that myosin can impede protrusion by bundling actin filaments at the boundary [18], we hypothesized that the actin growth rate is a decreasing function of local myosin density. Correspondingly, we use the following expression for *V*_{p} in the ZV model, *V*_{p} = *V*_{0}(*A*_{0}/*A*)(1+*M*)|_{∂Ω(t)}/M_{0})^{−1}−*K*(*A*−*A*_{0}), where *M*_{0} is a threshold beyond which myosin inhibits actin growth almost entirely. The expression has essentially the same dependence on cell area as in the ZS model, except that for the ZV model, *n* = 0 proved to be sufficient for preserving the target area. We discuss derivation of the mathematical expression for *V*_{p} from the force balance at the lamellipodial boundary and provide additional explanations in the Supplemental Material.

### Nondimensionalization

To nondimensionalize the model, we use the following set of units. The length unit *L* is defined as a characteristic linear size of the cell with a target area, (e.g., if this cell is a circle, *L* is its radius). We further choose *L*^{2}*D*^{−1} and *M*_{0} as the units of time and myosin concentration, respectively. Then the dimensionless variables, differential operator, and current and target cell areas are, respectively, *t* = *TDL*^{−2}, **x** = **X***L*^{−1}, **u** = **U***LD*^{−1}, *m* = *M*/*M*_{0}, ∇ ≡ ∇_{x} = *L*∇_{X}, |*ω*|=|*Ω*|*L*^{−2} and *a*_{0} = *A*_{0}*L*^{−2}. Correspondingly, Eq (1) takes the form,

(2)

where **u**_{eff} = **u**(1−*m*/*m*_{max,u}), if *m*<*m*_{max,u}, and zero otherwise, and *d*_{eff} = 1−*m*/*m*_{max,d} Eq (2) include two dimensionless parameters: the dimensionless viscosity-adhesion length parameter *α* = *ηL*^{−2}*ξ*^{−1} and the myosin contractility constant *β* = *σM*_{0}(*Dξ*)^{−1}. Note that the mechanical effect of localized myosin contraction spreads on the length scale , so the viscosity-adhesion length parameter *α* is the ratio of the length scale of the mechanical action to the cell size. The first of Eq (2), a diffusion-advection equation for myosin, is subject to the mass-conserving zero-flux boundary condition, **n** ⋅(−*d*_{eff}∇*m*+(**u**_{eff}−**v**_{f})*m*)|_{∂ω(t)} = 0, yielding an additional dimensionless parameter . The dimensionless boundary conditions for the force-balance equation (the second of Eq (2)) in the ZV and ZS models are **u**|_{∂ω(t)} = 0 and , respectively.

The dimensionless boundary velocity equation is **v**_{f} = *v*_{p}**n**+**n**|_{∂ω(t)}. In this equation, the dimensionless rate of membrane displacement caused by the actin polymerization and area preserving factors is *v*_{p} = *v*_{0}*a*_{0}/*a*−*k*(*a*−*a*_{0}(*a*_{0}/*a*)^{n}), where *v*_{0} = *V*_{0}*L*/*D* and *a* = |*ω*(*t*)|. For the ZS model, *n* = 2, whereas for the ZV model, *n* = 0 and there is the additional dependence on *m* in the first term: *v*_{p} = *v*_{0}*a*_{0}/(*a*(1+*m*|_{∂ω(t)}))−*k*(*a*−*a*_{0}).

Note that varying parameter *k* = *KL*^{3}*D*^{−1} is equivalent to rescaling the actin polymerization constant *v*_{0}. Also, because the myosin contractility constant *β* enters Eq (2) in combination with *m*, varying *β* is similar to rescaling *μ*_{tot}; in fact, *β* could be formally excluded from the ZS model by employing a different concentration unit, and the same is true for the ZV model defined in a fixed geometry, see section *Cell becomes motile when myosin contractility is higher than critical*. We therefore focus in our study on the role played by three essentially independent model parameters: *α*, *μ*_{tot} and *v*_{0}.

### Initial conditions

Steady dynamics of a motile cell were explored by solving Eq (2) in domains with free boundaries. Note that even though the force-balance equation does not involve time derivatives in and of itself, the coupled system (2) constitutes an initial-value problem and one must specify initial conditions for both variables and the domain, *m*(**x**,0), **u**(**x**,0), and *ω*(0).

To elucidate processes leading to instability of an initially symmetric stationary state of a motile cell and its transitioning to motility, we used initial conditions based on a stationary steady state of the ZV model in a circular geometry *ω*(0) = {(*x*,*y*):*x*^{2}+*y*^{2}≤1} (such that *a*_{0} = |*ω*(0)| = *π*): **u**(**x**,0) = 0, *m*(**x**,0) = (*μ*_{tot}/|*ω*|)(1−*gx*). Note the linear horizontal gradient, added to a steady-state uniform myosin distribution to probe stability of a stationary state; the gradient steepness *g* was assigned values from (0,1]. Note also that given the symmetry of *ω*(0), the definition of *m*(**x**,0) ensures that the solution has a prescribed *μ*_{tot}. The initial conditions specified above were used in solving both ZV and ZS models throughout this study.

### Numerical methods

Numerical solutions of the ZV and ZS models were obtained using a generalized version of a mass-conservative algorithm originally developed for solving parabolic equations in moving domains with known kinematics [30]. The method has been shown to converge in space with an order close to 2 in L^{2}-norm and ensures exact local mass conservation. The latter is achieved by employing finite-volume spatial discretization [31] and natural neighbor interpolation [32]. The algorithm was developed for modeling cell motility in *Virtual Cell* (VCell), a general-purpose computational framework for modeling cellular phenomena in realistic geometries [33].

To be applicable to a free-boundary problem with the models described above, the original method had to be augmented in several aspects. First, the boundary kinematics is generally not known *a priori* but rather needs to be computed on the basis of the rates that are functions of state variables—the actin velocities, in the ZS model, and the myosin concentrations, in the ZV model. To approximate values of the variables at the points of the boundary where the boundary velocities need to be evaluated, we used the second-order bilinear extrapolation. Once the boundary velocities are obtained, the cell boundary is advanced using a robust front-tracking technique implemented in FronTier, a freely available C++ library for tracking interfaces in two and three dimensions [34]. Accuracy of the algorithm coupled to FronTier was evaluated using several benchmark examples, one of which was based on the models of this study. The tests have shown that the accuracy of the original algorithm is preserved, if in addition to the second-order extrapolation, the front-tracking routine is also at least second-order accurate.

Second, the system (2), consisting of the coupled parabolic and elliptic equations, is nonlinear. Indeed, the equations are coupled via the advection term of the parabolic equation and myosin-dependent stress term in the elliptic equation, as well as through the boundary conditions at the moving boundary; also, the effective transport parameters are functions of the myosin concentration. To solve the system, we implemented a segregated solution strategy [31], in which equations are solved one at a time and nonlinear terms are treated by fixed-point iterations. One advantage of the segregated solver is that it prevents the matrix of a linearized system from becoming very large even with very fine computational grids. The system was advanced in time using an implicit backward Euler time discretization.

For each time step, the segregated method performs fixed-point iterations in two steps. First, we solve for actin velocities using fixed myosin concentrations from the previous iteration. The obtained velocities are then used as a fixed advection field at the next step, where we solve the linearized transport equation for myosin concentrations. Note that at this step, the values of the myosin concentrations in the discretized time derivative correspond to the previous time step, not to the previous fixed-point iteration. At the end of the iteration, maximum absolute differences of two consecutive iterates are checked for convergence. If they are within prescribed tolerances, the iterations stop and the solver reports the velocities and myosin concentrations as the current time step values, otherwise it proceeds to the next iteration and continues until the iterations converge or an imposed maximum for the number of iterations is exceeded. The algorithm is illustrated below for one time step by the mathematical pseudocode, where *m*^{k} and **u**^{k} are the variable values at the *k*th time step, and are the *n*th iterates for the (*k*+1)th time step, MaxNumIters is the maximum allowed number of iterations, and ‖⋅‖_{∞} denotes the L^{∞}-norm.

set and

for *n* = 1: MaxNumIters

– solve to get

– evaluate **u**_{eff} and *d*_{eff} using and

– solve to get

– calculate absolute errors and

– if solution converged, break the loop, else ,

end of segregated loop

if *n*< MaxNumIters , , else iterations are stagnant.

The segregated solver was validated against a coupled nonlinear solver implemented in COMSOL Multiphysics [35]. Good agreement was observed, with relative differences below 0.3%.

The computations were performed with the following solution parameters: the mesh sizes *h* varied between 0.05 and 0.16, whereas the time step was *Δt* = *ch* with *c* varying from 0.0002 to 0.025 (fast-moving cells required smaller mesh sizes and time steps), the tolerance for the differences of consecutive iterates was 1E-10, and the maximum allowed number of iterations, set at 35, was never reached.

## Results

The ZV and ZS models were used to simulate transitions of a motile cell from stationary to motile state. For this, as described in section *Initial conditions*, an initially radially symmetric cell was perturbed by superimposing a linear gradient over a uniform distribution of myosin. The emerged steady states fall into three asymptotically stable mechanical modes. For some parameter values, the cell, after a finite displacement, comes to a stop, with a final radially symmetric shape and myosin distribution, indicating the stability of the stationary state (Fig 1*A*). For other parameters, the cell irreversibly breaks symmetry, both in terms of its shape, distribution of myosin, and actin velocity field, and either acquires unidirectional motility (Fig 1*B*) or locks in a rotational mode (Fig 1*C*).

Fig 1.

Asymptotically stable mechanical states (in all cases, *α* = 0.5): (*a*) a stationary state of ZS model, (*v*_{0}, m_{tot}) = (2.5, 0.125*π*). (*b*) unidirectional translation in ZV model, (*v*_{0}, m_{tot}) = (2.5, 2*π*). (*c*) rotations in ZS model, (*v*0, mtot) = (2.5, 0.75*π*). Pseudo-colors depict distributions of myosin; arrows are actin velocities; a red dashed line/curve shows the trajectory of a centroid.

To analyze conditions for transitioning to different types of motility, we scanned the actin growth constant (*v*_{0}) and the contractility parameter (*μ*_{tot}) for two values of viscosity-adhesion length parameter, *α* = 0.5 and *α* = 1. The values of other model parameters, *β* = 5, *a*_{0} = *π*, *k* = 1.5, *m*_{max,u} = 15 and *m*_{max,d} = 125, were fixed in all simulations; the choice of these values ensures that the corresponding section of parameter space is representative of various states. The results of parameter scanning are presented in Fig 2 showing cell mechanical states as functions of the model parameters.

It should be noted that distinguishing between translational and rotational modes is sometimes ambiguous, particularly for the states near the borders between the corresponding regions in the parameters space. For example, some states of the ZS model shown in the diagrams of Fig 2 as translating were in fact only ‘piecewise unidirectional’, as the cell in those states would on occasion change its direction. Moreover, in some states in the ZS model, identified as translations, also neighboring the rotations in the diagrams of Fig 2, the cell actually changes its direction but very gradually, so the state could be a rotation with a very long radius. As a practical rule, we labeled states as rotations only if the radius of rotation of the cell centroid was comparable to, or less than the linear size of the cell.

Below we discuss in detail the conditions required for the straight and rotational motility in our models and the underlying mechanisms.

### Cell becomes motile when myosin contractility is higher than critical

The results in Fig 2 show that cells break symmetry and transition to motility when parameter*μ*_{tot} exceeds a threshold. The threshold behavior originates from a positive feedback between the actin flow and myosin gradients: the contractile forces, generating the centripetal flow of myosin, are proportional to the myosin gradients, which, in turn, are reinforced by the advection of myosin. This positive feedback results in steep myosin gradients and, potentially, symmetry breaking, but below a critical value of *μ*_{tot}, these gradients are prevented by dissipative viscous forces and myosin diffusion, and the cell remains stationary and radially symmetric. Above the critical value of *μ*_{tot}, steep gradients of myosin are developed and the radially symmetric stationary state becomes unstable. While kinematics of a free boundary plays a key role in the symmetry breaking and transitioning to motility in both models, the loss of symmetry in the ZV model also occurs in fixed domains. In the next section, we discuss effects of the boundary conditions for actin flow on solutions in domains with fixed and free boundaries.

Linear stability analysis can be used to estimate critical values of *μ*_{tot} in a simplified ZV model in a fixed domain. Consider a 1D ZV model on the fixed-length segment *∂ω* = {*x*∈(0,1)} in the limit *m*_{max,1}, *m*_{max,2} → ∞, so that *d*_{eff} = 1 and *u*_{eff} = *u*(*x*,*t*). Then, Eq (2) become *αu*_{xx} + *βm*_{x}−*u* = 0, *m*_{t} = *m*_{xx}−(*um*)_{x}. In this model, varying the contractility constant *β* is equivalent to rescaling *μ*_{tot}. Indeed, *β* could be excluded altogether by renormalizing *m*: , but in what follows parameter *β* is kept for generality. The symmetric steady state is characterized by uniform myosin distribution and absence of actin flow, *u* = 0, *m* = *μ*_{tot}. Its stability is probed by imposing small perturbations, *δu*(*x*,*t*) = *δu*_{0} exp(*λt*+*iqx*), *δm*(*x*,*t*) = *δm*_{0} exp(*λt*+*iqx*) with 0<*δu*_{0}<<1, 0<*δm*_{0}<<1 and *q* = *π*, 2*π*,…, so that *u* = *δu* and *m* = *μ*_{tot}+ *δm*. The perturbations satisfy the linearized system of differential equations, *αδu*_{xx}+*βδm*_{x}−*δu* = 0, *δm*_{t} = *δm*_{xx}−*μ*_{tot}(*δu*)_{x}, and the corresponding system of linear algebraic equations, −(1+*αq*^{2})*δu*_{0}+*iqβδm*_{0} = 0, *iqμ*_{tot}*δu*_{0}+(λ+*q*^{2})*δm*_{0} = 0, yields nontrivial solutions for *δu*_{0} and *δm*_{0}, if *λ*(*q*) = *q*^{2}(*βμ*_{tot}(1+*αq*^{2})^{−1}−1). The perturbations grow if *λ*(*q*)>0, with the fastest growing mode having the minimal wave number, *q*_{min} = *π*, and thus involving a large-scale redistribution of myosin. Thus, the symmetric state becomes unstable if *βμ*_{tot}>1+*π*^{2}*α* or, in the dimensional form, *σM*_{0}*πL*^{2}>*D*(*ξL*^{2}+*π*^{2}*η*).

The instability criterion predicts that the critical value of *βμ*_{tot} is an increasing function of *α*. The results of Fig 2 indicate that this prediction, while obtained by analyzing a ZV model in a fixed domain, holds for the free-boundary models as well. Indeed, the competition between the myosin contractile stress and dissipative processes, mathematically expressed in the instability condition, drives the initiation of motility in the free-boundary models. As described at the beginning of this section, the transition to motility occurs when the contractility, reinforced by the model positive feedback, prevails over the dissipation. For *α*≥1, the dimensional form of the instability criterion reduces to *σM*_{0}*L*^{2}>*πDη*: the total myosin stress needs to overcome the smoothing effects of actin viscosity and myosin diffusion, while the adhesion strength does not matter. In the limit of large values of *α*, the actin network is effectively stiff and thus does not allow for significant actin flows, which makes the cell more symmetric and as a consequence less motile and slower. Our simulations confirm this prediction (Fig 3A and 3B). In the opposite limit of highly deformable actin network, *α*<<1, the instability criterion reduces to *σM*_{0}>*ξD*, so the cell becomes motile if the characteristic myosin stress is able to generate actin flow that overcomes myosin diffusion, which requires the weakening of adhesions and strengthening of myosin, in agreement with the experiment [7].

Fig 3. Aspect ratios and translational or linear rotational speeds of steadily moving cells.

(*a*) Aspect ratio as a function of viscosity-adhesion length a; the results were obtained with (*v*_{0}, m_{tot}) = (2.5, 1.5p) for ZS model, and with (*v*_{0}, m_{tot}) = (5, 1.5p) for ZV model; aspect ratios were computed as ratios of the longest to shortest distances between cell boundary and cell centroid. (*b*) Dimensionless translational or linear rotational speed of a cell centroid as a function of viscosity-adhesion length a; model parameters are as in panel (*a*). (*c*) Aspect ratio as a function of *v*_{0} and m_{tot}; the results were obtained with a = 1 for ZV model and with α = 0.5 for ZS model. (*d*) Radius of rotation of a cell centroid as a function of *v*_{0} and m_{tot}, with values of a as in (c). (*e*) Dimensionless angular velocity of a cell centroid as a function of *v*_{0} and m_{tot}, with values of a as in (c).

Finally, it is worth noting that whereas the motility threshold in the ZS model is independent of *v*_{0}, the critical values of *μ*_{tot} in the ZV model, where the actin growth is affected by myosin, vary with *v*_{0} (Fig 2). Indeed, the cell described by the ZV model with *v*_{0} = 0 would not transition to motility for any *μ*_{tot}, because in this case, the myosin influence on the boundary is lost. Therefore, in this version of the model motility initiates only for finite values of *v*_{0}, which drop with the increase of *μ*_{tot}. In the ZS model, the cell with a sufficiently high *μ*_{tot} initiates motility even as *v*_{0}→0, because the asymmetric myosin pulls the boundary inward asymmetrically, and the area-preserving term causes the effective protrusion.

### Shape and speed of the motile cell

Similar to the 1D ZV model analyzed in the previous section, the symmetric state of the 2D version of the ZV model becomes unstable for sufficiently high *μ*_{tot} even in a fixed geometry, with myosin relocating to the cell boundary. Fig 4A illustrates the instability of the radially symmetric steady state, in which the maximum of myosin was slightly shifted to the left of the cell center. Qualitatively, because of the zero actin velocities at the boundary, a second, initially small, local maximum of myosin appears at the boundary point closest the main maximum due to slightly faster diffusion. The competition between the two maxima lowers the myosin gradients on the left side of the main one, resulting in a net force acting on it in the left direction. Hence, the relocation of myosin to the left segment of the boundary. In the cell with a free boundary, the redistribution of myosin is conferred to boundary velocity, resulting in slower outward and eventually inward movements of the part of the boundary that becomes the cell rear. The cell movement further skews the myosin towards the rear. For the small to moderate rate of actin growth and cell speeds, the cell maintains a convex shape and a steady unidirectional motion, with myosin forming a wide band at the rear edge (Fig 1*B* and S1 Movie).

Fig 4. Symmetry breaking in a fixed circle and in a free-boundary problem.

(*a*) Instability of a symmetric steady state of ZV model in a fixed circle: snapshots of dimensionless myosin density (pseudo-color) and actin velocities (arrows) at specified times *t* after myosin was slightly shifted left of center; computations were done for m_{tot} = 1.5p and α = 0.5. (*b*) Transition to unidirectional motility in ZS model; dimensionless myosin concentration (pseudo-colors) and boundary velocities (arrows) are shown for the solution obtained with α = 1, *v*_{0} = 5, and m_{tot} = 1.5p; the cell assumes steady unidirectional motility after *t* = 14 (S2 Movie).

In the ZS model of a fixed symmetric cell, an inward actin flow at the membrane prevents myosin from accumulating there. As a result, the symmetric solution with a myosin maximum at the center remains stable even for *μ*_{tot} above the threshold. Indeed, shifting the maximum from the center in this case increases the myosin gradients and centripetal forces on the ‘shorter’ side and decreases them on the ‘longer’ side, netting a stabilizing force. In the free-boundary problem, however, the symmetric solution, stable at low *μ*_{tot} (Fig 1*A*), becomes unstable above the motility threshold. Fig 4B and S2 Movie illustrate a transition to unidirectional motility that occurs in the ZS model with sufficiently large *μ*_{tot} and small to moderate *v*_{0}. As the myosin cluster shifts slightly from the cell center, the closer side is pulled inward faster and becomes the prospective cell rear, while the opposite side, where the protrusions are faster than retractions, becomes the cell front. In the ensued motility, myosin is pressed against the rear and in turn exerts a stronger inward force on the proximal portion of the boundary, developing a local concavity. If the movement is sufficiently fast, the myosin spreads along the portion of the boundary with negative curvature. This positive feedback between the myosin asymmetry, actin flow, and cell movement is the key to the stable motility.

Our simulations show that the aspect ratio of a steadily moving cell varies between 1 and 3 (Fig 3C), in agreement with experimental observations [5]. Note that in contrast to the ZV model, where the cell aspect ratio grows moderately with *v*_{0}, it becomes essentially independent of model parameters in the ZS model with *μ*_{tot}/*π*>1. This can be qualitatively understood by noting that myosin, which in a moving cell accumulates in the middle of the rear, exerts comparable forces on the front and side portions of the boundary. Then, because the myosin-generated flow decreases with distance at similar rates in all directions, the distances from the rear to the front and the sides should be on the same order, yielding the average aspect ratio ~ 2.

The ZS model generally predicts significantly higher cell speeds compared to those in the ZV model (Fig 5A). This is because in the ZS model, the fast centripetal flows generated by myosin at the rear boundary tend to decrease the cell area, leading to fast effective protrusion at the front, as actin can grow rapidly against the lowered membrane tension. As a result, the cell speed increases but the cell area decreases with total myosin (Fig 5B). Interestingly, the cell speed in the ZS model decreases slightly with the actin growth rate *v*_{0}, which can be understood by noting that the cell area in this model increases with *v*_{0}, thus mitigating the effect of myosin. In the ZV model, the cell speed is virtually insensitive to *μ*_{tot}, for *μ*_{tot}>*π*, because the term with *v*_{0} in the expression for *v*_{p} becomes inessential for *m*∼*μ*_{tot}/*π*>1. For the same reason, the cell area is also insensitive to *μ*_{tot}.

Fig 5. Cell speeds and areas as functions of *v*_{0} and m_{tot}.

(*a*) Dimensionless translational or linear rotational speed of a cell centroid were obtained with α = 1 for ZV model and α = 0.5 for ZS model. (*b*) Dimensionless steady-state areas for values of α as in panel (*a*). Insets in both panels: corresponding level sets.

Importantly, our model predicts that there are no short-wavelength instabilities in the cell shape (like fingering instabilities characteristic for some physical free-boundary models), which is supported by the experiment: there are small fluctuations on the experimentally observed cell boundaries, but they mostly do not grow.

### Mechanics of the straight and turning motility

The most nontrivial and important property of our models is that they predict rotational states with a radius of rotation comparable to the cell size in large regions of their parameter space (Fig 2). Note that both the radius of rotation and the angular velocity are not particularly sensitive to parameter values (Fig 3D and 3E).

Emergence of cell turning in the models can be qualitatively understood by analyzing the loss of stability of a planar axial symmetry characteristic of the straight moving cell. In the ZV model, the steady rotations are observed for *v*_{0} exceeding a threshold that is largely insensitive to either *α* or *μ*_{tot}Fig 6 and S3 Movie illustrate, for a particular parameter set, how rotations come about in the ZV model during a transient movement following a ‘nudge’ applied to a stationary cell in the form of an initial horizontal gradient of myosin. The initial convex cell shape is favorable for maintaining a unipolar axially symmetric myosin distribution, and the resulting motility is unidirectional. But due to sufficiently large *v*_{0}, fast boundary velocities at cell’s sides tend to elongate the cell in *y*-direction, making it prone to developing a concavity at the cell rear. In such a shape, the myosin spreads along the rear part of the boundary more uniformly (Fig 6*B*, *t* = 5), a distribution that is no longer stable. Indeed, even a slight asymmetry in the distribution of myosin, reinforced by the positive feedback from actin velocities and myosin accumulation due to faster movement of the corresponding portion of the rear boundary, breaks the axial symmetry (Fig 6*B*, *t* = 15). As a result, a stable asymmetric cell shape emerges as the cell locks in rotations (Fig 6*A*), with myosin aggregated at a high-curvature portion of the boundary.

Fig 6. Onset of steady rotations in ZV model, (*v*_{0}, m_{tot}, α) = (12.5, 2p, 1).

(*a*) Entire cell trajectory and cell centroid track (red dashed curve). (*b*) Snapshots of transient myosin distributions with individual color scales during a transient, and with white arrows representing actin velocities (S3 Movie).

While the same mechanisms underlie the turning behavior in the ZS model, the two models yield significantly different results for the parameter regions of rotations. This is due to the differences in boundary conditions that reflect the opposing assumptions about the strength of adhesions at the cell periphery, and in ways of conferring the myosin dynamics onto kinematics of the boundary. Unlike the ZV model, the *v*_{0} threshold for rotations in the ZS model strongly depends on *μ*_{tot} and *α* (Fig 2). In particular, the rotational states may exist for any *v*_{0}, if *α* is sufficiently small. Note also, that the concavity of the cell shape does not always destabilize unidirectional motility in the ZS model (Fig 4A).

Fig 7 and S4 Movie illustrate the onset of turning in ZS model. If the contractility due to myosin is strong, myosin forms a radially symmetric aggregate, which in a translating cell is skewed to the cell rear, pulling the rear boundary inward and maintaining the cell propulsion. When the myosin aggregate is sufficiently close to the rear boundary, it pulls the center of the cell rear inward stronger than the sides of the rear edge, creating a ‘dip’ at the center of the rear edge and giving the cell a characteristic keratocyte fan-like shape, in which the sides of the cell lag behind the center. For the parameters in the upper left corner of the parameter space (Fig 2), the cell motion is fast, and in the frame of the cell, myosin is effectively swept towards the rear and ‘pressed’ against the rear boundary; in these conditions, the translational motility remains stable (Fig 4A). For intermediate values of *μ*_{tot}, the cell moves slower and the myosin aggregate maintains its radial symmetry and remains close to the cell centroid (Fig 7*A*, *t* = 7). In this position, myosin is able to pull inward not only the rear but also the front of the cell, making the axial symmetry of the system unstable. Indeed, even a slight random asymmetry in either the myosin distribution or the cell shape induces and reinforces the asymmetry of the other. If, for example, the myosin aggregate becomes slightly closer to one side, this side is pulled inward faster than the other, which brings even more myosin to the side that is pulled inward, because the shift of that side effectively sweeps myosin towards it (Fig 7*A*, *t* = 9).

(*a*) Transient distributions of myosin (pseudo-colors) and actin velocities (arrows): *t* = 2, an initially symmetric cell with centroid at (*x*, *y*) = (0,0) self-polarizes and assumes fast unidirectional motility, myosin accumulates in a semi-circular band, pulling the rear inwards to form a ‘dip’; *t* = 7, the cell slows down and becomes unstable, as myosin is now close enough to cell front to be able to pull it in as well; *t* = 9, loss of axial symmetry, as the lower part of the cell with steeper myosin gradients is pulled inwards faster than the upper one; *t* = 23.5: emergence of stable asymmetric myosin distribution and cell shape, as the cell locks in rotations (see Fig 1*C*). (*b*) Cell shape and boundary velocities in steady rotations. Positions of the cell boundary and centroid at *t* = 23.5, 23.6, and 24 (solid, dashed, and dotted-dashed contours, respectively, and filled circles with larger size corresponding to later time). Faster boundary velocities (arrows) in the high curvature region, consistent with the location of steep myosin gradients (panel (*a*)), ensure rotational motility with a circular trajectory of the centroid (dotted arc), see also Fig 1*C*.

Once the axial symmetry of cell shape and the position of the myosin aggregate is broken (Fig 7*A*, *t* = 23.5), the boundary velocity field becomes asymmetric as well (Fig 7*B*). It is then clear that a steady movement of a cell with an asymmetric shape and asymmetric boundary velocities (where faster displacements occur at the location of higher myosin gradients) must involve rotations. Indeed, by connecting consecutively the ends of the arrows representing normal displacements of points of the boundary in Fig 7*B*, one recovers the same contour in a rotated position, as the centroids of the cell in different positions belong to the same circle (Fig 7*B*). We also note that the shape of the expanding portion of the cell boundary is reminiscent of spirals described by more abstract models of rotating free boundaries [36].

## Discussion

In this paper we systematically explored the ability of a minimal actin-myosin contractility model [7, 14] to reproduce observed mechanical states of the simplest motile cell. The model analysis has shown that the mechanical state of the cell critically depends on just three dimensionless parameters representing the myosin contractility, characteristic viscosity-adhesion length, and actin growth. For the large viscosity-adhesion length, the actin network becomes effectively stiff and does not allow for significant actin flows, which makes the cell more symmetric and as a consequence less motile and slower, and in the limit of very large values of viscosity-adhesion lengths, the cell is stationary. In the opposite limit of short viscosity-adhesion lengths, myosin forms a very small high-density aggregate, which affects the actin network only locally. In this regime, the steady motility is impossible, and the cell starts to pivot. Thus, an important conclusion is that to move straight and steadily, the cell has to keep the viscosity-adhesion length on the order of unity (to adjust the ratio of actin viscosity to adhesion strength so that it is on the order of the cell area). Interestingly, this conclusion is consistent with estimates based on the experimental data for keratocytes [5].

Intuitively, if myosin contractility is weak, the myosin spreads uniformly and the cell remains stationary and symmetric. Above a contractility threshold, the cells become motile. The mode of motility depends on the boundary conditions. For the zero actin velocity at the boundary and the sufficiently small actin growth constant and cell speed, the convex-shaped cell maintains unidirectional motility, with myosin concentrated in a band at the rear edge. With the rate of actin growth above a certain threshold, the increase of the cell speed is sufficient for the cell to lose its planar axial symmetry and start rotating. With the zero-stress boundary conditions, rotations occur for intermediate contractility strengths, whereas in the high contractility range, the fast cells stabilize their unidirectional movement, as myosin being effectively compressed into a long band at the rear edge. We found that both explored boundary conditions explain general features of the keratocyte motility, but there are interesting differences in the predicted behaviors, as discussed above.

The main finding of our study is that the contractile mechanism of motility results in a very robust turning behavior of the cell: in the models with both explored boundary conditions, the cell moving along a circular trajectory is not an anomaly but rather a solution that exists in a large region of the model parameter space. Broadly speaking, the cell starts turning in conditions of breaking the planar axial symmetry of its myosin distribution; in the ZV model the transition to rotation is controlled by the rate of actin growth, where in the ZS model–by all three independent model parameters. Turning motile behavior is an important part of the cell mechanical response in chemotaxis [37] and galvanotaxis [38], and our model generates intuition about the turning mechanism.

One important test of our model is that the solutions exhibit a characteristic fan-like keratocyte shape, with the side-to-side distance greater than the front-to-rear distance and aspect ratio between 1 and 3, in excellent agreement with the observations [5]. Moreover, the predicted aspect ratio in the ZS version is nontrivial and biphasic, reaching a maximum at intermediate myosin contractility and decreasing at very weak or strong contractility, indeed observed in [5]. Similarly, the model predicts that the lamellipodial area increases at higher adhesion and lower myosin contractility, and that the cell speed increases with myosin contractility, as observed [5]. Lastly, in agreement with the experiment [7], higher myosin contractility and/or lower adhesion strength are predicted to promote the cell polarization and motility initiation. One significant prediction of our model is that both self-polarization of the cell and its turning behavior can, in principle, occur, without complex adhesion dynamics. While it was observed that a nonlinear stick-slip adhesion behavior accompanied cell polarization [7], it remains an open question whether this nonlinear behavior is essential.

The model predicts that the stable motile behavior of the cell requires tight regulation of the total lamellipodial area. We hypothesized that this regulation is mechanical, through the membrane tension. Indeed, perturbations of the total membrane area and membrane tension were found to change the lamellipodial area in a predictable way, and drastic perturbations destabilized the cell [29]. We find that cell polarization may not depend on the cell ability to move: in the ZV model, myosin distribution and actin flow become asymmetric even in the stationary symmetric domain. However, for the motility initiation, protrusion of the boundary is obviously essential (note that while motility in the ZS free-boundary model with sufficiently high *μ*_{tot} can be initiated even with *v*_{0}→0, the area-preserving term of *V*_{p} in this limit effectively induces protrusion of the front, which is less affected by myosin, see sections *Model* and *Cell becomes motile when myosin contractility is higher than critical*).

Keratocyte motility and especially the peculiar and steady shape of the moving cell inspired a great deal of free-boundary modeling in the past decade. Our model is based on the well-justified assumption that the mechanical force balance determines cell shape and movements. Conceptually, our model is similar to the active gel models [24, 39], originating in the soft-matter physics. Some models were based on the viable idea that certain self-organized chemical patterns are upstream from the actin-myosin machinery [22, 40, 41], but majority of studies explored mechanical models [42–44]. A variety of numerical techniques–Potts models [40, 45], phase-field method [42, 44], immersed boundary method [43]–were used in respective simulations. Keratocyte polarization was modeled in [7, 46]. Alternative turning mechanisms, very different from the one predicted by our model, were computationally explored in [47, 48]. The fact that the majority of the models reproduce the keratocyte shapes and motile behavior corroborates the existing biological intuition about the keratocyte lamellipodium as the most basic, streamlined and robust actin-myosin motile structure [49]. Each of the cited studies added invaluable insights to understanding multifaceted aspects of cell motility; a relative advantage of our model is in that it is most easily connected to the experimentally observed biophysics of force balance and myosin transport in keratocytes [5, 7].

The minimal model we explored already predicts a wealth of motile behaviors (Fig 8). It is known, however, that even the cell as streamlined for locomotion as keratocyte has complexities that far exceed our minimal model. The two main aspects that need to be added to the model to make it more realistic are: spatially graded actin polymerization independent of myosin [10] and dynamic nonlinear adhesions. Complex effects of dynamic and non-homogeneous adhesions already attracted special attention and were simulated in [7, 46, 50]. It will also be interesting to explore how the predicted cell dynamics depend on actin density [51], more complex constitutive relations for the actin-myosin stress [52], membrane curvature [53, 54], elastic [55] and anisotropic [49] effects in the actin network. Our minimal free-boundary model might be useful for future modeling of other modes of cell motility [56] and collective cell movements [57, 58]. Lastly, for decades, research focused on understanding cell movements on flat 2D surfaces, and only recently exploration of cell crawling through three-dimensional (3D) matrices, more physiologically relevant, has begun experimentally [59] and theoretically [60–62]. Extension of our model to 3D will be a challenging, yet necessary, effort.

## Acknowledgments

M.N., I.L.N. and B.M.S. thank Leslie Loew for his continuing support and helpful discussions.

## References

- 1.

Ridley AJ, Schwartz MA, Burridge K, Firtel RA, Ginsberg MH, Borisy G, et al. Cell migration: integrating signals from front to back. Science 2003 Dec 5;302(5651):1704–9. pmid:14657486 - 2.

Lämmermann T, Sixt M. Mechanical modes of ‘amoeboid’ cell migration. Curr Opin Cell Biol. 2009 Oct;21(5):636–44. pmid:19523798 - 3.

Verkhovsky AB. The mechanisms of spatial and temporal patterning of cell-edge dynamics. Curr Opin Cell Biol. 2015 Oct;36: 113–21. pmid:26432504 - 4.

Keren K, Pincus Z, Allen GM, Barnhart EL, Marriott G, Mogilner A, et al. Mechanism of shape determination in motile cells. Nature 2008 May 22;453(7194):475–80. pmid:18497816 - 5.

Barnhart EL, Lee KC, Keren K, Mogilner A, Theriot JA. An adhesion-dependent switch between mechanisms that determine motile cell shape. PLoS Biol. 2011 May;9(5): e1001059. pmid:21559321 - 6.

Yam PT, Wilson CA, Ji L, Hebert B, Barnhart EL, Dye NA, et al. Actin-myosin network reorganization breaks symmetry at the cell rear to spontaneously initiate polarized cell motility. J Cell Biol. 2007 Sep 24;178(7):1207–21. pmid:17893245 - 7.

Barnhart EL, Lee KC, Allen GM, Theriot JA, Mogilner A. Balance between cell-substrate adhesion and myosin contraction determines the frequency of motility initiation in fish keratocytes. Proc Natl Acad Sci USA 2015 Apr 21;112(16):5045–50. pmid:25848042 - 8.

Ream RA, Theriot JA, Somero GN. 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. 2003 Dec;206(24):4539–51. - 9.

Verkhovsky AB, Svitkina TM, Borisy GG. Network contraction model for cell translocation and retrograde flow. Biochem Soc Symp. 1999;65: 207–22. pmid:10320940 - 10.

Ofer N, Mogilner A, Keren K. Actin disassembly clock determines shape and speed of lamellipodial fragments. Proc Natl Acad Sci USA 2011 Dec 20;108(51):20394–9. pmid:22159033 - 11.

Holmes WR, Edelstein-Keshet L. A comparison of computational models for eukaryotic cell shape and motility. PLoS Comput Biol. 2012;8(12): e1002793. pmid:23300403 - 12.

Gracheva M, Othmer H. A continuum model of motility in amoeboid cells. Bull Math Biol 2004;66: 167–193. pmid:14670535 - 13.

Carlsson AE. Mechanisms of Cell Propulsion by Active Stresses. New J Phys. 2011;13–073009. - 14.

Recho P, Putelat T, Truskinovsky L. Contraction-driven cell motility. Phys Rev Lett. 2013 Sep 6;111(10):108102. pmid:25166712 - 15.

Kimpton LS, Whiteley JP, Waters SL, King JR, Oliver JM. Multiple travelling-wave solutions in a minimal model for cell motility. Math Med Biol. 2013;30(3):241–72. pmid:22789545 - 16.

Satulovsky J, Lui R, Wang YL. Exploring the control circuit of cell migration by mathematical modeling. Biophys J. 2008 May 1;94(9):3671–83. pmid:18199677 - 17.

Hecht I, Skoge M, Charest P, Ben-Jacob E, Firtel R, Loomis WF, et al. Activated membrane patches guide chemotactic cell motility. PLoS Comput Biol 2011;7: e1002044. pmid:21738453 - 18.

Lomakin AJ, Lee KC, Han SJ, Bui DA, Davidson M, Mogilner A, et al. Competition for actin between two distinct F-actin networks defines a bistable switch for cell polarization. Nat Cell Biol. 2015 Nov;17(11):1435–45. pmid:26414403 - 19.

Raynaud F, Ambühl ME, Gabella C, Bornert A, Sbalzarini IF, Meister JJ, et al. Minimal model for spontaneous cell polarization and edge activity in oscillating, rotating and migrating cells. Nature Physics 2016 January;12: 367–373. - 20.

Lee J, Ishihara A, Theriot JA, Jacobson K. Principles of locomotion for simple-shaped cells. Nature. 1993 Mar 11;362(6416):167–71. pmid:8450887 - 21.

Rubinstein B, Jacobson K, Mogilner A. Multiscale two-dimensional modeling of a motile simple-shaped cell. Multiscale Model Sim 2005; 3:413–439. - 22.

Wolgemuth CW, Stajic J, Mogilner A. Redundant mechanisms for stable cell locomotion revealed by minimal models. Biophys J 2011;101: 545–553. pmid:21806922 - 23.

Rubinstein B, Fournier MF, Jacobson K, Verkhovsky AB, Mogilner A. Actin-Myosin Viscoelastic Flow in the Keratocyte Lamellipod, Biophys J. 2009 Oct 7;97(7):1853–1863. pmid:19804715 - 24.

Bois JS, Jülicher F, Grill SW. Pattern formation in active fluids. Phys Rev Lett. 2011 Jan 14;106(2):028103. pmid:21405254 - 25.

Foster PJ, Fürthauer S, Shelley MJ, Needleman DJ. Active contraction of microtubule networks. eLife. 2015 Dec 23;4: e10837. pmid:26701905 - 26.

Keller EF, Segel LA. Model for chemotaxis. J. Theor. Biol. 1971;30:225–234. pmid:4926701 - 27.

Leduc C, Padberg-Gehle K, Varga V, Helbing D, Diez S, Howard J. Molecular crowding creates traffic jams of kinesin motors on microtubules. Proc Natl Acad Sci USA 2012 109(16):6100–6105. pmid:22431622 - 28.

Bruggeman DAG. Calculation of various physical constants in heterogeneous substances. I. Dielectric constants and conductivity of composites from isotropic substances. Annalen der Physik 1935 24:636–664. - 29.

Lieber AD, Yehudai-Resheff S, Barnhart EL, Theriot JA, Keren K. Membrane tension in rapidly moving cells is determined by cytoskeletal forces. Curr Biol. 2013 Aug 5;23(15):1409–17. pmid:23831292 - 30.

Novak IL, Slepchenko BM. Conservative algorithm for parabolic problems in domains with moving boundaries. J Comput Phys. 2014 Mar 20; 270:203–13. pmid:25067852 - 31.

Ferziger JH, Peric M. Computational Methods for Fluid Dynamics, 3^{rd}ed. Springer; 2002. - 32.

Sibson R., A brief description of natural neighbor interpolation. In: Barnett V, editor. Interpreting Multivariate Data. Chichester: John Wiley & Sons; 1981. pp. 21–36. - 33.

Slepchenko BM, and LM Loew. Use of Virtual Cell in studies of cellular dynamics. Int. Rev. Cell Mol. Biol. 2010; 283:1–56. - 34.

Du J, Fix B, Glimm J, Jia X, Li X, Li Y, et al. A simple package for front tracking. J. Comput. Phys. 2006; 213:613–628. - 35.

COMSOL Multiphysics. Version 5.2 [software]. Stockholm, Sweden: COMSOL AB; 2015. Available from: www.comsol.com. - 36.

Altschuler DJ, Altschuler SJ, Angenent SB, Wu LF. The zoo of solitons for curve shortening. Nonlinearity 2013; 26:1189–1226. - 37.

Yang HW, Collins SR, Meyer T. Locally excitable Cdc42 signals steer cells during chemotaxis. Nat Cell Biol. 2016 Feb;18(2):191–201. pmid:26689677 - 38.

Allen GM, Mogilner A, Theriot JA. Electrophoresis of cellular membrane components creates the directional cue guiding keratocyte galvanotaxis. Curr Biol. 2013 Apr 8;23(7): 560–8. pmid:23541731 - 39.

Tjhung E, Marenduzzo D, Cates ME. Spontaneous symmetry breaking in active droplets provides a generic route to motility. Proc Natl Acad Sci USA 2012 Jul 31;109(31):12381–6. pmid:22797894 - 40.

Mareé AF, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L. Polarization and movement of keratocytes: a multiscale modelling approach. Bull Math Biol 2006;68(5):1169–1211. pmid:16794915 - 41.

MacDonald G, Mackenzie JA, Nolan M, Insall RH. A computational method for the coupled solution of reaction-diffusion equations on evolving domains and manifolds: Application to a model of cell migration and chemotaxis. J Comput Phys. 2016 Mar 15;309: 207–226. pmid:27330221 - 42.

Shao D, Rappel WJ, Levine H. Computational model for cell morphodynamics. Phys Rev Lett 2010;105: 108104. pmid:20867552 - 43.

Vanderlei B, Feng J, Edelstein-Keshet L. A computational model of cell polarization and motility coupling mechanics and biochemistry. Multiscale Model Simul 2011;9: 1420–1443. pmid:22904684 - 44.

Ziebert F, Swaminathan S, Aranson IS. Model for self-polarization and motility of keratocyte fragments. J R Soc Interface 2012 May 7;9(70):1084–92. pmid:22012972 - 45.

Nishimura SI, Sasai M. Modulation of the reaction rate of regulating protein induces large morphological and motional change of amoebic cell. J Theor Biol. 2007 Mar 21;245(2):230–7. pmid:17113108 - 46.

Ziebert F, Aranson IS. Effects of adhesion dynamics and substrate compliance on the shape and motility of crawling cells. PLoS One 2013 May 31;8(5): e64511. pmid:23741334 - 47.

Mogilner A, Rubinstein B. Actin disassembly ‘clock’ and membrane tension determine cell shape and turning: a mathematical model. J Phys Condens Matter. 2010 May 19;22(19):194118. pmid:20559462 - 48.

Camley BA, Zhao Y, Li B, Levine H, Rappel WJ. Crawling and turning in a minimal reaction-diffusion cell motility model: coupling cell shape and Biochemistry. Phys. Rev. E 2017;95: 012401. pmid:28208438 - 49.

Rafelski SM, Theriot JA. Crawling toward a unified model of cell mobility: spatial and temporal regulation of actin dynamics. Annu Rev Biochem. 2004;73: 209–39. pmid:15189141 - 50.

Shao D, Levine H, Rappel W. Coupling actin flow, adhesion, and morphology in a computational cell motility model. Proc Natl Acad Sci USA 2012;109: 6851–6856. pmid:22493219 - 51.

Kuusela E, Alt W. Continuum model of cell adhesion and migration. J Math Biol. 2009 Jan;58(1–2):135–61. pmid:18488227 - 52.

Linsmeier I, Banerjee S, Oakes PW, Jung W, Kim T, Murrell MP. Disordered actomyosin networks are sufficient to produce cooperative and telescopic contractility. Nat Commun. 2016 Aug 25;7: 12615. pmid:27558758 - 53.

Kabaso D, Shlomovitz R, Schloen K, Stradal T, Gov N. Theoretical model for cellular shapes driven by protrusive and adhesive forces. PLoS Comput Biol 2011;7: e1001127. pmid:21573201 - 54.

Du X, Doubrovinski K, Osterfield M. Self-organized cell motility from motor-filament interactions. Biophys J. 2012 Apr 18;102(8):1738–45. pmid:22768929 - 55.

Zemel A, Safran SA. Active self-polarization of contractile cells in asymmetrically shaped domains. Phys Rev E Stat Nonlin Soft Matter Phys. 2007 Aug;76(2 Pt 1):021905. - 56.

Strychalski W, Guy RD. Intracellular Pressure Dynamics in Blebbing Cells. Biophys J. 2016 Mar 8;110(5):1168–79. pmid:26958893 - 57.

Camley BA, Zhang Y, Zhao Y, Li B, Ben-Jacob E, Levine H, et al. Polarity mechanisms such as contact inhibition of locomotion regulate persistent rotational motion of mammalian cells on micropatterns. Proc Natl Acad Sci USA 2014 Oct 14;111(41):14770–5. pmid:25258412 - 58.

Löber J, Ziebert F, Aranson IS. Collisions of deformable cells lead to collective migration. Sci Rep. 2015 Mar 17;5: 9172. pmid:25779619 - 59.

Petrie RJ, Yamada KM. Multiple mechanisms of 3D migration: the origins of plasticity. Curr Opin Cell Biol. 2016 Oct;42: 7–12. pmid:27082869 - 60.

Herant M, Dembo M. Form and function in cell motility: from fibroblasts to keratocytes. Biophys J 2010;98: 1408–1417. pmid:20409459 - 61.

Tjhung E, Tiribocchi A, Marenduzzo D, Cates ME. A minimal physical model captures the shapes of crawling cells. Nat Commun. 2015 Jan 21;6: 5420. pmid:25607536 - 62.

Zhu J, Mogilner A. Comparison of cell migration mechanical strategies in three-dimensional matrices: a computational study. Interface Focus. 2016 Oct 6;6(5):20160040. pmid:27708764