A Finite Element Method for a Fourth Order Surface Equation With Application to the Onset of Cell Blebbing

A variational problem for a fourth order parabolic surface partial differential equation is discussed. It contains non-linear lower order terms, on which we only make abstract assumptions, and which need to be defined for specified problems. We derive a semi-discrete scheme based on the surface finite element method, show a-priori error estimates, and use the analytical results to prove well-posedness. Furthermore, we present a computational framework where specific problems can be conveniently implemented and, later on, altered with relative ease. It uses a domain specific language implemented in Python. The high level program control can also be done within the Python scripting environment. The computationally expensive step of evolving the solution over time is carried out by binding to an efficient C++ software back-end. The study is motivated by cell blebbing, which can be instrumental for cell migration. Starting with a force balance for the cell membrane, we derive a continuum model for some mechanical and geometrical aspects of the onset of blebbing in a form that fits into the abstract framework. It is flexible in that it allows for amending force contributions related to membrane tension or the presence of linker molecules between membrane and cell cortex. Cell membrane geometries given in terms of a parameterization or obtained from image data can be accounted for by the software. The use of a domain specific language to describe the model makes is straightforward to add additional effects such as reaction-diffusion equations modeling some biochemistry on the cell membrane. Some numerical simulations illustrate the approach.


INTRODUCTION
We present and analyse a finite element approximation to parabolic fourth order surface partial differential equations of the form for a vector field u = (u 1 , u 2 , u 3 ) : Ŵ 0 → R 3 on a surface Ŵ 0 that is the smooth boundary of a bounded domain in R 3 . Here, Ŵ 0 is the Laplace-Beltrami operator and ψ ′ and k are given Lipschitz continuous functions (see section 2.2 for precise definitions and assumptions). We also present a software framework that facilitates a convenient implementation of specific problems and report on numerical simulations on an application in cell biology. Our investigation is motivated by cell blebbing, which refers to the detachment of the plasma membrane from its actin cytoskeleton and the fast formation of a spherical protrusion. This is then followed by a slower reformation of the actin cortex close to the deformed part of the membrane and a retraction phase [1,2]. The phenomenon is observed in various processes including apoptosis, spreading, migration, division, embryonic development, and viral entry (see Ponuwei [3] for a recent overview), and is the subject of significant ongoing research. In Tyson et al. [4] and Collier et al. [5] a mechanical model based on ideas in Young and Mitran [6] and Strychalski and Guy [7] for curves in 2D is presented. It accounts for membrane tension, excess cell pressure, membrane bending resistance, forces due to molecules connecting the membrane with the cortex, and drag due to motion through the viscous ambient fluid, and it is used to study the influence of the initial geometry on the propensity to nucleate blebs. One of our objectives is to extend this model to surfaces in 3D and to perform numerical simulations of the onset of blebbing.
Computational approaches to cell blebbing are based on various methodology for fluid-structure interaction. In 2D, the membrane may be tracked by a closed curve. Forces due to its elastic properties can be computed by approximating the curve with a polygonal chain and using finite difference techniques as in Tyson et al. [4] and Collier et al. [5]. The approach can be interpreted as a finite element method and then analyzed in a variational framework. We aim for a generalization of this approach including its convergence analysis to surfaces in 3D that are approximated by triangulated surfaces. If (viscous) fluid flow inside and outside of the cell is accounted for then these membrane forces can be incorporated into the flow equations with a Dirac delta distribution. Smoothing the delta distribution according to certain principles and approximating the fluid flow equations on a regular bulk grid underpins the immersed boundary method, which is understood as both a mathematical formulation and a numerical approach [8]. Based on these ideas, in Young and Mitran [6], a vorticity-stream formulation for the Newtonian, viscous flow is used, and in Strychalski and Guy [7], a staggered grid finite difference method. As an alternative, there are boundary element formulations that are set up directly on the polygonal chain [9]. In 3D, the fluid flow becomes significantly more expensive, and tracking points on a surface whilst maintaining a good representation becomes much more involved. In Maxian et al. [10], a spectral method for the flow is coupled with two approaches for the surface, one based on a parameterization with spherical harmonics and one based on a piecewise linear representation. In Campbell and Bagchi [11,12] a surface finite element method for reaction-diffusion equations on the cell membrane is coupled with a projection method for the flow. In both approaches the immersed boundary method is used for the coupling.
Alternatively to tracking the membrane, interface capturing methods may be used. We are not aware of such an approach to cell blebbing but in Moure and Gomez [13] 3D simulations of a phase field model for moving cells are presented, which uses isogeometric analysis for the spatial approximation and a second order stable time discretisation involving a two-stage predictor-corrector scheme.
The coupling between the surface terms and the bulk flow may also be realized using variational approaches that are amenable to finite element techniques. In the last chapter of Barrett et al. [14] this is presented for vesicles formed by biomembranes that are governed by the Helfrich energy [15]. Such methods address the quality of the evolving surface mesh but a convergence analysis seems out of reach. However, accounting for the fluid flow with a simple drag force or just relaxing some elastic surface energy with a gradient flow dynamics leads to geometric evolution equations, and for the simplest one, the mean curvature flow, convergence of a surface finite element method has been proved recently [16]. Cell blebbing will lead to a more complicated problem but as our focus is on the onset of blebbing, which involves only small deformations, we can expect reasonable results by parameterizing the membrane position over a reference surface. We choose the initial surface as a reference surface and show that the emerging PDE problem can then be of the form (1). This motivates one of our main objectives, namely to prove weak wellposedness and convergence of a surface finite element method for such type of problems. To this end, we proceed essentially as in Elliott and Ranner [17], where the Cahn-Hilliard equation, a scalar parabolic fourth order equation, on an evolving surface is analyzed. The procedure involves splitting the fourth-order problem into two equation of second order by introducing a support field, namely w = − Ŵ 0 u, that can be interpreted as the curvature of the membrane. The essential challenge is the non-linear second order term in (1). To deal with it, we exploit the above linear equation to show strong convergence of the gradient of u.
The surface finite element method dates at least back to Dziuk [18] for the Laplace-Beltrami operator on stationary surfaces and has seen significant development since [19]. It is a versatile tool that can be used for cell motility modeled by a geometric evolution equation coupled with a systems of surface reactiondiffusion equations [20] and can be coupled with bulk finite elements in one of the adjacent domains [21,22] in a fitted approach in the sense that the surface mesh is the boundary of a bulk mesh. Also approaches for the unfitted case of a surface mesh intersecting a fixed bulk mesh rather arbitrarily have been developed, where we mention the trace finite element method [23] and the cut finite element method [24]. In both methods, PDEs are solved on the intersection of a surface with a, usually, Cartesian bulk mesh on which additional terms are set up for stabilization to deal with small intersection patches. For the case of a PDE on an evolving surface the trace finite element method has been analyzed in depth [25]. All these methods FIGURE 1 | Sketch of some objects defined in in sections 2.1, 4 for the abstract framework and the specific cell blebbing model, respectively. The cell membrane Ŵ(t) is parameterized by u over the initial surface Ŵ 0 which, unlike in this sketch, is closed and of spherical topology. Force densities arising from membrane tension, denoted f tension , point outwards in concave areas (as in the sketch) and inwards in convex areas. Linker molecules (indicated in orange) connect the membrane with the cortex, which is close to the initial surface and not explicitly indicated. They are modeled with a force density that points toward the initial surface when the molecules are stretched. The force vanishes if the linkers break. The pressure f pressure points in the direction of the unit normal ν Ŵ 0 of Ŵ 0 , which is the external unit normal of the cell in its initial shape.
have in common that they are usually not provided by PDE software packages out of the box but require expert knowledge to implement them. One of our objectives is therefore to make the methodology we have developed accessible by providing code that enables users to implement specific problems with relative ease.
Let us briefly summarize our objectives and achievements, and outline the paper: • The surface PDE (1) is analyzed for well-posedness using surface finite element techniques. Stability estimates for the semi-discrete scheme are derived and exploited to show convergence. Under slightly more restrictive assumptions on the quality of the solution, error estimates are shown. The suitable abstract variational problem is thoroughly formulated in section 2.2, and the finite element approach is presented and analyzed in section 3. • The force balance model in Tyson et al. [4] and Collier et al. [5], which is based on ideas in Young and Mitran [6] and Strychalski and Guy [7], is extended from curves in 2D to surfaces in 3D. More precisely, a continuum model is presented such that, when restricting the model to a curve in 2D and discretising the governing equations using standard finite difference methods, the original computational model is obtained. This specific model has been used for some numerical simulations and is derived in section 4. • A software framework for numerical simulations has been developed. It features a high-level interface to implement a problem in the Unified Form Language (UFL) [26], which enables a user to conveniently alter the variational problem. Whilst the overall program control and the time stepping are done at the high level, bindings to the Distributed Unified Numerics Environment (DUNE) [27,28] are used for efficiently discretising and solving the spatial problems, more precisely, the Python bindings to the DUNE-FEM module [29,30]. Section 5 contains the time discretisation, details on the implementation, and numerical simulations, which also demonstrate the extensibility of the software package.
The setup is illustrated in Figure 1, which also shows some forces acting on the surface that are explained in section 4. The dependence on t will usually be dropped in the following. We furthermore introduce the following notation: d signed distance to Ŵ 0 , well-defined in a thin layer around Ŵ 0 , convention: d < 0 inside of (0), extended to a thin layer around Ŵ 0 , surface area element when integrating over Ŵ 0 , We aim for approximating the PDE problem (1) using finite elements and thus require a variational formulation. We say that a function f ∈ L 1 (Ŵ 0 ) has a weak derivative holds true for all smooth functions ϕ with compact support. We also use ∇ Ŵ 0 to denote this weak derivative and write ∇ k Ŵ 0 , k ∈ N, for the k-th derivative. Sobolev spaces on Ŵ 0 are defined by H 0 (Ŵ 0 ) = L 2 (Ŵ 0 ) and For a function η ∈ H 1 and a vector field v ∈ (L 2 ) 3 the definition of the weak derivative yields that For the L 2 "mass" inner product and for the H 1 "stiffness" semi-inner product of vector valued functions we write Based on these Sobolev spaces we will consider the Bochner spaces

