The basic building block of biological membranes is the lipid bilayer. 1 The compositional asymmetry of such biomembranes, that is, the difference in lipid species found in the outer and inner leaflets, has been studied for half a century 2 and is widely conserved throughout Eukarya. [3][4][5] More recently, it has been shown that human erythrocyte plasma membranes (and likely many other mammalian plasma membranes) are highly asymmetric in terms of overall phospholipid abundance as well. 6 Going along with this, a different kind of asymmetry-that of a mechanical stress difference, or differential stress, between the two leaflets-has drawn increasing interest for its impact on membrane properties and potential biological implications. [7][8][9][10] Differential stress has, among other things, been proposed as a mechanism by which curvature torques originating from lipid shape preference can be canceled out, resulting in a flat membrane. 7 While it has been proposed that frequently flip-flopping species like cholesterol would act to nullify any such differential stress, 11 it was recently shown that this need not be the case; instead, cholesterol may well create differential stress due to preferential lipid interactions. 8 Thus, the combined influence of many distinct membrane asymmetries determines important properties of the bilayer, such as its equilibrium shape, its cholesterol distribution, and its elastic moduli. In this work, we will investigate some peculiar aspects of the simulation of such asymmetric systems.
A standard technique for in silico investigations of lipid bilayers is Molecular Dynamics (MD) simulations. In order to simulate quasi-infinite continuous membrane systems, a common choice of boundary conditions (BCs) is that of fully periodic boundary conditions (PBCs). This choice inadvertently comes with the side effect of enforced membrane flatness; the boundary conditions essentially pin the bilayer into a planar configuration regardless of its preferred curvature state. There are numerous reasons this could be undesirable, for instance if one is interested in curvature induction and/or sensing by proteins interacting with membranes in their elastic ground state. An alternative choice for periodically connecting a bilayer are so-called P2 1 boundary conditions. 12 These have recently been shown to offer some distinct advantages when dealing
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp with asymmetric bilayers, 13 but unfortunately they are not readily available in most MD simulation packages. With this in mind, we propose an alternative set of simulation conditions under which the membrane is allowed to relax its mean curvature. In order to achieve this, we break the periodicity of the membrane along one of the lateral dimensions, resulting in a semi-infinite membrane strip with free edges running parallel to the remaining direction of periodicity. Ordinarily, such open edges on a lipid bilayer would result in highly accelerated lipid flip-flop, yielding an on-average symmetric membrane with zero curvature preference. To circumvent this situation and maintain any conceivably imposed membrane asymmetry, we introduce nano-scale "sticky tapes" which adhere to the membrane open edges and prevent flip-flop over the edge defect.
The precise physical nature of the adhesive patches is not of particular importance, as their purpose is to facilitate novel simulation conditions, not to serve as a template for an experimentally realizable molecular system. In this work, we illustrate the design principles in the ultra-coarse-grained (CG) Cooke lipid model, but the idea is readily transferable to other, more finely resolved models.
Before we discuss the details of this new method, let us first set the stage and revisit in Sec. I A some basic elasticity background for asymmetric membranes and also explain in Sec. I B why periodicity is such a constraining condition.
To the lowest order, lipid membrane curvature elasticity is well modeled by the Helfrich energy functional, 14
In this expression, the integral is taken over a two-dimensional surface S representing the membrane shape. The constants κ and κ are the bending modulus and Gaussian curvature modulus, which respectively quantify the energetic penalty for inducing local curvature K and Gaussian curvature K G . The constant K 0 is the spontaneous curvature the membrane would prefer to have, and it is only nonzero in cases of broken up-down symmetry. In this work, we will be concerned with situations in which neither membrane topology nor the geodesic curvature of any open boundary is changing, and as such we can disregard the second term by invoking the Gauss-Bonnet Theorem. 15 If a finite patch of membrane is stretched or compressed such that its area A differs from its rest area A 0 , a Hookean contribution is added to the free energy,
in which KA is the area or stretching modulus.
The sum E H + EA can describe the elastic free energy of a lipid bilayer membrane, or each constituent monolayer of the bilayer. 16 The second approach allows one to quantify the moduli of each monolayer individually (indicated by a subscript "m," e.g., "κm"), which will for instance depend upon the lipid species present in each layer. The composite bilayer energy is then the sum of the individual monolayer terms, neglecting contributions due to inter-leaflet FIG. 1. Illustration of lipid bilayer geometry. The reference surface for each leaflet (dashed curves) is displaced away from the bilayer midsurface (solid curve) along its normal by distance z±.
coupling. In all that follows, we arbitrarily label one of the monolayers as the "upper" leaflet, and indicate its relevant quantities with a subscript "+," and similarly use a subscript "-" for quantities pertaining to the "lower" leaflet. Quantities with no subscript refer to the composite bilayer as a whole or its midsurface, as appropriate.
It must be noted that our expression for the free energy implicitly assumes that there is no curvature-area cross-coupling term proportional to (K -K 0 ) ⋅ (A -A 0 ), which should reasonably appear in a second-order expansion of the free energy in terms of K and A. The vanishing of this term implies that we take the neutral surface of each monolayer as our reference surface describing its geometry, as by definition this is the surface at which bending and stretching contribute to the free energy independently. 17 The locations of these reference surfaces will be assumed to be a constant distance z± away from the bilayer midplane (see Fig. 1). At times where it is necessary to distinguish the neutral surface from other possible reference surfaces, we will denote it by zn.
As has been shown previously, 7,10 the total bilayer elastic energy per unit area resulting from this description can be expressed in a succinct form in which both contributions to the energy resemble the Helfrich bending term,
plus terms of higher order in K. Here, κ is the bilayer bending modulus κ+ + κ-, and κ nl is a nonlocal "bending modulus" arising through stretching and compression of the leaflet areas, given by κ nl = z 2 + KA+ + z 2 -KA-. K is the average of the midsurface curvature over the whole membrane area. For surfaces of constant mean curvature, which will be our primary interest in this work, this distinction can be discarded. The quantities K 0b and K 0s define the optimal curvatures that minimize the parts of the free energy arising due to bending and stretching, respectively. They can be shown to be 7,18
and
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
In the above equations, K 0± and A 0± are the monolayer spontaneous curvatures and rest areas, respectively. It then follows that a general asymmetric membrane's equilibrium curvature preference can be found by minimizing Eq. ( 3),
This expression makes clear that the preferred bilayer curvature arises as a balance between the curvatures that optimize the two contributions to the free energy, weighted by their respective moduli.
The equilibrium curvature K ⋆ 0 just derived assumes that the bilayer is able to relax its curvature and area. If one simulates an asymmetric lipid membrane (differing in both lipid species and number in each leaflet) using the typical MD simulation setup of fully periodic BC, however, the membrane will almost invariably remain flat. 19,20 This is despite the fact that there will in general be a sizable differential stress present in the system under such conditions, even when the membrane is allowed to relax its area under conditions of zero net tension. 7,8 One might expect, intuitively, that the bilayer would tend to relieve the relative area strain by bending, similar to the way in which a bimetallic strip bends upon heating. This is not the case due to the restraining influence of the PBC.
At a hand-waving level, this can be understood rather straightforwardly: Any bending the membrane would undergo to relieve area strain differences between the two leaflets has to be undone somewhere else in the simulation box in order for the membrane to remain continuous across PBC. We can formalize this idea by directly calculating the area difference between the two monolayers' reference surfaces. Consider a square membrane patch prepared in a flat configuration under PBC with side lengths L.
Since the reference surfaces are parallel-displaced away from the bilayer midsurface, their area elements dA± can be related to the midsurface element dA via the parallel surface relation dA± = dA(1 ± z±K + z 2 ± K G ), 15 as illustrated in Fig. 1. We thus find
where we have again used the Gauss-Bonnet theorem to discard the integral of K G . Let us now consider membranes that can be parametrized in Monge gauge, meaning, by a height function h(x, y)
If the membrane is nearly flat, gradients are small, |∇h| j 1, and in this limit area element and curvature simplify to dA ≈ dx dy and K ≈ ∇ 2 h. In this case, it immediately follows that
At the second equality, we use the divergence theorem to transform the integral over the base plane into an integral along the square boundary with outward-pointing normal l. Under PBC, opposite sides of the simulation cell contribute equally (same shape) but with an opposite sign (direction of l flips), such that the entire expression integrates to zero. We thus see why membranes subject to PBC tend not to relax differential area strain between leaflets by bending into the third dimension: Doing so would not actually relax anything, but rather introduce curvature strain energy with no compensatory benefit. As such, membranes simulated subject to PBC eventually resort to alternative mechanisms to relieve (sufficiently large) differential stress, such as ejecting lipids from the compressed leaflet in the form of micellar buds, as seen in recent coarse-grained simulations. 20
We seek a protocol to simulate asymmetric lipid membranes such that the membrane is able to relax both its area and curvature simultaneously. As we have just elaborated, this is not possible for membranes subject to PBC in all directions. If one simply cuts the membrane along one of the lateral directions to break the periodicity, yielding a membrane strip of finite width and infinite (periodic) length in one direction, then the desired relaxation can occur. However, this solution is ultimately self-defeating, as membrane edges constitute defects along which lipid flip-flop is strongly accelerated. 21 Such a membrane would rapidly equilibrate lipid chemical potentials between the two leaflets, yielding a symmetric membrane with K ⋆ 0 = 0. However, there is a very easy fix to this problem: tape up the edge.
To rescue this free-curvature protocol, we introduce a new element to the simulation setup: adhesive "sticky tape" adhering to the membrane edges that blocks lipid flip-flop. Figure 2 shows the basic design implemented alongside our CG lipid model (described further below). The idea is straightforward, and we will describe it here in general terms applicable to MD lipid models of essentially any resolution. Implementation details that are specific to our CG lipid model are discussed in the supplementary material.
The sticky tapes have a height roughly equal to the hydrophobic thickness of the membrane, as they are designed to adhere to the lipid tails. Only one side of each tape structure is "sticky" (that is, has an attractive interaction with the lipid tails), because we do not wish lipids to "flow around" the tape. Due to some idiosyncrasy of our CG
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp model, the sticky side is subdivided into two regions corresponding to the upper and lower leaflets: The upper half only interacts favorably with lipids that have been designated as belonging to the top leaflet, and similarly the lower half only adheres to lower-leaflet lipids (for more on this, see the discussion of our model below). We do not believe this to be an essential feature of our method, though, especially not when replicated in an atomistic simulation.
The reverse side of the sticky tapes interact with other simulation particles through purely repulsive potentials, as this prevents the tape from being engulfed by the membrane. The tapes have length approximately equal to the length of the simulation box in the (now single) direction of membrane periodicity in order to completely cover the membrane open edges. The tapes should have a fairly high rigidity, such that the membrane edge remains straight and parallel to the direction of periodicity. One could even design the tape structures to be self-connected along the periodic direction, and additionally under tension in order to maintain their straight geometry, but we do not take this approach here.
To evaluate the sticky tape protocol, we employ MD simulations of the CG Cooke lipid model. 19,22,23 This model belongs to the class of very highly coarse-grained representations (just a few beads per lipid), which typically come with an artificially high flip-flop rate that impedes maintaining compositional or stress asymmetry. We therefore employ its recently developed "flip-fixed" variant that circumvents this limitation. What follows here is a very brief summary; for a detailed description, see Ref. 19.
The flip-fixed Cooke lipid model is an implicit-solvent model representing individual generic lipids as four CG beads in a row. One bead represents the head group (blue in Fig. 3), with the other three defining the tail region (yellow in Fig. 3). As there is no water in this model, the fluid phase is stabilized via a cohesive attraction between the lipid tails. Head beads interact with other beads through purely repulsive potentials. In the flip-fixed version of this model, lipids are additionally labeled according to the leaflet in which they are initially placed in order to penalize flip-flop, which is accomplished by disabling the attractive interaction between the two middle beads of lipid tails belonging to opposite leaflets. Observe that this leaflet-designation does not require lipids to be chemically distinct. Furthermore, we can exploit the existence of this label in the construction of the sticky tape: by having the adhesive side facing the upper leaflet only be adhesive to upper leaflet lipids and vice versa. More finely resolved models with intrinsically low flip-flop rates would not need to resort to such a leaflet-labeling trick, and so we would not have to make it part of the sticky tape construction either.
To illustrate the workings of our sticky tape setup with some nontrivial leaflet-based elastic asymmetry, we additionally introduce in this work a tapering angle α that determines the overall lipid shape, generating a set of lipids with differing intrinsic curvature preference (see Fig. 3). Changing the shape of the lipids alters more than just spontaneous curvature, and as such each lipid shape employed in this work (each corresponding to a particular value of α) was run through a series of benchmarks to determine relevant elastic parameters. These include area per lipid (APL) a ℓ , monolayer bending modulus κm, and monolayer stretching modulus KAm. The values of these parameters, as well as the procedures by which they were measured in simulation, are provided in the supplementary material. In order to ensure that all membranes remained in the fluid phase throughout our simulations, we simulated at a slightly higher temperature than originally proposed in Ref. 19, as noted below.
To validate our design principles, we carried out several simulation series in which the number and types of lipids in either leaflet are changed. This collection of simulations consists of three sequences, each comprising five simulations in which the lipid type remains fixed in both leaflets but their number is changed.
In the first sequence, both leaflets contain identical lipids with taper angle α = 0 (the "default" Cooke lipids). The initial simulation is fully symmetric, with 512 lipids in each leaflet (1024 in total). In each subsequent simulation, the number of lipids in the + leaflet (N+) is increased by >2% of the initial 512 (rounded to the nearest whole lipid), while the number in theleaflet (N-) is decreased by the same amount. This setup demonstrates the development of curvature preference purely due to leaflet area imbalance between the two leaflets. The second simulation sequence consists of bilayers in which the + leaflet contains negatively tapered Cooke lipids (α = -1 ○ ) while theleaflet is populated with positively tapered lipids (α = 0.5 ○ ). Each subsequent simulation is modified in the same manner as the first series; N+ is incremented and N-is decremented equally. The third sequence in this set has α+ = 0 and α-= -1.5 ○ . Unlike the previous two sequences, here we decrease N+ while increasing N-, though the magnitude of the change is the same as in the previous cases. This collection of simulations serves to verify that the sticky tapes are able to maintain stable asymmetric membranes with open edges across a variety of curvatures and states of differential stress.
Going beyond simply testing whether the new protocol functions on a basic level, we carried out a second collection of simulations inspired by a simplifying special case of Eq. ( 6). If the + andleaflets contain identical lipids (as in the case of the first simulation sequence described above), then K 0b = 0 and 6) then takes the form (see supplementary material)
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
where we have introduced the number asymmetry parameter δn = (N+ -N-)/(N+ + N-), which measures the bilayer's fractional deviation away from number-symmetry.
The moduli KA and κ are measurable in simulation through a variety of protocols, K ⋆ 0 can be determined from the resulting geometry by fitting the projected tail-bead positions to a circular segment, and δn is user-controlled. Thus, the only unknown in Eq. ( 9) is zn, the distance from the bilayer midplane to the neutral surface of each monolayer. If we run a sequence of simulations containing only one lipid type, over a range of number asymmetries δn and for taper angles α ∈ {-2 ○ , -1.5 ○ , . . . , 0.5 ○ }, we can then determine the neutral surface position zn and find how it depends on lipid shape.
All simulations presented in this work were carried out using version 4.1 of the ESPResSo MD package. 24 All sticky tape simulations as elaborated in this section were run under constant NVT conditions using a Langevin thermostat with k B T = 1.5ε and friction constant = 1τ -1 . The simulation cell dimensions were set to (Lx, Ly, Lz) = (60σ, 16σ, 60σ) and the integration time step was set to δt = 5 ⋅ 10 -3 τ. The membrane is initially placed in a flat configuration with its normal vector along the ẑ direction, spanning from x = 0 to x = 39σ, and continuous in the periodic y-direction. The sticky tape structures then placed on each membrane edge, with the "inward-facing" layers of adhesive CG beads placed at x = -1σ and x = 40σ. Figure 2(c) shows a snapshot from the equilibrium portion of one of these simulations, viewed looking in the +ŷ direction. Each system was run for a total of 3 ⋅ 10 5 τ. Analyses were carried out on the latter 2 ⋅ 10 5 τ of each simulation, discarding the initial equilibration period.
0 as a function of asymmetry δn within the bulk membrane region. Measurements of the average curvature K were taken at the bilayer midsurface during the equilibrium portion of each simulation. The standard error of the mean is smaller than the plotted points in all cases. Dashed lines are linear fits to the simulation data. × symbols indicate state points at which the upper and lower monolayer rest areas are equal. ☀ symbols indicate state points that yield on-average flat bilayers.
We found the sticky tape protocol to successfully maintain membrane asymmetry and differential stress in open-edge membranes for all the cases simulated. The membrane in each simulation is able to dynamically adjust both its area and curvature subject to free boundary conditions, allowing the elastic free energy to assume its minimum (although still subject to periodicity in one direction). Figure 4 shows the resulting average bilayer curvature K ⋆ 0 , measured at the midsurface, as a function of lipid number asymmetry for each simulation in our three series. Notably, the variation of curvature is found to be linear in δn over the entire range of asymmetries investigated, even nearing critical asymmetry values beyond which spontaneous breakdown of bilayer asymmetry is expected in our CG model. 19 This is true both in the case of lipidomic symmetry (blue circles in Fig. 4) as well as lipidomic asymmetry (orange triangles, green diamonds).
Some particular features of interest are highlighted in Fig. 4 by the × and ☀ symbols. A common choice for the construction of asymmetric membranes in MD simulations subject to PBC is to assemble the two monolayers such that the individual monolayer rest areas are equal. [25][26][27] One way to determine a δn value that ostensibly achieves this goal is to run symmetric bilayer simulations of the two respective individual leaflet types (or even compositions); this is sometimes called the "area per lipid (APL) protocol." For our three simulation sequences, the δn values resulting from such an approach are indicated in Fig. 4 by the × symbol. Notice that for the compositionally asymmetric systems (green and orange data), this results
Comparison of results for the location of the neutral surface and pivotal surface as a fraction of the monolayer thickness (determined by the mean position of the lipid head bead). Dark blue triangles are from the sticky tape curvature measurements, red points are from the lateral stretching modulus profile (described in more detail in the supplementary material), and the green crosses are measurements of the pivotal plane based on the method of Ref. 28. Error bars represent the error of the mean. The listed p-values are the probability of a zn difference zn between the two methods at least as big as the observed one occurring by chance, under the null hypothesis of an identical underlying distribution. Inset: zn as a function of taper angle; the blue dashed line shows the best fit to a constant systematic offset.
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp in membranes with sizable nonzero curvature. This means that such membranes simulated in flat configurations subject to PBC would be elastically strained, resulting in a (perhaps unexpected) residual differential stress. 7,13 If one prefers to simulate bilayers that voluntarily assume an on-average flat conformation, then the state points indicated by ☀ symbols in Fig. 4 give the corresponding δn values to use to set up such a simulation.
As explained in Sec. II, our second sequence of simulations comprises six sets of five simulations, each being analogous to the blue data in Fig. 4, but for Cooke lipids of varying taper angle α. For each set of fixed α simulations, fitting to the slope of K ⋆ 0 (δn) as given by Eq. ( 9) yields an inferred value for the neutral surface location zn for the given lipid type. The resulting zn(α) values are plotted in Fig. 5 as dark blue triangles. We find that as the lipid taper angle α is increased, zn also increases monotonically. That is, monolayers composed of lipids with more positive intrinsic curvature preference tend to have their neutral surfaces located farther away from the bilayer midplane.
Ideally, the only role of the sticky tape is to adhere to the membrane edges and prevent lipid flip-flop. However, the introduction of an attractive surface at the membrane edge somewhat predictably leads to changes in membrane properties near the edge, which decay away as one travels inward from the edge toward the bulk membrane phase. In order to quantify the range of influence of sticky tape boundary effects, we can examine the local area per lipid, a ℓ , and the hexatic order parameter |ψ 6 | as a function of distance from the central axis of the membrane. These are shown in Fig. 6, calculated from a sticky-taped simulation of a symmetric standard Cooke lipid membrane system (blue point in Fig. 4). This allows us to make direct comparison with the values for these quantities obtained from standard PBC simulations at zero tension, also shown in Fig. 6. A more detailed comparison of full |ψ 6 | order parameter distributions from both sticky-taped and PBC simulations is presented in the supplementary material.
In the center of the membrane, there is excellent agreement between the order parameters in the two systems. As one gets closer to the membrane edge, and therefore the sticky tape, membrane order increases, and area per lipid correspondingly decreases. The approximate cutoff between the bulk and edge phases is shown by the dotted line in the figure. All analyses pertaining to membrane properties, such as curvature and neutral surface, are performed using only information from the unperturbed bulk. As the number of lipids in the bulk portion of each monolayer is determined by the equilibration of lipid chemical potentials between the bulk and edge phases, the observed number asymmetry of the bulk phase can slightly differ from δn, the globally imposed asymmetry. We refer to this bulk asymmetry as δn b , and it is this asymmetry that is plotted on the horizontal axis of Fig. 4.
Relatedly, our sticky-taped membranes are somewhat reminiscent of scaffolded lipid nanodisks, 29 with the distinguishing feature of our protocol being the existence of a single direction of infinite FIG. 6. Area per lipid a ℓ (s) (top) and hexatic order parameter |ψ 6 | (bottom) as a function of distance s from the membrane arc midpoint (purple curves). The thickness of the curves corresponds to the error of the mean. The horizontal black line in each plot is the mean value determined from a flat PBC simulation. The vertical dashed line gives the approximate cutoff between the bulk and edge regions, showing that (in our setup) the influence of the sticky tape reaches about 5σ into the bulk. periodicity, contrasted with nanodisks' inherent finiteness. The fact that our sticky-taped membranes' bulk properties are roughly in line with their untaped counterparts may then be somewhat surprising, given that nanodisks often have bulk properties that differ from their native counterparts. 30,31 This could be entirely geometric, though: The open edges of our sticky-taped membranes are on average straight lines parallel to the direction of membrane periodicity. There is then no surface Laplace pressure arising due to the line tension γ of a curved boundary, Σ L = γ/R. The sticky tape geometry also allows the bilayer to freely adjust its area without bending or bulging due to confinement by the scaffold. Moreover, since the edges are straight, they have no geodesic curvature, no matter how much the membrane deforms in order to relax a bending torque. This ensures that no Gaussian curvature contribution enters via the edge-something we could not guarantee for the circular edge of a nanodisk. Interestingly, recent simulations of lipid bicelle systems, 32 which are in many ways similar to lipid nanodisks, do not seem to exhibit such noticeable deviations in bulk properties, though this could be due to particular simulation details, which we revisit later.
In Fig. 5, we present the result of our neutral surface measurements based on Eq. ( 9). The natural question to ask is how these results compare to other methods for determining the location of the monolayer neutral surface in simulation. The authors of the work of Campelo et al. 17 ARTICLE pubs.aip.org/aip/jcp of the monolayer lateral stress profile σ(z) with area strain and is defined as
This function can be approximated by calculating σ(z) for several small, flat, PBC simulations at a series of increasing area strains (as explained in the supplementary material). Figure 7 shows λ(z) found for a single-component Cooke lipid membrane with α = 0.5 ○ . The curvature-area cross-coupling modulus turns out to be given by the first moment of λ(z) with respect to the reference surface height z 0 . 17 By definition, this modulus vanishes for a reference surface at zn, implying
Here, z is measured from the bilayer midplane and the integral upper bound h is the total height of the monolayer, which in practice can be taken arbitrarily large since λ(z) → 0 rapidly once outside the membrane (see again Fig. 7). The results of this calculation for our Cooke lipid systems are also shown in Fig. 5 (red points). By eye, our new method for determining zn seems to be in fairly good agreement with the λ(z) profile method. Indeed, for each individual pair of measurements, we can calculate a two-tailed p-value under the assumption of Gaussian errors (reported in Fig. 5), which would suggest compatibility of the two methods. However, taken together, we see that all of the zn values calculated via λ(z) are less than those calculated from our analysis of K ⋆ 0 (δn), which we would only expect to occur >3% of the time by chance. The inset of Fig. 5 plots the differences zn between the two methods of calculation for each lipid type, along with a fit to a constant offset (blue dashed line), found to be 0.08 ± 0.04σ, which is about 2% of a lipid height, and two standard deviations away from zero. While very close, there does appear to be a small systematic difference between the two methods.
The method we have presented here does not require calculation of multiple lateral stress profiles and their connection to the overall moduli. However, it does require running multiple sticky tape simulations from which the mean curvature is measured. Regardless of which method is more efficient, it is reassuring to see that a measurement that relies on direct observation of the largescale geometric response of the membrane is in good agreement with micro-elastic considerations.
It is informative to compare our result for the neutral surface location zn with the location of another common monolayer reference surface, the pivotal plane, zp. The defining property of this surface is that it is the location of zero area strain upon pure membrane bending. We measured zp for all of the lipid shapes employed in this work using the method presented by Wang and Deserno, 28 which relies on counting the lipid imbalance between the two leaflets of a curved membrane buckle. The results of this analysis are presented alongside the neutral surface results in Fig. 5. The measured values for zp and zn are indistinguishable within error-a slightly surprising result, given that there is no reason to expect these two locations to coincide. Indeed, the values of κm and K 0m are generally dependent upon the choice of reference surface. This might reflect the inherent simplicity of our highly coarse-grained lipid model, as these two surfaces are generally not found to coincide experimentally. [33][34][35][36]
As alluded to in Sec. III A, there have been several protocols presented in the literature for how to assemble and carry out MD simulations of general asymmetric lipid membranes. These range from simply matching the total rest areas of the lipids on each side (as previously mentioned), [25][26][27] to positing that the two leaflets should be simultaneously tensionless, 37 to much more sophisticated schemes involving equilibrating chemical potentials of specific lipids between the two monolayers through the use of nontrivial boundary conditions. 13 Going beyond infinite periodic protocols, the previously mentioned bicelle setup presented in the work of Pöhnl et al. 32 has similar aims to our sticky tape method. Their protocol is to simulate specially restrained lipid bicelles, which are essentially nanodisks whose edges are stabilized by detergent molecules or short-tailed lipid species. 38 For properly tuned mixtures, the high-curvaturepreferring short-chained species localize at the disk rim, with the bilayer-forming lipids creating the core domain. Such systems have previously been re-created in simulation in order to, e.g., investigate peptide-induced membrane curvature. 39 The setup presented in the work of Pöhnl et al. 32 includes artificial restraining potentials that maintain selectively chosen lipids either within or outside a given cylindrical region to maintain separation between the bulk and rim phases. Interestingly, unlike other lipid disk protocols, their simulations do not exhibit noticeable deviations in lipid density, suggesting a beneficial influence of the external rim potential. While successful, it remains somewhat unclear how the membrane is affected by the fixed restraining potentials, which should in principle suppress membrane bending beyond certain thresholds. The presence of a curved interface at the membrane edge also raises concerns about the influence of the often-neglected boundary term in the Helfrich
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp energy (say, a κ-contribution via the Gauss-Bonnet theorem, coming from the boundary's geodesic curvature), as well as some form of radial compression via the Young-Laplace pressure, as discussed above.
All these protocols strive to realize certain "elastic ensembles" in which a particular set of extensive (like area) or intensive (like stress) thermodynamic variables are set. What we do not know, of course, is what the right ensemble would be in the first place. It may well depend on the situation whether a bilayer with vanishing differential stress is physically relevant, a bilayer with vanishing curvature torque, or yet some other condition. Ideally, one would know this from experiment, but this can be tricky, for instance because currently no method exists to measure the differential stress. One might have to indirectly infer the most appropriate ensemble, or simply make an executive decision in this matter. At any rate, in the present paper, we take no view on the "correct" boundary condition and merely wish to provide tools that help to enact or identify certain choices, which the user needs to justify by independent means.
Our sticky tape protocol allows for membranes to assume their preferred mean curvature conformation while simultaneously relaxing the overall bilayer area. If one prefers to simulate a flat system subject to full PBC, the sticky tape method is still potentially useful because it allows one to determine the composition and number asymmetry that renders the membrane voluntarily flat (as shown in Fig. 4). This can then be transferred to a fully periodic box for production simulation. Observe that such systems are generally under differential stress, but its origin is physically clear: We chose an ensemble in which the equilibrium shape is flat and hence satisfies the different mechanical equilibrium condition of zero torque-unlike setups such as the APL protocol that introduce hidden differential stress via the forced "unbending" of a curvature-preferring non-torque-free bilayer by the PBC.
We also wish to point out that the freely varying curvature conditions of sticky tape simulations are more than just a useful trick for finding a flat state. The sticky tape protocol also opens the possibility of simulating membrane-interacting proteins without the constraints imposed by full PBC. Such simulations could for instance provide insight into curvature sensing and/or induction.
We have presented the novel "sticky tape" protocol for the simulation of asymmetric lipid membranes under simultaneous freearea and free-curvature conditions. Its stability and robustness have been demonstrated in the context of the ultra-coarse-grained Cooke model, in which we find that asymmetric membranes maintain their asymmetry and relax to their preferred areas and curvatures, even when subject to sizable differential stress. We have also shown how the newly unlocked simulation ensemble can shed light on an essential elastic parameter of a lipid monolayer: the location of its neutral surface. So far, the method has only been implemented in the context of the Cooke lipid force-field. The design principles of the sticky tape protocol are, however, not specific to coarse-grained models and can be generalized to higher resolution systems.
While the portions of the membrane closest to the edges exhibit deviations from native behavior, the properties of the bulk phase are unperturbed by the presence of sticky tapes. It should also be emphasized that in this initial exploration, we made no effort to systematically tune the adhesive interaction potentials in a way that could minimize these edge deviations. As the structure of fluids is strongly influenced by the repulsive part of the pair potential, 40 softening the core part of the sticky tape-lipid attraction, or a weakening of the adhesion strength (depth of the potential), could potentially lessen the finite-size effects seen here.
The supplementary material for this article includes the following: a detailed description of the tapered Cooke lipid model; details of the implementation of the sticky tapes with this model; the derivation of Eq. ( 9); details of the computation of lateral stretching modulus profiles; plots of Cooke lipid hexatic order parameter distributions. An additional archive supp_code.zip is provided with Python scripts for running sticky tape membrane simulations using the ESPResSo MD package.
Published under an exclusive license by AIP Publishing
160, 064111-3Published under an exclusive license by AIP Publishing