Local sensitivity analysis of the “membrane shape equation” derived from the Helfrich energy

Abstract

The Helfrich energy is commonly used to model the elastic bending energy of lipid bilayers in membrane mechanics. The governing differential equations for certain geometric characteristics of the shape of the membrane can be obtained by applying variational methods (minimization principles) to the Helfrich energy functional and are well studied in the axisymmetric framework. However, the Helfrich energy functional and the resulting differential equations involve a number of parameters, and there is little explanation of the choice of parameters in the literature, particularly with respect to the choice of the “spontaneous curvature” term that appears in the functional. In this paper, we present a careful analytical and numerical study of certain aspects of parametric sensitivity of Helfrich’s model. Using simulations of specific model systems, we demonstrate the application of our scheme to the formation of spherical buds and pearled shapes in membrane vesicles.

Keywords

Helfrich model membrane curvature spontaneous curvature parametric sensitivity analysis

1. Introduction

The elastic behavior of lipid bilayers has been studied using mechanical models for 50 years or more [1, 2]. A pseudoelastic strain energy functional was first described by Canham [1] as a way to explain the biconcave shape of red blood cells and subsequently by Helfrich [2] to describe the mechanical behavior of lipid bilayers in various situations. This energy functional has now become the most accepted model for describing the mechanical properties of the cell membrane. This model has been used to study the shapes associated with whole cells, particularly that of the red blood cell [1, 3, 4]. Subsequently, the Helfrich energy has been used to study the different shapes associated with vesicles as a function of pressure, volume, and membrane composition [5].

While the original model proposed by Helfrich considered the lipid bilayer as a thin shell with negligible thickness, over the years there have been many mathematical developments to represent the different physical properties of biological membranes. Some of these include the area difference model [6] and the mattress model [7]. Among these, the spontaneous curvature model has been one of the most popularly used models to represent the asymmetry between the two leaflets of the bilayer (Figure 1(a)) [8]. The idea of spontaneous curvature to represent the asymmetry between the two leaflets of the lipid bilayer was first introduced by Helfrich in his seminal 1973 paper [2]. Subsequently, this term has been used to capture compositional asymmetry, protein-induced spontaneous curvature, coat proteins etc. (see [8] for a more detailed discussion).

For many problems of biophysical interest, it is necessary to consider the heterogeneity in composition across cell membranes [9]. In such cases, the spontaneous curvature function is often modeled as a spatially varying quantity rather than as a constant to represent different membrane domains without introducing a discontinuity for computational purposes [10, 11]. Alongside theoretical modeling efforts, there have been significant advances in computational methods for solving the partial differential equations resulting from the minimization procedure of the Helfrich energy [12, 13]. A vast majority of the simulations are implemented under the assumption of axisymmetry; this assumption enables us to transform the partial differential equations describing the shape of the membrane into a system of ordinary differential equations (ODEs), which are then equipped with appropriate boundary conditions to be solved [14] (Figure 1(b)). However, a major challenge associated with such simulations remains to be the number of free parameters associated with the spontaneous curvature function.

We and others have found that the specific choice of the spontaneous curvature function and the resulting parameters play an important role in determining the shape of the membrane [10, 15]. For certain parameters, such as membrane bending modulus, there exist sufficient experimental measurements to establish a range of physically relevant values [16]. However, for the spontaneous curvature function, such measurements are limited or do not exist in forms that are always amenable to modeling. As a result, throughout the literature, spontaneous curvature has been represented by different functions, using a wide range of parameter values. Thus, it seems that there is a need for a better understanding of the sensitivity of the spontaneous curvature model with respect to the various parameters involved.

To address these issues and gain some insight into the role of parameters, as a first step, in this work we study the local sensitivity of solutions of certain equations associated with the Helfrich energy model with spontaneous curvature. We note that while parametric sensitivity analysis methods are well documented for initial value problems (IVPs) [17], in this work the problems that we will consider are boundary value problems (BVPs). We will show that parametric sensitivity analysis can provide valuable insight into the impact of different parameters on energy minimization and a tool for better understanding critical phenomena associated with the shape of the membrane.

In what follows, we present a summary of the Helfrich model with spontaneous curvature and axisymmetric parametrization in Sections 2 and 3, development of the parametric sensitivity analysis method in Section 4, present some numerical results in Section 5, and end with our interpretations and conclusions in Section 6.

2. Overview of the Helfrich energy with spontaneous curvature

The Helfrich energy serves as the constitutive equation for the lipid bilayer. We use a modified version of the Helfrich energy that includes spatially varying spontaneous curvature C as opposed to a constant uniform value, as in [10, 18, 19]:

w = κ {[H - C]}^{2} + κ_{G} K,

(1)

where w is the energy per unit area, $κ$ is the bending modulus, H is the mean curvature, $κ_{G}$ is the Gaussian modulus, and K is the Gaussian curvature. This form of the energy density accommodates the local heterogeneity in the spontaneous curvature C. Note that w differs from the standard Helfrich energy [2] by a factor of 2, which is accounted for by setting the value of $κ$ to be twice that of the standard bending modulus typically encountered in the literature.

2.1. Equations of motion

We refer the interested reader to [18] for a detailed derivation of the governing equations. We make the following simplifying assumptions for simplicity in our model. We assume that the bending modulus and Gaussian modulus are uniform, the pressure difference across the membrane is zero, and there are no externally applied forces. Furthermore, we assume that the membrane is areally incompressible and introduce a Lagrange multiplier $λ$ to impose this constraint. It can be shown that among all areally incompressible surfaces, the one that minimizes the Helfrich energy functional has mean and Gaussian curvatures H and K satisfying

κ Δ [(H - C)] + 2 κ (H - C) (2 H^{2} - K) - 2 κ H {(H - C)}^{2} = 2 λ H

(2)

and

\frac{\partial λ}{\partial θ^{α}} = 2 κ (H - C) \frac{\partial C}{\partial θ^{α}},

(3)

where $θ^{α}$ ( $α = 1, 2$ ) denotes the surface coordinates and $Δ$ is the Laplace–Beltrami operator. Equation (2) gives the normal force balance on the membrane. Equation (3) represents the tangential force balance along the membrane. For a heterogeneous membrane (nonuniform or coordinate-dependent $C (θ^{α})$ ), this equation also represents the variation of $λ$ along the surface [18]. Detailed discussions of interpretation of $λ$ are given in [18, 20, 21].

2.2. Choice of spontaneous curvature function

In the existing literature, for cases where the situation under consideration can be classified as an axisymmetric problem, the spontaneous curvature function, $C (θ^{α})$ , is often represented using a hyperbolic tangent function to capture differences in the “coated” and “uncoated” regions. In this work, we use two forms of this function. Each form involves certain free parameters and we perform our parametric sensitivity analysis to understand the extent to which these parameters can affect the elastic bending energy and the shape of the membrane. These functions are given as

Type I : C (u) = 0.5 C_{0} [1 - \tanh [ξ (u - u_{0})]],

(4)

Type II : C (u) = - 0.5 C_{0} \frac{u - u_{0}}{u_{0}} [1 - \tanh [ξ (u - u_{0})]] .

(5)

In these formulas, we used the generic variable u instead of $θ^{α}$ . This variable can be viewed as one of the coordinates used in the parametrization of the surface of revolution representing the membrane, as will be described in the next section.

3. Axisymmetric parametrization

Under the assumption of axisymmetry, which in part implies that the surface representing the membrane is a surface of revolution obtained by rotating a regular curve about an axis, the governing equations shown in equations (2) and (3) can be recast as a system of ODEs. When equipped with appropriate boundary conditions, reflecting the geometric and physical constraints of the problem, the solution of the resulting boundary value problem will determine the shape of the membrane. In what follows, we first describe this viewpoint in more detail, and then we summarize the various ways in which our ODE system can be represented, depending on the choice of arc-length or area parametrization.

Under the assumption of axisymmetry, without loss of generality, we may assume that the surface of the membrane is generated by rotating a curve in the right half of the rz-plane, about the z-axis. We let $p : [0, S] \to R^{2}$ , defined by $s \mapsto (r (s), z (s))$ , denote the arc-length parametrization of the generating curve (or the profile curve). Then

r (s, φ) = (r (s) \cos φ, r (s) \sin φ, z (s))

(6)

is a parametrization of the membrane surface (here $0 \leq s \leq S$ and $0 \leq φ \leq 2 π$ ). Since $[r' (s {)]}^{2} + [z' (s {)]}^{2} = 1$ , there exists a function $Ψ (s)$ such that

r' (s) = \cos Ψ (s), z' (s) = \sin Ψ (s),

(7)

where $Ψ$ can be viewed as the angle made by the curve with the horizontal.

Using well-known formulas for the Gaussian and mean curvatures of a surface of revolution (see, for example, [22]), it is easy to show that if we orient the surface with the unit normal vector that points in the direction of

- \frac{\partial r}{\partial s} \times \frac{\partial r}{\partial φ},

then

r Ψ' = 2 rH - \sin Ψ,

(8)

K = H^{2} - {(H - \frac{\sin Ψ}{r})}^{2} .

(9)

To be able to write the final equations as first-order equations, we introduce the auxiliary variable L as

L = r [H - C]' .

(10)

Now, using these equalities and equation (2), it is easy to see that

H' (s) = r (s)^{- 1} L (s) + C' (s),

(11)

r (s)^{- 1} L' (s) = 2 H (s) [{(H (s) - C (s))}^{2} + \frac{λ (s)}{κ}] - 2 (H (s) - C (s)) [H {(s)}^{2} + {(H (s) - \frac{\sin Ψ (s)}{r (s)})}^{2}],

(12)

where we have used the fact that if f is any function defined on our surface of revolution whose values depend only on s, then

