Solvent-Free, Highly Coarse-Grained Models for Charged Lipid Systems.

We present a methodology to develop coarse-grained lipid models such that electrostatic interactions between the coarse-grained sites can be derived accurately from an all-atom molecular dynamics trajectory and expressed as an effective pairwise electrostatic potential with appropriate screening functions. The reference nonbonded forces from the all-atom trajectory are decomposed into separate electrostatic and van der Waals (vdW) forces, based on the multiscale coarse-graining method. The coarse-grained electrostatic potential is assumed to be a general function of unknown variables and the final site-site interactions are obtained variationally, where the residual of the electrostatic forces from the assumed field is minimized. The resulting electrostatic interactions are fitted to screened electrostatics functions, with a special treatment for distance-dependent dielectrics and screened dipole-dipole interactions. The vdW interactions are derived separately. The resulting charged hybrid coarse-graining method is applied to various solvent-free three-site models of anionic lipid systems.

Introduction

The lipid membrane provides remarkable structural stability and integrity for the cell. It acts as a semipermeable boundary that selectively controls the diffusion of ions, proteins, and other molecules into and out of the cell, based on the biological requirements.[1,2] Cell membranes also play a very active role in many biochemical processes, which involve a complicated interplay between the membrane and membrane-bound proteins.[3] In a eukaryotic cell, almost one-fifth of the lipids in the cytoplasmic membrane leaflet are anionic.[4−6] The net negative charge drives the recruitment of proteins to the membrane surface through long-range nonspecific electrostatic interactions. This fact has been demonstrated for proteins such as GRP1 and AKT1, which carry the positive-residue-rich membrane-targeting Pleckstrin Homology (PH) domain motif.[7−11] The binding affinity to the membrane for the PH domain is found to decrease substantially when no anionic lipids are present in the bilayer.[12] Multivalent anionic lipids such as phosphatidylinositol (4,5)-biphosphate (PIP2) and phosphatidylinositol (3,4,5)-triphosphate (PIP3) play special roles in signaling pathways and are known to act as membrane markers.[10,13] Peripheral proteins bind to PIP lipids with high specificity due to the specific conformation and stereochemistry of charged moieties of lipid molecules. Some examples where PIP lipids play an important role in the early stages of membrane association are the Gag proteins in the HIV-1 virus,[14] epsin proteins in membrane remodeling processes,[15] and the GRP1 protein PH-domain[12,16] in certain membrane-docking-related chemotaxis processes. In addition, net negative charge is key in binding of membrane-remodeling proteins that do not have affinity to specific lipids, such as BAR domain proteins.[17−20] Therefore, the relative contribution from the hydrophobic and the electrostatic interactions must be taken into account to fully describe the processes of membrane sensing and remodeling by proteins.

Changes in lipid composition are sometimes known to induce dramatic effects on biological function,[21,22] a mechanism poorly understood at the molecular scale. In particular, anionic lipids specifically interact with positively charged protein residues, such as arginine, lysine, and histidine. The heterogeneities and topological diversities characteristic to real membranes has limited our understanding of them because they are difficult to study in vitro and detailed quantitative information is difficult to obtain in vivo.[4,23] Computational and theoretical modeling[24−26] studies can help to bridge these gaps, especially with the development of more-accurate and transferable potential energy functions.

The relevant length scales and time scales in biological systems span many orders of magnitude and so multiscale computational methods are required to model these systems. Molecular dynamics (MD) simulations at atomistic resolution are generally carried out for systems that have computationally tractable length and time scales. However, when studying large-scale biological problems, coarse-grained (CG) MD simulations can be used to explore much larger length and time scales. In a CG model, an explicit atomistic system is replaced by a representative reduced resolution system with fewer degrees of freedom, which evolve on a CG effective potential. Several different techniques for coarse-graining molecular systems are available in the literature.[27−31] The level of chemical specificity and the resolution of the model may vary based on the problems that are addressed. With decreasing resolution of coarse-graining, the chemical specificity and molecular scale properties such as local density distribution and pair correlation function are increasingly difficult to reproduce. Developing a low-resolution CG model that reproduces the electrostatics of the finer system is equally challenging, since electrostatic interactions are averaged over and the CG model may not quantitatively capture the interactions due to the charge distribution. Recently, many lipid models have been proposed with varying degrees of complexity, which include the explicit representation of electrostatic interactions at the CG resolution.[32,33] The recently developed ELBA force field[32] for membranes is one where the lipid sites have point charges and dipoles, and the water molecules are represented explicitly as a point dipole. The MARTINI force field parameters for glycolipids including phosphatidylinositol (PI) and its phosphorylated forms (PIP and PIP2) were also recently published.[33] In the latter, clusters of four water solvent molecules are explicitly represented as CG particles while the electrostatic interactions are represented using the Columbic potential energy function with a constant dielectric screening.

