On differentiability of the membrane-mediated mechanical interaction energy of discrete–continuum membrane–particle models
We consider a discrete–continuum model of a biomembrane with embedded particles. While the membrane is represented by a continuous surface, embedded particles are described by rigid discrete objects which are free to move and rotate in lateral direction. For the membrane we consider a linearized Canham–Helfrich energy functional and height and slope boundary conditions imposed on the particle boundaries resulting in a coupled minimization problem for the membrane shape and particle positions.
When considering the energetically optimal membrane shape for each particle position we obtain a reduced energy functional that models the implicitly given interaction potential for the membrane-mediated mechanical particle–particle interactions. We show that this interaction potential is differentiable with respect to the particle positions and orientations. Furthermore we derive a fully practical representation of the derivative only in terms of well defined derivatives of the membrane. This opens the door for the application of minimization algorithms for the computation of minimizers of the coupled system and for further investigation of the interaction potential of membrane-mediated mechanical particle–particle interaction.
The results are illustrated with numerical examples comparing the explicit derivative formula with difference quotient approximations. We furthermore demonstrate the application of the derived formula to implement a gradient flow for the approximation of optimal particle configurations.
This research has been partially supported by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 ”Scaling Cascades in Complex Systems”, Project AP01 ”Particles in lipid bilayers”
Membrane proteins are crucial for various processes that involve the shaping of biological membranes. In some of these cases the proteins act as inclusions in the membrane and induce local deformations in the vicinity of the interface between the membrane and the proteins . As the membrane itself consists of a lipid bilayer, which in lateral direction can be seen as a fluid , particles are able to move easily within the membrane. Since the local protein-induced membrane deformation implicitly defines a membrane mediated mechanical protein–protein interaction proteins can tend to assemble in energetically preferable patterns .
Various research has been done in order to gain more insight into the shaping of membranes as well as membrane-mediated particle-particle interactions. A common description of biological membranes is based on a representation as continuous surface that minimizes an elastic energy [3, 13]. It was later shown that there are long-range interactions between particles that are predominantly membrane-mediated . Since then further work has been done to investigate particle interactions within elasticity models. Typically a flatness assumption on the membrane is made and particles are modeled as circular disks or points and their coupling to the membrane is prescribed by radially symmetric boundary conditions for which the interaction energy can either be computed analytically or approximately by asymptotic expansion [17, 21, 6, 23, 8]. However, it turns out that the shape of particles has a significant impact on their interaction . More recent work is also interested in numerical computations with pattern formation of many non-circular particles [12, 15], and attention was also given to situations where the flatness assumption is no longer fulfilled [19, 20]. Also more elaborate models for proteins in continuum elastic models have recently been considered in .
For efficient computations with moving particles it therefore is desirable to quantify the forces exerted on the particles by the membrane in a framework that is as widely applicable as possible. General results in this direction have been obtained based on arguments from differential geometry . The methods derived therein give insight into the qualitative behavior of particle interactions, but—to the best of the authors’ knowledge—they have not yet been made fully available for numerical computations.
In this paper we consider a discrete–continuum model where the membrane is modeled as a continuous graph minimizing a linearized Canham–Helfrich bending energy and where an arbitrary amount of particles are embedded into the membrane. These particles are modeled as discrete entities which are coupled to the membrane through certain boundary conditions. As particles are free to move in the membrane, those boundary conditions depend on each particle’s position. Consequently, the overall system’s energy given fixed boundary conditions and an optimal membrane shape can be written as a function of the particle positions, which we call the interaction energy.
In this setting we propose a method to prove differentiability of the interaction energy for arbitrary shapes and boundary conditions. Furthermore, we derive an expression for the derivative that can be evaluated numerically within a finite element scheme and where the evaluation error is bounded in terms of the discretization error of the finite element approximation. The proof is based on an application of the implicit function theorem and ideas from shape calculus [14, 4]. As such, the method is rather general and hence it naturally extends to a wider class of models that for example use nonlinear elastic energies or certain other membrane–particle couplings.
In the following we give an outline of this paper. In \crefsec:model we introduce the Canham–Helfrich energy in Monge-gauge as a model for the membrane and parametric boundary conditions for the coupling of the particles. \crefsec:interaction is then concerned with further mathematical notation that we use in order to define the interaction energy. There we also reformulate the parametric boundary conditions as linear constraints by using trace operators and appropriate projection operators. Afterwards, in \crefsec:differentiation, we prove differentiability of the interaction energy and derive a numerically feasible expression for the gradient. Finally, \crefsec:examples shows some example computations that illustrate that the derived formula can indeed be applied in a numerical scheme.
2. Membrane and particle model
To model the membrane itself we choose the well-established linearized Canham–Helfrich model in Monge-gauge. Given a 2-dimensional reference domain and a function , the membrane shape is described by the graph of . The bending energy of this membrane is approximated by
where and denote the bending rigidity and the surface tension, respectively. It is noted that this is already a linearized formulation of the bending energy that makes the assumption that the membrane is “rather flat” with respect to the reference domain .
In absence of particles that interact with the membrane, this model determines the stationary shape of the membrane solely by minimizing this energy. In the following we explain how the embedded particles are coupled to the membrane, before we state the model problem that is central to this paper.
For simplicity we first consider a single transmembrane particle that interacts with the membrane. Such a particle is not merely connected to the membrane but rather is included in it, comparable to a wedge in the membrane. We assume that the particle is rigid, which means that it stays constant in shape, and we approximate it by a rigid 3-dimensional shape .
Suppose furthermore that this particle has a hydrophobic belt, i. e. a region to which the membrane connects preferentially, and that the particle’s belt is approximated by a curve . We assume that is a simple closed curve that can be parameterized over the 2-dimensional Euclidean plane. This means that there exists a simple closed curve and a continuous function such that . This gives rise to the boundary condition which models that the membrane is connected to the particle at the interface . It is common to impose the additional constraint that the membrane attaches to the interface with a fixed slope. This is modeled via a function describing the slope and via the boundary condition Here is an oriented unit normal on and denotes the normal derivative of on .
Those constraints do not yet account for the fact that the particle is in principle free to move in space. To this end we parameterize the current position of the particle using translations along the -axes and rotations around the -axes. More precisely, let be a reference state of the particle that is centered in the origin, and let be the -rotation matrix around the -axis. Then we define the parameterized particle as
The curves and that are associated to could be parameterized like this, too. However, we decide to only approximate this parameterization instead. We make the simplifying assumption that the reference set is oriented in such a way that the belt is flat, by which we mean that is small. And, in the spirit of the Monge-gauge linearization, we also assume that and are small so that we can apply a small angle approximation in the rotation matrices. We define
Under these geometric assumptions and upon linearization it is justified to approximate the actual reparameterized curve by
The zeroth and first order boundary condition then essentially stem from a linear perturbation of the reference conditions. We use the shortcut and define the helper
Using these functions we get
Together with the boundary conditions then read
Since those boundary conditions are linear in , , and it is convenient to include the variability of these variables into the boundary conditions, which then become parametric boundary conditions:
Therefore, , , and are no longer relevant for describing the particle’s position and are now rather implicit to the boundary conditions. A particle in the model is then solely determined by its reference curve , its reference boundary conditions , and its position in the Euclidean plane.
In the case where multiple particles are present we state the constraints analogously by imposing the above constraints for each particle separately.
3. Interaction energy
Before we can formulate the final model problem, we need to introduce some notation. We also augment the parametric boundary conditions by zero boundary conditions on the outer boundary and we are going to reformulate these conditions as a simple linear constraint.
We consider particles with reference curves , height profiles and slopes . Given a particle configuration we define the curves
and . Their union is where we use and for the sake of a consistent notation. For we denote the set enclosed by by and the union of these is denoted by . We define the p-dependent reference domain as .
Based on this we define the interior of the set of feasible particle configurations as
and consequently we have .
Now suppose . Then we define the trace operators
Here is the unit outer normal on with respect to the domain . We define the joint trace operator by . We also define , and .
In the next step we prove a useful reformulation of the parametric boundary conditions as linear constraints. To this end, define for
and for let
We furthermore define for .
Lemma 3.1 (parametric boundary conditions as linear constraint).
Let . Then is well-defined and a projection operator such that for all
holds if and only if there exists a such that for all and
For equation (3.2) is equivalent to .
The statement for is trivial, hence let . First we show that is indeed well-defined. To that end note that the adjoint is given by
It is readily seen that the scalar product is positive definite on the space that is spanned by the . Moreover, the are linearly independent on and thus is symmetric positive definite and in particular invertible. This makes well-defined.
We define the shorthands and and note that and . Now, is obviously linear and
Therefore is a linear projection operator.
Let and . Suppose that (3.2) holds. If we define , then we get
This implies for almost-every with , , and
where we used the relation . Therefore (3.3) holds with and .
For notational convenience we define . The set of feasible membranes given the particle configuration p is defined by
We use the shorthand to define the interaction energy
where we use the convention .
We note that the interaction energy is well-defined for smooth enough. This statement can be proven analogously to [7, Theorem 1].
Altogether, our model problem reads
4. Differentiation of the reduced interaction energy
In this section we investigate the differentiability of on . First we derive a technical result which shows that the admissible membrane sets are isomorphic and which allows us to pose our problem locally over a fixed reference domain . Afterwards we apply the implicit function theorem to prove differentiability of the reduced interaction potential and use matrix calculus to derive an explicit and numerically feasible expression for the first order derivatives.
4.1. Trace-preserving diffeomorphisms between the reference domains
In this part we construct a local diffeomorphism between the domains that preserves the boundary conditions. The construction is based on ordinary differential equations (ODEs) and in particularly requires the following result from ODE theory.
Let be an open connected set, and let be Lipschitz-continuous. For and let be the unique solution of the ordinary differential equation
Then the map defined by fulfills and is an -diffeomorphism onto its image for all .
For all with as the unique solution of
holds . Also, for all with as the unique solution of
The global existence and uniqueness of is a consequence of the Lipschitz-continuity of and the well-known Picard–Lindelöf theorem. In particular, is well-defined.
The smoothness of and the characterization of its derivatives is a consequence of [11, Theorem 3.1, Theorem 4.1]
Concerning the inverse of , let be the unique solution of
Let . From
and the uniqueness of we infer , and in particular also . Therefore exists and is given by . Again, the smoothness of implies -smoothness of , and consequently is an -diffeomorphism. ∎
In the following we restrict ourselves to a special class of vector fields that is described in the result below. We show afterwards that the diffeomorphisms induced by such vector fields have a certain trace preserving property that again can be used to construct an isomorphism between the admissible membrane sets.
Let and . Then there exists an open neighborhood of and a Lipschitz-continuous map such that for all , , and holds
This is a consequence of the Whitney extension theorem, \crefthm:whitney. It uses the fact that the are pairwise disjoint and that the right hand sides in (4.3) smoothly extend to . ∎
Let be as in \crefthm:existenceCV for and let as in \crefthm:propertiesODE be induced by . Then , , and for all holds
From follows immediately that , i. e. .
We now prove (4.4). First we show that preserves the boundaries. From the properties (4.3) it follows that for and the value of can be computed explicitly. Recalling the definition of in (2.2) and the convention , this value is given by
Now, let . From (4.5) we infer for that
and by (4.9)
Altogether this proves equation (4.4).
In order to show the equality we define the set
Now let and assume that . By continuity of this would imply that there exists a such that and therefore, by definition of , there would exist a such that . As of this would be a contradiction to the uniqueness of . ∎
Let be as in \crefthm:XExistence. Then for all the map
is well-defined and an isomorphism.
From \crefthm:XExistence we know for every that the restriction is a -diffeomorphism onto . Because is compact we can assume without loss of generality that is uniformly bounded, i. e. that there exists a such that for all and holds . Elsewise we may replace by an appropriate sub-neighborhood. Hence, [1, Theorem 3.35] is applicable and the map
is well-defined and an isomorphism, and in particular also the restriction is well-defined and an isomorphism onto its image.
It remains to show that . Suppose and . Because of the trace preserving property (4.4) and by definition of and it follows that and , and so . ∎
In this part we use the maps from \crefthm:XExistence and from \crefthm:PhiExistence to transform the domain of definition for the functions from to . Afterwards we apply the implicit function theorem to derive a differentiability result.
For and the transformed energy is defined as
and the transformed reduced interaction energy is
We write the affine linear subspace as where is a linear subspace and is a function such that .
For notational convenience we define to be the differential with respect to the spatial coordinates.
Let and define
Equation (4.11) is a direct application of \crefthm:transformationDerivatives applied to .
Furthermore, for all we have
As of we know that and are both times continuously differentiable. Because the are diffeomorphisms with and because is continuous we have . Because is compact, we can assume without loss of generality that there exists a such that , else we replace by a suitable sub-neighborhood. Consequently, the integrand of is times differentiable with respect to q. Moreover, the integrand is even smooth with respect to and hence application of the dominated convergence theorem yields . ∎
There exists a neighborhood of such that and for all and multi-indices with holds
In particular, if and then