Δ f = \frac{[r (s) f' (s)]'}{r (s)} .

Finally, equation (3) gives

λ' (s) = 2 κ (H (s) - C (s)) C' (s) .

(13)

Remark 1. In what has been discussed so far, we have referred to the assumption of axisymmetry a number of times without clearly explaining what this assumption entails. Now we are at a position to give a careful description of this key assumption. In the present work, we say that our problem falls into the axisymmetric category if and only if:

The surface of the membrane is a surface of revolution parametrized by $r (s, φ)$ , as described already.

The values of the spontaneous curvature function C depend only on s.

3.1. Arc-length formulation

In this formulation, as described previously, s is the arc-length along the membrane and the unknown functions are

r (s), z (s), Ψ (s), H (s), L (s), λ (s) .

(14)

These functions must satisfy the following system of ODEs on the interval $[0, S]$ :

r' = \cos Ψ, z' = \sin Ψ, r Ψ' = 2 rH - \sin Ψ, H' = r^{- 1} L + C',

(15a)

r^{- 1} L' = 2 H [{(H - C)}^{2} + \frac{λ}{κ}] - 2 (H - C) [H^{2} + {(H - r^{- 1} \sin Ψ)}^{2}],

(15b)

λ' = 2 κ (H - C) C' .

(15c)

Before we attempt to solve this system of equations numerically, we need to provide the system with appropriate boundary conditions. Certain suitable classes of boundary conditions for this system will be introduced in Section 5.

3.2. Dimensionless arc-length formulation

In this formulation, we fix two positive constants $R_{0}$ and $κ_{0}$ that will be used to nondimensionalize the variables in the arc-length formulation by setting $t = s / R_{0}$ and $\tilde{κ} = κ / κ_{0}$ , and defining

x (t) = \frac{1}{R_{0}} r (R_{0} t), y (t) = \frac{1}{R_{0}} z (R_{0} t), ψ (t) = Ψ (R_{0} t), h (t) = R_{0} H (R_{0} t),

(16a)

c (t) = R_{0} C (R_{0} t), l (t) = R_{0} L (R_{0} t), \tilde{λ} (t) = \frac{R_{0}^{2}}{κ_{0}} λ (R_{0} t) .

(16b)

Note that $s \in [0, S]$ and $t \in [0, T]$ , where $S = R_{0} T$ .

A simple application of chain rule shows that the six dimensionless unknown functions $x (t)$ , $y (t)$ , $ψ (t)$ , $h (t)$ , $l (t)$ , and $\tilde{λ} (t)$ must satisfy the following system of ODEs on the interval $[0, T]$ :

\overset{\cdot}{x} = \cos ψ, \overset{\cdot}{y} = \sin ψ, x \overset{\cdot}{ψ} = 2 xh - \sin ψ, \overset{\cdot}{h} = x^{- 1} l + \overset{\cdot}{c},

(17a)

x^{- 1} \overset{\cdot}{l} = 2 h [{(h - c)}^{2} + \frac{\tilde{λ}}{\tilde{κ}}] - 2 (h - c) [h^{2} + {(h - x^{- 1} \sin ψ)}^{2}],

(17b)

\overset{\cdot}{\tilde{λ}} = 2 \tilde{κ} (h - c) \overset{\cdot}{c} .

(17c)

(Here, the dot denotes the derivative with respect to t.)

Remark 2. It is worth mentioning that if we assume that $C (s)$ has the form

C (s) = 0.5 C_{0} [1 - \tanh [ξ (s - s_{0})]],

(18)

where $C_{0}$ , $ξ$ , and $s_{0}$ are certain constants, then

c (t) = R_{0} C (R_{0} t) = 0.5 C_{0} R_{0} [1 - \tanh [ξ (R_{0} t - s_{0})]] = 0.5 C_{0} R_{0} [1 - \tanh [ξ R_{0} (t - \frac{s_{0}}{R_{0}})]] .

(19)

So if we let $γ = R_{0} ξ$ and $t_{0} = s_{0} / R_{0}$ , then

c (t) = 0.5 C_{0} R_{0} [1 - \tanh [γ (t - t_{0})]] .

(20)

Similar calculations show that if we assume $C (s)$ has the form

C (s) = (- 0.5 C_{0}) (\frac{s - s_{0}}{s_{0}}) [1 - \tanh [ξ (s - s_{0})]],

(21)

then

c (t) = (- 0.5 C_{0}) R_{0} (\frac{t - t_{0}}{t_{0}}) [1 - \tanh [γ (t - t_{0})]] .

(22)

where $γ = R_{0} ξ$ and $t_{0} = s_{0} / R_{0}$ .

Remark 3. In the following sections of this paper, we will be interested in the derivatives of the solution functions with respect to such input parameters as $C_{0}$ . These derivatives will be referred to as sensitivities. An application of the chain rule shows that sensitivities calculated using dimensionless variables are constant multiples of sensitivities computed using the original variables. For example,

\frac{\partial h}{\partial C_{0}} |_{t = \hat{t}} = R_{0} \frac{\partial H}{\partial C_{0}} |_{s = R_{0} \hat{t}},

(23)

\frac{\partial h}{\partial γ} |_{t = \hat{t}} = \frac{1}{R_{0}} \frac{\partial h}{\partial ξ} |_{t = \hat{t}} = \frac{1}{R_{0}} R_{0} \frac{\partial H}{\partial ξ} |_{s = R_{0} \hat{t}} = \frac{\partial H}{\partial ξ} |_{s = R_{0} \hat{t}},

(24)

\frac{\partial h}{\partial t_{0}} |_{t = \hat{t}} = \frac{1}{1 / R_{0}} \frac{\partial h}{\partial s_{0}} |_{t = \hat{t}} = \frac{1}{1 / R_{0}} R_{0} \frac{\partial H}{\partial s_{0}} |_{s = R_{0} \hat{t}} = R_{0}^{2} \frac{\partial H}{\partial s_{0}} |_{s = R_{0} \hat{t}} .

(25)

3.3. Area formulation

Here we introduce the new variable a as the area of the surface of revolution produced by rotating the segment of the curve spanned as the arc-length varies from 0 to s, that is,

a (s) = \int_{0}^{s} 2 π r (u) d u .

(26)

Since there is a one-to-one relationship between a and s, we can view the six unknown functions as functions of a rather than s. As pointed out in [10], this formulation has the advantage of prescribing the total area of the membrane as the domain size (rather than the corresponding arc-length), which is more physical and amenable to laboratory measurements.

A simple application of the chain rule shows that, if we denote the total area (of the membrane) by A, then the unknown functions

r (a), z (a), Ψ (a), H (a), L (a), λ (a)

(27)

must satisfy the following system of ODEs on the interval $[0, A]$ :

2 π rr' = \cos Ψ, 2 π rz' = \sin Ψ, 2 π r^{2} Ψ' = 2 rH - \sin Ψ, 2 π r^{2} H' = L + 2 π r^{2} C',

(28a)

2 π L' = 2 H [{(H - C)}^{2} + \frac{λ}{κ}] - 2 (H - C) [H^{2} + (H - r^{- 1} \sin Ψ)^{2}],

(28b)

λ' = 2 κ (H - C) C' .

(28c)

As noted, we will discuss suitable boundary conditions in Section 5.

Remark 4. We emphasize that, in these formulas, the functions r, z, $Ψ$ , H, L, $λ$ , and C, which were originally introduced as functions of s, are viewed as functions of the variable a; indeed, if we let $g (s) = \int_{0}^{s} 2 π r (u) d u$ , any function of s can be written as a function of a by replacing each instance of s in the expression of the function by $g^{- 1} (a)$ ; so if we wanted to be completely rigorous, in equations (27) and (28) instead of r, z, $Ψ$ , H, L, $λ$ , and C, we should have written $\hat{r}$ , $\hat{z}$ , $\hat{Ψ}$ , $\hat{H}$ , $\hat{L}$ , $\hat{λ}$ , and $\hat{C}$ , where $\hat{r} = r (g^{- 1} (a))$ , $\hat{z} = z (g^{- 1} (a))$ , etc. It is customary to abuse notation and denote these new variables by the same symbols as the original variables. In Section 6.1, we will make use of this very elementary observation to better understand the relationship between the simulation results obtained using the arc-length formulation and those obtained using the area formulation.

3.4. Dimensionless area formulation

In this formulation, we fix two positive constants, $R_{0}$ and $κ_{0}$ , which will be used to nondimensionalize the variables in the previous formulation by setting $α = a / 2 π R_{0}^{2}$ and $\tilde{κ} = κ / κ_{0}$ , and defining:

x (α) = \frac{1}{R_{0}} r (2 π R_{0}^{2} α), y (α) = \frac{1}{R_{0}} z (2 π R_{0}^{2} α), ψ (α) = Ψ (2 π R_{0}^{2} α), h (α) = R_{0} H (2 π R_{0}^{2} α),

(29a)

c (α) = R_{0} C (2 π R_{0}^{2} α), l (α) = R_{0} L (2 π R_{0}^{2} α), \tilde{λ} (α) = \frac{R_{0}^{2}}{κ_{0}} λ (2 π R_{0}^{2} α) .

(29b)

Note that $a \in [0, A]$ and $α \in [0, α_{max}]$ , where $A = 2 π R_{0}^{2} α_{max}$ .

A simple application of the chain rule shows that the six dimensionless unknown functions, $x (α)$ , $y (α)$ , $ψ (α)$ , $h (α)$ , $l (α)$ , and $\tilde{λ} (α)$ , must satisfy the following system of ODEs on the interval $[0, α_{max}]$ :

x \overset{\cdot}{x} = \cos ψ, x \overset{\cdot}{y} = \sin ψ, x^{2} \overset{\cdot}{ψ} = 2 xh - \sin ψ, x^{2} \overset{\cdot}{h} = l + x^{2} \overset{\cdot}{c},

(30a)

\overset{\cdot}{l} = 2 h [{(h - c)}^{2} + \frac{\tilde{λ}}{\tilde{κ}}] - 2 (h - c) [h^{2} + {(h - x^{- 1} \sin ψ)}^{2}],

(30b)

\overset{\cdot}{\tilde{λ}} = 2 \tilde{κ} (h - c) \overset{\cdot}{c} .

(30c)

Remark 5. It is worth mentioning that if we assume that $C (a)$ has the form

C (a) = (- 0.5 C_{0}) (\frac{a - a_{0}}{a_{0}}) [1 - \tanh [ξ (a - a_{0})]],

(31)

where $C_{0}$ , $ξ$ , and $a_{0}$ are certain constants, then

\begin{matrix} c (α) = R_{0} C (2 π R_{0}^{2} α) = - 0.5 C_{0} R_{0} (\frac{2 π R_{0}^{2} α - a_{0}}{a_{0}}) [1 - \tanh [ξ (2 π R_{0}^{2} α - a_{0})]] \\ = - 0.5 C_{0} R_{0} \frac{α - a_{0} / (2 π R_{0}^{2} α)}{a_{0} / (2 π R_{0}^{2} α)} [1 - \tanh [ξ 2 π R_{0}^{2} (α - \frac{a_{0}}{2 π R_{0}^{2}})]] . \end{matrix}

(32)

So if we let $γ = 2 π R_{0}^{2} ξ$ and $α_{0} = a_{0} / 2 π R_{0}^{2}$ , then

c (α) = - 0.5 C_{0} R_{0} (\frac{α - α_{0}}{α_{0}}) [1 - \tanh [γ (α - α_{0})]] .

(33)

Similarly, one can show that if $C (a)$ has the form

C (a) = (0.5 C_{0}) [1 - \tanh [ξ (a - a_{0})]],

(34)

then

c (α) = 0.5 C_{0} R_{0} [1 - \tanh [γ (α - α_{0})]],

(35)

where $γ = 2 π R_{0}^{2} ξ$ and $α_{0} = a_{0} / 2 π R_{0}^{2}$ .

Remark 6. The same argument as that discussed in Remark 3 shows that sensitivities calculated using dimensionless variables are constant multiples of sensitivities computed using the original variables.

We conclude this section by making a few comments about the choice of the positive constants $R_{0}$ and $κ_{0}$ that are used to nondimensionalize the variables. The following observations indeed show the existence of certain symmetries in the system of ODEs under consideration. Generally, it is true that symmetries of differential equations can be used to gain valuable information about the equations under consideration. In Section 6.2, we will use these observations to examine a conjecture about the solutions of our system of equations.

Let $η$ be any positive number. It is easy to see that $(x, y, ψ, h, l, \tilde{λ})$ is a solution of the ODE system using $(R_{0}, κ_{0})$ if and only if $(x, y, ψ, h, l, η \tilde{λ})$ is a solution of the system using $(R_{0}, η κ_{0})$ . That is, changing $κ_{0}$ merely amounts to scaling $\tilde{λ}$ . In our numerical experiments we choose $κ_{0} = κ$ so that $\tilde{κ} = 1$ .

For the dimensionless arc-length formulation, $(x (t), y (t), ψ (t), h (t), l (t), \tilde{λ} (t))$ is a solution of the ODE system with $c (t)$ defined using parameters $C_{0}, R_{0}, γ, t_{0}$ on the interval $[0, T]$ using $(R_{0}, κ_{0})$ if and only if

(\frac{1}{η} x (η t), \frac{1}{η} y (η t), ψ (η t), η h (η t), η l (η t), η^{2} \tilde{λ} (η t))

is a solution of the ODE system with $c (t)$ defined using parameters $C_{0}, η R_{0}, η γ, t_{0} / η$ on the interval $[0, T / η]$ using $(η R_{0}, κ_{0})$ .

For the dimensionless area formulation, $(x (α), y (α), ψ (α), h (α), l (α), \tilde{λ} (α))$ is a solution of the ODE system with $c (α)$ defined using parameters $C_{0}, R_{0}, γ, α_{0}$ on the interval $[0, α_{max}]$ using $(R_{0}, κ_{0})$ if and only if

(\frac{1}{η} x (η^{2} α), \frac{1}{η} y (η^{2} α), ψ (η^{2} α), η h (η^{2} α), η l (η^{2} α), η^{2} \tilde{λ} (η^{2} α))

is a solution of the ODE system with $c (α)$ defined using parameters $C_{0}, η R_{0}, η^{2} γ, α_{0} / η^{2}$ on the interval $[0, α_{max} / η^{2}]$ using $(η R_{0}, κ_{0})$ .

Finally, we remark that when solving the boundary value problem using different values of $R_{0}$ (or $κ_{0}$ , or both), in order to be able to compare the results in the sense described here, we also need to change the boundary conditions accordingly.

4. Mathematical framework for sensitivity analysis

In this section, we give a brief overview of the theoretical framework of sensitivity analysis for a system of ODEs dependent on parameters. The subject is well studied in the context of initial value problems [17] and, as illustrated in this and in the following sections, the same key ideas can be employed to analyze the sensitivity of boundary value problems such as those that will be considered in this paper.

Consider the following system of ODEs for the unknown $\vec{x} (u) = [x_{1} (u), \dots, x_{k} (u {)]}^{T}$ involving parameters $\vec{p} = [p_{1}, \dots, p_{m}]^{T}$ :

F (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) = [\begin{matrix} f_{1} (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) \\ ⋮ \\ f_{k} (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) \end{matrix}] = \vec{0},

(36)

with appropriate initial or boundary conditions. Two key objectives of sensitivity analysis would be to provide answers to the following questions:

How can we compute the rate of change of solution with respect to each of the parameters? That is, we are interested in computing

\frac{\partial x_{i}}{\partial p_{j}} \forall 1 \leq i \leq k \forall 1 \leq j \leq m .

How can we compute the rate of change of some functional $W (\vec{x}, \vec{p}) = \int_{0}^{u_{max}} w (\vec{x}, \vec{p}) d u$ of the solution with respect to each parameter? That is, we are interested in computing

\frac{\partial W}{\partial p_{j}} \forall 1 \leq j \leq m .

Throughout this document, we will refer to $\partial x_{i} / \partial p_{j}$ and $\partial W / \partial p_{j}$ as sensitivities of solution components and sensitivities of the functional, respectively. Notice that under appropriate smoothness assumptions on the functions involved, it follows from the chain rule that

\frac{\partial W}{\partial p_{j}} = \int_{0}^{u_{max}} \frac{\partial w}{\partial p_{j}} + (D_{\vec{x}} w) \frac{\partial \vec{x}}{\partial p_{j}} d u,

(37)

and so the answer to the first question can be used to answer the second question. However, as we shall see, if the ultimate objective is just to compute the sensitivity of a functional, there might be more efficient tools available for the job.

We begin by describing a standard method to answer the first question. Define the sensitivity vectors ${\vec{s}}_{1}, \dots, {\vec{s}}_{m}$ as follows:

{\vec{s}}_{1} = \frac{\partial \vec{x}}{\partial p_{1}} = [\begin{matrix} \frac{\partial x_{1}}{\partial p_{1}} \\ ⋮ \\ \frac{\partial x_{k}}{\partial p_{1}} \end{matrix}], \dots, {\vec{s}}_{m} = \frac{\partial \vec{x}}{\partial p_{m}} = [\begin{matrix} \frac{\partial x_{1}}{\partial p_{m}} \\ ⋮ \\ \frac{\partial x_{k}}{\partial p_{m}} \end{matrix}] .

(38)

Now note that

F (\vec{x}, \overset{\cdot}{\vec{x}}, p_{1}, \dots, p_{m}) = \vec{0} \Rightarrow (D_{\vec{x}} F) \frac{\partial \vec{x}}{\partial p_{j}} + (D_{\overset{\cdot}{\vec{x}}} F) \frac{\partial \overset{\cdot}{\vec{x}}}{\partial p_{j}} + \frac{\partial F}{\partial p_{j}} = \vec{0} .

(39)

That is,

(D_{\overset{\cdot}{\vec{x}}} F) {\overset{\cdot}{\vec{s}}}_{j} + (D_{\vec{x}} F) {\vec{s}}_{j} + \frac{\partial F}{\partial p_{j}} = \vec{0} .

(40)

Thus, the sensitivities ${\vec{s}}_{j}$ themselves satisfy a system of ODEs (consisting of $k \times m$ equations). The initial or boundary conditions for ${\vec{s}}_{j}$ are obtained by taking the partial derivatives of the initial and boundary conditions satisfied by $\vec{x}$ with respect to $p_{j}$ .

By appending equation (40) consisting of $k \times m$ equations to the original system (equation 36), we can construct a system of ODEs with $k \times (m + 1)$ equations:

F (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) = \vec{0},

(41)

(D_{\overset{\cdot}{\vec{x}}} F) {\overset{\cdot}{\vec{s}}}_{j} + (D_{\vec{x}} F) {\vec{s}}_{j} + \frac{\partial F}{\partial p_{j}} = \vec{0} .

(42)

Solving this system with a fixed set of values ${\tilde{p}}_{1}, \dots, {\tilde{p}}_{m}$ as parameters and with the appropriate initial and boundary conditions (as described) gives the original unknowns together with the sensitivities ${\vec{s}}_{j} = \partial \vec{x} / \partial p_{j}$ evaluated at $({\tilde{p}}_{1}, \dots, {\tilde{p}}_{m})$ . If we are only interested in computing sensitivities with respect to one particular parameter $p_{j}$ , it is enough to append the k equations corresponding to $p_{j}$ to the original system (and hence solve a system with $2 k$ equations).

Now let us focus on the second question. As mentioned, we can answer the second question by computing each individual sensitivity using the method explained previously (see equation (37)). However, there is at least one more approach that can be used to directly compute the sensitivity of a functional W with respect to the parameters. In this work, we are interested in both sensitivities of solution components and sensitivities of certain functionals of solutions; hence, we will only employ the first approach. Nevertheless, for the sake of completeness, here we briefly describe this alternative approach. To explain this second approach, which is sometimes referred to as the “adjoint sensitivity analysis” [23], we need to make a simple observation which we state as a proposition.

Proposition 1. Let $F$ and w be as before. Suppose that $\vec{v} (u) = [v_{1} (u), \dots, v_{k} (u {)]}^{T}$ solves the following system of ODEs:

\frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F] - {\vec{v}}^{T} D_{\vec{x}} F = - D_{\vec{x}} w .