Assumptions and Weak Formulation
In order to enable us to approximate (1) with linear finite elements, we split up the fourth order operator by using the curvature The function ψ : R 3×3 → [0, ∞] is assumed to be continuously differentiable with uniformly Lipschitz continuous partial derivative, i.e., denoting with ψ ′ this (3 × 3 tensor-valued) partial derivative we assume that there is a constant C ψ > 0 such This implies that |ψ ′ (A)| ≤ C ψ |A|+C for some constant C > 0. Moreover, we assume that ψ ′ is tangential in the following sense: Note that then ψ ′ (A)κ = 0 because κ points in the normal direction. The PDE (1) involves the term ∇ Ŵ 0 · ψ ′ (∇ Ŵ 0 u), and thus by (2) Ŵ 0 for sufficiently smooth functions v, z. With a slight abuse of notation we write The function k : Ŵ 0 × R 3 → R 3 is assumed to be bounded, measurable with respect to the first argument, and uniformly Lipschitz continuous in the second argument, i.e., there is some constant C k > 0 such that for all y ∈ Ŵ 0 This implies that |k(y, a)| ≤ C k |a| + C for some constant C > 0. The (weak) variational formulation of (1) reads: Problem 2.1. Find u, w ∈ L 2 (0, T; H 1 (Ŵ 0 )) with ∂ t u ∈ L 2 (0, T; L 2 (Ŵ 0 )) such that for all φ, η ∈ H 1 (Ŵ 0 ) and almost all t ∈ (0, T) and such that u(·, 0) = id Ŵ 0 .
Remark 2.2. One could consider a function ψ that, like k, depends on the position y on Ŵ 0 . For the finite element approximation discussed in the next section, we make use of an approximation k h of k (see around (10)). Something similar could be done for ψ as well. This would result in some additional terms that would need to be estimated using a consistency assumption similar to (10). As the procedure is similar to the one for k we omit the details for conciseness.

