PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Phys Chem B. Author manuscript; available in PMC 2010 April 2.
Published in final edited form as:
PMCID: PMC2666921
NIHMSID: NIHMS101065

A Solvent-Free Lipid Bilayer Model Using Multiscale Coarse-Graining

Abstract

The multiscale coarse-graining (MS-CG) approach developed in our previous work is extended here to model solvent-free lipid bilayers. The water (solvent) molecules are completely integrated out of the coarse-grained (CG) effective force field. The MS-CG force field, a sum of pairwise central forces, accurately approximates the many-body potential of mean force in the coarse-grained coordinates. It thus incorporates both energetic and entropic contributions. To improve the stability and elastic properties of the MS-CG simulated bilayer, an additional constraint was adopted: the partial virial associated with CG bilayer sites was matched to its corresponding atomistic value for each configuration of the system. The resulting solvent-free MS-CG model reproduces a liquid-state lipid bilayer with accurate structural and elastic properties. Finally, the solvent-free MS-CG model is used to simulate a very large, flat bilayer and two liposome geometries, demonstrating its greatly enhanced computational efficiency.

1. Introduction

Lipid bilayers are among the most important and diverse self-assembled supramolecular biological structures. The lipid membrane serves as an interface between the cell and its environment, so it is involved in many key cellular processes spanning an exceptionally wide range of length and time scales.1,2 Certain information on the properties and functions of realistic biological membranes can be obtained from experiment. As computational power increases, however, numerical simulation has become an important method to explore membrane dynamics and structure.313 Accurate all-atom simulations of molecular dynamics (MD), for example, have been successfully used to model local membrane structure and associated phenomena.

All-atom MD methods are still limited to fairly small membrane samples and short time scales (tens of nanometers wide and at most a few hundred nanoseconds).14 Such simulations are thus not able to access mesoscale phenomena such as membrane self-assembly, domain formation in multi-component membranes, and membrane fusion or rupture. Even some “macroscopic” elastic properties of the membrane such as the bending modulus are difficult to determine using all-atom MD simulations.10

Fortunately, atomistic resolution is not of principal importance when the goal is to gain insight into mesoscale membrane phenomena. Such phenomena can be simulated using “larger”, i.e., coarse-grained (CG) structural units representing groups of lipid atoms, whose dynamics are modeled through effective interactions. CG simulations are becoming increasingly popular1517 because they remain particle-based while greatly reducing the computational expense. Another advantage of CG methods is that the particle groups can be constructed at various resolutions, permitting the study of membranes at multiple scales. The computational efficiency of CG models can be further enhanced with a solvent-free description,14 meaning that the lipids are controlled by effective solvent-mediated interactions.

The development of solvent-free CG models for ordered phases of lipid bilayers (particularly the lamellar phases) has proven difficult, and the addition of cholesterol complicates matters further. The primary obstacle is the many-body character of the hydrophobic interaction, which is a driving force of phase separation in lipid systems (e.g. bilayer assembly). Even for simple hydrophobic solutes like methane, many-body effects are of importance to the hydrophobic interaction.

Quantitatively, the importance of many-body effects can be seen in how well the many-body potential of mean force (PMF), which describes the thermodynamics of solute association, can be approximated by a sum of pairwise PMFs.18 (The latter govern solute association in sufficiently dilute solutions.) The accuracy of any coarse-grained model is determined by how well the potential surface in CG coordinates represents the free energy surface of a fully atomistic simulation in the same coordinates. Previous methods have chosen to parameterize the CG interactions indirectly. In other words, a suitable (in most cases pairwise) form is chosen for the interaction then fit to a set of system observables such as structural, thermodynamic, or elastic properties.1927 However, tuning an ad hoc potential to fit even a broad range of macroscopic system properties does not guarantee that the result approaches the “true” effective interactions as determined by integrating out atomistic degrees of freedom in a thermodynamically consistent coarse-graining procedure. This limitation is manifest in the difficulty of casting relevant physics into CG models based on simple analytical (e.g. Lennard-Jones) potentials. For example, such models show an excessive tendency to crystallize.2830

Some past attempts have been made to obtain a stable, solvent-free, fluid phase bilayer in MD simulations. Some have even used analytical potentials with many-body (e.g. density-or orientation-dependent) terms.22,28,31 Importantly, it has recently been shown that a stable fluid bilayer can be obtained using pairwise potentials alone.22,23,3234 This result suggests that effective pairwise interactions can encompass many-body effects, at least to some extent. More systematic numerical approaches seen in the literature include Iterative Boltzmann Inversion24,3539 and Reverse Monte Carlo (RMC),2427 both of which attempt to model the radial distribution functions obtained in atomistic simulations. These methods rely on tabulated representations of the effective interactions, which are more flexible than analytic functions. However, as the complexity of the system increases the RMC methods can become difficult to converge.

The “multiscale coarse-graining” (MS-CG) approach uses a sum of pairwise CG potentials to approximate the atomistic free energy surface in CG phase space. It relies on explicit atomistic interactions, and has been proven to properly accommodate many-body effects.4048 In the MS-CG method the average generalized CG force field exerted on the CG coordinates is modeled by a conservative force field (i.e. for which a potential can be introduced) with a known functional form of the potential. The best-fit MS-CG effective potential, determined by least-squares approximation, may be considered an approximation to the corresponding many-body PMF for the CG sites. Coarse-grained coordinates are typically associated with the centers of mass of atomistic groups, and span the full coordinate space. The conservative vector force field chosen for our application is a sum of central, pairwise force terms. Each term’s magnitude is dependent only on the distance between two CG sites. Note that in this approximation, the individual terms contributing to the many-body PMF are quite different from the actual pairwise PMFs.46 The former can be regarded as an effective potential describing the coupling between a specific pair of coarse-grained particles in the context of a many-body system.

The effective pairwise CG potentials can be considered as a sum of potential and entropic terms, each likewise representing an approximation to the respective many-body quantity. In the case of solvated molecular structures where the many-body PMF coordinates are explicitly chosen to be the molecular CG degrees of freedom (e.g. the centers of mass of only the solute molecules), the MS-CG potentials are in fact the solvent-free effective interactions between pairs of solute molecules. As such, they must accommodate both energetic and entropic effects arising from the solvent subsystem, which has been formally removed from the system by the coarse-graining defined in this fashion.

In the present article, we present the solvent-free MS-CG model of a mixed DMPC/Cholesterol 1:1 bilayer. This model is similar to the full-solvent model reported in Ref.42

The remainder of this article is organized as follows. Sec. 2 briefly outlines the MS-CG method and links it to a pairwise decomposition of the many-body PMF in a CG coordinate space. Sec. 3 applies the method to the DMPC/Chol 1:1 system, starting with a presentation of the MS-CG fitting procedure. Sec. 4 goes on to discuss the MS-CG interactions, examines the quality of the MS-CG bilayer model, and finally proposes several other applications, including a large bilayer, self-assembly, a bilayer with open edges, and finally liposomes. Sec. 5 presents conclusions from this study.

2. The Multiscale Coarse-Graining Method

The algorithmic development of the MS-CG method is described in detail elsewhere.40,41 Refs.47,48 rigorously demonstrate the relationship between MS-CG potentials and a pairwise decomposition of the many-body PMF in the CG coordinate space. This section summarizes the main results and extends them to solvent-free CG modeling.

We introduce an n -particle Cartesian phase space (rn,pn) with Hamiltonian h(rn,pn). The n particles are partitioned into N CG groups as a canonical coordinate transformation to the new phase space (RN (rn,pn),PN (rn,pn),xP, pxP). The Cartesian CG subspace (RN (rn,pn),PN (rn,pn)) is spanned by the atomic coordinates of N CG sites, which in turn are specified by the linear mapping operations