(43)

Then

(D_{\vec{x}} w) \frac{\partial \vec{x}}{\partial p_{j}} = - {\vec{v}}^{T} \frac{\partial F}{\partial p_{j}} - \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}] .

(44)

Proof. We have

\begin{matrix} (D_{\vec{x}} w) \frac{\partial \vec{x}}{\partial p_{j}} \overset{(43)}{=} {\vec{v}}^{T} D_{\vec{x}} F \frac{\partial \vec{x}}{\partial p_{j}} - \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F] \frac{\partial \vec{x}}{\partial p_{j}} \\ = {\vec{v}}^{T} D_{\vec{x}} F \frac{\partial \vec{x}}{\partial p_{j}} - \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}] + {\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \overset{\cdot}{\vec{x}}}{\partial p_{j}} \\ \overset{(39)}{=} {\vec{v}}^{T} [- \frac{\partial F}{\partial p_{j}}] - \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}] . \end{matrix}

It follows that if $\vec{v} (u)$ satisfies equation (43) then, for each $1 \leq j \leq m$ ,

\begin{matrix} \frac{\partial W}{\partial p_{j}} = \int_{0}^{u_{max}} \frac{\partial w}{\partial p_{j}} + (D_{\vec{x}} w) \frac{\partial \vec{x}}{\partial p_{j}} d u \\ = \int_{0}^{u_{max}} \frac{\partial w}{\partial p_{j}} - {\vec{v}}^{T} \frac{\partial F}{\partial p_{j}} - \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}] d u \\ = \int_{0}^{u_{max}} (\frac{\partial w}{\partial p_{j}} - {\vec{v}}^{T} \frac{\partial F}{\partial p_{j}}) d u - {[{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}]}_{u = 0}^{u = u_{max}} . \end{matrix}

