{'fraction_non_alphanumeric': 0.0561070174370217, 'fraction_numerical': 0.02008606865192944, 'mean_word_length': 3.920936969333733, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Variational Approach to Particles in Lipid Membranes

A variety of models for the membrane-mediated interaction of particles in lipid membranes, mostly well-established in theoretical physics, is reviewed from a mathematical perspective. We provide mathematically consistent formulations in a variational framework, relate apparently different modelling approaches in terms of successive approximation, and investigate existence and uniqueness. Numerical computations illustrate that the new variational formulations are directly accessible to effective numerical methods.


Introduction
The interplay of proteins and curvature of lipid bilayers is well-known to regulate cell morphology and a variety of cellular functions, such as trafficking or signal detection. Here, interplay not only means that proteins can induce curvature by shaping and remodelling the membrane, but also that the membrane curvature plays an active role in creating functional membrane domains and organizing membrane proteins, including their conformation dynamics, [46,50,68]. Moreover, microscopic causes, such as hydrophobic mismatch of proteins and amphiphilic lipids, may have macroscopic effects, such as budding or fission. For example, the membrane remodelling during endocytosis is a consequence of the interplay between the elastic membrane and concerted actions of highly specialized membrane proteins that can both sense and create membrane curvature. In the case of clathrin-mediated endocytosis more than forty different proteins are involved, many of which are only transiently recruited to the plasma membrane [51].
A popular particle-based approach to the simulation of proteins in lipid membranes are coarse-grained molecular models in which the membrane constituents, i.e., lipids and proteins, are represented by short chains of beads, typically consisting of 3 -10 beads. Often the coarse-grained models are simulated using dissipative particle dynamics (DPD), with pairwise dissipative and random interactions between the beads that locally conserve momentum and yield the correct hydrodynamics. DPD simulations can reach time and length scales beyond those available by traditional molecular dynamics and thus allow for studying, e.g., cluster formation, budding or anomalous diffusion of transmembrane proteins [61,64,67]. Moreover, recent coarse-grained models incorporate structural and mechanical properties from experiments and all-atom bilayer simulations and claim to be 'semiquantitatively consistent' with experiments even on larger time and length scales [9,44,63,70].
On the macroscopic side of the model hierarchy, there are pure continuum models based on the fundamental Canham-Helfrich (CH) model of lipid membranes [12,34]. The associated Canham-Helfrich (CH) energy is obtained by expansion of the bending energy with respect to the principal curvatures up to second order. It is noteworthy that the CH energy is related to the much studied Willmore functional in differential geometry [43,73]. The CH model describes equilibrium and close-toequilibrium properties of biological membranes, such as the entropic repulsion of fluctuating membranes or the equilibrium shapes of lipid membranes. It has been modified so as to include the effect of different lipid components or large-scale protein assemblies. Lipids and proteins are then represented by areal concentrations. Their impact on the membrane morphology is modelled by concentration-dependent bending rigidities, phase-dependent spontaneous curvature, and a line energy associated with the phase boundaries [39,47]. The Euler-Lagrange equations associated with CH energy functionals usually are fourth-order, nonlinear (geometric) partial differential equations. These equations are typically coupled with nonlinear partial differential equations describing the evolution of lipids and proteins on the surface. The numerical solution of these coupled systems is an emerging field of current research (for an overview see, e.g., [15,21,22,28]). In numerical computations, binary and multicomponent models of membranes were able to reproduce the typical equilibrium shapes of vesicles, such as dumbbells, discocytes or starfishes (see, e.g., [6,20,24,25,36] and the references cited therein).
Hybrid models finally intend to bridge the gap between molecular dynamics based models that are expensive and have certain limitations in terms of the accessible time and length scales, and continuum models with continuous densities of proteins that do not incorporate the effect of small numbers of discrete particles. The basic idea is to treat proteins as discrete rigid particles coupled to the continuous membrane by suitable interface conditions. Aiming at moderate length scales, the continuous membrane is often described by a linearised CH model (Monge gauge). We emphasize that only mechanical curvature-induced interaction of particles and no chemical forces are considered in this context. However, additional chemical forces between particles could be incorporated by appropriate additional potentials and, under suitable conditions on these potentials, our mathematical analysis can be easily extended to this case. While the history of hybrid models can be traced back to the early nineties and meanwhile became a well-established field of research in theoretical physics (see, e.g. [17,18,30,33,37,40,48,55,58,60,71,76,77] only to mention a few of a multitude of references) the mathematical and numerical analysis of hybrid models is still in its infancy.
This paper is devoted to a mathematical consistent formulation of a variety of existing and some new models which reveals a kind of hierarchy of different modelling approaches in terms of successive approximation and is directly accessible to effective numerical methods. We start with, in a sense, most detailed hybrid models of proteins in lipid membranes that are based on the coupling of rigid particles of finite size to a linearised CH membrane by suitable conditions on the membrane displacement and its normal derivative (angle condition) along the boundaries of the particles [30,33,37,40,58,71]. The boundary values are determined by the specific interaction at the interface between the hydrophobic belt of transmembrane proteins and the surrounding lipid bilayer. We derive an equivalent reformulation of the corresponding fourth order boundary value problem by incorporating the boundary conditions in terms of constraints along the particle contour lines (curve constraints) in a similar way as in fictitious domain methods [29]. Variation of height and tilt is represented by additional degrees of freedom. This approach and most of our theoretical results also extend to rotation of particles which, however, will be considered elsewhere [75]. For fixed location of particles we prove existence and uniqueness of minimizers of the CH energy under curve constraints, together with convergence of a penalty formulation for vanishing penalty parameter. A gradient flow approach to varying locations of particles with soft-wall constraints is discussed only briefly as more detailed investigations are the subject of current research.
Approximating the conditions along the particle boundaries by their (constant) mean value, we introduce a novel class of finite-size particle models. The resulting averaged hybrid models no longer impose any conditions on the contour of the membrane along the particle boundaries and no longer provide a representation of tilt. By Green's formula, the angle conditions along the particle boundaries now take the form of averaged curvature conditions over the area of particles. In addition to existence and uniqueness for corresponding curvature-constrained minimization problems and its penalized counterparts, we also prove existence of global minimizers for varying locations of particles with hard-and soft-wall constraints.
It is sometimes convenient, both for analytical and numerical purposes, e.g., when proceeding to larger copy numbers or to larger spatial scales, to approximate averages of curvature over the whole area of finite-size particles by weighted point values of curvature, e.g., in the barycenters. The resulting point-like hybrid models have quite a history in theoretical physics, see, e.g. [7,17,18,40,48,53,55,72] and also [30,58]. However, as point curvature constraints are not well-defined for functions with only second order weak derivatives, the resulting problem is not well-posed from a mathematical point of view. While this issue is often addressed by suitable truncation of Fourier expansions in physics literature (see, e.g., [18]), we consider an extension of the linearised CH energy by regularizing third and fourth order terms (cf. e.g., [7]). From a physical point of view these additional terms could be justified in terms of higher moments in the expansion of the membrane bending energy [52]. Such kind of extended CH energy then allows for existence and uniqueness results for point-like particles that are in complete analogy to averaged hybrid models. We also provide a representation of global minimizers in terms of suitable Green's functions. For unbounded domains, we thus recover a representation of solutions that was first suggested by Bartolo and Fournier [7].
For particles interacting with the cytoskeleton, we consider two different classes of models, describing the interaction with the membrane by point values or by point forces [26,31,66]. As point values of functions with second-order weak derivatives are well-defined, the classical linearised CH membrane energy is used in both modelling approaches. We prove existence and uniqueness for prescribed point values at fixed locations and derive representations of the solutions in terms of suitable Green's functions. In addition to existence of global minimizers for varying locations, we prove that solutions for N particles are entirely determined by the two-particle problem considering only the particles with the largest and smallest prescribed point values. Again, we derive representations of the solutions of these problems in terms of Green's functions. Similar results are presented for prescribed point forces. Apart from existence and uniqueness results for fixed location of particles, we derive a representation of global minimizers in terms of Green's functions and prove clustering of point forces to one or two clustering points depending on the sign of point forces. In this way, computation of global minimizers with N point forces can be reduced to a problem with at most N = 2 point forces.
The paper concludes with some numerical experiments. First, we compute the interaction potential for two circular or elliptical particles as described by the finitesize hybrid model with coupling conditions on the membrane displacement and its normal derivative at the particle contour lines and compare the results with those obtained from a related point-like model. As a second scenario, we numerically investigate the interaction potential of two point forces and, in particular, illustrate how clustering depends on the ratio of bending rigidity and membrane tension.
The paper is organised as follows. After a short introduction to the biophysical background in Section 2, Section 3 summarizes the mathematical description of lipid membranes by the CH energy. Section 4 is devoted to finite-size hybrid models performing the coupling by conditions on the membrane displacement and its normal derivative (angle condition) at the particle boundaries while Section 5 concentrates on averaged hybrid models arising from approximations of these boundary conditions by mean values of curvature. Hybrid models with point-like particles are introduced by approximating those mean values by point curvature constraints and are analysed in Section 6. Models for particles interacting with the cytoskeleton by point values and point forces are considered in Section 7. The paper concludes with some numerical experiments collected in Section 8. The theoretical findings in Section 4 -7 are based on an abstract variational framework that is presented in the Appendix together with some required regularity results on Green's functions.

Particles in lipid membranes
A biomembrane is a thin layer that surrounds biological cells and cellular organelles acting as a barrier between the cell and its surroundings. It consists of a lipid bilayer with embedded and attached proteins. The lipid bilayer is composed of phospholipids, made up of a phosphate group and a diglyceride. The resulting form is that the molecules have a head, the phosphate group, and a tail comprised of the two fatty acid hydrocarbon chains that make up the diglyceride. The heads are hydrophilic and the tails hydrophobic. As a consequence, when placed in water these molecules form structures in which the heads point outwards and the tails inwards. There are a number of such formations, but the one of interest here is the bilayer sheet. This is when the phospholipids line up so that the heads form two distinct layers, sandwiching the tails between them as depicted in Figure 2.1 (taken from [50]).
The biological membrane is composed of different types of lipid molecules, which may differ in their head groups, in the length of their hydrocarbon chains, or in the number of unsaturated bonds within these chains, which is associated with different conformations (phases) leading to order (dense packing) or disorder (light packing). Since the different lipid types and phases tend to separate, the composition of such a multicomponent membrane can become laterally inhomogeneous.  Another crucial source of inhomogeneity are proteins that are responsible for the diverse functions performed by biological cells, which range from the transport of specific molecules across the membrane to the reception of chemical signals from the extracellular environment. Transmembrane or embedded proteins, at least partially, penetrate the membrane while peripheral proteins, (for example BAR domains [50]), are attached to it by electrostatic interaction or other weak binding forces. A pictorial description of a transmembrane protein is provided in the left picture of Figure 2.2 (taken from [50]). Interaction of embedded proteins with the membrane is characterised by the shape of their hydrophobic belt whose width might be different from the membrane thickness and might vary along the particle contour line. In contrast, peripheral proteins usually impose their shape and thus curvature to the lipid membrane like a scaffold. Both embedded and peripheral proteins may not only tilt and move up and down with the membrane, but also move laterally and eventually cluster according to mechanical forces induced by membrane curvature.
Cytoskeletal protein assemblies such as actin filaments or bundles of actin filaments locally impose membrane curvature by applying a point force or constraining the membrane to take a fixed height. A cartoon picture is given in Figure 2.2b (taken from [50]). Branching, bundling, and treadmilling of actin filaments is, e.g., responsible for the formation of filopodia which are slim, finger-like projections of a cell membrane that are formed by the bundling of actin filaments attached to the cell cytoskeleton at one end and pushing against the cell membrane at the other [49]. For further information on the mechanical interaction of particles in lipid bilayers we refer, e.g., to the overviews [50,62] and the literature cited therein.

Mathematical description of lipid membranes
3.1. Canham-Helfrich free energy. As the width of a bilayer (10 −9 m) is much smaller than its lateral extension, it is natural to model the membrane as a smooth, two-dimensional hypersurface M embedded in the three-dimensional Euclidean space R 3 . Note that this simplification neglects transversal stretching and transversal shearing as possible elastic deformations. Since the fluidity of the membrane excludes lateral shearing, the deformations of the membrane are caused by lateral stretching and bending.
The mathematical study of biomembranes principally concerns the minimization of the energy functional describing the energy associated to displacements of M. Fundamental to the macroscopic approach to modelling biomembranes is the Canham-Helfrich (CH) model [12,27,34] which is based on the expansion of the bending energy with respect to the principal curvatures up to second order. It describes equilibrium and close-to-equilibrium properties of biological membranes. The fundamental object of the CH model is the elastic bending energy J CH defined by Here H and K stand for mean and Gaussian curvature of the membrane M ⊂ R 3 while κ and κ G are the corresponding bending rigidities and dM is the surface element of M. The additional parameter c 0 is called spontaneous curvature and accounts for a possible asymmetry between the outer and inner layers in the otherwise flat reference configuration, e.g., due to different lipid compositions in the layers. The related energy supplements the bending energy with a surface energy´M σ dM that is associated with membrane tension σ ≥ 0. Here, the surface energy penalises area change and thus accounts for the incompressibility constraint of the fluid membrane in the lateral direction. These energies, depending on the type of problem, may be augmented with reduced volume or bilayer area difference constraints [65]. The mathematical derivation of Canham-Helfrich-type models from molecular descriptions of lipid bilayers by Γ-convergence techniques has been started only recently [8,59]. Note that for a general membrane M the Gaussian curvature term´M κ G KdM gives a non-constant contribution to J CHS (M). However, assuming that M is closed or that the geodesic curvature along ∂M-composed by the membranes outer boundary and possible particle contour lines-is fixed, this term becomes a topological invariant by the Gauss-Bonnet theorem. Hence, it can be ignored when computing equilibrium membrane shapes minimizing J CHS (M). For a detailed discussion of geodesic curvature along particle contour lines we refer to [75].
3.2. Monge gauge. We will now outline a geometrically linearised approximation of the Canham-Helfrich energy J CHS defined in (3.2). For simplicity, let us assume that spontaneous curvature c 0 is zero, so that we have In the Monge gauge, one assumes that the surface is nearly flat, so that the membrane surface can be parametrized as a graph over a two-dimensional reference domain Ω ⊂ R 2 . Then, the mean curvature H and the Gauß curvature K of the membrane M are given by A common approach to derive an approximate model is to assume that the displacement of the membrane from the x − y plane produced by the particles is small, i.e. |∇u| 1. In this case, it is sufficient to consider the geometric linearisation (3.5) which models perturbations from a flat surface. Inserting the geometric linearisation (3.5) of mean curvature, Gauss curvature, and the surface element dM = (1 + |∇u| 2 ) 1/2 dx into (3.3) yields, up to a constant term, the quadratic energy Ignoring Gaussian curvature in light of the Gauss-Bonnet theorem discussed above, a quadratic approximation of the energy J CHS from (3.3) finally takes the form We emphasize that Gaussian curvature may still vary across the membrane allowing for saddle-like deformations (cf., e.g., the numerical results in Subsection 8.1).
Observe that minimization of J is leading to fourth order plate-like equations.
3.3. Boundary conditions and coercivity. We consider the Canham-Helfrich free energy J (u) in the Monge gauge (3.7) with membrane displacement u defined on a reference domain Ω ⊂ R 2 . From now on, we assume that Ω ⊂ R 2 is a bounded, convex domain with a piecewise smooth Lipschitz boundary ∂Ω, e.g., a square. Only in Section 6.5, we shall consider Ω = R 2 . Deformations u of the membrane are taken from a closed subspace V ⊂ H 2 (Ω) satisfying suitable boundary conditions. We consider the following three cases , often referred to as Dirichlet, Navier, and periodic boundary conditions with zero mean, respectively. Note that is Ω-periodic and´∂ Ω v ds = 0}. is only defined for a rectangular domain Ω. The space V is equipped with the canonical norm · 2 = · H 2 (Ω) in H 2 (Ω). Throughout the following, we assume that κ > 0 as well σ ≥ 0 for all three choices of V .
associated with the energy functional J is continuous and coercive on V .
Proof. While continuity of a( · , · ) is obvious, we refer to [32] for a proof of coercivity.