In this work, we provide a systematic method to develop solvent-free three-site highly CG models of anionic phospholipids with an accurate representation of the effective electrostatic interactions. Formulating such highly CG solvent-free models is nontrivial, because the net electrostatic interactions at the CG level are not obvious and, moreover, they cannot be easily derived from elementary considerations. Therefore, the present CG modeling is based on the multiscale coarse graining (MS-CG) method,[34−37] where electrostatic forces from all-atom (AA) trajectories are used to derive effective charge–charge, charge–dipole, and dipole–dipole interactions in the solvent-free reduced representation. The CG electrostatic interactions are fitted to screening functions while the CG van der Waals (vdW) interactions are derived from the residual forces by subtracting the effective fitted electrostatic interactions from the total force. Following the hybrid coarse-graining (HCG) scheme, developed in our previous work,[38] analytical functions are used at short-range where the atomistic MD sampling is limited, although, in this work, we use more-flexible analytical functional forms. These models for charged lipids are called the charged hybrid coarse-grained (cHCG) models. The present method is applied to develop a mixture of 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and monovalent anionic 1,2-dioleoyl-sn-glycero-3-phospho-l-serine (DOPS) in 1:1 and 3:1 molar ratios. In addition, a CG model of the lipid PIP3 is developed as an example. The method is general enough to be extended to other compositions and resolutions in a straightforward way, and it is validated with a higher resolution two-site CG model for methanol. All of the data pertaining to the methanol model are discussed in the Supporting Information (SI).

The remainder of this article is organized as follows. The theoretical framework for determining the pairwise effective electrostatic and vdW interactions for the cHCG systems is presented in section 2, and then the steps involved in systematically developing the cHCG models are then enumerated, followed by the details of the reference AA MD simulations and the force-matching process. The application of the method to large lipid bilayer sheets, vesicles, and tubules containing anionic lipids is next discussed in section 3. The cHCG model is then used in a membrane-protein example, where the interaction of PH domain protein with a membrane is studied, focusing on the PIP3 lipid. Finally, section 4 presents conclusions and suggests future directions.

The Charged Hybrid Coarse-Grained (cHCG) Method

A new MS-CG based scheme[34−41] is proposed in this work to obtain accurate effective electrostatic interactions for the chosen CG representation. The new cHCG model has separate nonbonded interaction parameters for vdW interactions and electrostatic interactions, making the model more transferable. For all systems, both the pairwise vdW interactions and the effective electrostatic interactions for the cHCG systems are determined from the AA simulation using the force-matching scheme. The bonded interactions are obtained using the Boltzmann inversion method, which is described elsewhere.[38,41] The reference electrostatic force (Fcharge) is due to charge–charge interactions at the AA resolution. The cHCG effective electrostatic interactions are obtained from these reference forces using the MS-CG procedure and fitted to an appropriate screening function. The difference between the total reference force and the electrostatic CG force field is then used to derive the cHCG vdW interactions. The flowchart in Figure 1 provides a summary of the cHCG method, and details of the parametrization process and validation are discussed in the section as follows. The CG representations of lipids are discussed in section 2.1, while the process of obtaining the effective electrostatic interactions from the reference AA trajectory and the procedure of fitting these interactions to a screened electrostatics functional form is discussed in sections 2.2 and 2.3, respectively. The parametrization of the vdW interactions with the corrections at short range are discussed in section 2.4, and the full algorithm to construct the cHCG model is given in section 2.5. Details of the atomistic simulations and force matching are provided in sections 2.6 and 2.7, respectively.

Figure 1

Flowchart showing the steps involved in the charged hybrid coarse-grained (cHCG) modeling protocol.

CG Representation

The DOPC, DOPS, and PIP3 lipids are modeled with three highly CG sites (Figure 2). The composition of each lipid CG site is listed in Table 1, whereas the mass and the headgroup dipole moments are listed in Table 2. As shown in Figure 2, the three sites of the CG lipid are called “head”, “mid”, and “tail”, respectively. The mid and tail sites for both DOPC and DOPS lipids have the same composition and are indistinguishable in the cHCG model. The head sites for DOPC, DOPS, and PIP3 lipids are hereafter called “headPC”, “headPS”, and “headPIP”, respectively. The interactions involving water molecules and salt ions are implicit in the lipid CG site interactions.

Figure 2

Three-site cHCG models of (a) DOPC, (b) DOPS, and (c) PIP3 lipid molecules in this study.

Table 1

Atoms Represented by Each CG Modela

CG site	DOPC atoms	DOPS atoms	PIP3 atoms
head	1–44 (44)	1–37 (37)	1–40 (40)
mid	45–66, 92–113 (44)	38–59, 85–106 (44)	41–77, 106–120 (52)
tail	67–91, 114–138 (50)	60–84, 107–131 (50)	78–105, 121–152 (60)

The numbering scheme of DOPC and DOPC is based on ref (51), and that of PIP3 is taken from ref (52).

Table 2

Mass and Point Dipole Moment of CG Sites

CG Site	DOPC sites	DoPS sites	PIP3 sites
mass (amu)
head	196.160	182.137	518.053
mid	311.440	157.147	310.435
tail	278.524	140.270	291.543
dipole moment (debye)	24.26	21.43	53.46

In our cHCG representation of lipids, there are two bonds that connect the head groups to mid groups and the mid CG sites to tail CG sites. There is one angular potential energy function for the head–mid–tail angle. The effective bonded interactions for the DOPC and DOPS CG lipid systems are the same as those used in the HCG model developed previously.[38] New bonded interactions parameters were derived for PIP3 lipids. Since the head groups are different for all lipids, each set of pair types also has different vdW pair interactions. The electrostatic pairwise interactions are restricted to the head groups since the mid and tail sites have nominal charge–charge interactions. Besides the screened monopole electrostatics interactions, the CG lipid head groups also carry fixed effective dipoles. A special distant-angle dependent screening function for dipole–dipole interactions is used in this work. This screening function is obtained from the effective CG electrostatic interactions derived from the AA trajectory.

Effective Pairwise CG Electrostatic Interactions

