## Abstract

A very simple, one-dimensional, discrete, autonomous model of cell crawling is proposed; the model involves only three or four coupled first-order differential equations. This form is sufficient to describe many general features of cell migration, including both steady forward motion and oscillatory progress. Closed-form expressions for crawling speeds and internal forces are obtained in terms of dimensionless parameters that characterize active intracellular processes and the passive mechanical properties of the cell. Two versions of the model are described: a basic cell model with simple elastic coupling between front and rear, which exhibits stable, steady forward crawling after initial transient oscillations have decayed, and a poroelastic model, which can exhibit oscillatory crawling in the steady state.

## 1. Introduction

Cell crawling is important in biological processes as diverse as embryogenesis, wound healing and metastasis [1]. The detailed mechanical and biochemical processes involved in cell crawling remain active topics of research, although the basic components of the cycle are well known [2]. One puzzle that remains is why some cells exhibit steady motion (‘gliding’), and other cells exhibit oscillatory behaviour (‘stepping’) [3].

In a seminal paper on modelling cell migration, DiMilla *et al*. [4] describe a viscoelastic cell subject to internal cytoskeletal forces and dynamic interaction with the substrate. The study explicitly considers cells in mammalian tissue and blood, such as leukocytes or fibroblasts. Asymmetric adhesivity (differential adhesion at the front and rear of the cell) is an important feature of the DiMilla model. Numerical simulations of a six-compartment discretized representation of the cell were performed to obtain predictions of the effects of substrate adhesiveness and cell properties, such as contractile force, stiffness and viscosity, on cell speed. Although no time series are shown, the behaviour is described as cyclical, or oscillatory.

Other studies [3,5,6] have focused on the fish keratocyte, which may exhibit steady gliding motions. The recent modelling paper by Herant & Dembo [7] describes a poroelastic continuum mechanical model of cell locomotion that can exhibit impressive simulations of either oscillatory (‘fibroblast-like’) motion or smooth (‘keratocyte-like’) motion. In fact, keratocytes exhibit oscillatory, jerky motion under some conditions, as shown by Barnhart and colleagues [3]. A series of recent papers [8–11] have supported the notion that poroelasticity (the effects of moving fluid within a porous elastic solid) may be important in modulating the pressure distribution within a cell. Poroelasticity introduces time delays in transmission of mechanical signals between parts of the cell, which could lead to instability and oscillations.

The primary objective of the current study is to illuminate the factors that determine whether a cell exhibits steady or oscillatory crawling. A related objective is to clarify the essential features of autonomous cell crawling. We use a simple discrete representation of the cell, which moves in one dimension. Though the cell is a three-dimensional system, its properties are assumed to be uniform in the direction perpendicular to its motion and parallel to the substrate. This one-dimensional model retains the basic ideas of earlier models [3,4,7] but in a minimal form. The most important contributions of this study are: (i) a simple free-body analysis of the fundamental force balances involved in cell crawling; (ii) mathematical expressions for the effects of physical parameters on cell crawling speed and stability; (iii) a hypothesis for the mechanism of oscillatory cell locomotion characterized by alternating contraction at the rear of the cell and extension at the leading edge. In particular, this type of oscillatory locomotion can be explained by the effects of poroelasticity on the mechanism of cell crawling. Equations and graphical results are supplemented by animations of the crawling cell showing the motion of the leading and trailing sections, and the treadmilling of the actin network.

## 2. Overview of the mathematical model

### 2.1. Growth and contraction

Cell crawling fundamentally involves two active processes—actin polymerization in the leading edge of the lamellipodium (growth zone), and myosin-dependent contraction and actin depolymerization at the rear of the lamellipodium (the contractile zone) [4,12]. These processes are characterized by rates of growth at the front of the cell (*g*_{f}) or contraction at the back of the cell (*g*_{b}), respectively, and produce force between different parts of the cell. Let *v*_{f} be the (forward) speed of the leading edge, *v*_{b} the (forward) speed of the trailing cell body and *v*_{m} the (retrograde) speed of the actin network, all with respect to the laboratory frame and the substrate. (The subscripts f, b and m refer to *front*, *back* and *middle*, respectively.) If these variables are positive in the assumed directions, the growth rates satisfy the following kinematic relationships (figure 1*a*,*b*):
2.1and
2.2

