Directed cell migration poses a rich set of theoretical challenges. Broadly, these are concerned with (1) how cells sense external signal gradients and adapt; (2) how actin polymerisation is localised to drive the leading cell edge and Myosin-II molecular motors retract the cell rear; and (3) how the combined action of cellular forces and cell adhesion results in cell shape changes and net migration. Reaction–diffusion models for biological pattern formation going back to Turing have long been used to explain generic principles of gradient sensing and cell polarisation in simple, static geometries like a circle. In this minireview, we focus on recent research which aims at coupling the biochemistry with cellular mechanics and modelling cell shape changes. In particular, we want to contrast two principal modelling approaches: (1) interface tracking where the cell membrane, interfacing cell interior and exterior, is explicitly represented by a set of moving points in 2D or 3D space and (2) interface capturing. In interface capturing, the membrane is implicitly modelled analogously to a level line in a hilly landscape whose topology changes according to forces acting on the membrane. With the increased availability of high-quality 3D microscopy data of complex cell shapes, such methods will become increasingly important in data-driven, image-based modelling to better understand the mechanochemistry underpinning cell motion.

## Introduction

Directed migration of cells in response to extracellular signals is of fundamental importance during development [1], in the immune response [2], and in cancer [3]. It has been studied in great detail in experiments, but theoretical models helped tremendously to understand the regulatory principles of cell motility. Figure 1 gives an example of Dictyostelium amoeba responding to the inversion of a signal gradient of mechanical shear flow by localising activators of actin polymerisation to the cell front. There we also sketch one widely used framework for modelling gradient sensing, namely that of Turing-type reaction–diffusion models, and the principal steps of linking gradient sensing to cell polarisation, force generation, and net movement. Fully integrated models that include cell shape changes and actual movement are still scarce, but thanks to increased computational power and new mathematical methods, impressive progress has been made recently. Examples are models by Niculescu et al. [4] who showed that a smooth transition between amoeboid cells moving with multiple, round protrusions and keratocyte-like cells moving with a single flat lamella can be achieved by changing one parameter in a Cellular Potts model. In the Cellular Potts model, cell shape is modelled as a connected set of points on a grid. Cells can locally expand and shrink by incorporating new grid points or loosing points, depending on what is energetically more favourable. Camley et al. [5] investigated how cell turning can arise from a coupling between the reaction–diffusion mechanism and cell shape, and Alonso et al. [6] considered cell-to-cell variability in models of amoeboid motility. The latter two use what is called an interface capturing approach to model the evolution of the cell membrane, which has become increasingly popular, because of its mathematical rigor and the ability to efficiently model complex cell shapes. In the current minireview, we want to highlight recent developments in mathematical models for membrane evolution, with a specific focus on describing methods for interface tracking and capturing. More extended reviews of mathematical models of cell motility can be found in [7–10].

### Modelling the mechanochemistry of cell migration.

Figure 1.
Modelling the mechanochemistry of cell migration.

(A) A Dictyostelium cell labelled with a marker for F-actin (green) undergoes repolarization in response to inverting the direction of hydrodynamic shear flow ($σ=2.1Pa$) indicated by white arrows. $T=0s$ is the time of flow reversal. Bar: $10μm$. (B) Time-space plot of the actin fluorescence averaged over 18 responses from 14 cells as in (A) (red: high intensity, blue: low). Numbers along the cell contour in (A) indicate the position along the circumference (y-axis) in (B). (C) Spatial profiles of actin along the cell contour at selected time points. Green curves show model fits of a modified Meinhardt model [46] for cell reorientation. (A,B,C adapted from [47]). (D) A large number of models for gradient sensing, such as the Meinhardt model, are based on Turing-type reaction–diffusion models, with a caricature of a generalised 2-variable model sketched here. An external signal $S$ stimulates the production of an activator variable $A$, through a reaction $fA$. $fA$ is required to be non-linear, which can be met by assuming, for example, that $A$ stimulates its own production in an autocatalytic way. The production of $A$ will be counteracted by an inhibitor $I$ entering $fA$, either in a linear or non-linear way. Here, $A$ will stimulate the production of its own inhibitor through $fI$. This system of equations which describe the rate of change of $A$ and $I$ is extended by assuming diffusion of $A$ and $I$. In our example in (C) the actin concentration has been considered a read-out for the activator variable $A$. Diffusion has been considered to take place along the cell contour. Other models allow diffusion throughout the entire cell interior. In Turing-type models a shallow signal gradient can result in a sharp activator profile at the cell front. To obtain such a pattern the diffusion of the inhibitor is required to be larger than that of the activator. Many variations on this theme of models exist, for example to allow for perfect adaptation, stochasticity, limited numbers of molecules, molecules that can shuttle between active and inactive states, etc. (E) Gradient sensing as described in (D) has been shown to be independent of the actin system. However, the activator/inhibitor system can stimulate localised polymerization of actin at the cell front and Myosin-II recruitment to the cell rear, thus establishing cell polarity which defines the axis of cell migration (left). Middle: Diagram illustrating the combined action of actin, Myosin-II and differential substrate adhesion. Right: The net result is protrusive activity of the cell front and retraction of the cell rear. Details of internal membrane forces will be presented in the following figures.

Figure 1.
Modelling the mechanochemistry of cell migration.