No initial assumption is made for the functional form of the cHCG pairwise electrostatic interaction and it is assumed to be a general function of unknown variables. Using the MS-CG method, the forces from the CG force field are minimized, with respect to reference electrostatic forces obtained from the AA MD trajectory.[34−37] In short, the minimization of the variational residual χ leads to the effective cHCG electrostatic interaction and is given bywhere f(r) are the reference AA forces due to charge interactions projected on the CG sites, calculated from the atomistic MD simulations. A “pseudo” simulation is carried out with the available trajectory to extract the electrostatic component of the forces from the AA trajectories, and a “re-run” of the system with just the electrostatic force field outputs the electrostatic forces for the given coordinates of the existing trajectory. The FIelec-MS(R;ϕelec)) are the CG forces calculated based on the assumed general functions. Here, n is the number of atoms in the reference AA system and N is the number of CG sites. The vectors r = {r1, r2, ..., r} and R = {R1, R2, ..., R} denote Cartesian spatial coordinates for the AA and CG representations, respectively. The CG force field (called a basis set[36]) contains Nd parameters, ϕ = {ϕ1, ϕ2, ..., ϕ}, which are determined using variational minimization of the residual shown above.

The effective electrostatic interaction potential, Uelec (R), for a CG pair is obtained by numerically integrating the CG force results as given byThe cHCG-derived effective electrostatic interactions between CG sites between the headgroup sites in DOPC:DOPS bilayer system are shown in Figure 3a. We find that the pairwise effective electrostatic interactions exhibit screening effects and do not fit the standard Coulomb interaction functional forms. Indeed, the effective electrostatic interactions for the highly CG solvent-free lipid head-groups are much more complicated (as compared to two-site methanol model in the SI), because of the level of coarse-graining and the implicit inclusion of solvent molecules and counterions. In fact, it would be extremely difficult or even impossible to guess the behavior of these highly CG effective electrostatic interactions. For this reason, they are treated separately.

Figure 3

Pairwise effective (electrostatic and vdW) interactions using the cHCG method between the headgroup sites in a 1:1 DOPC:DOPS system: (a) PC–PC effective electrostatic interactions, (b) PC–PS effective electrostatic interactions, (c) PS–PS effective electrostatic interactions, and (d) nonbonded vdW interactions between the head groups sites in 1:1 DOPC:DOPS system (error bar is on the order of 0.1–0.2 kcal/mol).

Fitting the Effective cHCG Electrostatics Interactions

For the two-site CG model of methanol, which is discussed in the SI, a simple Yukawa screening function that is given byis sufficient to fit the derived CG electrostatic force field between two CG sites I and J. In the equation above, A is the magnitude of the scaling constant and κ is the screening length in units of reciprocal distance. The Yukawa potential behaves like a Coulomb potential with exponential screening, converging to zero at a shorter distance. In the case of the CG model for methanol, the parameters of the Yukawa potential are determined by fitting with the AA-derived effective electrostatic interactions. The figure with the fitting for methanol CG model and the parameters are given in the SI.

However, we find that the above simplifications cannot be made for the much-lower-resolution solvent-free three-site cHCG models of monovalent and multivalent anionic lipid molecules. As shown in Figure 3, the effective electrostatic interactions between headgroup CG sites for the DOPC:DOPS lipid do not follow a simple exponential decay but have medium-range repulsion and strong effective attraction at short range. This effective electrostatic interaction can be understood in the context of colloid chemistry[42,43] and from the theory behind the range and shielding of dipole–dipole interactions in phospholipid bilayers.[26] The long-range repulsion acts as a stabilizing factor against the strong short-range attraction. Although the scales of these interactions are very low, the nonmonotonic nature of the effective electrostatic interaction indicates significant dipole effects and anisotropic electrostatic interactions. Instead of fitting the effective electrostatic interaction to empirical forms, we have used standard electrostatic functions (screened-Yukawa and dipole–dipole interactions with smoothening functions) to fit the effective electrostatic profile between CG sites I and J. The fitting is shown by the red curves in Figure 3.

The functional form used to fit the effective electrostatic interactions is given below asHere, A0 and κcg are the magnitude of scaling constant and screening lengths for the Yukawa type monopole electrostatic interaction between the lipid headgroup sites and RC is the cutoff length (20 Å). The charge-dipole interaction term is not included in eq 4, because it would amount to overfitting. This is because the charge-dipole interaction decays with 1/r3 and it is considered implicit in the first term of dipole–dipole fitting that also decays by 1/r3. We note that the parameters resulting from the fitting procedure may not always identify the various physical components of the electrostatic interactions in the CG space, because of the very low resolution of the cHCG lipid model. By contract, we are able to extract the physically relevant parameters when the method is applied to the high-resolution 2-site methanol model (see the SI). The blurring or “abstraction” of physical information for the 3-site solvent-free cHCG lipid model is due to the aggressive nature of coarse-graining, wherein we are trying to capture the complex electrostatic interactions within a few CG sites. Therefore, our focus is to describe the effective interaction (Figure 3) with a physically motivated ansatz, since it provides a much better framework to fit the data than simply fitting it with ad-hoc parameters.

