|Home | About | Journals | Submit | Contact Us | Français|
Polarizable force fields for lipid and solvent environments are used for molecular dynamics simulations of a fully hydrated dipalmitoylphosphatidylcholine (DPPC) bilayer and gramicidin A (gA) dimer embedded in a dimyristoylphosphatidylcholine (DMPC) bilayer. The lipid bilayer is modelled using the CHARMM charge equilibration (CHEQ) polarizable force field for lipids and the TIP4P-FQ force field to represent solvent. For the DPPC bilayer system, results are compared to the same system simulated using the nonpolarizable CHARMM27r (C27r) force field and TIP3P water. Calculated atomic and electron density profiles, headgroup orientations as measured by the phosphorus-nitrogen vector orientation, and deuterium order parameters are found to be consistent with previous simulations and with experiment. The CHEQ model exhibits greater water penetration into the bilayer interior, as demonstrated by the potential of mean force calculated from the water density profile. This is a result of the variation of the water molecular dipole from 2.55 D in the bulk to 1.88 D in the interior. We discuss this finding in the context of previous studies (both simulation and experiment) that have investigated the extent of penetration of water into DPPC bilayers. We also discuss the effects of including explicit polarization on the water dipole moment variation as a function of distance from the bilayer. We show distributions of atomic charges over the course of the simulation, since the CHEQ model allows the charges to fluctuate. We have calculated the interfacial dipole potential, which the CHEQ model predicts to be 0.95 V compared to 0.86 V as predicted by the C27r model. We also discuss dielectric permittivity profiles and the differences arising between the two models. We obtain bulk values of 72.77 for the CHEQ model (TIP4P-FQ water) and 91.22 for C27r (TIP3P), and values approaching unity in the membrane interior. Finally, we present results of simulations of gA embedded in a DMPC bilayer using the CHEQ model and discuss structural properties.
Membranes and membrane-bound proteins are a vital part of biological systems. Indeed, membrane-bound proteins are involved in a number of physiological functions including but not limited to passive and active transport, signaling processes, and interfacial enzymatic processes1. Membrane-bound proteins have been found to make up approximately one third of the human genome2, a fact that further emphasizes their importance. Equally important is the lipid environment within which these proteins function. Recent studies have explored structual properties of membranes as well as electrostatic properties such as the dielectric variation within a bilayer3, the interfacial potential4, and the interactions of polar or charged amino acid side chains with hydrocarbon tails5. Experimental studies have also probed these properties in recent years. For example, structural properties of bilayers have been determined by x-ray and neutron scattering6,7 and NMR8, while the water penetration into the bilayer interior has been investigated via electron spin spectroscopy techniques9,10 as well as x-ray scattering11. However, even current state-of-the-art experimental measurements are not always able to provide the type of detailed atomic level resolution that would provide significant insights into the mechanisms of membrane systems. To this end, computational methods such as molecular dynamics and Monte Carlo simulations have been employed to study properties and processes in such systems at the atomic level1,5,12–23.
Molecular dynamics simulations require information about the atomic interactions within and between species. Quantum mechanical calculations would give the most accurate results; unfortunately, such computations are very time-consuming for all but the smallest systems. Because of this limitation, classical “force fields” have been developed which consist of potential energy functions chosen to model inter-species and intra-species interactions. In addition to models in which all atoms are modeled explicitly (i.e. all-atom models), there has been significant effort towards “coarse-grained” models24,25 in which entire small molecules or functional groups are represented as a single unit. Development and applications of such models continues today26–31. Implicit solvent and lipid models32–35, which represent solvent or lipid as a continuum, are also a viable alternative to all-atom models for simulating large biomolecular systems. However, all-atom simulations remain the standard to which such coarse-grained and implicit models are compared, so the need for more accurate atomistic models is clear.
Early, pioneering applications, along with today’s paradigmatic (state-of-the-art) all-atom simulation methods have employed fixed-charge representations of electrostatic interactions based on a Coulomb model. The shortcomings of such models, and the plausible importance of explicitly accounting for non-additive electrostatic and charge transfer effect, have been widely discussed in the literature36–39. The past decade has witnessed an increasing pace of development and application of polarizable force fields for a range of applications, though such models have not yet realized the popularity enjoyed by fixed-charge models. In order to begin to explore the effects of polarization in biological systems, the first step undoubtedly has to be the development of self-consistently parameterized models.
Polarizable interaction models that incorporate dipole induction effects have already proven an indispensable tool for obtaining an accurate theoretical estimation of solution structure and thermodynamics in interfacial systems such as aqueous solutions of inorganic salts. The development of non-additive, or polarizable, force fields for small molecular and larger biologically relevant macromolecules40–50 has attracted considerable interest and has resulted in the development of several conventional approaches for modeling atomic and molecular polarization. These approaches include point-dipole (and higher-order multipole) polarizable models41,51–53, Drude oscillator models44,54,55, and charge equilibration/fluctuating charge models42,43,56–61 in addition to fully ab initio based approaches such as Car-Parrinello molecular dynamics techniques.
The goals of this contribution are as follows. First, in Section II we discuss the charge equilibration (CHEQ) formalism, how it explicitly accounts for non-additive electrostatic effects, and our methodology for carrying out simulations using the CHEQ model. In Section III we discuss various properties of a model DPPC bilayer predicted using molecular dynamics simulations. In Section IV we present results of simulations of gramicidin A in a DMPC bilayer. Section V concludes our contribution and discusses continuing work in our laboratory.
The present simulations employ polarizable force fields based on the charge equilibration model for the entire lipid molecule and the solvent. The formalism employed to explicitly treat non-additive electrostatic effects is the charge equilibration model. In the following, we discuss the specifics of the formalism. The development of the CHEQ model for DPPC will be discussed further in the text.
Although applied here as in a classical potential60, the CHEQ formalism derives rigorously from the density functional theory of atoms in molecules62 based on Sanderson’s idea of electronegativity equalization63,64. Polarization is effected via the migration of charge density (condensed to a partial charge in the classical sense) between atomic species within a given molecule. The electronic density adjusts within the molecule so as to equalize the electrochemical potential (or equivalently, the electronegativity) at each point in the molecule. The direction and ease of flow are determined by physical properties of individual atoms as will be discussed. The reader is referred to the literature for more details42,43,56–60,62,65–67.
The electrostatic energy of a system of M molecules containing N atoms per molecule is:
where χ denotes the atomic electronegativities and η denotes the atomic hardnesses. The former quantity gives rise to a directionality of electron flow, while the latter represents a resistance, or hardness, to electron flow to or from the atom. The third term in Equation 1 is a standard Coulomb interaction between sites not involved in dihedral, angle, and bonded interactions with each other (the primed notation indicates a summation only over such sites). The second term represents the local charge transfer interaction generally restricted to within a molecule (no charge transfer) or some appropriate charge normalization unit. The last term is a Lagrange multiplier based constraint on total charge on a given normalization unit; this constraint helps to restrict charge equilibration (hence charge redistribution) over chemically relevant and distinct units68. We note that although the electronegativity and hardness follow exactly from the definitions of the electron affinity and ionization potential, they are considered here as empirical parameters to be determined as described below. Homogeneous hardness values (for each atom type) are parameterized as discussed in Patel and Brooks59. Heterogeneous elements (interaction elements between different atom types) are derived from the individual atom type values based on the combining rule58:
where Rij is the separation between atoms (or more generally sites) i and j. This local screened Coulomb potential has the correct limiting behavior as for separations greater than about 2.5 Å. This interaction is computed for 1–2, 1–3, and 1–4 sites (sites included in bonds, angles, and dihedrals). Sites in a molecule separated by 5 or more sites interact via a Coulomb interaction; in the case of interacting molecules, the interaction between sites on different molecules is again of the Coulomb form.
Regarding the polarizability, the charge equilibration model is indeed a polarizable model as the molecular polarizability, α, can be derived as follows:
where η̿ denotes the molecular hardness matrix, and and are the β and γ Cartesian components of the atomic position vector, respectively. A more detailed derivation can be found elsewhere68. The hardness matrix can be augmented to enforce charge constraints within a molecule68 for explicit calculations of the polarizability as were carried out in this study for the refinement of electrostatic parameters. In addition, the charge equilibration model, being an all-atom representation with partial charges assigned to all atomic species, contains all higher-order electrostatic multipole moments, in contrast to point dipole polarizable models69–71 and Drude oscillator models45,72,73. As such, the charge equilibration model incorporates higher-order electrostatic interactions explicitly.
The charge degrees of freedom are propagated via an extended Lagrangian formulation that imposes a molecular charge neutrality constraint, thus strictly enforcing electronegativity equalization at each dynamics step. The system Lagrangian is
where the first two terms represent the nuclear and charge kinetic energies, the third term is the total potential energy, and the fourth term is the molecular charge neutrality constraint with λi the Lagrange multiplier for each molecule, i. The fictitious charge dynamics, analogous to the fictitious wavefunction dynamics in Car-Parinello (CP) type methods74, are determined with a fictitious charge “mass” (adiabaticity parameter in CP dynamics). The units for this mass are . The charges are propagated based on forces arising from the difference between the average electronegativity of the molecule and the instantaneous electronegativity at an atomic site.
We will discuss the charge equilibration model for DPPC in further detail below. We comment here that the polarizable TIP4P-FQ water model is used to model solvent-solvent and solvent-solute interactions. The TIP4P-FQ water model is a 4-site model, based on the original TIP4P water model of Jorgensen et al75. The charges reside on the hydrogen atoms and a virtual site situated along the perpendicular bisector of the HOH angle 0.15 Å from the oxygen atom. The model has been characterized in previous studies, and the reader is referred to the relevant literature for further details.
Simulations were carried out in the constant pressure, surface area, and temperature (NPAT) ensemble using the CHARMM molecular modeling package19,76. The DPPC bilayer system consisted of 72 lipid molecules and 2511 molecules of TIP4P-FQ43 water. The total system size was initially 48.075 Å × 48.075 Å × 69.0 Å (based on a surface area per lipid of 64.2 Å2 from experimental measurements77); the starting geometry was obtained from a simulation using the C27r nonpolarizable force field. Dynamics were propagated using a Verlet leapfrog integrator78 with time steps of 0.5 femtoseconds (fs). The system temperature was maintained at 323 K using the Nosé-Hoover78,79 method with a thermal piston mass of 3000 . Pressure was maintained at 1 atmosphere using the Langevin piston method with a piston mass of 750 amu in the z-direction only (bilayer normal)78, resulting in constant surface area. Particle Mesh Ewald80,81 (PME) summation with screening parameter κ = 0.320 was used in all simulations to account for the conditionally convergent long-range electrostatic interactions; a grid spacing of 1 Å was used for the Fast Fourier transform (FFT) grid. Several replicate simulations of varying lengths were run, for a total simulation time of ~ 90 ns. In addition, several replicate simulations (a total of ~ 60 ns) were run on this system using the nonpolarizable C27r force field for lipids and TIP3P water75. Simulation parameters were the same as described above, except that a 1 fs timestep was used.
For the simulation of DMPC-bound gramicidin A, a hexagonal DMPC system with dimensions 29.7 Å × 29.7 Å × 76.8 Å was used. The system contained 20 DMPC molecules, 1070 TIP4P-FQ water molecules, the gA dimer, 18 K+ ions and 18 Cl− ions. The crystal structure of the gA dimer was obtained from the PDB (entry 1JNO). This structure was equilibrated using the nonpolarizable model to obtain a starting structure for the current simulations. The simulation was carried out in the constant pressure and temperature (NPT) ensemble with the temperature maintained at 298 K via a Nosé-Hoover thermostat with a piston mass of 100 . Pressure was maintained at 1 atm with a Langevin piston mass of 100 amu. To prevent the gA dimer from drifting inside the bilayer, a center of mass constraint with force constant 10 was applied to the protein only. Other simulation parameters were the same as described above for the DPPC system. Only the polarizable CHEQ model was used for the gA simulations. Total simulation time was ~ 4 ns.
Atomic charge degrees of freedom (the partial atomic charges) are propagated within an extended Lagrangian formulation. Since the Nosé-Hoover charge dynamics does not inherently enforce strict charge neutrality (or charge conservation in general), during each Nosé-Hoover iteration, we enforce charge neutrality for individual normalization units by subtracting out the average excess charge (excess relative to the required total charge constraint) from each atom. This approach, which serves as an efficient means to ensure strict charge neutrality during the course of the simulation, has been applied and validated previously82. Regarding computational cost, the overhead per time step is only 10% relative to the fixed charge models. However, due to the propagation of the light electronic degrees of freedom, timesteps on the order of 0.5 fs are generally applied, while timesteps twice as large (1 fs) are typically used with fixed charge models. Thus, the time required for a simulation using the polarizable model is roughly 2–4 times greater than using the nonpolarizable model
Finally, we will discuss the parameterization of the CHEQ model. The protein model used for gA was transferred directly from earlier work59. The lipid force field used for DPPC builds on an earlier model59. The dihedral, electrostatic, and Lennard-Jones (van der Waals) parameters were refined to better reproduce ab initio and experimental data for a series of model compounds. These model compounds were chosen to represent the various regions of a lipid molecule. For the phosphatidylcholine headgroup, tetramethylammonium (TMA) and dimethylphosphate (DMP) were used. The alkyl tails were represented by hexane. Parameters for the glycerol group (ester linkage) were transferred from earlier work59.
We will briefly outline the parameterization strategy here. Much greater detail can be found elsewhere82,83. For the electrostatic parameters (defined in Section II A), η was scaled to reproduce ab initio (for TMA and DMP) or experimental (for hexane) molecular polarizabilities, and χ was modified accordingly to reproduce ab initio atomic charges. The Lennard-Jones parameters are defined by the Lennard-Jones potential,
where ε and rmin are the adjustable parameters. ε represents the potential well depth while rmin is the distance at which the potential is at a minimum. For DMP, these values were adjusted to reproduce ab initio interaction energies and hydrogen bond distances with water calculated at the MP2/cc-pVTZ level of theory83. In the case of TMA, these values were found to be in satisfactory agreement with the reference data so those parameters were not modified. For hexane, the Lennard-Jones parameters were fit to reproduce the experimental liquid density and vaporization enthalpy. This was accomplished by adjusting the parameters for all relevant atom types at once and computing the properties from short simulations using various sets of adjusted parameters. From this data, linear regression was used to extrapolate a set of parameters which reproduced both properties almost exactly82. Finally, the dihedral parameters were refined to reproduce target ab initio conformational energies. The dihedral potential is
where K, n, and δ are the adjustable parameters. Based on the observation that the initial dihedral parameters would require only relatively small perturbations to come into agreement with the reference data, only the amplitude parameters, K, were modified. n and δ were not modified from their original values59,84. For the tailgroups, the relevant torsions are CH3–CH2–CH2–CH2 and CH2–CH2–CH2–CH2. Pentane and heptane were used in addition to hexane; the reference data was taken from a previous study which parameterized the nonpolarizable C27r force field for alkanes84. The DMP parameters were similarly fit to reference data calculated at the MP2/cc-pVTZ level of theory83. The relevant torsion in this case is C–O–P–O. Finally, o-phosphorylcholine85 was used to model the two torsions (N–C–C–O and C–C–O–P) connecting the choline to the phosphate in the phosphatidylcholine headgroup. Using the CMAP method86,87 in CHARMM, the energy surface of the two dihedrals was fit to a reference surface generated by the C27r force field, which was originally fit to ab initio data.
One of the most basic properties characterizing the equilibrium bilayer structure predicted by the simulation is the number density profiles of various atomic species as a function of distance along the bilayer normal. Figure 1 shows the component number density profiles for water and atomic species of the lipid (headgroup nitrogen, oxygens, and phosphorus). The CHEQ model is consistent with the nonpolarizable C27r model, and both are consistent with previous studies88–90, indicating stability in the bilayer structure. Both models predict gaussian distributions of the component atomic species, and it is this gaussian represenation that is often used to connect structural models to x-ray scattering experiments. We point to the subtle difference in the extent of water penetration into the membrane interior between the polarizable CHEQ and the nonpolarizable C27r models. The top panel of Figure 1 demonstrates the penetration of polarizable TIP4P-FQ water further into the bilayer (water density is higher than the carbonyl oxygen atom density for positions 10 Å or less from the bilayer center. For the C27r force field, we observe, as has been also observed in prior studies, that the carbonyl oxygen region is essentially the boundary beyond which very little nonpolarizable water resides. This phenomenon is not related to the larger surface area per headgroup for DPPC (64.2 Å2) nor the increased temperature of 323 K at which the present simulations were performed, since in that case, one would expect to see geometric/packing and thermal effects allowing for greater water penetration for both water models. We believe that this effect arises solely from the polarizability of the water (primary) and lipid (secondary) as the structural and dynamic properties of the lipid bilayer are similar for both models (and compare well to experiment). We will discuss this in the context of a potential of mean force (PMF) in Section III B.
We compute the electron density profile by summing the individual component density profiles, scaling each by the appropriate number of electrons for that atom type. Figure 2 shows the total and component electron density profiles for both polarizable and nonpolarizable models. Both are consistent with previous simulations19,91. Furthermore, the computed component profiles are in very good agreement with the those fit to experimental data77. The most noticeable differences between the polarizable and nonpolarizable models appears in the vicinity of 15–20 Å. The component profiles suggest that this is due to subtle differences in the carbonyl/glycerol (CG) and phosphate contributions between the two models. The C27r model underestimates the peak heights of these components relative to the polarizable model. However, both models underestimate the CG component relative to the experimental profile, while the phosphate component is overestimated. These differences result in an overall higher peak in the total density profile around 20 Å for both models. We also point out that the polarizable model performs slightly better in reproducing the shape of the experimental profile in the vicinity of the interface. However, it is worth noting that the nonpolarizable CHARMM model was recently refined to better reproduce the experimental electron density profile of DMPC91. It is important to emphasize that the electron density profile was not used as a fitting criterion for the CHEQ model presented here. Therefore, although the nonpolarizable model has been improved to better reproduce the electron density, it is promising that the polarizable CHEQ model is in such good agreement with experiment without that property being considered explicitly in the fitting. Another item of note is that the experimental profile consists of fitting functions based on contributions to the electron density from various regions of the bilayer. The fitting functions used indicate a vanishing water contribution going towards the center of the bilayer77. This suggests that very little or no water is present near the center of the bilayer, but this data alone is not sufficient to gauge the significance of the presence of water in the bilayer interior.
A recent study by Erilov et al9 used electron spin-echo envelope modulation (ESEEM) spectroscopy to measure the penetration of water into DPPC bilayers. Their calculations indicate that the concentration of water in the bilayer interior near the top of the alkyl chain (C4–C7) is ~ 2.9 M. A previous study that estimated the intramembrane water concentration based on more indirect measurements of hydrogen bonding of water to spin labels10 suggest the water concentration around C4 is ~ 1.0 M. From our simulations, based on the water density profiles shown in Figure 1 and the density profile for the C4 carbon on the alkyl chain (data not shown), the water concentration around C4 is ~ 0.06 M for the CHEQ model and ~ 0.04 M for the C27r model. One of the studies just mentioned10 also gave an estimation of the water concentration closer to the center of the bilayer, ~ 0.3 M. Based on the water density data from our simulations, we find the water concentration up to 5 Å from the center is approximately constant at 7.61 × 10−4 M for the CHEQ model and 5.75 × 10−5 M for the C27r model. The polarizable and nonpolarizable models give quite different values for the water concentration close to the bilayer center, an effect that will be discussed further in Section IIIB below. In addition, both models give concentrations much lower than those suggested by experiment. It is worth noting that the experimental studies9,10 used nitroxide-containing spin labels on different positions on the alkyl chain. It is possible that the presence of these spin labels, a polar group in an otherwise nonpolar environment, caused some artificial enhancement of the water concentration since water can hydrogen-bond with these spin labels. Granted, even with the polar spin labels the hydrophobic enviroment of the bilayer interior is still less favorable to the presence of water compared to the bulk solution. However it is possible that even the relatively low concentrations measured are higher than those of water in a membrane without such spin labels. We acknowledge that this experimental data indicates a significant role of water permeation into the bilayer, but we do not believe it can be accurately compared to our simulations.
An estimate of the free energy cost for water transfer across the membrane can be obtained from the potential of mean force (PMF). Taking the water density profile in Figure 1 as an approximation of the probability of finding a water molecule at a particular z-coordinate, the PMF relative to the bulk solution is given by:
where ρ0 is the bulk density. The calculated PMF profiles for the polarizable CHEQ and nonpolarizable C27r models are shown in Figure 3. The current estimate shows a free energy barrier of ~ 4.5 kcal/mol using the polarizable CHEQ model. The CHEQ model allows more water density near the center of the bilayer as a result of polarization effects between the water and lipid that are not captured by the nonpolarizable C27r model, which allows much fewer water molecules into the membrane during the timescale of the simulations resulting in a less precise free energy estimate of ~ 6 kcal/mol. It is worth noting that since the polarizable water model is able to respond to electrostatic effects of the local environment, it is able to adjust to a lower dipole moment in the bilayer interior (see Section III F 1). This suggests that desolvating a polarizable water molecule constitutes much less of a penalty than desolvating a nonpolarizable water molecule, which has a higher dipole moment that does not change with its local environment.
Interestingly, an early study of DPPC by Marrink and Berendsen92 suggested a free energy barrier on the order of 6 kcal/mol, similar to that just mentioned for the nonpolarizable C27r model. In addition, the shape of the C27r water PMF in Figure 3 is qualitatively similar to that obtained by Marrink and Berendsen, exhibiting a shallow minimum at the bilayer center. Further, recent studies11 have decomposed such a profile into a four-region kinetic model. We acknowledge that more rigorous calculations of free energy profiles of solutes are warranted, and work in our laboratory continues towards that end.
To further characterize the water permeation properties of our simulated lipid bilayer, we have calculated profiles of the free volume as a function of distance along the bilayer normal. Such free volume analyses have been used previously92–94 to obtain information on the ability of various solutes to “fit” into the membrane interior. For our purposes, this analysis may also provide further insight into the source of the differences in water penetration between the CHEQ and C27r models. We compute the free volume profile as described by Marrink et al93. Briefly, the simulation system is first mapped onto a 3D grid. All grid points that fall within the van der Waals radius of an atom are defined as occupied volume, while the remaining (empty) grid points are free volume. The ratio of the number of empty grid points to the total number of grid points in an interval is then the free volume fraction in that interval. To obtain a profile as a function of distance along the bilayer normal (z-axis), the system is divided into intervals of size Δz along the z-axis. Given a grid of size Gx × Gy × Gz and a system with dimensions X × Y × Z, the spacing between grid points is then in the x-direction and so on. We define our interval Δz to be equal to the grid spacing . In other words, each interval is only one grid spacing wide in the z-direction. Therefore, the total number of grid points contained in an interval is simply GxGy. To obtain the free volume fraction in a given z-interval, then, we simply count the number of grid points in that interval that fall inside an atom’s van der Waals radius and divide by the total number of grid points in that interval, GxGy.
The top panel of Figure 4 shows the empty free volume profile of the membrane as a function of distance along the bilayer normal using a grid spacing of 0.2 Å. Atomic van der Waals radii were derived from the force field Lennard-Jones parameters, described in more detail in Section III F 2. Since this analysis is more computationally expensive than the others, only the last 5 ns of the trajectory was used. The resulting profiles are qualitatively similar to those obtained previously for membrane bilayer systems92–94, however we note that those studies used a wider grid spacing of ~ 0.5 Å. There are two major observations to be made here. First, we note the “shoulder” apparent in the free volume profile obtained from the CHEQ model; this feature is not present in the C27r model. This “shoulder” feature is also present in the water density profile shown in the middle panel of Figure 4 (expanded from Figure 1). Alone, the CHEQ water density profile suggests an effect due only to the polarization of water compared to C27r. However, together with the free volume profiles, it suggests that there is a nontrivial difference in the bilayer structure around the headgroup region, ~ 20 Å from the center, between the two models. Second, both models show approximately the same fraction free volume at the center of the bilayer; however, in the region from 0 to 10 Å from the center, the CHEQ model shows less free volume than the C27r model. Further, this feature suggests a constriction in the membrane near the center, as the free volume decreases and then increases. This is opposite what would be expected based on the water density profiles in Figure 1 (reproduced in the bottom panel of Figure 4) the water density PMF in Figure 3, all of which show more water density in the bilayer interior with the CHEQ model. In other words, since there is more of a constriction present with the CHEQ model than the C27r model, one would expect less water to penetrate. The fact that this is not the case suggests that the increased water penetration with CHEQ arises from polarization effects and not structural differences, since the CHEQ model allows less empty free volume yet more water inside the bilayer.
The middle panel of Figure 4 shows the empty free volume profiles computed using experimentally derived van der Waals radii95 (set 1 in Table II). Using these values results in a higher overall free volume, but the profile is qualitatively the same, with the same features as noted above. Here, the same set of radii was used for both models. This is worth noting, since it suggests that the difference in free volume between the two models is a significant structural effect due to the force field and not due to small differences in the Lennard-Jones parameters between the two models. We also computed profiles using experimentally derived covalent atomic radii (set 2 in Table II). However, as those radii are much smaller, unrealistically high free volume fractions on the order of 0.8 resulted (data not shown). Finally, we considered the profile of accessible free volume as described by Marrink et al93. This analysis is identical to the empty free volume analysis except that the van der Waals radius of the “solute” or “penetrant molecule”, in this case water, is added to the radius of each atom in the system. This results in a profile of the fraction of free volume accessible to a solute molecule, or in other words the amount of volume where the solute can “fit”. For our DPPC system, we found very little water accessible free volume with either model (data not shown). In this respect, the empty free volume profiles are more informative, as they show a clear structural difference between the models as discussed above.
The orientation of the lipid headgroups in the bilayer can be analyzed by calculating the angle between the phosphorus-nitrogen vector (P) and the bilayer normal. The importance of the orientation of the P dipole in modulating the hydration of the interfacial region has been recently discussed23. The CHEQ model predicts an average headgroup orientation of 78.3° relative to the bilayer normal, while C27r predicts 73.3°. Both are consistent with previous studies of DPPC in which the headgroups favored an orientation of ~ 80–90° relative to the bilayer normal88,94,96. In addition, the distributions are equally as broad as those calculated in a recent study by Siu et al23 which used the generalized AMBER force field (GAFF)97 and a united-atom model with the Berger force field98,99 as well as the CHARMM force field. Finally, it is worth noting that the probability distribution obtained from the CHEQ model is narrower than that obtained from the nonpolarizable C27r model. This seems to suggest a more restricted range of rotation for the phosphate and/or choline groups. In the case of the phosphate group, this effect could be attributed to this higher cis energy barrier as discussed in our previous study of DMPC83. With a higher barrier to overcome, the range of likely orientations of the phosphate group is more restricted.
We also examined the orientation of the carbonyl groups by calculating the angle between the C=O vector and the bilayer normal. On average, the carbonyl groups on the sn-1 chain assume an orientation of 68.6° relative to the bilayer normal with the CHEQ model, while the C27r model predicts an average sn-1 carbonyl orientation of 76.3°. For the sn-2 chain, CHEQ predicts an orientation of 77.8° while C27r predicts 75.5° on average. An experimental study of carbonyl group orientations in DPPC bilayers100, indicated that on average the sn-1 carbonyl groups are oriented at 62 ± 2° and the sn-2 carbonyls are oriented at 66 ± 2° relative to the bilayer normal. As observed above with the P orientation, the nonpolarizable C27r model shows a broader distribution indicating a less restricted range of rotation compared to the CHEQ model. These observations suggest a need to refine the torsional parameters of the ester portion of our current force field, which were transferred from earlier work59.
The order and dynamics of the lipid tailgroups can be quantified by the variation of deuterium order parameters as a function of position along the alkyl chain. The deuterium order parameter, SCD, is calculated as:
where θ is the angle between a particular carbon-hydrogen vector and the bilayer normal, and
is the second Legendre polynomial. Quantitatively, SCD is equivalent to the SZZ component of the NMR quadrupolar splitting tensor23. Qualitatively, it is a measure of the order and orientational dynamics of the tailgroups that can be compared to experiment. Figure 5 shows the magnitude of the calculated order parameters as a function of position along the alkyl chain. Both the polarizable CHEQ model and the nonpolarizable C27r model are shown for comparison. The data are also presented in Table I. Both the CHEQ and C27r models reproduce experimental trends8,66,101,102, with the magnitude being slightly underestimated in most cases. Furthermore, the CHEQ model shows greater order approaching the ester region compared to C27r and experiment. This may be attributable to the more polar environment closer to the headgroup area, which using the CHEQ model results in a larger polarization effect on the carbons closest to the headgroup, particularly the C2 carbon. Further, we note that the current model is not able to discern between the two C2 hydrogen atoms that are steroscopically different. This data is included in Table I. The 2R and 2S order parameters are indistinguishable compared to experiment, and consistent with previous estimates of this property using empirical force fields23,84. As mentioned before, this may be attributable to a more polar environment as the C2 carbon is closer to the headgroup; however, this appears to be a general deficiency of classical force fields and further work continues to address the issue.
The conformational disorder of the tailgroups can also be measured by the fraction of gauche conformations. We define a torsion angle as being in a gauche conformation if −150° < < 150°, i.e. if is outside the range of trans conformations. We calculated the fraction of gauche conformations at each position along the alkyl chains. Specifically, for carbon Cn, we say there is a gauche conformation if the torsion angle Cn−2 − Cn−1 − Cn − Cn+1 is in a gauche conformation. It is therefore only meaningful to define gauche conformations at positions C4 through C15 since to go outside that range would either go beyond C16 (the terminal carbon) or include non-tailgroup atoms. Figure 6 shows the fraction of gauche conformations as a function of position along the alkyl chain. These are in qualitative agreement with a previous study by Klauda et al84; however, since that simulation was done under different conditions, most notably a lower temperature, we do not expect the trends to be reproduced exactly. We also show in Figure 6 experimental data103 for comparison. We note that the uncertainty in the experimental measurements is somewhat high (approximately 20%)104. The increase in disorder at the end of the tailgroup is in very good agreement with experiment, as is the gauche fraction at the C6 position. However, neither model reproduces the decrease in disorder at the C2 position. We note that the CHEQ model predicts slightly fewer gauche conformations at C2 than C27r, a result consistent with the increase in the deuterium order parameters predicted by CHEQ at that position compared to C27r (Figure 5).
The nature of water-lipid interactions can be observed through the radial distribution function (RDF) between relevant water and lipid atom pairs. We compute the RDF for two particular atom types by counting all atom pairs separated by a distance r over all frames in a trajectory. The resulting distribution is normalized by dividing by the number of frames and the average density of the integration shell, resulting in a function that approaches unity at large values of r. In addition, pairs on opposite sides of the bilayer are excluded to ensure that the peaks represent interacting pairs that are not separated by the bilayer. Figure 7 shows RDFs of water oxygens to various headgroup atoms. Both polarizable and nonpolarizable models show significant interactions. In the choline nitrogen to water oxygen RDF, the first peak around 4.5 Å indicates a significant solvation structure around the choline group, although the water oxygen does not interact directly with the nitrogen atom. The second, less prominent, peak around 7.5 Å defines a second solvation shell surrounding the choline group. These results are in very good agreement with a previous study by Lopez et al89 which looked at the interactions between water oxygens and headgroup atoms in a DMPC bilayer using the AMBER force field and SPC/E water. The water oxygen to phosphate phosphorus RDF indicates a difference between the polarizable and nonpolarizable models. The nonpolarizable model shows a slightly stronger first solvation shell around the phosphate group, a feature similar to that observed by Lopez et al89 and indicative of bridging interactions between waters solvating the headgroup. Another difference between the models can be seen in the water oxygen to phosphate non-bridging oxygen RDF. The nonpolarizable model, consistent with Lopez et al, shows a first solvation shell around 3 Å, a second solvation shell around 5 Å, and possibly a weak third hydration shell further out. The polarizable model shows the first and second hydration shells, but ~ 0.5 Å more distant than the nonpolarizable model. There is also virtually no indication of the weak third hydration shell. Furthermore, the second peak is more well-defined in the polarizable model, it does not appear to be made of two “lobes”, as in the nonpolarizable model. Lopez et al attribute the first lobe to a bridging water interaction and the second to the bulk second hydration shell. The fact that the closer lobe does not appear in the polarizable model suggests that the briding interaction is not as strong. Finally, the water oxygen to carbonyl oxygen RDF indicates a first and second hyrdation shell in both polarizable and nonpolarizable models. However, the first hydration shell is ~ 0.5 Å closer in the nonpolarizable model, and the second hydration shell is weaker in the polarizable model. In general, these RDF peaks represent the distances between interacting lipids and waters.
Polarizable force fields allow charges to respond to local electrostatic environment. One would expect, therefore, a variation in the average molecular dipole moment of water when moving along the bilayer normal from the membrane interior to the bulk solution. Figure 8 shows the molecular dipole moment distributions obtained by averaging the dipole moment of water molecules found in slabs of width 0.05 Å along the interface normal. In the bulk solution, the water dipole moment plateaus at a value of 2.55 Debye, ~ 0.8 Debye above the experimental gas phase value of 1.85 Debye95. This is consistent with the bulk TIP4P-FQ dipole moment. The average molecular dipole moment decreases to ~ 1.88 Debye in the membrane interior. It is worth noting that the polarizable water model is able to capture the effect of the bulk condensed phase environment on the water molecular dipole moment, namely the effect of an enhanced dipole relative to the gas phase. The value of 2.55 Debye is reasonable considering that there is no consensus value for the liquid phase dipole moment of water, with empirical and ab initio estimates ranging from 2.5–3.0 Debye105–109. We also observe that the dipole moment does not decrease exactly to the gas phase value at the center of the bilayer. It is fairly constant from the center of the bilayer until ~ 4 Å from the center, indicating that there is a significant region near the center of the bilayer in which the water dipole affected only by the local alkyl chain polarizability and is shielded from the headgroup region. Moving further out toward the interface, the average molecular dipole moment increases monotonically. Therefore, there is a significant region over which the dipole moment is enhanced in the lipidic environment. This suggests a nontrivial dielectric effect of the polarizable membrane; that is, the membrane has a dielectric constant different from unity. This effect arises not only from the local polarizability of the lipid chain but also from longer-ranged effects from the polar headgroup, which as mentioned before are negligible close to the center. In our previous study of DMPC83, we estimated a total interior dielectric constant of ~ 2 based on the Clausius-Mossotti relation and the calculated polarizability of tetradecane. This estimate is also appropriate for DPPC, and it is also consistent with the computed dielectric permittivity as will be discussed in Section III F 3. Finally, it is worth noting that as the water dipole moves away from the center of the bilayer, it begins to interact with the stronger electric fields generated by the phosphate, ester, and choline groups.
The membrane dipole potential is difficult to reproduce from atomistic simulations. In part, this difficulty is due to the rather ambiguous definitions of the interfacial potential as an experimentally measured physical quantity. For lipid bilayers based on phosphatidylcholine, such as DPPC, the potential has been measured indirectly using conductance measurements, resulting in values ranging from 220 to 280 mV110–112. Recently, novel cryo-EM techniques have been used to estimate membrane electrostatics with electrons as probes, measuring values on the order of 510 mV4; these studies concerned lipids containing both ester and ether chains.
The Galvani potential113 is calculated from simulations by twice integrating the charge density as a function of distance from the bilayer center:
where z = 0 is at the center of the bilayer, ε0 is the permittivity of vacuum, and ρ(z) is the charge density obtained by dividing the system into slabs along the bilayer normal, summing the charges in each slab, and dividing by the slab volume. For component contributions to the interfacial potential, charge densities are calculated for the individual molecular species, then twice integrated as in Equation 10. We note that for constructing the profile as a histogram and for numerical integration, we use a spacing (Δz) of 0.05 Å. Figure 9 shows the total and component contributions to the interfacial potential for both the polarizable CHEQ and the nonpolarizable C27r models. The CHEQ model predicts a potential drop of 0.9462 ± 0.001 V, while CHARMM27r predicts a drop of 0.8573 ± 0.002 V. Both models overestimate the magnitude of the potential drop compared to experiment, with the polarizable model displaying a drop that is ~ 0.1 V deeper. Regarding the contributions from lipid and water, Figure 9 shows that the polarizable model predicts lower individual absolute potential changes across the interface, though their directions are the same in both models (as they must be). Though the dipole potential computed in this manner appears to be in gross error relative to experiment, one need consider the variability in the experimental measurements, and their physical meaning. It has been argued114 that the potential drop across the bilayer-water interface is not a well-defined, experimental measureable. In fact what appears to be a more robust comparison to experimental observables is the difference in water-vapor and lipid monolayer-vapor interfaces114. Though for the present work we do not consider that specific system (this will be explored in a future work), we anticipate that the inclusion of polarizability should allow a better description of the difference in interfacial potential between water-air and monolayer-air.
It is also useful to decompose the water contribution to the interfacial potential into dipole moment and quadrupole moment contributions. Consistent with previous work115,116, we define the dipole moment contribution to the interfacial potential as
where the dipole moment density μz(z) is defined as
The integration in Equation 11 is taken from a point far into the bulk solution (z = z0) to the center of the bilayer (z = 0). The indices m and i in Equation 12 refer to a molecule and an atom in that molecule, respectively. The coordinate zim is taken relative to an arbitrary center which we choose to be the oxygen atom in the case of water. Similarly, then, we define the quadrupole moment density as
and the quadrupole moment contribution to the interfacial potential as
where the reference value is taken to be zero. Figure 10 shows the water dipole and quadrupole moment contributions to the interfacial potential calculated using both CHEQ and C27r models. For CHEQ, the dipole and quadrupole moment contributions are 0.423 V and 1.06 V, respectively, while C27r predicts 1.84 V for the dipole contribution and 0.758 V for the quadrupole contribution. In both cases, the sum of these contributions is equal to the total water component (1.48 V for CHEQ and 2.60 V for C27r). Comparing to the profile of the total water component in Figure 9, we note the small peak in the area of the interface. It is clear from Figure 10 that this feature arises from the dipole moment contributions for both models. We also note that the CHEQ model predicts a larger contribution from the quadrupole component while C27r predicts a larger contribution from the dipole component. This can be attributed to the explicit polarization accounted for by the CHEQ model. The TIP3P water model has a fixed geometry and charges, and therefore a fixed dipole moment as shown in Figure 8. Since the polarizable model is able to adjust to a lower dipole moment in the membrane interior, the contributions from water molecules in the interior will be smaller and therefore the total contribution will be smaller.
Before moving to a discussion of dielectric permittivity profiles through the bilayer-water interface, we consider the effect of alternate representations of the charge density in defining the interfacial dipole potential from solution of Poisson equation. Physically, charge distributions described as point charges are artificial constructs to represent a localized charge density; consider the various approaches for condensing information about the electronic wavefunction to a single metric embodied by the charge. An obvious alternative is to diffuse the partial charge density over a suitable spatial region. In the following, we consider an alternative method for computing the interfacial potential in which the atomic charges are represented as gaussian distributions. This approach has been used previously94. We represent the atomic charge distribution for an atom i as
where σ is the width of the gaussian, zi is the atomic z-coordinate (i.e. the center of the distribution), and is the total atomic charge on atom i. The factor normalizes the gaussian distribution such that the integral over all z is equal to . The width of the gaussian, σ, is determined as one third of the van der Waals radius, συdW. This is to ensure that most of the charge is contained within the van der Waals radius of the atom, since ~ 99% of a normal distribution is contained within ±3σ. The van der Waals radii συdW are related to the force field nonbonded parameters rmin by the relation
which is derived from the Lennard-Jones potential and the definitions of rmin and συdW as the separation at which the interatomic potential is minimized and the separation at which the interatomic potential equals zero, respectively.
The interfacial potential is computed as in Equation 10 except with this alternative definition of the charge density. Since this computation is more time-consuming than the previous analyses, only the last 5 ns of the trajectory is considered here. In addition to using radii derived from the force field Lennard-Jones parameters, we also test the effect of using experimentally derived atomic radii95 to determine the width of the distributions. Set 1 contains nonbonded contact van der Waals radii refined with crystallographic data, and Set 2 contains covalent contact radii used to determine whether atoms in a crystal are covalently bonded. These values for relevant atoms are shown in Table II. Figure 11 shows the profiles computed using these radii for both CHEQ and CHARMM27r models. Profiles calculated using point charge density (first method) are also shown. We note that as the radii get larger overall, resulting in broader distributions, the magnitude of the potential drop decreases. We confirmed this effect by artifically scaling the radii by several factors (data not shown) and also found that the profile matches very closely with that calculated from the point charge density when the distributions are made very narrow. In this sense, the point charges could be viewed as infinitely narrow gaussian distributions. As a final note, broader distributions seem to have a stronger effect on the CHARMM27r model than the CHEQ model. Any discrepancy in the trend between the force field radii and the experimental values between the CHEQ and CHARMM27r models is likely due to the difference in force field parameters between the two models. We do observe a smaller drop using Set 1 than Set 2, consistent with the observations mentioned previously, since the radii in Set 1 are larger than those in Set 2.
Next, we consider the application of the approach developed by Stern and Feller3 for computing the profile of the parallel (in-plane) component of the dielectric permittivity as a function of distance along the bilayer normal. We compute profiles of the parallel component of the dielectric permittivity using Equations 71 and 26 of Stern and Feller3 for tin-foil boundary conditions:
where (z) is the local polarization density and is the total system dipole moment. a(z) represents the contribution from explicit polarization in Stern and Feller’s formulation. However, with the CHEQ model, that contribution is self-consistently included in the first term of Equation 18. Briefly, the polarization density is calculated using an approach similar to that of Stern and Feller, in which a set of bond-charge increments (bjk) are defined for each bond such that each atom (j or k) associated with a bond jk receives an amount of charge ±bjk. The total charge on an atom is then the sum over all the bond-charge increments to which the atom belongs:
Given the atomic charges, we obtain the bond-charge increments bjk by inverting the matrix via singular value decomposition or via straightforward inversion for well-conditioned matrices. This inverse is computed once for a given molecule using the minimal topology description of bond-charge increments (i.e. chemical bonds) and reused for analysis of the trajectories. Given the bond-charge increments, we compute the polarization density (z) as a sum of local bond dipole moments in a bin of width dz at a position z. Following the approach of Stern and Feller, we use:
where A is the cross-sectional area (constant in our simulations) and the bond dipole jk is determined from the bond-charge increments bjk and the position vectors j and k simply as:
Figure 12 shows the parallel component of the dielectric constant as a function of distance along the bilayer normal (z-axis) for the polarizable CHEQ and nonpolarizable C27r models. The bulk water values are consistent with the properties of the pure solvent models: 72.77 ± 5.8 for TIP4P-FQ (closer to experiment) and 91.22 ± 2.3 for TIP3P. In the bilayer interior, the CHEQ model exhibits a dielectric of 1.31 ± 0.6 and CHARMM27r exhibits a dielectric of 1.10 ± 0.5. Both models approach a value of 1 in the interior, which is expected for alkyl species. The headgroup regions display larger dielectric constants, which has been attributed to the large magnitude of the dipoles (19–25 D)3 in these regions. In addition, the errors are largest in this region, especially for the polarizable model. Such large errors in the headgroup dielectric are consistent with a study by Nymeyer and Zhou117 in which the dielectric profile for a bilayer was computed by free energy and Poisson calculations. The fact that the headgroup dielectric is so much higher with the polarizable model can be attributed to the effects of polarization implicitly included in the contribution from the total system dipole. For a fixed-charge model, polarization would need to be included explicitly via some effective polarizability as a function of distance from the bilayer center. However, like Stern and Feller we have not done so here.
Finally, we apply our polarizable force field to simulations of gramicidin A embedded in a DMPC bilayer. Gramicidin A is a membrane-bound ion channel which has been studied extensively through both experiment and simulation1,12,14,118. As mentioned in Section II B, a center of mass constraint was applied to the gA dimer to prevent drift of the channel within the bilayer. However, another measure of the stability of the channel structure is the tilt of the dimer axis (defined as the axis connecting the monomer centers of mass) relative to the bilayer normal. Previously, Allen et al found the tilt to be ~ 12° on average14. In the top panel of Figure 13, we plot the time profile of the tilt angle. On average, we obtain a value of 9.4 ± 2.1°, reasonably consistent with the simulations of Allen et al.
The stability of the channel structure can also be measured by the root mean square deviation (RMSD) of the backbone atoms from the experimental crystal structure (PDB entry 1JNO). The middle panel of figure 13 shows the time profile of RMSD from the crystal structure over the simulation. On average, we observe a RMSD of 0.68 ± 0.05 Å relative to the crystal structure, which indicates very little structural deviation, especially considering that at least some deviation from the crystal structure is to be expected due to different environments.
Another measure of the stability of both the bilayer and the gA channel is the variation in surface area over the course of the simulation. In the bottom panel of Figure 13, we show the time profile of the surface area. On average, we observe a surface area of 753.6 ± 11.5 Å2. For comparison, we use a reference value based on data derived from experiment. Woolf and Roux119 estimated a value of 250 Å2 for the surface area of the gA channel. A experimental study by Kučerka et al6 obtained a value of 60.6 Å2 per lipid for DMPC. Our hexagonal system consists of 20 DMPC molecules, 10 per leaflet. Therefore, based on this data, we obtain a value of 856 Å2 for a system of 10 DMPC molecules per leaflet plus the gA dimer. This may seem inaccurate, but we note that 250 Å2 is still only an estimate of the gA channel surface area; there are estimates ranging from 110 to 305 Å2 derived from experiment119.
As discussed in Section III F 1, the CHEQ model is able to capture the variation in the water molecular dipole moment going from bulk water into the bilayer interior. In the case of the gA channel, one would expect water molecules interacting with the channel to have an enhanced dipole moment relative to water molecules in the bilayer far from the channel. We find that near the bilayer center, water molecules far from the channel (> 8 Å from the dimer axis) have a dipole moment of ~ 1.85 D, consistent with the water dipole moment profile for the DPPC bilayer shown in Figure 8. However, within the channel (≤ 8 Å from the dimer axis), the average water molecular dipole moment is found to decrease to ~ 2.45 D compared to the bulk value of 2.55 D.
We have presented results of all-atom molecular dynamics simulations of a DPPC bilayer using the recently refined CHEQ model for lipid bilayers83. We have shown that the charge equilibration method is capable of modeling induction effects in such systems. Component number and electron density profiles as well as deuterium order parameters are reasonably reproduced. The CHEQ method is able to model the variation in the molecular dipole moment of water moving from the bulk into the bilayer interior, something inaccessible to fixed charge models. Furthermore, the CHEQ model is able to capture the electrostatic response of various atomic species as they go from the gas phase to the condensed phase, as demonstrated by the distribution of atomic charges throughout the simulation. Most atom types take on charges different from the gas phase values to which they were parameterized; this is indicative of a response to the bulk environment of the lipid which is far different from a gas phase environment.
Consistent with the assertion that the lack of explicit polarization may be a nontrivial simplification of current non-polarizable force fields, the CHEQ model predicts subtle but nontrivial differences in several properties relative to the nonpolarizable C27r model. Water density profiles show that the polarizable model allows for the presence of more water near the bilayer center, an observation consistent with the ability of the polarizable model to account for the decreasing water dipole from the bulk to the interior. Based on an analysis of the fraction of free volume available to both polarizable and non-polarizable water molecules, the presence of polarizable water molecules appears to be driven by a smaller electrostatic free energy penalty resulting from the ability of the polarizable water to reduce its dipole moment. Ongoing studies detailing more specifically the free energetics of water and permeant transfer through the interfacial regions using polarizable models will be presented in the future.
We have also demonstrated the suitability of the present polarizable force fields for modeling membrane-water-protein systems with short simulations (4 ns) of the gramicidin A ion channel embedded in a DMPC bilayer solvated using TIP4P-FQ water molecules. The current work sets the stage for considering the effects of protein, solvent, and lipid polarizability on permeation free energetics in narrow ion channels such as gA; ongoing work in our laboratory will report on the influence of polarization in defining the barriers to permeation in this classic model system.
The authors gratefully acknowledge support from the National Institute of Health sponsored COBRE (Center of Biomedical Research) Grant Numbers P20-RR017716 at the University of Delaware (Department of Chemistry and Biochemistry) and P20-RR015588 (Department of Chemical Engineering).