RI(rn)=i=1,sIcIiri,PI(pn)=i=1,sIpi,i=1,sIcIi=1,forI=1,,N,
(1)

where sI is the number of atoms in the Ith grouping. The CG coordinates RI (rn) and PI (pn) are canonical coordinates due to the last condition in Eq. (1). RI (rn) and PI (pn) are the position and momentum of the center of mass (CM) of the Ith group provided that cIi = mIi/MI, where mIi is mass of the ith atom in the group and MI is the net mass of the group. Other choices of cIi are possible. For example, setting cIi = 1/sI is equivalent to the assumption that all atoms within the group have equal mass. In this case the CG “particle” is located at the geometrical center (GC) of the group. In either case, the remaining canonical degrees of freedom (xP, pxP), which are formally integrated out in the MS-CG model, can be chosen so that N + P = n.

The local free energy at a given value RN = RN and temperature T is given by an N-body potential of mean force

W(RN)=kBTlnPR(RN),
(2)

where the probability PR(RN) of finding the system in the state with coordinates RN is

PR(RN)=δ(RN(rn,pn))RN).
(3)

The brackets denote an ensemble average over atomic variables with the Hamiltonian h (rn, pn). The Helmholtz free energy change ΔRN (t) W (here we assume that the ensemble is constant NVT, but the formalism should remain valid for other ensembles) along a path RN(t) is given by a line integral

ΔRN(t)W=RN(t)dRNFRN(RN).
(4)

The 3N- dimensional vector FRN = (FI, I = 1, N) represents the generalized forces FI (which are three-dimensional vectors) acting at the Ith generalized CG coordinates RI introduced in Eq. (1), and are directly related to the potential field W(RN) as follows

FI(RN)=RIW(RN)RN=RN.
(5)