Curve constraints
4.1. Particles interacting with lipid membranes. Transmembrane proteins are interacting with the membrane curvature by the shape of the hydrophobic belt of single molecules or by oligomerisation, i.e., macromolecular clustering [50]. Other particles, such as peripheral proteins, may be attached to the membrane surface, acting as active or passive scaffolds [50], others may be partially wrapped due to adhesion energy [4,42]. All these phenomena can be captured by the same type of mathematical model that treats each particle as a finite-sized, rigid body B that interacts with the membrane by suitable conditions on the common interface Γ.

4.2.
Boundary values and curve constraints. We assume that the membrane occupies a graph over a subset of the given reference domain Ω ⊂ R 2 with non-empty, open subsets B i ⊂ Ω representing i = 1, . . . , N particles with diameter 2r i and centres of mass X i included in the membrane. For example, the particles B i could be (but do not have to be) circles with radii r i and midpoints X i . We assume that the particles have C 3boundaries Γ i = ∂B i . We also assume that the particles do not overlap and do not touch the boundary ∂Ω of Ω in the sense that B i ∩ B j = B i ∩ ∂Ω = ∅ for all i = j = 1, . . . N . In the sequel, we will consider varying the height and location of particles. We will also consider linearized variation of tilt (by neglecting the corresponding variation of B i ), but exclude rotation for simplicity. We assume that the deformation u defined on Ω B satisfies (artificial) boundary conditions on ∂Ω as specified implicitly in the preceding section (via the function spaces in which u lies) and additional boundary conditions on the particle boundaries Γ i that determine the membrane-particle coupling. For a (a) transmembrane protein, (b) scaffold, and (c) partly wrapped particle this is illustrated in Figure 4.1. For a transmembrane protein, the particle boundary Γ of B (dashed line) is the projection of the protein-membrane contact line γ ⊂ R 3 (only shown as cut with the mid-surface of the protein) to R 2 , the height h is the (spatially varying) distance from γ to R 2 , and −s describes the (spatially varying) angle between the bilayer mid-plane and the horizontal outward normal of B. Similarly, for a peripheral protein, B stands for the projection of the particle-membrane contact area, while this set is augmented by part of the membrane thickness for partly wrapped particles.  Assuming that height, tilt, and location of particles are fixed, these additional essential (Dirichlet) boundary conditions are given by Here, and throughout the paper, n denotes the outward normal to Ω B (while the inward normal was more convenient in Figure 4.1). Furthermore, We assume that the reference particles B 0 i ⊂ R 2 are open sets with C 3 -boundaries, X 0 i = 0 ∈ B 0 i , and thus X i ∈ B i . The essential boundary conditions (4.1) are incorporated in the affine subspace It is convenient to introduce the vector-valued functions h = (h i ) and s = (s i ).
The deformation u of the membrane with fixed height, tilt, and location of particles is the solution of the following minimization problem.
The minimization problem 4.1 is well-known to be equivalent to find u ∈ V h,s In order to show existence and uniqueness of a solution to Problem 4.1, we now derive its reformulation in terms of functions defined on the whole domain Ω.
To this end, it is convenient to introduce the particle trace operators T i , Proof. Let i = 1, . . . , N be fixed and recall that the boundary Γ i of B i is C 3 . As B i ∩ ∂Ω = ∅, a compactness argument shows that dist(B i , ∂Ω) > 0. Hence, there is a subsetΩ ⊂ Ω with a C 3 -boundary ∂Ω and B i ∩∂Ω = ∅. In particular,Ω i =Ω\B i then has a C 3 -boundary and classical trace theorems [74,Theorems 8.7 and 8.8] imply that the trace operator is linear, continuous, and surjective. As ∂Ω i = Γ i ∪ ∂Ω, this provides the assertion.
Note that denseness of C ∞ (Ω) in H 2 (Ω) implies that T i remains the same if constructed from the inside of B i instead of the outside, i.e., as trace operator for Proof. Lemma 4.1 provides linearity and continuity of T .
As above, classical trace theorems [74,Theorems 8.7 and 8.8] imply that there is a w ∈ V such that (w| Γi , ∂ ∂n w| Γi ) = v and thus surjectivity of T .
Note that the condition B i ∩ B j = ∅ for i = j is necessary for surjectivity of T . Now, we introduce the affine subspace Homogeneous data h i = s i = 0 provide the linear space V 0,0 ⊂ V . In this way, the boundary conditions (4.1) now become constraints on the curves Γ i ⊂ Ω. The equivalent variational formulation of Problem 4.2 is to find u ∈ V h,s such that We now show existence and uniqueness of solutions to Problem 4.2. Proof. By surjectivity, linearity, and continuity of T we know that V h,s is a nonempty, affine, closed, subspace of V . In addition, a( · , · ) is bounded and coercive on V by Lemma 3.1. Hence, the Lax-Milgram lemma provides existence and uniqueness.
We now clarify the relation of Problems 4.2 and 4.1.
Proof. By definition of V Ω B there must be an extensionṽ  to Ω by u i ∈ H 2 (B i ), i = 1, . . . , N , then satisfies u ∈ V h,s . Furthermore, for each v ∈ V 0,0 we have v B = v| Ω B ∈ V 0,0 Ω B and v i = v| Bi ∈ H 2 0 (B i ). As u B satisfies (4.4) and u i satisfies (4.8), this leads to Adding these equations, we obtain a(u, v) = 0. Hence, u is the solution of Problem 4.2.
Existence and uniqueness of a solution to Problem 4.1 is now an immediate consequence of Proposition 4.1.

4.3.
Parametrized curve constraints. We consider curve constraints of the form (4.1) which now allow for varying heights γ = (γ i ) ∈ R N of the particles B i , i = 1, . . . , N . More precisely, we impose the conditions Find (u, γ) ∈ V × R N minimizing the energy J subject to curve constraints (4.9).
We consider existence and uniqueness of a solution together with an equivalent variational reformulation of Problem 4.3.
and for each fixed i = 1, . . . , N there is an η i ∈ V , such that Moreover, v ∈ V satisfies curve constraints (4.9) with some γ i ∈ R if and only if Proof. The first two assertions follow from Lemma 4.2. Assume that v ∈ V satisfies the curve constraints (4.9) with some γ i ∈ R. Then, by definition of η i , we have and thus v ∈ η 0 + span{η i | i = 1, . . . , N } + V 0,0 . The converse is obvious. Proof. By Lemma 4.4, the Problem 4.3 is equivalent to minimizing J on the nonempty, affine, and closed subspace Hence, the Lax-Milgram lemma provides existence, uniqueness, and the variational formulation (4.13).
We now explicitly decouple the computation of u and γ using the orthogonal projections (4.14) with respect to the L 2 -scalar product ( · , · ) Γi = ( · , · ) L 2 (Γi)×L 2 (Γi) . The kernel of P Γi is given by ker P Γi = span{(1, 0)} ⊂ X i . The induced projection P = (P Γi ) : X → X on the product space X = N i=1 X i now provides a decoupling of u and γ.
Here, the dot stands for the Euclidean scalar product in R 2 and the two components of the varying angles can be equivalently written as Utilizing suitable orthogonal projections, the computation of u ∈ V and of the coefficients γ i , α i,j ∈ R can be decoupled in analogy to Problem 4.4 and Proposition 4.3. We refer to [75] for details.