(A) A Dictyostelium cell labelled with a marker for F-actin (green) undergoes repolarization in response to inverting the direction of hydrodynamic shear flow ($σ=2.1Pa$) indicated by white arrows. $T=0s$ is the time of flow reversal. Bar: $10μm$. (B) Time-space plot of the actin fluorescence averaged over 18 responses from 14 cells as in (A) (red: high intensity, blue: low). Numbers along the cell contour in (A) indicate the position along the circumference (y-axis) in (B). (C) Spatial profiles of actin along the cell contour at selected time points. Green curves show model fits of a modified Meinhardt model [46] for cell reorientation. (A,B,C adapted from [47]). (D) A large number of models for gradient sensing, such as the Meinhardt model, are based on Turing-type reaction–diffusion models, with a caricature of a generalised 2-variable model sketched here. An external signal $S$ stimulates the production of an activator variable $A$, through a reaction $fA$. $fA$ is required to be non-linear, which can be met by assuming, for example, that $A$ stimulates its own production in an autocatalytic way. The production of $A$ will be counteracted by an inhibitor $I$ entering $fA$, either in a linear or non-linear way. Here, $A$ will stimulate the production of its own inhibitor through $fI$. This system of equations which describe the rate of change of $A$ and $I$ is extended by assuming diffusion of $A$ and $I$. In our example in (C) the actin concentration has been considered a read-out for the activator variable $A$. Diffusion has been considered to take place along the cell contour. Other models allow diffusion throughout the entire cell interior. In Turing-type models a shallow signal gradient can result in a sharp activator profile at the cell front. To obtain such a pattern the diffusion of the inhibitor is required to be larger than that of the activator. Many variations on this theme of models exist, for example to allow for perfect adaptation, stochasticity, limited numbers of molecules, molecules that can shuttle between active and inactive states, etc. (E) Gradient sensing as described in (D) has been shown to be independent of the actin system. However, the activator/inhibitor system can stimulate localised polymerization of actin at the cell front and Myosin-II recruitment to the cell rear, thus establishing cell polarity which defines the axis of cell migration (left). Middle: Diagram illustrating the combined action of actin, Myosin-II and differential substrate adhesion. Right: The net result is protrusive activity of the cell front and retraction of the cell rear. Details of internal membrane forces will be presented in the following figures.

## Modelling paradigms for membrane movements

The cell membrane has a unique role in cell migration. It is the interface upon which intracellular forces act in the form of pushing, through the growth of actin filaments against the membrane, and pulling, through Myosin-II mediated contraction of the cortical actin network attached to the membrane. Furthermore, it will exchange forces with the environment, for example, through sites of attachment to the extracellular matrix, or experience drag forces by moving through the extracellular fluid. Membrane tension and resistance to bending are additional (internal) forces at play, determined by the physical properties of the membrane.

Treating the cell as a continuous domain with the membrane as domain boundary, a set of partial differential equations can be employed to describe the movement of the membrane under the influence of all external and internal forces. The exact geometry of the cell boundary at a certain point in time is then found by solving the system of equations, and mathematically we term this a free boundary problem. Closely related problems describing either evolving free surfaces, or interfaces of two fluids or a fluid with a solid, appear in many applications of fluid dynamics. In this community, two main paradigms have evolved that underpin modelling and computational approaches: interface tracking approaches and interface capturing approaches.

## Interface tracking

Interfaces such as the cell boundary may be modelled by evolving surfaces (in 3D) or curves (in 2D) if the scales allow to neglect their thickness, which for a lipid bilayer of ∼10 nm is true when compared with typical cell sizes in the range of tens of micrometres. Interface tracking approaches are based on explicitly describing the position and movement of the interface, typically by inferring it from knowledge about its local velocity $v→$ (Figure 2). Numerical methods based on this paradigm then compute trajectories of points with co-ordinates $x→(t)$ on the interface at timepoints $t$ by integrating their velocity, i.e. solving $∂tx→(t)=v→(x→(t),t)$. Note that points on the interface may move in the normal direction and in the tangential direction, within the plane of the interface. Only movement of points in the normal direction results in a change of the shape itself. The tangential motion describes transport along the membrane surface but has not impact on the geometry.

### Modelling cell shape changes using interface tracking.

Figure 2.
Modelling cell shape changes using interface tracking.

(A) Snapshots of a surface finite element simulation from [12] of a model as depicted in (C). The reaction–diffusion system is of Turing-type (see Figure 1D). Colours: concentration of an activator variable A (red: high, blue: low). The deformation along the surface normal is directly coupled to the activator concentration. (B) Surface mesh obtained from 3D image data [48]. Numerical methods such as the finite element method can be set up on such triangulations to solve models as in (C). (C) Model of the form as used for the simulation in (A); All forces $f→$ acting on each point of the surface sum up to zero. The force balance yields an equation for the velocity $v→$ in the direction of the surface normal $n→$. It is coupled with a reaction–diffusion system from [49] on the evolving membrane. The chemical species in the reaction–diffusion system can be expressed in vector form, for example $c=(A,I)$ would include the activator and inhibitor of a Turing-type system in a single vector. The individual forces are as follows: $f→drag=−μv→$ accounts for friction the surface experiences when moving through the extracellular medium ($μ>0$ is a friction parameter), $f→pressure=pn→$: hydrostatic pressure which prevents the surface from collapsing, $f→coupling=g(c)n→$: describes how c deforms the surface by an activator dependent force $g(c)$ pushing in direction of the surface normal $n→$. $f→tension=σκ→$: surface tension, resisting stretching of the model membrane with $κ→$ being the curvature and $σ$ a material tension parameter. Additional forces may be due to resistance to bending etc. The terms in the reaction–diffusion system can be interpreted as follows: $(∂t+v→⋅∇)c=∂tv→c$: material derivative of $c$, describing how the chemical species contained in $c$ change at a given node of the surface mesh in time; $c∇Γ⋅v→$: this term, with $∇Γ$ being the gradient of the surface $Γ$, models changes to the densities in $c$ due to stretching or compressing of the membrane during the movement; $DΔΓc$ : Diffusion of chemical species $c$ within the surface, with $ΔΓ$ being the Laplace–Beltrami operator and $D$ the diffusion constant; $r(c)$ reaction term describing production and decay of $c$.

