Sandeep Patel (Assistant Professor, Chemistry & Biochemistry)
Integral membrane proteins have become an arena of intense effort of fundamental research due to the ubiquitous nature of these macromolecules. Comprising roughly one-third of the human genome, they are implicated in myriad physiological function and, unfortunately, dysfunction. For instance, for normal physiological functioning, integral membrane proteins are involved in signaling processes, passive and active transport, and interfacial enzymatic processes [2]. Of course, one cannot speak about integral membrane proteins independent of the lipidic context in which they function. Lipid membranes in their own right have garnered much attention as well, with particular focus on membrane properties and behaviors including structural deformation (in the presence of small molecules as well as integral membrane proteins), and electrostatic properties such as the interfacial dipole (or total) potential [3] and dielectric constant variation with location in a bilayer [4]. Of particular recent interest have been protein-lipid interactions, and specifically, the interactions of charged and/or polar amino acid residues as they pertain to an understanding of the thermodynamics of structural and energetic stability of integral membrane proteins upon desolvation of such systems.
Complementing the enormous volume of experimental effort to understand lipid bilayers and integral membrane proteins, computational approaches based on the analytics of statistical mechanics (molecular dynamics and Monte Carlo simulations) have been an indispensable tool to understanding properties and processes in these systems at the atomic and molecular level, and providing insights at a resolution inaccessible to current state of the art experimental methods [2, 5-17]. These methods require information about the forces between interacting species; ideally, the necessary information would be contained in a quantum mechanical potential energy surface, which could then be directly applied to generate forces. However, simplifying assumptions must be made in order to arrive at tractable functional forms of a classical nature. Current state-of-the-art modeling approaches, thus, employ empirical potentials, or force fields to model inter-species interactions.
All-atom simulation methods invariably employ 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 literature [18-21]. 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. Though several models have appeared discussing application to proteins and nucleic acids, polarizable models for membrane systems have not yet been reported.
Under the support of the current COBRE, we have been able to develop and apply a prototypical lipid interaction model for DMPC that is fully polarizable. In the following, we discuss properties of model lipid bilayers from molecular dynamics simulations using these models.
Time profiles of the predicted surface area per headgroup are shown in Figure 1. Molecular dynamics simulations under constant pressure (1 atmosphere) and temperature (303 K) conditions maintained using a Nose-Hoover thermostat and Langevin piston barostat were performed using a Verlet leap-frog integrator; three independent trajectories were generated for sampling. The average headgroup area is predicted to be 55.8+1.1A2, lower than recent experimental data of 60.6+0.52 [22]and the fixed-charge prediction of 58.3+1.6A2. We observe for one of the trajectories rather large fluctuations of the surface area, an effect of small system size observed in earlier simulations on systems of comparable size[5]; nevertheless, quite satisfyingly, the final equilibrium values of the surface area converge for all three simulations. Ongoing refinement, as is currently pursued for most major force fields [23]will improve this system property.
One measure of the stability of the bilayer is the number density of various atomic species as a function of distance along the bilayer normal. Figure 2 shows the component density profiles for water and atomic species of the lipid (headgroup phosphorous, oxygens, and nitrogen). These are consistent with previous studies [24-26] as well as the non-polarizable CHARMM22 model. We note the subtle difference in the extent of water penetration into the membrane interior between the polarizable and non-polarizable force fields. Polarizability of both solvent and lipid tends to accommodate, free energetically, the presence of waters in the low-dielectric environment of the lipid interior. The potential of mean force for the water molecules derived from the computed density profiles shows a “slight” stabilizing effect introduced by the polarizability of the lipid and water as compared to the simulations without polarizability. This behavior has direct bearing on the thermodynamics associated with integral membrane proteins. The current results, though far from definitive, suggest the role of non-additive electrostatic effects (polarization) in contributing to the stability of lipid-exposed residues of integral membrane proteins. Continuing work in our laboratory will address the free energetics associated with transfer of polarizable amino acids into such systems much in the spirit of recent studies [8].
The dynamics of the phospholipid tailgroups are probed by the variation of deuterium order parameters with respect to the position along the alkyl chain; these are directly measurable via experiment. The deuterium order parameter, SCD, is obtained as:
where
is the angle between a particular CH vector and the bilayer normal, and
,
the second Legendre polynomial. This is equivalently the SZZ component of the NMR quadrupolar splitting tensor [17]. Figure 7 shows the magnitude of the calculated order parameters as a function of position along the alkyl chain. Also included for comparison are values from the nonpolarizable CHARMM force field and experimental values for the sn-2 chain [1]. Both polarizable and nonpolarizable models reproduce experimental trends [1, 27, 28] with the revised polarizable model showing slightly better agreement. The polarizable model displays less order in the bilayer center (more closely in agreement with experiment). As discussed earlier, this difference would also contribute to the enhanced water penetration into the bilayer core. We note that the current polarizable model is not able to discern between the two C2 hydrogen atoms that are stereoscopically different. The differences between the 2R and 2S order parameters are indistinguishable compared to experiment, and consistent with previous estimates of this property using force field based methods[17]. Furthermore, neither model reproduces the relative degree of order between the sn-1 and sn-2 chains. This general behavior is attributable to a more fundamental flaw in a classical representation of the local valence interactions rather than electrostatic or polarization effects. This appears to be a general deficiency for classical force fields [17] and further work continues to address this issue.
Polarizable force fields allow for a response to local electrostatic environment. As such, one anticipates a variation in the average dipole moment of water molecules when moving along the interface normal from bulk solution to membrane interior. Figure 10 shows the water molecular dipole moment distributions obtained by averaging the dipole moment of water molecules found in slabs of width 0.25 A along the interface normal. In bulk solution, the water dipole moment plateaus at a value of 2.60 Debye, 0.8 Debye above the gas-phase value of 1.85 Debye; this is consistent with the bulk TIP4P-FQ dipole moment. The profile exhibits a monotonic decrease to a value of approximately 1.9 Debye within the membrane interior. There are two items of note. First, the polarizable water model captures the condensed-phase environment effect on the local water molecular electrostatics via an enhanced dipole moment; we note that there still is no single consensus value of the liquid phase dipole moment of water, with values ranging from 2.5 – 3.0 Debye [29-33]. Second, we observe that the dipole moment does not fall exactly to the gas-phase value at the center of the membrane. Moving toward the interface from the center, the average molecular dipole moment increases monotonically; thus, there is a significant interior region over which the water dipole moment, though not at the gas-phase value, exhibits an enhanced electrostatic moment in the lipidic environment. This suggests a non-trivial dielectric effect exerted by the polarizable membrane, i.e., the membrane possesses a dielectric constant different from unity.
The membrane dipole potential still remains an elusive physical quantity to reproduce based on atomistic simulations. In part, the difficulty is based on the still rather ambiguous definitions of the interfacial potential as it is experimentally measured; this is a particularly sinister effect for assessing the accuracy of predicted liquid water contributions to the neat water and lipid system interfacial regions. For phosphatidylcholine-based lipid bilayers, such as DMPC, experimental values range from 220 to 280 mV; these are determined indirectly using conductance measurements [34-36]. More recently, novel cryo-EM approaches have been used to estimate membrane electrostatics with electrons as probes[3] measuring values on the order of 510 mV (0.51V).
The surface potential is calculated from simulations by twice integrating the charge density as a function of distance from the center of the bilayer:

where z=0 is at the center of the bilayer,
is the permittivity of vacuum, and
is the charge density obtained by dividing the system into slabs along the bilayer normal (z-axis) and summing the charges in each slab. Individual molecular species charge densities are twice integrated to obtain component electrostatic potential contributions. Figure 11 shows the total and component contributions to the interfacial potential for both polarizable and non-polarizable models. There is a slightly deeper potential drop for the polarizable model, on the order of 0.1 V, with both models overestimating the magnitude of the experimental range of values, 220 – 280 mV. We note here that the difference is within the range of fluctuations observed based on all-atom simulations using a variety of force fields.