4.4.
Varying the location of particles.

Hard and soft wall constraints. For varying locations
i we want to enforce that particles do not overlap, touch the boundary ∂Ω or even escape from the considered domain Ω. To this end, we constrain the locations X = (X i ) to the subset of all positions with non-overlapping particles contained in the domain. Recall that the particles B i (X) ⊂ R 2 are open sets with C 3 -boundaries and X i ∈ B i (X).
is open. Hence, ω can be written as As Ω is bounded, ω ⊂ Ω N is also bounded and therefore compact.
Later on, we will also require that particles do not touch neither each other nor the boundary. It can be shown by elementary arguments that the set of such particle positions is given by the interior of ω which can be characterized by The constraints X = (X i ) ∈ ω are enforced by additional terms in the energy. We first consider so-called soft-wall constraints V soft = V 1 + V 2 consisting of a Lennard-Jones potential This term accounts for the repulsion and attraction of particles. We also define This term is accounting for escaping particles. For circular particles we have The energy V soft associated with soft-wall constraints can be regarded as a penalization of hard-wall constraints V hard =V 1 +V 2 with Note that ω is the domain of V hard .

4.4.2.
A Gradient flow approach. With given reference particles B 0 i located at 0 and data (h i , s i ) prescribed on Γ i (X) = X i + Γ 0 i according to (4.2) (see Section 4.2), we now allow for varying locations X = (X i ) of particles B i (X) = X i + B 0 i , and fix height and tilt, for simplicity. We impose soft wall constraints V soft according to (4.18), (4.19) .

Problem 4.5 (Varying locations of particles with soft wall constraints).
Find (u, X) ∈ V ×int ω minimizing the energy J +V soft subject to curve constraints We describe a gradient flow approach to the iterative solution of Problem 4.5. For t ≥ 0 let X(t) be a trajectory of locations of particles B i (X(t)) with corresponding boundaries Γ i (X(t)), i = 1, . . . , N .
∂Xi,2 , we consider the gradient flow with suitable initial locations X 0 ∈ ω. Then, by construction, we have monotonically decreasing energy J (u(X(t)))+V soft (X(t)) along a solution X(t) of (4.21) and we may expect that (u(X(t)), X(t)) tends to a solution of Problem 4.5 for t → ∞. While the gradient ∇ X V soft (X) can be obtained by straightforward calculus, the derivation of ∇ X J (u(X)) requires some care. In particular, in light of a generic lack of smoothness of the minimizers u(X(t)) ∈ V h,s across Γ i , it is convenient to consider their restriction to Ω B . Utilizing basic techniques from shape calculus, we then obtain the following representation of ∇ Xi J (u(X)) (see [75] for details).

4.5.1.
Penalization of curve constraints. The numerical approximation of the curve constrained Problem 4.2 leads to saddle point problems involving suitable approximating spaces for primal and dual variables. As an alternative, we present an adaptation of the boundary penalty approach, [3,5] to essential boundary conditions, see also Nitsches approach [57]. More precisely, we penalize deviation from curve constraints (4.1) by the additional energy term with given penalty parameter ε > 0, the trace operators T i defined in (4.5), and Γi is the L 2 -norm on the trace space X i . Note that the resulting mathematical problem may be considered as a model in its own right by interpreting the resulting boundary conditions as physical contact conditions and as a modelling parameter. The counterpart of Problem 4.2 with soft constraints then reads as follows.
An equivalent variational formulation of Problem 4.6 is to find u ε ∈ V such that Proof. By Theorem 4.1 Problem 4.2 admits a unique solution u ∈ V h,s = {v ∈ V | T v = (h, s)}, and the trace operator T introduced in Lemma 4.2 is a bounded, linear map into X . Hence, the assertion follows from Proposition 9.3 in the appendix.

4.5.2.
Parametrized soft curve constraints. We consider a soft version of the parametrized curve constraints (4.9) involving variable heights γ i ∈ R of particles B i . The corresponding additional energy term reads with given penalty parameter ε > 0, given functions h i , s i as in (4.9), and unknown heights γ i ∈ R, i = 1, . . . , N . The corresponding minimization problem reads as follows.
Again, we use the projections P Γi defined in (4.14) to decouple the computation of u ε and γ ε . In particular, orthogonality of P Γi and (γ ε,i , 0) ∈ ker P Γi provides the reformulation Γi of the penalty term (4.24). Here, the second summand on the right hand side is vanishing for Proposition 4.7. The pair (u ε , γ ε ) ∈ V × R N solves Problem 4.7, if and only if u ε ∈ V is minimizing the energy and The assertion follows from a general result stated in Proposition 9.2. Proposition 4.7 provides an explicit representation of the heights γ ε = (γ ε,i ) in terms of the deformation u ε that in turn can be computed from a fully decoupled minimization problem. The variational form of the minimization problem is to find are continuous, the Lax-Milgram Lemma provides the assertion.
Proof. The same arguments as in the proof of Proposition 4.7 show that the assumptions of Proposition 9.3 in the appendix are satisfied and the latter provides the assertion.
} where the additional two functions describe tilt as in Remark 4.2. (4.2). In light of Proposition 4.7 and Remark 4.3, we assume that height and tilt of particles are fixed, for simplicity.

Varying the location of particles. We consider varying locations
We first consider hard-wall constraints V hard as defined in (4.20). This means that we constrain the locations X = (X i ) to the domain ω of V hard as given in (4.16).
Problem 4.8 (Soft curve constraints with varying locations and hard-wall constraints). Find (u ε , X ε ) ∈ V × ω minimizing the energy Here, the trace operators T (X) i are defined according to (4.5) but with Γ i replaced by Γ i (X). We set T (X) = (T (X) i ) in accordance with Lemma 4.2.
Now we augment the energy by the additional term V soft (X) associated with soft wall constraints that penalize overlapping and escaping particles (cf. Subsection 4.4.1). Problem 4.9 (Soft curve constraints with varying locations and soft wall constraints). Find (u ε , X ε ) ∈ V × int ω minimizing the energy We briefly describe a straightforward gradient flow approach to the iterative solution of Problem 4.9. As in Subsection 4.4.2, we consider a trajectory X(t) ∈ int ω of locations of particles with corresponding boundaries Γ i (X(t)), i = 1, . . . , N . Let u(X) ∈ V denote the unique minimizer of E(·, X) on V for given X = X(t) (cf. Proposition 4.5). Again denoting ∇ X = (∇ X1 , . . . , ∇ X N ) and ∇ Xi = ∂ ∂Xi,1 , ∂ ∂Xi,2 , we might consider the gradient flow with u = u(X) which, however, might not be applicable for lack of smoothness. Existence and approximation of global minimizers for varying locations will be considered in a separate paper.
4.6. Discussion. We have analysed various descriptions of proteins in lipid membranes in terms of rigid particles of finite size, interacting with a linearised Canham-Helfrich membrane model by suitable conditions on the membrane displacement and its normal derivative (angle condition) at the particle boundaries. The history of such kind of hybrid models can be traced back to the early nineties, see, e.g., [30,33,37,40,58,71]. For example, Problems 4.1 and 4.3 set out in Section 4.2 and 4.3 are equivalent to what is referred to as the 'strong coupling regime' in [30] or the 'microscopic model' in [58]. See, e.g., [56] for a careful discussion of the boundary conditions (4.1). Investigations particularly focus on the effect of the mechanical properties of the membrane on the interaction between different particles in interdependency with corresponding membrane deformations. Such kind of membrane-mediated interaction was studied for particles with parametrized height and tilt angle, circular or non-circular cross-section, respecting or breaking reflection symmetry [30,41,58]. Varying particle positions as treated in Section 4.4 have also been considered in this context, see, e.g., [40]. The notion of 'soft inclusions' as introduced from a physical point of view, also referred to as 'perturbative regime' [30] or 'phenomenological model' [58], does not completely agree with the notion of 'soft constraints' as considered in Section 4.5. Soft inclusions were originally introduced to study the influence of regions with excess concentration of lipids on fluctuating membranes, where it is appropriate to assume the bending rigidity of these regions to be close to the value of the surrounding membrane. Note, that regions of excess lipid concentration usually respect reflection symmetry (across the membrane) and thus, do not impose any constraints on the curvature. However, applying this approach to model 'soft', i.e., non-rigid proteins that break reflection symmetry, due to spontaneous curvature terms c i , leads to an energy which takes the form of soft point curvature constraints as considered in Section 6, when formally passing to the limit of point-like particles.

Averaged curve constraints
5.1. Fixed heights and locations. In order to derive approximations of boundary conditions or equivalent curve constraints (4.1) on Γ i , we introduce the averaging functionals f i , g i : V → R, The functionals f i , g i , are linear and bounded on V , because each T i is a linear and bounded map into Remark 5.1. Utilizing Green's formula, the functionals g i , i = 1, . . . , N , can be expressed as averaged mean curvature according to where the sign results from the fact that n is an inward normal to B i . Conditions on g i are therefore called mean curvature constraints in the sequel.
We introduce the averaged data with h i , s i given according to (4.2) and the notation As an approximation of the boundary value Problem 4.1 or its equivalent fixeddomain formulation Problem 4.2, we now consider the following minimization problem with average constraints.

Problem 5.1 (Averaged curve constraints).
Find u ∈ V minimizing the energy J subject to the average constraints Note that, using the componentwise integral operator´: X → R N ×2 defined by the constraints (5.5) can be rewritten as ffl i.e., as an average of the constraints in Problem 4.2 (cf. Theorem 4.1). Obviously, averaged constraints (5.5) or (5.6) are weaker than curve constraints (4.1).
Proof. The trace operator T = (T i ) is linear, continuous, and surjective by Lemma 4.2. Continuity of averaging´implies continuity of´T : V → R N ×2 . Therefore, the Lax-Milgram lemma provides existence and uniqueness in the non-empty, affine, closed, subspace See also Proposition 9.1 in the appendix.

5.2.
Averaged mean curvature constraints. We now investigate averages of parametrized curve constraints (4.9) of the form with h i , s i given according to (4.2) and freely varying heights γ = (γ i ) ∈ R N . In analogy to (5.6) the parametrized constraints (5.7) can be rewritten as i.e. as averages of parametrized curve constraints (4.9) occurring in Problem 4.3.
Utilizing the orthogonal projection P = (P Γi ) : X → X with P Γi introduced in (4.14), the constraints (5.8) can be decoupled according to Now observe that by definition of P = (P Γi ) the constraints hold for all v ∈ V . The remaining constraints in (5.9) and thus the parametrized constraints (5.7) take the form In light of g X (u) = ( ffl T u) i,2 , Remark 5.1, ands = ffl s, minimization of J subject to averaged parametrized constraints (5.7) finally amounts to the following averaged version of Problem 4.4.
Problem 5.2 (Averaged mean curvature constraints). Find u ∈ V minimizing the energy J subject to the constraints (5.10) g X (u) =s.
According to (5.9), the optimal heights γ = (γ i ) ∈ R N are obtained from once the solution u ∈ V of Problem 5.2 is available. Note that, by the representation (5.2), one can construct a function w ∈ V with g X (w) =s, even without the assumption B i ∩ B j = ∅. Hence, the following existence result is a consequence of the Lax-Milgram lemma.
Proposition 5.2. There exists a unique solution u ∈ V to Problem 5.2.