During the fitting process, the initial values for A0 and κcg are taken as −1.0 kcal mol–1 Å–1 and 0.15 Å–1, respectively. The final values are shown in Table 3 for all the headgroup pairs. The vectors μ̂ and μ̂ are the effective dipole vector on each CG head site and is given a magnitude of 1. Unit dipole vectors are used for computational convenience. The actual magnitude of the dipoles are determined from coefficients A1 and A2. The initial magnitude of the dipole for head groups sites belonging to DOPC, DOPS, and PIP3 are obtained by calculating the average dipole moment from the AA trajectory using bare charges in the AA representation and is found to be 21.43, 24.26, and 53.46 D, respectively. The initial values of A1 and A2 are estimated as products of the dipole moments for a given pair. The final fitted values are shown in Table 3. The normal and in-plane dielectric screening parameters for the dipole–dipole interactions are denoted as ϵ⊥ and ϵ∥, respectively. The relative dielectric constant of water is ∼80, and that of the interior of the membrane is between 2 and 4.[44] However, to our knowledge, the relative dielectric permittivity at the interface (headgroup region) is not known experimentally.[25,26,45] Also, theoretical studies have shown that the normal component of dipole–dipole interaction is repulsive, whereas the in-plane dipole–dipole interaction is attractive, and both of them are screened by very different dielectric permittivities. The repulsive part dominates at larger lengths and the in-plane attractive part dominates at the shorter distances. To fit the effective electrostatic interactions, the initial value of ϵ⊥ and ϵ∥ is assumed to be 10 and 80, respectively. The final optimized parameters in eq 4 for different pairs are listed in Table 3.

Table 3

Optimized Parameters for the Effective Electrostatic Interaction in eq 4

	headPC–headPC	headPC–headPS	headPS–headPS
A0 (kcal mol–1 Å–1)	–3.555	–3.467	–4.258
κcg (Å–1)	0.1198	0.1282	0.1485
A1	588.214	519.892	459.014
A2	346388.27	270287.48	210905.87
ϵ⊥	8.04	9.61	8.78
ϵ∥	76.81	78.23	77.36

Anisotropy in the dipole moment is also introduced by biasing the dipole orientation, using a loose harmonic restraint on the angle that the dipole vector makes with the head–tail vector of the lipid. This angle is denoted with the symbol ω below. The biasing potential[32] is given bywhere ω is calculated as the angle between the dipole vector and the vector R⃗ht that connects the head CG site of the lipid with the tail CG site:Kω is the rigidity constant. A loose rigidity constant of 5 kcal/(mol degree2) is used for all the head groups and the mean orientations (ω0) for headgroup belonging to DOPC, DOPS, and PIP3 are given the values of 17°, 25°, and 40°, respectively. We used the distribution of the angle between the vector joining the phosphorus and nitrogen in the phospholipid headgroup, with respect to outward normal of the membrane to calculate ω0 for the PC and PS head groups. However, in the absence of sufficient sampling for the PIP lipids, we estimated the angle distribution data from the simulation work done by others on PIP2 and PIP3 lipids in bilayer.[46]

The vdW Interaction Parameters of cHCG Systems

The vdW interaction parameters for the cHCG systems are derived from the AA trajectory using the nonelectrostatic reference forces.[38] MS-CG procedure is applied on these forces to obtain the cHCG vdW pairwise interactions. The vdW interactions parameters obtained from the MS-CG process complete the derivation of the nonbonded interactions. However, the MS-CG procedure is unable to provide a CG potential energy function of the short-range interactions, because of statistically insufficient sampling at short ranges. This limitation is discussed in detail in the original HCG paper.[38] Therefore, analytical functions are used to account for poorly sampled regions of the configurational space. The analytical function selected for the short-range interaction in the current lipid models is of the formwhere ε, σ, m, n, and RSR are the five parameters that define the analytical short-range interaction. The initial values of ε, σ, m, n, and RSR are given as 1 kcal/mol, 6.0 Å, 12, 6, and 8.0 Å, respectively. The parameters m and n are restricted to integer values. Final values of the first four parameters are obtained iteratively through the fitting process. User adjustment is made to the fifth parameter, RSR, which determines the range of analytical function form for the vdW interaction. This is because the reference AA trajectory is unable to fully sample the interactions at short range for the highly CG representation. (If the sampling was more extensive, no adjustment would be needed.) The final fitted values are incorporated in tabular form, and all the tables for the vdW interactions between all pairs are provided in the SI. Although the analytical functional form for the short-range interaction in the cHCG model is different than the one used in the earlier HCG model[38] (where a simple 12–6 LJ interaction was used), the fitting procedure is same as that described in the original paper.[38] The new analytical function form was chosen because it provides greater flexibility in the fitting of the short-range vdW interaction.

The CG force field for bonded interactions is obtained using the Inverse Boltzmann method.[47,48] The bonded and angular interactions of the DOPC and DOPS are the same as in the HCG model.[38] The bond and angle distributions for the PIP3 lipids are given in the SI.

Building the Solvent-Free cHCG Model

The protocol involved in developing the cHCG model for molecular systems are enumerated below:

(1) Generate a short AA reference trajectory for the chosen system.

(2) Choose the CG resolution and the composition of each CG site. The reference CG coordinates and forces for the given representations are mapped from the corresponding AA trajectories. Perform additional calculations to store the forces due to electrostatic interactions from the AA simulations. Separate modules in LAMMPS[49] MD package were written for this purpose. This force is used to obtain the effective electrostatic interactions in the CG systems.

(3) Use the Inverse Boltzmann method to obtain the bond and angle-bending interactions of the given systems, using the probability distributions from the AA trajectory.

(4) Use the reference coordinates and forces to obtain the nonbonded pairwise interactions for the CG systems as described in detail above. Broadly, this is carried out in three steps:

Using the force-matching optimization algorithm, obtain the effective electrostatic interactions due to reference electrostatic forces.

Fit the effective electrostatic interactions to screened electrostatics interaction functions.