The polymerization process produces a force *f*_{p} pushing the cell membrane forward from the front of the F-actin network, and the contraction process produces a tensile force *f*_{c} pulling the cell body towards the rear of the actin network.

### 2.2. Adhesion

Force is transmitted to the substrate via adhesive connections from the F-actin network and cell body. Adhesive force is roughly proportional to, and opposite in direction to, the velocity of the cytoskeletal network or cell body relative to the substrate [13]. Recall that behind the zone of polymerization and in front of the contractile zone, the actin network moves in a retrograde direction. Let *f*_{m} represent the forward force on the retrograde-moving actin meshwork and *f*_{b} the backward force on the forward-moving cell body (figure 1*b*). Any force on the cell in front of the polymerization zone is several orders of magnitude smaller than these [14] and is neglected. Force equilibrium for the cell requires that the vector sum of external forces is zero:
2.3

The assumed force–velocity equations are 2.4and 2.5

### 2.3. Pressure feedback

Mechanical feedback is assumed to couple the active growth and contractile processes. Contraction in the rear will produce increased hydrostatic pressure, which is transmitted to the front of the cell. This pressure acting on the cell membrane at the leading edge will decrease the force opposing actin polymerization, which will increase the rate of polymerization. In the mathematical model, hydrostatic pressure is represented by the compressive force in a spring that is shortened when contraction exceeds growth (figure 1) 2.6

Here *x*_{b} and *x*_{f} are the positions of the trailing cell body and the leading edge of the lamellipodium, respectively. Including the pressure force, the requirements for equilibrium of internal forces between sections in this model are (figure 1*c*)
2.7and
2.8

The rates of growth and contraction are assumed to be governed by simple first-order equations with time constants *t*_{f} (growth) and *t*_{b} (contraction). The model for polymerization-driven growth at the front of the cell is a coarse-grained, linearized version of the ‘Brownian ratchet’-based models [15–17]. In the current model, the rate of polymerization *g*_{f} will approach a target rate, *g*_{f0}, which depends on the polymerization force *f*_{p} [15–18]. At zero polymerization force, the target growth rate is the ‘no-load rate’, *g*_{f00}. At the ‘stall force’, , the target growth rate is zero. In between the ‘no-load’ and ‘stall’ conditions, the rate is linearly dependent on force (figure 2). Analogous equations govern the rate of contraction at the rear of the cell [19]. Physically, the force–velocity relationship in the contractile zone reflects the dependence of myosin step rate or step size on force. Direct evidence of an inverse relationship between myosin contractile speed and force is given by Aratyn-Schaus *et al*. [20]; previous studies [21,22] have shown that the contractile mechanism slows when adhesiveness of the substrate is increased. Accordingly, the contraction rate, *g*_{b}, at the rear of the cell and the growth rate, *g*_{f}, at the front are assumed to be governed by
2.9and
2.10where
2.11and
2.12

Pressure provides the basic feedback mechanism: if the contraction speed at the back of the cell is greater than the speed of growth at the front (*g*_{b} > *g*_{f}), the length of the cell will decrease, the pressure (*f*_{k}) will increase and the force resisting polymerization (*f*_{p}) will decrease. The decreased resistance will allow polymerization to proceed faster, tending to correct the imbalance between the speeds of the back and front of the cell.

## 3. Basic model: cell crawling regulated by immediate pressure feedback from back to front

The earlier-mentioned model can be formulated as a three-state system of ordinary differential equations. Substituting the kinematic relationships, equations (2.1) and (2.2), into the expression for the pressure force equation (2.6), we obtain 3.1 3.2and 3.3

The expressions for target contraction and polymerization rates are governed by linear relationships (as long as the system is not near no-load or stall conditions): 3.4and 3.5

Additional equations are obtained from satisfying force equilibrium conditions (internal and external). Velocities of the different parts of the cell are obtained in terms of growth rates and the non-dimensional ratio, *α*, between the drag coefficients for the actin network and the cell body:
3.6
3.7
3.8and
3.9

## 4. Poroelastic model: cell crawling regulated by delayed pressure feedback from back to front