Figure 2.
Modelling cell shape changes using interface tracking.

(A) Snapshots of a surface finite element simulation from [12] of a model as depicted in (C). The reaction–diffusion system is of Turing-type (see Figure 1D). Colours: concentration of an activator variable A (red: high, blue: low). The deformation along the surface normal is directly coupled to the activator concentration. (B) Surface mesh obtained from 3D image data [48]. Numerical methods such as the finite element method can be set up on such triangulations to solve models as in (C). (C) Model of the form as used for the simulation in (A); All forces $f→$ acting on each point of the surface sum up to zero. The force balance yields an equation for the velocity $v→$ in the direction of the surface normal $n→$. It is coupled with a reaction–diffusion system from [49] on the evolving membrane. The chemical species in the reaction–diffusion system can be expressed in vector form, for example $c=(A,I)$ would include the activator and inhibitor of a Turing-type system in a single vector. The individual forces are as follows: $f→drag=−μv→$ accounts for friction the surface experiences when moving through the extracellular medium ($μ>0$ is a friction parameter), $f→pressure=pn→$: hydrostatic pressure which prevents the surface from collapsing, $f→coupling=g(c)n→$: describes how c deforms the surface by an activator dependent force $g(c)$ pushing in direction of the surface normal $n→$. $f→tension=σκ→$: surface tension, resisting stretching of the model membrane with $κ→$ being the curvature and $σ$ a material tension parameter. Additional forces may be due to resistance to bending etc. The terms in the reaction–diffusion system can be interpreted as follows: $(∂t+v→⋅∇)c=∂tv→c$: material derivative of $c$, describing how the chemical species contained in $c$ change at a given node of the surface mesh in time; $c∇Γ⋅v→$: this term, with $∇Γ$ being the gradient of the surface $Γ$, models changes to the densities in $c$ due to stretching or compressing of the membrane during the movement; $DΔΓc$ : Diffusion of chemical species $c$ within the surface, with $ΔΓ$ being the Laplace–Beltrami operator and $D$ the diffusion constant; $r(c)$ reaction term describing production and decay of $c$.

A common way of finding the velocity is to formulate a force balance for the cell boundary. This modelling ansatz is illustrated in Figure 2C. Forces arise from the viscoelastic properties of the cell membrane, such as resistance to stretching (tension) and bending. A drag force proportional to the membrane velocity may be postulated due to the fact that the ambient (viscous) fluid has to move during protrusion or retraction. This approach has been used in [11,12] for modelling and simulations of cells moving on substrates in 2D, and in [13,14] for the nucleation and evolution of blebs, which are rapid hemispherical membrane protrusions driven by intracellular pressure (also in 2D). The portability to surfaces in 3D is illustrated in [12], a figure of which is reproduced in Figure 2A. Further contributions to the cell membrane's force balance may arise from other cytoskeletal components such as microtubules, which can be modelled as point forces. In [15], a discrete model for microtubules with a simplified description of the membrane is used for computing membrane tethers in 2D. Such tethers have also been computed in 3D with a continuum model for vesicles formed by biomembranes and prescribed point forces [16].

One approach to describe how actin polymerisation pushes the cell membrane outwards is to make the forces acting in the normal direction of a model membrane directly dependent on the concentration of the activator variable in Turing-type models for front activation. A generalised formulation of such Turing models is outlined in Figure 1D [17,18]. In [17], pulling forces resulting in rear retraction (in real cells this is linked to high Myosin-II activity) were made dependent on the curvature of the cell contour, assuming that regions which show no protrusive activity tend to be round and retract.

On a moving interface with associated velocity field $v→$, it is convenient to take a Lagrangian view to derive balance equations for interface resident species (molecules), for example, receptors embedded in the membrane or activators of actin polymerisation. The time derivative $∂t$ is then replaced by the material derivative$∂tv→$ which, in addition to the time derivative, contains a term to account for transport with the velocity field of the interface. For any interfacial field $c$ of, for example, concentrations of molecules, the material derivative is its instantaneous change of its value with respect to time in a material point $x→(t)$, i.e. $∂tv→c(x→(t),t)=ddtc(x→(t),t)$. A typical system of reaction–diffusion equations for a vector of species densities $c$ on the evolving cell membrane $Γ$ reads as in Figure 2C, with a detailed description of individual terms given in the legend.

