INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

One of the current challenges in computational biology is the modeling of biological structures at scales approaching the cellular level. Appropriate description of solvent is crucial for such an endeavor; whereas bulk properties of water can be approximated by continuum electrostatic theories, it is generally accepted that an explicit solvent description is necessary at interfaces and molecular boundaries (Levy and Gallicchio, 1998; Makarov et al., 1998; Boresch et al., 1999; Lounnas et al., 1999; Simonson, 2000; Im and Roux, 2001). Among such interfaces, membrane bilayers are of particular interest due to their importance in a variety of biological phenomena, including molecular recognition of many signal transduction processes in biological cells (Hurley and Misra, 2000). The membrane potential is known to have a profound influence on the permeation of small molecules and ions (Roux, 1997), the insertion of proteins or small peptides (Ladokhin and White, 2001), and the adsorption state of various membrane binding proteins (Murray et al., 1999). Three major sources contribute to the existence of the membrane potential: the membrane voltage, the membrane surface potential, and the membrane dipole potential. The membrane voltage, or the transmembrane potential, arises from the imbalance of ions at the two sides of membranes. The membrane surface potential is due to the charged lipids mixed in the membranes. The origin of the membrane dipole potential is more subtle and typically attributed to the anisotropic distribution of the lipid head group orientation and the "overcompensation" of the dipolar responses from the membrane-covering water molecules.

The current work examines the applicability of continuum electrostatics theories, such as the Poisson equation, to the calculation of the membrane dipole potential across pamlitoyloleoylphosphatidylcholine (POPC) lipid bilayers.

Early experiments (Liberman and Topaly, 1969; Haydon and Myers, 1973) have shown that the passive permeability of zwitterionic bilayers is orders of magnitude greater for hydrophobic anions than for cations, leading to the hypothesis that the membrane dipole potential is responsible for the large disparity in membrane permeability of ionic species of opposite charge. More specifically, it was supposed that the internal membrane dipole potential in the hydrophobic core of the bilayer must be positive relative to the aqueous phase, leading to the discrimination against cationic species.

This positive dipole potential has been supported by many explicit solvent molecular dynamics (MD) simulations (Zhou and Schulten, 1995; Feller et al., 1996; Essmann and Berkowitz, 1999; Mashl et al., 2001; Saiz and Klein, 2002). Implicit solvent approaches, such as Poisson-Boltzmann theory, have the advantage of reducing the dimensionality of the simulation by approximating the electrostatic properties of the discrete solvent as a dielectric continuum. However, previous Poisson-Boltzmann calculations of zwitterionic bilayers failed to obtain the correct sign of the membrane potential in the hydrophobic core (Zheng and Vanderkooi, 1992). Attempts have been made to explain such failings by addressing the role of explicit water molecules using short-time (15 ps) restrained MD simulations (Gabdoulline et al., 1996). The present work illustrates how a consistent membrane potential profile can be obtained from implicit solvent calculations using a limited number of explicit solvent molecules at membrane surface.

	COMPUTATIONAL DETAILS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

MD simulations of the POPC bilayer

The molecular dynamics simulation of the POPC bilayer was conducted using the SANDER modules in AMBER 7.0 (Case et al., 2002). The POPC bilayer membrane consisted of 100 × 2 lipids and 19,577 TIP3P water molecules (Jorgensen et al., 1983) with equilibrated lateral dimensions of roughly 74.5 Å × 90.5 Å and with dimension of 118.1 Å in the membrane normal direction. The average surface area per lipid is 67.4 Å². The force field parameters were primarily adapted from previous parameterization for a DPPC bilayer (Smondyrev and Berkowitz, 1999), whereas the partial charges were recalculated using the RESP scheme (Bayly et al., 1993). Electrostatic interactions within the MD simulation were calculated using the particle-mesh Ewald method (Essmann et al., 1995). Bonds between hydrogens and heavy atoms were constrained with the SHAKE algorithm (Ryckaert et al., 1977), permitting the use of 2-fs time steps to integrate the equations of motion. The weak-coupling method (Berendsen et al., 1984) was used to couple the system to a thermal bath of 300 K and a barostat of 1 bar with coupling constants of 0.2 and 1.0 ps, respectively. The system center of mass motion was removed at every picosecond to avoid the "cold solute/hot solvent" problem (Harvey et al., 1998). The solute temperature and the solvent temperature were also monitored continuously, and no significant deviation from the whole system temperature was detected. The relative errors of both solute and solvent temperatures are within 0.01. Anisotropic scaling for the pressure regulation was used, and the surface tension of the membrane was observed throughout the simulation; the average surface tension was zero. It should be noted that anisotropic scaling for pressure regulation is required for controlling the surface tension of membrane simulations. Isotropic scaling for the pressure regulation will fail to achieve a tension-free membrane system if the initial box dimensions are not assigned properly, which is often the case. The entire system, consisting of 69,131 atoms, was simulated for a total duration of 2 ns.