Calculate the vdW interaction using force-matching procedure on reference vdW force, which are obtained by subtracting the forces due to fitted screened electrostatics from the total reference forces calculated in step 2. For highly CG lipid models with insufficient short-range sampling, supplement the vdW interaction with an analytical function, as discussed in section 2.4.

Reference AA MD Simulations

Three different AA reference bilayer systems were simulated. System 1 consisted of a bilayer patch composed of 36 DOPC and 36 DOPS lipid molecules. System 2 contained 54 DOPC and 18 DOPS lipids. In both lipid systems, the monolayers had an equal composition of lipids. In system 3, one DOPS lipid in each leaflet from system 2 was replaced by PIP3. Initial configurations of the mixed DOPC–DOPS bilayer systems were generated with the CHARMM-GUI.[50] The CHARMM36 force-field was used for the DOPC and DOPS lipids.[51] The parameters for PIP3 were used from Lupyan and co-workers.[52] The lipid system was then solvated with 2100 TIP3P[53] water molecules and ionized with Na+ and Cl– ions with a salt concentration of 150 mM. A cutoff of 12 Å was used for the nonbonded vdW forces with a LJ switching function from 8 Å to 12 Å, without long-range corrections. Particle mesh Ewald summation[54] was used to calculate long-range electrostatic interactions. All C–H bonds were constrained using the SHAKE algorithm.[55] An integration time step of 2 fs was used. After the initial energy minimization, the fully solvated mixed lipid bilayer system was simulated for 40 ns under constant NPT conditions using the Nosé–Hoover thermostat at 310 K and a Langevin piston[56] at 1.01 bar.

The bilayer were validated against reference data[51] and the average area per lipid (A) in the last 10 ns of the NPT run was found to be 71.82 Å2 and 71.91 Å2 for 1:1 and 3:1 (DOPC:DOPS) systems, respectively. Following the NPT simulations, all the systems were run in the constant NVT ensemble for 10 ns. The last 4 ns were used to sample reference forces and coordinates for the CG model development. System minimization, equilibration, and dynamics were performed using the LAMMPS[49] molecular dynamics software package. Visualization was performed with VMD.[57] Analyses were performed using LAMMPS and GROMACS.[58]

Details of the Force Matching

For the chosen resolution, the reference coordinates (CG centroid positions) and the average reference forces for each configuration were obtained from the AA trajectory. The forces due to water and ions were included in the effective force calculation of the lipid CG sites. The MS-CG least-squares optimization was carried out using the normal matrix equations and B-Spline basis functions were used to interpolate the MS-CG forces. Details about the B-Spline-based parametrization can be obtained in refs (41 and 59). The final CG force–displacement results for the electrostatic interactions and the vdW interaction were integrated to obtain the cHCG potential energy functions, as discussed above.

Results and Discussion

Here, we discuss the results from the cHCG models for DOPC–DOPS and DOPC–DOPS–PIP3 lipid systems. Finally, we apply the electrostatic bilayer model on a membrane-protein system and show that the model can be used for further studies to understand the role of anionic lipids in recruiting proteins to the membrane surface.

The structural properties of lipid bilayers, such as area per lipid, z-density profile, bilayer thickness, area compressibility modulus, bending modulus, and lateral pressure profile, are calculated and compared to available reference data. To show the scope and fidelity of the model for long length scales and time scales, cHCG simulation results from large-scale simulations on tubules and vesicles are also shown. Simulations of large patches with a 3:1 DOPC:DOPS bilayer are also carried out and compared with the 1:1 DOPC:DOPS bilayer.

DOPC–DOPS Bilayer

Three sets of system sizes were considered for the 1:1 DOPC:DOPS CG flat bilayer simulations. The three systems have 72, 1250, and 5000 lipids with equal number of lipids in each leaflet. The smallest system size of 72 lipids is equivalent to the AA system used to derive the cHCG interaction parameters. The intermediate-sized 1250-lipid system (∼100 Å × ∼100 Å) and a large patch with 5000 lipids (∼400 Å × ∼400 Å) are used to show the transferability of the model for large systems. Another CG simulation with a 5000-lipid DOPC:DOPS bilayer system with 3:1 composition is also carried out and the trajectory is used to analyze the composition-dependent property changes in bilayer systems. The membrane simulations were run in the constant NPT ensemble for the flat sheets. The z-dimension was fixed and the membrane was allowed to breathe in the xy dimensions to simulate fixed surface tension boundary conditions. The 5000 CG lipid bilayer sheet was simulated for 40 ns. The area per lipid (AL), bilayer thickness (dHH), z-density profile, area compressibility modulus (KA), lateral pressure profile, and bending modulus (κc) for each system were calculated and compared with available reference data. The structural parameters AL, dHH, KA, and κc are reported in Table 4.

Table 4

Area Per Lipid, Membrane Thickness, Area Compressibility Modulus, and Bending Modulus of cHCG Models of DOPC:DOPS Patches with (a) Different Sizes of the 1:1 System and (b) Different Compositions with 5000 Lipids

DOPC:DOPS	AL (Å2)	ZHH (nm)	KA (dyn/cm)	κc (kBT)
(a) Parameters Determined Using Different Sizes of the 1:1 System
72	70.92	47.3	146.01	18.23
1024	70.63	48.1	132.12	17.89
5000	70.48	48.7	97.09	15.54
refs (61, 67, and 76)	71.82	40–50	150–300	12–25
(b) Parameters Determined Using Different Compositions with 5000 Lipids
1:1	70.48	48.7	97.09	15.54
3:1	71.02	49.1	121.46	15.87

z-Density and Lateral Distribution