The essence of numerical methods for interface tracking is that the position of the cell boundary is represented by marker points that are advanced by a time-stepping scheme according to the forces acting on them. A key difficulty is finding a suitable discretisation of the moving cell boundary $Γ$, usually using a triangulation as depicted in Figure 2B. In particular, the resolution, both of the geometry and for any fields on the cell boundary, must be sufficiently high in regions of high curvature of the boundary, or steep concentration gradients of a chemical species. Parametric approaches attempt to reformulate the problem on a reference domain. For example, 2D models can be conveniently based on closed curves for the cell boundary as used in [13,14], and for 3D problems the corresponding reference domain can be is a sphere as used for modelling blebbing by Woolley [19]. Note that the definition of a reference domain that is topologically equivalent to the cell surface is only needed for obtaining numerical solutions. The actual cell surface can be arbitrarily complex (see Figure 2B).

When considering the interaction of the cell membrane with the extracellular fluid, the movement of the cell boundary will transfer momentum onto the extracellular fluid creating a flow, which can be determined from the momentum equation. Assuming viscous, Newtonian fluids for which Stokes flow is a suitable model the equation of fluid motion can be written in the form
$0=∇⋅Σ+δfmemb,$
1
where $Σ$ is the stress tensor that incorporates viscous friction between the extracellular fluid and the cell surface, as well as the pressure of the fluid resting on the surface, $fmemb$ are all the force densities located at the membrane, and $δ$ is the surface Dirac delta function (which lets each point of the surface mesh used to model the cell boundary behave like a point force). Equation (1) underpins immersed boundary methods for fluid–structure interaction problems, and we refer to [20–23] for the specific case of cell boundaries immersed in surrounding fluid.

Moving mesh methods that monitor and attempt to keep a good mesh quality have been developed for 2D domains [24]. For 3D modelling where the reference domain is a sphere, numerical schemes are less straightforward to implement. Schemes where the reference domain is explicitly used are termed extrinsic. In intrinsic schemes, there might be a notion of a reference domain, but all computations are performed on the physical domain. For instance, the scheme in [12] uses finite elements on a triangulated surface to compute the forces and to advance the vertices. Figure 2B gives an impression of such a triangulated surface. Combining this with a system of reaction–diffusion equations from [17], 3D simulations as in Figure 2A have been performed.

Let us consider now the case of bulk-surface coupled problems, which means that model equations have to be solved inside or outside of the cell and are coupled to models on the cell membrane. The development of computational tools for such problems is very challenging. The discretisation has to enable transfer model data from the surface to the equations governing the bulk fields and vice versa. With regards to mesh-based numerical methods, there is a distinction between fitted and unfitted methods. A mesh for the interface can be obtained by taking the boundary of a bulk mesh, so the two meshes perfectly fit and exchange of data is relatively straightforward. In [18], a reaction–diffusion model for cell polarisation with equations inside and on the boundary of a (stationary) cell is solved based on triangulating the 3D cell domain and using the triangular boundary faces for a triangulation of the 2D membrane surface. Similarly, in [24], the outside of a moving 2D cell domain is triangulated so that the moving boundary edges approximate the cell boundary.

Such fitted methods involving moving meshes are not easy to implement though some software packages already provide functionality to change the mesh geometry and even extract submeshes (DUNE 2016). In turn, the idea behind unfitted approaches is to us a fixed, standard bulk mesh that is not adjusted to the geometry of the problem, and the interface mesh then moves through the bulk mesh. Here, the challenge is to determine intersections of elements of the interface domain with elements in the bulk domains to compute the contribution of interface terms to the bulk equations, or contributions of bulk terms to the interface equations. Such functionality usually does not come as a standard in software packages. But approaches have been developed and implemented for a variety of problems including multi-phase flow [25] and transferred to cell motility. We refer to [26] for simulations of vesicles formed by biomembranes, which feature mechanical properties similar to cells, and to [22] for simulations of cell motility.

## Interface capturing

In interface tracking, surfaces are explicitly defined, that is, co-ordinates of points on the surface are known. In contrary, in interface capturing, the interface representing the cell boundary is implicitly given. To explain this concept, consider a closed-smooth surface enclosing a domain which here is the cell volume. One can then define a signed distance function close to the surface, i.e. a function that is positive inside the domain, negative outside, that values 0 on the surface, and its absolute value is the distance to the surface. Inversely, if the distance function is known then the surface can be recovered.

This idea can be generalised beyond distance functions. The level sets of any smooth scalar bulk function $ϕ$ are surfaces in general. An evolving surface is obtained if the bulk function changes in time. The velocity of a level set in the direction of $∇ϕ$, which is normal to the level set, is given by $−∂tϕ/|∇ϕ|$. If now a velocity field $v→$, by which material points on the membrane move, is known (or can be computed) then we obtain the identity $∂tϕ+v→⋅∇ϕ=0$.

Level set methods (see [27,28] for overviews) are based on numerically solving this transport equation. Regarding cell motility, the membrane forces that determine the membrane's velocity $v→$ have to be expressed in terms of the level function $ϕ$ in a suitable way. The normalised gradient $∇ϕ/|∇ϕ|$ is the unit normal of the level set. Also, other geometric fields of the surface such as the mean curvature can be expressed in terms of $ϕ$. In particular, the force balance in Figure 2C can be written as a partial differential equation for $ϕ$. A common way to treat the Dirac delta function in (1) is by smoothing it so that it is defined in a thin layer around the interface [29]. Regarding applications to cell motility we refer to [30,31] where simulations in two spatial dimensions have been performed on uniform meshes. In [32], 2D simulations of the level set equation have been performed with additional equations to model adhesion. Neilson et al. [17] use a level set method for the membrane of a migrating cell in 2D, too. In the spirit of a hybrid approach, they then construct a surface mesh from the zero level set of the discrete solution in each time step and then use a surface finite element method to solve a reaction–diffusion system on the thus computed evolving membrane.