Remark 5.2. As a consequence of
the tilt of particles is no longer represented by parametrized constraints of the form (4.15) after averaging.

5.3.
Varying the location of particles. We now consider particles with varying location and hard-and soft-wall constraints as defined in Section 4.4.1. In case of hard-wall constraints, we allow the locations X = (X i ) to vary only in the domain ω of V hard as given in (4.16).
Problem 5.3 (Averaged mean curvature constraints with varying locations and hard-wall constraints). Find (u, X) ∈ V × ω minimizing the energy J (u) subject to the constraint Notice that Problem 5.3 is equivalent to minimizing J (u) + V hard (X) over V × ω under the constraint (5.12). For soft-wall constraints, we obtain the following related problem.
Problem 5.4 (Averaged mean curvature constraints with varying locations and soft-wall constraints). Find (u, X) ∈ V × ω minimizing the energy J (u) + V soft (X) subject to the constraint g X (u) =s.
The key ingredient to show existence of solutions for Problems 5.3 and 5.4 is the following lemma.
Lemma 5.1. Let M ⊂ R n be bounded and measurable and p > 1. Then the operator Proof. Assume that v ∈ L p (Ω) is extended by zero to R n and q = p p−1 is the dual exponent. Now let x, y ∈ Ω. Then Hölder's inequality gives Proof. We want to apply the general result of Proposition 9.4 in the appendix. To this end we first select an By Lemma 4.5, the set ω is compact. Furthermore it is easily checked that the functional V hard : R N ×2 → R ∪ {∞} is lower semi-continuous. In order to apply Proposition 9.4, it remains to show that X → g X ( · ) ∈ V is continuous. This follows from the representation of g X given in Remark 5.1, Lemma 5.1 with p = 2, and the continuity of ∆ : V → L 2 (Ω).
Under a slightly stronger assumption we also get existence for soft-wall constraints.
Proof. We want to apply Proposition 9.4 in the appendix. To this end, we now fix X 0 ∈ int ω = ∅. Then, equation (4.17) implies that we have dist(B i (X 0 ), B j (X 0 )) > 0 and dist(B i (X 0 ), ∂Ω) > 0 for i = j and thus V soft (X 0 ) < ∞. Now we can literally proceed as in the proof of Proposition 5.3, but with V hard replaced by V soft . This provides a minimizer (u, X) ∈ V × ω. Since V soft is infinite on ω \ int ω, we even have X ∈ int ω.
Note that we cannot expect uniqueness, e.g., for reasons of symmetry. 5.4. Soft averaged mean curvature constraints.

Fixed locations of particles.
We skip penalty approximations of the fixedheight Problem 5.1 and directly concentrate on averaged mean curvature constraints (5.10) and varying heights at fixed location of particles. We introduce the penalty term 2) and penalty parameter ε > 0. The corresponding minimization problem reads: Problem 5.5 (Soft averaged mean curvature constraints).
Find u ε ∈ V minimizing J (u) + 1 2ε g X (u) −s 2 R N . An equivalent variational formulation amounts to find u ε ∈ V such that (5.14) a(u ε , v) The Lax-Milgram lemma yields existence and uniqueness.
Proposition 5.5. There exists a unique solution u ε ∈ V to Problem 5.5.
Problem 5.5 can be regarded as an approximation of Problem 5.2.
Proposition 5.6. Let u denote the solution of Problem 5.2 and u ε denote the solution of Problem 5.5 for fixed ε > 0. Then we have Proof. The assertion is a direct consequence of Proposition 9.3 in the appendix.

5.4.2.
Varying the locations of particles. We first consider hard wall constraints V hard as defined in (4.20). This means that we allow the locations X = (X i ) to vary only in the domain ω of V hard as defined in (4.16). Recall that ω is compact in R N ×2 .
Proof. In the proof of Proposition 5.3 we have shown that Problem 5.3 satisfies the assumptions of Proposition 9.4 in the appendix. Under the same assumptions Proposition 9.5 in the appendix provides existence of a solution of its penalized analogue, i.e., Problem 5.6. Now we consider soft wall constraints associated with the additional energy term V soft (X) that penalize overlapping and escaping particles (cf. Subsection 4.4.1).
Problem 5.7 (Soft averaged mean curvature constraints with varying locations and soft wall constraints). Find (u ε , X ε ) ∈ V × int ω minimizing the energy Proposition 5.8. Assume that int ω = ∅. Then there exists a solution (u ε , X ε ) ∈ V × int ω to Problem 5.7.
Proof. Again we only need to note that the existence result of Proposition 9.5 for penalized problems is valid under the same assumptions as the one for the nonpenalized analogue in Proposition 9.4, and that we already verified these assumptions in the proof of Proposition 5.4. Hence we have a solution (u ε , X ε ) ∈ V ×ω that must also satisfy X ε ∈ int ω because we would otherwise have V soft (X ε ) = ∞.

5.5.
Discussion. To our knowledge, averaged hybrid models were not yet considered in the existing literature. They could be regarded as an intermediate approximation step from finite-size hybrid models, performing the coupling by boundary conditions at the particle boundaries (cf. Section 4), to well-established pointparticle models to be considered below (cf. Section 6). In particular, the transition from particle boundary conditions on the membrane displacement and its normal derivative (angle condition) to curvature constraints turns out to be equivalent to simple averaging of boundary conditions (see, e.g., [18] for another motivation in terms of finite differences). Further transition to point-particle models amounts to approximations of mean values by point values. Averaged curve constraint models no longer represent tilt (cf. Remark 4.2), but still preserve some information on the shape of particles which is no longer present in point-particle models. Averaged hybrid models could as well be regarded as regularizations of, in the first instance ill-posed, point-particle models. The size of particles then acts as a regularization parameter which, in contrast to previous approaches (cf., e.g., [53,54]), is physically meaningful and independent of discretization.
6. Point curvature constraints 6.1. Point approximation of mean values. When the particles have a small diameter with respect to the diameter of the domain then it is of interest to consider modelling the particles as points. One approach to obtaining such models is to replace integrals by point evaluations. That is, the mean value ffl Bi ∆u dx is naturally approximated by ∆u(X i ) by sending the diameter of the particle B i to zero. This may be understood in a different way as approximating the integrals in the averages (5.1) by a first-order Gauss formula.
denotes the Dirac functional. In order to state a well-posed minimization problem on a smaller solution space of sufficiently regular functions, we augment the Canham-Helfrich energy J defined in (3.7) by additional higher order terms to obtain with some given regularization parameters κ 8 , κ 6 > 0. This artificial extension could be replaced by a more realistic fourth-order expansion of the bending energy with respect to principal curvatures [38,52]. The strict positivity of κ 8 guarantees that functions which have bounded energy lie in H 4 (Ω). In turn, the continuous embedding H 4 (Ω) ⊂ C 2 (Ω) implies that the G i = δ Xi (∆ · ) are bounded linear functionals on H 4 (Ω). Note that the functionals G i are linearly independent for distinct locations X i , i = 1, . . . , N . This also holds for point values of second-order derivatives δ Xi (∂ xx · ), δ Xi (∂ xy · ), and δ Xi (∂ yy · ).
Differentiation ofJ yields the associated bilinear form The higher order terms inJ give rise to additional boundary conditions defining a suitable closed solution spaceṼ ⊂ H 4 (Ω). For example, we might choosẽ is Ω-periodic and´∂ Ω v ds = 0}.
For the final two cases we consider only rectangular domains Ω. Each choice for V also provides complementary natural boundary conditions for solutions to variational problems posed in that space. Observe thatã( · , · ) is bounded and symmetric onṼ . It is also coercive for any κ 8 > 0 and κ 6 , κ, σ ≥ 0, we refer to [32,35] for a proof. For a more detailed discussion of boundary conditions see [35].
6.3. Fixed locations of particles. We consider the following version of Problem 5.2 with hard point constraints.
Proposition 6.1. There exists a unique solution u ∈Ṽ to Problem 6.1.
Proof. The particles B i are disjoint, so that we have X i = X j , for i = j = 1, . . . , N . Hence, the functionals G i are linearly independent. Now the assertion follows from Proposition 9.7 in the appendix.
Possible anisotropies can be represented by the geometry of particles B i and boundary conditions in Problem 4.1. These are lost completely in the approximation by point mean curvature constraints. Accounting for anisotropies we now prescribe different curvatures Problem 6.2 (Point curvature constraints).
Find u ∈Ṽ minimizing the energyJ (u) onṼ subject to the constraints The functionals G i,j defined in (6.5) are linearly independent for distinct locations X i . Hence, existence and uniqueness again follows from Proposition 9.7. Proposition 6.2. There exists a unique solution u ∈Ṽ to Problem 6.2.
Note that existence and uniqueness of solutions φ i,j of (6.7) follows from the Lax-Milgram lemma, becauseã(·, ·) is bounded and coercive and the linear functionals G i,j are bounded onṼ ⊂ H 4 (Ω).
Proof. The assertion follows from the abstract Proposition 9.7 as applied to the special case = 0 and thus ϕ = 0.
Remark 6.1. Note that Proposition 9.7 also provides a corresponding representation of the solution of Problem 5.2.
We now consider soft curvature constraints. To this end we augment the energỹ J by the penalty term with some small penalty parameter ε > 0 and the Frobenius norm · R N ×3 .
While existence and uniqueness follows from the Lax-Milgram lemma, Proposition 9.3 implies convergence to the hard-constrained solution. 6.4. Varying the locations of particles. We now seek a global minimizer over prescribed curvatures in the sense that we allow for varying locations X = (X i ) of particles. To emphasize that the functionals G = (G i,j ) defined in (6.5) depend on the locations X i , we introduce the notation (6.10) First we consider hard-wall constraints, i.e., we restrict the locations X to ω defined in (4.16) so that the particles would not overlap.
Problem 6.4 (Point curvature constraints with varying locations). Find (u, X) ∈Ṽ × ω minimizing the energyJ (u) onṼ subject to the constraint that there is a X = (X i ) ∈ ω such that G X u = r holds with given r ∈ R N ×3 .
Lemma 6.1. For any X ∈ ω the family (G X,i,j ) ∈ (Ṽ ) N ×3 of functionals G X,i,j ∈ V is linearly independent and G X :Ṽ → R N ×3 is surjective.
Proof. Let X ∈ ω. First we note that surjectivity of G X :Ṽ → R N ×3 is equivalent to linear independence of the family (G X,i,j ). By the assumption 0 ∈ B 0 i we have X i ∈ B i (X) and hence, by (4.16), X i = X j for i = j. Hence we can construct smooth functions v with G X v = r for any r ∈ R N ×3 . Proposition 6.5. Assume that ω = ∅. Then there exists a solution (u, X) ∈Ṽ × ω to Problem 6.4.
Proof. In order to apply Proposition 9.4, we first note that, by Lemma 6.1, for any Y ∈ ω there is a v ∈Ṽ with G Y v = r, i.e., the feasible set is non-empty.
We now provide a reformulation of Problem 6.4 in terms of suitable Green's functions.
Proposition 6.6. Assume that ω = ∅. For given Y ∈ ω, let the Green's functions φ Y,i,j ∈Ṽ and the matrix A Y =ã(φ Y,i,j , φ Y,k,l ) ∈ R (N ×3)×(N ×3) be defined as in (6.7) and Proposition 6.3, respectively. Then each solution (u, X) of Problem 6.4 has the representation Proof. Note that Lemma 6.1 implies that the family (G Y,i,j ) is linearly independent for all Y ∈ ω. Hence Proposition 9.8 can be applied for the special case = 0 (and thus φ 0 = 0). Now we consider soft-wall constraints by augmenting the energy functional with the term V soft (X).