Implicit solvent electrostatic calculations

The canonical equation for implicit solvent electrostatics is the Poisson-Boltzmann equation (Davis and McCammon, 1990; Honig and Nicholls, 1995), a mean field approximation that expresses the electrostatic potential due to charges and mobile counterions in a dielectric medium in terms of the solution to a nonlinear partial differential equation. However, due to the electroneutrality of the POPC system, mobile counterions were not included in the present work and the electrostatic potential calculated from solutions to the related Poisson equation:

(1)

(2)

was solved in a three-dimensional domain  in conjunction with a Dirichlet boundary condition (Eq. 2) on , the boundary of . In Eq. 1, (r) is the position-dependent dielectric constant, u(r) is a dimensionless electrostatic potential (scaled by e_c), e_c is the electron charge,  = (k_BT)¹ is the inverse thermal energy, k_B is the Boltzmann constant, T is the absolute temperature, N is the number of particles, z_i is the atomic partial charge, and r_i is the charge location. It should be noted that the current computational setup corresponds to experiments of salt-free systems. The condition in the Poisson calculations is also consistent with that of the MD simulation. Due to the generous "padding" of the electrostatic calculations with periodic copies of the MD simulation box (see below), the zero boundary condition (g(r) = 0) was used. Any artifacts due to the zero boundary condition (rather than periodic boundary conditions in the MD simulation) are offset by the explicit presence of additional periodic copies in the continuum calculations and the averaging used to calculate the electrostatic potential results (see below).

The Poisson equation was solved using the parallel focusing methods available in the Adaptive Poisson-Boltzmann Solver software package (Baker et al., 2001), described in more detail at http://mccammon.ucsd.edu/apbs. In the Poisson calculations, the basic cell of the MD simulation was replicated in the membrane plane direction to make a 3 × 3 × 1 supercell, as shown in Fig. 1. Atomic partial charges and radii were adopted from the force field values used in the MD simulation. The volume occupied by explicit atoms (either lipid or water) was assigned a dielectric value of 1, whereas the remainder of the domain was assigned the bulk water value of 78.54. The equation was solved on a 161³ grid focusing from a uniform 1.55-Å coarse mesh spacing to 0.466 Å on the finest level. This focusing resulted in a highly accurate solution in the central cell of the 3 × 3 × 1 supercell (see Fig. 1); the remaining eight cell copies were simply used to provide a degree of periodicity to the system, thereby mimicking the simulation conditions.

View larger version (104K): [in this window] [in a new window]   FIGURE 1   System used for electrostatic calculations. The basic MD cell was replicated in the membrane plane to form a 3 × 3 × 1 supercell. Quantitative measurements were only obtained from 75 Å × 75 Å framed region that contained the highest accuracy potential from the electrostatic focusing calculations.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Structural properties of the interfacial water

Fig. 2 a depicts the electron density profile of each hydration layer plotted along the membrane normal. The hydration shells were defined in terms of their contacts: water molecules in the first shell have direct van der Waals contacts with the lipid molecules, molecules in the second hydration layer directly contact the first layer, and so on. Increased electron density in the first hydration layer is due to the penetration of water molecules into the head group region of the bilayer; the wide dispersion of the lipid head-group region enhances the electron density profile of the first hydration layer. According to the above definition of hydration layers, these lipid molecules do not have direct contacts with the second hydration layer. It has been shown previously that very stable hydrogen bond and salt bridge patterns exist in the head group region of a zwitterionic bilayer (Pasenkiewicz-Gierula et al., 1999). The average numbers of water molecules per lipid, n_w, for these seven layers are 13.1, 20.9, 27.7, 34.1, 40.4, 46.5, and 52.6, respectively. As shown in Fig. 2 a, electron density distributions beyond the third layer are Gaussian, whereas asymmetry can be observed in the first and the second shells.

