A PDE model for bleb formation and interaction with linker proteins

The aim of this paper is to further develop mathematical models for bleb formation in cells, including cell-membrane interactions with linker proteins. This leads to nonlinear reaction-diffusion equations on a surface coupled to fluid dynamics in the bulk. We provide a detailed mathematical analysis and investigate some singular limits of the model, connecting it to previous literature. Moreover, we provide numerical simulations in different scenarios, confirming that the model can reproduce experimental results on bleb initation.


Introduction
Bleb formation or "blebbing" is a biological process during which the cell membrane of an eucaryotic cell is disconnected from the cell cortex. The inflowing cytosol pushes out the free membrane part which builds a protrusion called a bleb. Due to regeneration mechanisms of the cell, the membrane connection to the cortex is restored and the protrusion is healed after some time. Despite these phases seem to be well established, the particular cause for the transition from one phase to the other, or the initialisation of the whole process are still subject to debate. Blebbing has been related to many interesting biological processes such as mitosis [7], cell spreading [8], and apoptosis [31]. It has also been noticed as migration mechanism, especially in embryonic [23], [12] and cancer cells. Nevertheless, according to [10], the lamellipodia migration mechanism had received a lot more interest and the authors found it necessary to emphasise the importance of deeper investigations into bleb formation. Following their promotion, effort has been made to derive bio-physical models and develop an understanding of this phenomenon by numerical means focusing on fluid-membrane interaction [36], [37], [39] on membrane dynamics including linker influence [27], or on linker kinetics [1]. There has also been effort in enhanced mechanical modelling [38] and on developing models of a full bleb life cycle in three dimensions [28].
In this work we propose a model for bleb formation in three dimensions taking into account the following aspects: • Elastic properties of the membrane described via Canham-Helfrich energies.
• Forces of the actin cortex exerted on the membrane via linker proteins.
• Activation and deactivation of linker proteins as well as possible movement of the proteins.
• Intracellular fluid-dynamics and the corresponding fluid-structure interaction.
A key issue of the cell blebbing phenomenon is the fact that there are actually two free boundaries, the membrane and the cortex that strongly interact via linker proteins. We describe the membrane by a height relative to the cell cortex, which is subject to forces exerted by a surrounding fluid as well as proteins that connect the cell membrane with the cell cortex and act like springs (cf. [1]). We take into account that these proteins may disconnect by introducing a ripping rate function whose steepness is controlled by a small parameter ϑ. This way we rediscover the model of [27] by ignoring movement of the linker proteins and passing to the limit ϑ 0. Besides the detailed mathematical modelling we provide a detailed analysis of the model in the case of small deformations of the membrane relative to the cortex, where a linearization of the mechanics applies. The key nonlinearities we focus on are hence due to the presence of the linker proteins. Besides well-posedness of the time-dependent model we prove the existence of stationary solution and show that some critical pressure (e.g. arising from cortex contraction) is necessary and sufficient to form a bleb, which we define as a deformation above a critical height of the membrane relative to the cortex at which linkers are ripping off. As mentioned above we study singular limits and show that the models of [1] and [27] arise as special cases respectively scaling limits of our model. Moreover, we provide a numerical study of the bleb initiation by a critical pressure.
The paper is organised as follows: After having fixed some basic notational conventions in Section 2, we derive a fourth order evolution system of equations in Section 3 starting from first principles. The main result of Section 4 is a global-in-time existence theorem for solutions of this system. The following Section 5 is devoted to proving existence of stationary solutions and studying stability of a particular subclass of stationary solutions by means of nonlinear semigroups. The analytical part of this paper ends with Section 6, in which we pass to the limit ϑ 0 in the parameter controlling the disconnection rate. Finally, we illustrate some properties of our model numerically by presenting simulations of different biological situations in Section 7 and conclude in Section 8.

Preliminaries
General By I X we denote the identity on the set X. We denote the minimum of two values x, y by x ∧ y. Let F : X → Y for Banach spaces X and Y ; the Gateaux derivative of F in direction v ∈ X at x ∈ X, if it exists, written as d d v (F ) x . The space of Radon measures or regular, countably additive measures on a measurable space Ω ⊆ R n that are absolutely continuous with respect to the Lebesgue measure is rca (Ω).
Differential operators and vectorial Sobolev spaces Gradients ∇u (in the weak and strong sense and independent of whether they are on open sets or manifolds) of scalar functions u are column vectors. The Jacobian ∇v of a vector-valued function v is the matrix whose lines are the transposed gradients of the components of v. We write J (v) = ∇v + (∇v) where M ⊆ R 3 is some two-dimensional manifold with outer unit normal field ν M and σ 2 the Hausdorff measure with Hausdorff dimension two. Shape derivatives For denoting the shape derivative of a functional W in directions θ, ϑ, we following the notation in [14] using d (W; θ) for the first and d 2 (W; θ, ϑ) for the second shape derivative. We denote by for a hypersurface M , δ ∈ I ⊆ R, I being an interval, and v : M → R 3 , the perturbation of M by v. In this particular case, the shape derivative of the functional W being defined on Differential geometry The notation of this paragraph follows [4] and [6]. With Diff q (X, Y ) we denote the set of all C q -diffeomorphisms mapping the Banach space X into the Banach space Y . The set of functions where J is a real interval, is denoted by Diff k,q (J × X, Y ). We call an n-dimensional manifold Γ ⊆ R n+k a real submanifold. The tangent space of a real submanifold at a point p ∈ Γ is denoted by T p Γ ⊆ R n+k . Let (Γ s ) s∈J be a family of real submanifolds in R n with a mapping X : J ×Γ → s∈J Γ s such that X s is a global parametrisation of Γ s on a reference manifold Γ. We call the set G (Γ s ; s ∈ J) = s∈J {s}×Γ s an evolving manifold. If Γ s are hypersurfaces, we use the term evolving hypersurface. For functions f s : Γ s → N , where N is a set, we define the function An evolving manifold G (Γ s ; s ∈ J) is smooth if T p G (Γ s ; s ∈ J) = {0} × R n for all p ∈ G (Γ s ; s ∈ J). For such a smooth evolving manifold, we denote the velocities ∂ s (X s )| r • (X r ) −1 by V r . If the manifolds are orientable, i. e., there exist smooth outer unit normal fields ν s : Γ s → R 3 , the velocities can be decomposed in their normal V s ν = V s ν ν s , V s ν = V s ν · ν s , and tangential V r τ components. We define a differential operator ∂ • s (f )| r = ∂ s (f •(s, θ) → (s, X s (θ)))| r • (X r ) −1 called the material derivative (of f with respect to X).-The parametrisation mapping shall always be given by the context if not stated explicitly. For any real submanifold Γ ⊆ R n , the tangential gradient ∇ Γ of a function f : Γ → R at p ∈ Γ is given by ∇ Γ (f )| p = P TpΓ ∇ f p , wheref is any differentiable extension of f to an open neighbourhood of p in R n and P TpΓ is the orthogonal projection onto the tangent space T p Γ of Γ at p. The projection matrix I − ν Γ (p) ⊗ ν Γ (p) and the projection P TpΓ are identified. This way, we also acquire the partial tangential derivatives For a differentiable, orientable real submanifold Γ ⊆ R n with normal field ν, we denote the Weingarten map by H Γ : Γ → R (n,n) , p → ∇ Γ (ν)| p . Then, the mean curvature of Γ is H Γ = tr H Γ = ∇ Γ ·ν and the Gaussian curvature is K Γ = det H Γ .
Physical dimension and units The sets R n , n ∈ N, are identified with the product sets R n ×D×U, where D is the set of all physical dimensions D = {T, L, M, N, . . . } (meaning time, length, mass, etc.) and U the set of all physical units. For x ∈ R n , we denote by x its second component (called the physical dimension of x). When we write x, we always refer to the first component.