Problem 6.5 (Point curvature constraints with varying locations and soft-wall constraints).
Find (u, X) ∈Ṽ × int ω minimizing the energyJ (u) + V soft (X) subject to the constraint G X u = r for given r ∈ R N ×3 . Proposition 6.7. Assume that int ω = ∅. Then there exists a solution (u, X) ∈ V × int ω to Problem 6.5.
Proof. Noting that V soft (Y ) < ∞ for all Y ∈ int ω, we can prove existence of solutions (u, X) ∈Ṽ × ω as in the proof of Proposition 6.5. Then X ∈ int ω follows from the definition of V soft . Proposition 6.8. Assume that ω = ∅. Then for each solution (u, X) of Problem 6.5 we can represent u as in Proposition 6.6 where where X is now a minimizer of the mapping Proof. The proof can be carried out literally as for Proposition 6.6. Remark 6.3. The results stated in the Propositions 6.5, 6.6, 6.7, and 6.8 also hold for G X = (G X,i ) ∈ (Ṽ ) N with functionals G X,i = δ Xi (∆ · ) defined according to (6.2). Remark 6.4. A similar representation to (6.11) can be derived for the solutions of the average constrained Problems 5.3 and 5.4 with hard-and soft-wall constraints.

Unbounded domains.
In the preceding sections, we have focussed on problems with particles on a membrane that is parametrized over a bounded domain Ω ⊂ R 2 . However, our variational approach is not limited to this case. As an example, we will now consider a physical model as suggested by Bartolo and Fournier [7] with an unbounded membrane parametrized over R 2 . We will formulate this model in terms of our variational framework and then use our general theory to recover some, but not all results that were obtained for bounded domains. In order to facilitate the transfer between [7] and our presentation, we will mostly adopt new notation to [7]. 6.5.1. Hard and soft point constraints in R 2 . Following Bartolo and Fournier [7], we consider an extension of the energy functionalJ defined in (6.3) by an additional confining potential γu 2 with a given constant γ > 0 (there is no danger of confusion with the height γ used elsewhere in this paper). The coupling of the membrane with N pointwise isotropic particles at pairwise distinct locations X i ∈ R 2 is represented by the interaction term the Dirac functionals G X,i,j defined in (6.10), i = 1, . . . , N , and given data C = (C i,j ) ∈ R N ×3 , Γ ≥ 0, and > −1/2, such that N is symmetric positive definite. Now the model considered by Bartolo and Fournier [7] reads as follows. Problem 6.6 (Soft point curvature constraints in R 2 ). Find u ∈ H 4 (R 2 ) minimizing the energy Proposition 6.9. There exists a unique solution to Problem 6.6.
Proof. First we note that the bilinear form associated with the energy functional F m is bounded on H 4 (R 2 ). By partial integration and the fact that the C ∞ -functions with compact support are dense a m ( · , · ) is also coercive on this space. Furthermore the Sobolev embedding H 4 (R 2 ) → C 2 (R 2 ) (see, e.g., [1, Theorem 4.12]), guarantees continuity of each G X,i,j . Since N is symmetric and positive definite for > −1/2 this implies that the bilinear form associated with the energy functional in (6.15) is symmetric and H 4 (R 2 )-elliptic. Hence, the assertion follows from the Lax-Milgram lemma.
In order to identify the hard constrained version of Problem 6.6, we reformulate the interaction term defined in (6.13) according to For increasing penalty parameter Γ → ∞, we therefore obtain the following hard-constrained problem. Problem 6.7 (Point curvature constraints in R 2 ). Find u ∈ H 4 (R 2 ) minimizing the energy F m (u) subject to the constraints Proposition 6.10. There exists a unique solution u ∈ H 4 (R 2 ) to Problem 6.7. Moreover, denoting the solution of Problem 6.6 for fixed Γ ≥ 0 by u Γ , we have Proof. First recall that the bilinear formã(u, v) + γ´R 2 uv dx is bounded, symmetric, and coercive on H 4 (R 2 ) and that the operator G X is bounded. Furthermore by Lemma 6.1 the functionals G X,i,j are linearly independent for pairwise distinct locations X i . Hence existence and uniqueness follows from Proposition 9.7. Finally, since Q is regular, the constraint (6.18) is equivalent to Q(G X,i u − C i ) = 0, i = 1, . . . , N such that the convergence (6.19) is a consequence of Proposition 9.3.
We now derive a representation of the stationary energy F m (u) of Problem 6.7 which will finally turn out to be a hard constrained limit version of the stationary energy obtained by Bartolo and Fournier [7]. Since we would like to emphasize the dependence of this energy on the locations X = (X i ), X i ∈ R 2 , we will denote the solution of Problem 6.7 for fixed X by u X from now on. Then Proposition 9.7 provides the representation with Green's functions φ X,i,j and the matrix A X ∈ R (N ×3)×(N ×3) given by and A X = (a m (φ X,i,j , φ X,k,l )), respectively. Lemma 9.3 implies that the energy at the minimizer is given by We conclude this section with an explicit representation of the entries of A X as appearing in [7]. To this end let G ∈ H 4 (R 2 ) denote the Green's function given by ∀v ∈ H 4 (R 2 ) (6.22) and define the differential operators ∂ (1) = ∂ xx , ∂ (2) = ∂ xy , and ∂ (3) = ∂ yy . Then G can be related to the Green's functions φ X,i,j in the following way. Proposition 6.11. For j = 1, 2, 3 we have ∂ (j) G( · − X i ) = φ X,i,j .
Proof. Let w ∈ C ∞ (R 2 ) and set v = ∂ (j) w( · + X i ) ∈ C ∞ (R 2 ). Inserting this in (6.22) gives By the regularity result given in Lemma 9.4 in the appendix we have G ∈ H 6 (R 2 ). Hence we have ∂ (j) G ∈ H 2 (R 2 ) and partial integration and translation invariance of integrals over R 2 yields Since C ∞ (R 2 ) is dense in H 4 (R 2 ) we have shown that ∂ (j) G( · − X i ) coincides with the unique solution φ X,i,j of (6.20).

6.5.2.
Varying the location of particles. In analogy to Problem 6.4 with locations of particles varying in the compact set Ω, we now consider varying locations in R 2 . We first fix some notation for the R 2 analogue of the set ω ⊂ Ω N .
Recall G X = (G X,i,j ) with G X,i,j defined in (6.10) and let C = (C i,j ) ∈ R N ×3 be given.

Problem 6.8 (Point curvature constraints with varying locations in R 2 ).
Find (u, X) ∈ H 4 (R 2 ) ×ω such that u minimizes the energy F m (u) subject to the constraint Note that Proposition 9.4 can be no longer applied to prove existence, becausẽ ω is not compact. While we could still show thatω is closed, it fails to be bounded here. Furthermore, we should not expect a solution to Problem 6.8 to exist. To see this, consider the case of N = 2 particles at varying locations X = (X 1 , X 2 ) ∈ω and choose C = (C 1 , C 2 ) in such a way that their interaction energy decreases with separation. As an example, one might chose identical isotropic particles as considered, e.g., by Bartolo and Fournier [7]. Let u X be the unique solution of the corresponding Problem 6.7 with fixed locations X. Then F m (u X ) is precisely the interaction energy and is well-known to strictly decrease as the points X 1 and X 2 are moved apart. As the distance between X 1 and X 2 can become arbitrary large, there can be no minimal set of locations and thus no solution of the corresponding Problem 6.8. 6.6. Discussion. The constrained minimization Problems 6.1 and 6.2 stem from the models discussed, e.g., in [7,17,18,40,48,53,55,72]. See also [30,58]. The addition of the higher order terms to the energy functional (see (6.3)) to ensure well posedness is done in [7]. Prior to this the ill posedness is dealt with by only studying large separation distances between particles [48] or by truncating the Fourier expansion of the solution, termed the high wave-vector cutoff in [18]. When put in our framework, these papers study an energy functional of the form given in equation (6.15).
The hard inclusions limit as described in [7] corresponds to the limit Γ → ∞, here we are able to rigorously understand this limit. By what we have established in Proposition 6.10, the hard inclusions limit is equivalent to the quadratic minimization Problem 6.7.
This limit problem with anisotropic particles prescribing curvatures as in (6.5) is studied in [18]. The elastic interaction energy is calculated and a thermal equilibrium is approximated using a Monte-Carlo algorithm. In the equilibrium configuration proteins aggregate into one region and form an egg carton type structure with the anisotropic particles located at the saddle points of this structure. This equilibrium is analogous to the global minimizer in Problem 6.8.

Point value constraints and point forces
7.1. Particles interacting with the cytoskeleton. In the preceding sections we have studied a variety of models describing the interaction of lipid membranes with embedded particles, scaffolds or wrapped particles (cf. Figure 4.1). We now consider interactions of the membrane with thin actin filaments that are anchored to the cytoskeleton and may prescribe displacements of the membrane or with particles that apply forces that may be due to either entropic or direct chemical interactions (see, e.g., [26,31]). In the mathematical models to be considered in this section, these effects are represented by point value constraints or point forces.

Point value constraints at fixed locations.
Prescribed point values at given locations X = (X i ) ∈ Ω N are represented by the constraints as illustrated in Figure 7.1. Note that δ Xi ∈ V and thus F X ∈ (V ) N due to the continuous embedding V ⊂ H 2 (Ω) ⊂ C(Ω). Hence, from a mathematical point of view, minimization problems with point value constraints (7.1) share their basic properties with minimization problems with point curvature constraints (6.1) as considered above. We first consider prescribed point values at fixed locations.

Problem 7.1 (Point value constraints).
Find u ∈ V minimising the energy J on V subject to the constraints (7.1). In order to avoid possible conflicts of point constraints (7.1) with boundary conditions (3.8), we exclude X i ∈ ∂Ω.
Proposition 7.1. For distinct locations X 1 , . . . , X N ⊂ Ω there exists a unique solution u ∈ V to Problem 7.1.
Proof. As the locations X 1 , . . . , X N are distinct and contained in Ω the functionals δ Xi are linearly independent for all three choices (3.8) of V . Hence, the assertion follows from Proposition 9.7.
Remark 7.1. In applying Proposition 9.7, we can also derive a representation of the solution in terms of Green's functions φ i ∈ V , defined by This can be viewed as finding the optimal locations for filaments which prescribe a particular set of displacements.