View larger version (30K): [in this window] [in a new window]   FIGURE 2   Location and orientational polarization of solvent hydration layers averaged over 2-ns MD simulation: layer 1 (), layer 2 ( ), layer 3 (  ·), layer 4 (· · ), layer 5 ( ), layer 6 (· · · ·), layer 7 (). (a) Electron density profile of hydration shells at the POPC membrane/water interface. (b) Polarization profile of water molecules in each hydration shell.

To investigate the effect of interfacial water structure on continuum electrostatics approximations, the position-dependent polarization for each hydration layer can be calculated and compared with the bulk dielectric typically used in implicit solvent simulations. The mean polarization (normal to the membrane surface) of hydration layer h at position z is defined by the ensemble average:

(3)

in which the sum runs over the N_h water molecules in layer h, µ_i is the dipole moment of the i-th water molecule, n is the membrane normal, and   denotes an average over all frames of the trajectory. As Fig. 2 b indicates, water molecules in the first three hydration layers demonstrate a higher degree of polarization with respect to the water in the more remote hydration shells. This suggests that water in the immediate hydration layers shows significantly different structure than "bulk" water and likely exhibits different electrostatic behavior than is included in typical continuum approximations, a result that is in good agreement with recent inelastic incoherent neutron scattering experiments (Ruffle et al., 2002).

Dielectric properties of the interfacial water

To further examine the differences between interfacial and bulk solvent, the dielectric properties of the system were calculated as a function of distance along the membrane normal. As mentioned above, the dielectric function (r) plays a central role in the Poisson and Poisson-Boltzmann equations. Most continuum electrostatics calculations assume  is a scalar-valued function of position, which assumes a bulk value in the solvent and a smaller value in the solute;  is usually constant throughout the solute and solvent and often changes discontinuously at the boundary between these two media. The use of a scalar-valued dielectric constant assumes polarization of the dielectric medium is orthogonal (no off-diagonal terms in the dielectric tensor) and isotropic (the three diagonal terms are the same). However, in the most general case, the dielectric should be represented by a tensor-valued function that includes the possibility of anisotropic and nonorthogonal polarizability. Additionally, nonlocal contributions can be included by integral equation approaches, e.g., Beglov and Roux (1996).

The normalized total dipole moment fluctuation tensor  is directly related to the dielectric tensor by linear response (de Leeuw et al., 1980) and fluctuation-dissipation (Kubo, 1966) theories. To assess the validity of the assumptions implicit in the typical implicit solvent simulation assignment of isotropic "bulk" dielectric values to all solvent-accessible regions,  was calculated for a series of 5 Å cubic boxes placed along the membrane normal:

(4)

in which the dipole moment components the N water molecules in each box were summed to obtain the total dipole moment components: (M =  µ_,i). Similar calculations were performed on a separate simulation of pure TIP3P water to obtain "bulk" values () for the fluctuation tensor.

Fig. 3 shows the fluctuation tensor profile of the POPC bilayer normalized by the bulk values (/^bulk, ^bulk =  =  = ). Like the polarization, the fluctuation tensor shows significant deviations from bulk behavior for hydration layers near the solvent-membrane interface, which typically vanish by the third or fourth hydration layer (±25 Å). It should be noted that the usual Kirkwood relation (de Leeuw et al., 1980) would not apply to such small, cubic systems (as used in above sampling); therefore it is not straightforward to relate this fluctuation tensor to the dielectric tensor. Because all of the diagonal components of fluctuation tensor converge to the bulk value when it is beyond 20 Å with respect to the bilayer center, the dielectric tensor is expected to exhibit a qualitatively similar behavior.

View larger version (17K): [in this window] [in a new window]   FIGURE 3   The total dipole fluctuation tensor () for interfacial solvent normalized by the "bulk" solvent values obtained from the MD simulation. Components are depicted as: xx/bulk (), yy/bulk ( ), zz/bulk (· · · ·), xy/bulk ( · ·), yz/bulk (·  ·), xz/bulk (). Hydration layers are approximately located at ±15, ±22, ±25, ±27, ±30, ±33, and ±37 Å from the membrane center.

Electrostatic potentials with varying amounts of explicit interfacial water

Fig. 4 a shows the electrostatic potential of the POPC as obtained by solution of the Poisson equation in the presence of varying numbers of explicit hydration layers. Of physical interest is the membrane dipole potential:

(5)

