Computational modelling of cell motility on substrates is a formidable challenge; regulatory pathways are intertwined and forces that influence cell motion are not fully quantified. Additional challenges arise from the need to describe a moving deformable cell boundary. Here, we present a simple mathematical model coupling cell shape dynamics, treated by the phase-field approach, to a vector field describing the mean orientation (polarization) of the actin filament network. The model successfully reproduces the primary phenomenology of cell motility: discontinuous onset of motion, diversity of cell shapes and shape oscillations. The results are in qualitative agreement with recent experiments on motility of keratocyte cells and cell fragments. The asymmetry of the shapes is captured to a large extent in this simple model, which may prove useful for the interpretation of experiments.
Crawling motion of eukaryotic cells on substrates is caused by polymerization–depolymerization cycles in the actin filament network. The propulsion forces exerted by the cell on the substrate are determined by a balance of actin polymerization at the leading edge, intermittent adhesion of the cell to the substrate and detachment of the back of the cell; all of these processes are further mediated by a complex protein regulatory system [1–3]. Detailed modelling of cell motion that includes the entire complexity of all regulatory pathways is a formidable challenge; many pathways are intertwined and poorly understood, and forces that give rise to motion are not fully quantified . Moreover, cell motion is a complex dynamic process associated with cell shape deformation and complicated flow geometries. Direct numerical modelling requires the inclusion of movable deformable interfaces or thin-film theories of multi-phase reactive flow, which introduce additional computational challenges [5–11].
The accepted, though simplified, picture of cell motility is the following: actin filament polymerization occurs predominantly at the front of the cell, where Wiskott–Aldrich syndrome proteins nucleate new branches of actin filaments by activating the Arp2/3 complex located close to the cell membrane . Actin polymerization, via adhesion complexes, produces forces on the substrate and, consequently, induces cell movement. At the same time, conventional myosin molecular motors are expressed at the rear of the cell, triggering formation of actin filament bundles. The bundles exert contractile stresses at the rear and are responsible for the overall anisotropy and crescent-like shape of the moving cell. Actin filaments are concurrently disassembled into G-actin monomers in the bulk of the cell. Rapidly diffusing G-actin is then recycled at the cell membrane allowing for polymerization of new actin filaments.
In the proposed model, we use a phase-field approach based on ideas developed for solidification , crack propagation , fluidization transition in granular media [15,16] and for actin gel-bead motility . In the context of vesicles, this approach is also known as the advected-field approach . Recently, in a similar framework, Shao et al.  suggested a model for the cell shape dynamics during motion. The main advantage of this approach is that no explicit tracking of the cell's boundary is required as the auxiliary phase-field distinguishes the interior of the cell from the exterior. Shao et al.  applied this method to the motion of epithelial keratocytes , crescent-shape cells that extend a thin lamellipodium at the front and sides. While the model reproduced cell motion and shapes, it used two scalar fields, cross-linked actin filaments and actin bundles, instead of a vector field describing the orientation and intrinsic anisotropy of actin. The latter one is believed to be important for motility and has been characterized experimentally .
In this work, we propose a simple phase-field model describing the cell shape coupled to the orientation (or polarization) of the actin filament network. Two fundamental issues distinguish our work from Shao et al. : first, we include the actin filament polarization in our modelling; this is more realistic and yields additional information, e.g. of the traction on the substrate. Second, we do not need to explicitly use two separate fields (hence separate forces) for the protrusion and the retraction to sustain cell motion as was the case in the study of Shao et al. .
Our model reproduces the primary phenomenology of cell motility: we find a discontinuous onset of cell motion, as observed experimentally for cell fragments . To date, this is the first model describing this transition without pre-imposing the shape. We also obtain correct stationary crescent-like shapes of moving cells as well as periodic cell shape oscillations in the course of motion. Finally, the outcome of our modelling can be directly compared with recent experiments on actin network polarization distribution, distribution of traction on the substrate and cell shape [12,20–23].
2. Phase-field model
Our approach is based on two fields: first, the phase-field variable ρ(x, y; t) is an auxiliary field separating the interior of the cell from the exterior: ρ has a value of 1 inside the cell, 0 outside the cell and varies smoothly at the interface. This diffuse interface is interpreted as the location of the cell membrane, see below. Second, the vector field p(x, y; t) describes the polar orientation (or polarization) of the actin filaments. It points in the direction of the local mean orientation of filaments and its absolute value is a measure of the degree of ordering. The overall description is depth-averaged, i.e. we do not resolve the actual height of the cell and the model is effectively two-dimensional.
Before proceeding to the derivation of the model, let us discuss in detail the interpretation of the phase-field variable ρ. As we stated above, the phase-field approach became a legitimate mathematical tool and effective numerical alternative to the sharp interface limit for systems undergoing first-order phase transition [13–15]. It was successfully applied to a broad range of phenomena, from solidification to granular fluidization. While it was originally designed as an auxiliary field, in some cases the phase field can be associated with relevant physical variables, such as the number of broken atomic bonds in the case of crack propagation  or the ratio of static versus sliding contacts in dense granular flows [15,16]. In the context of vesicles, the membrane is not a boundary between true phases; hence, the advected-field interpretation (the membrane is advected by flows in the two domains) has been developed . It is in this spirit that we interpret the phase-field approach as the advection of the membrane (given by the diffuse interface) by forces created in the cell's interior. We should stress that the phase-field variable does not capture any of the structure of the cytoplasm; it is rather the orientation field p that resolves this missing piece by describing the average orientation and the degree of orientation of the filaments that polymerize and push against the membrane. Below (see §5), we show that a simplified sharp interface limit of our model indeed yields the correct kinematic condition, i.e. the steady-state velocity of the cell is proportional to the protrusive force diminished by a surface tension contribution owing to the local two-dimensional curvature of the membrane.
We start by writing the evolution equation for the phase-field ρ, 2.1
The parameter Dρ (with units of length squared as time is already rescaled) has two meanings: first, it characterizes the width of the phase-field interface. Second, it can be shown to be the ratio of the surface tension of the cell membrane to the friction with the substrate . As the width of the interface is inconsequential, we do not introduce two separate parameters for these two effects. It is known how to generalize this approach , a step which is needed to perform a strict sharp-interface limit (where the width goes to zero while the surface tension remains finite), on which we will report elsewhere. In the following, we also neglect the membrane's bending rigidity, which would correspond to higher order differentials in ρ [18,19]. While of possible importance in the unresolved z-direction (where the radius of curvature is 100 nm), the bending rigidity can be estimated to be negligible for typical curvature radii of the order of 10 µm involved in the two-dimensional crawling motion of a cell.
The parameter δ is introduced in equation (2.1) to control the motion of the cell boundary and to ensure the approximate (effective two-dimensional) volume or area conservation of the cell (cf. [15,16]). Equation (2.1) has three fixed points resulting from the local nonlinearity (1 − ρ)(δ − ρ)ρ: ρ = 0, ρ = 1 and ρ = δ. For 0 < δ < 1, the two fixed points ρ = 0,1 are stable and separated by the unstable fixed point ρ = δ. Correspondingly, in spatially extended systems for δ = δ0 = 1/2 the interface between the two ‘phases’ ρ = 1 and ρ = 0 is stationary. This is the so-called Maxwell rule. Increasing δ above 1/2, on the other hand, leads to an advance of the phase with ρ = 0 into the one with ρ = 1, and vice versa for δ < 1/2. For the present situation, we assume δ to be of the form 2.2
The rational for this special choice in equation (2.2) is the following: the parameter δ0 characterizes the ‘internal’ hydrostatic pressure of the cell. For δ0 < 1/2, the cell tends to expand (positive pressure), and for δ0 = 1/2, the cell shape is quasi-stationary (no pressure). The second term, ∫dx dy ρ(x,y) − V0, is an implementation of the constraint of total area conservation: if the cell's total area ∫dx dy ρ(x,y) exceeds the chosen initial area V0, the parameter δ increases, leading to overall contraction of the cell or an advance of the phase with ρ = 0 into the one with ρ = 1. The parameter μ characterizes the stiffness of this constraint, and may depend on membrane surface tension and elasticity. However, the overall behaviour is not very sensitive to variations of μ. The last term ∼|p|2 in equation (2.2) models a contractile stress due to motor-induced contraction of actin filament bundles with the dimensionless model parameter σ being a ratio of maximum tensile stress to maximum internal pressure. The bundles, comprising either parallel- or anti-parallel-oriented filaments , can be described by a nematic tensor Sij, which is related to the orientation vector p according to 2.3This equation describes the following behaviour: in the absence of polarization (p = 0), the nematic tensor relaxes towards an equilibrium value Sij = Sij0 given by filaments that are not in motor contact but still may display a certain amount of nematic ordering. τ is the relaxation time for the nematic field and Sij0 is a constant traceless tensor, which for simplicity we set to zero here. If there is finite motor-induced polarization p, the value of the nematic order parameter Sij increases. In the following, we assume that the relaxation time τ is small and that the nematic field is enslaved to the polarization , i.e. we neglect the time derivative in equation (2.3). The nematic tensor is then of the form 2.4
It is known from the theory of ‘active gels’  that the stress tensor due to motor contraction of bundles is given by σij ∼ Sij. Correspondingly, the tensor's invariant, |p|2, describes the magnitude of contractile deformations and couples cell dynamics and mechanics.
The full coupled model equations for the phase-field ρ and the orientation vector p read 2.5and 2.6The two fields ρ, p are dimensionless and 0 < ρ, |p| < 1, i.e. the orientation field is normalized by its maximum value. In equation (2.5), α describes the advection of ρ along the orientation vector p due to polymerization of actin filaments against the membrane. Its relation to the cell's velocity via a kinematic condition is discussed in §5. Also note that here cell adhesion is tacitly assumed: for the growing filaments to push, they have to transfer momentum to the substrate, see also the discussion of traction forces below.
The orientation field p is non-conserved, as it can be created by polymerization: this is described by the β-term in equation (2.6), stating that polarization is created by ∇ρ. Hence, at the location of the membrane, given in the diffuse interface description by the maximum of ∇ρ, orientation is created perpendicular to the interface, owing to polymerization of new actin filaments. On the other hand, the τ1-term in equation (2.6) describes degradation of orientation by diffusion, depolymerization, etc. The τ2-term dampens the orientation field p outside the cell while not changing the behaviour inside where ρ = 1. Note that the cell does not advect the filaments (and hence p) considered here; rather, since the polymerization exerts a force on the membrane, the respective filaments must be connected to the substrate via adhesion sites. Hence, we do not advect the orientation field but rather suppress the variable p outside the cell. We verified that this changes neither the shapes nor the dynamics. The term DpΔp describes diffusion and stiffness of the polar orientation field.
Finally, the γ-term breaks the reflection symmetry of the orientation p. This can be motivated as being an effect of another (overdamped) field, namely the concentration of myosin motors m. For the latter, we can write the following equation: 2.7where Dm, Vm and τm are the motors' diffusion, speed along the filaments and relaxation time, respectively, and m0 is the steady-state motor concentration. The motors are suppressed near the front of the cell with a rate . For small τm, we can neglect the terms DmΔm,Vmp · ∇m and write . Hence, the motor density is higher at the rear, and the corresponding term in the equation for p effectively decreases polarization as filaments are absorbed into (anti-parallel) bundles at the rear of the cell.
We solved equations (2.5) and (2.6) on a periodic two-dimensional square-shaped domain using a numerical, quasi-spectral split-step method. Typically, 256 × 256 fast Fourier transformation harmonics were used. The integration domain was chosen large enough to ensure negligible contributions from the boundaries (typically 100 × 100 dimensionless units). The initial conditions were generally a circular domain for the phase-field variable, i.e. ρ = 1 inside a circle of radius r0 (hence V0 = πr02 is the cell's area) and ρ = 0 otherwise. The orientation field was set to a constant value p = const inside the circle and pointing in the x-direction. Depending on the model parameters, the initial conditions either restored a circular symmetric shape and became stationary, or they evolved into stable polarized moving crescent-shape patterns, closely reminiscent of forms observed experimentally for keratocyte fragments [20,21].
Before discussing moving cells, a brief comparison between the obtained static cell shape and nematic droplets is in order: a similar shape (called radial structure) is found for nematic droplets with perpendicular surface anchoring (homeotropic alignment) [27,28]. Despite the similarity in overall conformation and alignment orientation, there are several important differences due to the activity of the cytoskeleton: first, a nematic system is not described by a vector field p like the polar cytoskeletal filaments (the modulus of which describes the amount of ordering ), but rather by a director field with no direction specified, , and where the extent of ordering is given by the second moment, i.e. the nematic tensor Sij. Second, in a nematic droplet the active terms, i.e. the advection term (proportional to α) and the motor-related term (proportional to γ) do not exist. Finally, there is no continuous creation of new ordered filaments at the interface (the β-term), but only passive surface anchoring.
Select shapes of moving cells obtained by solving equations (2.5) and (2.6) are shown in figure 1. While the model has a substantial number of parameters, the overall outcome of the numerical modelling is the following: self-sustained moving cells emerge from initial conditions for sufficiently large values of the parameters α, β and γ (characterizing the polymerization force, the polymerization rate and the motor-induced asymmetry, respectively). In turn, increase of the parameter σ accounting for the contribution from tensile stresses generated by actin–myosin bundles results in an increase of the aspect ratio of the cells (figure 1a–c).
Surprisingly, if the parameters α and γ are moderate but of similar order, we observed a transition from steady motion to a state displaying motion with periodic oscillations in both the velocity and the shape of the cell. Typical shapes for this oscillating moving state are shown in figure 1d–f. A time-dependent characterization is displayed in figure 2. Interestingly, the oscillations of the velocity V and the cell's aspect ratio and anisotropy are not in phase: there is approximately a quarter-period shift between the maximum velocity and the maximum of the aspect ratio. While no conclusive data are available on shape oscillations of keratocytes in vivo, periodic lamellipodium contractions of crawling fibroblasts  were reported as well as periodic motion and dynamics of filopodia on compliant substrates . Usually, a stick-slip mechanism is used to rationalize these findings. However, our simple model shows that the nonlinear coupling between the filament polarization and the cell shape is sufficient to allow for oscillatory states.
To characterize the onset of motion and the minimal propagation velocity, we performed simulations for various parameter ranges. Results are shown in figure 3 as a function of the parameter γ describing the asymmetry induced by the molecular motor distribution. We started with a stationary moving structure for a large value of γ and then slowly ramped down γ. As one can see from the figure, the cell's motion stops abruptly below a critical value γc ≈ 0.42, at a minimal velocity of Vmin ≈ 0.62, which of course both depend on the other model parameters. The velocity V, the aspect ratio h − 1 and the asymmetry measures ζ, η (see §4 for a detailed definition) drop to zero for γ < γc. Remarkably, one of the asymmetry measures, ζ, increases approaching the transition point, while the second, η, decreases monotonically.
In figure 3, we showed that a moving cell stops abruptly when the parameter γ is reduced. On the other hand, when a stationary spherical fragment is perturbed by noise it remains stationary for the parameters shown in that figure. One needs a nonlinear perturbation to trigger motion, e.g. by choosing a region inside the cell and adding a constant p inside that region. This overall behaviour shows that the onset of cellular motion is a subcritical instability, in agreement with experiments on fragments : there fragments had to be perturbed to initiate motion, while unperturbed or slightly perturbed fragments stayed immobile. This subcriticality is caused by the motor-related symmetry-breaking term (proportional to γ), and reflects the fact that a perturbation has to be large enough for the motors to establish a significant asymmetry in their distribution across the cell.
To further characterize the shapes, we have investigated the dependence of the cell's aspect ratio h − 1 and the asymmetry measures characterizing deviations from reflection symmetry, ζ, η, on the parameters and the cell's velocity V. The aspect ratio, in particular, has been used experimentally to characterize moving cells and has been related to the cell's speed. Figure 4 illustrates the dependence of the cell's aspect ratio h − 1 and the asymmetry measures ζ, η on the cell's velocity V. The simulation conditions were as shown in figure 3. The aspect ratio and η decrease with decreasing velocity, while the other asymmetry measure ζ increases.
Figure 5 features the cell's velocity V and the aspect ratio h − 1 as a function of the cell's overall area V0. While there is a monotonous and only slight increase of the velocity with cell size, the variation for the aspect ratio is high and exhibits a surprising non-monotonous behaviour with a maximum at some ‘optimal area’. Namely, as the parameter σ describing the contractility was kept constant, cells of small area attain a nearly elliptical shape. Increasing the area leads to more asymmetric and extended shapes, as seen in figure 1. For even larger volume the shape becomes extended in the direction of motion, with a comparatively small propulsive front and a large tail shaped as a conical frustum dragged behind.
Thus, the relation between the cell's velocity and its shape, characterized here both by the aspect ratio and the asymmetry measures, is far more complicated than recently suggested [19,20]. The influence of the overall volume, in particular, must be accounted for (when experimentally addressing this issue, cells with similar area should be compared).
4. Material and methods
The velocity, aspect ratio and asymmetry measures of oscillating cells have been shown in figure 2. To extract the velocity, we determined the centre-of-mass position xc according to 4.1The aspect ratio was determined via the corresponding 2 × 2 variance matrix Iij 4.2Its eigenvalues λ1,2 were calculated and their ratio 4.3is a measure for the aspect ratio of the cell. For a cell moving in x-direction (i = 1), for the off-diagonal elements follows I12 = I21 = 0 and the aspect ratio is simply given by . As for an axisymmetric shape the aspect ratio parameter is h = 1, in the figures we trace h − 1.
We also characterized the asymmetry of moving cells, namely the deviation from reflection symmetry. For this purpose, we calculated the skewness tensor Gijk 4.4Obviously, for an ellipse Gijk = 0. For an asymmetric, e.g. crescent-shape, cell moving in the x-direction (i, j, k = 1), only four elements of the skewness tensor are non-zero: G111 ≠ 0, and G122 = G212 = G221 ≠ 0. Hence we can define the following relative asymmetry measures 4.5and 4.6For an ellipse, ζ, η = 0; for the asymmetric moving cells obtained by our simulations one has ζ, η ≠ 0.
5.1. Interface dynamics
Our description is depth-averaged and does not explicitly account for the substrate. However, the propulsion of the phase-field ρ by polymerization obviously implies that filaments are attached to the substrate or to the rest of the actin gel. Hence, equation (2.5) for the phase field can be related to a simple kinematic condition at the interface. Let us consider a constant velocity solution to equation (2.5), namely ρ(x, y, t) = ρ0 (x − Vt, y), where V is the velocity. Substituting into equation (2.5), in the local coordinate system, we get for the front position 5.1where n is the coordinate normal to the front and κ is the local two-dimensional curvature. Thus, to leading order the front dynamics is described by 5.2If we multiply equation (5.2) by a viscous friction coefficient ξ, we rewrite it as a local force balance equation (cf. ): Ffr + Fpol + Ften = 0, where Ffr = ξV is a viscous friction force, Fpol = ξαpn is the force due to actin polymerization and Ften = ξDρκ is the force due to interface tension. The force balance allows us to evaluate model parameters from experimental data (see below) and to interpret the corresponding terms in equation (2.5) as force densities.
5.2. Traction distribution
The polymerization force ∼αp is balanced by the viscous drag force proportional to the cell's velocity V. The resulting traction T on the surface will be the sum of these two forces (neglecting additional elastic bulk forces and other complications). Hence, T ∼ p − νV, where ν is an effective friction coefficient. Note that the latter can be determined self-consistently from the condition that the total force applied to the surface must be zero (cell is a force dipole). In conclusion, the approximate traction field can be determined from the simulation as follows: 5.3Here we introduced the common factor ρ to ensure that the traction vanishes outside of the cell. The obtained traction distribution is shown in figure 6b. It is in qualitative agreement with experimental measurements .
5.3. Retrograde flow
A so-called retrograde flow of actin filaments from the leading edge towards the rear is observed in the course of motion for most cell types; for keratocytes, this was not experimentally observed until recently  due to a rather strong interaction between the actin network and the substrate. It seems like the retrograde flow is not of crucial importance for motility in this case, although this is at present still debated. To approximately reconstruct the retrograde flow v inside a moving cell, we apply Darcy's law known from porous media, assuming that filaments flow through meshes in the actin network. Hence, 5.4where ψ is transport coefficient. Here we assumed that the flow is generated by an effective pressure (due to polymerization at the front and depolymerization at the rear) determined by the absolute value of the orientation vector. The image shown in figure 6c is again in qualitative agreement with experimental data : the flow is retrograde and slightly increased at the sides of the cell. The flow is of the order of 10 per cent of the cell's velocity.
5.4. Estimation of model parameters
Some parameters in our model can be estimated from available experimental values: the maximum actin polymerization velocity (α, β ≈ 0.1 µm s−1), the actin depolymerization rate (1/τ1 ≈ 0.02–0.1 s−1, depending on physiological conditions of the cell) and the diffusion coefficient of actin filaments Dp ≈ 0.1–0.2 µm2 s−1 . As in our model we set 1/τ1 = 0.1, the unit of time is set approximately to 10 s. From the diffusion coefficient of p used, Dp = 0.2, we obtain the typical length scale of 1 µm, which is of the order of an actin filament. Correspondingly, the range for the dimensionless values of the parameters α, β in equations (2.5) and (2.6) is 0.2–2. The stiffness of the interface is given by the ratio of surface tension (typically 1 pN) and friction coefficient (typically several pNs µm−2)  which yields the estimate Dρ ≈ 1. The motor-related parameter γ was varied from zero to a value γ ≈ 0.1 µm s−1, which is of the order of a typical motor speed, implying that it takes at least 100 s to establish a broken symmetry across a 10 µm cell. A list of the model parameters and typical values used in the simulations are summarized in table 1.
In this work, we proposed and studied a nonlinear dynamical model for cell motility, which is computationally effective and gives insights into the dynamics and shapes of moving cells or cell fragments. The results are in overall qualitative agreement with available experimental data. Our model yields some non-trivial predictions worthy of experimental verification. In particular, our results indicate that there is no unique relation between the cell's aspect ratio and speed; the mutual dependence is influenced by the total area/volume of the cell, as well as by other model parameters, and it could be even non-monotonous as shown in figure 5. This multiparametric dependence may explain significant scatter of the data on aspect ratio versus speed plots measured experimentally (cf. fig. 4b in ).
The proposed model can be extended in many directions. For example, the effects of mechanical stresses (elasticity or viscoelasticity of the actin gel) or intracellular flow (which we discussed here only very briefly in the Darcy limit) could be incorporated into the model as additional fields directly coupled to the polarization and phase field. Another important avenue for the development is a more realistic implementation of the interaction with the substrate. The substrate's elastic properties (as shown in Chan & Oded ) strongly affect the character and dynamics of cell motion, an effect that can be possibly captured in a refined model.
Finally, in addition to insights into the biological problem of cell motility, our results could prove useful for the design of synthetic systems, such as self-propelled polymer functionalized microcapsules [31,32] that are currently investigated for the design of self-healing materials with artificial microvascular networks: their mobile healing agents (e.g. glue-filled microcapsules ) are transported along the microvascular network to damaged areas to make repairs .
We thank J. Oliver, P. Sens, E. Raphaël, J. Prost, J.-F. Joanny, F. Jülicher and K. Kruse for stimulating discussions. We also would like to thank one of the referees to point out the connection to nematic droplets. F.Z. thanks the DFG for partial funding via IRTG 1642 Soft Matter Science. S.S. acknowledges support from KAUST, the University of Oxford and the DDR&E and AFOSR under Award No. FA9550-10-1-0167. I.S.A. thanks the ESPCI for hospitality and the CNRS for support during his stay. The work of I.S.A. was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering, under Contract DE-AC02-06CH11357.
- Received July 4, 2011.
- Accepted September 29, 2011.
- This journal is © 2011 The Royal Society