Problem 7.2 (Point value constraints with varying locations).
Find (u, X) ∈ V × Ω N such that u is minimising the energy J on V subject to the constraint F X u = α with given α ∈ R N . Lemma 7.1. The mapping Ω X i → δ Xi ∈ (V ) and thus the mapping Ω N X → F X ∈ (V ) N is continuous.
Proof. Since the Sobolev embedding theorem provides continuity of the injection V ⊂ H 2 (Ω) → C 0,λ (Ω) for any Hölder-exponent 0 < λ < 1 (see, e.g., [1, Theorem 4.12]), the estimate Proof. It is well-known that Ω N ∈ R N is compact. Furthermore, the mapping Now the assertion follows from Proposition 9.4. Proof. The solution of Problem 7.2 is equivalent to solve and to take X ∈ Ω N such that F X u = α. Hence, it is sufficient to show that To this end let v ∈ V α1,1 and X N ∈ Ω such that v(X N ) = α N . Then, by continuity of v on Ω, for all three choices (3.8) of V there is X 0 ∈ Ω with v(X 0 ) = 0. Now, by continuity of v and convexity of Ω, the intermediate value theorem implies that v attains each value α i ∈ co{0, α N } at some point X i ∈ Ω and hence v ∈ V α,N .
We now move on to the case where the prescribed point values α i do not have the same sign. Similarly to the previous case, the behaviour is governed by the extreme values of α, in this case the greatest and least. We have thus reduced N to one or two constraints, based upon the signs of the prescribed point values. We first concentrate on the case of point values with the same sign and reformulate Problem 7.2 in terms of the Green's function φ x ∈ V , defined by By definition, a(φ x , φ x ) = φ x (x) holds for all x ∈ Ω. In order to exclude the degenerate case φ x = 0 we will constrain x to the set Notice that Ω V depends on the choices (3.8) of the boundary conditions incorporated in V : For V = H 2 0 (Ω) and V = H 2 (Ω) ∩ H 1 0 (Ω) we have Ω V = Ω whereas V = H 2 p,0 (Ω) allows for Ω V = Ω. Proposition 7.5. Assume that the prescribed point values have the same sign and let 0 ≤ |α 1 | ≤ · · · ≤ |α N |. Then the solution of Problem 7.2 is given by Proof. First we note that (7.4) is well-defined because we have φ x (x) = 0 for x ∈ Ω V . Proposition 7.3 implies that the solution u of Problem 7.2 is the minimizer of J subject to the constraint that δ X N u = α N holds with some X N ∈ Ω. For α N = 0 we only have the trivial minimizer u = 0 which is in accordance with (7.4).
Hence we must have X N ∈ Ω V for all solutions (u, X N ) and the representation (7.4) follows from Proposition 9.8 for N = 1, M = Ω, M = Ω V , and T 1 (y) = δ y for y ∈ M.
Note that in the generic case α N = 0, the optimal location X N is independent of α N and u depends linearly on α N .
If the prescribed point values do not have the same sign, then Problem 7.2 can be reformulated in terms of two Green's functions φ Y1 , φ Y2 ∈ V defined by and the associated Gramian matrix A Y = (a(φ Yj , φ Yi )) = (φ Yi (Y j )).
If the assumption α 1 < 0 < α N is not fulfilled, then all α i have the same sign. In this case we can use the representation given by Proposition 7.5. 7.4. Point forces at fixed locations. We now consider forces exerted to the membrane that are localized to certain points X i ∈ Ω, i = 1, . . . , N . These forces are perpendicular to Ω with positive or negative direction and give rise to the additional term in the energy functional to be minimized. Here, β i ∈ R \ {0} are given constants representing the magnitude of point forces at the locations X i . We set with the closed subspace V ⊂ H 2 (Ω) defined in (3.8) and consider the following minimization problem.

Problem 7.3 (Point forces at fixed locations).
For given X = (X i ) ∈ Ω N find u ∈ V minimising the energy E(u, X) on V .
Existence and uniqueness of a solution u ∈ V of Problem 7.3 follows from the Lax-Milgram lemma. It is characterized by the variational equality The solution can be represented by Green's functions φ x as defined in (7.3).
Lemma 7.2. For given X ∈ Ω N the solution u ∈ V of Problem 7.3 is given by Proof. The assertion follows directly from the linear representation (7.7) of X by the functionals δ Xi . 7.5. Point forces at varying locations.
7.5.1. Existence of global minimizers. We now seek a global minimizer over prescribed point forces, in the sense that we allow the point forces to be applied at varying locations X = (X i ) ∈ Ω N .
Problem 7.4 (Point forces at varying locations).
Proof. In light of the continuity of Ω X i → δ Xi ∈ V as stated in Lemma 7.1, the assertion follows from Proposition 9.6.
Using the representation for fixed X given in Lemma 7.2, we will also construct a representation of solutions to Problem 7.4. To this end, we from now on denote by u Y ∈ V the unique minimizer of E( · , Y ) for given Y ∈ Ω N . As a first step, we compute the energy of such minimizers. Lemma 7.3. Let Y ∈ Ω N be given and A Y = (a(φ Yi , φ Yj )) ∈ R N ×N . Then Proof. After inserting v = u Y into the variational equality (7.8) for u Y , we use the definition of Y , and the representation of u Y as given in Lemma 7.2 to obtain . This completes the proof.
As a direct consequence we get Clustering. Having established the existence of global minimizers we will now explore the properties of minimizers for particular combinations of the parameters β i , i = 1, . . . , N . Of particular interest will be exhibiting cases, where optimal locations of point forces lie inside Ω. Of course, this is of no interest under periodic boundary conditions. As such, we will assume V = H 2 0 (Ω) or V = H 2 (Ω) ∩ H 1 0 (Ω) for the rest of this section, but will remark where this assumption plays a role. We will also show clustering behaviour for larger numbers of point forces and that opposite point forces do not annihilate each other.
We first show that point forces do not cluster on the boundary ∂Ω of Ω. Proof. Assume that (u, X) ∈ V × (∂Ω) N solves Problem 7.4. Then X (v) = 0 holds for all v ∈ V and therefore u = u X = 0. Hence, we have . . , N . This contradicts optimality of (u, X).
The following lemma quantifies the change of energy that is caused by moving a single point force. This is the key ingredient to prove clustering of point forces later on.
The representation of energy in Lemma 7.3 and the binomial formula provide the estimate . Inserting these identities into (7.10) we obtain the assertion.
In the forthcoming clustering analysis we will make use of the equivalence relation Recall that we have δ x = 0 on V for all x ∈ ∂Ω. Hence, x = y holds, if and only if x = y or x, y ∈ ∂Ω. By definition, the locations of a solution can be replaced by equivalent locations, i.e., if (u, X) solves Problem 7.4 and Y i = X i holds for all i = 1, . . . , N , then (u, Y ) does also solve Problem 7.4. Now we are ready to prove clustering of point forces.
Proposition 7.9. Assume that (u, X) ∈ V × Ω N is a solution to Problem 7.4.
Proof. Let (u X , X) ∈ V × Ω be a solution of Problem 7.4. Then we have Recall that Lemma 7.4 implies X / ∈ (∂Ω) N . Hence, at least one point force must be located in Ω. Let X j ∈ Ω be arbitrarily chosen. Then it is sufficient to show that X i = X j must hold for all i ∈ I j = {l = 1, . . . , N | sgn(β l ) = sgn(β j )}.
Without loss of generality assume that β j > 0 and that j is selected such that δ Xj (u X ) ≥ δ Xi (u X ) holds for all other X i ∈ Ω with i ∈ I j . Then the same estimate is valid for all i ∈ I j , because (7.13) yields δ Xj (u X ) ≥ 0 = δ Xi (u X ) for all X i ∈ ∂Ω.
In contradiction to the assertion, we now assume that X k = X j holds for some k ∈ I j . Application of Lemma 7.5 with Y i = X i , i = k and Y k = X j , together with . Now we have either X k ∈ ∂Ω, and therefore φ X k = 0, or X k ∈ Ω, and therefore that φ X k and φ Xj are linearly independent. In both cases a(φ Xj − φ X k , φ Xj − φ X k ) > 0 providing E(u Y , Y ) < E(u X , X). This contradicts optimality of (u X , X). Remark 7.4. As a consequence of Proposition 7.9, the point forces with positive (negative) sign either cluster in one point X + ∈ Ω (X − ∈ Ω) or are all located on the boundary ∂Ω. By Lemma 7.4, not all N point forces can be located at the boundary. Hence, point forces of a solution (u, X) of Problem 7.4 are clustering in exactly one of the following three ways.
(i) X i = X + ∈ Ω for all i with β i > 0 and X i = X − ∈ Ω for all i with β i < 0, (ii) X i = X + ∈ Ω for all i with β i > 0 and X i ∈ ∂Ω for all i with β i < 0, (iii) X i = X − ∈ Ω for all i with β i < 0 and X i ∈ ∂Ω for all i with β i > 0. We may regard the occurrence of one of these three cases as a property of a( · , · ), the parameters β i , and Ω.
As another consequence of Proposition 7.9 we can characterize the solutions to Problem 7.4 with N forces in terms of an equivalent problem with at most two forces.
Then (u, X) ∈ V × Ω N is a solution of Problem 7.4 with (X + , X − ) ∈ Ω 2 satisfying (7.12), if and only if (u, (X + , X − )) ∈ V × Ω 2 is a solution of Problem 7.4 with point forces Proof. Proposition 7.9 implies that the locations X of all solutions to Problem 7.4 are contained in the subset Hence minimizing E(u, X) over V × Ω N is equivalent to minimization over V × M .
By definition of M , we can identify M with Ω 2 by the condition (7.12) up to componentwise equivalence in the sense of (7.11). Now, let X ∈ M be identified with X 0 = (X + , X − ) ∈ Ω 2 in this way. As a consequence of (7.12), we then have X = X0 and thus u X = u X0 . In light of Lemma 7.3, this leads to Therefore, minimization of E(u, X) over V × M is equivalent to minimization of the energy E 0 (u X0 , X 0 ) over V × Ω 2 . This concludes the proof. By Proposition 7.9 at least one of the clustering points X + , X − ∈ Ω must be contained in Ω. Utilizing the values of β + and β − defined in (7.14), we can often exclude one of the three cases in Remark 7.4. Proposition 7.10. Let (u, X) be a solution of Problem 7.4. If |β + | > |β − |, then X + ∈ Ω and if |β + | < |β − |, then X − ∈ Ω.
We now assume that all forces point in the same direction. In this case, the solutions of Problem 7.4 can be obtained by solving Problem 7.4 with a single force.
Corollary 7.2. Assume that all of the coefficients β i have the same sign. Then (u, X) ∈ V × Ω N is a solution to Problem 7.4, if and only if X 1 = · · · = X N ∈ Ω and (u, X 1 ) is a solution of Problem 7.4 with one point force Proof. By Lemma 7.4 there must be at least one X k ∈ Ω. Then, Proposition 7.9 provides X 1 = · · · = X k = · · · = X N and the assertion follows from Corollary 7.1.
We have thus characterised the behaviour of systems with forces that are all pointing in one direction: The global minimizer is simply when all of the particles lie at the same point and that point is a global minimizer for only one point force. There is still no uniqueness however, as global minimizers for the one point force problem are not unique in general. The uniqueness for the one point force problem may be regarded as a property of the domain Ω.
Remark 7.5. All results given above can be extended to the case V = H 2 p,0 (Ω) by replacing all occurrences of Ω by Ω and dropping all cases where ∂Ω shows up. 7.6. Discussion. The models formulated in Problem 7.1 and Problem 7.2 describe the optimal shape of a membrane under point constraints and the optimal location of such constraints. This approach could be used to describe the action of actin filaments bound to the membrane. Such kind of problems also occur in the study of thin plates. For example, Problem 7.1 is the central object in the study of thin plate splines and Problem 7.2 is analysed in [11] which studies support points of a plate, producing this problem with homogeneous data α = 0.
The model set out in Problem 7.3 and further extended and analysed in Section 7.5 is motivated by the general approach in [26] where protein membrane interaction is described by an additional term in the membrane energy functional representing the work done by the pressure exerted by proteins. In [26] the particles are assumed to have a positive diameter and are bound to membrane. We have adapted this model to particles anchored to the cytoskeleton applying point forces. Note that the results on clustering of point forces, derived above, do not agree with the interaction of finite sized particles, as investigated in [26]. The key difference between the two models is that the point forces in Problem 7.4 do apply a net force to the membrane which is not the case for the interactions studied in [26]. The action of protrusive forces on a membrane is discussed in [31,69]. The variational framework we have introduced may be also applied in this case and used to analyse the membrane mediated interactions between particles.

