Reversible elastic phase field approach and application to cell monolayers

Motion and generation of forces by single cells and cell collectives are essential elements of many biological processes, including development, wound healing and cancer cell migration. Quantitative wound healing assays have demonstrated that cell monolayers can be both dynamic and elastic at the same time. However, it is very challenging to model this combination with conventional approaches. Here we introduce an elastic phase field approach that allows us to predict the dynamics of elastic sheets under the action of active stresses and localized forces, e.g. from leader cells. Our method ensures elastic reversibility after release of forces. We demonstrate its potential by studying several paradigmatic situations and geometries relevant for single cells and cell monolayers, including elastic bars, contractile discs and expanding monolayers with leader cells.


Introduction
Cell and tissue mechanics is an essential element of many physiological processes, including development, tissue homeostasis and wound healing [1,2]. Both single cells and cell collectives are highly dynamic. For animal cells, fluorescence-based experiments have shown that subcellular structures like the actomyosin cortex, lamellipodia and adhesion complexes turn over on the timescale of minutes, despite their function to provide mechanical stability to cells and tissues [3][4][5]. In most developing and even in some homeostatic tissues (notably skin and intestine), there exists a constant flow of cells [6]. Together, these observations suggest that biological systems should be viscous rather than elastic on large time scales, at least in the absence of extracellular matrix [7].
Surprisingly, recent experiments with epithelial monolayers did reveal elastic signatures despite the high cellular and subcellular dynamics. The standard setup in this context is the so-called wound healing assay. In these experiments, cell monolayers of well-defined geometry migrate into free space created by the removal of a straight barrier. This setup has been used to quantify cellular velocity fields, traction forces and intramonolayer tension [8][9][10]. Several new effects have been discovered, including plithotaxis [11] and collective durotaxis [12].
Very important in our context, it has been shown that one can extract a linear relation between stress and strain, thus defining an elastic modulus [13]. This agrees with the a e-mail: f.ziebert@thphys.uni-heidelberg.de results of experiments that stretch free-standing monolayers without adhesion, from which a well-defined elastic modulus can be extracted [14,15]. Later it has been argued that the seemingly elastic signature in expanding monolayers can also be explained as an emergent property of an active fluid with a purely viscous material law [16].
Three phenomena that have received special attention in the context of the wound healing assay are mechanical waves, interface protrusions with leader cells and monolayer organization at boundaries. Mechanical waves in expanding epithelial monolayers were discovered with traction force microscopy and initially explained by repeated cycles of cytoskeleton fluidization and reinforcement [17]. Later it was shown that both an elastic material law [18] or a viscous material law [19] can explain their origin. A similar situation exists for the finger-like protrusions that are often observed to form at the wound margin [8]. Reminscent of fingering in flow cells, these protrusions are often explained by viscous theories that include mechanisms for wavelength selection [20]. Recently, however, it has been shown that these protrusions tend to have a characteristic distance between each other that can be explained by the elastic properties of the monolayer [21]. Experiments with circular wounds and monolayer flow around obstacles have also provided evidence for both elastic and viscous processes in cell monolayers. While flow around a circular obstacle has been shown to be described best by the Maxwell model for viscoelastic fluids [22], the mechanical properties around a gap seem to correspond more to those of the Kelvin-Voigt model for a viscoelastic solid [23][24][25].
Taken together, a growing body of experimental and modelling results suggest that cell monolayers are highly dynamic and reconcile both viscous and elastic signatures. However, the identification of the appropriate material law is often not clear and might strongly depend on the context. Rather than extending the ongoing discussion of viscous versus elastic laws in a continuum mechanics framework, here we suggest a new approach that brings together these two aspects in a different mathematical framework, namely the phase field method. It is nowadays widely used, e.g., for vesicle dynamics [26], actin gel growth [27] or for the crawling motility of single [28,29] or collectives [30] of cells, but to date has not been applied to the (bulk) elastic aspects of cells and cell monolayers. For these elastic aspects, we resort to standard elasticity theory, which here we turn into a dynamic description by coupling it to the phase field. In fact the phase field method has already a tradition for problems involving elasticity, especially for fracture mechanics [31,32] and stress-induced instabilities [33,34]. However, because they address irreversible problems like fracture, these existing implementations are typically not reversible under a removal of the forces and stresses. They hence are not adapted to the biological situation described above, which requires the combination of dynamics and reversible elasticity. Reversibility has been shown e.g. by optogenetics, when cells return to a homeostatic level of contraction after transient stimulation of motor activity [35].
The aim of this work is to define and test a novel elastic phase field approach that is reversible under release of forces and stresses, and to apply it to few simple model problems suggested by single cell and cell monolayer experiments. It is organized as follows. First, in sect. 2 we briefly review the phase field method, how elastic stresses have to be defined, and how previous studies incorporated elastic effects. We show that reversibility, one of the hallmarks of elasticity and the fact that after a release of the applied forces the system should go back to its initial state, is not recovered in the existing approaches. Therefore, in sect. 3 we propose a novel way to couple elasticity to the phase field via an imbalance of forces at short times that allows to drive the phase field interfaces. After testing the approach in one dimension, in sect. 4 we extend the approach to two-dimensional sheets and show that both reversibility and the Poisson effect, another hallmark of elasticity, are captured. Finally, in sect. 5 we discuss several applications of the method in the context of cells and cell monolayers, namely: a contractile cell (or cell sheet) adhering to a substrate, a hole in a contractile cell monolayer, a contractile cell pinned via strong focal adhesions, and the formation of a finger-like protrusion driven by a localized force at the boundary, cf. sketches (a)-(d) in fig. 1. We conclude by discussions and conclusions in sect. 6.

