article: 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.
  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.
  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, L=A0/π (e. g. , if this cell is a circle, L is its radius). We further choose L2D−1 and M0 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 = XL−1, u = ULD−1, m = M/M0, ∇ ≡ ∇x = L∇X, |ω|=|Ω|L−2 and a0 = A0L−2. Correspondingly, Eqs (1) takes the form, ∂tm=∇⋅ (deff∇m−ueffm) αΔu+β∇m−u=0, (2) where ueff = u (1−m/mmax, u), if m<mmax, u, and zero otherwise, and deff = 1−m/mmax, d. Eqs (2) include two dimensionless parameters: the dimensionless viscosity-adhesion length parameter α = ηL−2ξ−1 and the myosin contractility constant β = σM0 (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 Eqs (2), a diffusion-advection equation for myosin, is subject to the mass-conserving zero-flux boundary condition, n ⋅ (−deff∇m+ (ueff−vf) m) |∂ω (t) = 0, yielding an additional dimensionless parameter μtot=∬ω (t) m⋅d|ω|. The dimensionless boundary conditions for the force-balance equation (the second of Eqs (2) ) in the ZV and ZS models are u|∂ω (t) = 0 and n⋅ (α∇u+βmI^) |∂ω (t) =0, respectively. The dimensionless boundary velocity equation is vf = vpn+u|∂ω (t). In this equation, the dimensionless rate of membrane displacement caused by the actin polymerization and area preserving factors is vp = v0a0/a−k (a−a0 (a0/a) n), where v0 = V0L/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: vp = v0a0/ (a (1+m|∂ω (t) ) ) −k (a−a0). Note that varying parameter k = KL3D−1 is equivalent to rescaling the actin polymerization constant v0. 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 v0. Steady dynamics of a motile cell were explored by solving Eqs (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): x2+y2≤1} (such that a0 = |ω (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 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 L2-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 mk and uk are the variable values at the kth time step, mnk+1 and unk+1 are the nth iterates for the (k+1) th time step, MaxNumIters is the maximum allowed number of iterations, and ‖⋅‖∞ denotes the L∞-norm. set m1k+1=mk and u1k+1=uk for n = 1: MaxNumIters - solve αΔun+1k+1+β∇mnk+1−un+1k+1=0 to get un+1k+1 - evaluate ueff and deff using un+1k+1 and mnk+1 - solve (mn+1k+1−mk) /Δt=∇⋅ (deff∇mn+1k+1−ueffmn+1k+1) to get mn+1k+1 - calculate absolute errors ‖un+1k+1−unk+1‖∞ and ‖mn+1k+1−mnk+1‖∞ - if solution converged, break the loop, else mnk+1=mn+1k+1, unk+1=un+1k+1 end of segregated loop if n< MaxNumIters mk+1=mn+1k+1, uk+1=un+1k+1, 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.
  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 mmax, 1, mmax, 2 → ∞, so that deff = 1 and ueff = u (x, t). Then, Eqs (2) become αuxx + βmx−u = 0, mt = mxx− (um) x. In this model, varying the contractility constant β is equivalent to rescaling μtot. Indeed, β could be excluded altogether by renormalizing m: m˜=β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) = δu0 exp (λt+iqx), δm (x, t) = δm0 exp (λt+iqx) with 0<δu0<<1,0<δm0<<1 and q = π, 2π, …, so that u = δu and m = μtot+ δm. The perturbations satisfy the linearized system of differential equations, αδuxx+βδmx−δu = 0, δmt = δmxx−μtot (δu) x, and the corresponding system of linear algebraic equations, − (1+αq2) δu0+iqβδm0 = 0, iqμtotδu0+ (λ+q2) δm0 = 0, yields nontrivial solutions for δu0 and δm0, if λ (q) = q2 (βμtot (1+αq2) −1−1). The perturbations grow if λ (q) >0, with the fastest growing mode having the minimal wave number, qmin = π, and thus involving a large-scale redistribution of myosin. Thus, the symmetric state becomes unstable if βμtot>1+π2α or, in the dimensional form, σM0πL2>D (ξL2+π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 σM0L2>π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 σM0>ξ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]. Finally, it is worth noting that whereas the motility threshold in the ZS model is independent of v0, the critical values of μtot in the ZV model, where the actin growth is affected by myosin, vary with v0 (Fig 2). Indeed, the cell described by the ZV model with v0 = 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 v0, which drop with the increase of μtot. In the ZS model, the cell with a sufficiently high μtot initiates motility even as v0→0, because the asymmetric myosin pulls the boundary inward asymmetrically, and the area-preserving term causes the effective protrusion. 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 1B and S1 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 1A), 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 v0. 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 v0, 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 v0, which can be understood by noting that the cell area in this model increases with v0, thus mitigating the effect of myosin. In the ZV model, the cell speed is virtually insensitive to μtot, for μtot>π, because the term with v0 in the expression for vp becomes inessential for m∼μtot/π>1. For the same reason, the cell area is also insensitive to μtot. 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. 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 v0 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 v0, 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 6B, 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 6B, t = 15). As a result, a stable asymmetric cell shape emerges as the cell locks in rotations (Fig 6A), with myosin aggregated at a high-curvature portion of the boundary. 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 v0 threshold for rotations in the ZS model strongly depends on μtot and α (Fig 2). In particular, the rotational states may exist for any v0, if α is sufficiently small. Note also, that the concavity of the cell shape does not always destabilize unidirectional motility in the ZS model (Fig 4B). 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 4B). 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 7A, 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 7A, t = 9). Once the axial symmetry of cell shape and the position of the myosin aggregate is broken (Fig 7A, t = 23. 5), the boundary velocity field becomes asymmetric as well (Fig 7B). 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 7B, one recovers the same contour in a rotated position, as the centroids of the cell in different positions belong to the same circle (Fig 7B). 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].
  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, whereas 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 v0→0, the area-preserving term of Vp 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.
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 account 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.
keywords: velocity
  cell motility
  classical mechanics
  built structures
  engineering and technology
  cell processes
  geometry
  aspect ratio
  simulation and modeling
  molecular motors
  actin motors
  mathematics
  membrane structures
  motor proteins
  research and analysis methods
  contractile proteins
  actins
  actin polymerization
  membrane technology
  proteins
  physics
  biochemistry
  cytoskeletal proteins
  cell biology
  myosins
  biology and life sciences
  structural engineering
  physical sciences
  motion
title: A free-boundary model of a motile cell explains turning behavior