Numerical experiments
In this section, we present various numerical computations with hybrid models as introduced above. The variational setting of these models naturally suggests discrete finite element counterparts which will be specified and analysed in more detail in future publications. We first concentrate on a comparison of finite-size and point-like particle descriptions (cf. Section 4 and 6). In particular, we investigate the effect of the arrangement and the separation distance of particles on the membrane-mediated interaction energy. Later, we complement our theoretical investigations of clustering of point forces (cf. Section 7) by numerical computations. Our numerical results are discussed in comparison with qualitative and quantitative results in the existing literature (see, e.g., [7,26,30,40,71]). 8.1. Finite size particles and point curvature constraints. We consider a membrane decorated with N = 2 particles with circular or elliptical cross-sections and equal or opposite orientations. Here, we focus on two different types of models, treating particles as rigid objects with finite size (cf. Section 4) or as described by point curvature constraints (cf. Section 6), respectively. 8.1.1. Finite size particles. We consider a parametrized curve constrained problem on V = H 2 0 (Ω) (cf. Section 3.3) that involves both varying height and tilt angle of the particle as described in Remark 4.2. However, the coupling to the membrane is now performed by boundary conditions on ∂Ω B \ ∂Ω instead of curve constraints on Γ i ⊂ Ω (cf. Section 4). More precisely, we aim to approximate (u, γ, α) ∈ V Ω B × R 2 × R 2·2 , minimizing the energy J (u) subject to the boundary conditions (4.15) with unknown height γ = (γ 1 , γ 2 ) ∈ R 2 and tilt α = (α 1 , α 2 ) ∈ R 2·2 of particles. In the elliptical case, the particle rotation is an additional parameter that is fixed in each computation.
The discretization is performed by (non-conforming) Morley finite elements (see [13], [45]). Straightforward (elementwise) application of the bilinear form (3.10) to the Morley finite element space yields an indefinite stiffness matrix. We therefore utilize the following equivalent reformulation of (3.10), depending on a parameter c ∈ (0, κ] where · , · ∂Ω B stands for the dual pairing of H − 1 2 (∂Ω B ) and H 1 2 (∂Ω B ). For sufficiently smooth functions · , · ∂Ω B is just the L 2 -inner product on ∂Ω B . For details on the numerical treatment of this boundary term, we refer to [75]. We use the parameter value c = κ in all subsequent computations.
We compute approximate displacements u of a membrane over the domain Ω = (−3, 3) 2 , generated by two particles that are either circular with radius r = 0.1 or elliptical with major and minor half-axis a = 0.14 and b = 0.06. Note that the unit of length is 10r, i.e., length is measured in terms of the particle size. We choose height functions h 1 = h 2 = 0, and normal derivatives with equal s 1 = s 2 = 1 (left) or opposite orientation s 1 = −s 2 = 1 (right), respectively. We select the material parameter κ = 1 and also fix the unit of energy in this way. Figures 8.1 and 8.2 show the clipping of the results to (−1, 1) 2 for σ = 0. The elliptical particles are arranged collinearly ( ) for equal orientation (cf. Figure 8.2a) or parallel ( ) for opposite orientation (cf. Figure 8.2b). As illustrated in Figure 8.3, these configurations seem to be energetically optimal, in the sense of having lower energies than other configurations. In these figures, the approximate interaction potential, i.e., the approximate minimal energy J (u), of collinear and parallel arrangements is plotted over the separation distance R between the two particles, here defined as the    shortest distance between the particle's boundaries. While collinear particles are energetically preferred for equal orientation (cf. Figure 8.3a), parallel arrangement is preferred for oppositely oriented particles (cf. Figure 8.3b). For other arrangements we found energy curves lying in between these two extremal configurations. We depict the interaction potential for an orthogonal arrangement ( ) as an example. These numerical results confirm related theoretical considerations in [41].  We now study the membrane-mediated interaction potential as a function of the separation distance R between two particles which are arranged as in Figure 8.1 and 8.2. For ease of comparison with point-like particle descriptions (cf. Subsection 8.1.2), the separation distance R from now on stands for the distance between the midpoints of particles. In order to reduce the possible influence of the (artificial) boundary ∂Ω, we select separation distances R ≤ 1. Figure 8.4 shows the interaction potential over the separation distance R for two particles with circular cross-section and various values of σ. We choose σ = 0 and σ = 1, 4, 9, 16, 25 which, for a particle size of 3nm, corresponds to usual interaction lengths κ/σ ranging from 6nm to 30 nm [26]. We say that the interaction is repulsive or attractive if the derivative of the energy with respect to the separation distance is negative or positive, respectively. In the case of equally oriented particles, i.e., for contact angles s 1 = s 2 = 1, we observe repulsion at all separation distances R and for all values of σ. (cf. Figure 8.4a). In contrast, for oppositely oriented particles, i.e., for contact angels s 1 = −s 2 = 1, the interaction is repulsive for small and attractive for larger distances R provided that σ > 0 (cf. Figure 8.1b). The preferred separation distance R * > 0, i.e., the distance resulting in the lowest energy, decreases with increasing membrane tension. In both cases, the interaction potential becomes larger for growing σ. For vanishing membrane tension, the interaction is always repulsive.
In the case of particles with elliptical cross-section, we obtain a fundamentally different picture: Figure 8.5 shows the interaction potential over the separation distance for various elliptical shapes in the case σ = 0. Both for equal (left) and opposite orientation (right), the interaction is repulsive for small and attractive for larger separation distances. The resulting preferred separation distance R * > 0 decreases with increasing ratio a/b, i.e. for slimming particles. Increasing membrane tension σ also yields a decrease in preferred distance R * and an increase in interaction potential, similar to the circular case.  The numerical approximation utilizes a splitting of the eighth order Problem 6.3 into an equivalent system of three second order equations for the unknown functions u, w = ∆u and z = ∆w. For this purpose, we formally rewrite the energyJ (u) defined in (6.3) in terms of w: Note that the corresponding Euler-Lagrange equation is fourth order in w. We impose essential boundary conditions u = 0, w = ∆u = 0, and z = ∆ 2 u = 0 on ∂Ω. For fixed p ∈ (2, ∞) and q ∈ (1, 2) such that 1/p + 1/q = 1, we consider the variational problem to find (u, w, z) ∈ H 1 0 (Ω) × W 1,p 0 (Ω) × W 1,q 0 (Ω) such that the three equations (8.2)ˆΩ −κ 8 ∇z · ∇v + κ 6 ∇w · ∇v + κwv − σuv dx hold for all v ∈ W 1,p 0 (Ω). Regularity of solutions to Problem 6.3 implies that u ∈Ṽ solves Problem 6.3, if and only if (u, ∆u, ∆ 2 u) solves the system (8.2). For details we refer to [35]. The system (8.2) is finally discretized by P 1 finite elements. The penalty parameter is taken as ε = 1 × 10 −8 in our computations. For the material parameters we use κ 8 = 1.23 × 10 −6 and κ 6 = 1.11 × 10 −3 or, equivalently, fix the ratios (κ 6 /κ) 1/2 = (κ 8 /κ) 1/4 = 1/30. Selecting κ = 1 and the same physical length scale as in the   previous finite size particle case, this means that both ratios correspond to 1nm which is of order of the thickness of the membrane as suggested in [7, section 2]. We also choose the computational domain to be Ω = (−3, 3) 2 . The values r k = ±20 for the point constraints associated with circular particles of radius 0.1 are obtained from the approximation (6.1). Figures 8.6a and 8.6b plot the clipping to (−1, 1) 2 of the approximate membrane displacement obtained for equal and opposite curvature constraints, respectively, which correspond to equal and opposite orientations of particles. Investigating their interaction potential for various values of membrane tension σ in analogy to the previous section, Figure 8.7a shows that equally oriented particles repel each other for separation distances R > 0.2. The strength depends upon membrane tension σ. For distances R < 0.2 we observe attraction. For oppositely oriented particles the interaction is repulsive for small and attractive for larger separations as depicted in Figure 8.7b. The strength of the attraction increases with σ.

Discussion.
For the point inclusion model, separation distances shorter than R = 0.2 are physically impossible, as they would represent situations where the particles with radius r = 0.1 would overlap. This explains and contextualises the observed unphysical interactions for short distances R < 0.2. For comparison with the results obtained for the finite-sized circular particles, it is therefore sufficient to consider separation distances R ≥ 0.2. The major difference is that the point model does not reproduce the strong repulsion at the limiting separation R = 0.2. This is unsurprising as the point model introduces a cutoff length similar to the length the strong repulsion acts over. For separation distances above this cutoff i.e., for R > 0.3, the observed interactions agree qualitatively and quantitatively as depicted in Figures 8.7 and 8.4. A point inclusion counterpart to the finite-sized model with elliptical particles would require anisotropic functionals of the form (6.5) which are not considered here.
The interaction potential between two particles with circular cross-section has been intensively studied for more than 20 years (see, e.g., [7,18,19,30,40,48,58,71]). For finite-size particles Weikl, Kozlov & Helfrich [71] found by analytical considerations that in case of positive membrane tension, i.e., σ > 0, the interaction depends on the relative orientation of the two particles: Equally oriented particles repel each other at all separation distances, whereas for particles with opposite orientation the interaction is repulsive at small and attractive at larger distances. For σ = 0 the interaction is found to be repulsive at all distances independent of the particles' orientation, confirming earlier results for finite-size particles in [30,40,58]. Similar results have been obtained for point-like particles [7,18,48]. These theoretical findings are in accordance with our numerical computations.
Investigations of membrane-mediated interaction of particles with non-circular cross-sections are rare. In [58] and [41] the effect of the particle shape on the character of interaction is studied. The authors of [41] consider a pair of identical particles whose elliptical cross-section is a small perturbation of the unit circle. They obtain that introducing the horizontal orientation of the particle as an additional degree of freedom, qualitatively changes the asymptotic interaction character from a repulsion to an attraction, which is in agreement with a previous result in [58] and with our numerical computations.
Note that well-known interaction laws of the form 1/R 4 for circular and 1/R 2 for non-circular particles are based on large distance asymptotics assuming an unbounded asymptotically flat membrane with particles separated by distances which are large compared to their size, i.e. r R. Although our numerical experiments cover the complementary situation r ≈ R, it turns out that our numerical results quantitatively reproduce the asymptotical results given in [71] even for small separation distances.