Surface Triangulations and Finite Elements
The membrane Ŵ 0 is approximated by a family of polyhedral surfaces {Ŵ 0 h } h , each one being of the form where the E are closed, flat non-degenerate triangles whose pairwise intersection is a complete edge, a single point, or empty.
For each E belonging to the set T h of triangles we denote by h(E) = diam(E) its diameter and then identify h = max E∈T h h(E) with the maximal edge length of the whole triangulation. We assume that the vertices of Ŵ 0 h belong to Ŵ 0 so that Ŵ 0 h is a piecewise linear interpolation of Ŵ 0 . We also assume that h is small enough so that Ŵ 0 h lies in the thin layer around Ŵ 0 in which the signed distance function d is well-defined. Furthermore, we assume that Ŵ 0 h is the boundary of a domain that approximates (0) and denote the external unit normal, which is defined on the triangles and thus piecewise constant, with we denote the projection to the tangent space in points on Ŵ 0 h where it exists (i.e., in the interiors of the triangles E ∈ T h ). This gives rise to the piecewise (i.e., triangle by triangle) definition of a surface gradient ∇ Ŵ 0 h on Ŵ 0 h . The same notation ∇ Ŵ 0 h is used again for the weak derivative. We write dσ h for the surface area element when integrating functions on Ŵ 0 h . For the error analysis we have to measure the distance of functions such as u on Ŵ 0 to functions such as the finite element solution on Ŵ 0 h . For this purpose, we assume that for each y h ∈ Ŵ 0 h there is a unique point y ∈ Ŵ 0 such that This bijection gives rise to the lift of any function η : Ŵ 0 h → R to Ŵ 0 defined by Writing µ h for the local change of the surface area element, i.e., dσ h = µ h dσ , integrals transform as A straightforward calculation show that in points where both η and η ℓ are differentiable The following two lemmas on the errors due to the approximation of the surface and on the stability of the lift are due to Dziuk [18], Dziuk and Elliott [31].
Lemma 3.1. The following estimates hold true for some constant C > 0 independent of h: Let also E ∈ T h and E ℓ = {y ∈ Ŵ 0 | y h ∈ E} with y as in (7) for a given y h . The following estimates hold true with a constant C > 0 independent of h and the element E: .
These inequalities generalize to the whole surfaces by summing over the elements.
The standard finite element space used throughout is Note that the identic map of Ŵ 0 h belongs to S 3 h . Bilinear forms corresponding to m and s are defined for finite element functions R h , Z h ∈ S 3 h on the triangulation by and we will also use again the notation s h (ψ ′ ; For the discrepancy to the forms on Ŵ 0 we note the following result: We define the Ritz projection h : H 1 (Ŵ 0 ) → S h by It's lift is denoted by π h (ξ ) = h (ξ ) ℓ and has the following approximation properties: where C > 0 is a constant independent of h and ξ . The projection and the convergence results extend to functions in L 2 H 1 with a pointwise (in time) definition of the projection and with the norms · H k replaced by · L 2 H k , k = 0, 1, 2.

Semi-discrete Problem
In applications, we may only have access to a triangulated surface Ŵ 0 h but not Ŵ 0 , for instance, when Ŵ 0 h is computed from image data. In such cases we may also know functions such as k only approximately. For instance, the specific choice (55) of k for our numerical simulations involves the unit normal ν Ŵ 0 and the cortex position u c , both of which may not be known exactly if we only have Ŵ 0 h rather than Ŵ 0 . However, we have the approximations For the analysis of the abstract model we therefore assume that k is approximated by some function (properly, a h family of functions) k h : Ŵ 0 h × R 3 → R 3 that has the same regularity properties as k. In particular, k h is Lipschitz continuous in the second argument with the same Lipschitz constant C k > 0 independently of h. We define its lift k ℓ h : Ŵ 0 × R 3 → R 3 by k ℓ h (y, a) = k h (y h , a), a ∈ R 3 , with y and y h related as defined around (7). We assume that k h is an approximation of k in the following sense: There is a constant C > 0 independent of h such that for all a ∈ R 3 In the next subsection we will show the following main result: Theorem 3.6. The semi-discrete problems 3.5 are well-posed for all h > 0 small enough.
converge to some functions (u, w) that uniquely solve the abstract variational problem 2.1 and satisfy with some C > 0 that depends on data only.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org

Proof of Theorem 3.6
We generally follow the procedure in Elliott and Ranner [17]. Essential differences consist in the approximation of the data k by k h and the non-linear function ψ ′ of the gradient. To deal with the former, the consistency assumption (10) will turn out sufficient, whilst for the latter we will exploit the relations (6) and (12) to show strong convergence of the gradient of the deformation. Short time existence for (11), (12) is straightforward to show. Estimates are now derived that are, at first, only valid at times of existence but then in the usual way can be used to show existence over the whole time interval by a continuation argument. We therefore state these estimates directly on the whole time interval. We also use the standard notion of C > 0 for a generic constant that depends on the problem data but not on any solution, and which may change from line to line.
Testing with h = U h in (11) and H h = W h in (12) and subtracting these identities yields that Here and in the following we use the Lipschitz continuity of ψ ′ and k h , which implies linear growth [see (3), (4) and the comments after]. Choosing H h = U h in (12) and applying Young's inequality we see that forε > 0, and choosingε small enough we thus obtain from (14) that A Gronwall argument therefore yields the estimate Testing with h = W h in (11) and (12) and then adding these equations yields that With H h = W h in (12) and using Young's inequality again we get for any smallε > 0 that Using the Lipschitz continuity of ψ ′ and k h again we thus can conclude that Choosingε small enough and then applying (15) and a Gronwall argument we obtain the estimate Taking the time derivative of (12) (note that, using that ∂ t U h exists, this equation can be used to show that Noting that and using the Lipschitz continuity of k h again we see that Therefore, with (15) we obtain the estimate These estimates (15)-(17) are now lifted from Ŵ 0 h to Ŵ 0 . We can then apply compactness arguments to deduce the existence of limits (u, w), which we will show to satisfy Problem 2.1. As a first step, the stability estimate (13) will be derived. Using Lemma 3.2 the lifted solutions satisfy the estimates Hence, there are functions u ∈ L 2 H 1 with ∂ t ∈ L 2 L 2 and w ∈ L 2 H 1 such that for a subsequence as h → 0 and these limits also satisfy (18) and (19) and, thus, the stability estimate (13).
Let us now show that (u, w) satisfies (6). For any η ∈ L 2 By the properties of the Ritz projection (Lemma 3.4) we have Thanks to (20) we thus have that J 1 → 0 as h → 0, and similarly J 4 → 0 thanks to (21). Lemma 3.3 together with the estimates (16) and (17) ensures that J 2 → 0 and J 3 → 0 as h → 0. Therefore, (u, w) satisfies the following identity, which implies (6): Next, we show strong convergence of ∇ Ŵ 0 u h . We note that Using again Lemma 3.4 we see that π h (u) → u in L 2 H 1 , and with (20) this implies that K 1 → 0. Regarding the second term we note that thanks to (22) and (12) As both π h (u) → u and u h → u by (21) With w h ⇀ w in the same space we obtain that K 21 → 0. From the definition and properties of the Ritz projection it easily follows that h (ξ ) H 1 ≤ C ξ H 1 with some C > 0 independent of h and ξ ∈ H 1 (Ŵ 0 ). The stability estimate (13), which is already proved, and the estimates (15) and (16) therefore yield that h (u) − U h H 1 is uniformly bounded in h. Using (15) and (16) again for W h and Lemma 3.3 we obtain that K 22 → 0 and K 23 → 0 as h → 0. This finally shows that To conclude the proof of Theorem 3.6 we need to show that (u, w) satisfies (5). For any φ ∈ L 2 Thanks to (21) we have that L 1 → 0 as h → 0. Lemma 3.4 on the Ritz projection ensures that L 2 → 0. It also ensures that h is uniformly bounded in h, and with Lemma 3.3 and (17) we obtain that Analogously one can show that Next, we can write Thanks to (23) and the Lipschitz continuity of ψ ′ we have that thanks to the estimate (16) and Lemma 3.4. In the last term we lift the second integral to Ŵ 0 (recall (8) and (9) for the transformation of the derivative): We can now apply the Lipschitz continuity of ψ ′ and the geometric error estimates in Lemma 3.1 (which imply that Using estimate (18) and that also φ h L 2 is uniformly bounded (follows from Lemma 3.4) we see that M 3 → 0, and we can conclude that For the last term in (5) we note that Thanks to (21) and the Lipschitz continuity of k we have that k(u h ) → k(u) in L 2 L 2 and almost everywhere, so that N 1 → 0 as h → 0. The second term converges to zero thanks to φ h → φ in L 2 L 2 . Regarding N 3 , we lift the second integral to Ŵ 0 : Using now the consistency (10) of the approximation of k by k h , the Lipschitz continuity of k h , and the geometric error estimates in Lemma (3.1) we obtain that using estimate (18) and that also φ h L 2 The convergence results (25), (26), (30), and (34) show that (u, w) satisfies (5), which is the limit of (11) as h → 0.
In the next section we show error estimates. These techniques can also be used to show uniqueness of the solution (u, w) to Problem 2.1. We therefore omit the details for brevity. This concludes the proof of Theorem 3.6.

Error Estimates
Deriving error estimates is possible when assuming higher regularity of the solution, henceforth: We will derive error estimates on the triangulated surfaces and for this purpose us the bijection (7) to define the inverse lift of the solution (u, w) to Ŵ 0 h : We use the Ritz projection to split the errors into a projection error ρ and a discrete error θ : Thanks to the regularity assumption (35), the properties of the Ritz projection (Lemma 3.4), and the properties of the lift (Lemma 3.2) error bounds for the projection errors are straightforward: To estimate the discrete errors let h ∈ S 3 h , test (11) with − h and test (5) with −φ h , which is the lift of h . We get that Frontiers in Applied Mathematics and Statistics | www.frontiersin.org

by the definition of the Ritz projection this yields that
Proceeding similarly with (6) and (12) for any H h ∈ S 3 h with lift η h we obtain that The error terms satisfy the following estimates: Proof: To show the estimates, we will frequently apply Lemma 3.3 on the approximation of the bilinear forms, Lemma 3.2 on the stability of the lift, and Lemma 3.4 on the Ritz projection without explicitly pointing it out for conciseness. The Ritz projection commutes with the time derivative thanks to the regularity of u. Therefore The termL 3 is similar to L 3 in (24) but without the time integral and with π h (∂ t u) instead of u h and thus can also be estimated similarly: Furthermore, which altogether yields (40).
We can also split up E ψ : The first termM 3 is similar to the term M 3 in (27), without the time integral and with π h (u) instead of u h . Following the lines of (28) and (29) we obtain that Using that ψ ′ is Lipschitz, the other term is estimated as which together shows (41).
For the third estimate we use the splitting Noting and exploiting the similarity ofÑ 3 with N 3 in (27) we proceed as in (32) and (33) to obtain that Furthermore, which finally yields the estimate (42). The last estimate (43) can be proved analogously to (40). This concludes the proof of Lemma 3.7.
With these estimates we can derive the following estimates for the error: with a constant C > 0 independent of h.
Proof: We proceed as for deriving (15) but start with (38), where we test with h = θ (u) , and with (39), where we choose H h = θ (w) . Taking the difference we obtain that In Lemma 3.7 we absorb the norm of u into C to obtain that and similarly for the other two errors. Using that ψ ′ and k h are Lipschitz we then get that Substituting H h = θ (u) in (39) gives for anyε > 0 that We can thus estimate the terms involving ∇ Ŵ 0 h θ (u) 2 L 2 (Ŵ 0 h ) on the right-hand-side of (44) by terms involvingε θ (w) 2 Choosing nowε > 0 small enough, these terms involving θ (w) 2 L 2 (Ŵ 0 h ) can then be absorbed in the left-hand-side to that altogether By standard interpolation theory (recall that the identic map of the triangulated surface Ŵ 0 h linearly interpolates the identic map of Ŵ 0 ) the initial error satisfies Applying Gronwall therefore yields that From (45) we now see that also Together with (36) and (37) these two estimates conclude the proof of Theorem 3.8.

Force Contributions
Mathematical models that aim to provide insight into the mechanisms that control cell blebbing require an approach to describe the evolving geometry and have to account for various force contributions acting on the plasma membrane. Based on previous ideas [4][5][6][7] we postulate that a force balance of the form governs the cell membrane's shape. The contributions are force densities on the cell membrane and are modeled in a form such that a surface partial differential equation of the form (1) is obtained. Figure 1 gives an idea of the setup. We specifically aim for generalizing the model in Tyson et al. [4] and Collier et al. [5] for curves in 2D to surfaces in 3D but also put the force contributions into the wider literature context.
• Pressure: Building up internal pressure, for instance, by actinmyosin contraction in the cortex, is essential for blebbing (see Charras et al. [32] and references). We write the corresponding force density as where V(u) = max{ Ŵ 0 1 3 u · ν Ŵ 0 dσ , 0} is an approximation of the volume of and p 0 is a pressure coefficient so that p 0 /|V(u)| is the pressure difference between interior and exterior of the cell. This simple assumption of a constant (in space) pressure difference should be sufficient to study the influence of the initial shape on the blebbing propensity. However, for the dynamics of blebs the pressure distribution will be of importance and requires to model the ambient fluid as in Strychalski and Guy [33] and Fang et al. [34] and probably also the actin-myosin biochemistry.
• Coupling between membrane and cortex: Forces arise due to molecules connecting the membrane with the actin cortex, and when the membrane detaches during the blebbing process these linkers break. Figure 1 gives an impression of such a force (denoted by f coupling ). An actin scar is left behind the bleb [4,7] that disintegrates in the longer run, and the cortex reassembles close to the new membrane position. However, as we are interested in short times and the onset of blebbing we assume the cortex to be stationary and positioned a small distance l 0 away from the initial membrane. Connection points of linkers in the cortex are given by u c = id Ŵ 0 − l 0 ν Ŵ 0 . The linker molecules can be modeled as the density of simple springs with parameter k l and assumed to be initially at rest. In a continuum setting this can be modeled with an energy density e coupling = k l 2 (|u − u c | − l 0 ) 2 as long as they are intact. As a critical length u B is exceeded they break. Moreover, when they get closer than a distance u R to the cortex then the repulsion force is increased to prevent any intersection. A model for the force density thus reads with k coupling (y) = k l 1 + k L H(u R − y) H(u B − y) with some constant k L > 0 and the Heaviside function H(r) = 1 if r ≥ 0 and H(r) = 0 otherwise. The assumption of a constant k l , which corresponds to a homogeneous linker, is a simplification. The cell usually is able to control the density and, thus, the local effective linker strength by its biochemistry [35]. We briefly report on a way to account for membranebound biochemistry in section 5.5. • Tension: Lateral tension in curved membranes leads to net forces in normal direction that, in the case of isotropic materials, are proportional to the curvature. The distribution of tension is of relevance in moving cells [36]. It has been noted in Tyson et al. [4] and further investigated in Collier et al. [5] that the propensity of blebs in concave regions of the cell membrane is higher. In such areas, the force due to the membrane tension points outwards (Figure 1 gives an impression, the force is denoted by f tension ) and then adds to the pressure whilst in convex regions it opposes the pressure. We consider an energy density of the form where k ψ > 0 and x 0 ∈ [0, 1] are parameters. Note that |∇ Ŵ 0 u(·, 0)| = |∇ Ŵ 0 id Ŵ 0 | = |P| = √ 2 so that in the case x 0 = 1 the membrane initially is at rest. If Ŵ 0 was a curve in 2D then we would obtain the model in Tyson et al. [4] and Collier et al. [5], which is motivated by chains of linear springs with spring constant k ψ and resting length x 0 . More sophisticated elastic energies for membranes than (50) can be derived from discrete models for meshes formed by springs [37,38] and then usually also account for resistance to bending (discussed further below around (52)). Let us also note that biomembranes, the basic component of cell membranes, rupture when stretched beyond a few percentages. Whilst this seems satisfied by the small deformations during the onset of blebbing that we study, the area increase during full bleb formation can be more significant and require the supply of membrane area [39], which then is likely to affect the tension, too. The (tension) force density is given by minus the variation of our energy (50), • Regularization: The membrane resists bending, though much less than stretching. The corresponding elastic energy may be modeled as in Helfrich [15] and Woolley et al. [40], see [41] for a discussion of minimal approaches. The impact of the bending force on the blebbing site selection and its shape has been found to be significantly smaller than that of the tension [4]. However, we suspect that it is of relevance in the area where the bleb is connected with the cell as there the curvature and its derivative can be high. Therefore, we want methodology and software capable of addressing this aspect for future studies. We here choose a simple linear bending model that may be considered as a regularization with the energy density where k b is a (small) bending resistance coefficient. Minus its variation yields the regularization force density • Viscous drag: Blebbing approaches often explicitly account for the fluid in the interior and exterior of the cell, which then is assumed Newtonian and viscous [6,7,9]. Our focus is on the onset of blebbing rather than on the time evolution and the growth of blebs. We therefore only postulate a viscous drag force density that opposes any membranes movement: where ω is an effective material parameter related to the viscosity of the ambient fluid. Though this approach might yield a good approximation to the dynamics, at least at short time scales, it should rather be considered as another regularization term and a mean to decrease membrane energy in a controlled way. It may also enable us to compute bleb shapes, noting that the drag force is zero if the membrane is at rest, i.e., if ∂ t u = 0. To address questions such as the origin of the fluid in the bleb (from outside of the cell via pores in the membrane [42] or from inside through the cortex [43]) our approach will be insufficient, of course.
With these choices, the force balance (46) yields the PDE.

Non-dimensionalization and Regularization
The PDE (53) has been non-dimensionalized by choosing a length scale U and by using k ψ as an energy density scale. Choosing the time scale = U 2 ω/k ψ then eliminates the viscosity parameter ω. We furthermore define the nondimensional parameters λ b = k b /(U 2 k ψ ), λ l = k l U 2 /k ψ , and λ p = p 0 /(U 2 k ψ ) and note that x 0 and k L are non-dimensional already. Writing again Ŵ 0 , u, u c , u B , u R , l 0 , and V(u) for the respective non-dimensional objects, Equation (53) in nondimensional form and rearranged reads Note that the model in Tyson et al. [4] and Collier et al. [5] is obtained by reducing the dimension of this Equation (54) (i.e., Ŵ 0 is a curve in 2D). The curve then is parameterized by arc-length, and their computational model is obtained by using standard finite difference techniques. Equation (54) can be cast in the form (1) (without loss of generality we can assume that λ b = 1 as we may divide by λ b , rescale in time, and absorb the 1/λ b term into k and ψ defined below) by defining Note neither k nor ψ ′ satisfy the Lipschitz continuity conditions in section 2.2. We still used them for the simulations in section 5 but also have redone some simulations with regularized versions of k and ψ ′ . To ensure that the assumptions in section 2.2 are satisfied, smoothing the Heaviside function and ensuring that the denominators do not degenerate is sufficient. With a small parameter ε > 0 the choice λ coupling,ε (y) = λ l 1 + k L 1 + exp(2(y − u R )/ε) 1 1 + exp(2(y − u B )/ε) and then satisfies the assumptions around (4). The approximation satisfies the consistency assumption (10). Similarly, the choice satisfies the assumptions around (3).

Time Discretisation
We performed some numerical simulations for (54) to illustrate the capability of the theoretical framework that has been presented and analyzed. The variational form with operator splitting in Problem 2.1 is discretised in time with a simple semiimplicit first order scheme as follows: We split the time interval [0, T] into M ∈ N equal parts of size τ = T/M, denote the time steps with t (m) = mτ , and write f (m) = f (t (m) ) for any time dependent fields or functions. with

Implementation
We have solved the above problem using the Python bindings from the DUNE-FEM module [29], which is based on the Distributed and Unified Numerics Environment (DUNE) [28].
DUNE is an open source C++ environment that uses a static polymorphic interfaces to describe grid based numerical schemes. The package provides a large number of realizations of these interfaces including many finite element spaces on structured and unstructured grids. This approach allows for the efficient and flexible simulation of a large variety of mathematical models based on partial differential equations. The Python bindings described in Dedner and Nolte [44] simplify the rapid prototyping of new schemes and models, while maintaining the efficiency and flexibility of the DUNE framework. This is achieved by using the domain specific language UFL [26] to describe the mathematical model and implementing the high level program control within Python. All computationally critical parts of the simulation are carried out in C++ using just in time compilation of the required DUNE components. Consequently, the assembly of the bilinear forms and solving of the linear and non-linear problems is implemented in C++ while the time loop and the input and output of data is carried out using the Python scripting language. More information on the concepts can be found in Dedner et al. [30].
Meshes can be provided using a GMsh file or, as done for this work, by using the internal Dune Grid Format (DGF). All simulations reported on in this paper were performed using a first order Lagrange space over a simplicial, locally adaptive, distributed grid, which can be used for both bulk and surface domains [45]. Bindings for a number of different solver packages are available through DUNE-FEM including the iterative solvers from DUNE-ISTL [46] (used for this work), direct solvers from the SuiteSparse package [47], and a number of solvers and preconditioners from the PetSc package [48]. The simulation results were exported using VTK and visualized using ParaView [49].
In the following we show how to setup the grid and how some parts of the mathematical model are defined within UFL. The full code needed to perform the simulations shown in this paper is available (see the Data Availability Statement).
The first listing shows how to read in a grid for a cell obtained from experimental data that was used for the simulations in In each time step a saddle point problem is solved using a Uzawatype algorithm where a CG method is used to invert the Schur complement as described, for example, in Braess [50]. The main algorithm is implemented in Python calling C++ routines to compute the matrix-vector operations and to solve the inner problem. The time loop with the solver is fairly long and doesn't seem worthwhile listing here but, as stated above already, the whole code is publicly available. A number of tests have been performed for problems with known solutions (u, w) to validate the convergence (rates) of Theorems 3.6 and 3.8. Recall that the choices of the tension term ψ (56), and the coupling term k (55), in the specific model (54) do not satisfy the requirements of the analysis. However, in our simulations, the denominators in these terms did not become very small. Comparative simulations with the regularized choices (58) and (57) with ε = 10 −5 did not reveal any essential difference. More details on computations to support the theoretical results and validation of the code can be found in Nixon [51]. For conciseness, we don't report on these here but focus on an investigation of the parameter space instead.

Brief Parameter Study
One question of interest has been whether surface tension and pressure are sufficient to initiate blebbing without any weakening of the cortex, as found in Collier et al. [5] in 2D. We also further study the parameter space but remark that the simulation results are at a qualitative level. An in-depth discussion involving quantitative information is beyond the scope of this article and left for future investigations. The parameters are steered into the software at the high-level python interface which, in principle, can be conveniently automated at that level for more substantial parameter studies.
We consider an initial shape Ŵ 0 obtained by deforming a sphere of radius one by (all lengths in 10 −6 m ) y = (y 1 , y 2 , y 3 ) → (4y 1 , 4y 2 ,ỹ 3 ), with r = (4y 1 ) 2 + (4y 2 ) 2 . This yields a shape similar to a discocyte (or red blood cell, see Figure 2) with a volume of about V(id Ŵ 0 ) ≈ 1.5 × 10 −16 m 3 and a largest distance of 4.0 × 10 −6 m from the center. Let us briefly motivate this choice for the initial shape. In the concave parts the force due to the membrane tension points outwards and in the convex parts it points inwards. When the cell membrane moves away from the cortex then the coupling force, which is due to stretched linker molecules as long as these are not broken, points inwards. Therefore, if pressure and tension are able to overcome the coupling force and initialize a bleb we will expect this to happen first (and possibly only) in the concave parts of the initial shape.
Parameters for the various force densities vary in the literature, not least due to differing cell types and differences in the models.  [5]. The parameters u R = 7.5 × 10 −9 m and k L = 500.0 were chosen ad hoc but repeating some simulations with k L = 0 (particularly those with higher tension so that the membrane got closer to the cortex) didn't reveal any visual difference. As discussed earlier around (47), pressure inside of the cell has to be higher than outside of it to initiate blebbing. In the literature we find values between 10Pa [40] and 81Pa [4]. We found that a value of p 0 /|V(u(0))| ≈ 2.25Pa was already sufficient in our simulations to break the linker molecules and, thus, to generate a bleb site. With a length scale of U = 1.0 × 10 −6 m the set of non-dimensional parameters is given in Table 1 and was used for simulations unless stated otherwise. A triangulation Ŵ 0 h is obtained by starting with a cube with vertices on the unit-sphere, then diagonally cutting the square faces into triangles, and then bisecting all triangles 14 times such that the longest edge is halved. New vertices are projected to the unit-sphere after each refinement step. After, the above map (61) is applied to the 196608 vertices. Figure 2 gives an impression of a mesh thus obtained but with ten refinements only. The time step size was set to τ = 0.0025. As mentioned, the model serves to study the onset of blebbing and we therefore ended the simulation at T = 2. At that time the final shapes usually weren't at rest yet but the deformations were sufficient to break linker molecules and thus generate a blebbing site. Figure 3 gives an overview of some shapes at the final time for the data set in Table 1 and some variants (see Figure caption for details). The initial shape Ŵ 0 is axisymmetric around its central axis. The initial shape Ŵ 0 h for the computations is not. However, on visual inspection, also the computed shapes look approximately axisymmetric. For this reasons, we compare cuts through the centers for more insight. Figure 4 displays the slices through the initial and the final shape that is visible in Figure 3A. The color code from Figure 3 is used again so that parts of the membrane with broken linkers are colored red. Differences are predominant in the concave part of the initial shape, where the membrane has moved outwards and detached. The tension force in such concave parts points outwards and, together with the pressure, initiates a bleb without requiring any weakening of the cortex. This simulation thus supports the finding in Collier et al. [5].
In Figure 5, we compare the final shapes for different parameters of the linker strength λ l , more precisely, slices of the shapes in Figures 3A,B. Note that the color code is different (see caption of Figure 5). The deformation isn't much stronger as, once the membrane is detached, the linker term doesn't influence the evolution any further. But a weaker linker strength λ l and, thus, less resistance to breaking leads to a wider bleb site.
A smaller resting length parameter x 0 increases the surface tension, which leads to a faster evolution and a stronger final deformation. This is visible in Figure 6 where we compare the slices of the shapes in Figures 3A,C, and the (red) curve for the smaller x 0 indicates that the membrane has moved further away from the initial shape.
The impact of a higher pressure is illustrated in Figure 7 where slices through the shapes in Figure 3A (blue) and Figure 3D (red) are overlayed. The effect resembles a bit that of a smaller linker strength in that the deformation isn't much different and in that the bleb site is much bigger. The pressure term doesn't break down after detachment and continues to push outwards, though, so that the membrane has moved a bit further throughout the bleb site.  4 | Slices through the central axis of shapes obtained in simulations that we report on in section 5.3. A magnified image of the black box is presented on the right. The initial (dark gray) and the final shape (red/light blue) for the computation with the data from Table 1 are displayed. The color code is as in Figure 3. Strongly deformed parts of the membrane, where the linkers are broken, are red and predominant in the concave part of the initial shape.  Table 1 (blue), and one with a smaller linker strength parameter, namely λ l = 12 (red). In the latter case, the forces that keep the membrane attached to the cortex are smaller (see Equation (48), (49), λ l corresponds to k l in the dimensional model). The area where the membrane has detached from the cortex is bigger, so the blebbing region is wider.  Table 1 (blue), and one with x 0 = 0.85 (red). The latter effectively increases the membrane tension but also reduces the area at which the membrane is at rest (see Equations (50), (51) and in between). In the concave region of the initial surface the membrane thus moves farther and its area shrinks further than before.  Table 1 (blue), and one with a higher pressure parameter, namely p 0 = 30 (red). In the latter case, the cell membrane is further pushed outside, and the area where the membrane has detached from the cortex is wider. TABLE 2 | Non-dimensional parameters for the simulation with an initial surface obtained from image data, which we report on in section 5.4.

Application to Experimental Data
Apart from given, "in-vitro" geometries and their influence on blebbing, users may also be interested in studying the effect of "in-vivo" geometries that are obtained from experimental data. The image postprocessing outlined in Du et al. [53] enables users to extract triangulated surfaces representing the cell membrane from 3D images of cells, which then can be steered into the software framework. This was done with data of a Dictyostelium cell (also used in Du et al. [53]) moving by actin-driven pseudopods without any blebbing. However, the purpose is again to showcase the capability of the software framework rather than to extract any quantitative information, which is left for future investigations. We used the non-dimensional parameters in Table 2, T = 20, and τ = 0.02. Figure 8, left, shows the triangulated surface Ŵ 0 h that has been obtained from the image data. On the right of Figure 8 the final shape is displayed where the same color code as in Figure 3 for the deformation strength is used. As in the simulations before we observe that blebs form in concave regions. We also see some deformations at the sides where small protrusions become quite spiky. Both tension and resistance to bending are expected to prevent any singularities to occur, however the geometry seems under-resolved by the mesh in these areas.

Extensibility
Due to the use of the Python scripting language for the high level control of the simulation, it is fairly easy to extend our mathematical model to include additional effects. So for example we can conveniently add surface reaction-diffusion equations to model some biochemistry occurring on the membrane, i.e., adding a system of the general form for some reactants c : Ŵ 0 → R r , where the function R : R r → R r describes the reactions. The evolution of c can depend on the evolution of the membrane, i.e., both the diffusion tensor D and the source and sink term R may depend on the deformation u. The reactants c may just be passive in the sense that they do not influence the evolution of the membrane but, in general, they will enter the evolution equation for the surface deformation. For example, they may change some material parameters and, thus, the forces acting on the membrane. As a proof of concept, we investigate the effect of adding a chemical signal (scalar, c maps to R) to the model to contain the spread of the bleb. In the current model a bleb will, in general, not remain localized but the membrane will detach  Table 2. The color scheme is as in Figure 3. See section 5.4 for further details.
from large parts of the cortex. In experiments, bleb protrusion stops, forming a fairly spherical cap attached to the part where membrane and cortex still seem connected [2]. To reproduce this phenomenon at least qualitatively, we add a scalar signal c that is initially zero but is subsequently produced where the linkers are broken, and it then diffuses along the membrane. In regions where the linker molecules are still intact the signal leads to an increase of the linker strength. This can prevent further detachment of the cortex and effectively stop any further protrusion of the bleb. We can localize this effect further by making the diffusion tensor D in the reaction-diffusion equation for the signal c depend on the cortex detachment from the membrane. More precisely, we may assume a large diffusion coefficient in regions where |u − u c | > u B and a small diffusion coefficient elsewhere.
The variational form of the equation for the signal reads for test functions z with initial condition c(·, 0) = 0 and augments Problem 2.1. Here, D c > 0 is a diffusivity parameter and χ(u) = H(|u − u c | − u B ) (recall that H denotes the Heaviside function) is the characteristic function of the region where the cortex is detached. In the diffusion term, d f is a reduction factor. If d f = 1 then the diffusion is a uniform constant D c while if d f < 1 then the diffusion is reduced by this factor in regions where the cortex is still attached. The signal is produced at a rate r c > 0 as soon as the cortex detaches up to a maximum value of l c . We discretise (63) in space using the surface finite element method introduced in section 3. In time we treat the diffusion term implicitly and the reaction term explicitly, resulting in a linear problem in each time step. The model for the signal evolution can be conveniently added to the code with a couple of lines: = ( i n n e r ( c , z )+ t a u * s i g n a l T e r m s _ i m ) * dx == \ ( i n n e r ( s i g n a l _ n , z )− t a u * s i g n a l T e r m s _ e x ) * dx As stated above, we assume that the signal is not passive but influences the linker strength, i.e., λ l in (55) is replaced by a c dependent function. We want the linker strength to increase from the original value λ l to some value λ L > λ l when c increases from some value c b > 0 to some value c B > c b . For this purpose we define where ξ is a piecewise defined monotone C 1 weight function with ξ Note that we treat the impact of the signal on the surface evolution explicitly in time to simplify the presentation. Consequently, we can solve for the new membrane deformation and the new signal concentration in separate steps within the time loop. Details are available in the published code (see the Data Availability Statement). In summary, the signal evolution and its interaction with the surface evolution is characterized by seven additional constant parameters, D c , l c , r c , d f for the signal and λ L , c b , c B for the surface.
We test the effect of the signal on the evolution of blebs with an artificial surface geometry but with less symmetry than in section 5.3. Figure 9, left, gives an impression of Ŵ 0 . In particular, the surface is not rotationally symmetric. Nondimensional parameters that remain fixed for all simulations that we report on further below are given on the left of Table 3. They coincide with the parameters used in section 5.3 (see Table 1) except for a smaller x 0 (this increases the tension force) and a larger λ p (this increases the pressure difference). All the simulations reported on in this subsection were carried out on a triangulated surface with around 40 000 elements with a FIGURE 9 | Surface evolution with parameters given by Table 3 for the extended model in section 5.5 but with r c = 0 so that the signal c remains zero. The left figure shows the initial surface. The white curve indicates a plane along which we cut the shapes for visualization purposes in Figures 10,11. In the middle and on the right, both the initial surface and the deformed surface are displayed at the final time T = 5. The initial surface is fully displayed with the white line in the middle and cut on the right to illustrate the cutting procedure. The coloring of the original surface indicates the magnitude of the deformation as in Figure 3.  We first checked the case r c = 0, for which c = 0 solves (62), i.e., there is no signal. Figure 9, middle and right, displays the surface at the final time T = 5. We spot that the cortex has detached from the membrane everywhere. We will show now that this can be prevented from happening by accounting for a signal and a strengthening of the linker terms.
We next consider the case d f = 1 so that the diffusion parameter D c is constant. The additional parameters required for the signal and the coupling are given on the right of Table 3. The evolution of the surface with signal is shown in Figure 10. We see the signal being generated where the deformation is strong, i.e., where membrane moves away from the initial shape and thus from the cortex. However, at time t = 5 the membrane is still close to the initial shape in regions away from the concave parts. Comparing this with the simulation without the signal (see Figure 9, right) it is clearly visible that the increase of the linker strength due to the diffusing signal leads to a localization of the blebs. We also see in Figure 10 that the shape is not stationary at time t = 5 but the deformation continues to grow, which we now want to prevent.
We reduce the diffusion coefficient in regions where the linkers are not broken by setting d f = 0.01 whilst keeping the other parameters as in the previous simulation. The expectation is that the strengthening of the linkers is more focused on the boundary of the bleb than before and increases there faster as the signal is not transported away by the diffusion. Thus, also the linker strength should increase faster at the boundary of the bleb. The simulation result is demonstrated in Figure 11. Now, the simulation reaches a stationary state at time t ≈ 2.5 with a small bleb that has developed in the concave region at the top of the surface. At the beginning of the simulation, secondary protrusions also start to develop at the three saddle-shaped sides of the surface. Subsequently, the resulting production of the signal c leads to a strengthening of the linkers, which then pull the cortex back so that the protrusion disappears again.

CONCLUSION
A general modeling framework for the onset of blebbing has been presented and analyzed. It is formulated in terms of partial differential equations on the initial membrane, which is considered as a hypersurface. Various forces acting on the plasma membrane due to its elastic properties, linker molecules coupling it to the cell cortex, and pressure difference across the membrane are accounted for. Fluid flow within and outside of the cell is essentially neglected modulo a drag force but may be considered in future studies.
The general framework is particularly flexible with regards to membrane tension and the coupling forces. A convergence analysis of a surface finite element discretisation shows its robustness to model alterations within not too restrictive limits.
There are some open questions with regards to the discretisation in time, and as blebs are local events, spatial mesh adaptivity may be beneficial.
Software for a specific instance of the general model is provided and has been used to perform some numerical simulations. A convenient high-level interface in Python allows for directly implementing the model in its variational form and solving it by an efficient software back-end. Standard software usually does not provide functionality for numerically solving problems on moving domains or hypersurfaces in 3D out of the box but requires a substantial amount of coding. We hope that our approach will address this issue and simplify the  Table 3. The diffusion of the signal is a uniform constant. The color from blue via green and yellow to red corresponds to the signal strength: Dark red indicates c > c B , i.e., regions where the linkers (if not broken yet) are maximally strengthened by the signal c, dark blue indicates c < c b and thus regions where the linker strength is not altered by c.  Table 3 but with d f = 0.01, i.e., a reduction of the diffusivity of the signal to 1% of its original value in the region where the linkers are not broken. The color corresponds to the signal strength and is the same as in Figure 10. The shape at the final time T = 5 is a stationary state. The figure on the very right displays slices through the shapes at the final time T = 5 for two different grid resolutions and time steps. The yellow curve is the stationary state corresponding to the left three figures, i.e., using around 40 000 element and a time step of 0.001. The green curve corresponds to the final state on 160 000 elements and with a time step of 0.0005. As before, the white curve shows the initial surface. The final deformed surface with the coarser triangulation also serves to give an impression of the grid resolution compared to the size of the developed bleb.
implementation of such moving boundary problems. Future work on the software will include improvements of the efficiency of the implementation, by investigating, for example, the use of well-established techniques like local grid adaptivity and parallelization (which are available through the grid manager used for the simulations shown here [45]) as well as incorporating more sophisticated preconditioning methods in the linear solver implementation. This will be an important step toward extending the model to include, for example, fluid flow inside and outside of the cell.

DATA AVAILABILITY STATEMENT
The Python scripts used to obtain the results reported on here are available in a git repository hosted on the DUNE gitlab server: https://gitlab.dune-project.org/bjorn.stinner/sfem_blebs.
The main scripts are blebbing_artgeom.py, blebbing_imgdata.py, and blebbing_signal.py used for the results from sections 5.3, 5.4, and 5.5, respectively. The setup of the model, the time loop, and the solver used are contained in blebbing_compute.py. Some auxiliary functions can be found in blebbing_tools.py.
The experiments reported on in this paper, can be reproduced using the DUNE-FEM docker container. A script "startdune.sh" is available in the git repository to download and start the container. This requires the "docker" software to be available on the system. It can be downloaded for different platforms including Linux, MacOS, and the latest Windows version. More information is available under https://dune-project.org/sphinx/ content/sphinx/dune-fem/installation.html.

AUTHOR CONTRIBUTIONS
BS provided context and background, significantly contributed to the model and the numerical analysis, and set most of the paper. AD was core developer of the Python bindings and of the DUNE software framework and set parts of the paper. AN contributed to the model and the numerical analysis, performed the simulations, and set parts of the paper.

FUNDING
This project was supported by the Engineering and Physical Sciences Research Council (EPSRC, United Kingdom), grant numbers EP/K032208/1 and EP/H023364/1.