Modelling
As mentioned in the introduction, many details of the process of bleb formation are still subject to research and not fully understood. Therefore, we aim at a rather abstract model following the general description of the process in [10]: The main parts of an eucaryotic cell that are involed in bleb formation are the cell membrane, which is basically a bilayer of lipid molecules, the cell cortex, which is a network of actin fibres, and elastic proteins which connect the cell cortex to the cell membrane. These linker proteins are only stretchable to a certain length above which they disconnect from the membrane. Inside the cell there is a fluid that is called the cytosol and the cell is itself swimming in an extracellular fluid. Caused by mechanisms which have not completely been understood yet, a certain patch of the cell cortex contracts and raises the pressure on the membrane locally. This way, the corresponding membrane patch is pushed so far away from the cortex that most of the linker proteins disconnect. The cytosol that pushes against the free membrane patch now causes the formation of a protrusion which is called a bleb. Over time, the protein linkers are reconnected to the cell membrane causing the membrane patch to be fixed to the cortex again and the bleb to vanish.

The fluid system
Let D ⊆ R 3 be a bounded, connected, and open set with sufficiently regular boundary. We require this set to be partitioned into the open connected set Ω ext 0 , modelling a reference region exterior to the cell, the open connected set Ω int 0 , modelling a reference region interior to the cell, and the boundary M 0 of Ω int 0 , which shall be a two-dimensional orientable C 2 -manifold, modelling the membrane in its initial state; its unit normal field is denoted by ν M0 : M 0 → R 3 . The region which is occupied by Consequently, the exterior region is Ω ext (t) = Φ (t, Ω ext 0 ) and the whole domain is denoted by Ω(t) = Ω ext (t) ∪ Ω int (t) and M(t) = Φ (t, M 0 ). Furthermore, the cell cortex is denoted by C ⊆ Ω int 0 and is modelled as a sphere which is fixed in time. Figure 1 illustrates a typical geometry compatible with the previous description.
Both the fluid in the inner region, representing the cytosol, and in the outer region wtih pressures p i : Ω i T → R and velocities u i : Ω i T → R 3 , i ∈ {int, ext}, are described by incompressible stationary Stokes equations in Eulerian coordinates on Ω i T = t∈[0,T ] {t}×Ω i (t) with final time T > 0 (time dependency will be brought into the system by boundary conditions). This is justified by small length scales as in [36]. For any functions g i : Ω i T → X, i ∈ {int, ext}, into a vector space X, we associate the function We further pose a homogeneous Dirichlet boundary condition at the exterior boundary which does not change over time, , be a force we will specify below. The stress at the interface M(t) is subject to a Neumann-type boundary condition: where T = µJ (u) − pI is the Cauchy stress tensor of the interior and exterior fluid and we define To assure well-posedness of the problem, we further require Taking an energetic point of view, we consider a variational formulation of the Stokes equations (1) with boundary conditions (2) (the no jump condition is encoded in the solution space): for all ϕ ∈ H 1 0 D, R 3 and q ∈ L 2 (D).

Remark 1.
Another formulation of the Stokes equations includes the term (∇u, ∇ϕ) L 2 (D,R (3,3) ) instead of 1 2 (J (u) , J (ϕ)) L 2 (D,R (3,3) ) . We observe that, because of u being solenoidal, so both expressions are equal. The motivation to use the latter is related to the structure of the employed Neumann boundary conditions.
We introduce a function h : [0, T ] × C → R, h = L, which is intended to give the membrane's height relative to the cortex in normal direction ν C , i. e., We further require M(t) lying in a sufficiently small tubular neighbourhood of C (this approach is analogous to [17]): exists.
Yet there is no guarantee that the mapping in (4) is bijective between C and M(t) as multiple points in M(t) may have the same projection point rendering the existence of a function h impossible. Therefore, we also pose an invertibility condition for the parametrisation of the membrane over the cortex: for y ∈ M 0 , x = P C (y).
Potential energy In view of Condition 1, we define a rescaled height h = δ h and we may regard M(t) as small perturbation of C by hν C with order of magnitude δ. In particular, we write In order to derive a mathematical model, let us turn to the bio-physical properties of the membranecortex system just described: The membrane shall consist of lipid molecules arranged in two layers. On the cortex, there are proteins connected to the membrane, therefore called linkers proteins. They are considered stretchable and shall obey Hooke's law for springs, so we can assign the potential energy density functional is a function playing the role of a spring constant in every spatial point and ρ a ∈ L 2 (C), ρ a = L −2 , is the density of active linkers (further explanation below, see Section 3.4). There are several models (cf. [35]) for the surface energy of membranes. A widely used example is the Helfrich energy (cf. [22], [40]): where H M : M → R, H 0 M ∈ R, and K M : M → R are the mean curvature, spontaneous mean curvature, and Gaussian curvature of M, respectively, with bending rigidity κ and Gaussian bending rigidity κ G .
The second order expansion of the total energy density of the membrane-cortex system at time where we have used that the derivatives of W C δ h(t, ·)ν C with respect to δ at 0 are equal to the shape derivatives of W in directionĥ(t, ·)ν C at zero. If the cortex was an equilibrium shape of Helfrich's energy, we would have d W; hν C = 0. But this may be especially not true if the cortex is contracted due to myosin motor activity. Therefore, we model d W; hν C = δp 0 , h for p 0 ∈ L ∞ ([0, T ], L ∞ (C)) interpreting p 0 as a stress that is exerted on the membrane due to the cortex contraction and transmitted by the fluid. The energy functional up to second order therefore is introduces a mechanism by which an initially flat membrane in a resting fluid may be deformed after all: Not considering d W; hν C as a parameter, but instead taking the terms that come out of a computation of this shape derivative would only change the coefficients of the ∆ C h and h terms in the variational principle we will derive below. The resulting equation is homogeneous and therefore does not show any deforming behaviour in case the membrane is initially flat and the fluid velocity zero. The more physical but also rather complex approach for introducing this mechanism would be to relate the pressure p 0 to shape deformations of the cortex and then describe the influence of p 0 on the fluid introducing another surface-bulk coupling this way.

Connecting the fluid and the membrane model
The fluid system is not closed, but subject to external forces f t . This is exactly where the membranecortex system comes into play: The potential energy of this system is considered to be the source of forces acting on the fluid and therefore being transformed into kinetic energy of the fluid. We also take a damping effect due to friction between the fluid particles and the cortex into account with a linear friction model with friction constant c having dimension c = MT −1 L −2 . In order to enforce Condition 2, these forces are all directed normally to the cortex, so the tangential part of f t is zero: Recalling (5), we have the fluid particles at the membrane moving in the direction of the cortex normal.
Taking the length of the velocity vector to be the change of the membrane's height in time, we have specified the Dirichlet boundary of u at M(t), and we therefore may express where [ϕ] X = ϕ • X(t, ·) −1 for X(t, ·) := I C + h(t, ·)ν C and any function ϕ with domain C and is the Dirichlet-to-Neumann operator of the Stokes problem Problem 1. A definition and references to important properties of this operator is given in Appendix A. By combination of both descriptions of the Neumann data, we are going to derive a PDE model:

