Abstract
We present a mathematical and a computational framework for the modelling of cell motility. The cell membrane is represented by an evolving surface, with the movement of the cell determined by the interaction of various forces that act normal to the surface. We consider external forces such as those that may arise owing to inhomogeneities in the medium and a pressure that constrains the enclosed volume, as well as internal forces that arise from the reaction of the cells' surface to stretching and bending. We also consider a protrusive force associated with a reaction–diffusion system (RDS) posed on the cell membrane, with cell polarization modelled by this surface RDS. The computational method is based on an evolving surface finiteelement method. The general method can account for the large deformations that arise in cell motility and allows the simulation of cell migration in three dimensions. We illustrate applications of the proposed modelling framework and numerical method by reporting on numerical simulations of a model for eukaryotic chemotaxis and a model for the persistent movement of keratocytes in two and three space dimensions. Movies of the simulated cells can be obtained from http://homepages.warwick.ac.uk/∼maskae/CV_Warwick/Chemotaxis.html.
1. Introduction
Modelling the directional motility of cells is of great importance especially because of the central roledirected cell migration plays in several biological phenomena, such as embryonic development, cancer, tissue development and immune responses [1]. Broadly speaking, the motile cycle of a cell consists of the following processes: polarization where the cell develops a front and a back through the redistribution of proteins and lipids within the cell, protrusion at the leading edge of the cell pushing the front of the cell outwards, and retraction of the rear of the cell towards the leading edge [2]. Although the main aspects of the motile cycle appear deceptively simple, as further details are added to the modelling, various complexities arise. For example, in the case of chemotactic eukaryotic cells, the molecular mechanisms that govern gradient sensing and cell polarization are still not fully understood [3]. Furthermore, it is difficult to quantify the forces associated with motility and only recently has experimental progress been made in this direction [4,5]. Direct numerical simulation of cell motility necessitates the consideration of deformable surfaces [6,7] or multiphase flow models [8,9], both of which are computationally challenging. Finally, the deformation of the cell surface is linked to the dynamics of actin and other cellresident proteins, and these dynamics must be coupled in a consistent way with the evolution of the cell surface.
In this work, we present a mathematical framework for the modelling of cell motility and a numerical method for the simulation of such models. The approach we propose uses ideas of existing twodimensional models but generalizes these and extends the modelling to the threedimensional setting. It consists of partial differential equations, specifically those of reaction–diffusion type, posed on the cell boundary coupled with an evolution law for the cell membrane. Further, we discuss the inclusion of external forces and illustrate this with a phenomenological model for the interaction between a cell and a obstacle. We present a numerical method, based on evolving triangulated surfaces, that consists of an evolving surface finiteelement method [10] for the approximation of the surface partial differential equation and a parametrized finiteelement method [11] for the approximation of the surface evolution law.
It is our hope that the parametric approach we employ will be more efficient than other standard approaches such as phase field [6] or level set methods [12]. The reasoning behind this statement is that our methodology based on triangulated surfaces formulates the problem in one dimension less than other approaches, in which the equations are discretized in the embedding space. See, for example [13–16] and references therein for further discussion.
We consider two specific models for cell motility. First, a model for eukaryotic chemotaxis. Aspects of chemotaxis, such as changes in direction owing to splitting and biased generation of pseudopods as well as response to a changing chemotactic signal [17], are present in this surfacebased model. We also present simulations of pseudopoddriven migration of a cell in three dimensions. We next consider a surfacebased model for the persistent migration of fish keratocytes, presenting numerical simulations in two and three dimensions. The surfacebased model is able to capture the different shapes that characterize migrating keratocytes and the correlation between aspect ratio and cell speed [18].
A summary of the contributions of our study is as follows. We derive a rigorous mathematical framework for the modelling of cell motility and chemotaxis in two or three space dimensions, our modelling includes both surface tension and bending rigidity with volume conservation and allows the inclusion of external forces. We present a numerical method for the simulation of the model. Equations for and on the surface of the cell are discretized on a discrete surface. The efficacy of our methodology is illustrated by computer simulations of pseudopoddriven migration and persistent migration in three space dimensions and the simulation of cell migration in the presence of obstacles.
While a major part of this work is the investigation of modelling generalizations through numerical simulation, especially with respect to cell motility in three dimensions, our intention is to present a general modelling framework and numerical method that will be a potentially useful methodology for experimentalists and theoreticians alike in future studies of cell motility.
The remainder of our discussion proceeds as follows. In §2, we present our general modelling framework and our modelling assumptions. In §3, we present a numerical method for the approximation of surface evolution laws coupled with surface partial differential equations. In §4, we report the results of numerical simulations of a model for chemotaxis. In §5, we report the results of numerical simulations of a surfacebased model for the persistent motion of cells such as fish keratocytes. In §6, we discuss the implications of our findings in the study of cell motility and possible applications of our methodologies in future studies. We provide the technical details of our modelling and the numerical methods we employ in the electronic supplementary material.
2. A surfacebased model for cell polarization and movement
We consider models for cell motility and chemotaxis that consist of a geometric evolution law for a hypersurface representing the cell boundary coupled to a spatial pattern generator on the evolving surface describing polarization of the cell. The particular form of the spatial patterning mechanism we shall investigate is a Turing pattern generator, i.e. a semilinear reaction–diffusion system (RDS). The use of Turing type systems to model biological pattern formation phenomena is widespread (see Murray [19] for a review), and recent numerical studies of Turing type systems on evolving surfaces show that while the key features of Turing mechanisms persist, such as spontaneous pattern formation and bifurcations owing to surface evolution, the geometry of the evolving surface strongly influences the patterns expressed [20,21]. We stress that our general modelling strategy and the numerical methods we employ can be generalized to other possible polarization mechanisms, such as gradientbased models or excitable network and wavebased models all of which effectively couple surface partial differential equations to a surface evolution law [3].
2.1. Geometric evolution model
The cell membrane is represented by an evolving hypersurface, with the movement of the cell determined by the interaction of various forces that act normal to the cell membrane. We consider external forces such as a protrusive force associated with the RDS species and a pressure that constrains the enclosed volume, as well as internal forces that arise from the reaction of the cells' surface to stretching and bending. We use the following force balance on the membrane, where we neglect the inertia of the membrane: 2.1where denotes the outward pointing unit normal to the surface Γ. We account for the following force contributions appearing in equation (2.1).
— A protrusive force depending on the densities of chemical species resident on the membrane (cf. equation (2.9)) is denoted by 2.2In the subsequent numerical simulations, we make the phenomenological modelling assumption that the force is proportional to the species densities and given by The sign of the component of the vector governs whether the ith species promotes protrusion (positive) or retraction (negative) of the cell membrane. For such a force may model the protrusive force generated by crosslinked filamentary actin in the vicinity of the cells' surface, while could correspond to the contraction force generated by actin bundles [22].
— Experimental studies suggest that while the cell surface area may exhibit variability during movement, the enclosed volume is relatively constant [18,23]. We take this fact into account as a hard constraint, which means that the cell is able to immediately counterbalance small volume changes on the time scale of the RDS and the boundary evolution. In the following, a corresponding Lagrange multiplier will be denoted by . It can be interpreted as a pressure difference between the interior and exterior of the cell. The corresponding force simply reads
Note that the Lagrange multiplier λ is spatially constant and therefore models a spatially uniform force such that the enclosed volume is conserved.
— We include a viscous force that opposes motion: 2.3where is a kinetic coefficient and V is such that , where is the material velocity of the cell boundary (which we assume to be normal to the cell membrane). In the twodimensional case, adhesion and deadhesion may be modelled as an effective friction, i.e. a force, of the form (2.3), proportional to the local speed [22,24]. In the threedimensional case, the situation is more complicated and we intend to address this subject in future work.
— We write for any other external normal forces acting on the cell boundary, where we have interaction with the medium in mind. As a concrete example and in order to illustrate the versatility of the proposed approach for cell motility in the complicated environments encountered in vitro, in §4.1.3, we present a simple phenomenological model for the movement of cells in the presence of obstacle particles that the cell cannot invade but that it may push out of its way.
— Resistance of the cell boundary to stretching may be incorporated by means of a surface energy of the form: 2.4where can be interpreted as a surface tension. The variation of the area functional is the mean curvature (we refer to Deckelnick et al. [25] for details of a derivation), hence the force arising from the surface energy is given by: 2.5where H is the mean curvature of Γ.
— The lipid bilayer forming the basic component of the cell membranes also resists bending. We consider the established model of Helfrich [26] for the bending energy: 2.6where is the bending rigidity. The variation of the bending energy yields the force contribution 2.7where and Δ_{Γ} denote the surface gradient and Laplace–Beltrami operator, respectively (electronic supplementary material). We refer to Willmore [27] for a derivation.
Summing up the forces with their specific choices except for the external force, we obtain the following equation for the evolution of the cell boundary: 2.8The variational formulation of the evolution law (2.8), which we use to construct a finiteelement discretization, is given in the electronic supplementary material.
2.2. Cell polarization model
We consider an RDS posed on an evolving surface : 2.9where , is the number of chemical species involved, denotes the density of the ith chemical species, is the material velocity of the surface (cf. (2.3)), 2.10is the material derivative with respect to the velocity , is a diagonal matrix of positive diffusion coefficients and is the reaction. For details of the derivation we refer, for example, to Dziuk & Elliott [10] and Barreira et al. [21]. We assume in the following that the evolving hypersurface is closed so that no boundary condition is required. For the initial condition, we write 2.11
3. Discretization
Here, we describe the numerical methods we shall employ for the approximation of models for cell motility of the form described in §2. We keep the exposition relatively nontechnical referring the interested reader to the electronic supplementary material for the technical details. We decouple the approximation of the surface evolution and the RDS by treating the RDS concentration explicitly in the surface update step. The numerical method is based on approximating the surface Γ(t) with a triangulated surface Γ_{h}(t), which stems from the method described in the seminal paper of Dziuk [28].
3.1. Definition (triangulated surface)
A triangulated surface Γ_{h} is a polygon or a polyhedron for d = 2 or 3, respectively, with linear edges for d = 2 and planar faces for d = 3, such that 3.1where consists of a finite number of closed intervals (d = 2) or a finite number of closed nondegenerate triangles (d = 3). For the simulations on surfaces (d = 3), we make use of a quadratic triangulated surface. That is, a surface that consists of curvilinear triangles each of which is the image of a reference triangle under a quadratic map as illustrated in figure 1. We will use to denote both a triangulated surface and a quadratic triangulated surface interchangeably. So as no confusion arises, we stress that for the approximation of smooth curves, we consider triangulated (polygonal) discrete surfaces and for the approximation of smooth surfaces (d = 3) we consider quadratic triangulated surfaces. For details on triangulated surfaces, quadratic triangulated surfaces and approximation results, we refer to the earlier studies [10,29–32].
The evolution law (2.8) is discretized using a surface finite element described in detail in the electronic supplementary material. The method is based on the parametric finiteelement methods for fourthorder geometric evolution equations derived in Dziuk [33] and Barrett et al. [11]. Under the proposed method, the movement of the nodes of the triangulation satisfies the evolution law in the normal direction and includes a tangential velocity (that leaves the evolution law unchanged), which gives highly desirable meshproperties in practice. Figure 2 illustrates an example of the robustness of the proposed scheme in the approximation of large deformations over schemes where nodes are moved solely in the normal direction.
To approximate the RDS posed on the evolving surface, we employ a surface finiteelement method based on the evolving surface finiteelement method proposed by Dziuk & Elliot [10], where we account for the tangential velocity induced by the surface update scheme.
We also describe, in the electronic supplementary material, a framework for the inclusion of stochastic external signals into the model. In particular for the subsequent simulations of chemotaxis, we include a signal that is modelled by independent stochastic differential equations posed in each element (interval or triangle) which we approximate using the Euler–Maruyama method.
4. Modelling pseudopoddriven chemotaxis
We investigate a model for eukaryotic chemotaxis originally proposed by Meinhardt [36]. The original model was posed at the discrete level and consisted of a three species RDS with a spatially varying local activator, a spatially varying local inhibitor, and a spatially constant global inhibitor. Meinhardt proposed that such a model could account for the polarization of chemotactic cells and the subsequent relocation and splitting of activator peaks in response to changing external signals. He did not, however, consider the mechanical aspects of the evolution of the cell membrane. Neilson et al. [7] investigated a continuous form of the model where the three species were all spatially dependent, approximating the model equations with a surface finiteelement method for the RDS approximation and a level set method for the surface update. They have conducted detailed comparisons of their simulations using a level set method with experiments [37], as well as some preliminary investigations into robust computational methods, specifically shorttime simulations using a surface finiteelement method [16]. All their modelling and simulation was conducted in two dimensions, the model we consider extends the previous work by increasing the dimension, accounting for the bending energy and modelling obstacles.
In the original model posed by Meinhardt, the spatial independence of the global inhibitor is used in the derivation of the model. Since the global inhibitor is spatially constant, its concentration can be obtained by averaging, i.e. the use of a nonlocal term (the mean value of the local activator). We, therefore, consider the following transformation of Meinhardt's model from the spatially discrete fixed surface setting to a continuous evolving surface:
4.1
Here the r_{i}, s_{i} and b_{i}'s are material parameters, γ is a scaling parameter that governs the overall timescale of the reaction rate and the term s(x,t) models both the underlying noise from the external media and a noisy chemotactic signal. The term involving the chemotactic signal feeds multiplicatively into the autocatalytic term in the activator equation. Weak signals are amplified owing to autocatalysis, thus the model provides a mechanism for gradient sensing of small chemotactic gradients. The RDS is coupled to the surface evolution law solely through the activator concentration a_{1}, which promotes protrusion of the cell membrane. Adopting the notation of §2, the contribution of the RDS species to the evolution law (cf. (2.2)) reads 4.2
In fact, the activator and inhibitor are in phase and this is one potential drawback of the Meinhardt model as it presents no obvious mechanism for coupling surface concentrations to retraction of the cell membrane. One can show that the model (4.1) is equivalent to the model considered by Neilson et al. if the diffusivity of the global inhibitor in the Neilson model is sufficiently large relatively to the diffusivity of the other two species [38].
We model the stochastic term s(x,t) in equation (4.1) as the sum of an underlying noise term owing to, say, heterogeneity in the medium and a term that models the cells' sensing of the chemotactic signal, i.e. . The underlying noise term (present in all the simulations) satisfies the following Ornstein–Uhlenbeck (meanreversion) stochastic process: 4.3where W^{t} denotes the Wiener process. We discretize in space by assuming is constant over each finite element, and the Euler Maruyama method is used to approximate the solution numerically (for details see the electronic supplementary method).
4.1. Results
For the results on curves, we took the unit circle as the initial steady state and used the reaction kinetic parameter values given in table 1. The parameter values are those Meinhardt used in his original study rescaled such that the diffusivity of the activator is 1.0, the only parameter we have tuned is the saturation of the activator concentration s_{1}, which is smaller to account for dilution in the activator concentration owing to domain growth. For all the simulations on curves, we used the same equidistributed initial mesh with 1024 d.f. (further refinement did not change the solutions qualitatively) and a timestep of 10^{−5}. The central processing unit (CPU) times of all the calculations on curves is in the order of minutes, for example, the two simulations reported in figure 3 had CPU times of just over 2000 s.
4.1.1. Random migration
Figure 3 shows the centroid trajectories of five cells migrating under two different geometric evolutions, surface tension evolution () and combined surface tensionelastic evolution ), with no chemotactic signal and different realizations of the noise term. In both cases, the cells initially polarize and then migrate, with roughly two pseudopods at the front of the cell at any given time. The cells change direction via splitting and decay of pseudopods with one pseudopod splitting and persisting while the other decays. This leads to the characteristic linear motion over short times interspersed with sharp changes in direction corresponding to a pseudopod splitting/decay event, similar to that observed by Neilson et al. [7]. Under surface tension evolution, the cells maintain a characteristic shape (two pseudopods at the front of the cell with a valley between them), while the introduction of bending rigidity generates a greater variety of cell shapes and the cell no longer has a characteristic shape with varied bananalike shapes evident (the cell shapes resemble those in the simulations presented in figure 5). From the simulations, we observe that the diluting effect of protrusion plays an important role in destabilizing activator peaks and pseudopod splitting. By this, we refer to the fact that the formation of an activator peak results in protrusion of the cell membrane, which in turn leads locally to an decrease in activator density (as protrusion may be viewed as volumetric expansion). As the local maxima corresponding to the activator peak is reduced most at the tip of the peak, where protrusion is largest, this has the effect of increasing the propensity of activator peaks and hence pseudopods to split.
We may proceed to estimate some of the parameter values using available experimental data, which is readily available for dictyostelium cells [2]. The typical radius of a cell cross section is 4 µm, which sets the length scale for the computations. The maximum actin polymerization velocity (which is related to the nondimensional parameter ) is approximately 0.1 m s^{−1}, thus the value of together with the maximum density of a^{1} in the simulations (approx. 30) sets the timescale for the simulations. Typical values of the surface tension are 10 pN µm^{−1}, assuming a cell height of 0.1 µm, this sets a kinetic scale for the simulation. The remaining physically elevant parameters may thus be estimated and are given in table 2. Note that the timescale we refer to in table 2 corresponds to one unit of computational time. The length of the simulations in figure 3 is one computational time unit or 30 min in actual time and corresponds to roughly 20 pseudopod lifetimes (each change in direction in figure 3 represents a pseudopod splitting/decay event), thus the timescale of an individual pseudopod is around 90 s.
We note that other choices of the material parameters, specifically weaker surface tension, gives cells with more elongated shapes, larger protrusions, and cell bodies that appear less rounded.
4.1.2. Migration in the presence of a chemoattractant
We now include a chemoattractant in the model. We use the stochastic receptor model proposed by Neilson et al. [7] to model the noisy chemotactic signalling. For completeness, we state the essential details. At time , they model the cells sensing of the chemotactic signal R^{t} with an Ornstein–Uhlenbeck stochastic process of the form: 4.4where W^{t} denotes the Wiener process, models the strength of the chemotactic signal, the rate of reversion to the mean μ, and the variance. The mean μ is local and prescribed by the model, while the rate of reversion to the mean and variance is local too as they are chosen such that θ = 1/(1 − μ) and σ = cμ^{1/2} (for details, see Neilson et al. [7, §6.2.2]). To compare with Neilson et al. [7, §6.2.5], we model the signal such that if a chemoattractant is present, the mean μ varies from the base signal strength (at the back of the cell) of 0.5 to 0.5 + ρ at the front, where ρ > 0 represents the signal strength. For a given signal direction , the position of the rear of the cell is such that . We have also conducted experiments where , where denotes the location of a static point source of chemoattractant and observe similar results. We discretize (in space) by assuming the mean μ is constant over each finite element, and hence θ and σ are also constant on each finite element.
Figure 4 shows the trajectories of the centroids of five cells migrating leftwards in a linear chemotactic gradient of varying strength under conserved surface tension evolution (, . The results are similar under the other geometric evolution considered () and are not illustrated. The migration of the cells with only the base signal and the signal strength set to zero is reported in figure 5a. We observe no clear directional preference similar to the migration observed in the absence of a chemoattractant. As the signal strength is increased (at a signal strength of 0.04 around 8% of the base signal), the cells start to exhibit a clear directional preference and successfully navigate up the chemotactic gradient. In table 3, we report on chemotaxis measures of 100 cells migrating under the six different signal strengths shown in figure 4. We state the average value over the 100 simulations (for each signal strength) of the following quantities, all evaluated at t = 0.5.

— The chemotactic index (CI), defined as the cosine of the angle between a line connecting the present position of the cell centroid to the starting point and a line directly up the chemotactic gradient [39].

— The persistence length (PL) of the centroid trajectory in the x and ydirections. The persistence length is taken to be the displacement in the chosen direction divided by the total length of the trajectory of the cell centroid [40].

— The squared displacement of the cell centroid from its initial position.

— The speed of the cell.
The data suggest that as the strength of the chemotactic signal is increased, the cells exhibit greater propensity for persistent migration up the chemical gradient with chemotactic indices similar to those observed experimentally in the case of Dictyostelium cells [39; 0.71–0.94]. The persistence length in the xdirection (up the chemotactic gradient) also increases with the signal strength, while the persistence length in the ydirection is reduced. We also note increasing the signal strength leads to larger displacements. The results suggest that for values of ρ ≥ 0.06, the cell is able to migrate successfully up the chemical gradient with all the reported statistics converging to similar means for further increases in the signal strength. The (physical) cell speeds are similar to those observed in migrating leukocytes [41, table 1] and Dictyostelium cells [39] (in both cases reported in µm min^{−1}).
We now investigate the ability of this model to capture the ability of a cell to respond to a changing chemotactic signal. We use the same stochastic receptor model for the chemotactic signalling but now we change the direction of the signal at various stages in the evolution. Figure 5a,b shows snapshots of the cells shaded by activator concentration under the two different geometric evolutions (, and ) in response to a changing chemotactic signal. Initially, we include only the base signal with noise, i.e. the signal strength is set to zero. At the times in the evolution at which the arrows appear in the figure, we change the direction of the signal, with signal strength ρ = 0.1, to the direction indicated by the arrows. We see that under both geometric evolutions, the cell successfully responds to the changing signal, exhibiting a clear directional preference for movement in the direction of higher chemoattractant concentration. As a final example of response to a changing signal, we consider the case where the signal direction is changed by 180°. The results of such a simulation are shown in figure 6. We observe the cell successfully responds to the change in signal direction and does so via turning gradually through 180°. This corresponds to the socalled ‘hops’ (consecutive right/right or left/left splitting of pseudopods) that are an important mechanism for the reorientation of Dictyostelium cells moving in a direction more than 90° off the chemotactic gradient [42, fig. 4]. Under this model, we have however, thusfar, not observed the formation of de novo pseudopods towards the direction of increasing chemoattractant, which are another significant mechanism for major directional corrections [42].
4.1.3. Migration in the presence of obstacles
We now include an external force in the evolution law, which arises from a model for the migration of cells through a field of obstacles. We model the obstacle particles as rigid spherical bodies. The obstacle–cell interaction is described by a phenomenological repulsive force that points in the direction normal to the cell membrane with no tangential component. Unlike the models proposed by Hecht et al. [43] and Grima [44], the obstacles and the cell both move owing to mechanical interactions.
Let denote the number of obstacles with centres and radii . For the force acting on a point on the cell boundary owing to the interaction with obstacle i we postulate 4.5where is a material coefficient and is a thickness parameter: the force is zero if the distance between the cell membrane and the obstacle boundary is bigger than . The force becomes infinite as this distance approaches zero and then dominates any other forces on the cell membrane, thus preventing intersection of the cell and the obstacle. The external force acting on the cell boundary is given by: 4.6
For the obstacle particles, we postulate a viscous motion law, too, where the reaction forces from the cell boundary and obstacle–obstacle interactions are taken into account. We postulate 4.7
Here, the are positive kinetic coefficients related to the mass of the particle, the first term on the righthand side modelling the cell–obstacle interaction is 4.8and is the force from particle j exerted on particle i for which we postulate 4.9where the are material coefficients. Note that in the absence of the cell, the initial location of the obstacles is such that the sum of the forces yields zero so that the particles do not move. Moreover, we have the following balance of forces exerted by the cell on the obstacles and forces on the cell membrane owing to the obstacles
Figures 7 and 8 show a series of snapshots of cell migration through a field of obstacles, with parameter values as in table 1 and the two previously considered geometric evolutions. Our numerical experience suggests that under the simple model of cell obstacle interactions we have employed, the increase in computational time, even with a large number of obstacles, from the case of no obstacles is negligible. We include the forcing terms in the evolution law for the cell membrane and the obstacle centres given by equations (4.6)–(4.9), with parameter values for all i, j and . We observe that the cell successfully migrates through the field of obstacles maintaining the characteristic shape as it deflects the obstacles. Our numerical experiments suggest that this behaviour is sensitive to the parameter values chosen in the repulsive potential (4.5). In particular, if we set the kinetic coefficient related to the mass of the obstacles cf. equation (4.8) to be comparable in magnitude to the kinetic coefficient related to the mass of the cell (1 by assumption), which means the obstacles inhibit more strongly the protrusion of pseudopods, then pseudopod splitting no longer occurs and the cell exhibits persistent motion in the direction of an obstacle (not reported).
4.1.4. Migration of cells in three space dimensions
We now present results for the motion of threedimensional cells in the absence of a chemoattractant. We took the unit sphere as the initial steady state and used the reaction kinetic parameter values given in table 4. We selected a timestep of 10^{−5} and used the adaptive strategy described in the electronic supplementary material with parameters N_{H} = 0.5, N_{h} = 0.75, M_{H} = 0.25 and M_{h} = 0.5. We considered an evolution law of the form (2.8) with parameters . As in table 2, we give a physical interpretation of the parameter values in table 5, assuming the radius of the spherical cell at rest is 1.17 µm. The model for the random signalling was the same as the model used in the twodimensional case. Snapshots of the cell surface shaded by activator concentration are reported in figure 9. We see qualitatively similar behaviour to the case of curves with protrusion of activator peaks leading to pseudopod formation and the rest of the cell retracting behind the pseudopods. We also observe pseudopod splitting as the cell changes direction via biased generation and retraction of existing pseudopods. The simulation is considerably more challenging than the curve case considered previously and the total CPU time of the simulation was just over 27 h.
5. Modelling the persistent motion of keratocytes
We present a model based on the general modelling framework described in §2 that seeks to capture the persistent motion of fish keratocytes. The cells deform rapidly into a temporally persistent shape and once in this shape move at a constant speed without changes in direction. Keren et al. [18] conducted an analysis of the shapes taken by moving keratocyte cells and propose a simple phenomenological model to account for the observed movement and cell shapes. Their results suggest that the steadystate shapes of the cell are broadly described by two modes and that cell shape, specifically the aspect ratio (length/width), is strongly correlated with the speed of motion. They also examined the actin distribution within the cell. Branched actin filaments promoting protrusion are concentrated at the fastmoving front and retraction promoting actin bundles are concentrated at the rear. The steady state appears stable to perturbations and if movement of the cell is disrupted, the cell rapidly regains its previous shape and speed of movement, usually moving in a new direction.
The observed behaviour of spontaneous polarization and subsequent development of a steady state stable to perturbations suggests a Turing type mechanism coupled to a surface evolution law could accurately capture the observed dynamics. Shao et al. [22] considered a membrane subject to surface tension, bending rigidity and forcing with volume conservation. The forcing strength was dependent on the concentrations of a two component RDS posed in the bulk of the cell. They present computational results, for twodimensional cells with weak volume conservation (enforced via penalization), based on a phasefield method. Ziebert et al. [6] present a model for keratocyte movement, again based on a phasefield method, where they couple the surface evolution to a vector field that seeks to describe the polarization of the actin network. Studies suggest that branched filamentary actin and actin bundles are concentrated primarily near the cell membrane near areas of protrusion and retraction, respectively, while away from the cell surface the actin is in a remodelling phase between that of branched and bundled actin [18,45,46]. This suggests a surface model where the pattern formation process occurring on the cell membrane itself may be appropriate. We propose the activatordepleted substrate model [47]: 5.1
We first present results for curves, with material and RDS parameters given in table 6. We considered an initially circular cell with radius 1 centred at the origin. The initial condition for the RDS was taken as the linearly stable steady state with a symmetry breaking perturbation of the form added to the initial condition of the a_{2} species. The specific form of the initial condition leads to cells that migrate only along the xaxis (we verified that the choice of other initial conditions only changed the direction of the movement). The hypothesis of Keren et al. [18] is that variability in the actin dynamics is the major factor governing the observed variations in shape and speed. To investigate this hypothesis, we propose that the a_{1} species in the RDS (5.1) corresponds to the density of retraction promoting actin bundles, while the a_{2} species corresponds to the density of protrusion promoting actin filaments, which is similar to the model considered in Shao et al. [22]. We can model variable actin dynamics by changing the constant k_{2} which can be interpreted as the growth rate of actin filaments. Increasing k_{2} leads to higher concentrations of a_{2} relative to a_{1}, and thus should lead to faster moving cells with stronger forcing at the front.
In all the simulations on curves, we used an initially equidistributed mesh with 1024 d.f. and a fixed timestep of 10^{−3}. The CPU times were on the order of seconds with a typical simulation taking approximately 200 s. Figure 10a,b shows, for different values of k_{2}, the initial position of the cells at time 0 and the cell positions and surface RDS concentrations at time 5 (by which time all the cells have reached a steady state with constant speed and time independent RDS concentrations). We see faster cell speeds and larger aspect ratios for increased values of k_{2}, similar to the models where the RDS is posed in the bulk [22]. The shapes of the cells at steady state also resemble those observed experimentally [18]. In particular, we see the rounded ‘D’ shape in the right most cell corresponding to k_{2} = 0.6, and the much more elongated ‘canoe’ shape in the left most cell corresponding to k_{2} = 1.8. We report on the aspect ratio , as considered in Ziebert et al. [6], where for i = 1,2, is as follows (the 's are the eigenvalues of the diagonal 2 × 2 variance matrix of the cells centroid): where is the ith coordinate of the cells centroid. We also report on the deviation from reflection symmetry of the migrating cells as considered in Ziebert et al. [6]. This is measured by the following quantities (the nonzero components of the skewness tensor of the cells centroid scaled by a constant factor):
Figure 11 shows plots of the speed of the cell centroids and the aspect ratio of the cells both against time. We clearly see the positive relationship between aspect ratio and cell speed evident in the experimental studies. In physical units, the range of the speed at steady state of the cells shown in figure 10 is 0.178–0.445 µm s^{−1}}. Both the speed and aspect ratios are similar to those observed in the experimental results reported in Keren et al. [18, fig. 4b]. Figure 12 shows values of the asymmetry measures against time. We see that the cells travelling at a slower speed exhibit larger deviations from reflection symmetry. This is in contrast to the results obtained under the model considered in Zibert et al. [6, fig. 4], where does increase as the speed of the cells decreases, but is positively correlated with speed. We note that for the first two simulations and , is negatively correlated with speed. As an ellipse satisfies , it is not clear that the faster moving cells with larger aspect ratios, but more elliptical profiles, should exhibit larger deviations from reflection symmetry. We also observe that after a brief initial stage in which all the reported values oscillate, the values converge to a steady state as the cells travel in a persistent fashion with a fixed shape and a constant speed.
We also report on simulations of threedimensional keratocyte motion. We took the unit sphere as the initial cell shape, the same initial conditions for the RDS concentrations as in the curve case, a fixed timestep of 10^{−4} and remaining parameter values for both the surface evolution and adaptive strategy as given in table 8. The CPU times were on the order of minutes with a typical simulation taking approximately 2000 s. Proceeding as in §4, we give physical interpretations of the parameter values for curves and surfaces in tables 7 and 9, respectively. Figure 13a,b shows a similar experiment to the one carried out for curves now on surfaces, specifically we report for different values of k_{2}, the initial position of the cells at time 0 and the cell positions and surface RDS concentrations at time 5 (by which time all the cells have reached a steady state with constant speed and timeindependent RDS concentrations). The gross behaviour is the same as the curve case, in that as the parameter k_{2} is increased from 0.6 to 1, the cells move faster at steady state and appear more elongated. Figure 14 shows plots of the speed of the cell centroids and the surface area of the cells both against time. We plot the surface area as it is proportional to the two, roughly equal, aspect ratios in the (x, y)and (x, z)directions. We observe the same positive relationship as in the curve case, with both surface area and speed converging to steady states. We have also verified that the aspect ratios converge to steady states with the aspect ratio in the (y, z)direction approaching 1. Note for larger values of k_{2}, the cells developed a selfintersection which is inadmissible under our modelling as it would correspond to a change in topology in the physical setting. Scenarios where one wishes to consider a topological change, or respectively, methods that avoid topological change, are a subject of our current research.
6. Conclusion
In this work, we have presented a computational framework for the modelling of cell motility. We propose a simple and consistent means of coupling cell movement with gradient sensing, polarization using surface partial differential equations (PDEs), and external forces. Our methods can be generalized to the modelling of more complex phenomena such as adhesion and crawling on a substrate or cell–cell interactions and we illustrate one such generalization with a concrete example of migration in the presence of obstacles.
A contribution of our study is the description of a numerical method for the simulation of cell movement that can account for the large deformations that arise in simulations of cell motility and that can be applied to the study of threedimensional cell migration. The model equations consist of PDEs for and on surfaces and the numerical method seeks to approximate these equations on a discrete surface. Thus both the continuous and discrete problems are posed in one dimension less than the underlying spatial dimension, which in the case of the discrete problem typically means fewer degrees of freedom are needed than would be the case for embedded methods [10,13–16]. On curves, our experience is that the method maintains a mesh suitable for computation without the use of remeshing or adaptive mesh refinement. For the simulations of cell migration in three dimensions, we occasionally observed deterioration in the mesh quality, even with the redistribution of the vertices implicit in our numerical scheme, necessitating spatial adaptivity. An area of our ongoing research is the investigation of numerical methods robust to large deformations in the cell surface.
We consider a pseudopodcentred model for chemotaxis similar in form to that considered in Neilson et al. [37]. Unlike compass models [39], which are reasonable for cells with flexible polarity where a large gradient may induce pseudopods at any position on the cell membrane, pseudopodcentred models [17] are suitable for strongly polarized cells, where pseudopods are generated preferentially at the front with directional bias, owing to a chemotactic gradient, restricted primarily to small changes in direction [42]. The major contributions and imports of our study are the inclusion of bending rigidity, the inclusion of external forces, the observation that the gross behaviour, of pseudopod splitting, observed in Neilson et al. [37] for twodimensional cells persists in threedimensional simulations, and that the model remains qualitatively unchanged when one considers a twocomponent RDS, with a spatially constant global inhibitor, rather than a three component RDS with a biologically implausible nonlocal term. Our computational method based on surface finite elements extends the method in Neilson et al. [16] and is an alternative to the level set method considered in Neilson et al. [7,37]. The simulations illustrate that the model is capable of reproducing aspects of pseudopoddriven cell migration, described in Insall [17], in both two and three space dimensions. We report on many widely used chemotaxis measures and observe values similar to experimental observations. We also note that the simulations exhibit a dilution effect at the tip of a pseudopod, where the local maxima corresponding to an activator peak is reduced. This suggests experimental investigation of the relative importance of mechanical effects of membrane protrusion on the distribution of cellresident proteins is warranted.
We also investigated a model for the motion of fish keratocytes. The model appears to reproduce some experimental observations of the shapes of motile keratocyte cells and the experimental observation of the correlation between cell shape and speed [18]. The computational model in Keren et al. [18] reproduces the velocity–aspect ratio relationship. However, unlike our model, both polarization and cell shapes are not explicitly modelled, with a parabolic actin profile at the leading edge assumed and the shape of the cell rear neglected. Studies [6,22] propose models where polarization is modelled by equations in the bulk of the cell which are coupled to an evolution law for the cell surface. The import of our study is to show that a surface RDS coupled to a surface evolution law gives qualitatively similar results. A further contribution is the use of surface finite elements rather than the phasefield method considered in Ziebert et al. [6] and Shao et al. [22]. This allows simulation of threedimensional keratocyte migration, in which studies [6,22] both note is computationally expensive with the phasefield methodology. We do observe minor differences from Ziebert et al. [6], for example, in the measures of deviation from reflection symmetry.
Our numerical experience suggests that some aspects of cell migration and chemotaxis can be captured by the Schnakenberg RDS, equation (5.1). This RDS is considerably simpler from a mathematical analysis viewpoint than say the Meinhardt model, equation (4.1). One can show that the model is well posed on evolving (planar) domains [48], which is an open question even on fixed domains for the Meinhardt model. As the two components are out of phase, the model lends itself naturally to the case of a species that promotes protrusion (e.g. actin) and another that promotes retraction (e.g. myosin).
In this work, it is our intention to present a framework for future modelling rather than suggest any definitive models for cell migration. We hope that future studies will employ the framework we have set out to refine existing models for cell motility and make predictions based on numerical simulations that can be used to direct and inform experimental studies.
Acknowledgements
This research has been supported by the UK Engineering and Physical Sciences Research Council (EPSRC), grant no. EP/G010404.
 Received April 2, 2012.
 Accepted May 8, 2012.
 This journal is © 2012 The Royal Society