Point Forces.
8.2.1. Interaction potential for fixed locations. We consider Problem 7.3 with point forces at fixed locations X = (X i ) ∈ Ω with the solution space V = H 2 (Ω) ∩ H 1 0 (Ω). Similar to Subsection 8.1.2, our numerical approach is based on a splitting of the fourth order problem (7.8) into two second order problems for the unknown functions u and w = κ∆u − σu. For fixed p ∈ (2, ∞) and q ∈ (1, 2) such that 1/p+1/q = 1, we consider the variational problem to find (u, w) ∈ H 1 0 (Ω)×W 1,q 0 (Ω) such that  [35]. We discretize the system (8.3) using P 1 finite elements.
Recall that point forces of equal sign cluster to one point (see Corollary 7.1). Thus, we restrict the investigation of membrane-mediated interactions to the case N = 2 and α 1 < 0 < α 2 . We choose the same physical length scale as in Section 8.1, fix the same scaling of the energy by selecting κ = 1, and consider the domain Ω = {x ∈ R 2 | |x| < 1} for all subsequent computations. We also choose the forces α 1 = −10 and α 2 = 10. For a motivation of this choice, let us consider the length scale κ/|α i | and the typical value κ ≈ 20 κ B T with the Boltzmann constant κ B and absolute temperature T ≈ 300K [26]. Then our selection corresponds to the strength |α i | ≈ 27 pN . This is reasonable as forces applied by actin polymerisation are of the order 10pN [2].
To study the interaction potential between these two opposite forces, we fix X 1 = (0, 0), allow X 2 to vary along the abscissa and compute the resulting approximate minimal energy J as a function of the separation distance R. This is done for a variety of values of σ corresponding to interaction lengths κ/σ ranging from 1 nm to 30 nm [26]. The results depicted in Figure 8.8 show that opposite forces repel each other and that this repulsion depends upon the ratio κ/σ. Increasing σ, i.e. decreasing the ratio κ/σ, yields a decrease in the distance over which the repulsive interaction plays a role. For σ = 1 the repulsion persists close to the boundary, whereas for higher values of σ the repulsion has a shorter length scale and is, from a certain distance R * on, dominated by the inwards force applied at the boundary. Note that this inwards force is a consequence of the (artificial) boundary condition u = 0. 8.2.2. Discussion. The above numerical findings could be related to the theoretical results derived in Section 7.5. According to Remark 7.4, there are essentially two types of global minimizers solving Problem 7.4. A type 1 global minimizer is characterized by X + , X − ∈ Ω (case (i)) and type 2 means that either X + or X − is located on ∂Ω (cases (ii) and (iii)). The numerical results shown in Figure 8.8 indicate that type 1 or type 2 minimizers occur for sufficiently small or large ratio κ/σ, respectively.  κ = 1, σ = 100. Recall that type 2 global minimizers only occur due to influence of the domain boundary. Thus, it is not surprising that there is no such result in the existing literature that mostly concentrates on unbounded asymptotically flat membranes. However, dependence of the length scale of repulsive interaction between particles which apply forces to the membrane upon the ratio κ/σ is well known and discussed, e.g., in [26].
9. Appendix 9.1. Constrained minimization. In this section we collect some results on linearly constrained minimization problems and their approximation by penalization. Throughout this section, J denotes the quadratic functional , v ∈ V, defined on a Hilbert space V with norm · . We assume that the bilinear form a( · , · ) : V × V → R is symmetric and elliptic in the sense that holds with positive constants γ, Γ ∈ R, and ∈ V . We consider minimization of J on the affine subspace V b ⊂ V , given by Here, X 0 ⊂ X is a complete subspace of a pre-Hilbert space X with scalar product ( · , · ) X and norm · X , T : V → X is a bounded linear operator, and b ∈ X. Note that where P : X → X is the orthogonal projection with ker P = X 0 .
It is well-known that Problem 9.1 is equivalent to find u ∈ V b such that where V 0 = {v ∈ V | T (v) ∈ X 0 } is obtained for b = 0. As T is continuous and X 0 is closed, the affine subspace V b is also closed. Hence, the Lax-Milgram lemma provides the following existence and uniqueness result. In order to derive a penalty approximation of the constraints T v ∈ X 0 + b occurring in (9.2), we consider the following minimization problem. Proof. As v ∈ V b ⇔ T v = y + b with some y ∈ X 0 , the sets of admissible v are the same for both minimization problems.
A penalty approximation of Problem 9.1 now reads as follows.

Problem 9.3 (Penalized parametrized minimization problem).
Find (u ε , x ε ) ∈ V × X 0 minimizing the energy The associated variational formulation of Problem 9.3 is given by The orthogonal projection P : X → X with ker P = X 0 was used to characterize the solution space V b according to (9.3). Now we can use the projection P to characterize the solution of Problem 9.3 as follows.
Proposition 9.3. Assume that V b = ∅ and let u, u ε denote the solutions of Problem 9.1 and 9.3, respectively. Then (u ε , x ε ) → (u, x) converges in V × X as ε → 0. Moreover, P T (u ε − u) 2 X ≤ Cε holds with a positive constant C, depending only on γ, Γ, and .
for k → ∞. Since this holds for any weakly convergent subsequence, we have shown u ε → u in V . Now continuity of P and T provides x ε → x in X.
Note, that constraints of the form T u = b are covered by the choice X 0 = {0}. 9.2. Global minimizers. We now consider a generalization of Problem 9.7 by parametrizing the functional T over a suitable set M. More precisely, for a given metric space M, and a given mapping T : M → L(V, X), with L(V, X) denoting the space of all bounded linear mappings from X to V , we define the admissible set W = {(v, y) ∈ V × M | v ∈ V b,y }, V b,y = {v ∈ V | P (T (y)v − b) = 0} (9.11) with P : X → X denoting the orthogonal projection with ker P = X 0 as in the preceding section. Furthermore we allow for an additional term V : M → R ∪ {∞} to model certain preferences of y ∈ M. Note that ξ is well-defined by Proposition 9.1. Proof. We will show that ξ is lower semi-continuous. Then ξ + V : M → R ∪ {∞} is also lower semi-continuous on the compact set M which provides existence of a minimizer x ∈ M [16, (12.7.9)].
Let (x n ) ∈ M with lim n→∞ x n =x ∈ M. We have to show m = lim inf n→∞ ξ(x n ) ≥ ξ(x).
Since this is trivial for m = ∞ we now assume m < ∞. Coercivity of J implies m > −∞. Furthermore, there is a subsequence, still denoted by (x n ), such that lim n→∞ ξ(x n ) = m and ξ(x n ) < ∞ for all n ∈ N. As a consequence of ξ(x n ) < ∞, we have V b,xn = ∅ so that u n = arg min is well-defined. Boundedness of J(u n ) = ξ(x n ) < ∞ and coercivity of J imply that u n is bounded. Hence there is aū ∈ V and another subsequence, still denoted by (u n ), such that u n ū for n → ∞.
Now continuity T (x n ) → T (x) yields P T (x n )u n P T (x)ū. Hence,ū ∈ V b,x and (ū,x) ∈ W. As J is convex and continuous, J must be weakly lower semi-continuous (see, e.g. [23][Chapter I, Corollary 2.2]) which gives This concludes the proof.
In order to show existence of minimizers for a penalized version of Problem 9.4 we will use a general result on parametrized problems without constraints. In the following we will identify continuous bilinear forms a x ( · , · ) : V × V → R with their operator representations a x : L(V, V ). Lemma 9.2. Assume that M is compact, that the functions M y → a y ( · , · ) ∈ L(V, V ), M y → y ( · ) ∈ V are continuous, that a y ( · , · ) is symmetric and uniformly coercive for all y ∈ M, that V : M → R ∪ {∞} is lower semi-continuous, and that there is Proof. For any y ∈ M we define u y = arg min v∈V J y (v). We will show that the function η : M → R defined by First we note that (v, y) → J y (v) is continuous on V ×M, because y → a y ( · , · ) ∈ L(V, V ) and y → y ( · ) ∈ V are continuous. Now let (y n ) ∈ M with lim n→∞ y n = y ∈ M. Then the first lemma of Strang [14, Theorem 4.1.1] implies u y − u yn ≤ (1 + a y ) sup w∈V |a y (u y , w) − a yn (u y , w)| w + sup w∈V | y (w) − yn (w)| w ≤ (1 + a y ) a y − a yn u y + y − yn −→ n→∞ 0.
After these preparations, we consider a penalized version of Problem 9.4.
Now continuity of T implies continuity of y → (P T (y)) * P T (y) ∈ L(V, V ) and y → (P T (y)) * P b ∈ V and thus of y → a y ( · , · ) and y → y ( · ). Furthermore a(v, v, ) ≤ a y (v, v) implies uniform coercivity of a y ( · , · ) with respect to y ∈ M. Hence Lemma 9.2 provides the assertion.
To conclude this subsection we consider a minimization problem with parametrized source term M y → y ( · ) ∈ V . Problem 9.6 (Parametrized source term). Find (u, x) ∈ V × M minimizing the energy J x (u) = 1 2 a(u, u) − x (u) on V × M.
Existence is an immediate consequence of Lemma 9.2 with a y ( · , · ) = a( · , · ) and V = 0. Proposition 9.6. Assume that M is compact and that M y → y ( · ) ∈ V is continuous. Then there is a solution (u, x) ∈ V × M of Problem 9.6. 9.3. Equivalent characterization of minimizers. We first consider Problem 9.1 in the special case (9.14) Problem 9.7 (Finite dimensionally constrained minimization).
Find u ∈ V minimizing the energy J subject to the constraints (9.15) T i u = b i , i = 1, . . . , N.
While Proposition 9.1 provides existence and uniqueness of solution to Problem 9.7, if and only if V b = ∅, we will now derive an equivalent characterization of the solution in terms of suitable basis functions. Utilizing the Lax-Milgram lemma, we define φ 0 ∈ V to be the unique solution to We also define φ i ∈ V , i = 0, . . . , N , such that and the associated Gramian matrix A = (A ij ) ∈ R N ×N with A ij = a(φ i , φ j ).
Proposition 9.7. Assume that Problem 9.7 is non-degenerate in the sense that the functionals T i , i = 1, . . . , N , are linearly independent. Then there is a unique solution of Problem 9.7 given by Proof. Linear independence of the functionals T i , i = 1, . . . , N , implies that T is surjective. Hence, V b = ∅ and Problem 9.7 has a unique solution by Proposition 9.1. Furthermore the functions φ i defined in (9.16) are linearly independent and the associated Gramian matrix A is regular. Let u be given by (9.17). Then u ∈ V b , as for any i = 1, . . . , N we have In addition, for all v ∈ V 0 , we have Hence, u is the unique solution of Problem 9.7.
We now consider Problem 9.4 in the related special case Recall that W is defined in (9.11) and V b,y = {v ∈ V | T i (y)v−b i = 0, i = 1, . . . , N }.
While existence of solutions is considered in Proposition 9.4, we now derive an equivalent characterization in terms of suitable basis functions. For all y ∈ M, we introduce y-dependent counterparts φ i (y) ∈ V of φ i by replacing T i by T i (y) in (9.16). The associated Gramian matrix is denoted by A(y). We also define M = {y ∈ M | T i (y), i = 1, . . . , N, are linearly independent}. as obtained from Proposition 9.7, we exploit (9.20), the definition of A(y) and U (y), to obtain J(u(y)) = 1 2 U (y) A(y)U (y) − 1 2 a(φ 0 , φ 0 ) Now the following characterization of solutions to Problem 9.8 is an immediate consequence of Lemma 9.3 and Proposition 9.7. Proposition 9.8. Assume that Problem 9.8 is non-degenerate in the sense that all solutions (u, x) satisfy x ∈ M . Then (u, x) ∈ V × M solves Problem 9.4, if and only if x ∈ M is a minimizer of ξ + V on M with ξ given in (9.19) and

9.4.
Regularity of Green's functions. We provide the following regularity result for Green's functions of eighth order problems on unbounded domains as exploited in Section 6.5.