The z-density profile is in good agreement with the corresponding distribution from the AA trajectory (Figure 4a) The bilayer thickness is taken as the average distance between the head groups in each leaflet and was found to be ∼50 Å. Moreover, the overlaps in the distribution of different CG sites indicate that the bilayer is fluid, which is further validated by mean-squared displacement (MSD) and lateral pressure profile calculations.

Figure 4

(a) z-density profile for 1:1 DOPC:DOPS three CG sites from the atomistic simulation (black) and the CG simulation (red). (b) The accumulated average number of lipids in one of the membrane leaflets with total 5000 lipids. The figure shows the lateral organization of the lipids for 1:1 DOPC:DOPS and 3:1 DOPS system. The y-axis denotes the number of lipids around a central lipid and the results indicate almost perfect mixing. The dashed black line overlaps the black line (complete mixing), and the blue and red lines (1:1 systems) also overlap, since PC and PS lipids have almost identical lateral organization.

The mixing of lipids is quantified with a nearest-neighbor analysis[60] by calculating the accumulated average number of lipids in one of the membrane leaflets.[41,59] The increase in the accumulated average of both lipids, as a function of radial distance (Figure 4b), is directly proportional to the area, which rules out any strong hexagonal correlations or domain formation. The nearest-neighbor analysis yields values of 0.504 and 0.496 for DOPC and DOPS lipids, respectively, in the 1:1 mixture, whereas it gives 0.753 and 0.247 in the 3:1 mixture system. The values obtained are proportional to their mole fractions, indicating a near-ideal mixing. Such behavior is expected for bilayers containing only DOPC and DOPS lipids.

Area per Lipid and Area Compressibility Modulus

The average Al for 1:1 and 3:1 DOPC:DOPS bilayer system are found to be 70.48 Å2 and 71.02 Å2, respectively, in good agreement with values obtained in AA simulations (71.82 Å2 and 71.91 Å2). For comparison, Al of pure DOPC system was found to be 67.4 Å2 in CG simulations.[38] The net decrease in Al can be attributed to the softer CG potential energy landscape. We did not see any perceptible change in Al in systems containing one PIP2 or PIP3 lipid (see Table 4).

The area compressibility modulus (KA) of the CG lipid system is estimated from the fluctuations in Al (see Figure 5). KA of the 5000-lipid system is found to be 97.09 dyn/cm (1:1 composition), comparable to recent micropipette-aspiration experiments[61] (67 dyn/cm) (Table 4b). Pure DOPC bilayer is reported to have much higher KA (188–256 dyn/cm).[62,63] Since the area per lipid of unary and binary systems does not appreciably change, such a remarkable change in area compressibility modulus is primarily due to the significant increase in membrane area fluctuations. These fluctuations are caused by the disparate chemical and steric nature of head groups between the two types of lipids. In addition, KA seems to decrease with the increased size of the simulated bilayer patch (Table 4).

Figure 5

Fluctuations in area per lipid with time in CG simulations with the 1:1 DOPC:DOPS cHCG model. The average area per lipid is shown as dashed lines, and the compressibility modulus is calculated using the average area per lipid and its fluctuations.

Lipid Diffusion

The two-dimensional MSD is calculated to estimate the lipid diffusion constant and thus evaluate the fluidity of the bilayer model. As expected, the MSD of both CG lipid systems is linear with time (see Figure 6) and indicates the long-time-scale liquid-like behavior in the lateral plane. The diffusion constant of lipids in both 1:1 and 3:1 DOPC:DOPS mixture is 2.5 × 10–7 cm2/s. The diffusion constant is an order of magnitude higher than experimental measurements (∼10–8 cm2/s), indicating that the model effectively simulates much larger time scales than predicted from AA MD. Such behavior is consistent with all CG models, considering that the reduced degrees of freedom, no friction terms in the dynamics,[64] and the softer free-energy profile facilitate faster dynamics. It seems that the diffusion constant does not depend on the relative ratios of the two lipids, confirming the expected high miscibility of the simulated composition.

Figure 6

MSD of DOPC:DOPS bilayer system with 5000 lipids calculated using center of mass diffusion. The MSD for 1:1 DOPC:DOPS is denoted in black, and the MSD for 3:1 DOPC:DOPS is shown in red.

Lateral Pressure Profile and Surface Tension

At any plane on the bilayer, lateral pressure is defined as the difference between the in-plane and normal component of the pressure tensor at that plane. For a bilayer with surface normal parallel to the z-axis, if the diagonal elements of the pressure tensor are given by P, P, and P, then the lateral pressure across the membrane thickness (P(z)), is equal to (1/2)(P + P) – P at a given membrane depth. The lateral pressure profile across the thickness of the bilayer is indicative of the fluid nature of the membrane and is related to the larger-scale properties, such as surface tension, spontaneous curvature, and bending modulus.[65−67] The details of calculating the P(z) can be found in an earlier work.[38] The P(z) of the 1:1 5000-lipid DOPC:DOPS system is shown in Figure 7a. The qualitative shape of the P(z) function is consistent with liquid-phase lipid systems. The pressure has a net positive value near the headgroup and is a consequence of effective repulsion between the head groups. The large drop in lateral pressure immediately below the headgroup is characteristic of the hydrophobic interactions. The positive pressure at the center of the molecule is thought to originate from the loss of conformational freedom of the tails in the center of the bilayer, which leads to an effective repulsion between the tail sites. A plot of the surface tension of the bilayer with time (Figure 7b) shows that the repulsive forces near the headgroup and the center of the bilayer balance the attractive interactions in the hydrophobic regime of the membrane. The surface tension (γ) of the system is given bywhere dHH is the membrane thickness. The convergence of surface tension to zero indicates that the bilayer has achieved an optimal packaging configuration and the free energy is minimized with respect to the area, which is an observation that is corroborated with an excellent convergence in area per lipid.