PDE description of the height function
Approximation of the Dirichlet-to-Neumann operator For sufficiently regular Stokes flow velocity, we can make use of the small height condition Condition 1 to approximate the time-dependent Dirichlet-to-Neumann operator (see Appendix C) with its stationary version on C. This way, we arrive at a gradient-flow structure with L = cI H 1 (C,R 3 ) + DN 0 and ϕ ∈ H 2 (C).
Calculating the variation of the potential energy To derive a full PDE description for h, we have to calculate the variation of the potential energy functional. So first, we calculate the first and second shape derivatives of W. Recall, where H δ is the mean curvature of C δ ĥ ν C Observe that the integral over the Gaussian curvature is constant in δ due to the Gauss-Bonnet theorem (C δ ĥ ν C0 is homeomorphic to C) and therefore vanishes when differentiated in δ. The neccessary caclulations have been carried out before for the Willmore energy, e. g. in [17, p. 7]: The additional calculations use the same techniques and the interested reader may consult the appendix (Corollary 3) and Corollary 5) for details: Spherical cortex shape Significant simplification of these terms is achieved by considering C to be a sphere with radius R: (cf. [17]) and (cf. Corollary 4) and Force density equation All together, we arrive at the following expression for (7): where Remark 3. The form a (·, ·) may not be coercive on H 2 (C) nor may it be non-negative. In order to assure at least non-negativity, we make the following considerations: In case H 0 ≥ 0, we follow [17] and derive a Poincaré-type inequality from Courant's min-max principleˆC where 2 R 2 is the second eigenvalue of the Laplace-Beltrami operator, on span {1} ⊥ , i. e. for all functions with mean zero. As H 2 (C) = span{1} ⊕ span{1} ⊥ , for every u ∈ H 2 (C) there is a constant m (u's mean value) and u 0 (being mean-value-free) such that u = m + u 0 . We observe In case H 0 < 0, we need a compatibility condition on H 0 and R. We require 2