in which u_n(z) is the value of the potential calculated with n explicit water layers present, z is the distance along the membrane normal (z = 0 is the center of the bilayer), and Z_n is the location where the electron density of the n-th explicit water layer (or membrane, for n = 0) vanishes and the dielectric constant assumes the "bulk" value. This somewhat unusual definition for the dipole potential was motivated by the Dirichlet boundary conditions (u = 0) imposed in the Poisson calculations. As evidenced in Fig. 4 a, this condition causes the electrostatic potential to decay from a nonzero value at the explicit-implicit interface to zero at the domain boundary. This difference can be attributed to the nonzero constant background potential implicit in molecular dynamics simulations and the presence of periodic boundary conditions (Figueirido et al., 1995). This difference does not affect the calculated membrane dipole potential as this quantity is independent of the chosen background potential. If the Poisson calculations were performed with periodicity in the x and y direction, it is expected that the above decay would no longer exist, and u_n(Z_n) would be the same as u_n().

View larger version (24K): [in this window] [in a new window]   FIGURE 4   Membrane electrostatic potential and dipole potential calculated by Adaptive Poisson-Boltzmann Solver with varying numbers n of explicit water layers. (a) Total electrostatic potential with the first 0 (), 1 ( ), 2 (· ), 3 (·  ·), 4 ( · ·), 5 ( ), 6 (· · · ·), 7 () explicit hydration layers included. (b) Membrane dipole potential calculated with varying numbers (n) of explicit water layers; see Eq. 5 for definition of un.

Fig. 4 b depicts the behavior of the membrane dipole potential (u_n) with increasing numbers of explicit water layers. The membrane dipole changes by only 2% upon the addition of the fourth explicit hydration layer to obtain nearly 99% of its asymptotic value (38.3 k_BT/e_c or 993.8 mV). This convergence of the membrane potential at the inclusion four hydration layers is consistent with the results of previous sections and the neutron scattering data described above (Ruffle et al., 2002). As with other explicit solvent molecular dynamics simulations (Zhou and Schulten, 1995; Feller et al., 1996; Essmann and Berkowitz, 1999), the observed membrane potential obtained is larger than experimental values (Flewelling and Hubbel, 1986; Gawrisch et al., 1992). This discrepancy between simulation and experiment is well known and has been discussed in a variety of references (Zhou and Schulten, 1995; Feller et al., 1996; Tobias et al., 1997; Essmann and Berkowitz, 1999; Mashl et al., 2001; Saiz and Klein, 2002). To rectify this discrepancy, several general advances in computational biology would be necessary, including more accurate parameterization of force fields, better water models, inclusion of atomic polarizability, and larger scale simulations to account for the membrane undulation effects. However, this work is aimed at examining the relative contributions of explicit water and this difference does not affect the conclusion that four hydration layers are required for consistent descriptions of the membrane dipole potential in implicit solvent simulations.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

In conclusion, these calculations have shown that water at the membranes surface is highly ordered and has a profound influence on the electrostatic potential of the system. When the whole solvent is treated merely as a continuum dielectric medium, even the sign of the membrane dipole potential will be incorrect. Specifically, implicit solvent electrostatics calculations for membrane systems must include three to four layers of explicit water at the membrane surface to reproduce the membrane dipole potential observed in MD simulations. This observation is consistent both with dipole orientation and fluctuation calculations obtained from the MD data, as well as inelastic neutron scattering experiments (Ruffle et al., 2002).

This work suggests several methods for improving implicit solvent simulations of membrane bilayer. First, it may be possible to incorporate some aspects of the high polarization of interfacial solvent by appropriate modification of the local dielectric constants near the bilayer surface. However, given the linear response assumptions implicit to the Poisson equation, this modification will not correctly model solvent behavior under conditions of dielectric saturation (Beglov and Roux, 1994). Alternatively, the necessary hydration layers could be explicitly included in the model in a hybrid implicit-explicit solvent as previously used in protein dynamics (Lounnas et al., 1999; Im and Roux, 2001) and ion channel (Im et al., 2000). This method offers the advantage of offering a more accurate condition of solvent behavior in the presence of high electric field but suffers from the need to explicit average over solvent coordinates in simulations.

It is clear from these results that although membrane simulations have benefited from the use of implicit solvent methods, simple implementations of such techniques do not accurately describe the electrostatics of the membrane-solvent system. Given the ubiquitous nature of biomembranes in biology, future work must aim to enhance simple continuum electrostatics methods to correctly model the behavior of interfacial solvent and improve agreement with both experiment and explicit solvent simulations.