Determining Structural and Mechanical Properties from Molecular Dynamics Simulations of Lipid Vesicles.

We have developed an algorithm to determine membrane structure, area per lipid, and bending rigidity from molecular dynamics simulations of lipid vesicles. Current methods to extract structure from vesicle simulations define densities relative to the global center of mass of the vesicle. This approach ignores the long-wavelength fluctuations (undulations) that develop across the sphere and broaden the underlying structure. Our method establishes a local reference frame by defining a radially undulating reference surface (URS) and thereby removes the broadening effect of the undulations. Using an arc-length low-pass filter, we render the URS by defining the bilayer midplane on an equi-angular θ, ϕ-grid (colatitude, longitude). This surface is then expanded onto a truncated series of spherical harmonics. The spherical harmonic coefficients characterize the long-wavelength fluctuations that define both the local reference frame-used to determine the bilayer's structure-and the area per lipid (AL) along the undulating surface. Additionally, the resulting power spectrum of spherical harmonic coefficients can be fit to a Helfrich continuum model for membrane bending in spherical geometry to extract bending rigidity (kc). kc values determined for both DMPC and DMPC + cholesterol (30 mol %) vesicles are consistent with values from corresponding flat-patch systems determined using an independent, previously published spectral method. These new tools to accurately extract structure, AL, and kc should prove invaluable in evaluating the construction and equilibration of lipid vesicle simulations.

Molecular dynamics (MD) simulations are commonly used to study the structure (number density ρ(z) and area per lipid AL) and bending rigidity (kc) of lipid bilayer systems.[1−6] Recently, MD simulations of lipid vesicles of biologically significant size (radius >20 nm) have become possible due to more reliable coarse-grained force fields,[7−13] but analysis tools for extracting structure and rigidity are lacking. In simulations of flat-patch bilayers (<200-lipids), the membrane structure normal to the bilayer (z-dimension) can be determined by binning the relative z-position for each bead in the system relative to the global bilayer center of mass (COM).[3] A similar approach has been employed for simulations of vesicles, where the global COM of the vesicle defines the radial reference frame for determining the radial membrane structure profile.[8,9,11,14]

However, in simulations of large bilayers—both flat-patch and vesicles—long-wavelength fluctuations convolve a smoothing function with the intrinsic structure profile. This results in a broadened structure when the calculation relies upon the global COM. We previously developed a method to remove this effect in the flat-patch geometry that defines a local reference frame using an undulating reference surface (URS).[3] Our approach—surface referencing with undulation correction—characterizes the long-wavelength fluctuations and removes them from the bilayer structure calculation. This paradigm relies on the notion that undulations introduce fluctuations in the lipid bilayer’s local normal, but do not alter the underlying membrane structure.[15] Because the undulation correction isolates the local normal vector, the resulting number density profile is system-size independent.[3]

A bilayer’s number density profile, ρ(z), is a one-dimensional measure of its structure. The area per lipid (AL) is directly coupled to ρ(z), and is of fundamental importance in describing changes in membrane structure.[5,16−19]AL is commonly determined as a projected area, defined as AL = NLipids/2LL, where NLipids is the number of lipids and L, L are lateral periodic box dimensions.[17−19] A more appropriate AL metric determines the area along the URS for the system, as we showed that the correction accounts for a systematic ∼1% change in the calculated AL.[3]

One important question that arose in our previous study, and which is further complicated in the case of vesicles, is what should be the gold-standard against which extracted structures should be compared? This issue is confounded by the fact that experimental structure determination for bilayers (e.g., from X-ray scattering) often use simulations or modeling as a guide.[3,18−21] For our previous work, we used a flat-patch system small enough to not develop any undulations (Nlipids < 200) as the comparison point for structure profiles extracted from much larger systems.[3]

With vesicle simulations, the appropriate gold-standard is not as clear. Lipid asymmetry may exist in a radius-dependent manner. If a vesicle is large enough, it could be compared to the flat-patch. However, at smaller length scales, local curvature alter the distributions. Furthermore, the methods for building vesicle starting configurations are still evolving. It is difficult to know exactly how many waters to include in the interior of the vesicle as well as the appropriate ratio of lipids on the inner and outer monolayer. It is also unclear how long simulations need to be in order to ensure the appropriate lipid distribution is met in the two monolayers (e.g., lipid flip-flop). It has not been our goal to address these important complexities, for which an extensive effort is needed but for which the appropriate measures for evaluation is lacking. Indeed, our new algorithm should provide a more quantitatively reliable framework for rationally building and testing physical characteristics (e.g., water density) and convergence.