Standard elastic phase field approach
Partial differential equations with moving boundaries are difficult to solve numerically since the position of the latter has to be tracked to impose the boundary conditions at ev- ery time point. An alternative, efficient solution strategy is the phase field method, originally developed for solidification and crystal growth [36,37]. In this method, the sharp boundaries are replaced by diffusive ones at the cost of introducing an additional auxiliary field, the so-called phase field ρ(r, t). It distinguishes two bulk "phases", assigning to each one a constant value, with smooth transition regions in between, identified with the boundaries or interfaces. This method has been now widely applied to many physical systems, such as fracture mechanics [31], vesicles in flow [26], actin gel growth [27] and cell motility [29]. Detailed descriptions of the application of this approach to a variety of problems can be found in specialized reviews, e.g. for microstructure evolution [38], solidification [39], soft matter systems [40] and biological cells [41]. The method associates a free energy functional with free energy density f (ρ) to the phase field where D describes the cost of creating interfaces and sets the width of the diffusive interfaces (∝ √ D), and g(ρ) = ρ 2 (1 − ρ) 2 is a double-well potential with minima at ρ = 0 and ρ = 1. The choice of the double-well potential is not unique, but the two phases should be minima and the given choice is the simplest one. Assuming relaxational dynamics, the dynamic equation for the phase field is then given by the functional derivative of the above energy functional Let us now assume that both phases, described by the regions where ρ = 0 and ρ = 1, respectively, are elastic materials. In order to couple the phase field to the elastic variables, first an expression for the stress field in the whole domain is needed. Assuming both materials are linearly elastic and isotropic, for each phase one has Hooke's law [42], with the strain field ij and Lamé coefficients μ and λ that can be different in the two materials, and the respective elastic energy density To define the stress tensor in the whole domain, one interpolates by writing with σ 0 ij and σ 1 ij the stresses in the materials described by ρ = 0 and 1, respectively, and an interpolation function h(ρ) = ρ 2 (3 − 2ρ) [27,34]. This function is, again, not unique, but it should have values 0 and 1 for the phases 0 and 1, respectively. The given choice has in addition a minimum for 0 and a maximum for 1.
Mechanical equilibrium is described by ∇ · σ + f = 0 [42], with f an external force density. The same relation should hold for the stress in the phase field sense, i.e.
∇ · Σ + f = 0, (6) with Σ being defined in eq. (5). Note that the divergence operator then generates several terms, including derivatives of the interpolation field h(ρ) and hence of ρ. Similar to the stress field, the elastic energy can be interpolated by replacing σ ij by Σ ij in eq. (4), resulting in The standard approach to couple elasticity to the phase field dynamics is to add the interpolated (and hence phase field-dependent) elastic energy to the free energy density of the phase field defined in eq. (1), f tot = f (ρ)+f el (ρ) [27,[32][33][34][43][44][45]. Again performing the functional derivative, this results in where the new term can be interpreted as a driving force for the phase field due to elasticity. Since ∂h(ρ)/∂ρ is a positive function peaked at the phase field interface, the phase field (associated to phase 1) will advance if f 1 el > f 0 el and retract otherwise. For example, consider the stressinduced surface instability (Grinfeld instability) [33,34], where phase 0 is the outside, non-material phase (hence f 0 el = 0) and phase 1 is under stress (hence f 1 el > 0). Consequently, phase 1 will grow -in fact, because of incompressibility, by undulating its surface-to release the stress. Note that eq. (8) has to be solved together with eq. (6) describing mechanical equilibrium.
The standard elastic phase field approach describes situations in which elastic deformations drive changes in the position of a domain. However, it does not describe any elastic relaxation back to the original configuration. The reason is that in eq. (8) only the elastic energy enters, which is quadratic in deformation (or, equivalently, in stress). Hence a sign change of the applied force f does not lead to a sign change in the driving force term (last term of that equation) and the interface does not go back. We should note that for the problems treated with eq. (8) so far, this "non-reversibility" was not a problem: in fracture mechanics problems of brittle materials, when the material is broken, it does not close under release of the force [32,43]. Similarly, in ref. [27] the growth of an actin gel was modeled, but the healing/depolymerisation of the gel was not considered. However, for many other situations of interest, and in particular for cell monolayers as extensively discussed in the introduction, a phase field approach with elastic reversibility would be highly desirable.

Reversible elastic phase field approach
We now propose an alternative approach to couple the elasticity to the phase field dynamics, which is able to recover the reversibility of linear elasticity. To make it more transparent, we first describe it in one spatial dimension (1D). The generalization to two dimensions (2D) is then described in the next section.
The idea is as follows: the phase field should be regarded as a numerical method to describe deforming or moving boundaries. As such, although this may be convenient, its dynamics does not need to be the functional derivative of a potential. We propose to keep the phase field potential defined in eq. (1), but to implement the coupling to elasticity via forces instead of energy. We hence write for the phase field dynamics The last term should be interpreted as follows: the term in brackets is the sum of the internal elastic force (1D divergence of stress) and the external force density F (we use F here to avoid confusion with the free energy density). At mechanical equilibrium, this sum is zero and hence the whole term vanishes and the interface described by the phase field is stationary. If there is an imbalance of forces, however, the interface will be advected in the direction of this force imbalance due to the term ∂ρ/∂x. This proceeds until the phase field has attained a new shape fulfilling force balance, where mechanical equilibrium again holds and the movement stops. α is a coupling parameter (a mobility or inverse friction), setting the characteristic velocity of this movement.
The question now arises how to implement the elastic movement during the times of force imbalance. One could think of various choices, depending on the system to be modeled. In the soft matter and biology context, one would argue that motion should be overdamped. The simplest assumption then is a relaxational (overdamped) dynamics for the displacement field, where ξ is a friction coefficient setting the relaxation timescale into mechanical equilibrium. The dynamics of eqs. (9) and (10) then is as follows: let us assume we start with a stationary phase field (for the standard phase field in 1D, this corresponds to a tanh(x)profile) and vanishing displacement field u(x) = 0. If we apply a force F (specified in detail below) at the phase field boundary, as the stress in eq. (10) is zero, u will increase with a certain timescale proportional to ξ. This will lead to a build-up of stress until mechanical equilibrium is reached via eq. (10) and at the same time the phase field interface moves via eq. (9) because of the transient force imbalance.
In order to ensure the independence of the phase field dynamics on ξ, the coupling parameter α should be proportional to the inverse of ξ. This becomes clear by recognizing that the right hand side of eq. (10) is the driving term of eq. (9). We chose α = 1/ξ since this yields indeed consistent results, meaning that the phase field interface moves as far as the displacement field at the initial boundary indicates, independent of ξ. In turn, ξ should be chosen as small as possible to minimize computational costs, although for a specific system of interest, the exact choice should depend on microscopic details. We use an implicit scheme (Crank-Nicolson) to solve eq. (10), while the phase field dynamics is solved by the Fourier pseudospectral method.
We now test the proposed modeling approach with an instructive example, namely an elastic bar. In 1D, with the ρ = 1 phase describing an elastic material with Young's modulus E and the ρ = 0 phase being empty space, the stress is simply given by The externally applied force density in eq. (10) has to be modeled in the phase field sense. That means, if one wants to apply a force at the boundary, which is diffuse, also the force has to be smeared out. As it should converge to a delta-function in the sharp-interface limit, we choose a Gaussian force profile with an amplitude F max and a width given by the variance υ 2 . For consistency, the latter should be chosen in a similar range as the width of the phase field interface.
As the force should be always located at the boundary, its position is phase field-dependent. We hence assume the mean x 0 (ρ) to be always located at the position where ∂ρ/∂x is maximum, which is the simplest identification of the boundary position in the phase field approach. Finally, there is one more thing to take into account: if the interface moves backwards, i.e. if the domain shrinks, one has to assure that the displacement field is suppressed in what will become the outside of the domain, where no material and hence no displacement field exists. This problem is easy to remedy by adding a suppression term to eq. (10) for the displacement field One sees that the new suppression term is active only in the outside, is linear in u and has a rate γ.
Since the suppression term is active also in the interface region, as a consequence the total force actually applied to the domain is not the integral F (x)dx of the Gaussian force. Rather, the integral of the last two terms in (12), should be regarded as the total applied force, which we will call "effective force" in the following. A numerical implementation of our scheme shows that it is indeed able to describe the reversible dynamics of a 1D domain after the release of an applied force as predicted by linear elasticity. This is demonstrated in figs. 2 and 3, where eqs. (9) and (12), i.e.
were solved numerically for a 1D domain of length 2L 0 under both extension and compression, i.e. with a force couple F applied (with opposite sign) on both sides. Figure 2 shows the reversibility: the black symbols and solid black curve show the phase field before and after the action of the force. Because they lie on top of each other, our model is reversible. The dashed curve shows the phase field during the action of the force, after its movement has stopped and before relaxation by removal of the force. The corresponding displacement fields are shown in red: before application of the force (symbols), during application of the force (dashed, linear profile in bulk) and after force removal (solid flat line). Again it behaves exactly as expected for a reversible model. Figure 3 compares the obtained stationary, mechanical equilibrium state to the analytical solution of the elastic problem. Here the red solid curve is the numerically obtained displacement field. The symbols show the analytical solution of the underlying elasticity problem, which is u(x) = F ext x/E, where F ext is given by the effectively applied force, eq. (13), obtained using the numerical solution without any fitting parameter. In the bulk, the numerical displacement is indeed linear with the correct slope, until it rapidly transits to zero at the boundary (since it has to vanish in the outside where there is no material). The blue curve shows the shape of the integrand of eq. (13), normalized by the total effective force F ef f . The peak at the interface corresponds to the applied force density and the dip outside to the local suppression of u. Obviously the dip cannot be neglected and is important to obtain the excellent agreement with the analytical solution.
Apart from applying forces at the boundaries, we also checked the approach under the action of a constant homogeneous internal stress inside the material. This also works and is explained in more detail for the 2D case of a contractile cell sheet in sect. 5.1.

Two-dimensional elastic sheets
Having detailed the model for the case of a 1D elastic domain, it is not difficult to generalize to higher dimensions and to study more complex geometries. As we are aiming at quasi-two-dimensional cell monolayers and cell sheets, the generalization to 2D involves the fields ρ = ρ(x, y, t) and u = (u x (x, y, t), u y (x, y, t)) and reads Here a new term appears that is related to κ, the 2D curvature of the interface. The curvature can be calculated from the phase field via κ = −∇ · (∇ρ/|∇ρ|). In practise, its calculation can be restricted to the region around the interface, because the prefactor rapidly goes to zero away from the interface. This term is needed to correct for the known effect of the used phase field potential to shrink curved interfaces [46][47][48]. Otherwise eq. (17) is formally identical to the 1D version. We note that the stress as defined by eqs. (5), (3) now involves the 2D plane stress Lamé coefficients where d is the thickness of the layer, assumed to be much smaller than the in-plane dimensions.
As in 1D, we again suppress displacement outside the moving elastic domain. We now write the damping term in eq. (17) as γ(F ) = γ 0 + γ 1 (F ). The use of two terms make the procedure more flexible and in particular allows us to adapt to the Poisson effect, which leads to movement at free boundaries that are not directly pulled by an external force.   The results depend on the choice of the damping coefficients for the displacement outside the elastic domain. Figure 4(a) shows the result for γ 1 = 1 (suppression at regions with large force) and γ 0 = 0 (no global suppression). One can see that the Poisson effect is indeed very well captured. Figure 4(b) shows the result including a small global suppression, γ 0 = 0.35 and γ 1 = 0.65. The Poisson effect is still captured well, but the displacement in the normal direction is not perfectly linear (probably also due to edge effects) and slightly underestimated close to the boundary. Despite this slight disadvantage, the second version is favorable in regard to reversibility. If the force is released, for the suppression parameters as in (b) the system is completely reversible, while for the version shown in (a) the phase field retracts back but leaves some displacement artefacts, that can even accumulate to a peak in the outside region. In conclusion, our approach yields reversibility, but it comes at the cost of resolving the Poisson effect less accurately. Obviously, if needed, the latter accuracy can still be improved by decreasing D and by using a finer grid resolution. merical realization in 2D plane stress, Fmax already takes into account the layer height d. Hence to compare the numerical and analytical results, F ef f has to be devided by d.  Figure 5 shows the full 2D displacement field, u = (u x , u y ), as arrows. In addition, the domain's initial interface (defined by the ρ = 0.5-isocurve) is shown in blue and the domain's interface after application of the force and reaching mechanical equilibrium in red. Clearly, beyond the mid sections shown in fig. 4, the behavior is well captured: The left and right boundaries move under the acting force and the upper and lower boundaries, where no force is acting, retract due to the Poisson effect. Only the corners lightly round up, which however is a known feature for the phase field method, related to the wall energy being removed only to leading order by the correction term (∝ κ) in eq. (16). In general, as most geometries of interest in soft matter and biological systems do not show corners, this is not a strong weakness.
Finally, fig. 6 investigates different values for the Poisson ratio. From the figure (obtained for the case γ 0 = 0) one sees that for lower values of the Poisson ratio, the simulations agree less well with the analytical results. Again this is not a strong weakness as in soft matter and biological systems, we are usually interested in close-toincompressible materials. Moreover the overall trends are captured well for lower Poisson ratios as well. For ν = 0 the method experiences problems due to insufficient resolution and, probably, edge effects. As expected, for the case with global suppression (γ = 0) the agreement worsens. Nevertheless, alternative approaches are either not reversible at all (as discussed above) or only work for the incompressible case [49].

Application to cells and cell monolayers
Having validated the approach for passive 2D elastic domains deformed by external forces, we will now discuss its usefulness to describe biological systems that deform themselves due to active internal forces. More specifically, we will now use the developed elastic phase field approach to investigate several archetypical geometries, as aleady sketched in fig. 1. Cells and cell monolayers continuously generate internal forces to probe the mechanics of their environment, to stabilize their interactions with the environment and to control their shape and mechanics. In order to address this biological situation, we have to include two important aspects in our modelling approach, namely internal contractility and adhesion to the underlying surface. Regarding contractility, single cells and cell sheets are known to be contractile due to myosin II motors actively contracting the actin cytoskeleton [50]. For epithelial monolayers, this contraction is coherent over very large distances (hundreds of micrometers) due to the strong cohesion provided by the cadherin-based adherens junctions [51,21]. Such a global contractility can be modeled on a coarse scale by adding an isotropic contractile stress σ 0 < 0 to the constitutive relation from eq. (3): For simplicity, we here consider only a homogeneous, timeindependent contractile stress, but clearly the method is also applicable for time-dependent and spatially inhomogeneous active stresses σ 0 (x, y, t).
Regarding adhesion, cells and cell monolayers are connected to the substrate through a layer of integrin-mediated adhesions. This adhesion layer allows for exchange of information and for mechanical coupling between cells and their substrates [52]. Cell-matrix adhesion to the substrate has been extensively characterized with traction force microscopy [53,54]. A simple approach to model an adhesive soft interface that transmits force is the so-called elastic foundation [55][56][57][58]. Here one assumes a homogeneous surface coverage with springs of spring constant density Y (measured in N/m 3 ) positioned between a stiff substrate and the cell. This leads to an additional contribution to the force balance, eq. (6), reading In the dynamics, this term simply enters the r.h.s. of eq. (17). Finally, the case of very strong adhesions can also be modeled by defining regions A where both the phase field and the displacement field are pinned, i.e. via the boundary condition With these rather simple ingredients added, we are now prepared to model the specific geometries relevant for cell and cell monolayer statics and dynamics, as sketched in fig. 1.

Contractile cell adhering to a substrate
We first address the contractile elastic disk on an elastic foundation, cf. fig. 1(a), which constitutes a simple model for an adherent cell or cell monolayer [55][56][57][58]. This model has been used before to explain the characteristic localization of traction forces at the cell or cell layer periphery. Important in our context, an analytical solution is known [55] and the radial displacement field is given by Here l = Ed/Y (1 − ν 2 ) is the localization length, which combines the moduli and dimension of the disk with the substrate stiffness, i.e. the spring stiffness density, Y . I 0 and I 1 are the modified Bessel functions of the first kind. Figure 7 shows results of phase field simulations for an elastic disk of radius r 0 = 25 under the action of a strong homogeneous contractile stress and for different values of the substrate spring stiffness density Y . The agreement of the displacement field, which by symmetry is only radial, with the analytical solution is very good inside the bulk material. At the interface, the displacement smoothly crosses over to zero in the no-material region.
Note that to make this comparison, in eq. (22) one has to account for two effects: first, the initial disk radius r 0 has to be replaced by r c , the current radius, since the displacement field moves with the phase field. Second, in analogy to sect. 3 where the damping γ reduced the applied force, here the contractile stress is reduced. Hence similar to the effective force given in eq. (13), here an effective (radial) stress is acting that is given by where r c is again the current radius of the disk. Note that to get this good agreement, the contribution of the substrate spring stiffness density Y in the interface region has also to be accounted for 2 . The inset of fig. 7 investigates the reversibility. It should be noted that the used stress value is rather high, so deformations are only small for large stiffnesses Y . Nevertheless, the reversibility is well captured down to Y = 250 implying displacements of order 10%. We nevertheless should give the warning that for too large displacements (i.e. large σ 0 and/or small Y ), reversibility may be no longer complete, because the phase field model includes non-linear effects associated to the large scale motion of its interface that go beyond linear elasticity.
Let us briefly discuss the orders of magnitude of the parameters used. The equations we solve have been made dimensionless by assuming a typical length scale of 1 μm (for cell layers we use 10 μm) and a force scale of 1 pN. In the example shown in fig. 7, this corresponds to a strongly spreading cell of diameter 50 μm, thickness 1 μm and modulus E = 5 kPa creating a contractile stress σ 0 ∼ 3.15 kPa. The values for the spring stiffness densities Y are in the range 0.25-2 nN/μm 3 . This demonstrates that the method is working well in the range of parameters relevant for the aimed at biological problems. The actual values for the following examples are in a similar range.
If adhesion to the substrate is not considered the analytical solution for a hole of radius r 0 in an infinite monolayer under an isotropic stress σ 0 is given by [64] u r (r) = 1 + ν E yielding a long-ranged 1/r-decay. In the numerical implementation we have to apply periodic boundary conditions, corresponding to a regular array of holes. The result for the radial displacement field for a hole of r 0 = 10 is shown in fig. 8 for three different system sizes L = 200, 300, 400 (colored curves) and compared to the analytical solution (dashed). One can see that the -most interestingpart of the displacement close to the hole is well captured, while the long-ranged decay suffers from finite size effects that become smaller for larger system sizes. The inset of fig. 8 again confirms the reversibility after the release of the stress. Note that to make this comparison -as already explained in the previous section-in eq. (24) one again has to use the current radius r c instead of r 0 and σ 0,ef f as given by eq. (23), where Y = 0 in the present case.

Contractile cell pinned at focal adhesions
Studying cells on micropatterned adhesive substrates has a long tradition in cell biology and biophysics [65][66][67]. Today this approach is used on a routine level to mimic the behaviour of different cell types in their physiological environment, which is more structured than a glass or plastic dish [52]. We envision that our new model approach can be used in the future to simulate mechanosensitive cell behaviour in such environments.
With this motivation, we next studied another test case as sketched in fig. 1(c): A cell was pinned to the corners of a (25 × 25) square by strong focal adhesions. The latter were implemented as being centered around the four corners, where we draw circles of radius r A = 3 wherein the rigid boundary condition, eq. (21), was applied. We then allowed the cell to contract under an isotropic stress (here σ 0 = 400).
In mechanical equilibrium the cell's boundaries displayed invaginated arcs as shown in fig. 9(a) in red (corresponding to the ρ = 1/2−isocurve). The displacement field is shown as arrows and the adhesion sites, where the cell is pinned, as blue circles. To get a scalar quantification of the stress inside the cell, we calculated the von Mises stress, which is defined (in plane stress) as Note that in the phase field sense, σ ij was replaced by the interpolated stress Σ ij . The result for the cell shape shown in fig. 9(a) is given in fig. 9(b), where one can nicely see the expected stress accummulation close to the pinned corners.
We also checked that taking a small adhesion to the substrate (with a spring stiffness density Y ) into account -in addition to the pinning at the corners-does not change the overall qualitative picture (not shown). In contrast, when we implemented only substrate adhesion via a spring stiffness density Y , but no pinning, the squareshaped cell did just contract homogeneously (in an affine fashion), i.e. without displaying invaginations, and the stress was highest in the center, as expected (not shown).

Finger formation in wound healing
As a last test case, we study finger formation at the boundary of cell monolayers, cf. the sketch in fig. 1(d). Such finger-like protrusions emerge spontaneously in wound healing assays [8,10] and have attracted much attention recently, both experimentally and from the modeling side [68][69][70]. Especially, it has been shown recently [21] that the elasticity of the cell sheet is a determinant of the emergence and the spacing of the fingers.
To study the emergence of fingers in our new modeling framework, we initialize a semi-infinite cell monolayer with a straight boundary (at a fixed x 0 ). We then apply a localized force at and perpendicular to the boundary, given in generalization of eq. (11) by Note that again x 0 = x 0 (ρ), i.e. the position where the force is acting moves with the phase field and hence is always acting at the boundary. This force mimics the effect of a cell having decided to spontaneously migrate to invade the free space not belonging to the epithelial monolayer, dragging followers with it.
Our numerical results are shown in fig. 10. Panel (a) shows the boundary of the monolayer in red and the displacement field as arrows. Panel (b) is the calculated von Mises stress, showing stress localization in the layer close to the finger tip. Note that we added a small contractile stress and a weak adhesion to the substrate here, otherwise the displacement propagates even further inside the layer. As in the former cases, for small forces the finger formation is completely reversible.

Discussion and conclusions
Cells and cell monolayers are both elastic and dynamic at the same time, making it very challenging to develop appropriate mathematical models. The phase field approach is very suitable for describing moving interfaces and versions accounting for elasticity have already been developed. However, existing elastic phase field approaches are not reversible under release of forces. As this is crucial in the biological context, e.g. when a protrusion first forms and then relaxes again in a wound healing assay, here we introduced a novel approach. It is based not on a total phase field energy that includes the elastic energy, but rather implements elasticity on the level of forces.
In order to validate the new approach, we performed several tests in 1D and 2D. The Poisson effect is captured best for close to incompressible materials. Care has to be taken when using strong damping outside of the domain (to ensure good reversibility) and when implementing an elastic foundation, since both involve a rescaling of the forces/stresses for finite widths of the phase field interface. Accounting for this, we have shown that the method is in agreement with all tests against analytical results performed and is completely reversible for not too large forces and stresses. Reversibility may become only partial (i.e. the system does not go back completely to its initial state) for higher forces or stresses. This is to be expected since the phase field moving under the action of elastic forces is an effect going beyond linear elasticity, and the more so, the further the phase field boundaries move.
We applied the method to several standard situations that are often studied experimentally for both single cells and cell monolayers. Important features of biological systems, namely active stresses generated inside the layer as well as both weak adhesion and strong pinning to an underlying substrate can be integrated easily into the method. All tests worked well, including reversibility, and several features observed experimentally were well captured, such as the appearance of invaginated arcs and stress focusing for strongly pinned contractile cells.
In the future, the method should prove very useful to investigate dynamically self-organized forces and stresses, especially concerning contractile cells pinned at focal adhesions or finger formation and dynamics at monolayer boundaries. Building on existing cellular phase field models, as additional features relevant to cell dynamics one could implement actin filament orientation [29], concentration fields [28], the effects of biochemical signaling like the Rho-pathway [71], as well as adhesion and traction force dynamics. The latter have been modeled previously within the phase field approach [72,73] by introducing reaction-diffusion kinetics for the engaged adhesive bonds, transmitting traction forces on an elastic substrate, while elasticity of the cell was disregarded. In the framework proposed here, the simplest approach would be to let the distribution of engaged adhesive bonds modulate the substrate's spring stiffness density. Furtheron, in view of the viscoelastic flow behavior of monolayers migrating around an obstacle [22], a generalization of the approach to different viscoelastic models [49] would be highly interesting. We also envision extending our approach to multicellular situations in which single cell resolution is required, by using different phase fields for different cells [30,[74][75][76][77]. The approach could be generalized to three dimensions, to model e.g. cell spheroids or cells moving in strong confinement [78]. In summary, the new method of reversible elastic phase fields introduced here is expected to find many interesting applications in modelling biological systems.
Open Access funding enabled and organized by Projekt DEAL. We thank Chaouqi Misbah, Karin John and Igor Aranson for useful discussions. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2181/1 -390900948 (the Heidelberg STRUCTURES Excellence Cluster). USS is a member of the Interdisciplinary Center for Scientific Computing (IWR).

Author contribution statement
USS and FZ designed the research. All authors were involved in developing the method and its applications. RC developed the numerical code and analyzed the data. FZ supervised the research. All authors wrote the paper together and approved the final manuscript.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.