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.

## 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\u2192$ (Figure 2). Numerical methods based on this paradigm then compute trajectories of points with co-ordinates $x\u2192(t)$ on the interface at timepoints $t$ by integrating their velocity, i.e. solving $\u2202tx\u2192(t)=v\u2192(x\u2192(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.

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\u2192$, 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 $\u2202t$ is then replaced by the *material derivative*$\u2202tv\u2192$ 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\u2192(t)$, i.e. $\u2202tv\u2192c(x\u2192(t),t)=ddtc(x\u2192(t),t)$. A typical system of reaction–diffusion equations for a vector of species densities $c$ on the evolving cell membrane $\Gamma $ 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 $\Gamma $, 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).

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 $\varphi $ 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 $\u2207\varphi $, which is normal to the level set, is given by $\u2212\u2202t\varphi /|\u2207\varphi |$. If now a velocity field $v\u2192$, by which material points on the membrane move, is known (or can be computed) then we obtain the identity $\u2202t\varphi +v\u2192\u22c5\u2207\varphi =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\u2192$ have to be expressed in terms of the level function $\varphi $ in a suitable way. The normalised gradient $\u2207\varphi /|\u2207\varphi |$ 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 $\varphi $. In particular, the force balance in Figure 2C can be written as a partial differential equation for $\varphi $. 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 $\u2202t\varphi +v\u2192\u22c5\u2207\varphi =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.

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 $\phi $ 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 $\epsilon $ 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 $\phi $, and the corresponding force approximates $\sigma \kappa \u2192$, 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 $\phi $. The force balance as in Figure 2C thus can be approximated by the bulk partial differential equation for $\phi $ in Figure 3B for which again standard numerical approaches can be used.

*.*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 $\phi $ 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

*Dictyostelium*