Regarding bending rigidity, methods to extract kc from periodic, flat-patch lipid bilayers rely on Helfrich-like continuum models. In this well-known model, the continuum behavior of the bilayer manifests as a wave-vector to the inverse fourth power for long wavelengths (undulations).[1,2,22] Similar theoretical treatments have been developed for spherical vesicle geometries[23−26] and are commonly used in interpreting fluorescent and neutron spin–echo experiments to extract bending rigidity from lipid vesicles.[27−29] Several algorithms have been developed for extracting kc from flat-patch bilayer simulations (e.g., spectral fluctuations analysis, simulated buckling, and lipid tilt modulus).[2,30−32] However, implementation of any of these approaches in spherical geometry has not yet been developed within the MD simulation literature, making it critical that we develop appropriate tools to handle the more complex geometry.

In this study, we develop a novel algorithm to extract the undulation-corrected number density profile (ρ(r)), area per lipid (AL) and bending rigidity (kc) from vesicle simulations by eliminating the smoothing effect inherent in undulating bilayer systems. Translating current algorithms for extracting these bilayer properties from lipid vesicles is far from trivial because of the loss of periodicity and closed membrane envelope. Our algorithm is based upon a spectral method that characterizes the radial undulating membrane fluctuations using spherical harmonics analysis (SPHA), the Fourier analysis corollary for spherical geometries. We thereby provide a seamless link to our previous work in flat-patch systems.[2]

METHODS and ALGORITHM DEVELOPMENT

Defining the Vesicle Surface, r(θ, ϕ)

In order to characterize the radial undulations from a simulated vesicle, we start by recasting the positions of the bilayer beads into an equi-angular discrete surface representation using a θ, ϕ-grid (colatitude and longitude respectively), where θ ∈ [0, π] and ϕ ∈ [0, 2π], with θ = 0 at the northern pole. The angular resolution is defined as where λS is the arc-length where the bilayer fluctuations transition between a continuum mode to molecular (λS ≈ 4 nm), r0 is the average vesicle radius, and ns is the number of sampling points per λS with ns ≥ 2 to satisfy Nyquist sampling theorem.

To render the surface, we define the origin at the COM of the vesicle (lipids + membrane inclusions) and then transform the system into spherical coordinates (x, y, z ⇒ r,θ, ϕ, where r,θ, ϕ corresponds to radius, colatitude, and longitude for each bead i, respectively). The direct use of an equi-angular -grid to bin bead positions introduces discontinuities at the poles. We mitigate these discontinuities by implementing an arc-length low-pass filter with filter cutoff qarc = 2π/λS. Our previous work with undulation analysis on flat patch systems identified a cutoff wavenumber, q0 = 1.5 nm–1, at which the long-wavelength undulations transition into molecular structure fluctuations.[2,3] We use this q0 as a guide in defining an appropriate qarc for the vesicle systems.

Figure 1A shows a snapshot from an equilibrated, 34 nm DMPC vesicle, and Figure 1B presents a cutaway view, emphasizing the vesicle’s diameter relative to its thickness. We parse all lipids into the inner and outer monolayer using the average vector between the acyl-chain terminal beads and the phosphate bead, where the radial component defines which monolayer the lipid belongs to (i.e., r < 0, inner monolayer; r > 0, outer monolayer). Using an arc-length filter with qarc = 2.5 nm–1, we define the local θ, ϕ-membrane surface for both inner and outer monolayers, rin(θ, ϕ) and rout(θ, ϕ), respectively. Figure 1C schematizes the first stages of our method, with the inner (red) and outer (blue) monolayer beads shown for numerous θ, ϕ-positons. On the left hemisphere of Figure 1C, regions of rin(θ, ϕ) and rout(θ, ϕ) are presented for illustration. The average of the two surfaces define the undulating radial surface as

Figure 1