In the MS-CG method,4048 the generalized forces FI(RN) are usually approximated by a conservative vector field (i.e., for which a potential can be introduced) GI (R1,…, RN,Ω) whose functional form depends on a set of parameters Ω. If the potential Wg ((R1,…, RN,Ω) of the model field GI is known, then the integral in Eq. (4) can be approximated by the difference in Wg at the initial and final states of the path (RN) (t). The actual forces (FI(RN) are approximated by choosing parameters Ω that minimize the residual function

χ2(Ω)=dRNPR(RN)I=1,NFI(RN)GI(RN,Ω)2,
(6)

where PR(RN) is given by Eq. (3). This form has the property that the exact form of FI(RN) are typically not known in the CG coordinate space.

Next, we assume that GI (R1,…, RN,Ω) is a linear functional of the parameters Ω. It can be shown that this property allows us to find the optimal parameters Ω0 by minimizing a new functional [chi]2(Ω), which relies on the difference between instantaneous and model forces.47 For a discrete sampling of the atomistic phase space, [chi]2(Ω) can be written as

χ¯2(Ω)=1LNl=1,LI=1,NFI({RN}l,{xP}l)GI({RN}l,Ω)2.
(7)

Here FI({RN}l,{xP}l) represents the instantaneous net force acting on the CM of the Ith CG group in the lth configuration (e.g., within an atomistic MD simulation at the lth sampled moment of time). This term can be easily evaluated from the atomistic forces fi as follows:

FI({RN}l,{xP}l)=i=1,sIfi.
(8)

The model usually chosen for GI ({RN}l,Ω) is a sum of pairwise, central forces depending only on the interparticle distances

GαI({RN}l,Ω)=βΓJ=1,MβΓgαβ(RαI,βJ,Ω)nαI,βJ,
(9)

where αI now denotes the Ith CG site of type α(in a context where sites of different types interact differently). The index Γ denotes the set of all types, so MαΓ is the total number of sites of type α. The interparticle separation is represented as a distance RαI,βJ = |RαI,βJ| = |RαIRβJ| and the unit vector nαI,βJ=RαI,βJ/RαI,βJ, where both CG sites belong to the lth configuration. This force field is conservative, and its potential is readily evaluated as

Wg(RN)=αIβJwαβg(RαI,βJ,Ω)+const,
(10)

where wαβg(R) is an integral of the gαβ(R)function in Eq. (9). The CG effective potential Wg(RN) is an approximation to the N-body PMF W(RN) used in Eq. (2). It must be emphasized that the pairwise functions wαβg(RαI,βJ,Ω), which depend only on the positions of two CG sites (αI and βJ), are not two-body PMFs. Rather, they represent an approximate pairwise decomposition of the many-body PMF.

The next stage of the present solvent-free CG procedure, determining gαβ(RαI,βJ,Ω), essentially coincides with that previously adopted for the MS-CG method.40,41 Each term gαβ(R,Ω) is represented by, e.g., a cubic spline: a piecewise cubic polynomial g(R,{Rαβ,γ,k},{gαβ,γ,k}, { gαβ,γ,k}) constructed on a mesh {Rαβ,γ,k}. The free parameters Ω to be fit are its values gαβ,γ,k and second derivatives gαβ,γ,k at the mesh points Rαβ,γ,k We have also introduced an index γ, which distinguishes between the nonbonded (γ = nb) and bonded (γ = b) interactions (defined in accordance with the connectivity between CG sites).41 Other piecewise (e.g., linear) functions can be used, as long as they are linear in their parameters.

The problem of minimizing the residual function [chi]2(Ω)in Eq. (7) can now be reduced to an overdetermined system of linear equations

γ=nb,bβΓJ=1,MβΓ(g(RαIl,βJl,{Rαβ,γ,k},{gαβ,γ,k},{gαβ,γ,k})nαIl,βJl=FαIl(CG),I=1,,MαΓ,αΓ,l=1,,L.
(11)

The solution is obtained by equating the force Fα(Γ)Il(CG) [which is FI({RN}l,{xP}l) in Eq. (7)] to the corresponding model force [Eq. (9)] calculated at each site, for all L configurations. When using the pairwise representation Eq. (9), we assume that contributions only come from those sites which belong the to set of types Γ. This system of equations (11) can be augmented by a set of constraints, as long as they are linear with respect to {gαβ,γ,k, gαβ,γ,k}.41

Because the number of configurations L should be large enough to produce reliable statistics, it is not feasible (and might be not even necessary) in practice to solve Eq. (11) simultaneously for all configurations. To make the problem tractable, the configurations are partitioned into equal blocks; L in Eq. (11) can now denote the size of each block. The size of L has to be large enough that Eq. (11) still overdetermines the force parameters. Eq. (11) is then solved for each block separately, and the resulting solutions are averaged. The final result is consistent with Eq. (6), in that FI(RN)has the meaning of a “mean” force. Block averaging is equivalent to determining a “global” mean force from the mean forces acting within each block of configurations [that is, the solutions of Eq. (11) for individual blocks] depending on the degree of correlation.47 The mean force field from block averaging might not be a global minimum of [chi]2(Ω), but it is typically much smoother than solutions obtained by direct minimization of [chi]2(Ω) (e.g. by solving the full system of equations). In this sense, block averaging may be considered a way of regularizing the solution.

Importantly, in the procedure outlined above, if Γ in Eq. (11) denotes only solute sites, the MS-CG procedure results in solvent-free MS-CG potentials between solute sites, in which the solvent effects are effectively built in to the CG potentials.

3. Building the Solvent-Free MS-CG Model

3.1. Reference All-Atom MD Simulation

A MS-CG model of the interactions in a DMPC/Chol mixture in a particular phase or thermodynamic state should be matched to an atomistic MD simulation carried out under similar conditions. For ordered systems (e.g. lamellar), such MS-CG interactions will be phase-dependent. Furthermore, for configurations beyond the reference structure there will be substantial uncertainty in the effective interactions due to unsampled MS-CG site-site separations. However, such interactions might well be relevant to other phases. In the present study, two simulations were therefore used as the source of reference atomistic forces and trajectories: one of a solvated bilayer and the second of randomly dispersed DMPC and Chol molecules in solvent. The latter (random) phase yielded a better sampling of intersite separations, but the MS-CG model matched to a bilayer simulation is clearly preferable for lamellar structures.

The atomistic lipid bilayer simulation contained 32 DMPC/Chol pairs of molecules and 1312 water molecules, while the random (solution) phase was represented by a same number of DMPC/Chol pairs solvated in 2688 water molecules, and placed in a periodic cell. For the random phase, to improve sampling, and because the DMPC and Chol molecules tended to segregate, four simulations were carried out with different initial configurations. The DMPC molecules were modeled using a united atom force field, while the Chol molecules were described by a modified AMBER force field.8 The rigid TIP3P model49 was used for the water solvent molecules. Electrostatic interactions were calculated via the particle-mesh Ewald (PME) summation,50 and van der Waals interactions were cut off at 0.7 nm.

The bilayer simulation’s initial configuration was taken from Ref.51 and additionally equilibrated for 20 nanoseconds in the constant NPT ensemble. Melchionna’s modification of the Hoover algorithm for coupling a Nosè-Hoover thermostat to an isotropic barostat52 was used to implement an NPT ensemble with temperature 308 K and zero pressure. The ratio of supercell dimensions, which is maintained constant in the simulation with an isotropic barostat, corresponded to a tensionless state of the bilayer. The random (solution) systems were pre-equilibrated for about 1 ns in a constant NPT ensemble implemented with the isotropic barostat at the same pressure an temperature as for bilayer system.

Reference trajectories and forces for the bilayer were then collected from a 50 ns simulation in the constant NVT ensemble. The solution systems were sampled for 10 ns each under the same conditions. The cell dimensions in the NVT simulation of the bilayer were set to 4.07 × 3.05 × 7.42 nm, close to equilibrium. This geometry will be used repeatedly to generate initial bilayer configurations for the MS-CG simulation, and referred to as a base (1 × 1) cell.

3.2. Coarse-Graining and the Details of Force-Matching

The DMPC and Chol molecules were mapped onto 11 different types of CG sites, as shown in Fig. 1. (The same scheme is given in Ref.40,42) This mapping results in a total of 66 distinct, nonbonded interactions. The nearest-neighbor head group sites are connected by bonds, as are all nearest-neighbor and next-nearest-neighbor sites involving either a DMPC glycerol backbone (GL), ester (E1,2), alkyl (SM, ST) groups or Chol groups as shown in Fig. 1. Representing the DMPC alkyl chain with a (CH2)3 group as a single CG site reduces the degrees of freedom available for chain disorder, but may also introduce uncertainty with respect to its shape (the site may be in the trans- shape or gauche- shape, for example). In general, this ambiguity would complicate the determination of effective interactions for the alkyl sites. Fortunately, this was not an issue for our study; at 308 K, the DMPC/Chol (1 × 1) bilayer occupies a liquid-ordered (Lo) phase where all alkyl chains are in the trans- conformation. The alkyl chains may undergo thermal fluctuations within the rotamer. However, associated fluctuations in the effective interactions between CG groups induced by such fluctuations are rather small and are expected to be captured in the MS-CG force averaging procedure.

Figure 1
Coarse-grained representation of the lipid bilayer molecules: (a) DMPC and (b) Cholesterol. The MS-CG bonded interactions were explicitly matched for the connected sites shown.

As with our previous MS-CG model of a full solvent mixed bilayer,42 the next-nearest-neighbor bonds were brought in to emulate an angle-dependent interaction (as the current MS-CG procedure can only treat two-site interactions). These bonds are referred to as angle bonds. Angle bonds were used for pairs involving the acyl tail sites (E1, E2, SM, ST). As head groups are more rigid, only nearest-neighbor bonds were introduced. This explicit separation of bonds in the MS-CG procedure increases the number of CG site-site interactions to 79.

For the purposes of force-matching, atomic trajectories, velocities, and forces were sampled at 5 ps intervals. For the bilayer system, this resulted in 10,000 stored configurations. For each configuration, the trajectories and forces (ri, fi) for all DMPC and Chol atoms were combined into corresponding data ( RI(CG),FI(CG)) for the CG sites as in Eq. (1), (8). The data (total forces and coordinates) associated with water atoms were not included from the calculation and thereby folded into the other CG interactions. MS-CG interactions for the bilayer geometry were derived using constraints on a virial associated with the bilayer sites; the goal is to match the partial pressure of the bilayer as seen in the atomistic MD. The atomistic and CG partial virials associated with the bilayer were evaluated as

Watm=13ifiri,WCG=13IFI(CG)RI(CG),
(12)

after subtracting the “drift” force Fdrift = Σifi/Nbil from the atomistic forces fi. Nbil is the number of bilayer atoms. For each configuration we require that WCG be equal to Watm+Ekin, where ΔEkin is the difference between the instantaneous kinetic energy of all bilayer atoms and the translational kinetic energy of all bilayer CG groups.

After the CG forces were calculated, the MS-CG procedure was applied as given in Eq. (11) with an explicit separation of bonded and nonbonded interactions.41 The mesh spacings used for the force splines were approximately 0.005 nm for nonbonded pairs and 0.0025 nm for bonded pairs. In the case of nonbonded forces, the mesh extended to 1 5. nm. The size of the blocks used in Eq. (11) was L = 10. The MS-CG bonded interactions were approximated by harmonic potentials, as will be discussed in the next section. To elucidate the effect of the solvent, the solvent-free and full-solvent CG interactions in a bilayer parameterized without any virial constraint and using a coarser distance mesh were compared. The MS-CG interactions from solution were derived with no virial constraint and used to simulate aggregation of the DMPC/Chol mixture from a randomly dispersed solution. The bilayer MS-CG interactions are difficult to extend to random configurations.

The MS-CG forces, especially for hydrophilic groups, appeared to converge slowly at large separations. The final force profiles had substantial noise but apparently zero background. The noise affected the quality of simulations under constant pressure conditions; for example, the equilibrium of the MS-CG bilayer area was slightly off. As a remedy, we reduced the cutoff in bilayer MS-CG interactions; this change eliminated some noise and at the same time preserved the correct equilibrium area under zero tension. The forces were smoothly shifted to zero at the cutoff to avoid discontinuities. A cutoff of 1.2 nm appeared to be appropriate.

4. Results and Discussion

4.1. MS-CG Interactions

Solvent-free and full-solvent effective potentials were tabulated as integrals of the MS-CG forces, and defined to be zero at the cutoff radius. In Fig. 2, selected solvent-free potentials for the bilayer geometry are plotted along with those derived with full-solvent interactions. This comparison provides insight into water’s stabilizing effect on interactions within the bilayer. For hydrophilic groups the solvent-free potentials are also expected to be more attractive and less intense due to embedded solvent screening effects. However, the differences between the two potentials are rather mixed. For instance, many of the head group sites interacted more strongly in the solvent-free formalism.

Figure 2
Comparison of selected solvent-free (solid) and full solvent (dashed) MS-CG nonbonded potentials as a function of the CG inter-site separation. See Fig. 1 for CG site definitions.

It is possible that behavior such as that described above arises from peculiarities of the MS-CG methodology, which relies on correlations between the total atomistic forces associated with CG groupings and configuration changes to determine the effective interactions. Thus, any significant correlation between forces and trajectories at a particular separation may be interpreted by the MS-CG procedure as an effective interaction. The interaction will be deemed stronger for more correlated data. Correlation can be introduced by cooperativity between atomic groups, even if the groups do not interact directly. In condensed phases, such cooperativity might develop from environment-mediated interactions. Solvent-mediated interactions (induced by hydrophobic or screening effects) and water bridging between hydrophilic groups are both examples of such interactions.

In the solvent-free MS-CG model, interactions are affected not only by the cooperative dynamics of solvent degrees of freedom (which are integrated out) but also by cooperativity in the bilayer’s explicit degrees of freedom. For example, while solvent screening softens the solvent-free MS-CG potentials between head groups, this effect can be masked by the contributions from cooperative bilayer modes (e.g. undulation). This observation might explain why some hydrophilic sites interacted more strongly in the solvent-free model than in the explicit solvent model, as well as the noticeable shift in interaction energies for sites buried in the bilayer. The solvent-free potentials appeared to be shifted to higher energies for many site pairs; most belong to the DMPC and Chol head groups (e.g. CH-CA, PH-CA, PH-PH), but a few belong to the hydrophobic region (e.g. SM-CC). For sites positioned deeper in the bilayer most potentials were shifted to lower energies. The E1-CA and CB-CB pairs were among those sites, which showed the greatest increase in attraction (> 10 kJ/mol). The solvent-free interactions involving lipid alkyl sites and Chol CC, CD sites include smaller contributions from the water subsystem (< 4 kJ/mol), but their net effect on solvent-free bilayer stabilization could be substantial because there are a great many lipid tail sites. For this reason the E1-SM interactions, which became attractive by about 6 kJ/mol in the solvent-free system, likely contributing significantly to bilayer stability.

That DMPC molecules were attracted most strongly through CH-PH interaction is perhaps anticipated, as the CH and PH groups bear the largest partial charges of opposite sign. The DMPC and Chol molecules were mainly attracted through E1-CA interactions, and sterol molecules were mainly attracted through CB-CB interactions. All three interactions have strengths of about −20 kJ/mol. The atomistic description of the E1-CA bond that gives rise to DMPC-Chol coupling is a hydrogen bond between Chol OH hydrogen and ester carbonyl oxygen. This picture agrees with observations from both simulation and experiment that at higher cholesterol concentrations the hydroxyl groups of Chol molecules tend to remain close to the sn-1 carbonyl groups.11,5355

In addition to hydrogen bonding, the E1-CA interaction accommodates contributions from water-bridged oxygen pairs. In descending order of strength, we see interactions between Chol CA sites and the following DMPC sites: E1 (−21 kJ/mol at the potential well minimum), PH (−9 kJ/mol), CH (−4 kJ/mol), and E2 (−1.4 kJ/mol). These interactions, which in the atomistic description were due to the hydrogen bonding and water bridging, are all very important in the formation of DMPC-Chol complexes. The alteration of CH-CA with E1-CA or PH-CA bonds may lead to oligomerization, creating a linear network of bonded molecules (DMPC-Chol-DMPC-Chol…).56 The interaction between nonpolar parts of the Chol molecule and alkyl chains reduces the probability of trans- conformations, which may lead to a “condensing effect”. (Specifically, it would create large negative deviations in the additivity of molecular areas for unlike molecules in the bilayer.)

The interaction of hydrophobic sterol sites CB, CC, CD with the alkyl SM, ST sites were weakly attractive (<−2 kJ/mol), but sufficiently strong to keep the alkyl chains in trans- conformations. The CA-SM interaction (not shown), however, was repulsive. The coupling of sterol molecules through CB groups competed with their repulsion by CA groups. In fact, at distances of 0.6–1.0 nm the repulsive CA-CA potential strongly resembled a weaker version (by ≈7 kJ/mol) of the inverted CB-CB potential. Adding both interactions to the mildly repulsive CA-CB potential (not shown) results in an extended region of weak attraction between Chol molecules. Considering their steep repulsion at separations of < 0.5 nm, we might expect this weak attraction to bias the lateral organization of Chol molecules at high concentrations into quasiperiodic superlattices with a lateral separation varying roughly from 0.5–1.0 nm.57 This kind of organization can also be inferred from the CA-CA distribution function shown in Fig. 4 (the CB-CB and CC-CC distributions qualitatively are similar to the CA-CA distribution)

Figure 4
Comparison of selected CG site-site RDFs from atomistic (solid) and solvent-free MS-CG (dashed) simulations. See Fig. 1 for CG site definitions.

Force profiles obtained from the bilayer and randomly dispersed phases are only qualitatively similar. A major cause of the observed differences is methodological: the MS-CG minimization procedure targets the total force experienced by a group in the atomistic simulation. The resulting pairwise decomposition is therefore sensitive to variations in sampling.

Bonded interactions exhibited much less difference under the solvent-free MS-CG treatment; Fig. 3 compares some of the bonded force profiles in detail. The largest shift in energy (less than 2 kJ/mol) was for the CH-PH bond. Less solvated bonds were affected even less by the solvent-free treatment. The bonded interactions deviated from being harmonic, with softer parts of the molecule experiencing the strongest deviations. Rigid bonds such as those between CA, CB, and CC sites of the Chol molecule were almost perfectly harmonic (e.g., the CA-CB bond shown in Fig. 3). Bonds within flexible parts of the molecules (mostly bonds within the DMPC acyl tails and bonds involving aliphatic CC or CD sites of the Chol), however, were essentially anharmonic. The angle bond forces for acyl chains were always repulsive, while the bonds between nearest-neighbor sites had one equilibrium state. In the solvent-free MS-CG model, the bonded force profiles were approximately linear (i.e., harmonic) in a least-squared sense around the equilibrium position, as shown in Fig. 3. The highest-frequency bond in the DMPC molecule (GL-E1) was about 310 cm−1, that is only slightly larger than the GL-E2 frequency. In the Chol molecule, the CA-CB bond had the highest frequency at about 230 cm−1. These low-lying CG bond spectra permit a substantially larger (about an order of magnitude) time step in the MS-CG simulations compared to atomistic MD. The properties of DMPC/Chol mixtures in the MS-CG simulations also appear insensitive to variations in the bonded force constants, at least within the uncertainty of our meshed bonded force profiles.

Figure 3
Effective MS-CG bonded forces between selected CG sites in the solvent-free (solid) and full solvent (dashed) models, as a function of inter-site separation. The straight (dot-dashed) lines show harmonic fits used in the solvent-free MS-CG simulations. ...

4.2. Model Validation

The performance of the solvent-free MS-CG model was verified by evolving a system with the same geometry and thermodynamic conditions as the reference atomistic simulation. The timestep of the MS-CG MD simulation could be increased to as high as 0.01 ps (we used 0.008 ps). The CG image of the final configuration of the reference atomistic MD simulation with the water subsystem removed was chosen as the initial bilayer configuration for the CG MD simulation.

First, the CG bilayer was simulated in a constant NVT ensemble for 50 ns after 20 ns of equilibration. Site-site RDFs can be regarded as the most important test of model validity, as their accurate reproduction implies that the atomistic free energy landscape has been captured sufficiently well. The solvent-free MS-CG bilayer was stable, and its site-site structure agreed well with the reference all-atom MD simulation. This result is demonstrated in Fig. 4, which compiles some typical RDFs. The hydrophilic parts of the MS-CG bilayer exhibited less than perfect agreement, owing to the many-body character of water-mediated interactions. This difference can be attributed partly to our use of pairwise MS-CG potentials, and partly to large fluctuations in the atomistic forces associated with strongly solvated groups. The CA-CA structure has been somewhat enhanced relative to the atomistic model (the same is true of the CB-CB structure, which is similar in behaviour and therefore not shown), implying a more ordered lateral arrangement of the sterol hydroxyl and aromatic groups in the MS-CG simulation. The MS-CG structure of Chol aliphatic sites (as seen in the CD-CD structure) was in good agreement with the atomistic simulation. The RDFs between Chol sites of the same type exhibited a characteristic two-hump shape. The ratio of their two favored positions resembles that of a closely packed 2D array, suggesting the possible lateral ordering of Chol molecules.

Fig. 5 compares several density profiles in the bilayer normal direction, taken from the solvent-free MS-CG and atomistic MD simulations. Bilayer thickness is defined as the normal distance between peaks in the PH (phosphate group) distribution. In the solvent-free and atomistic simulations, these thicknesses are 4.16 nm and 4.26 nm respectively. According to experiment, the barycenter of Chol ring A is vertically positioned at ≈1.5 nm from the bilayer center. This distance is nearly invariant with temperature (within 0.1 nm).53 The CA density maximum was 1.46 nm in the MS-CG simulation, agreeing well with the above result (the CA is close to center of Chol ring A).

Figure 5
Number density profiles in the bilayer normal direction at selected CG sites, comparing the solvent-free MS-CG (solid) and atomistic (dashed) simulations. See Fig. 1 for CG site definitions.

The dynamics of the solvent-free bilayer was accelerated compared to the reference all-atom MD simulation as expected. The lateral MSDs of molecular CMs obtained from 50 ns long simulations, are compared in Fig. 6. The coefficients of lipid lateral diffusion, however, were DL = 2.6 × 10−12m2/s in the atomistic MD simulation and DL = 2.0 × 10−11 m2/s in the MS-CG simulation. The atomistic value agrees favorably with experiment (pfg-NMR method),5860 which yields DL = 3.0 × 10−12 m2/s. The lateral diffusions of DMPC and Chol molecules in the atomistic MD simulation were essentially equal, in agreement with previous simulations of high-Chol content bilayers.55 As seen in Fig. 5, the MS-CG model reproduced this feature well. The close similarity of DMPC and Chol lateral diffusion constants can probably be attributed to a high probability of coupling between Chol and DMPC molecules, which occurs through bonding of a Chol CA site to the DMPC E1, PH, or CH sites.

Figure 6
Mean-squared center-of-mass displacements of DMPC (solid) and Chol (dashed) molecules, comparing the solvent-free MS-CG and atomistic simulations.

We next simulated the MS-CG bilayer under tensionless conditions for 100 ns. The zero-tension state was emulated by coupling bilayer in-plane (xy) cell dimensions to an anisotropic barostat; the z dimension was kept fixed. The bilayer remained stable, but the ratio of xy cell dimensions relaxed from a value of ≈0.9 (in the reference atomistic MD simulation) to ≈0.75 in the tensionless state. The bilayer area was 0.772 nm2 per DMPC/Chol pair in the MS-CG simulation, and 0.776 nm2 in the atomistic MD simulation.

The average configurational MS-CG energy per DMPC/Chol pair was −215 kJ/mol. This value includes the entropic contribution from missing atomistic degrees of freedom in the MS-CG description.

4.3. Simulation of a Large Bilayer

The simulations described here model a large MS-CG bilayer, which in most cases is a 15×15 replica of the reference system (7200 DMPC/Chol pairs) under zero-tension conditions. Where explicitly stated, a 2×2 (128 DMPC/Chol pairs) replica is sometimes used as well. All the atomistic simulations referred to in this section were carried out under constant NPT conditions using an isotropic barostat. In addition to the original reference system, 2×2 and 9×5 (1440 DMPC/Chol pairs) atomistic replicas were simulated for 30 ns.

Note that MS-CG simulation of a 1 × 1 solvent-free bilayer was about three orders of magnitude faster than the corresponding fully atomistic MD simulation. In lager systems, the gain in computation time will be even higher.

A. Large Bilayer Structure

Fig. 7a shows a side view of the MS-CG and atomistic bilayers in the CG representation, with molecules projected onto their centers-of-mass (CMs) top. Views of the upper leaflets are given in Fig. 8. Lateral diffusion in the MS-CG bilayer was by about an order faster than the reference atomistic 1×1 value (reported in Sec. 4.2) and its structure displayed no long-range lateral molecular order. Together with the fact that all the alkyl chains were visibly in trans- configurations, these observations allowed us to definitively identify the phase as liquid-ordered (Lo). Undulations in the bilayer can be discerned in Fig. 7c. Fig. 8a, a top-down view of the monolayer with molecules projected onto their CMs, clearly shows that the DMPC and Chol formed thread-like complexes. Most were of the A-B-A- B-… type, but pure chains of DMPC or Chol molecules could also be found. The threads tended to be oriented along a longer side of the cell lateral plane. A similar pattern can be seen in Fig. 8b, a snapshot of an atomistic 2×2 bilayer in the CM representation. Lateral oligomerization in PC bilayers with high Chol content tends to create networks of such complexes. This process has been suggested to occur through alternating hydrogen bond OH···O (between a hydrogen of the Chol hydroxyl group and preferably a carbonyl oxygen of the DMPC molecule’s ester groups) and hydrogen bond CH···OH (between the CH3 in the choline of a PC molecule and an oxygen in a Chol hydroxyl group) or CA-E1,2 and CH-CA bonds in the CG description, and a “condensing interaction” involving DMPC alkyl chains and nonpolar parts of the Chol molecules; the latter contributes to the stability of such complexes.56,61 The Chol molecules are likely to form a lateral superlattice with local quasi-order resembling a closely packed 2D structure, as seen in Fig. 8c. As previously discussed, this phenomenon is also manifested in Chol CA-CA (see Fig. 4), CB-CB, and CC-CC RDFs. An indication of the formation of a Chol superlattice is a vacancy-like mechanism of Chol diffusion often observed in the MS-CG simulation as demonstrated in Fig. 9.

Figure 7
Snapshots of a membrane simulated using the solvent-free MS-CG (15 × 15) [Panels (a), (c)] and reference atomistic simulations [Panels (b), (d)]. Panels (c) and (d) show the DMPC (green/large balls) and Chol (grey/small balls) molecular CMs; Panels ...
Figure 8
Snapshots of a membrane leaflet (top view, down the + z axis) modeled using solvent-free MS-CG [Panel (a)] and atomistic (2 × 2) [Panel (b)] simulations with DMPC (green/large balls) and Chol (grey/small balls) molecules projected onto their centers-of-masses ...
Figure 9
Bilayer-projected CM trajectories of two Chol molecules over a period of 15 ns of MS-CG MD simulation, demonstrating the vacancy assisted character of Chol lateral diffusion.

According to experiment, a DMPC/Chol 1:1 bilayer does not form a gel phase.62 However, an energetically optimal structure of the DMPC/Chol 1:1 bilayer has been studied.63 To further explore these results, we first simulated a MS-CG bilayer in the 1×1 arrangement at T = 250 K in a tensionless state. After an ordered structure was detected in this first monolayer, the system was replicated to create a 2×2 bilayer and the simulation was continued. Snapshots from the latter are given in Fig. 10. The 2×2 bilayer exhibited lateral order, which however was more pronounced in the one leaflet. The other 2×2 bilayer’s leaflet apparently had some difficulty finding a minimum in the energy landscape created by the more ordered monolayer. The dynamics of the more ordered monolayer was crystalline; even after a long simulation the periodicity inherited from the 1×1 structure could be still seen (Fig. 10a). The less ordered monolayer exhibited faster dynamics and appeared to be in the fluid state. The packing in the more ordered monolayer resembled that in structure A from Ref.63, with two DMPC-Chol dimers per asymmetric unit (similar to the tetramer squared in Fig. 10a). Structure A possesses cm planar symmetry with a centered primitive cell. In such geometries the nearest neighbor is always the same type of molecule.

Figure 10
Snapshots from a 2×2 solvent-free bilayer simulation carried out at 250 K and zero tension. Panel (a): Packing of DMPC (green/large balls) and Chol (grey/small balls) molecules is shown by projecting their CMs onto the monolayer (top view, down ...

Many distorted, hydrogen-bonded tetramers appeared in the MS-CG fluid membrane (see Fig. 8a), contributing to its thread-like structure. This form of packing agrees with the observation that these superlattices have a centered rectangular symmetry.64 Such superlattices may form in systems with both long-range attraction and short-range repulsion. In the MS-CG description, a kind of Chol superlattice is formed by the weak but relatively long-range (0.5–1.0 nm) attraction of sterols, as discussed in Sec. 4.1, augmented by the PC-sterol directional coupling. This long-range attraction competes with a short-range repulsion induced in part by steric strain in the acyl chain matrix.

Several theoretical and experimental studies have suggested that sterols (e.g., cholesterol and ergosterol) tend to adopt regular, superlattice-like distributions in lipid bilayers, and PC bilayers in particular.57,64,65 It has been proposed that cholesterol imposes steric strain in the acyl chain matrix, which drives cholesterol molecules towards a regular distribution.64 The same study provided evidence that some of the putative superlattices have a centered rectangular, rather than hexagonal, symmetry. It should be noted that a substantial number of experimental studies have found deviations in some physical properties of the membrane at Chol fractions such as 1:5, 1:3 and 1:1. These “magic” Chol concentrations have often been interpreted as evidence of specific Chol-PC complexes. It has also been reported that the effect of cholesterol on the thermotropic behavior of PC/Chol bilayers is remarkably well explained by the cholesterol superlattice model.

B. Bilayer Elastic Properties

The good accuracy of the MS-CG bilayer model in a zero-tension simulation might be taken as evidence of reasonable elastic properties such as area compressibility, bending modulus, and bulk modulus. For this model, the atomistic properties were obtained from a 9×5 replica. In addition, elastic moduli were evaluated for the full-solvent MS-CG model42 of a bilayer of the same size. The area compressibility KA was obtained using a fluctuation formula, KA = kB TA/< δA2 >, where A is bilayer area. The bending modulus kb was evaluated using a continuum limit for the membrane’s structure factor S(q) in a tensionless state at small wavelengths: S(q) = kBT/(Akbq4). In this case, the structure factor showed q−4 behavior at q≤ 0.3 nm−1.

To calculate a bulk modulus, the barostat was coupled only to the z -dimension of the cell (equivalent to the application of pressure in that direction only) while keeping the other dimensions unchanged. The average bilayer volumes V1 and V2 were then measured at pressures P1 = 0 and P2 = 0.1 kbar. The bilayer volume was calculated as the product of bilayer area and thickness, where the latter was assumed to be the separation between maxima in z-density at a PH site (see Fig. 5). The bulk modulus E was then evaluated as the inverse of compressibility: 1/E = β= ln(V1/V 2)/(P2−P1).

The elastic moduli for the MS-CG models are compared to their corresponding atomistic MD and experimental values in Table I. The atomistic moduli are well reproduced by the solvent-free and full-solvent MS-CG models. The atomistic bending and bulk moduli agree reasonably well with experiment. The atomistic and MS-CG model values of area compressibility were substantially higher than the experimental values, however. This discrepancy may be attributed in part to the inaccuracy of the fluctuation formula.10 More accurate methods are available to estimate the area compressibility. However, in the present study our main goal was a comparative study of the bilayer properties in the atomistic and CG representations, when the property is measured using the same methodology for both models.

TABLE I
Elastic properties of the different membrane models: E, bulk modulus; KA, area compressibility; kb, bending modulus.

C. Bilayer Stress Profile

Next we compared the local stress profiles P(z) for solvent-free MS-CG and atomistic MD systems. P(z) was calculated according to the scheme of Lindahl et al.,66 as modified by Gullingsrud et al.,67 with a bin size of 0.05 nm. Both atomistic and MS-CG simulations were carried out for 30 ns on a 2×2 replica of the reference system. In the atomistic simulation, the electrostatic contribution to local stress was calculated using a cutoff of 1.8 nm. The local pressure Ploc (z) = (Pxx (z)+Pyy (z)+Pzz(z))/3 and lateral pressure Plat (z) = Pzz (z)−(Pxx (z)+Pyy(z))/2 profiles are shown in Fig. 11.

Figure 11
Analysis of membrane pressure profiles obtained from the solvent-free MS-CG (red) and atomistic (black) simulations. Panel (a): Local pressure, Ploc (z); Panel (b): Lateral pressure, Plat (z).

The shape of the atomistic Plat (z) profile agreed with that obtained for the DPPC/Chol 1:1 bilayer,68 but the amplitude of the DMPC/Chol Plat (z) was higher (3.6 kbar vs. 2.5 kbar for repulsive regions, and −4.0 kbar vs. −3.5 kbar for attractive regions). Interestingly, the atomistic Plat (z) has a strong compression region (≈1.0 nm) located close to a maximum in the average density of alkyl (SM) sites (see Fig. 5), which in turn was close to the center of the region spanned by the Chol (which is close to maximum in Chol CB density). This structure is very different from that of a pure DMPC bilayer, where a strong compression region separates the hydrophilic and hydrophobic layers.67,69 With respect to the atomistic Ploc (z), a strong negative region developed in between two different hydrophobic regions (≈1.5 nm) that is close to maximum in Chol CA density. A shift of the main compression region toward lipid tail sites might be explained by a condensing effect of Chol on DMPC chains, which would increase the rigidity of the bilayer’s hydrophobic core.

Qualitatively, the MS-CG and atomistic Ploc (z) and Plat (z) profiles were similar. In particular, for both representations positions of minima(maxima) in the Ploc (z) profile corresponds to a position of maxima(minima) in the Plat(z) profile. In the MS-CG simulation, the region of negative tension in Plat (z) occurred farther out. The position of the outermost negative Plat (z) region (≈ 2.0 nm), which probably developed due to strong solvation of the CH and PH groups (see Fig. 5), was reproduced sufficiently accurately by the MS-CG model. However, this feature was much shallower than in the atomistic model. The amplitude of the MS-CG Ploc (z) profile was close to the atomistic data, but Plat (z) was less intense overall. Possible reasons for this discrepancy in the MS-CG include an absence of local stresses due to the water subsystem and bilayer intra-group (bonded and nonbonded) interactions, softer MS-CG bonded interactions, and our harmonic approximation to the bonded forces. (With respect to the last approximation, there is indeed a pronounced difference in the stress profiles calculated using the atomistic sequence of configurations and the MS-CG force field assuming harmonic bonds.)

The accurate reproduction of local stress at every location is not necessarily of great importance. More important are the integral characteristics of stress taken over the whole bilayer (pressure and tension) or key regions thereof (head groups and the hydrophobic core). The first and second moments of Plat (z) are also important, because they determine the bending properties of the bilayer. This property is demonstrated below.

Let us calculate the excess free energy of a small spherical deformation in a fluid state membrane. The deformation has curvature c = 1/R, where R is its radius. The excess free energy is Kc2/2, where K = kb + k/2; kb is the bilayer bending constant, and k is a saddle-splay constant.70 The elastic constants kb and k can be expressed using moments of Plat (z) as

kbc0=zPlat(z)dz,k¯=z2Plat(z)dz,
(13)

where c0 is the bilayer’s spontaneous curvature and the second integral is centered on the neutral plane (which for a symmetric bilayer is its mid-plane). The same expressions hold for leaflets, except that the integrals are narrowed to a half-cell region in the z-direction. For either a leaflet or an asymmetric bilayer, c0 can be nonzero. Correct reproduction of elastic constants and the spontaneous curative of bilayer leaflets are important criteria if we are to study phenomena such as membrane-protein interaction, where the constants enter into the chemical potentials of embedded proteins.71 Positive values of K predict a stable bilayer, while negative values suggest the spontaneous formation of spherical aggregates such as liposomes. Using values of kb from Table I (assuming a bending rigidity of kb/2), the spontaneous curvatures of the MS-CG and atomistic MD monolayers are −0.30 nm−1 and −0.27 nm−1 respectively. The MS-CG and atomistic values of k for bilayers were 0.9× 10−19 J and 1.0 × 10−19 J respectively. For both MS-CG and atomistic bilayers the K modulus was positive, pointing toward stability of the bilayer when subjected to positively curved perturbations.

4.4. Lipid Aggregation in Solution

In this simulation, evolution of the system from a random dispersion phase was studied using MS-CG interactions matched to the atomistic MD simulation of a dispersed solution. The initial system consisted of 800 pairs of DMPC/Chol molecules distributed randomly in a cubic supercell of length 36 nm. The system rapidly (in about 10 ns) formed either DMPC/Chol bilayered structures of quasidiscotic shape or aggregates composed mostly of Chol molecules. Typical DMPC/Chol aggregates are shown in Fig. 12. The bilayered aggregates were about 6–7 nm in diameter, but were not exactly “bicelles” as the hydrophobic core was exposed at the boundaries. After this initial rapid formation of aggregates, the system showed little change. It is therefore unlikely that these models are completely transferable to a bilayer configuration without additional modifications in the CG interactions.

Figure 12
Typical self-assembled structures of DMPC and Chol molecules in solution (not in a pre-assembled bilayer configuration), from a solvent-free MS-CG MD simulation.

The differences between the MS-CG models fitted to the ordered lamellar and disordered phases highlights the phase dependence of the effective MS-CG interactions, which should be taken into account when a model is constructed to study a particular phenomenon. For instance, if the main goal is a n accuate reproduction of the self-assembled ordered phase, then the best choice is to use the MS-CG interactions fitted to the ordered state. At early stages of self-assembly, when the system is mostly disordered, such a model may produce a bias toward excessive segregation. The simulation of self-assembly might be improved by using the MS-CG interactions fitted to disordered state at early stages of self-assembly and, then, as the system evolves toward the ordered state gradually “switching” them to the interactions derived from the lamellar phase. This, and others approaches to enhance the transferability of MS-CG models between thermodynamic conditions are currently being explored in our group.

4.5. Lipid Bilayer with Open Edges

A solvated membrane with free boundaries was simulated as a 15 × 15 bilayer in free space. Fig. 13 shows a snapshot of the membrane after 30 ns of simulation. The membrane appeared to be rather stable, but resisted curving. Its deviation from planar geometry was rather limited, indicating a rigidity of the bilayer structure. Molecules at the membrane edges rearranged themselves to minimize exposure of the hydrophobic core, resulting in a bicelle-like structure. Overall, the system evolved into a disk-like bicelle as expected for lamellar structures smaller than the size needed for vesicle formation. This process was quite slow, however; evidence of the initial square shape persisted throughout the simulation (several tens of nanoseconds), suggesting low line tension along the edges.

Figure 13
Snapshot of a lipid bilayer with free edge boundaries [Panel (a)] and the edge of the membrane [Panel (b)], from the MS-CG MD simulation.

4.6. Liposome Simulation

Our final illustrative application of the solvent-free MS-CG method is the liposome geometry. Liposomes of two sizes were simulated using MS-CG interactions fitted to the reference atomistic bilayer simulation. One liposome was 40 nm in diameter and was composed of 11,960 DMPC/Chol pairs. The other was 100 nm in diameter and formed of 79,282 DMPC/Chol pairs. Including solvent (assuming box size equal to liposome outer diameter), these effectively represent approximately 6 4×106 and 108 total atoms, respectively. In both liposomes, the inner and outer layers initially had identical compositions. In the outer layer, the available area per DMPC/Chol dimer was close to that of a tensionless bilayer. The initial states are shown in Figs. 14a and 14c. Snapshots of the liposomes after 20 ns, with molecules projected onto their CMs, are shown in Figs. 14b and 14d. The small liposome was in a high-tension regime, and quickly relaxed to non-spherical form as shown in Fig. 14b with its shape fluctuations destroyed. Several areas of high and low curvature developed similar to those seen in Fig. 14b. The regions with high curvature had larger projected areas, and came close to developing a rupture. Molecules from high-curvature areas then migrated into adjacent low-curvature areas, which became flatter. The large liposome was stable, and showed more pronounced shape fluctuations. These simulations reflect the remarkable increase in computational efficiency afforded by the solvent-free MS-CG model.

Figure 14
Solvent free MS-CG simulation of small liposome, 40 nm in diameter [Panels (a) and (b)] and a larger liposome 100 nm in diameter [Panels (c) and (d)]. The green and black dots are the CMs of DMPC and Chol molecules, respectively. The left-hand panels ...

5. Conclusions

The MS-CG method4048 was extended in this paper to derive a solvent-free model of the mixed DMPC/Chol lipid bilayer. As shown in this and previous work, the MS-CG sum of pairwise potentials approximates the coarse-grained, many-body PMF, and thus accounts for the missing interactions and entropic contributions due to removed degrees of freedom. To improve the MS-CG bilayer’s elastic properties, the MS-CG interactions were constrained to match the instantaneous partial virial associated with bilayer atoms in an equivalent atomistic simulation. For purposes of comparison, a similar model was obtained by matching the forces to an atomistic simulation of randomly dispersed DMPC/Chol 1:1 solution without any pressure-related constraints.

The solvent-free MS-CG model performed well in simulations of lamellar structures. A large bilayer at 308 K showed good elastic properties and local stress when simulated using the MS-CG model matched to the bilayer structure. The simulated bilayer lacked long-range order and exhibited accelerated lateral diffusion compared to the reference atomistic MD bilayer, indicating a fluid state. Its structure suggested the presence of hydrogen-bonded complexes, often with a vacancy-like mechanism for Chol diffusion.

Various applications of the computationally efficient solvent-free MS-CG methodology for lipid bilayers are currently being developed in our group.

Acknowledgments

This research was supported by the National Institutes of Health (R01-GM063796). Computational resources were provided by the National Science Foundation through TeraGrid computing resources, specifically San Diego Supercomputer Center.

References

1. Gennis RB. Biomembranes: Molecular Structure and Function. Springer; New York: 1989.
2. Katsaras J, Gutberlet T. Lipid Bilayers: Structure and Interactions. Springer; Berlin: 2001.
3. Feller SE, Pastor RW. Biophys J. 1996;71:1350. [PubMed]
4. Tieleman DP, Marrink SJ, Berendsen HJC. Biophys J. 1997;1331:235. [PubMed]
5. Feller SE, Yin D, Pastor RW, MacKerell JAD. Biophys J. 1997;73:2269. [PubMed]
6. Feller SE, Venable RM, Pastor RW. Langmuir. 1997;13:6555.
7. Tu K, Klein ML, Tobias DJ. Biophys J. 1998;75:2147. [PubMed]
8. Smondyrev AM, Berkowitz ML. Biophys J. 1999;77:2075. [PubMed]
9. Smondyrev AM, Berkowitz ML. J Chem Phys. 1999;111:9864.
10. Feller SE, Pastor RW. J Chem Phys. 1999;111:1281.
11. Pasenkiewicz-Gierula M, Róg T, Kitamura K, Kusumi A. Biophys J. 2000;78:1376. [PubMed]
12. Smondyrev AM, Berkowitz ML. J Comp Chem. 1999;20:531.
13. Feller SE. Curr Opin Colloid Interface Sci. 2002;5:217.
14. Brannigan G, Lin LCL, Frank LH. Eur Biophys J. 2006;35:104. [PubMed]
15. Ayton GS, Noid WG, Voth GA. Curr Opin Struct Biol. 2007;17:192. [PubMed]
16. Tozzini V. Curr Opin Struct Biol. 2005;15:144. [PubMed]
17. Voth GA, editor. Coarse-graining of condensed phase and biomolecular systems. Vol. 455 CRC Press/Taylor and Francis Group; Boca Raton: 2009.
18. Nèmethy G, Scheraga HA. J Phys Chem. 1962;66:1773.
19. Shelley JC, Shelley MY, Reeder RC, Bandyopadhyay S, Klein ML. J Phys Chem B. 2001;105:4464.
20. Marrink SJ, Mark AE. J Am Chem Soc. 2003;125:15233. [PubMed]
21. Marrink SJ, de Vries AH, Mark AE. J Phys Chem B. 2004;108:750.
22. Brannigan G, Brown FLH. J Chem Phys. 2004;120:1059. [PubMed]
23. Brannigan G, Philips PF, Brown FLH. Phys Rev E. 2005;72:011915. [PubMed]
24. Soper AK. Chem Phys. 1996;202:295.
25. Lyubartsev AP, Laaksonen A. Phys Rev E. 1995;52:3730. [PubMed]
26. Garde S, Ashbaugh HS. J Chem Phys. 2001;115:977.
27. Murtola T, Falck E, Patra M, Karttunen M, Vattulainena I. J Chem Phys. 2004;121:9156. [PubMed]
28. Drouffe JM, Maggs AC, Leibler S. Science. 1991;254:1353. [PubMed]
29. Noguchi H, Takasu M. Phys Rev E. 2001;64:41913.
30. Cooke IR, Deserno M. J Chem Phys. 2005;123:224710. [PubMed]
31. Noguchi H, Takasu M. Phys Rev E. 2002;65:51907. [PubMed]
32. Farago O. J Chem Phys. 2003;119:596.
33. Cooke IR, Kremer K, Deserno M. Phys Rev E. 2005;72:011506. [PubMed]
34. Lyubartsev AP. Eur Biophys J. 2005;35:53. [PubMed]
35. Müller-Plathe F, Florian M. ChemPhysChem. 2002;3:754.
36. Reith D, Putz M, Müller-Plathe F. J Comput Chem. 2003;24:1624. [PubMed]
37. Sun Q, Faller R. Macromolecules. 2006;39:812.
38. Sun Q, Faller R. J Chem Theory Comput. 2006;2:607.
39. Chen LJ, Qian HJ, Lu ZY, Li ZS, Sun CC. J Phys Chem B. 2006;110:24093. [PubMed]
40. Izvekov S, Voth GA. J Phys Chem B. 2005;109:2469. [PubMed]
41. Izvekov S, Voth GA. J Chem Phys. 2005;123:134105. [PubMed]
42. Izvekov S, Voth GA. J Chem Theory Comput. 2006;2:637.
43. Chu JW, Ayton GS, Izvekov S, Voth GA. Mol Phys. 2007;105:167.
44. Zhou J, Thorpe IF, Izvekov S, Voth GA. Biophys J. 2007;92:4289. [PubMed]
45. Liu P, Izvekov S, Voth GA. J Phys Chem B. 2007;111:11566. [PubMed]
46. Noid WG, Chu JW, Ayton GS, Voth GA. J Phys Chem B. 2007;111:4116. [PMC free article] [PubMed]
47. Noid WG, Chu JW, Ayton GS, Krishna V, Izvekov S, Voth GA, Das A, Andersen HC. J Chem Phys. 2008;128:244114. [PubMed]
48. Noid WG, Liu P, Wang Y, Chu JW, Ayton GS, Izvekov S, Andersen HC, Voth GA. J Chem Phys. 2008;128:244115. [PubMed]
49. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J Chem Phys. 1983;79:926.
50. Darden T, York D, Pedersen L. J Chem Phys. 1993;98:10089.
51. Ayton G, Smondyrev AM, Bardenhagen SG, McMurtry P, Voth GA. Biophys J. 2002;83:1026. [PubMed]
52. Melchionna S, Ciccotti G, Holian BL. Mol Phys. 1993;78:533.
53. Leonard A, Escrive C, Laguerre M, Pebay-Peyroula E, Neri W, Pott T, Katsaras J, Dufourc EJ. Langmuir. 2001;17:2019.
54. Guo W, Kurze V, Huber T, Afdhal NH, Beyer K, Hamilton JA. Biophys J. 2002;83:1465. [PubMed]
55. Hofsäβ C, Lindahl E, Edholm O. Biophys J. 2003;84:2192. [PubMed]
56. Pandit SA, Bostick D, Berkowitz ML. Biophys J. 2004;86:1345. [PubMed]
57. Somerharju P, Virtanen JA, Cheng KH. Biochim Biophys Acta. 1999;1440:32. [PubMed]
58. Filippov A, Oradd G, Lindblom G. Biophys J. 2003;84:3079. [PubMed]
59. Scherfeld D, Kahya N, Schwille P. Biophys J. 2003;85:3758. [PubMed]
60. Almeida PFF, Vaz WLC, Thompson TE. Biophys J. 2005;88:4434. [PubMed]
61. Chiu SW, Jakobsson E, Mashl RJ, Scott HL. Biophys J. 2002;83:1842. [PubMed]
62. Karmakar S, Raghunathan VA. Phys Rev E. 2005;71:061924. [PubMed]
63. Vanderkooi G. Biophys J. 1994;66:1457. [PubMed]
64. Chong PL. Proc Natl Acad Sci. 1994;91:10069. [PubMed]
65. Sengupta P, Singh RRP, Cox DL, Slepoy A. Phys Rev E. 2004;70:021902. [PubMed]
66. Lindahl E, Edholm O. J Chem Phys. 2000;113:3882.
67. Gullingsrud J, Schulten K. Biophys J. 2004;86:3496. [PubMed]
68. Patra M. Eur Biophys J With Biophys Lett. 2005;35:79. [PubMed]
69. Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, deVries AH. J Phys Chem B. 2007;111:7812. [PubMed]
70. Szleifer I, Kramer D, Ben-Shaul A, Gelbart WM, Safran SA. J Chem Phys. 1990;92:6800.
71. Marsh D. Biophys J. 2007;93:3884. [PubMed]
72. Evans E, Needham D. J Phys Chem. 1987;91:4219.
73. Meleard P, Gerbeaud C, Pott T, Fernandez-Puente L, Bivas I, Mitov MD, Dufourcq J, Bothorel P. Biophys J. 1997;72:2616. [PubMed]