The cytosol plus actin network may behave as a poroelastic material [8–10]. In such a material, a pressure gradient may exist between the rear and front of the actin network, and there will be an accompanying delay in transmitting force (fluid pressure) to leading edge. This is illustrated in figure 3 and described by equations (4.1)–(4.7) below. The ‘pressure’ force *f*_{k1} acts on the membrane at the leading edge, and *f*_{k2} is produced at the rear of the actin network; the ‘pressure gradient’ is the difference, *f*_{k2} − *f*_{k1}. In figure 3, the cause of the pressure gradient is a viscous force proportional to the difference in the velocity of the midpoint of the springs, *v*_{d}, and , a weighted average speed of the cell:
4.1

Here (the sum of the weighting coefficients *a*_{1} + *a*_{2} = 1, and the difference *a*_{2} − *a*_{1} = *a*; *a*_{2} > *a*_{1}). These coefficients influence the magnitude of the pressure gradient, as seen later (equation (4.8) and following); their values do not correspond directly to physical parameters but can be varied to explore the behaviour of the model. The dependence of ‘pressure gradient’ on relative fluid velocity in equation (4.1) is a lumped-parameter representation of Darcy's law [11]. Very roughly, *v*_{d} represents the fluid velocity and the solid velocity of the poroelastic medium. These assumed relationships provide an additional poroelastic time constant: *t*_{p} = *c*_{d}/*k*.

The state equations of this model are 4.2 4.3 4.4and 4.5

The equations for the target contraction and growth rates (far from no-load or stall conditions) become: 4.6and 4.7

The auxiliary equations (3.6)–(3.9) remain the same, and govern the speeds in the poroelastic model as well as the basic model. If the viscous coefficient *c*_{d} is zero, the two forces *f*_{k1} and *f*_{k2} would be identical, and the basic model would be recovered.

Note that if *g*_{b} and *g*_{f} were specified, the dynamics of the pressure gradient (Δ*f* = *f*_{k2} − *f*_{k1}) would be captured by a simple first-order equation:
4.8

In the case of an instantaneous contraction: , the change in length would be , where *δ* (*t*) is the Dirac delta function, and *H*(*t*) is the Heaviside, or step, function. In response, the pressure gradient would exhibit an instantaneous rise, then exponential decay:
4.9

## 5. Dimensionless form of the basic model

Dimensionless combinations of parameters are useful in understanding the behaviour of the model, and extending its general applicability. Dimensionless growth and contraction rates are obtained by dividing the physical rates by a characteristic rate: the no-load polymerization rate: 5.1and 5.2

Dimensionless forces are defined with respect to a reference force such as the stall force for contraction. For example 5.3