Therefore, if we can choose side conditions for $\vec{v}$ such that the boundary term

{[{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}]}_{u = 0}^{u = u_{max}}

vanishes, then we can follow this two-step process to compute the sensitivities for the functional W:

Step 1. Solve the following problem (with appropriate side conditions as described previously) for the $2 k$ unknowns $\vec{x} (u) = [x_{1} (u), \dots, x_{k} (u {)]}^{T}$ and $\vec{v} (u) = [v_{1} (u), \dots, v_{k} (u {)]}^{T}$ :

(\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) = [\begin{matrix} f_{1} (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) \\ ⋮ \\ f_{k} (\vec{x}, \overset{\cdot}{\vec{x}}, \vec{p}) \end{matrix}] = \vec{0}, \\ \frac{d}{d u} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F] - {\vec{v}}^{T} D_{\vec{x}} F = - D_{\vec{x}} w . \end{matrix}

Step 2. For each $1 \leq j \leq m$ , compute

\frac{\partial W}{\partial p_{j}} = \int_{0}^{u_{max}} (\frac{\partial w}{\partial p_{j}} - {\vec{v}}^{T} \frac{\partial F}{\partial p_{j}}) d u .

5. Numerical sensitivity analysis of the Helfrich model

In this section, we illustrate the application of the theoretical considerations of sensitivity analysis by applying them to the ODEs of our problem of interest and by conducting numerical experiments for certain parameter choices. In Section 5.1, we set up a framework for the numerical sensitivity analysis of our boundary value problem. In Sections 5.2 to 5.5, we will discuss some of our numerical results obtained using the dimensionless arc-length parametrization. The results obtained using the dimensionless area parametrization, which were consistent with what was observed using the arc-length parametrization, are presented in Appendix A.

5.1. Framework for numerical sensitivity analysis of the system

To facilitate applying the theoretical considerations in Section 4 to our system of ODEs, we relabel the variables as follows:

x_{1} = x, x_{2} = y, x_{3} = ψ,

(45a)

x_{4} = h, x_{5} = l, x_{6} = \tilde{λ} .

(45b)

Using these new labels, we may rewrite our boundary value problem and set up the equations for direct (or adjoint) sensitivity analysis. For each of our numerical experiments, in addition to the choice of parametrization (arc-length parametrization or area parametrization), we had to make several other choices, including:

(a) Choice of spontaneous curvature function (explained in the following);

(b) Choice of boundary conditions (explained in the following).

Once we set up the boundary value problem, the MATLAB BVP solver “bvp4c” was used to solve the system numerically. Roughly speaking, the domain is partitioned into subintervals and on each subinterval the solution functions are approximated by polynomials of degree at most three. Note that a cubic polynomial has four coefficients, so to find the approximate solutions, the solver needs to find four coefficients for each unknown function on each subinterval. The equations needed to solve for the unknown coefficients are obtained by requiring that each approximate solution must be continuous over the entire interval, and also requiring that the differential equation hold at certain points on each subinterval (so the method used by “bvp4c” is, in essence, a collocation method).

5.1.1. Case 1: Dimensionless arc-length parametrization

Our original boundary value problem can be written as

F (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, t_{0}) = \vec{0}

(46)