Figure 7

(a) Lateral pressure profile for the 1:1 DOPC:DOPS bilayer system with 5000 lipids. (b) Surface tension evolution with time for the same system.

Bending Modulus

Bending modulus is one of the most important material properties of the membrane and has a direct bearing on biological processes, such as curvature changes in cell membrane, membrane remodeling, cell division, and fusion.[17,68−70] One of the quickest ways to estimate the bending modulus of a lipid bilayer is by assuming that it is a thin elastic sheet and relating the bending modulus to the area compressibility modulus and membrane thickness. The relationship is given by κc = KAdHH2/α, where κc is the bending modulus. The dimensionless factor α is generally assumed to be 24.[68]

The κc value was found to be 18, 17, and 16 kBT for 1:1 bilayer systems with 72, 1250, and 5000 lipids, respectively, where kB is the Boltzmann’s constant and T is the temperature (in Kelvin). The size dependence of κc is attributed to increasing undulation with increasing patch size.[71−73] Also, because of undulations, the effective dHH is numerically estimated to be higher for bigger patches, since it is calculated across the entire system. The κc of 3:1 system with 5000 lipids is calculated to be 15 kBT. The bending modulus was also calculated from the long wavelength undulation spectrum using the Helfrich theory.[74] From this, κc was estimated to be 15.5 kBT and 15.9 kBT for 1:1 and 3:1 DOPC:DOPS bilayer with 5000 lipids. The 3:1 DOPC:DOPS bilayer is slightly stiffer than 1:1, which is expected since the bending modulus for the pure DOPC was estimated to be ∼18–19 kBT from the earlier HCG model. The values of κc from our cHCG model are also in the range of experimental observations. For similar lipid systems, the following values have been reported:[61,67,76] 13 kBT, from pipet aspiration;[62] 16 kBT, from X-ray data;[75] 30 kBT, from thermally excited shape fluctuations; and 32 kBT, from optical measurements.

Simulation of Vesicles and Tubules

Figure 8 shows snapshots and the time evolution of the trajectory of a simulation of a vesicle with high curvature (1/20 nm–1). The system contains 19 450 lipids and has an equal number of DOPC and DOPS lipids. The simulations for vesicles and tubules were carried out in the constant NVT ensemble with excess vacuum space around the lipid system. Time step of 20 fs was chosen and the simulation was run for 5 million time steps, which is equivalent to 100 ns of simulation time. During the course of the simulation, some pores initially form on the vesicle, compensating for the low area per lipid in the initial configuration. The long-time stability of small vesicles using the cHCG model demonstrates the fidelity of the cHCG model for high curvature systems.

Figure 8

Snapshots of the CG simulations on a vesicle with 1:1 DOPC:DOPS composition (the vesicle size is ∼40 nm diameter (∼20 000 lipids) and simulated for 100 ns): (a) initial system with a cross-section across the center; (b) the configuration after 5 ns; (c) cross-sectional view after ∼50 ns; and (d) final converged simulation after the length of the simulation. The test was performed to examine the resilience of the model to high-curvature systems.

Another large-scale simulation was carried out on a membrane tubule with a height of 50 nm (Figure 9a) and a diameter of 30 nm (Figure 9b). The tubule was composed of 39 125 cHCG lipids, with an equal number of DOPC and DOPS. In the absence of any membrane remodeling proteins or external stress, the tubular bilayer system should converge to the spherical system to minimize the surface tension. We observe that, after ∼150 ns of the simulation time, the tubule transforms to a vesicle. At the start of the simulation (Figure 9c), the two open ends of the cylinder start closing-in, followed by a transformation into an energetically more favorable spherical structure. Both the vesicle and tubule simulation were carried out on a 512-processor cluster computer system. Full convergence (100 ns runs) of the ∼20 000 cHCG lipid vesicle system was achieved within ∼10 h, while the transition from tubule to vesicle (∼180 ns runs) for an ∼40 000 lipid system was achieved within 24 h of simulation, emphasizing the significant acceleration that the cHCG model can achieve.

Figure 9

Snapshot of CG simulations on tubule with ∼40 000 lipids and simulated for 150 ns: (a) initial configuration showing the front view; (b) initial configuration showing the top view; (c) snapshot after 20 ns, showing the closing in of the open-end tops; and (d) fully converged sphere after 150 ns of simulation time.

Membrane–Protein CG System with Electrostatic Interactions

Since only one PIP3 lipid per monolayer is used with 1:1 DOPC-DOPS lipids, there is no remarkable effect on large-scale physical properties such as area per lipid and compressibility modulus, and they are not reported. However, the dipole orientations of the PIP3 lipid headgroup were monitored and compared with reference AA data as a part of the model validation. Since the current lipid cHCG model is developed with the objective of studying the role of anionic lipids (monovalent/multivalent) in membrane–protein interactions, a particular example of membrane-bound protein (PH domain)[12] is chosen to highlight the capability of the model to study such systems at large length and time scales.