A great advantage in comparison with interface tracking approaches is that the level set equation $∂tϕ+v→⋅∇ϕ=0$ is a transport problem for which there is a variety of numerical methods such as finite difference methods, finite volume methods, discontinuous Galerkin methods, or (pseudo-)spectral (Fourier) methods. As the numerical method and the underlying mesh are independent of the geometry of the problem but rather aligned with the view of an external, fixed observer such approaches are termed Eulerian. They are in contrast with Lagrangian approaches discussed before where moving marker points indicate the position of the interface. In Lagrangian approaches, changes to the topology and, consequently, to the mesh topology are challenging to detect and implement. In turn, they are naturally captured in Eulerian, implicit approaches. Figure 3D shows a simulation of phase separation on a biomembrane forming a vesicle where the topology of the phases (red versus blue) changes. The black region indicates the position of the phase interface, which is captured using a phase field method (explained below). Note that the shape of the vesicle changes concomitantly with phase separation, because here it is assumed that each of the phases has a preferred curvature. A typical equilibrium shape with different parameters is depicted in Figure 3E. The evolution of the shape has been solved using an interface tracking method discussed in Figure 2.

### Interface capture approaches for modelling cell shape.

Figure 3.
Interface capture approaches for modelling cell shape.

(A) Left: profile of a phase field function across an interface. Middle: diffuse interface. Right: smoothed surface Dirac delta function $δε$. (B) Force balance and reaction–diffusion system from Figure 2D, smoothed by the phase field method. $ε$ is a small parameter, $φ$ is the phase field variable, $w′$ is the derivative of a double well potential with minima that correspond to the interior and exterior of the cell (±1 in (A), $cε$ are chemical species which are smoothed across the interface. The other parameters ($μ$, $σ$, $p$, $D$) and functions ($g$ and $r$) are explained in the legend to Figure 2. (C) Finite element mesh for solving phase field models, which consist of bulk partial differential equations. The solutions feature thin layers. Adaptive mesh refinement can improve computational efficiency. (D) Interface capturing and tracking methods can be combined. Here, a surface finite element method is used for tracking a deforming surface with phase separation (phases are red, $φ=1$, and blue, $φ=−1$) for which a phase field equation is solved. The topology of the phases changes during the evolution due to differing material properties, which generally can implemented with ease in capturing approaches. First published in [50] by Global Science Press. (E) Simulated equilibrium shape, using the same method but different parameters as in (D). First published in [50] by Global Science Press.

Figure 3.
Interface capture approaches for modelling cell shape.

(A) Left: profile of a phase field function across an interface. Middle: diffuse interface. Right: smoothed surface Dirac delta function $δε$. (B) Force balance and reaction–diffusion system from Figure 2D, smoothed by the phase field method. $ε$ is a small parameter, $φ$ is the phase field variable, $w′$ is the derivative of a double well potential with minima that correspond to the interior and exterior of the cell (±1 in (A), $cε$ are chemical species which are smoothed across the interface. The other parameters ($μ$, $σ$, $p$, $D$) and functions ($g$ and $r$) are explained in the legend to Figure 2. (C) Finite element mesh for solving phase field models, which consist of bulk partial differential equations. The solutions feature thin layers. Adaptive mesh refinement can improve computational efficiency. (D) Interface capturing and tracking methods can be combined. Here, a surface finite element method is used for tracking a deforming surface with phase separation (phases are red, $φ=1$, and blue, $φ=−1$) for which a phase field equation is solved. The topology of the phases changes during the evolution due to differing material properties, which generally can implemented with ease in capturing approaches. First published in [50] by Global Science Press. (E) Simulated equilibrium shape, using the same method but different parameters as in (D). First published in [50] by Global Science Press.

Another interface capturing approach is the phase field method. Interfaces are represented by thin layers, giving rise to the notion of diffuse interface approaches. An order parameter, or phase field variable, denoted by $φ$ here, is introduced that is constant within the cell and outside of the cell away from the cell boundary, and smoothly but quickly changes its value across the interfacial layer representing the cell boundary, see Figure 3A for an illustration. A small parameter $ε$ determines its thickness. Phase field models have been derived for systems where the interface carries some sort of energy, for example, a potential energy resulting from stretching that is proportional to the area of the membrane. This energy can be written in terms of $φ$, and the corresponding force approximates $σκ→$, which is the tension force introduced in Figure 2C. Similarly, the velocity of the interface and the other forces can be expressed in terms of $φ$. The force balance as in Figure 2C thus can be approximated by the bulk partial differential equation for $φ$ in Figure 3B for which again standard numerical approaches can be used.

Any equations on the cell membrane for the biochemistry also have to be reformulated to account for their representation by a thin layer. Theory and techniques have been developed to render this process relatively straightforward. A reaction–diffusion system as in Figure 2C can be written in the following form, which again involves the Dirac delta function that was used for the membrane forces in (1):
$∂t(δc)+∇⋅(δcv→)=D∇⋅(δ∇Γc)+δr(c)$
2
By smoothing the delta function, this system can be incorporated into a level set or a phase field model. See Figure 3A for a sketch of a profile of such a smoothed $δε$ across the interfacial layer and Figure 3B for the smoothed system. With regards to numerical approaches, the reaction–diffusion system now lives in the bulk domain and, as for the level set equation or the phase field equation, one can attempt to treat it with standard numerical methods. One difficulty, however, consists in the degeneracy of $δε$. Without special consideration, it leads to unstable schemes. Techniques to address this problem are presented in [33,34] for a generic surface partial differential equation, and in [35] in the context of cell motility.

Similar ideas can be used to account for partial differential equations within or outside of the cell's domain. For instance, see Ref. [36] for a phase field model that accounts for actin filament and myosin concentrations in the cell. Adhesion to a substrate can also be incorporated, which is the focus of the recent work in [37] in 2D. Simulations in 3D for a cell moving in a fibrous environment have been performed by Moure and Gomez [35,38].

In relation to level set methods, one could consider the zero level set of $φ$ as a representation of a sharp interface. But as an important difference, the derivation of phase field models is based on energy or entropy functionals, which can be further exploited for deriving robust and stable numerical schemes. For instance, numerical level set methods occasionally require re-initialisation of the level set function for accurate simulations [27,28], which is not the case for numerical phase field methods. We refer to [39] for a recent review of computational methods for phase field equations, where also the relation to level set approaches is discussed.

To obtain accurate results in numerical simulations, the thin layer around the interface has to be resolved by the numerical method. Mesh-based methods then become quickly very expensive if the mesh is refined to the required level everywhere on the computational domain. One way to make simulations cheaper is using finite element methods with adaptive refinement as discussed in [39]. Figure 3C gives an impression of a mesh that is refined only in the interfacial layer for an ellipsoidal droplet in 2D. Regarding cell motility, we refer to [40,41] for the use of adaptive finite element methods for simulations in 2D, where the latter work is on vesicles formed by biomembranes.

In comparison with the interface tracking paradigm, computational methods based on interface capturing are generally easier to analyse, to stabilise, and to implement within standard software packages (Matlab, The MathWorks Inc.; COMSOL, The COMSOL Group; DUNE [42]) that nowadays provide parallel mesh management, discretisation tools for assembly, and solvers. The natural way to capture topological changes might also make them beneficial for problems involving cell division or similar. With regards to efficiency, comparative studies of numerical schemes for cell motility do not seem to exist. But accounting for fluid flow can be expected to be a (if not the most) significant contribution to the computational cost. Therefore, comparative studies on the two-phase flow [43,44] give some insight and suggest that the effort of develop tracking methods pays off. But we also note that the parallel test and benchmark computations of 3D problems are missing, and capturing methods involving only bulk differential equations seem more straightforward to parallelise and load-balance. Although generic software packages exist, and highly customised code was developed by experts in numerical methods, to date there are no packages that are easily available to biologists, where models can be adapted and conveniently implemented with ease. Ideally, advanced numerical tools for dealing with changing geometries would be included in tools such as virtual cell [45], which provide a high-level interface for implementing and sharing models.

## Perspectives

• Novel and modern methods to tackle free boundary problems for modelling cell motility are required if we want to account for realistic complex cell shape changes.

• Two main paradigms currently exist, which rely on either tracking the cell boundary explicitly or reformulating the problem using additional fields in which the cell boundary is embedded. Both approaches have been successfully used in 2D. The 3D case is more challenging because no straightforward parametrisation of the free boundary is at hand, and secondly because computations in three dimensions are significantly more costly, particularly if the fluid flow is involved.

• An important future application is the inference of parameters from experimental data which requires numerous simulations and thus highly efficient solvers. Although significant progress has been made, it is still a long way to go before biologists could routinely use such modelling tools to explore different hypotheses in an interactive manner.

## Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

## Author Contributions

B.S. contributed his expertise on numerical methods for interface tracking and capturing, T.B. provided expertise in modelling cell motility.

## Funding

TB received funding from BBSRC grant BB/R004579/1.

## References

References
1
Szabó
,
A.
and
Mayor
,
R.
(
2018
)
Mechanisms of neural crest migration
.
Annu. Rev. Genet.
52
,
43
63
2
Lämmermann
,
T.
and
Kastenmüller
,
W.
(
2019
)
Concepts of GPCR-controlled navigation in the immune system
.
Immunol. Rev.
289
,
205
231
3
Stuelten
,
C.H.
,
Parent
,
C.A.
and
Montell
,
D.J.
(
2018
)
Cell motility in cancer invasion and metastasis: insights from simple model organisms
.
Nat. Rev. Cancer
18
,
296
312
4
Niculescu
,
I.
,
Textor
,
J.
and
de Boer
,
R.J.
(
2015
)
Crawling and gliding: a computational model for shape-driven cell migration
.
PLoS Comput. Biol.
11
,
e1004280
5
Camley
,
B.A.
,
Zhao
,
Y.
,
Li
,
B.
,
Levine
,
H.
and
Rappel
,
W.-J.
(
2017
)
Crawling and turning in a minimal reaction–diffusion cell motility model: coupling cell shape and biochemistry
.
Phys. Rev. E
95
,
012401
6
Alonso
,
S.
,
Stange
,
M.
and
Beta
,
C.
(
2018
)
Modeling random crawling, membrane deformation and intracellular polarity of motile amoeboid cells
.
PLoS One
13
,
e0201977
7
Bretschneider
,
T.
,
Othmer
,
H.G.G.
and
Weijer
,
C.J.J.
(
2016
)
Progress and perspectives in signal transduction, actin dynamics, and movement at the cell and tissue level: lessons from Dictyostelium
.
Interface Focus
6
,
20160047
8
Othmer
,
H.G.
(
2019
)
Eukaryotic cell dynamics from crawlers to swimmers
.
Wiley Interdiscip. Rev. Comput. Mol. Sci.
9
,
e1376
9
Mogilner
,
A.
,
Barnhart
,
E.L.
and
Keren
,
K.
(
2020
)
Experiment, theory, and the keratocyte: an ode to a simple model for cell motility
.
Semin. Cell Dev. Biol.
100
,
143
151
10
Holmes
,
W.R.
and
Edelstein-Keshet
,
L.
(
2012
)
A comparison of computational models for eukaryotic cell shape and motility
.
PLoS Comput. Biol.
8
,
e1002793
11
Neilson
,
M.P.
,
MacKenzie
,
J.A.
,
Webb
,
S.D.
and
Insall
,
R.H.
(
2010
)
Use of the parameterised finite element method to robustly and efficiently evolve the edge of a moving cell
.
Integ. Biol.
2
,
687
695
12
Elliott
,
C.M.
,
Stinner
,
B.
and
Venkataraman
,
C.
(
2012
)
Modelling cell motility and chemotaxis with evolving surface finite elements
.
J. R. Soc. Interface
9
,
3027
3044
13
Tyson
,
R.A.
,
Zatulovskiy
,
E.
,
Kay
,
R.R.
and
Bretschneider
,
T.
(
2014
)
How blebs and pseudopods cooperate during chemotaxis
.
111
,
11703
11708
14
Collier
,
S.
,
Paschke
,
P.
,
Kay
,
R.R.
and
Bretschneider
,
T.
(
2017
)
Image based modeling of bleb site selection
.
Sci. Rep.
7
,
6692
15
Oelz
,
D.B.
,
del Castillo
,
U.
,
Gelfand
,
V.I.
and
Mogilner
,
A.
(
2018
)
Microtubule dynamics, kinesin-1 sliding, and dynein action drive growth of cell processes
.
Biophys. J.
115
,
1614
1624
16
Ma
,
L.
and
Klug
,
W.S.
(
2008
)
Viscous regularization and R-adaptive remeshing for finite element analysis of lipid membrane mechanics
.
J. Comput. Phys.
227
,
5816
5835
17
Neilson
,
M.P.
,
Veltman
,
D.M.
,
Van Haastert
,
P.J.M.
,
Webb
,
S.D.
,
John
,
A.
and
Insall
,
R.H.
(
2011
)
Chemotaxis : a feedback-based computational model robustly predicts multiple aspects of real cell behaviour
.
PLoS Biol.
9
,
1
11
18
Cusseddu
,
D.
,
Edelstein-Keshet
,
L.
,
Mackenzie
,
J.A.
,
Portet
,
S.
and
,
A.
(
2018
)
A coupled bulk-surface model for cell polarisation
.
J. Theor. Biol.
481
,
119
135
19
Woolley
,
T.E.
(
2017
)
Random blebbing motion: a simple model linking cell structural properties to migration characteristics
.
Phys. Rev. E
96
,
012409
20
Lim
,
F.Y.
,
Koon
,
Y.L.
and
Chiam
,
K.-H.
(
2013
)
A computational model of amoeboid cell migration
.
Comput. Methods Biomech. Biomed. Engin.
16
,
1085
1095
21
Taloni
,
A.
,
Kardash
,
E.
,
Salman
,
O.U.
,
Truskinovsky
,
L.
,
Zapperi
,
S.
and
La
,
P.C.
(
2015
)
Volume changes during active shape fluctuations in cells
.
Phys. Rev. Lett.
114
,
1
5
22
Campbell
,
E.J.
and
Bagchi
,
P.
(
2017
)
A computational model of amoeboid cell swimming
.
Phys. Fluids
29
,
101902
23
Wu
,
H.
,
Avila
,
M.
and
De
,
L.P.
(
2018
)
Getting in shape and swimming: the role of cortical forces and membrane heterogeneity in eukaryotic cells
.
J. Math. Biol.
77
,
595
626
24
MacDonald
,
G.
,
Mackenzie
,
J.A.
,
Nolan
,
M.
and
Insall
,
R.H.
(
2016
)
A computational method for the coupled solution of reaction–diffusion equations on evolving domains and manifolds: application to a model of cell migration and chemotaxis
.
J. Comput. Phys.
309
,
207
226
25
Barrett,
J.W
,
Garcke,
H.
and
Nürnberg
,
R
. (
2020
) Chapter 4—Parametric finite element approximations of curvature-driven interface evolutions. In
Geometric partial differential equations—part I
(
Bonito
,
A.
,
Nochetto
,
R.H.
, eds), vol.
21
, pp.
275
423
,
Elsevier
,
Amsterdam
26
Barrett
,
J.W.
,
Garcke
,
H.
and
Nürnberg
,
R.
(
2015
)
Numerical computations of the dynamics of fluidic membranes and vesicles
.
Phys. Rev. E
92
,
52704
27
Osher
,
S.
and
Fedkiw
,
R.P.
(
2001
)
Level set methods: an overview and some recent results
.
J. Comput. Phys.
169
,
463
502
28
Sethian
,
J.A.
and
Smereka
,
P.
(
2003
)
Level set methods for fluid interfaces
.
Annu. Rev. Fluid Mech.
35
,
341
372
29
Chang
,
Y.C.
,
Hou
,
T.Y.
,
Merriman
,
B.
and
Osher
,
S.
(
1996
)
A level set formulation of Eulerian interface capturing methods for incompressible fluid flows
.
J. Comput. Phys.
124
,
449
464
30
Yang
,
L.
,
Effler
,
J.C.
,
Kutscher
,
B.L.
,
Sullivan
,
S.E.
and
Robinson
,
D.N.
and
Iglesias
,
P.A.
(
2008
)
Modeling cellular deformations using the level set formalism
.
BMC Syst. Biol.
2
,
68
31
Shi
,
C.
,
Huang
,
C.H.
,
Devreotes
,
P.N.
and
Iglesias
,
P.A.
(
2013
)
Interaction of motility, directional sensing, and polarity modules recreates the behaviors of chemotaxing cells
.
PLoS Comput. Biol.
9
,
e1003122
32
Kuusela
,
E.
and
Alt
,
W.
(
2009
)
Continuum model of cell adhesion and migration
.
J. Math. Biol.
58
,
135
161
33
Rätz
,
A.
and
Voigt
,
A.
(
2006
)
PDE's on surfaces—a diffuse interface approach
.
Commun. Math. Sci.
4
,
575
590
34
Elliott
,
C.M.
,
Stinner
,
B.
,
Styles
,
V.
and
Welford
,
R.
(
2011
)
Numerical computation of advection and diffusion on evolving diffuse interfaces
.
IMA J. Numer. Anal.
31
,
786
812
35
Moure
,
A.
and
Gomez
,
H.
(
2017
)
Phase-field model of cellular migration: three-dimensional simulations in fibrous networks
.
Comput. Methods Appl. Mech. Eng.
320
,
162
197
36
Shao
,
D.
,
Levine
,
H.
and
Rappel
,
W.J.
(
2012
)
Coupling actin flow, adhesion, and morphology in a computational cell motility model
.
109
,
6851
6856
37
Cao
,
Y.
,
Karmakar
,
R.
,
Ghabache
,
E.
,
Gutierrez
,
E.
,
Zhao
,
Y.
,
Groisman
,
A.
et al (
2019
)
Cell motility dependence on adhesive wetting
.
Soft Matter
15
,
2043
2050
38
Moure
,
A.
and
Gomez
,
H.
(
2018
)
Three-dimensional simulation of obstacle-mediated chemotaxis
.
Biomech. Model Mechanobiol.
17
,
1243
1268
39
Du
,
Q.
and
Feng
,
X.
(
2020
)
The phase field method for geometric moving interfaces and their numerical approximations
.
Handb. Numer. Anal.
,
425
508
40
Marth
,
W.
and
Voigt
,
A.
(
2014
)
Signaling networks and cell motility: a computational approach using a phase field description
.
J. Math. Biol.
69
,
91
112
41
Aland
,
S.
,
Egerer
,
S.
and
Lowengrub
,
J.
(
2014
)
Voigt A. diffuse interface models of locally inextensible vesicles in a viscous fluid
.
J. Comput. Phys.
277
,
32
47
42
Blatt
,
M.
,
Burchardt
,
A.
,
Dedner
,
A.
,
Engwer
,
C.
,
Fahlke
,
J.
,
Flemisch
,
B.
et al (
2016
)
The distributed and unified numerics environment, version 2.4
.
Arch. Numer. Softw.
4
,
13
29
43
Aland
,
S.
,
Boden
,
S.
,
Hahn
,
A.
,
Klingbeil
,
F.
,
Weismann
,
M.
and
Weller
,
S.
(
2013
)
Quantitative comparison of Taylor flow simulations based on sharp-interface and diffuse-interface models
.
Int. J. Numer. Methods Fluids
73
,
344
361
44
Barrett
,
J.W.
,
Garcke
,
H.
and
Nürnberg
,
R.
(
2014
)
Phase field models versus parametric front tracking methods: are they accurate and computationally efficient?
Commun. Comput. Phys.
15
,
506
555
45
Resasco
,
D.C.
,
Gao
,
F.
,
Morgan
,
F.
,
Novak
,
I.L.
,
Schaff
,
J.C.
and
Slepchenko
,
B.M.
(
2012
)
Virtual cell: computational tools for modeling in cell biology
.
Wiley Interdiscip. Rev. Syst. Biol. Med.
4
,
129
140
46
Meinhardt
,
H.
(
1999
)
Orientation of chemotactic cells and growth cones: Models and mechanisms
.
J. Cell Sci.
112
,
2867
2874
PMID:
[PubMed]
47
Lockley
,
R.
,
,
G.
and
Bretschneider
,
T.
(
2015
)
Image based validation of dynamical models for cell reorientation
.
Cytometry A
87
,
471
480
48
Du
,
C.-J.
,
Hawkins
,
P.T.
,
Stephens
,
L.R.
and
Bretschneider
,
T.
(
2013
)
3D time series analysis of cell shape using Laplacian approaches
.
BMC Bioinformatics
14
,
296
49
Neilson
,
M.P.
,
Mackenzie
,
J.A.
,
Webb
,
S.D.
and
Insall
,
R.H.
(
2011
)
Modeling cell movement and chemotaxis using pseudopod-based feedback
.
SIAM J. Sci. Comput.
33
,
1035
1057
50
Elliott
,
C.M.
and
Stinner
,
B.
(
2013
)
Computation of two-phase biomembranes with phase dependent material parameters using surface finite elements
.
Commun. Comput. Phys.
13
,
325
360