Dimensionless time, *T*, is defined with respect to a characteristic time for the system; we use the time constant for contraction kinetics, so *T* = *t*/*t*_{b}.

If we write the equations of motion in terms of these quantities, the following dimensionless ratios of parameters appear: 5.4

The ratios *τ*, *γ* and *ϕ* describe the parameters (time constant, no-load rate and stall force) of the polymerization process relative to contraction/depolymerization. The parameter *β* describes the ratio of maximum adhesion force at the back of the cell to maximum contractile force. The parameter *λ* captures the relative importance of the mechanical (pressure) feedback term.

In terms of this set of dimensionless parameters, the equations for the three-state variables: *p*, *G*_{b} and *G*_{f} become
5.5
5.6and
5.7where the instantaneous target rates are
5.8and
5.9

In the linear regime (i.e. not near no-load or stall conditions), the dynamic equations can be put in matrix form (primes denote derivatives with respect to ): 5.10

## 6. Dimensionless form of the poroelastic model

For the poroelastic model, there are two dimensionless forces (*P*_{1} at the front, and *P*_{2} at the rear), two feedback ratios, *λ*_{i} (one for each spring), and a new dimensionless time constant (*τ*_{p}, the ratio of the poroelastic time constant to the time constant of depolymerization).
6.1

Incorporating these and other dimensionless parameters into equations (4.1)–(4.7), the following matrix form of the dynamical system is obtained: 6.2

## 7. Results

### 7.1. Behaviour of the basic model

Steady-state solutions of equations (5.5)–(5.7) (or equivalently equation (5.10)) can be found by setting the time derivatives to zero and solving the resulting set of algebraic equations simultaneously, for a given set of parameters. These solutions provide the steady-state values of crawling speed and internal state (contraction rate, polymerization rate and pressure ‘force’) of the cell. The steady-state response depends only on the dimensionless parameters *α*, *γ*, *ϕ* and *β*.
7.1
7.2and
7.3

Local asymptotic stability of these solutions is guaranteed if all eigenvalues of the state transition matrix have a negative real part. This can be verified directly by calculating the eigenvalues, or by using the Routh–Hurwitz criterion.

Solutions are shown as a function of the dimensionless adhesiveness ratio, *β*, in figure 4, under the assumption that the other dimensionless ratios (*α*, *γ*, *ϕ*) are constant. The monotonic decrease in crawling speed with increased adhesiveness seen in figure 4 is *not* consistent with the biphasic relationship observed in experiments [22]: low crawling speeds are observed experimentally at both low adhesiveness (small *β*) and high adhesiveness (large *β*). This behaviour can be produced by the model if polymerization rate is allowed to vary with adhesiveness as suggested by Gupton & Waterman-Storer [23]. A simple way to model this is with a sigmoidal function
7.4in which *γ*_{0} is the ‘baseline’, or maximum, dimensionless polymerization rate. The effects of adhesiveness on cell behaviour with this relationship are shown in figure 5.

The explanation for the bell-shaped curves in figure 5 can be summarized as follows. When adhesiveness is low, the coupling between adhesion and polymerization [23] leads to a low maximum polymerization rate, and thus low crawling speed. When adhesiveness is high relative to the maximum contractile force that can be exerted, the cell has difficulty releasing adhesions [22], and progress is also slow. At intermediate values, neither constraint is active and the progress of the cell forward and the F-actin network backward occurs at maximal rates.

Coupling between adhesiveness and polymerization rate also increases internal pressure for low adhesiveness (small *β*) values, because the baseline polymerization rate at the front of the cell is low. Contraction in the back, combined with initially slow growth in the front, shortens the cell and creates positive pressure. This increased positive pressure stimulates an increased growth rate at the front and slows contraction at the back, equilibrating their speeds and sustaining forward progress.

In the three-state system, the Routh–Hurwitz criterion establishes that steady-state solutions are always stable for positive (i.e. physically reasonable) values of parameters. However, the transient behaviour can include damped oscillations. Simulations of the system, performed with Runge–Kutta fourth–fifth-order integration (the *ode45* function in Matlab, The Mathworks, Natick, MA, USA), are illustrated in figure 6. The simulation is started with the cell at rest with a small positive value of contractile growth *G*_{b}. Results can also be visualized by animation of the cell (electronic supplementary material, videos S1 and S2).

Figure 6 illustrates the basic mechanism of speed regulation and coordination. Contraction at the rear is the initiating event. Contraction causes an increase in pressure, which is transmitted immediately to the front of the cell. The pressure on the cell membrane increases the rate of polymerization, causing the cell front to advance. The internal pressure decreases in response to the increased growth at the front, affecting the rate of contraction at the back. Pressure feedback ultimately brings the rates of contraction and polymerization into concordance, and the cell advances steadily. This model predicts that even if oscillations are not sustained, cells typically exhibit an initial oscillatory transient when the cell starts from rest, i.e. in response to a stimulus. Oscillatory transients are predicted over a broad range of parameters; stronger force feedback and greater adhesiveness tend to produce longer, less-damped, oscillations. (The effects of different dimensionless ratios on the speed of cell advance are the same as for the poroelastic model, and are discussed later.)

### 7.2. Behaviour of the poroelastic model

Steady-state solutions can be found by setting the time derivatives of equations (6.2) to zero, and solving the resulting linear algebraic equations. Steady-state solutions for the poroelastic model are identical to those of the basic model (the two pressure forces *P*_{1} and *P*_{2} equilibrate at the same value of *p*); however, some of these equilibria may be unstable (figure 7). Stability of these solutions can again be determined from the eigenvalues of the local Jacobian matrix, or by the Routh–Hurwitz criterion. In the poroelastic model, the stability of the steady state is not guaranteed. The eigenvalues of state transition matrix are complex with positive real parts under certain conditions (figure 8). Under these conditions, oscillations arise, leading to behaviour that looks like alternating extension–contraction. Because forces do not exceed the stall force, and polymerization and contraction saturate at their no-load values, the oscillations do not grow indefinitely, but settle into a limit cycle (figure 9).

The basic physical mechanism of speed regulation and coordination is the same for the poroelastic model as it is in the basic model, except that there is a delay between the effect of pressure feedback on contraction at the rear of the cell and the effect of pressure on growth at the front of the cell. This delay causes instability and oscillation. For example, growth at the front of the cell may be decreasing even though pressure at the back has increased. By the time the pressure has changed at the front, the growth rate may have overshot, amplifying the (delayed) pressure feedback to the rear. Instability in response to delayed feedback is a classic feature of automatic control systems [24].

The effects of different parameters on steady-state solutions of the poroelastic model and their stability are shown in figures 10–12. Figure 10 illustrates the effect of adhesiveness magnitude and distribution. As in the basic model, the cell crawling speed is a biphasic function of the adhesiveness relative to the maximum contractile force, which is characterized by *β*. As in the basic model, when adhesiveness is low, polymerization rate (and thus cell crawling speed) is reduced because of the modelled relationship between adhesion and polymerization. At high adhesiveness ratios, the cell's progress is limited by its ability to drag its rear adhesions from the substrate.

Bell-shaped solution curves analogous to those in figure 7 are shown for different values of the adhesiveness distribution ratio, *α*, in figure 10. The curves in figure 10*a*,*c* show that increasing the adhesiveness of the actin network towards the front of the cell relative to adhesiveness of the rear (larger *α*) increases the forward speed, *V*_{f}, of the cell, and decreases the retrograde speed, *V*_{m}, of the network. In effect, for large *α* the ‘clutch’ slips less. The mean growth/contraction rate, *G*, also decreases with increasing *α*, while the pressure feedback is almost unaffected. Steady solutions are unstable at higher levels of adhesiveness (larger *β*); stability is insensitive to the distribution of adhesiveness (*α*).

Increasing the baseline polymerization rate (characterized by *γ*_{0}) predictably increases the cell crawling speed (*V*_{f}, figure 11*a*) and treadmilling speed (*V*_{m}, figure 11*c*), and decreases the average internal pressure, *P* (figure 11*d*). Increasing *γ* also appears to de-stabilize crawling behaviour; at larger values of *γ*_{0}, unstable solutions appear over a greater range of adhesiveness.

Changing *ϕ*, the ratio of polymerization stall force to contractile stall force, has little effect on the average crawling speed in the range we studied. Decreasing *ϕ* tends to de-stabilize steady crawling solutions and favour oscillatory solutions.

The effects of the various parameters on stability of steady crawling are shown in the contour plots of figure 12. These boundaries define the regions in which any eigenvalue of the state transition matrix is positive. Contours are obtained for values of *ϕ* = 0.1, 0.2, … , 1.0. Unstable (U) regions are uniformly larger for smaller values of *ϕ*, suggesting that stronger contractility (or weaker polymerization) may make cells more prone to oscillatory motion. The plots in figure 12*a* confirm that higher values of adhesiveness (*β*) also favour instability, as does stronger pressure feedback (*λ*_{1}, figure 12*b*), and higher polymerization rate (*γ*_{0}, figure 12*c*). These parameters all tend to increase the feedback ‘gain’, leading to larger responses to internal signals.

Physical parameters that may be varied directly in experiments include polymerization rate (*g*_{f00}), maximum contractile force () or contraction rate (*g*_{b00}), and adhesiveness (*c*_{b}). In the scenario that *g*_{b00} is fixed, the effects of contractile stall force, adhesive force and polymerization rate are shown in figure 13. Stability boundaries as a function of physical adhesivity (*c*_{b} *g*_{b00}, with units of pN) and *γ* (the ratio of polymerization rate to baseline contraction speed) are shown in figure 13. Adhesivity can be roughly interpreted as the force needed to detach the cell body from the substrate at its nominal speed. Each boundary corresponds to a different value of contractile force, . Instability is favoured by high adhesive force, high polymerization rate and high contractile force.

## 8. Discussion

A simple autonomous model of cell crawling is presented in this paper. This highly simplified model is based on a discrete representation of the cell; adhesion is approximated by velocity-dependent frictional forces between the cell and substrate. The model consists of a set of ordinary differential equations describing polymerization at the leading edge, contraction and depolymerization at the rear of the actin network, and pressure feedback coupling these two active processes. Steady-state solutions of the model correspond to steady forward crawling with retrograde treadmilling of the actin network. These solutions, obtained in closed form, describe the well-known biphasic dependence of crawling speed on adhesiveness.

The model is augmented to model the effects of poroelasticity, as advocated by Charras, Mitchison, Mahadevan and co-workers [8–10]. In a poroelastic cell, changes that occur on a timescale *τ* will lead to significant differences in pressure between points in the cell separated by the length scale , where is the diffusion constant of the system, *k* is the hydraulic permeability, *K* is the elastic bulk modulus of the solid network and *ϕ* is the volume fraction of fluid. Charras *et al*. [9] estimated and , corresponding to a diffusion constant m^{2} s^{−1} and a length scale for pressure variations of . This length scale is comparable to cell length, indicating that pressure differences across the cell will indeed be significant. The magnitude of pressure variation is expected to be greater than 100 Pa, which is comparable to the magnitudes of the elastic and viscous stresses in the cell.

Inclusion of poroelastic effects leads to oscillatory behaviour with alternating, out-of-phase contraction at the back of the cell and extension of the leading edge. The general dependence of crawling speed on adhesiveness is preserved. The effects of the distribution of adhesiveness, the rate of polymerization, the strength of the pressure feedback term and the poroelastic time constant are readily identified from the dimensionless parameters used in the model.

Despite its coarseness, this model is consistent with a number of experimental observations. It is important to note that crawling speeds are not imposed (i.e. as an input parameter), but are obtained by solving the model equations. (i) The biphasic dependence of crawling speed on adhesiveness is preserved. (ii) For keratocytes on a given substrate, unstable (oscillatory) behaviour occurs in faster crawling cells, and oscillation frequency is correlated with cell crawling speed. (iii) Faster physical speeds are associated with high values of baseline contraction rate *g*_{b00}, which increases the non-dimensional parameters *β* (adhesive force) and *λ*_{1} (pressure feedback). Both tend to cause instability, and increasing *λ*_{1} increases oscillation frequency. (iv) Cells in which myosin activity is inhibited by blebbistatin, for example, crawl more slowly (but still crawl). This is consistent with the effects of reduced contraction rate, *g*_{b00}, in the current model.

The parameters explored in this theoretical study are mostly non-dimensional ratios of quantities that are not known precisely. Estimates of the underlying physical parameters of the current models are included in table 1. To obtain these lumped parameters, quantities reported as distributed parameters, such as stress or force per unit length, have been multiplied by estimates of the appropriate area or length. The values and interpretations of dimensionless parameters derived from physical values are summarized in table 2.

Why is a model like this useful, when more realistic and sophisticated models [4,28] are available? Representing the effects of well-understood processes into a few parameters, and making clear assumptions about others, allows hypotheses about the basic character of cell crawling to be tested. The most important prediction of this model is that the time delay in transmitting hydrostatic pressure through a poroelastic medium can produce alternating contraction/extension rather than steady crawling. The fundamental processes in steady crawling are the same as those in oscillatory crawling: polymerization, contraction and adhesion. Higher adhesiveness (relative to the maximum contractile force) higher polymerization rate (relative to contraction/depolymerization rate), lower polymerization force (relative to contraction force) and stronger pressure feedback all tend to produce oscillatory behaviour in the model. Finally, this discrete representation of the cell can be used to model and explore more complicated constitutive relationships, such as nonlinear adhesive force–velocity relationships or more sophisticated polymerization models.

## Acknowledgements

Support from NSF grant DMS-0540701 and NIH grants R01 NS070918 and R01 GM086882 is gratefully acknowledged.

- Received September 14, 2011.
- Accepted October 6, 2011.

- This journal is © 2011 The Royal Society