(Boundary conditions I) (\begin{matrix} x_{1} {(0}^{+}) = 0, x_{3} {(0}^{+}) = 0, x_{5} {(0}^{+}) = 0 \\ x_{2} (T) = 0, x_{3} (T) = 0, x_{6} (T) = {\tilde{λ}}_{0} \end{matrix},

(47)

(Boundary conditions II) (\begin{matrix} x_{1} {(0}^{+}) = \sin θ, x_{3} {(0}^{+}) = θ \\ x_{2} (T) = 0, x_{3} (T) = 0, x_{5} (T) = 0, x_{6} (T) = {\tilde{λ}}_{0} \end{matrix},

(48)

where $θ$ is a fixed angle, F is

\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, t_{0}) = [\begin{matrix} f_{1} (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, t_{0}) \\ ⋮ \\ f_{6} (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, t_{0}) \end{matrix}] \\ = [\begin{matrix} {\overset{\cdot}{x}}_{1} - \cos x_{3} \\ {\overset{\cdot}{x}}_{2} - \sin x_{3} \\ {\overset{\cdot}{x}}_{3} - 2 x_{4} + \frac{\sin x_{3}}{x_{1}} \\ {\overset{\cdot}{x}}_{4} - \frac{x_{5}}{x_{1}} - \overset{\cdot}{c} \\ {\overset{\cdot}{x}}_{5} - 2 x_{1} x_{4} [(x_{4} - c)^{2} + \frac{x_{6}}{\tilde{κ}}] + 2 x_{1} (x_{4} - c) [x_{4}^{2} + {(x_{4} - \frac{\sin x_{3}}{x_{1}})}^{2}] \\ {\overset{\cdot}{x}}_{6} - 2 \tilde{κ} \overset{\cdot}{c} x_{4} + 2 \tilde{κ} c \overset{\cdot}{c} \end{matrix}], \end{matrix}

(49)

and the function $c (t)$ represents the spontaneous curvature, which, in the existing literature, as mentioned before, is chosen to have one of the following forms:

Type I : c (t) = 0.5 R_{0} C_{0} [1 - \tanh [γ (t - t_{0})]],

(50)

Type II : c (t) = - 0.5 R_{0} C_{0} \frac{t - t_{0}}{t_{0}} [1 - \tanh [γ (t - t_{0})]] .

(51)

Remark 7. Note that the graph of $\tanh$ has two horizontal asymptotes. Roughly,

\tanh u \approx (\begin{matrix} 1 & u \geq 3 \\ - 1 & u \leq - 3 \end{matrix},

(52)

and so

1 - \tanh [γ (t - t_{0})] \approx (\begin{matrix} 0 & t \geq t_{0} + \frac{3}{γ} \\ 2 & t \leq t_{0} - \frac{3}{γ} \end{matrix},

(53)

Because of this property, we may choose $γ$ and $t_{0}$ so that the graph of $1 - \tanh [γ (t - t_{0})]$ possesses an abrupt jump near $t_{0}$ . However, if the parameters $γ$ and $t_{0}$ are chosen such that there is no “abrupt” jump in the graph of $1 - \tanh [γ (t - t_{0})]$ for positive values of t (in particular, this happens if $t_{0} >> 0$ and $γ t_{0} < 3$ ), we say our spontaneous curvature function has a mollifying property. As we shall see, whether or not the spontaneous curvature function is mollifying can affect the solution sensitivities.

Remark 8. An explanation as to why the first set of boundary conditions (equation (47)) physically makes sense can be found in [10]. Here is one way to interpret the second set of boundary conditions (equation (48)). We assume that the protein coat consists of two segments; a segment with constant curvature $C_{0}$ and another segment over which the curvature of the coat decreases from $C_{0}$ to 0. Clearly, we do not need to consider the constant curvature portion as part of our domain (because the shape of this part is already known). The prescribed value of $x_{1} (0^{+})$ should be thought of as the distance from the $x_{2}$ -axis marking the end of the a priori known spherical part of the coated membrane (with constant curvature $C_{0}$ ) and the start of the part of the membrane whose shape we want to determine (the unknown part of the membrane). By prescribing a value for $x_{3} (0^{+})$ we set the “direction” of the generating curve at the starting point of the unknown part of the membrane. When we use Type II boundary conditions, our goal is to determine the shape of the part of the membrane that we do not know a priori.

Remark 9. Although both types of boundary condition, equations (47) and (48), are taken from the existing literature, it seems to us that each type comes with certain deficiencies that we believe should be clearly discussed.

In the first set of boundary conditions, $x_{1} (0^{+}) = 0$ and $x_{5} (0^{+}) = 0$ correspond to $r (0^{+}) = 0$ and $L (0^{+}) = 0$ , respectively. However, $L (s) = r (s) [H (s) - C (s)]'$ . So it seems that these equalities are not independent. In practice, for all of our simulations the value of $r (0^{+})$ is set to be a number close to zero rather than zero, and it is scaled appropriately for the dimensionless parametrizations. This is done primarily to avoid numerical difficulties that may arise as a result of the presence of the term $1 / r (s)$ in the equations. As a byproduct, this will also take care of the dependency described previously.

In the second set of boundary conditions, the side condition $x_{5} (0^{+}) = 0$ is replaced by $x_{5} (T) = 0$ . Hence, the aforementioned dependency does not exist anymore. Nevertheless, it seems to us that the prescribed value of $x_{1} (0^{+})$ is rather arbitrary. As explained before, here the goal was to have a nonzero value for $x_{1} (0^{+})$ marking the end of the spherical part of the coat; however, we could not find any convincing geometric justification as to why $\sin θ$ is a good choice for this nonzero value.

Although we were aware of these issues, in order to make our results relatable to the existing literature on the subject, we decided not to depart from the standard boundary conditions used in the literature.

The energy functional of interest, that is, the total elastic bending energy of the membrane can be represented by

W (\vec{x}, C_{0}, γ, α_{0}) ~ \int_{0}^{T} [x_{4} (t) - c (t {)]}^{2} x_{1} (t) d t .

(54)

Indeed, if we use the change of variable $s = R_{0} t$ , then

\begin{matrix} \int_{0}^{T} [h (t) - c (t {)]}^{2} x (t) d t = \int_{0}^{T} R_{0}^{2} [H (R_{0} t) - C (R_{0} t {)]}^{2} \frac{r (R_{0} t)}{R_{0}} d t \\ = \int_{0}^{R_{0} T} R_{0} [H (s) - C (s {)]}^{2} r (s) \frac{d s}{R_{0}} \\ = \int_{0}^{S} [H (s) - C (s {)]}^{2} r (s) d s ~ W . \end{matrix}

(55)

To ensure that our notation is consistent with that discussed in Section 4, we also introduce $p_{1}$ , $p_{2}$ , and $p_{3}$ as follows:

p_{1} = C_{0}, p_{2} = γ, p_{3} = t_{0} .

(56)

Following the discussion in Section 4, for each $1 \leq j \leq 3$ , to find the sensitivities ${\vec{s}}_{j} = \partial \vec{x} / \partial p_{j}$ , it is enough to solve the following system of ODEs for the unknowns $\vec{x} (t)$ and ${\vec{s}}_{j} (t)$ :

(\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, p_{1}, p_{2}, p_{3}) = \vec{0} & (six scalar ODEs) \\ (D_{\overset{\cdot}{\vec{x}}} F) {\overset{\cdot}{\vec{s}}}_{j} + (D_{\vec{x}} F) {\vec{s}}_{j} + \frac{\partial F}{\partial p_{j}} = \vec{0} & (six scalar ODEs) . \end{matrix}

(57)

The boundary conditions for the sensitivities can be obtained by taking the derivative of the boundary conditions for the original unknowns $x_{1}, \dots, x_{6}$ . For example, if we decide to use the first type of boundary conditions, then we will equip this system of ODEs with the following boundary conditions:

(\begin{matrix} x_{1} {(0}^{+}) = 0, x_{3} {(0}^{+}) = 0, x_{5} {(0}^{+}) = 0, \\ x_{2} (T) = 0, x_{3} (T) = 0, x_{6} (T) = {\tilde{λ}}_{0}, \\ s_{j 1} {(0}^{+}) = 0, s_{j 3} {(0}^{+}) = 0, s_{j 5} {(0}^{+}) = 0, \\ s_{j 2} (T) = 0, s_{j 3} (T) = 0, s_{j 6} (T) = 0 . \end{matrix}

(58)

Here, $s_{ji}$ denotes the ith component of the vector ${\vec{s}}_{j}$ , that is, $s_{ji} (t) = \partial x_{i} (t) / \partial p_{j}$ . Once we find each ${\vec{s}}_{j}$ , we can use the following formula to compute the sensitivities of our energy functional:

\forall 1 \leq j \leq 3 \frac{\partial W}{\partial p_{j}} ~ \int_{0}^{T} \frac{\partial w}{\partial p_{j}} + (D_{\vec{x}} w) {\vec{s}}_{j} d t,

(59)

where $w (\vec{x}, p_{1}, p_{2}, p_{3}) = [x_{4} (t) - c (t {)]}^{2} x_{1} (t)$ . Sometimes we denote the integrand in this integral by $Dw / D p_{j}$ .

5.1.2. Case 2: Dimensionless area parametrization

Our original boundary value problem can be written as

F (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, α_{0}) = \vec{0}

(60)

(Boundary conditions I) (\begin{matrix} x_{1} {(0}^{+}) = 0, x_{3} {(0}^{+}) = 0, x_{5} {(0}^{+}) = 0 \\ x_{2} (α_{max}) = 0, x_{3} (α_{max}) = 0, x_{6} (α_{max}) = {\tilde{λ}}_{0} \end{matrix},

(61)

$or$

(Boundary conditions II) (\begin{matrix} x_{1} {(0}^{+}) = \sin θ, x_{3} {(0}^{+}) = θ \\ x_{2} (α_{max}) = 0, x_{3} (α_{max}) = 0, x_{5} (α_{max}) = 0, x_{6} (α_{max}) = {\tilde{λ}}_{0} \end{matrix},

(62)

where $θ$ is a fixed angle and F is

\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, α_{0}) = [\begin{matrix} f_{1} (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, α_{0}) \\ ⋮ \\ f_{6} (\vec{x}, \overset{\cdot}{\vec{x}}, C_{0}, γ, α_{0}) \end{matrix}] \\ = [\begin{matrix} {\overset{\cdot}{x}}_{1} - \frac{\cos x_{3}}{x_{1}} \\ {\overset{\cdot}{x}}_{2} - \frac{\sin x_{3}}{x_{1}} \\ {\overset{\cdot}{x}}_{3} - 2 \frac{x_{4}}{x_{1}} + \frac{\sin x_{3}}{x_{1}^{2}} \\ {\overset{\cdot}{x}}_{4} - \frac{x_{5}}{x_{1}^{2}} - \overset{\cdot}{c} \\ {\overset{\cdot}{x}}_{5} - 2 x_{4} [(x_{4} - c)^{2} + \frac{x_{6}}{\tilde{κ}}] + 2 (x_{4} - c) [x_{4}^{2} + {(x_{4} - \frac{\sin x_{3}}{x_{1}})}^{2}] \\ {\overset{\cdot}{x}}_{6} - 2 \tilde{κ} \overset{\cdot}{c} x_{4} + 2 \tilde{κ} c \overset{\cdot}{c} \end{matrix}] . \end{matrix}

(63)

Here, the function $c (α)$ represents the spontaneous curvature, which, as discussed in previous sections, is chosen to have one of the following forms:

Type I : c (α) = 0.5 R_{0} C_{0} [1 - \tanh [γ (α - α_{0})]],

(64)

Type II : c (α) = - 0.5 R_{0} C_{0} \frac{α - α_{0}}{α_{0}} [1 - \tanh [γ (α - α_{0})]] .

(65)

The energy functional of interest, that is, the total elastic bending energy of the membrane, can be represented by

W (\vec{x}, C_{0}, γ, α_{0}) ~ \int_{0}^{α_{max}} [x_{4} (α) - c (α {)]}^{2} d α .

(66)

Indeed, if we use the change of variable $a = 2 π R_{0}^{2} α$ , then

\begin{matrix} \int_{0}^{α_{max}} [h (α) - c (α {)]}^{2} d α = \int_{0}^{α_{max}} R_{0}^{2} [H (2 π R_{0}^{2} α) - C (2 π R_{0}^{2} α {)]}^{2} d α \\ = \int_{0}^{2 π R_{0}^{2} α_{max}} R_{0}^{2} [H (a) - C (a {)]}^{2} \frac{d a}{2 π R_{0}^{2}} \\ = \frac{1}{2 π} \int_{0}^{A} [H (a) - C (a {)]}^{2} d a ~ W . \end{matrix}

(67)

Again, to ensure that our notation is consistent with that discussed in Section 4, we also introduce $p_{1}$ , $p_{2}$ , and $p_{3}$ , as

p_{1} = C_{0}, p_{2} = γ, p_{3} = α_{0} .

(68)

(\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, p_{1}, p_{2}, p_{3}) = \vec{0} & (six scalar ODEs) \\ (D_{\overset{\cdot}{\vec{x}}} F) {\overset{\cdot}{\vec{s}}}_{j} + (D_{\vec{x}} F) {\vec{s}}_{j} + \frac{\partial F}{\partial p_{j}} = \vec{0} & (six scalar ODEs) \end{matrix} .

(69)

The boundary conditions for the sensitivities can be obtained by taking the derivative of the boundary conditions for the original unknowns $x_{1}, \dots, x_{6}$ . For example, if we decide to use the second type of boundary conditions, then we will equip this system of ODEs with the following boundary conditions:

(\begin{matrix} x_{1} {(0}^{+}) = \sin θ, x_{3} {(0}^{+}) = θ, \\ x_{2} (α_{max}) = 0, x_{3} (α_{max}) = 0, x_{5} (α_{max}) = 0, x_{6} (α_{max}) = {\tilde{λ}}_{0}, \\ s_{j 1} {(0}^{+}) = 0, s_{j 3} {(0}^{+}) = 0, \\ s_{j 2} (α_{max}) = 0, s_{j 3} (α_{max}) = 0, s_{j 5} (α_{max}) = 0, s_{j 6} (α_{max}) = 0 . \end{matrix}

(70)

As before $s_{ji}$ denotes the ith component of the vector ${\vec{s}}_{j}$ , that is, $s_{ji} (α) = \partial x_{i} (α) / \partial p_{j}$ . Once we find each ${\vec{s}}_{j}$ , we can use the following formula to compute the sensitivities of our energy functional:

\forall 1 \leq j \leq 3 \frac{\partial W}{\partial p_{j}} ~ \int_{0}^{α_{max}} \frac{\partial w}{\partial p_{j}} + (D_{\vec{x}} w) {\vec{s}}_{j} d α,

(71)

where $w (\vec{x}, p_{1}, p_{2}, p_{3}) = [x_{4} (α) - c (α {)]}^{2}$ . The integrand in this integral is sometimes denoted by $D w / D p_{j}$ .

Remark 10. Alternatively, according to the adjoint method, to compute the sensitivities of the energy functional W, we may use the following formula:

\frac{\partial W}{\partial p_{j}} ~ \int_{0}^{u_{max}} (\frac{\partial w}{\partial p_{j}} - {\vec{v}}^{T} \frac{\partial F}{\partial p_{j}}) d u,

(72)

where, for example, in the case where we use the dimensionless area parametrization with boundary conditions of Type I, $u_{max}$ and u in this formula stand for $α$ and $α_{max}$ , respectively, and the vector function

\vec{v} = [\begin{matrix} v_{1} \\ ⋮ \\ v_{6} \end{matrix}]

can be computed by finding a solution of the following system of ODEs (with 12 scalar unknown functions):

(\begin{matrix} F (\vec{x}, \overset{\cdot}{\vec{x}}, p_{1}, p_{2}, p_{3}) = \vec{0}, \\ \frac{d}{d α} [{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F] - {\vec{v}}^{T} D_{\vec{x}} F = - D_{\vec{x}} w, \\ x_{1} {(0}^{+}) = 0, x_{3} {(0}^{+}) = 0, x_{5} {(0}^{+}) = 0, \\ x_{2} (α_{max}) = 0, x_{3} (α_{max}) = 0, x_{6} (α_{max}) = {\tilde{λ}}_{0}, \\ v_{2} {(0}^{+}) = 0, v_{4} {(0}^{+}) = 0, v_{6} {(0}^{+}) = 0, \\ v_{1} (α_{max}) = 0, v_{4} (α_{max}) = 0, v_{5} (α_{max}) = 0 . \end{matrix}

The following proposition justifies the particular choice of boundary conditions for certain components of $\vec{v}$ in this boundary value problem.

Proposition 2. Let $\vec{x}$ and $F$ be as defined previously. Suppose that

\begin{matrix} x_{1} (0^{+}) = 0, & x_{3} (0^{+}) = 0, & x_{5} (0^{+}) = 0, \\ x_{2} (α_{max}) = 0, & x_{3} (α_{max}) = 0, & x_{6} (α_{max}) = {\tilde{λ}}_{0} . \end{matrix}

\vec{v} (α) = [\begin{matrix} v_{1} (α) \\ ⋮ \\ v_{6} (α) \end{matrix}]

satisfies

\begin{matrix} v_{2} (0^{+}) = 0, & v_{4} (0^{+}) = 0, & v_{6} (0^{+}) = 0, \\ v_{1} (α_{max}) = 0, & v_{4} (α_{max}) = 0, & v_{5} (α_{max}) = 0, \end{matrix}

then

{[{\vec{v}}^{T} D_{\overset{\cdot}{\vec{x}}} F \frac{\partial \vec{x}}{\partial p_{j}}]}_{α = 0}^{α = α_{max}} = 0 .

Proof. At $α = 0^{+}$ , we have

\begin{matrix} {\vec{v}}^{T} (0^{+}) D_{\overset{\cdot}{\vec{x}}} F |_{α {= 0}^{+}} \frac{\partial \vec{x} {(0}^{+})}{\partial p_{j}} = [\begin{matrix} v_{1} {(0}^{+}) & \dots & v_{6} {(0}^{+}) \end{matrix}] [\begin{matrix} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{1}} |_{α {= 0}^{+}} & \dots & \frac{\partial F}{\partial {\overset{\cdot}{x}}_{6}} |_{α {= 0}^{+}} \end{matrix}] [\begin{matrix} 0 \\ \frac{x_{2} {(0}^{+})}{\partial p_{j}} \\ 0 \\ \frac{x_{4} {(0}^{+})}{\partial p_{j}} \\ 0 \\ \frac{x_{6} {(0}^{+})}{\partial p_{j}} \end{matrix}] \\ = [\begin{matrix} v_{1} {(0}^{+}) & \dots & v_{6} {(0}^{+}) \end{matrix}] (\frac{\partial x_{2} {(0}^{+})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{2}} |_{α {= 0}^{+}} + \frac{\partial x_{4} {(0}^{+})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{4}} |_{α {= 0}^{+}} + \frac{\partial x_{6} {(0}^{+})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{6}} |_{α {= 0}^{+}}) \\ = \frac{\partial x_{2} {(0}^{+})}{\partial p_{j}} v_{2} (0^{+}) + \frac{\partial x_{4} {(0}^{+})}{\partial p_{j}} v_{4} (0^{+}) + \frac{\partial x_{6} {(0}^{+})}{\partial p_{j}} v_{6} (0^{+}) = 0 . \end{matrix}

Here, we used $v_{2} (0^{+}) = v_{4} (0^{+}) = v_{6} (0^{+}) = 0$ and the fact that $D_{\overset{\cdot}{\vec{x}}} F$ is the identity matrix.

Similarly, at $α = α_{max}$ , we have

\begin{matrix} {\vec{v}}^{T} (α_{max}) D_{\overset{\cdot}{\vec{x}}} F |_{α = α_{max}} \frac{\partial \vec{x} (α_{max})}{\partial p_{j}} \\ = [\begin{matrix} v_{1} (α_{max}) & \dots & v_{6} (α_{max}) \end{matrix}] [\begin{matrix} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{1}} |_{α = α_{max}} & \dots & \frac{\partial F}{\partial {\overset{\cdot}{x}}_{6}} |_{α = α_{max}} \end{matrix}] [\begin{matrix} \frac{x_{1} (α_{max})}{\partial p_{j}} \\ 0 \\ 0 \\ \frac{x_{4} (α_{max})}{\partial p_{j}} \\ \frac{x_{5} (α_{max})}{\partial p_{j}} \\ 0 \end{matrix}] \\ = [\begin{matrix} v_{1} (α_{max}) & \dots & v_{6} (α_{max}) \end{matrix}] (\frac{\partial x_{1} (α_{max})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{1}} |_{α = α_{max}} + \frac{\partial x_{4} (α_{max})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{4}} |_{α = α_{max}} + \frac{\partial x_{5} (α_{max})}{\partial p_{j}} \frac{\partial F}{\partial {\overset{\cdot}{x}}_{5}} |_{α = α_{max}}) \\ = \frac{\partial x_{1} (α_{max})}{\partial p_{j}} v_{1} (α_{max}) + \frac{\partial x_{4} (α_{max})}{\partial p_{j}} v_{4} (α_{max}) + \frac{\partial x_{5} (α_{max})}{\partial p_{j}} v_{5} (α_{max}) = 0 . \end{matrix}

Here, we used $v_{1} (α_{max}) = v_{4} (α_{max}) = v_{5} (α_{max}) = 0$ and the fact that $D_{\overset{\cdot}{\vec{x}}} F$ is the identity matrix.

5.2. Dimensionless arc-length formulation—Type I spontaneous curvature

The spontaneous curvature function used is

c (t) = 0.5 C_{0} R_{0} [1 - \tanh [γ (t - t 0)]] .

(73)

Note that, for each choice of $C_{0}$ , $c (t)$ is a constant multiple of $F (t) = 1 - \tanh [γ (t - t_{0})]$ , so $F (t)$ is the function that determines the shape of the spontaneous curvature. The graph of $F (t)$ is depicted in Figure 2.

Boundary conditions used for the original unknowns are:

(\begin{matrix} x_{1} {(0}^{+}) = 0, x_{3} {(0}^{+}) = 0, x_{5} {(0}^{+}) = 0, \\ x_{2} (T) = 0, x_{3} (T) = 0, x_{6} (T) = {\tilde{λ}}_{0} . \end{matrix}

(74)

The boundary conditions for the corresponding sensitivities are all set to be zero.

Our choices for the input parameters are given in Table 1. The value of $R_{0}$ is chosen such that $t_{0} = 1$ . A brief explanation of the choice $s_{0} = 80$ and $T = 5$ is given in Section 6. The value of $γ$ is chosen such that there will be a sharp transition at $t_{0} = 1$ , while $κ_{0}$ is chosen to be equal to $κ$ . (A range of acceptable values for $κ$ and $λ_{0}$ can be found in [10].)

Table 1.

Parameters used in the model. The corresponding diagrams are depicted in Figure 2.

Parameter	Input used
$R_{0}$ (nm)	80
$C_{0}$ (nm⁻¹)	Various values in [0,0.02]
$γ = ξ * R_{0}$	20
$t_{0} = s_{0} / R_{0}$	1
${\tilde{λ}}_{0} = λ_{0} R_{0}^{2} / κ_{0}$	12.8
Domain $= [0, T]$	[0,5]

We organize our results as follows: graphs of energy sensitivities with respect to the parameters $C_{0}$ , $λ$ , and $t_{0}$ are displayed on the left of Figure 2(b) to (d), respectively; the graphs of the dimensionless curvature sensitivities with respect to the parameters involved are depicted on the right of these panels.

Here we make two key observations. First, we notice that for our choice of the type of boundary conditions (Type I, equation (47)) and spontaneous curvature function, the energy sensitivity graphs do not intersect the horizontal axis (that is, there is no critical point). As we shall see later, this seems to be in correlation with certain interesting properties of the final shape of the membrane. Second, notice that for our choice of boundary conditions (Type I, equation (47)) and spontaneous curvature function (which possesses a sharp transition), there is an abrupt change in curvature sensitivity near the location where the sharp transition in spontaneous curvature function occurs, which, of course, is expected. As we shall soon see, the graphs of curvature sensitivities will smear out if we smooth the transition in the spontaneous curvature function.

5.3. Dimensionless arc-length formulation—mollifying Type I spontaneous curvature

All parameters are exactly the same as Section 5.2 except $γ$ . This time we choose $γ$ to be equal to $2.5$ so that $s_{0} γ < 3$ . This choice of $γ$ and $s_{0}$ smooths the transition in the spontaneous curvature function (see Figure 3(a)). As a result, we will see that the graphs of curvature sensitivities (the right diagrams in Figure 3(b) to (d)) are smeared out (compared with the observations in Section 5.2). Notice that this particular change in the shape of the spontaneous curvature function does not have any effect on the number of horizontal intercepts of the energy sensitivities.

5.4. Dimensionless arc-length formulation—Type II spontaneous curvature

The spontaneous curvature function used is

c (t) = - 0.5 C_{0} R_{0} \frac{t - t_{0}}{t_{0}} [1 - \tanh [γ (t - t 0)]] .

(75)

Note that for each choice of $C_{0}$ , $c (t)$ is a constant multiple of

F (t) = - \frac{t - t_{0}}{t_{0}} [1 - \tanh [γ (t - t_{0})]] .

The graph of $F (t)$ is depicted in Figure 4(a). The boundary conditions used for the original unknowns are:

(\begin{matrix} x_{1} {(0}^{+}) = \sin (0.9 π), x_{3} {(0}^{+}) = 0.9 π, \\ x_{2} (T) = 0, x_{3} (T) = 0, x_{5} (T) = 0, x_{6} (T) = 0 . \end{matrix}

(76)

The boundary conditions for the corresponding sensitivities are all set to be zero.

Our choices for the input parameters are given in Table 2. This particular choice of parameters is in agreement with [24].

Table 2.

Parameters used in the model. The corresponding diagrams are depicted in Figure 4.

Parameter	Input used
$R_{0}$ (nm)	200
$C_{0}$ (nm⁻¹)	Various values in $[0, \frac{4.5}{1000}]$
$γ = ξ R_{0}$	30
$t_{0} = \frac{s_{0}}{R_{0}}$	30
Domain $= [0, T]$	[0,100]

We organize our results as follows: graphs of energy sensitivities with respect to parameters $C_{0}$ , $λ$ , and $t_{0}$ are displayed in the left panels of Figure 4(b) to (d), respectively; the graphs of the dimensionless curvature sensitivities with respect to parameters involved are depicted on the right of these panels. We make two observations. First, notice that, for our choice of the type of boundary conditions (Type II, equation (48)) and spontaneous curvature function, the energy sensitivity graph with respect to $C_{0}$ indeed intersects the horizontal axis. As we shall see, this seems to be in correlation with certain interesting properties of the final shape of the membrane. Second, notice that for our choice of boundary conditions and spontaneous curvature function (which possesses a sharp transition), the graphs of curvature sensitivities have sharp bends near the location where the sharp transition in spontaneous curvature function occurs. As we shall soon see, the graphs of curvature sensitivities will be smoother near $t_{0}$ if we smooth the sharp transition in the spontaneous curvature function.

5.5. Dimensionless arc-length formulation—mollifying Type II spontaneous curvature

Almost all parameters are the same as those used in Section 5.4, except for the value chosen for $γ$ . This time we choose $γ$ to be equal to $0.08$ such that $s_{0} γ < 3$ . This choice of $γ$ and $s_{0}$ further smooths the transition in the spontaneous curvature function (see Figure 5(a)). As a result, we will see that the graphs of curvature sensitivities (diagrams on the right of Figure 5(b) to (d)) are smoother near $t_{0} = 30$ (compared with the observation in Section 5.4).

6. Results and conclusions

6.1. Insights obtained from sensitivity analyses on membrane shapes

6.1.1. Remarks on the size of the domain

As discussed earlier, we use a MATLAB BVP solver to approximate the solution to our system of ODEs. For this reason, the dimensionless size of the domain (that is, the interval over which we want to solve the system), the mesh size, and the initial guess for the solution (which is given to the solver as an input) can play key roles in the performance of the method and the accuracy of the results. In each of our numerical experiments, we normally start by dividing the interval into about 100 subintervals and then use a finer mesh if we run into trouble. As expected, we observed that, for larger domain size, finer partitions are needed for a good performance. In our simulations we adhered to the following rules with regard to the total size of the domain:

Dimensionless arc-length formulation. The total size is taken to be 2–5 times the size of the coated region.

Dimensionless area formulation. The total size is taken to be 4–25 times the size of the coated region.

It is important to mention that there are experimental measurements that can provide an estimate of the size of the coat [25]. In some of our simulations, we used $a_{0} = 20106$ nm² as an estimate of the area of the coat [10]. The corresponding dimensionless area is obtained by first choosing a number $R_{0}$ and then computing $a_{0} / 2 π R_{0}^{2}$ . We mention in passing that, although theoretically the value chosen for $R_{0}$ should not have any effect on the solution of the system, in practice, when using a numerical method to solve the problem, the value of $R_{0}$ can affect the performance of the numerical solver. Our chosen values for $R_{0}$ can be found in Tables 1 to 4, displayed in Section 5 and Appendix A.

We know that the area of a spherical cap corresponding to arc length s is approximately equal to $π s^{2}$ (assuming that s is much smaller than the radius of the sphere). For this reason, it is not completely unreasonable to use the equation $π s_{0}^{2} = 20106$ to obtain an estimate of the length corresponding to the coated region (which gives $s_{0} \approx 80$ nm).

6.1.2. Should we expect the area formulation and the arc-length formulation to produce the exact same results?

Let us assume that we have fixed the values of $a_{0}$ , $s_{0}$ , $C_{0}$ , and $ξ$ . An important question is the following: Is it reasonable to expect that we should obtain the exact same answers using the arc-length formulation and the area formulation? A little bit of deliberation tells us that the answer is no! For one thing, what the exact relationship between $a_{0}$ and $s_{0}$ should be is not known in advance. More importantly, although the spontaneous curvature function used for the area formulation has the same form as the spontaneous curvature function used in the arc-length formulation, the area spontaneous curvature used in our simulations is not exactly the function obtained by transforming $C (s) = 0.5 C_{0} F (s)$ into a function of a using the relationship $a = g (s) = 2 π \int_{0}^{s} r (u) d u$ . As discussed in Remark 4, the area spontaneous curvature function that would have produced the same answer would have been $\hat{C} (a) = 0.5 C_{0} F (g^{- 1} (a))$ . For example, in the case where $F (s) = [1 - \tanh [ξ (s - s_{0})]]$ , we would obtain $\hat{C} (a) = 0.5 C_{0} [1 - \tanh [ξ (g^{- 1} (a) - s_{0})]]$ , which, partly owing to the presence of the function $g^{- 1}$ , does not necessarily have the same form as the form that is actually prescribed for $C (a)$ in our simulations.

6.1.3. On the importance of boundary conditions and the value of ${\tilde{λ}}_{0}$

As noted in previous sections, we ran our simulations with two different sets of boundary conditions, representing distinct physical constraints or assumptions. Unfortunately, as opposed to the case of initial value problems, there is no general mathematical theory that can be used to ensure the existence and uniqueness of solutions to boundary value problems such as the one studied in this work. One thing that became clear to us was that the performance of the MATLAB BVP solver “bvp4c” was highly sensitive to the chosen boundary conditions, in particular to the value of ${\tilde{λ}}_{0}$ . In fact, in some cases, a small change in the value of ${\tilde{λ}}_{0}$ can result in the appearance of linear algebraic systems with singular coefficient matrices in the process of numerically solving the ODE system. At this point, it is not clear whether this is a theoretical issue related to the existence of solutions, or whether this should be explained by exploring the stability properties of the numerical method used to solve the system.

6.1.4. A cleverly chosen spontaneous curvature can smear out solution sensitivities

As we mentioned before, whether or not the spontaneous curvature function is mollifying can affect the solution sensitivities. Our simulations provide numerical evidence for the conjecture that the smoother the transition between the nonzero part of $F (u)$ and the zero part of $F (u)$ is, the less likely it will be to have abrupt jumps in the graphs of solution sensitivities. Notice that the shape of the function F over the transition region can be tuned by cleverly choosing the parameters $γ$ and $t_{0}$ (or $α_{0}$ ).

6.1.5. Is there any correlation between sensitivity diagrams and the final shape of the membrane?

One motivation of this work was to gain more insight into the circumstances that would result in bud-shaped membranes, as opposed to those that would give pearl-shaped membranes. To that end, we performed various numerical experiments to see whether we could find evidence indicating positive or negative correlation between the input data of the problem (such as the type of the spontaneous curvature, the parameters used in the expression of the spontaneous curvature, and the type of the boundary conditions) and the final shape of the membrane.

The following observations or conjectures are in agreement with all of our results, parts of which are depicted in Figure 6.

Our results provide numerical evidence for the conjecture that there might be a connection between the behavior of the energy sensitivity with respect to $C_{0}$ and the appearance of turning points in the generating curve of the surface of the membrane and the formation of pearls (loops). In particular, as becomes evident by comparing panels (a), (b), (c), (d), and (g) of Figure 6 with panels (e), (f), and (h), we observed that when W did not have any critical point (as a function of $C_{0}$ ), that is, when the energy sensitivity function with respect to $C_{0}$ did not intersect the horizontal axis, no loops were formed.

Previously we observed that using a mollifying spontaneous curvature function can smear out the solution sensitivities. However, the type of spontaneous curvature function (that is, Type I or Type II) or whether the function is mollifying or not, does not appear to preclude the possibility of formation of pearls.

For a fixed domain size and coat size, we are more likely to see an energy sensitivity function with oscillatory behavior about the horizontal axis (indicating existence of zeros) when we use the second type of boundary conditions. Indeed, in our experiments, every time we used the first set of boundary conditions, we noticed that the resulting energy sensitivity function (as a function of $C_{0}$ ) will not intersect the horizontal axis, and subsequently no pearls were formed.

6.2. Examining a conjecture related to the “pearling” transition

In our discussions with researchers working on the subject, we noticed that some implicitly believe in the conjecture that “The slope of the transition part of the spontaneous curvature function is a key factor in whether or not pearls will be formed; in particular, there is a correlation between formation of pearls and not having sharp transitions (big slopes) in the spontaneous curvature function used in the numerical solution of the boundary value problem.”

The purpose of this section is to provide a simple argument that disproves the validity of this conjecture in the generality stated. Indeed, in what follows, we will show that there exist spontaneous curvature functions with very large slopes on their transition regions that ultimately result in the formation of pearls. To be concrete, we focus on the dimensionless arc-length parametrization; however, an analogous argument can be applied to the dimensionless area parametrization.

Let $(C_{0}, R_{0}, γ, t_{0})$ be a set of parameters with $γ t_{0} > 3$ that results in the formation of pearls in the domain $[0, T]$ . Note that since $γ t_{0} > 3$ , the value of the spontaneous curvature function at 0 is approximately $C_{0} R_{0}$ and it decreases to near zero over the interval $[0, t_{0} + 3 / γ]$ . As discussed in Section 3, if $x (t)$ and $y (t)$ are the first two components of the corresponding solution (note that $R_{0} x (t)$ and $R_{0} y (t)$ are the coordinates of the generating curve of the surface of the membrane), then for any constant $η > 0$ , the functions $x (η t) / η$ and $y (η t) / η$ will be the first two components of the solution in $[0, T / η]$ using parameters $(C_{0}, η R_{0}, η γ, t_{0} / η)$ . In particular, the shape of the solution curve (and the corresponding surface) using these new parameters will be the same as the shape of the solution curve (and the corresponding surface) using the original parameters. If one has pearls, then the other will also have pearls. Now note that $(η γ) t_{0} / η = γ t_{0} > 3$ and so the value of the corresponding spontaneous curvature at 0 is approximately $(C_{0}) (η R_{0})$ and the spontaneous curvature function decreases to nearly zero over the interval $[0, t_{0} / η + 3 / η γ]$ . We can set the number $η$ to be as large as we want. The larger $η$ , the smaller the interval of transition will be and the larger the starting value $η C_{0} R_{0}$ will be. That is, for very large $η$ , the corresponding spontaneous curvature function must decrease from the huge value $η C_{0} R_{0}$ to zero over the very small interval $[0, t_{0} / η + 3 / η γ]$ , which means that it will have a huge slope. Nevertheless, since the original parameters resulted in the formation of pearls, this new spontaneous curvature function will also result in the formation of pearls.

6.3. Concluding remarks

In this study, we conducted sensitivity analysis on the well-known Helfrich model for lipid bilayer bending in the context of the spontaneous curvature function. We observed some interesting phenomena in our numerical experiments that relate the sensitivity of the energy with respect to the free parameters in the spontaneous curvature function to the shape of the membrane. Given the wide usage of the Helfrich model for simulating membrane-bending phenomena that have been reported experimentally, our approach of using sensitivity analysis can provide some insight into how one can design input functions and parameters for this system of ODEs. Our experiments also led us to make certain new, arguably nonobvious, conjectures about the behavior of the solution. Clearly, even if results of numerical experiments using a million different sets of inputs display positive correlation between certain output variables, that does not necessarily mean that theoretically there should be positive correlation between the output variables, regardless of what the inputs are. Nevertheless, such results may help us make conjectures that were not immediately obvious from the outset, and this is, in fact, a major way that science progresses. By no means should this work be viewed as the final word on the parametric sensitivity analysis of the shape equations in the Helfrich energy model. Rather, we suggest that it is merely a first step toward better understanding the role of parameters involved in shape equations, particularly in the context of numerical simulations. Indeed, future directions include identifying suitable “bump” functions for the spontaneous curvature function and extension of these methods to solutions in general coordinates.

Footnotes

A. Some numerical results for dimensionless area parametrization

In this appendix we will present some of our numerical results obtained using the dimensionless area parametrization.

Acknowledgements

The authors thank Haleh Alimohamadi, Jennifer Fromm, Christopher Lee, Can Uysalel, Ritvik Vasan, Cuncheng Zhu, and other members of the Rangamani group for critical discussions and feedback on the manuscript.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by NSF DMS/CM (grant number 1620366), NSF DMS/MB (grant number 1934411), the National Institutes of Health (grant number R01GM132106), and the Air Force Office of Scientific Research (grant number FA9550-18-1-0051).

ORCID iD

P Rangamani

References

Canham

PB.

The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J Theor Biol 1970; 26(1): 61–81.

Helfrich

Elastic properties of lipid bilayers: Theory and possible experiments. Z Naturforsch, C: Biosci 1973; 28(11–12): 693–703.

Alimohamadi

Smith

Nowak

, et al. Non-uniform distribution of myosin-mediated forces governs red blood cell membrane curvature through tension modulation. bioRxiv 2019; 668582.

Turlier

Fedosov

Audoly

, et al. Equilibrium physics breakdown reveals the active nature of red blood cell flickering. Nat Phys 2016; 12(5): 513–519.

Seifert

Configurations of fluid membranes and vesicles. Adv Phys 1997; 46(1): 13–137.

Miao

Seifert

Wortis

, et al. Budding transitions of fluid-bilayer vesicles: The effect of area-difference elasticity. Phys Rev E 1994; 49(6): 5389.

Mouritsen

Bloom

Mattress model of lipid–protein interactions in membranes. Biophys J 1984; 46(2): 141–153.

Alimohamadi

Rangamani

Modeling membrane curvature generation due to membrane–protein interactions. Biomolecules 2018; 8(4): 120.

Kozlov

Campelo

Liska

, et al. Mechanisms shaping cell membranes. Curr Opin Cell Biol 2014; 29: 53–60.

10.

Hassinger

Oster

Drubin

, et al. Design principles for robust vesiculation in clathrin-mediated endocytosis. Proc Natl Acad Sci, USA 2017; 114(7): E1118–E1127.

11.

Agrawal

Steigmann

DJ.

Modeling protein-mediated morphology in biomembranes. Biomech Model Mechanobiol 2009; 8(5): 371–379.

12.

Vasan

Rudraraju

Akamatsu

, et al. A mechanical model reveals that non-axisymmetric buckling lowers the energy barrier associated with membrane neck constriction. Soft Matter 2020; 16: 784–797.

13.

Sauer

Duong

Corbett

CJ.

A computational formulation for constrained solid and liquid membranes considering isogeometric finite elements. Comput Methods Appl Mech Eng 2014; 271: 48–68.

14.

Luke

JC.

A method for the calculation of vesicle shapes. SIAM J Appl Math 1982; 42(2): 333–345.

15.

Walani

Torres

Agrawal

Endocytic proteins drive vesicle growth via instability in high membrane tension environment. Proc Natl Acad Sci, USA 2015; 112(12): E1423–E1432.

16.

Waugh

RE.

Surface viscosity measurements from large bilayer vesicle tether formation. II. Experiments. Biophys J 1982; 38(1): 29–37.

17.

Varma

Morbidelli

Parametric sensitivity in chemical systems. Cambridge: Cambridge University Press, 1999.

18.

Steigmann

Fluid films with curvature elasticity. Arch Ration Mech Anal 1999; 150(2): 127–152.

19.

Agrawal

Steigmann

DJ.

Boundary-value problems in the theory of lipid membranes. Continuum Mech Thermodyn 2009; 21(1): 57–82.

20.

Rangamani

Agrawal

Mandadapu

, et al. Interaction between surface shape and intra-surface viscous flow on lipid membranes. Biomech Model Mechanobiol 2013; 12: 833–845.

21.

Rangamani

Mandadap

Oster

Protein-induced membrane curvature alters local membrane tension. Biophys J 2014; 107(3): 751–762.

22.

Gray

Abbena

Salamon

Modern differential geometry of curves and surfaces with Mathematica. 3rd ed. Boca Raton, FL: Chapman and Hall/CRC, 2006.

23.

Petzold

Cao

, et al. Sensitivity analysis of differential-algebraic equations and partial differential equations. Comput Chem Eng 2006; 30: 1553–1559.

24.

Yuan

Alimohamadi

Bakka

, et al. Membrane bending by protein phase separation. bioRxiv 2020; 10.1101/2020.05.21.109751.

25.

Nickels

Chatterjee

Stanley

, et al. The in vivo structure of biological membranes and evidence for lipid domains. PLoS Biol 2017; 15(5): e2002214.

Local sensitivity analysis of the “membrane shape equation” derived from the Helfrich energy

Abstract

Keywords

1. Introduction

2. Overview of the Helfrich energy with spontaneous curvature

2.1. Equations of motion

2.2. Choice of spontaneous curvature function

3. Axisymmetric parametrization

3.1. Arc-length formulation

3.2. Dimensionless arc-length formulation

3.3. Area formulation

3.4. Dimensionless area formulation

4. Mathematical framework for sensitivity analysis

5. Numerical sensitivity analysis of the Helfrich model

5.1. Framework for numerical sensitivity analysis of the system

5.1.1. Case 1: Dimensionless arc-length parametrization

5.1.2. Case 2: Dimensionless area parametrization

5.2. Dimensionless arc-length formulation—Type I spontaneous curvature

5.3. Dimensionless arc-length formulation—mollifying Type I spontaneous curvature

5.4. Dimensionless arc-length formulation—Type II spontaneous curvature

5.5. Dimensionless arc-length formulation—mollifying Type II spontaneous curvature

6. Results and conclusions

6.1. Insights obtained from sensitivity analyses on membrane shapes

6.1.1. Remarks on the size of the domain

6.1.2. Should we expect the area formulation and the arc-length formulation to produce the exact same results?

6.1.3. On the importance of boundary conditions and the value of λ ~ 0

6.1.4. A cleverly chosen spontaneous curvature can smear out solution sensitivities

6.1.5. Is there any correlation between sensitivity diagrams and the final shape of the membrane?

6.2. Examining a conjecture related to the “pearling” transition

6.3. Concluding remarks

Footnotes

A. Some numerical results for dimensionless area parametrization

Acknowledgements

Funding

ORCID iD

References

6.1.3. On the importance of boundary conditions and the value of ${\tilde{λ}}_{0}$