Protein linkers
The quantity ρ a has been mentioned before in modelling the potential energy of the membrane-cortex system. It models the density of linkers that are connected to the membrane. We also take linkers into account that are disconnected and whose density is denoted ρ i . Both active and inactive linkers are considered to be mobile species diffusing on the cortex. Moreover, they are transformed into each other as result of overstretching above a critical height h * ∈ C (C, [0, ∞)), which causes active linkers to disconnect, or regeneration mechanisms connecting inactive linkers to the membrane again. A reaction-diffusion-kind of system may be used to model these processes: where η a , η i ∈ [0, ∞) are the active and inactive linker diffusivities, k ∈ [0, ∞) a regeneration rate, and r : R → [0, ∞) a disconnection rate being Lipschitz continuous and r (−∞,0) = 0 (a typical example is the non-negative part). The disconnection of linkers from the membrane is considered to be a fast process. To account for this, a (small) parameter ϑ ∈ (0, ∞) is used for rescaling the argument of r.
(In Section 6, we analyse the solution's behaviour when ϑ 0.) For the sake of readability, we may use the abbreviation in the following.
Remark 4. Approaching the linker movement by a reaction-diffusion model is motivated by the work of [1]. They consider the following equation for the membrane height h and linker density ρ a (we adopt the notation of this work for their parameters and quanitites): with a maximal linker density ρ 0 and a disconnection rate k off . But there is an important new aspect to the model presented here: The concept of inactive linkers is not present in (10), but a gauge protein density ρ 0 is assumed of which a part is connected ρ a and ρ 0 − ρ a is disconnected. As consequence of this condition, their approach is limited to scenarios where the cortex is intact. Nevertheless, it has also been observed [10] that bleb formation may be triggered by cortex disruption (leading to a hole in the cortex). This case is contained in our active-inactive linker setting with ρ a (0, where D is the area of the hole in the cortex (cf. Section 7.3). Indeed, (10) is a specialisation of our model: Set η a = η i = η and ρ a (0, ·) + ρ i (0, ·) ≡ ρ 0 . Adding (9a) and (9b), we get ∂ t (ρ a + ρ i ) + η∆ C (ρ a + ρ i ) = 0, which is solved by ρ a + ρ i ≡ ρ 0 . This way, we can express ρ i = ρ 0 − ρ a giving (10).
After having derived a PDE model for the blebbing phenomenon, we will deal analytically with the following issues in the next sections: 1. global-in-time existence of weak solutions (Section 4), 2. existence of stationary solutions and their stability (Section 5), 3. convergence of stationary solutions to a singular limit when ϑ 0, 4. and rediscovering the model for bleb formation proposed in [27] (Section 6).

Time-dependent solutions
In the following three chapters, we analyse a variatonal formulation of (8), (9a), and (9b) having the following strong equivalent for sufficiently regular h, ρ a , ρ i : where we write r ϑ (h) = r h−h * ϑ . Let us summarise the properties of the parameters involved: • κ and λ are non-negative constants, whereas γ is also a constant but not necessarily non-negative.
• The operator L in front of the time-derivative is the sum of the identity and the Dirichlet-to-Neumann operator of the Stokes problem.
• The function ξ is in L ∞ ([0, T ], L ∞ (C)) as well as the pressure p 0 . Also, ξ is assumed to be non-negative a. e.
• The repairing rate k is a non-negative constant.
• The diffusivities η a , η i are taken to be positive.
• The disconnection rate r is assumed to be non-negative and Lipschitz; the corresponding steepness parameter ϑ shall be non-negative as well.
For better readability, we introduce the following forms: The antisymmetric structure of the linker equations allows for a simple conclusion: Lemma 1 (Mass Conservation). Let ρ a , ρ i be parts of a solution to Problem 2 in the strong sense.
Proof. Add (11b) and (11c), integrated over C, to achievê Then, with the Divergence Theorem on closed manifolds, we get The integral and the weak differential operator commute, so´C ρ a +ρ i dx is constant almost everywhere in time. Considering the non-negativity of the initial values, the claim follows.

Global-in-time existence
In the following, we refer to • the coefficients κ, γ, λ, the pressure p 0 , the critical height h * , the linkers spring constant ξ, the function r with its Lipschitz constant L r , the parameter ϑ, the repairing rate k, • the positive definite, self-adjoint operator L = S 2 , and the constant Ξ > 0 such that Ξ u 2 as the data of Problem 2.
We set and define the operator and Remark 5 (Well-defined). The unique existence of a solution to (14) follows by standard parabolic PDE theory. Existence and uniqueness of solutions to (13) may be shown by employing a Petrov-Galerkin-type approximation argument. (We refer the interested reader to Appendix B.) Hence, F T is well-defined.
We aim at a Banach-type fixed point argument. To this end, the following a priori estimates are derived, which will give us Lipschitz continuity of the map F T . We will then show that F T is contractive in a rescaled topology of X T and Y T by introducing the rescaled L ∞ norms for a constant β > 0 that we may choose with respect to the problem data and T . In the following, X T and Y T are equipped with their corresponding rescaled norms. We will not state the rescaling constant explicitely since, from now on, the rescaling factor of every rescaled norm we deal with shall be e −βt if not stated otherwise.
Remark 6. We will make use of the following interpolation embedding, which is a specialisation of [3], Theorem 3.1 for θ ∈ [0, 1] and m ∈ [0, ∞): isometrically. Choosing θ = 1 4 and m = 1, we have A priori estimates There will be a lot of constants in the following estimates, so for ease of notation and in the sake of readability, we will slightly abuse notation and write C(ξ, Ξ, . . . ) or similar for an expression that only depends on the problem data; it does not necessarily denote the same expression in every occurance. For convenience, we abbreviate = ffl C h dx. Lemma 2. Letρ a ∈ Y T ,ρ a ≥ 0 a. e. and T > 0. The following bound holds for h = H T (ρ a ): Proof. (i) We observe that because the mean value of h is constant in time, we have = 0 and where we dropped the term − (ρ a ξh, h) L 2 (C) because of its non-positivity. We note that With this simple observation, the remaining right hand side term is bounded by employing the Cauchy-Schwartz and then the Young inequality: The Grönwall inequality then gives us the bound: With the positive definiteness of L, we even have This finishes the proof.
The following a priori bounds hold for (ρ a , ρ i ) = G T (h): .
Proof. (i) We test (14a) by σ a = ρ a : and (14b) by σ i = ρ i : and add the equations leaving out the non-positive terms on the right hand sides (r ϑ ≥ 0 by assumption): The Grönwall inequality now implies thus giving the claimed L ∞ − L 2 bound.
(ii) We start again with (19), drop η i ∇ C ρ i 2 L 2 (C) on the left, integrate in time, and then also drop Together with the previous estimate, we obtain the claimed bound.
(iii) We test (14a) by σ a ∈ H 1 (C) with σ a H 1 (C) ≤ 1 such that ρ a , σ a H −1 (C) ≥ 0 w. l. o. g., shift the gradient term to the right, and use the Hölder inequality on the right hand side terms: Squaring the inequality and integrating in time, we obtain Lemma 4. Given a time T > 0, there exist constants M 1 > 0, M 2 > 0, depending on the data of Problem 2 such that the operator F T maps the set We choose M 1 as the expression on the right hand side of the equation in Lemma 2. The bound M 2 then directly follows with Lemma 3 by setting We may consider (14a), (14b) as reaction-diffusion system with right hand side According to [29], Lemma 1.1, quasi-positivity of f is sufficient to guarantee preservation of nonnegativity. We see immediately that for r 1 , r 2 ≥ 0, f 1 (0, r 2 ) = kr 2 ≥ 0 and f 2 (r 1 , 0) = r ϑ (h)r 1 ≥ 0, so ρ a and ρ i are non-negative since the initial data are non-negative. Therefore, To apply Banach's fixed point theorem, we have to show that there is a β > 0 such that F T is a contraction for an arbitrarily chosen T > 0. To this purpose, we choose an arbitrary pair of arguments where L ∈ [0, 1). For ease of notation, we abbreviate There is a β (depending only on problem data) such that the following estimate holds: for some L h ∈ [0, 1).
Proof. Subtracting (13) for the two different arguments one achieves (using the lineariy of L): (i) We then test by ϕ = h ∆ and insert a suitable zero expression on the right hand side. (Note that h ∆ is mean-value-free as The last term on the right hand side is non-positive, so we drop it. We apply the Hlder and then Young inequality, and use the positive definiteness of L: With the Grönwall inequality, we obtain By multiplying the inequality with e −βt and inserting e −βs e βs under the right hand side time integral, we arrive at: so by choosing β appropriately large, the claimed contraction estimate follows.
Proof. (i) Subtracting (14b) for the two different arguments and inserting a suitable zero expression on the right hand side, we obtain By choosing σ i = ρ i∆ , we find We apply the Hölder inequality on the right hand side and then use Young's inequality with a parameter ε so small that ε ρ i∆ 2 L 4 (C) can be absorbed by the H 1 norm on the left: where α is some positive constant. Next, we subtract (14a) for the two different arguments and insert a suitable zero expression on the right hand side to obtain Testing with σ a = ρ a∆ leads to Applying the Hölder and then the Young inequality on the right hand side terms, we obtain We add (25) and (24) leaving out all non-negative terms on the left hand sides: The Grönwall inequality now implies We multiply the inequality by e −βt and insert e −βs e βs under the right hand side time integral: Due to the a priori bounds on ρ 1 a and the interpolation mentioned in Remark 6, we know that ˆt is also bounded. Hence, choosing β large enough (depending only on problem data and T ), we obtain the claimed contractive estimate.
Theorem 1. (i) Given any T > 0, Problem 2 has a unique weak solution in K T .
(ii) Therefore, given the fact that the bounds in K T include only problem data, Problem 2 is wellposed in terms of the definition of Hadamard.
Proof. Let T > 0 arbitrarily chosen and define K T as above. It is clear that The fixed point operator F T maps K T into itself according to Lemma 4. Furthermore, F T is a contraction due to Lemma 5 and Lemma 6, so the Banach fixed point theorem applies on F T and guarantees the existence of (h, ρ a , ρ i ) such that F T (h, ρ a , ρ i ) = (h, ρ a , ρ i ) what immediately implies that (h, ρ a , ρ i ) is a weak solution of Problem 2. Banach's fixed point theorem also guarantees uniqueness of such a fixed point, so the weak solution is unique.

Stationary solutions
In this section, we deal with solutions (h, ρ a , ρ i ) of Problem 2 with ∂ t h a. e.
= 0, which we call stationary solutions. From now on, we only consider the case where a (·, ·) is coercive (cf. Remark 3) and ξ, h * , and p 0 shall be time-independent. As pointed out in Remark 3, the coercivity requirement puts restrictions on the spontaneous mean curvature, i. e., H 0 ∈ (−∞, − 4 R ) ∪ (0, ∞). Physically speaking, this means that we eiher consider a membrane whose natural tendency is to form a concave shape (negative curvature) with a curvature of an absolute value of at least 4 R , or a strictly convex shape (positive curvature). Apart from this, all assumptions on the parameters are the same as in the previous section. The associated stationary problem reads: for all ϕ ∈ H 2 (C) , σ a , σ i ∈ H 1 (C).

Basic properties
Lemma 7. Let h and ρ a be parts of a solution to Problem 3. If ρ a ≥ 0 a. e., h H 2 (C) is bounded by a constant depending only on κ, γ, λ, p 0 , and the domain C.
Proof. Testing (26a) with h, we obtain Due to the coercivity assumption on a (·, ·) we stated at the beginning of the section, we further have for some α > 0. Choosing ε small enough (depending on C and κ, γ, and λ) such that ε h 2 L 2 (C) may be absorbed on the left hand side, we derivẽ withα > 0. Since α also depends only on C, κ, γ, and λ, the claim now directly follows.
The special structure of (26b) and (26c) also gives: Lemma 8. Let ρ a , ρ i ∈ H 1 (C 0 ) be parts of a solution to Problem 3. Then Proof. Testing (26b) and (26c) with the same σ ∈ C ∞ c (C) and adding both, we get On closed Riemannian manifolds, all weak solutions of this problem are smooth and only differ up to a constant. Therefore η a ρ a + η i ρ i = ρ 0 a. e.
With the previous results, we are in the position to state the following observation: Proof. Choose an arbitrary function p 0 ∈ L 2 (C, [0, ∞)) and consider the problem of finding h s ∈ H 2 (C), for s ∈ [0, ∞), such that for all ϕ ∈ H 2 (C), where ρ 0 is the constant of Lemma 8. We note, h s a. e. = sh 1 , which directly implies Assume h 1 ≤ 0. Testing (27) with h 1 , we find h 1 H 2 (C) ≤ 0, so h 1 a. e. = 0, which is no solution to (27).
Since h s is continuous, there exists x * ∈ C such that h s (x * ) = sup x∈C h s (x). Choose s * large enough such that h s * (x * ) > h * ; then, we have a ball B δ (x * ) ⊆ C (in the induced subtopology of C), δ > 0, where h s * (x) > h * , x ∈ B δ (x * ), which is a set of non-zero two-dimensional Hausdorff measure. Now assume h being part of a stationary solution to Problem 3 with pressure s * p 0 and being lower or equal h * almost everywhere. Then (26c) becomes = ρ 0 for some ρ 0 ≥ 0 (cf. Lemma 8). Therefore, h fulfils (27) for s * . But this contradicts the observation h B δ (x * ) > h * made above and the claim must be true.
Elimination and reconstruction of the inactive linker density Motivated by Lemma 8, we introduce two auxiliary variational problems, which are parametrised by m 0 (recall Lemma 1) or ρ 0 , respectively, such that for every stationary solutions of Problem 2 there is an auxiliary problem that is fulfilled by it and whose solutions allow for construction of a solution of Problem 3. In case η a ≥ η i , we will employ Problem 4 (Auxiliary problem). Let m 0 ∈ [0, ∞). Find h ∈ H 2 (C) and ρ a ∈ H 1 (C) such that for all ϕ ∈ H 2 (C) and σ ∈ H 1 (C).
In case η i > η a , we will use Problem 5 (Auxiliary problem). Let ρ 0 ∈ R. Find h ∈ H 2 (C) and ρ a ∈ H 1 (C) such that for all ϕ ∈ H 2 (C) and σ ∈ H 1 (C). (ii) In case η i ≤ η a , all solutions (h, ρ a ) of Problem 4 (with parameter m 0 ) can be extended to a solution (h, ρ a , ρ i ) of Problem 3 such that η a ρ a + η i ρ i ≡ const and´C ρ a + ρ i dx = m 0 . If ρ a ≥ 0 a. e., ρ i ≥ 0 a. e.
(iii) In case η i > η a , all solutions (h, ρ a ) of Problem 5 (with parameter ρ 0 ) can be extended to a solution (h, ρ a , ρ i ) of Problem 3 such that η a ρ a + η i ρ i ≡ ρ 0 If ρ a ≥ 0 a. e., ρ i ≥ 0 a. e.
Proof. (i) Let h, ρ a , and ρ i be stationary solutions of Problem 2. According to Lemma 1, it holds Additionally, due to Lemma 8, we have Inserting (30) into the integrated version of (31) gives an equation for ρ 0 : Inserting this expression and ρ i = 1 ηi (ρ 0 − η a ρ a ) into (26b), we get (ii) Now let (h, ρ a ) be a solution of Problem 4. Choose We know and ρ a satisfies (26b). We further calculate: and ρ i fulfills (26c). A small computation showŝ Furthermore, if ρ a ≥ 0 a. e., a standard maximum principle guarantees non-negativity of ρ i . (iii) The reconstruction of ρ i is just the same as in (ii) with ρ 0 being directly given.

Remark 7.
In the case η i > η a , the total linkers' mass is given by It is not hard to see that (the latter inequality being due to´C ρ a dx ≤ 1 ηa´C ρ 0 dx) implying that there are stationary solutions with arbitrary small (but non-negative) or arbitrary large mass. We conjecture existence of stationary solutions for all non-negative masses m 0 . However, it is not clear how a surjective map from ρ 0 to m 0 can be defined-not even a continuous map, so the mean value theorem is not directly applicable-; hence, a rigorous argument is still missing. Fix m 0 ≥ 0 and let

Fixed point argument
where C is a constant depending on the constant in Lemma 7 and the embedding constant of H 2 (C) cont → L ∞ (C). K is clearly convex and bounded. K is also closed: Let (h n , ρ n a ) n∈N be a sequence in K that converges in L ∞ (C) × L 1 (C). Hence, for all n ∈ N, C ρ n a dx ≤ m 0 , and, as L 1 (C)-convergence implies convergence of the integrals´C ρ n a dx n→∞ −→´C ρ a dx, the inequality it preserved in the limit. Furthermore, for all x ∈ C and all ε > 0, ρ a dξ ≥ 0 which finally gives (using the Lebesgue point property) ρ a ≥ 0 a. e.
We now define for all ϕ ∈ H 2 (C) and all σ ∈ H 1 (C). It holds F (K) ⊆ K: In order to show non-negativity of ρ a , we test (32b) with (ρ a ) − : By assumption η a ≥ η i and the right hand side is always ≤ 0, which implies (ρ a ) − = 0 a. e. The appropriate a priori bound of ρ a is obtained by testing (32b) with 1 and leaving out the gradient term on the left hand side: Dividing both sides by k ηa ηi leads to the claimed bound. The bound of h follows directly with the non-negativity ofρ a and Lemma 7.
F is continuous: Let hn ,ρ a n n∈N be a sequence that converges in K. Take an arbitrary subsequence F hnk ,ρ a n k = (h n k , ρ n k a ). We have the a priori bound immediately by Lemma 7. We also observe We can choose ε small enough such that ε |C| < kηa ηi , so ε ffl C (ρ n k a ) 2 dx can be absorbed on the left hand side of (32b), which gives a uniform bound of ρ n k a H 1 (C) . So by weak compactness, there are subsequences h n k , ρ n k a ,h n k ,ρ a n k such that the integrals in (32) converge and we have for the limits Due to unique solvability of (32),h andρ a are the same limits for all subsequencesh n k l ,ρ a n k l , respectively, so h n , ρ n a converge to F (h,ρ a ). But this implies that F is continuous. F (K) is relatively compact: Choose a sequence (h n , ρ n a ) = F hn ,ρ a n . From the previous calculations, we have a bound of ρ n a H 1 (C) and know that h n H 2 (C) is bounded a priori. The compact impeddings H 1 (C) comp → L 1 (C) and H 2 (C) comp → L ∞ (C) directly give us a convergent subsequence (h n k , ρ n k a ) in L ∞ (C) × L 1 (C). With these properites of K and F , Schauder's fixed point theorem applies, so F has a fixed point in K. This fixed point (h, ρ a ) is a solution of Problem 4, which then may be extended to a solution (h, ρ a , ρ i ) of Problem 3 due to Lemma 10.

Local exponential stability of non-critical stationary solutions
In this section, we will be concerned with the local stability of stationary solutions for a specific ripping interpolation function r ϑ (h) = h−h * ϑ + for solutions below the critical height h * . We are going to make a linearised stability argument using nonlinear semigroup theory (for an introduction see, e. g., [5]). Note that, due to Lemma 7, there are stationary solutions h < h * a. e. for sufficiently small p 0 . Without loss of generality, we set ξ a. e.

= 1. Observe that the operator
is invertible as it is linear and its kernel is zero-dimensional (to see this, one may test with x and use the positive-definiteness of L).

Remark 1.
Well-posedness of this problem follows with sufficient regularity of initial data and the inhomogeneity almost by standard techniques as presented, e.g. in [19,Chapter 7]: Using a Galerkin approach, we take existing solutions h, ρ a , ρ i of Problem 2 and express them with eigenfunctions (w k ) k∈N of ∆ C constituting a Schauder basis of L 2 (C). We obtain the projections h m (t, for ϕ ∈ span{w 1 , . . . , w m }, where π m is the projection into span{w 1 , . . . , w m }. We differentiate in time and test by ∂ t h m . Integrating in time then gives an energy estimate of the L 2 [0, T ], H 2 (C) norm of ∂ t h m (the nonlinearity on the right hand side is uniformly controlled due to the already established a priori bounds of the solutions). L 2 (C) regularity for ∂ t ρ a , ∂ t ρ i , is established by testing the projected equations with ∂ t ρ m a , ∂ t ρ m i , respectively. Control of the nonlinearities also helps with establishing increased spatial regularity.
This result makes us confident that existence theory of (11) might as well be developed using nonlinear semigroups.
Lemma 11. A stationary solution to Problem 6 is locally exponentially stable iff it is an exponentially stable stationary solution to the corresponding linearised system (see below).
Proof. Due to [15], Theorem 2.1, the claim follows if we can show differentiability in H 2 (C) × L 2 (C) 2 of the nonlinear C 0 semigroup generated by A + F . Applying Theorem 3.3 in [25], it is sufficient to show Fréchet differentiability of F in H 2 (C) × L 2 (C) × L 2 (C) on a sufficiently small ball around a stationary solution h, ρ a , ρ i with h < h * a. e. and the Lipschitz continuity of the derivative therein. Indeed, its Fréchet derivative is for (h, ρ a , ρ i ) ∈ B r h, ρ a , ρ i , r > 0, where r is chosen sufficiently small such that h < h * a. e. (This is Note, as L −1 is linear and bounded, we may drop it in the following calculations without loss of generality. Consider for (h ∆ , ρ a∆ , ρ i∆ ) the difference quotient After straightforward simplifications, we obtain (Since h < h * a. e., r ϑ (h) vanishes a. e.) Again, we use the imbedding H 2 (C) cont → L ∞ (C) to conclude that if h ∆ is small enough, h + h ∆ < h * a. e., so the last term vanishes and all together the difference quotient goes to zero as (h ∆ , ρ a∆ , ρ i∆ ) For a stationary solution h, ρ a , ρ i of Problem 6 with ρ a , ρ i ≥ 0 a. e. and with mass´C ρ a +ρ i dx = m 0 andh < h * a. e., the linearised system of Problem 6 is declared as such that ρ 0 a , ρ 0 i ≥ 0.
We start by showing that the inactive linkers decay exponentially: Lemma 12. Let ρ i∆ ∈ H 2 (C) be part of a triple solving Problem 7. Then, for ω ∈ (0, ∞).
Despite equations (36b) and (36c) looking very symmetric, their decaying behaviour is not. To show exponential decay for the acctive linkers, we require an additional condition: The initial values (37) to be chosen such that´C ρ 0 a + ρ 0 i dx = m 0 , so´C ρ a∆ + ρ i∆ dx = 0, which we call the mass conservation property. where ω > 0 is the same as in Lemma 12 and D ∈ (0, ∞).
For ρ a∆ (t, ·) L 2 (C) , we have at least a time-uniform bound: Lemma 14. Let ρ a∆ be part of a triple solving Problem 7. Then ρ a∆ is bounded in · L 2 (C) uniformly for all t ∈ [0, ∞).
Proof. The operator B : H 2 (C) → L 2 (C) , ρ a∆ → −η a ∆ C ρ a∆ is self-adjoint and it's resolvent R λ (B) is bounded by 1 λ for λ > 0. Therefore, σ (B) ⊆ (−∞, 0] and so it generates a strongly continuous (even analytic) semigroup bounded by 1. As ρ a∆ solves (36b), it has representation as where T is the semigroup generated by B. Consequently, Combining these lemmas, we even find exponential decay for the active linkers: Lemma 15. Let ρ a∆ be part of a triple solving Problem 7. Then for α a , E a ∈ (0, ∞).
Proof. First, we observe that for stationary solutions with h < h * a. e., there are no inactive linkers, i. e., ρ i a. e. = 0. (This can can be seen in the last paragraph of the proof to Lemma 9: Since r ϑ (h) = 0, the right hand side of the inactive linkers equation is zero and therefore testing with ρ i shows that ρ i 2 H 1 (C) = 0.) This has the critical implication that ρ a is constant in space (see Lemma 8). Second, recall that the operator L is positive and self-adjoint, so we have positive square root S 2 = L. We apply L on both sides of (36a) and test with h ∆ . Then, we make use of Young's inequality and leave out the terms of the bilinear form on the left hand side

Application of Grönwall's inequality leads to
As in the proof of Lemma 15, we choose ε such that −2α a < β(ε) < 0. Going back to (36a), we find We see that the right hand side consists only of exponentially decaying terms, which gives the decay rate for h ∆ in · H 2 (C) .
Theorem 3. Every stationary solution h, ρ a , ρ i of Problem 6 with h < h * a. e. is locally exponentially stable under disturbance that fulfills the mass equality condition, i. e., a solution (h, ρ a , ρ i ) of Problem 6 in a sufficiently small neighbourhood of h, ρ a , ρ i with´C ρ 0 a + ρ 0 i dx =´C ρ a + ρ i dx converges exponentially fast in time to h, ρ a , ρ i .
Proof. The linearised problem of Problem 6 in a sufficiently small neighbourhood around h, ρ a , ρ i is Problem 7. Due to Lemma 11, we only need to show exponential stability of Problem 7 in zero. But this is follows from the results in Lemma 12, Lemma 15, Lemma 16 and we are finished.

Singular Limits
We are going to have a closer look at stationary solutions in the limit ϑ 0. This way, we also rediscover the model of [27] as specialisation of our model. Henceforth, we restrict to the special disconnection rate r ϑ (h) = h−h * ϑ + . with lim n→∞ ϑ n = 0, there exists a subsequence (ϑ n k ) k∈N with h ϑn k H 2 (C) h 0 , ρ ϑn k a H 1 (C) ρ 0 a , and ρ

Singular limit of the stationary system
for all ϕ ∈ H 2 (C), σ a , σ i ∈ C c (C), where r 0 ∈ rca (C). (We identify a Radon measure and its density function w. r. t. the Lebesgue measure in the following.) Moreover, a. e.
Proof. Let (ϑ n ) n∈N be a zero sequence. independent of ϑ n . By testing (26b) with ρ ϑn a , we derive an H 1 (C) bound on ρ ϑn a uniformly w. r. t. ϑ n (the critical term including r ϑn can be dropped due to its non-positivity). We employ Lemma 8 again; this time to see that η a ∇ C ρ ϑn a = −η i ∇ C ρ ϑn i and conclude that the H 1 (C) norm of ρ ϑn i is as well bounded uniformly w. r. t. ϑ n . Lemma 7 assures the boundedness of h ϑn H 2 (C) independently of ϑ n . In H 1 (C) 2 , the theorem of Banach-Alaoglu gives us a weakly convergent subsequence (ρ αn a , ρ αn i ) n∈N to ρ 0 a , ρ 0 i ∈ H 1 (C) 2 as well as going to another subsequence h βn n∈N of (h αn ) n∈N does in H 2 (C). We use the Rellich-Kondrachov theorem to single out another subsequence (ρ γn a ) n∈N converging in L 2 (C) to ρ 0 a (since weak convergence in H 1 (C) implies weak convergence in L 2 (C) to the same limit). Analogously, we have convergence of a subsequence h δn n∈N to h 0 in L 2 (C), which includes another subsequence h ζn n∈N converging pointwise to h 0 . Due to weak convergence, a h ζn , ϕ n→∞ −→ a h 0 , ϕ for any ϕ ∈ H 2 (C). Moreover, To retrieve (39b) and (39c), test (26b) with a smoothed signum S ε ∈ H 1 (C) of ρ a with S ε H 1 (C) −→ sgn • ρ ζn a . We then have . Therefore, for a constant C 1 > 0. So in the limit ε → 0, we get the estimate r ζn h ζn ρ ζn a L 1 (C) ≤ C 2 for a constant C 2 > 0 due to the a priori bound on ρ a and ρ i . This implies weak--convergence in rca (C) of a subsequence r ζn k h ζn k ρ ζn k a to a Radon measure r 0 ∈ rca (C). As h ζn k and h * are in ∈ rca (C) and h * − h ζn k + r ζn k h ζn k ρ ζn k a weak--converges to h * − h 0 + r 0 in rca (C). We observe that and h * − h 0 + r 0 = 0.

Model without diffusion
We now turn to the system (26) with η a = η i = 0, which will lead us to a model of [27] in the Γ-limit ϑ → 0. Recall, that a sequence of functionals F n : X → R, n ∈ N, defined on a topological space Γ-converges to a functional F : Vanishing diffusivities have not been treated in our previous existence proofs for the time-dependent or stationary case. We attack this issue by reducing the system (26) to minimising an energy functional. In this order we choose a function ρ ∈ H 1 (C) , ρ ≥ 0 a. e., and formally substitute ρ i = ρ − ρ a into (26b). This way, we obtain: where . Inserting into (26a), we have only one equation left: whose solutions are the critical points of the following energy functional Lemma 17. There exists a minimiser of J ϑ .
Proof. This functional is coercive in the H 2 (C) norm: C 1 > 0, and we choose ε smaller enough for ε h 2 L 2 (C) to be absorbed by Note that A is weakly lower semicontinuous in H 2 (C). We show that B is weakly continuous in H 2 (C), so A + B is weakly lower semicontinuous in H 2 (C).
B is continuous in C b (C): Take a sequence (h n ) n∈N converging in C b (C) to h. Now observe that for all δ > 0, there exists an N ∈ N such that for all n ≥ N , it holds Since H 2 (C) comp → C b (C), for every sequence (h n ) n∈N converging weakly in H 2 (C) to h, we may take from every subsequence a subsubsequence (h n k ) k∈N that converges in C b (C) to h . It holds h a. e.
= h as weak convergence in H 2 (C) implies weak convergence in L 2 (C) and convergence in C b (C) implies (strong) convergence in L 2 (C). So with the previous result, With this lemma we have proven that for every ρ ∈ H 1 (C) , ρ ≥ 0 a. e., we find a minimiser h of J ϑ with which we may then define ρ a = g ϑ • h, and further ρ i = ρ − ρ a such that (h, ρ a , ρ i ) solves (26) with vanishing diffusivities. It is easily checked that the constructed linker densities are non-negative.
Lemma 18. For ϑ 0, the energy functional J ϑ Γ-converges to where g 0 (x) = ρ 2 min{h(x) 2 , (h * ) 2 }. Proof. As the other terms in J ϑ are independent of ϑ, we only need to consider the functional It is first shown that from any (ϑ n ) n∈N → 0 and any (h n ) n∈N that converges in H 2 (C) to h, a subsequence such that G n (h n ) → G 0 (h) can be singled out. (For the sake of notational simplicity, we abbreviate G ϑn = G n and g ϑn = g n .) In this order, rewritê We take a subsequence such that h n k converges in C b (C) to h. For any δ > 0 consider a sufficiently large k such that This term converges to zero for δ → 0 since g n k L ∞ (C) is uniformly (in k) bounded. Now consider Due to the monotonicity of r, g n (s) monotonically decreases to zero and g n ≤ g 1 ∈ L 1 (C), the monotone convergence theorem applies, so the last term converges to zero for n → ∞. Therefore, Take any clustering point of G n (h n ) and a subsequence G n k (h n k ) converging to it. Using the previous result, we single out another subsequence that converges to G 0 (h), which has to be the clustering point, so G 0 (h) ≤ lim inf n→∞ G ϑn (h n ).
Choosing an arbitrary h and taking the sequence h n = h converging to it, the previous considerations also lead to G 0 (h) ≥ lim sup n→∞ G n (h n ) and the claimed Γ-limit is shown.
Lemma 19. The Euler-Lagrange equation of J 0 is given by being the Heaviside function.
Proof. Let h be a stationary point of J 0 , i. e., d d ε (J 0 (h + εv)) 0 = 0 for all v ∈ H 2 (C). Take an arbitrary v ∈ H 2 (C) and consider ε > 0 small enough such that h < h * a. e. implies h + εv ≤ h * a. e. Then calculate The other derivatives are standard and the claim follows.
We may summarise the results of this section as follows: Theorem 5. i) For a total mass density ρ ≥ 0, and a ripping parameter ϑ > 0, any minimiser h of J ϑ constitutes a (variational) solution to Problem 3 for η a = η i = 0 together with some ρ a , ρ i ∈ H 1 (C). ii) Sending ϑ to zero, minimisers of J ϑ converge to variational solutions of the model of [27], cf. p. 2, Equation (2) therein.
Proof. i) Take any ρ and ϑ. Existence of minimisers of J ϑ is guaranteed by Lemma 17. As critical points, the minimisers solve (42). Defining ρ a according to (41) and ρ i = ρ − ρ a , we get solutions of Problem 3 for η a = η i = 0.
ii) follows directly from Lemma 18 and Lemma 19.

Numerical examples
In this section, we discuss results from numerical simulations of the parabolic PDE system we analysed in the previous sections. The predicted bleb size after a typical bleb formation time is compared against heights observed by biologists. Furthermore, we investigate the role of the critical pressure defined in [27] for static systems in the case of our time-dependent PDE system. As the tabular Table 1 suggests, quite a lot parameters of our model have already been assessed and discussed in the literature. The interested reader may consult the given sources and the references therein as a thorough treatment of these parameters is not in the scope of this work. Nevertheless, the choice of γ requires a comment: We decided to follow [27], [2] and set H 0 = 0. Despite their choice of κ = 2 · 10 −19 J and γ = −2 · 10 −6 , we stay consistent with our computations in Section 3.3, p. 9 and take γ ∼ − κ R 2 . In correspondence to [2], we choose the linker density ρ a (0, ·) to be 10 14 m −2 and start with ρ i (0, ·) = 0, accordingly. homogeneous Neumann boundary conditions. This domain approximation is a simple approach towards the singularities at the poles of ω 2 R and there are other more elaborate methods like surface finite elements (cf. [16]) which can handle such problems. However, we will show that we can in fact recover established results and conclude that this approach is sufficient for studying the bleb height near the equator. The fluid influence is neglected for simplicity setting L = cI.

Discretisation
In order to avoid a H 2 (D) trial space, but to stick with conformal Galerkin methods, we use a standard splitting of the ∆ 2 operator with boundary conditions for h and ∆h by introducing w = −∆h. This leads to the following system: In the Galerkin approximation, we employ Lagrangian P 1 finite elements conforming with the trial space H 1 (D) for h, w, ρ a , and ρ i . The mesh width is denoted by ∆x. In this configuration, we obtain a semidiscretisation of the form For implementation of this spatial discretisation, we used the FEM solver Netgen/NGSolve (https: //ngsolve.org, [33]). Discretisation in time is achieved by applying a semi-implicit Euler scheme with time step size τ > 0 and time points 0 = t 1 < t 2 < · · · < t k < t k+1 < · · · < T , n ∈ N: To cope with the implicit terms and the nonlinearities, we use Newton's fixed point iteration.

Scenarios
For all the simulations we will discuss in the following, the parameters of the PDE are as in Table 1. The typical expansion time for a bleb to nucleate is about 30 s [11], so we normalised the simulation time with respect to this reference quantity. The cortex is modelled (as before in the analytic part) as a sphere with radius R = 10 −5 m [11].
Nucleation of a bleb We are interested in the (maximal) height of the bleb that is nucleated after this time. It has been observed that the typical height of blebs is about 2 µ m, cf. [11].
We prescribe a pressure as the function x → 10 3 e − d(x,m) 2r 2 Pa, where d is the geodesic distance between the argument x ∈ S 2 and a midpoint m ∈ S 2 and r controls the width of the pressure pulse. The pressure is constant in time. The scaling ϑ of the disconnection rate r ϑ (h) is still free. A parameter study shows that with ϑ = 1.1 · 10 −12 , we may achieve a bleb height of 1.57 · 10 −7 , which is about one tenth of the experimentally observed bleb height. With a more elaborate view on the protein distribution at the membrane and especially their behaviour after the critical height h * is passed, one may achieve better results. For some time points, we plotted the membrane height as a heat map on the cortex, see Figure 2.
Remark 8. Considering the difference between the maximal height of every time step (see Figure 3), there seems to be numerical indication of a stability property like that rigorously shown in Section 5.3, which does not apply in this case as H 0 = 0 and a (·, ·) is not coercive.
Critical pressure [27] consider a bleb to form when the the membrane height reaches above the critical height h * (on a certain interval I ⊆ R) and linker bonds are broken in response. According to this notion of a bleb in their static model, they define the critical pressure to be the greatest pressure below which the membrane height is beneath the critical height everywhere. We adopt this notion to our (dynamic) model in the sense that the critical pressure p * 0 is the greatest value below which the maximal height of the membrane is beneath the critical height after the nucleation phase of 30 s which is driven by the pressure function Pa.
Passing the critical height triggers the linker disconnection process. By applying the pressure as    previously described and increasing the pressure in steps of 2 Pa, we plotted the maximal height after the nucleation phase against the pressure, see Figure 4. Considering a linear function through the first two data samples gives evidence that the data samples do not grow linearly at least not on the whole pressure range; up to a pressure of about 60 Pa linear growth seems to be a good description. We hypothesise that there is a change in the growing behaviour and that this change occurs at the critical pressure. To give some evidence for this, we successively fitted linear functions (with the gnuplot implementation of the Marquardt-Levenberg algorithm) to the datasets of the pressure intervals [0 : q], q ∈ {2, . . . , 120}, and assessed the root mean square residuum (RMS), see Table 2. We notice that until 16 Pa there is very little change in the RMS. From 16 Pa to 18 Pa there is a specifically large increase of the RMS (three orders of magnitude). Afterwards, the RMS increases readily. Taking a look at, Figure 4, we see that 16 Pa happens to be the critical pressure which substantiates our hypothesis that reaching the critical pressure triggers a major change in the bleb growth behaviour.

Conclusion
We derived a PDE model by balancing the bending, stretching and linker forces coming from the variational derivative of an extended Helfrich energy functional with the stress at the interface between the cytosol and the extracellular fluid. Based on the restriction of only small membrane displacement normal to the cortex, we could derive a gradient flow describing the membrane height normal to the cell cortex. Additionally, linker kinetics were incorporated with reaction-diffusion equations where we also introduced the concept of inactive linkers to include the phenomenon of cortex disruption. To our knowledge, this effect cannot be modelled with any other model. For the resulting system, we established global-in-time existence and uniqueness of weak solutions by applying the Banach fixed point theorem. The stationary case can also be treated with a fixed point argument, which uses the Schauder fixed point theorem, to establish existence. However, we do not have results about uniqueness or at least classification of stationary solutions.
The a priori estimates used in the stationary solution existence proof could be exploited for passing to the limit in the rescaling parameter of the disconnection rate, ϑ 0. We observed that the model of [27] could be rediscovered this way. The existence of a singular limit in the time-dependent case remains an interesting open question for future research.
Finally, let us mention that so far the model is purely mechanistic and ignores any interaction with external and cell-internal signalling. In particular, the interaction of bleb formation with polarization of protein distributions influencing the mechanical properties will become important in order to fully understand cell migration by blebbing.
Proof. We argue by a Petrov-Galerkin-type approximation: 1) Let (ϕ i ) i∈N be an orthonormal Schauder basis of L 2 (C) with eigenvalues (λ i ) i∈N (sorted ascendingly) consisting of eigenfunctions of the Laplace-Beltrami operator ∆ C (due to the divergence theorem and ∂C = ∅, the eigenfunctions ϕ i , i ≥ 2, are mean value free; for a spectral theorem on Riemannian manifolds cf. [26], Theorem 4.3.1). For m ∈ N, we set h m (t, . Formally inserting into (47) and testing with ϕ j , j ∈ {1, . . . , m}, we obtain the finite-dimensional system for the coefficient vectors h m = (h i ) i∈{2,...,m} and right hand side where A = diag (λ 2 , . . . , λ m ) and M = (L (ϕ i ν C ) , ϕ j ν C ) L 2 (C) j∈{1,...,m},i∈{2,...,m} . Due to the symmetry and positive definiteness of L((·)ν C ) · ν C in L 2 (C) (see Lemma 20), Leaving out a (h m , h m ) (non-negative term, cf. Remark 3), and applying the Cauchy-Schwartz inequality on the right, we arrive at with L = S 2 . We recall the continuity of L (see Lemma 20) and observe (L(h m ), h m ) L 2 (C) = With Young's inequality, integration in time, and the previous bounds, we obtain a bound of h m in L 2 [0, T ], H 1 (C) uniformly w. r. t. m. All together there is a subsequence with indices m k such that h m k weakly converges in H 1 (C) to h ∈ H 1 mvf (C), L(h m k ) converges weakly in L 2 (C) to L(h ), h m k weakly converges in H 2 (C) to h ∈ H 2 mvf (C) and f m k weakly in L 2 (C) to f for almost every t ∈ [0, T ). Uniqueness follows with the linearity of the equation.
(ii) The mean value ffl C h dx is constant in time.
Then consider the solutionh of (47) with initial datah 0 and right hand sidef . Set h =h + ffl C h 0 dx and observe (L (∂ t hν C ) , ϕν C ) L 2 (C) + a (h, ϕ) = L ∂ th ν C , ϕν C so h is a solution as claimed.
(ii) ffl C h 0 dx is constant in time and this is by construction the mean value of h.

D Differential geometry
This section contains basic differential geometric formulae which are eventually used for calculating the shape derivatives in Section 3.3. Some simple algebraic observations: Lemma 23. Let Γ be an orientable differentiable real submanifold with normal field ν and f : Γ → R m a function with existing tangential Jacobian. Further, P Γ = I − ν ⊗ ν. It holds, An analogue of the Schwarz theorem for tangential gradients: Lemma 24 ([6, Lemma 15]). Let Γ ⊆ R n be a differentiable, orientable real submanifold with normal field ν and f : Γ → R a function with existing tangential derivatives up to second order.

D.1 Commutator rules
Lemma 25 (Commutator rule for the tangential gradient, [6,Lemma 38]). Let X : I × R n → R n+k , for an interval I ⊆ R, be parametrisations of n-dimensional manifolds Γ s with an associated material derivative ∂ • s and velocity fields V r = ∂ s (X (s, ·))| r • (X (r, ·)) −1 . Let f : G (Γ s ; s ∈ I) → R be sufficiently regular. It holds A similar rule exists for the tangential divergence: Lemma 26 (Commutator rule for the tangential divergence, [6,Lemma 38]). Let X : I × R n → R n+k , for an interval I ⊆ R, be parametrisations of n-dimensional manifolds Γ s with an associated material derivative ∂ • s and velocity fields V r = ∂ s (X (s, ·))| r • (X (r, ·)) −1 . Let f : G (Γ s ; s ∈ I) → R n be sufficiently regular. It holds We can now derive a commutator rule for the Laplace-Beltrami operator: Lemma 27 (Commutator rule for the Laplace-Beltrami operator). Let G (Γ s ; s ∈ I), I ⊆ R, be a smooth evolving manifold and f : G (Γ s ; s ∈ I) → R a function with tangential derivatives up to order two whose material derivative exists and let f be material differentiable itself. It holds, T · ∇ Γr f + 2 ∇ Γr V r − ∇ Γr V r T − P Γr ∇ Γr V r : ∇ 2 Γr f + ∇ Γr V r : ν r ⊗ H r ∇ Γr f . ∇ Γr V r − 2J Γr (V r ) : ν r ⊗ H r ∇ Γr f = ∇ Γr V r − P Γr ∇ Γr V r − ∇ Γr V r T P Γr : ν r ⊗ H r ∇ Γr f = ∇ Γr V r : ν r ⊗ H r ∇ Γr f by the same arguments as in (5).
E Derivatives E.1 Material derivative of the mean curvature Lemma 28 ([6,Lemma 39]). Let G (Γ s ; s ∈ I) be an orientable, smooth evolving manifold. We then have E.2 Second order derivative of the integral of the mean curvature Lemma 29. Let G (Γ s ; s ∈ I), I ⊆ R, be a smooth evolving hypersurface in n dimensions. It holds, Proof.  (2): as V s • X s is independent of s and (note X = I Γ ) due to the symmetry of ∇ Γ ν.