The electron-rich inositol rings in the PIP3 has an orientation of ∼40°, with respect to the bilayer normal,[46] and the orientation is known to change upon association with proteins.[12,16] In this paper, we provide a simple example (schematics shown in Figure 10) where we study the response of PIP3 dipole moment vector as a function of PH domain approach to the membrane surface. A 22-site CG model of the PH domain was developed using the heterogeneous elastic network model (heteroENM) of proteins.[77,78] The AA trajectory for the parametrization was generated using the high-resolution structural data of GRP1 PH domain [Protein Data Bank ID: 1FGY]. The hetroENM parameters are provided in the SI. The effective dipole moment of each CG site on the protein was also calculated using VMD dipole package and superimposed on each CG site of the protein. They are also listed in the SI. The cHCG simulations were carried out under fixed NPT boundary conditions for 40 ns.

Figure 10

(a) PH-domain protein with CG sites in black and dipole moment vectors on CG sites shown as small arrows; the overall dipole orientation of the protein is shown with the long gray arrow. (b) PIP3 lipid, with the arrow depicting the headgroup dipole orientation (the headgroup orientation is defined as in ref (77)). (c) Snapshots of CG simulation of PH domain on DOPC:DOPS–PIP3 membrane. The protein CG sites are shown in gray; the DOPC, DOPS, and PIP3 headgroup sites on the top leaflet are shown in blue, red, and green, respectively. The orange arrow indicates the dipole vector direction on the PIP3 headgroup as defined in ref (75).

In our previous studies with the AA simulations of the same system,[12] it was shown that the PH domain has an overall dipole (facing the membrane surface), which helps in the electrostatic steering of the protein toward anionic surfaces. In this cHCG simulation, we see that the PIP3 lipid responds to the approach of the PH domain and aligns itself favorably to optimize the electrostatic interactions between the membrane and the protein. The dipole orientation of PIP3 in the bottom leaflet does not show any realignment with time (Figure 11a), and the distribution of PIP3 dipole vector is found to be normal and agrees relatively well with the observations made in ref (46). On the other hand, the dipole vector of the PIP3 lipid on the top surface, which faces the PH domain, aligns with the protein dipole and constantly realigns toward favorable interaction energy (Figure 11b). The above example was carried out to highlight the capability of the cHCG model and to show the promise that this model holds for exploring the role of anionic lipids in large length and time scale protein–membrane association processes.

Figure 11

(a) Distribution of the PIP3 lipid dipole angle in the native bilayer environment on the bottom leaflet. (b) Time evolution of the angle between the protein and PIP3 dipole vectors on the top leaflet.

Concluding Remarks

A new force-matching based charged hybrid coarse-grained (cHCG) method is developed with the ability to explicitly represent screened electrostatic interactions between the coarse-grained (CG) sites. The method is applied to develop very-low-resolution three-site, solvent-free models of anionic lipids. Since the electrostatic interactions at the CG level are derived from the AA trajectories, they provide a realistic picture of the charge-screening environment of different CG sites. We show that the different CG pairs have different screening lengths and nonlinear distant-dependent dielectric behavior is also observed in certain cases. When applied to develop solvent-free three-site model of lipids, the model had to be extended (as compared to methanol model) to incorporate the dipole effects and nonlinear nonmonotonic dielectric screening behavior. Also, the insufficient AA sampling at short range was compensated with an MS-CG informed analytical functional form, similar to that used in the original HCG model. There is always a possibility that the ion distributions in the bulk water may not be sampled properly due to the small number of ions present in the reference bilayer system, which may affect the derivation of the effective electrostatic interactions. Increasing the number of ions (while maintaining physiological salt concentrations) also results in intractable size of the reference system for parametrization purposes. To test the convergence with respect to the electrostatics sampling, we calculated the electrostatic membrane potential for the lipid bilayer with 40-ns and 50-ns trajectories (see Figure S6 in the SI) and our data compare well with the benchmark results.[51] The precise reproduction of the experimental potential drop for bilayer in simulation is still an open area of research (and somewhat controversial), because the simulated values are very sensitive to the definition of correct ion distribution and also sensitive to the distributions of partial charges on the lipid molecules. We found that, generally, 40–50 ns of AA simulation time for a lipid bilayer system with 72 lipids provides sufficient sampling for the derivation of effective electrostatic interactions from our reference trajectories.

Relevant biophysical membrane properties such as area per lipid, area compressibility modulus, bending modulus, and the lateral stress profile were found to be in agreement with results from the literature. The diffusion analysis confirmed the two-dimensional fluid nature of the bilayer. No domain separation was observed and the mixed DOPC:DOPS system exhibited ideal lateral mixing for all compositions. The cHCG simulation of PH-domain protein on DOPC:DOPS:PIP3 bilayer was used as an example to illustrate the ability of the model to study membrane-sensing and membrane-binding processes at long length and time scales. Mechanistic studies electrostatically driven recruitment of proteins by lipids could not be studied by the earlier HCG model. However, the new cHCG model can be used to carry out long-time-scale simulations of large membrane–protein systems (for example, the PIP–PH domain studied in this work) to better understand the role of electrostatics in lipid recognition by various membrane-binding proteins. The dipole orientation of the PIP3 headgroup, which interacts with the PH domain protein, was tracked over time and the dipole orientation on the lipid is shown to evolve with the approach of the protein on the surface. The lipid dipole finally aligns such that interaction between the protein and the lipid is most favorable for the binding. In the future, an expanded coarse-graining of individual proteins will be carried out to quantify the effective electrostatic interaction between the cHCG lipid model and the proteins at a highly CG level more accurately. The resulting CG protein–cHCG lipid model can then be used to perform long-time-scale simulations on large membrane–protein systems to understand the role of electrostatics in lipid recognition by various membrane-binding proteins.