A. Snapshot of a 34 nm DMPC-lipid vesicle. B. Cutaway of the vesicle. C. An arc-length filter is used to define a surface for both inner and outer monolayers (blue and red selection beads respectively). D. rund(θ, ϕ) is determined as the average of the inner and outer monolayers (color-map indicates fluctuations about the average radius, r0, (units in nm).

Figure 1D shows rund(θ,ϕ) calculated from one frame of the DMPC trajectory. The color-map displays the deviations about the average vesicle radius, r0′ (units in nm), which is typically smaller than the ideal sphere radius r0. From rund(θ, ϕ), we then define the normalized radial fluctuations as

Spherical Harmonics Analysis (SPHA)

Next, spectral decomposition of rund(θ, ϕ) is accomplished using spherical harmonics analysis (SPHA). Spherical harmonics are standing waves on a sphere. Given that rund(θ, ϕ) is a discrete surface defined on a spherical manifold, we can represent it as a linear combination of spherical harmonics with degree (l) and order (m), corresponding to the number of waves in the θ, ϕ-dimensions, respectively. Figure 2 illustrates three standing wave patterns as calculated from our simulations.

Figure 2

SPHA decomposes rund(θ, ϕ) into a series of standing waves on a sphere with degree l and order m.

Helfrich’s formulism for SPHA expands the normalized radial fluctuations in spherical harmonics aswhere a are the spherical harmonic coefficients of degree l and order m (l ∈ 0,1···,lmax and m ∈ −l,...,l) with Lmax defined by the number of colatitude grid points.[23]Y are the spherical harmonic basis functionswhere P̅(θ) are the fully normalized associated Legendre polynomials with the normalization factor, N, defined as

Because f(θ, ϕ) is a real-valued function a = a*, the complex conjugate of a. We thereby redefine a such thatallowing us to writewhere now m ∈ 0,...,l.

We generate the matrix with dimensions (2N * Lmax) × N2 for a given θ, ϕ-distribution,[33] where N is the number of colatitude parallels, such that for each l ∈ 0,...,Lmax, the corresponding 2m + 1 rows are defined asfrom which the transformation matrix follows as

From the matrix we can write the spherical harmonic forward transform aswhere a1 is a recasting of the spherical harmonic coefficients a with dimension N2 × 1 corresponding to the row construction of , and f1 is a matrix of positions with dimension (2N * Lmax) × 1.

Whereas the forward transform is exact, the inverse transform (i.e., SPHA) is overdetermined. Using the FACTORIZE[34] package in Matlab, we compute the pseudoinverse of using QR-factorization and apply a least-squares approximation to determine a1 as

Implementation of the SPHA algorithm results in error propagation from numerous sources (i.e., truncation error, round-off error, and least-squares approximation). We evaluated the magnitude of these effects by calculating the root-mean-squared-difference (RMSD) between rund(θ, ϕ) and the transformed r̆und(θ, ϕ), where

RMSD calculated for the DMPC system is 3.5 × 10–3 nm and for the DMPC + cholesterol system is 1.8 × 10–3 nm. These errors are 2 orders of magnitude lower than the magnitude of the radial fluctuations, O (0.5 × 10–1 nm), allaying concerns regarding the inherent propagating errors.

Determining Flat-Patch Structure (ρ(z) and AL)

Membrane structure (ρ(z) and AL) for flat-patch membrane systems was determined using our previously published method.[3] Briefly the URS was determined using the phosphate atoms and subsequently filtered with a 1.5 nm–1 cutoff. The filtered surface was then used to reference every atom to the local bilayer midplane. These distances were subsequently binned and volume normalized to extract the number density profile.

Determining Vesicle Structure (ρ(r) and AL)

The standard method for determining the radial membrane number density profile, ρrbin(r), references every bead/atom in the system relative to the COM of the vesicle in spherical coordinates, averaging across all θ, ϕ-angles. We have previously shown that the membrane structure profile is smoothed out in systems where long-wavelength undulations develop.[3] By isolating the long wavelength undulations and referencing every bead/atom to the local undulating reference surface we can mitigate this broadening effect.

Following the same principles for our method developed for flat-patch systems, we apply a low-pass ideal filter to the spherical harmonic coefficients with an order cutoff Lcut = qarcr0′ – 0.5 such thatwhere ã1 are the filtered coefficients. For the DMPC vesicle Lcut = 26, whereas for the slightly larger DMPC + cholesterol vesicle, Lcut = 29. We next apply an the inverse spherical harmonic transform to resolve a filtered radial-undulation surface r̃und(θ, ϕ) as

With the filtered surface we resolve the vesicle’s local membrane structure, ρuc(r), by referencing every bead/atom i relative to it is position on the surface, such thatwhere r̃und(θ, ϕ)nearest corresponds to the closest equal-angular grid point on r̃und(θ, ϕ). The grid definition samples the corresponding filter wavelength at least 2 times (satisfying Nyquist sampling theorem). As a default we employ 4-samples per wavelength (at qarc = 2.5 nm–1 the arc-length resolution is ∼0.4 nm), providing a close approximation to the unique r̃und, (θ, ϕ). The resulting r̃(θ, ϕ) is then binned with a bin-width of dr = 0.1 nm and subsequently normalized by the differential volume shellwith r0—the radius of the ideal sphere with equal volume as r̃und(θ, ϕ) defined asand adjusted for every bin index b by the corresponding distance from the surface.[23]

In addition to number density profiles, we can identify the vesicle’s area per lipid, AL, both as a whole and independently for both monolayers. Instead of defining r0′ with all vesicle beads/atoms, we parse them independently as inner and outer monolayers. Then using eq 18, we obtain two additional ideal sphere radii (three in total), r0, r0, and r0. The corresponding AL for each ideal radius is simplywhere nL, inner and nL, outer refer to the number of lipids in the respective monolayer.

Determining Bending Rigidity (kc)

From a1, we obtain the undulation power spectrum by binning the modulus of the spherical harmonic coefficients across degree l. The resulting profile can be interpreted according to the Helfrich continuum model for undulations on a sphere with vanishing spontaneous curvature aswhere T is temperature and l ∈ 2,...,Lmax.[23]

In the flat-patch spectral method, the undulations power spectrum, ⟨|u(q)|2⟩, determined by direct Fourier transformation of the bilayer selection atoms can be modeled aswhere N is the number of selection atoms, a is the projected area per lipid, and Sρ(q) is the in-plane structure factor.[2] Sρ(q) is subtracted from N⟨|u(q)|2⟩ to provide a broader range of modes within the q–4 regime. kc is then determined by fitting the low-q, long-wavelength modes.

Arc-Length Filter Cutoff Effects

Using a filter to define the URS attenuates the spectral intensity of the undulations. This attenuation can span a broad frequency bandwidth, bleeding through into the desired frequencies of the signal. We explored a range of arc-length filter cutoffs to characterize the filter’s frequency response, with the goal of identifying a cutoff where signal attenuation is limited to frequencies above the crossover wavenumber, q0 = 1.5 nm–1. Figure 3 presents the power spectra for the DMPC vesicle system using four different arc-length filters with wavenumber cutoff, qcut, ranging from 0.5 nm–1 to 2.5 nm–1.

Figure 3

Undulation spectra for a range of arc-length cutoff wavenumbers, (qcut = 0.5 nm–1: blue; 1.5 nm–1: green; 2.0 nm–1: red; and 2.5 nm–1: cyan) for the DMPC vesicle system. The crossover wavenumber q0 = 1.5 nm–1 is denoted by the black dashed line. kc values for each qcut are detailed in the legend.

The frequency response of the arc-length filter is far from ideal, with significant bleed-through extending below the desired cutoff wavenumber. This is most noticeable when comparing the crossover wavenumber q0 = 1.5 nm–1 (black dashed-line)—the wavenumber where continuum undulations transition into molecular fluctuations[2]—to the 1.5 nm–1 filter cutoff (green triangles). There is significant loss of undulation intensity for degrees 12–25, all before the 1.5 nm–1 crossover wavenumber. This loss of intensity skews the kc-fit, resulting in a larger kc. Furthermore, increasing the filter cutoff increases the number of degrees that comprise the linear region below q0, thereby improving the fit to eq 22. Bending rigidities for each filter cutoff are listed in the Figure 3 legend. With qcut = 2.5 nm–1, all 26 degrees (spanning the full range, q0 ≤ 1.5 nm–1) can be used to determine kc.

Molecular Dynamics Simulations

All simulations were run using the MARTINI force field and using the GROMACS v4.5.3 program.[10,35−38] All systems were prepared and run in the isothermal–isobaric (NPT) ensemble at constant pressure and temperature (1 bar and 303 K, respectively) using either a 25 fs time step for pure lipid or 10 fs time step for cholesterol systems. Pressure coupling was isotropic for the vesicle systems and semi-isotropic for flat-patch, with independent xy- and z-barostats, resulting in a tensionless bilayer. Initial equilibration for vesicle systems included 100 000 steps of steepest descent minimization followed by 500 ns of dynamics using the velocity-rescaling thermostat and a Berendson barostat (flat-patch systems underwent 10 000 steps steepest descent minimization and 100 ns velocity-rescaling dynamics).[37,38] All production runs were simulated using the Nosé–Hoover thermostat and the Parrinello–Rahman barostat with a time constant of 2.5 and 250 ps, respectively.[39,40] Pressure coupling was applied isotropically for vesicle simulations and semi-isotropically for flat-patch systems. Vesicle production simulations were 2.5 μs (10 μs scaled simulation time), flat-patch simulations were 5 μs (20 μs scaled time) sampled every 1 ns.

The starting configurations for the DMPC vesicle (11 126-DMPC, 1 123 315-coarse grained water + antifreeze, “high-water”, 1 113 742-cg water + antifreeze “low-water”) and DMPC + cholesterol (12 771-DMPC + 5473-cholesterol, 1 313 698-cg water + antifreeze) were constructed by randomly seeding two opposing monolayers with the appropriate number of lipids (DMPC and/or cholesterol) based on the area per lipid of each species, the surface area of the monolayer shell, and the mole fraction of the mixture. A second DMPC vesicle was simulated using the 3 μs frame as a starting configuration with 10% of the internal water beads removed to explore effects of system construction on structure and bending rigidity. Flat-patch bilayer systems were constructed with 3200-DMPC lipids and 2240-DMPC + 960-cholesterol, randomly seeded in flat monolayers with 70,400 coarse grained water + antifreeze particles.

During the early stages of equilibration, rapid pore formation occurred throughout the vesicles. These pores coalesced and closed quickly—within 70 ns—providing an opportunity for both lipid flip-flop and water exchange across the vesicle to equilibrate the system. For the low-water system, the vesicle collapsed into an ellipsoidal shape within 100 ns, and maintained an ellipsoidal character throughout the 3 μs production run.

Data Analysis

Trajectories were manipulated and processed using both the GROMACS v.4.5.3 simulation package[37,38] and the MDAnalysis python library.[41] Further data analysis and figure rendering was performed using MATLAB (v.R2012a) with use of the FACTORIZE library.[34]

Results

Radial Membrane Structure

Having established the framework to define rund(θ, ϕ) for vesicle simulations, it now becomes possible to extract the underlying membrane structure from the fluctuating vesicle’s trajectory. Figure 4 presents the number density profiles for the DMPC vesicle, ρrbin(r) and ρuc(r), and the DMPC flat-patch, ρuc(z), systems. Vesicle profiles are illustrated as dashed lines, flat-patch profiles as filled distributions, with the component groups color-coded as described in the caption. Negative radial positions correspond to the inner monolayer, positive values correspond to the outer monolayer.

Figure 4

Comparisons of component number densities for DMPC vesicle and flat-patch systems. Vesicle profiles are dashed-lines. Flat-patch profiles are filled-distributions with headgroup (black/gray), carbonyl-glycerol (red/pink), and acyl-chain (blue/cyan). (A), (B) ρrbin for the high-water and low-water DMPC system, respectively. (C), (D) ρuc for high-water and low-water DMPC system. The same flat-patch profile has been included in all panels.

Figure 4A compares the high-water vesicle profile, ρrbin(r), to the flat-patch profile, ρuc(z). The ρuc(z) profile was translated by the average vesicle radius, r0′, to allow direct comparison to the vesicle profiles. As expected, using a global reference frame with a fluctuating bilayer results in a broadening of the component distributions (RMSD between ρrbin(r) and ρuc(z) is 3.08 nm–3). This broadening results in a loss of structural resolution in the vesicle profile, most noticeably in the acyl-chain distributions, where the typically pronounced terminal methyl trough is absent. Figure 4B compares ρrbin(r) from the ellipsoidal, low-water DMPC system to the flat-patch profile. As fluctuations increase, the structural broadening effect intensifies. The low-water ρrbin(r) displays significant distortion (RMSD is 7.12 nm–3) across all lipid-component distributions, highlighting the problem with using a global COM reference frame to define membrane structure.

Figure 4C,D present comparisons of our new method, ρuc(r), to the flat-patch profile for both the high- and low-water DMPC systems, respectively. In both cases, the resulting number density more closely resembles the flat-patch result (RMSD between ρuc(r) and ρuc(z) is 1.66 nm–3 for low-water DMPC and 2.34 nm–3 for high-water). Each component distribution is tightened up relative to ρrbin(r), and key structural features (e.g., the terminal methyl trough) are properly characterized. In the high-water system (Figure 4C), the headgroup distributions are centered at the appropriate position (i.e., the bilayer thickness matches that of the flat-patch). However, the amplitude and width of the headgroup distributions do not fully agree with the flat-patch result: the amplitudes are too high and the widths to narrow. This small discrepancy suggests that the lipids are under tension, either due to water-density or lipid density asymmetry across the bilayer.

Comparison of the low-water profiles (Figure 4B,D) demonstrates the dramatic improvement in calculated membrane structure when undulations are isolated and characterized during structure determination. Even with a dramatic ellipsoidal geometry, the ρuc(r) method extracts the underlying membrane structure. The low-water profile amplitudes agree with the flat-patch profile, more so than the high-water result. However, the membrane appears thicker and the component distributions are slightly broader than the flat-patch result.

Comparing Figure 4C with 4D (ρuc(r) for the spherical and ellipsoidal DMPC vesicle systems), we see a dramatic change in the structural profile of the membrane. As expected, varying the number of waters inside the vesicle can lead to different membrane tensions and corresponding structures. With our new algorithm, this type of comparison can serve as a guide in the initial construction of vesicle simulations, evaluating the membrane structure for imbalances in water density or lipid density across the bilayer.

Analysis of the DMPC + cholesterol system highlights additional complexities that can develop due to the different system geometries. Supplemental Figure 1 presents ρrbin(r) and ρuc(r) for the DMPC + cholesterol system. The undulation correction still sharpens the component distributions, albeit to a lesser extent than the pure DMPC system. Cholesterol’s facile ability to flip-flop during the time-scale of the simulation results in an asymmetric cholesterol distribution across the bilayer in the vesicle geometry. The development of an asymmetric cholesterol distribution is not surprising for the vesicle as the spontaneous curvature difference between inner and outer monolayers will induce an asymmetric partitioning of lipids with non-neutral spontaneous curvature.[42]

Close inspection of all vesicle number density profiles highlight an asymmetry (particularly in the headgroup distributions) where the inner monolayer has a higher density relative to the outer monolayer (Figure 4 and Supplemental Figure 1). We cannot know for certain whether this asymmetry is due to a redistribution of lipid components or an artifact from a poorly equilibrated vesicle.

Vesicle structure does not change across the last 1.5 μs for both the DMPC and DMPC + cholesterol systems. However, the differences in structure between high-water and low-water DMPC vesicle systems highlight the effect of water density imbalances across the bilayer. Furthermore, in all simulations, lipid flip-flop was only observed for cholesterol and never for DMPC. Significantly longer simulation times (>100s of μs) of standard MD would be necessary to allow for both water equilibration and sufficient lipid flip-flop events to fully equilibrate the vesicle.

Our focus here has been on the development of a method to calculate structure from vesicle simulations and not to perfectly refine these particular simulations. Nevertheless, the problem of vesicle construction and equilibration is an important challenge that must be answered to facilitate further study of these more complex systems. Risselada et al. employed an artificial pore allow for both water and lipid exchange across the vesicle prior to the production MD simulation.[9] For simple lipid compositions, this approach is preferable over an iterative refinement of the lipid and water distributions (altering them by hand on the basis of structural profile imbalance) and far more computationally efficient than extending the simulation until sufficient lipid flip-flop is observed.

Area Per Lipid, AL

Figure 5 presents the area per lipid (AL) time-series for the DMPC vesicle systems. We compare values determined along the URS for the full vesicle (ALL, blue) and both inner and outer monolayers (green and red, respectively) with those of the flat-patch AL (cyan). A summary of the average AL is presented in Table 1. In the high-water DMPC system (Figure 5A) there are three distinct differences between the vesicle and flat-patch result: (1) the inner and outer monolayer have different AL, with inner monolayer lipids occupying more lateral area than those in the outer leaflet; (2) all measures of AL are greater than the flat-patch AL; and (3) the fluctuations in AL are significantly decreased in the high-water vesicle system versus the flat-patch. The imbalance between inner and outer monolayers correlates with the structural asymmetry in the headgroup ρuc(r).

Figure 5

(A) AL trajectory for the high-water DMPC vesicle system with total vesicle (blue), inner monolayer (green), outer monolayer (red), and flat-patch system (cyan). (B) AL trajectory for the low-water DMPC system.

Table 1

AL (or AUC) for Flat-Patch and Vesicle Systems

geometry	flat-patch	vesicle
subset	all	all	inner	outer
DMP (high-water)	0.601	0.622	0.623	0.608
DMP (low-water)		0.653	0.652	0.641
DMP + cholestrol	0.644	0.672	0.652	0.684

We presume that this reflects that the inner monolayer is under greater tension than the outer, constraining the inner headgroups in the radial dimension and resulting in a tighter position distribution (recall that the thickness is the same, Figure 3B). There are two potential causes of this tension, either an imbalance in the water density or the lipid density across the vesicle. The area per unit cell, AUC = 2 ASurface/NDMPC, for the DMPC + cholesterol system shows a similar trend as the DMPC vesicle (see Supplemental Figure 2).

Figure 5B illustrates the same AL analysis for the ellipsoidal, low-water DMPC system. It is immediately evident that the water content dramatically affects the vesicle’s structure and dynamics. The ideal sphere AL is increased for all three AL metrics, even though the average vesicle radius is smaller. There is a partial recovery in the AL symmetry across the monolayers (i.e., a reduction in ΔAL = AL,out – AL,in from 0.015 nm2 to 0.011 nm2). The residual discrepancy between AL,out and AL,in likely stems from a lipid density imbalance across monolayers, which is still evident in Figure 4D. Although the structure from the low-water system more closely matches the flat-patch profile, it is the AL of the high-water system that is in better agreement. The number density profile calculation isolates the increased fluctuations in the low-water system, whereas the AL does not.

Bending Rigidity

Figure 6 presents the power spectra and subsequent model fits for DMPC high-water (black) and DMPC + cholesterol (red) in both vesicle (Figure 6A-C) and flat-patch (Figure 6D) geometries. In each case the spectra was determined over the final 1.5 μs of the trajectory. As expected, in both the flat-patch and vesicle simulations cholesterol reduces the magnitude of the undulation intensity, reflecting the sterol’s well-know rigidifying effect.[5,21,43,44]kc values were determined by fitting the linear regime across low degrees corresponding to q0 ≤ 1.5 nm–1 to the appropriate Helfrich continuum model for vesicle or flat-patch geometries (l ≤ 26 for DMPC and l ≤ 29 for the slightly larger DMPC + cholesterol vesicle). As shown in Table 2, we observe a similar cholesterol-induced increase in kc in the vesicle- and flat-patch geometry (the latter of which was obtained via our the previously published method).[2]

Figure 6

Power spectra for both vesicle systems (A) with corresponding kc fit for DMPC (B) and DMPC + cholesterol (C). Undulation power spectrum and kc-fits for flat-patch systems. For all panels, DMPC (black) and DMPC + cholesterol (red).

Table 2

kc Comparison Across System Geometry and Compositiona

kc[J]
system/geometry	flat-patch (3200 lipids)	vesicle (∼12 200 lipids)	flat-patch (30 000 lipids)*	experiment**
DMP (high-water)	2.1 × 10–19	1.3 × 10–19	1.5 × 10–19	5.8 × 10–20
DMP (low-water)		6.1 × 10–20
DMP + cholesterol	3.8 × 10–19	2.4 × 10–19		2.7 × 10–19
ratio	1.8	1.9 (hw); 3.8 (hw)		4.7

The * denotes values from Brandt et al.;[2] experimental values (**) from Pan et al.[45].

The high-water (hw) DMPC vesicle kc is in better agreement with our previously published kc from a large flat patch, where a 30 000 lipid DMPC coarse-grained flat-patch system had kc = 1.5 × 10–19 J.[2] To ensure the use of an equi-angular grid did not introduce any artifacts into the vesicle analysis, we repeated the SPHA on the high-water DMPC vesicle in a 90° rotated reference frame. The results were identical to the nonrotated system (see Supplemental Figure 3).

The increased kc from our smaller 3200 lipid DMPC flat-patch highlights potential system-size effects with spectral methods to extract kc. The low-water (lw) system shows a reduced kc, relative to the high-water kc, agreeing with the increased AL fluctuations (see Figure 5B, Table 2, and Supplemental Figure 4). As we discuss below, other factors may also contribute to the differences in kc across the two system geometries, including hydration levels. For example, in the low-water system, we observed a dramatic increase in fluctuations and concomitant decrease in kc.

The increase in rigidity for the DMPC + cholesterol system, determined as the ratio of kc between DMPC + cholesterol and DMPC, shows excellent agreement for the high-water (hw) case, 1.8 for flat-patch versus 1.9 for vesicles. For both cases (pure lipid and + cholesterol), the vesicle’s kc is smaller than the corresponding flat-patch (1.3 × 10–19 J versus 2.1 × 10–19 J for DMPC and 2.4 × 10–19 J versus 3.8 × 10–19 J for DMPC + cholesterol).

In both the flat-patch and vesicle systems, the calculated kc value for DMPC is an order of magnitude greater than that determined experimentally (0.58 × 10–20 J, determined from diffuse X-ray scattering on oriented bilayer stacks[45]). This overestimate of kc is known for the MARTINI force field with DMPC.[2] However, given that our goal is to establish a reliable and robust (i.e., force-field independent) algorithm, we consider the agreement between the established flat-patch method and our new vesicle method as the relevant measure of success, rather than close agreement between our result and experiment. That said, we do note that there is much better agreement between the kc calculated from simulation and the experimental values[45] for both DMPC low-water (6.1 × 10–20 J versus 5.8 × 10–20 J) and the DMPC + cholesterol systems (2.4 × 10–19 J versus 2.73 × 10–19 J, respectively). Achieving better agreement with experiment in these values remains a challenge for force-field development, though it remains to be seen whether any force field (even fully atomistic ones) will be able to accurately capture kc across all lipid types.

Because of the sensitivity of the spectrum to the arc-length filter cutoff (see Figure 3), it was important to evaluate the sensitivity of the kc fit to the loss of undulation intensity. We did so by varying the number of degrees used in the data-fitting, testing both the 1.5 nm–1 and 2.5 nm–1 filters. Supplemental Figure 5 illustrates the sensitivity of kc fits and time course of kc for both DMPC and DMPC + cholesterol vesicle systems. The increased linear region for the 2.5 nm–1 spectra results in a much broader span of converged kc fits for both vesicle systems (Supplemental Figure 5B). This added robustness confirms our choice of filter parameters, specifically the 2.5 nm–1 arc-length filter (for rendering the initial surface) and the 1.5 nm–1 ideal filter (to define the URS and truncate the SPH coefficients for the kc-fit at the transition from continuum to molecular length scales). The evolution of kc shows very small changes across the duration of the simulations for both DMPC and DMPC + cholesterol systems (Supplemental Figure 5C).

CONCLUDING REMARKS

We have developed a method to determine membrane structure (number density and area per lipid) and bending rigidity from MD simulations of lipid vesicles. Comparisons of number density profiles determined from vesicle simulations using either a global or a local reference frame highlight the broadening effect introduced by local bilayer fluctuations, similar to what was previously observed in large flat-patch simulations.[3]

Using a URS removes the fluctuation-induced broadening effect and recovers the underlying structure profile, even under extreme vesicle deformations (e.g., elliptical vesicles). Our new SPHA method can be applied to any MD simulation where the topography is such that radial lines intersect the topography only once (e.g., vesicles, bubbles, micelles). In systems where the topography is more tortuous (e.g., where radial lines intersect the topography multiple times) the SPHA method breaks down. Nevertheless, the local structure profile can still be resolved via the URS, as long as the method used to define the URS is modified to accommodate the more complex topography.

Our choice of using an equi-angular grid in order to define the URS as well as in implementing the SPHA was done to simplify the calculation of the SPH coefficients via matrix transformation (eq 12). Alternative surface definitions (e.g., geodesic gridding) are equally viable; however, they may introduce increased complexity and computational cost in determining the transformation matrix. Any predefined θ, ϕ surface reference frame requires either a real-space filter (e.g., arc-length filter) or a direct spherical harmonic transformation that requires treating every membrane bead as a radial delta function to obtain the URS. Although the latter is the direct corollary to our flat-patch direct Fourier method,[2] the increased computational cost of the explicit SPH transform makes that approach prohibitive.

The spherical harmonics analysis of the URS allows us to determine the AL for the equivalent ideal sphere that samples the same area as the undulating bilayer. We define a unique AL for the vesicle as well as both monolayers. These three AL metrics describe structural imbalances that exist due to initial vesicle construction (e.g., water or lipid imbalances across the bilayer). The convergence of these three metrics may serve as a simple readout for vesicle equilibration.

This paper lays the groundwork for an iterative process to improve vesicle construction. With flat-patch membranes, numerous studies have been successful in matching experimental structure profiles by tuning the AL through altering the periodic box dimensions.[5,18,20] The vesicle membrane geometry raises a new and more difficult challenge: obtaining AL, structure, and bending rigidity correct by varying the lipid distribution in the two monolayers and the water density inside the vesicle.

Our algorithm can extract these structural and mechanical properties to guide the construction and equilibration of these complex vesicle systems. Comparison of these structural and mechanical results with experimental measurables is the ultimate goal. Although we are currently limited by the accuracy of the simulation force fields, changes in structure that correlate with changes in rigidity may provide insight into the interpretation of experimental results.