INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Computational methods that allow a priori estimation of solute-free energies in various bulk solvents or in more complex heterogeneous environments such as biomembranes have many potential applications. Of particular interest to us are computational approaches for predicting drug transport across lipid membranes for application in the rational design and selection of new drug candidates. The semi-empirical functional group-contribution approach has been widely used by a number of investigators to relate partition coefficients, thermodynamic activities, and permeability coefficients to molecular structure with a reasonable degree of success (Davis et al., 1974; Hansch, 1971; Leo et al., 1971; Leo, 1993; Mayer et al., 2000; Xiang and Anderson, 1994b). The group-contribution concept was originally introduced by Langmuir (1925) as the "principle of independent surface action" and subsequently refined and verified by Butler (1937). Implicit in the approach is the assumption that the free energy of solvation is additively composed of independent contributions from the constituent functional groups, though various schemes have also been developed to account for the differences in the properties of isomers and possible couplings between neighboring substituents (Abildskov et al., 1996).

The contribution of a functional group X to the change in standard free energy accompanying the transfer of a molecule from one thermodynamic phase to another can be defined as the difference in the corresponding free energy changes for a substituted (RX) and unsubstituted (RH) compound,

(1)

The second expression highlights the fact that the substituent contribution is commonly obtained from the ratio (F_X) of the partition coefficients of a substituted and parent compound. Providing that the environment adjacent to substituent X is similar from molecule to molecule (e.g., the substituent is well isolated from other substituents on the same molecule with which it might interact), (G°)_X is expected to remain constant for a given X.

We have recently used functional group contributions to the free energy of transfer of small-molecule permeants from water into the rate-determining barrier domain of lipid bilayer membranes generated using relative permeability coefficients of substituted versus unsubstituted small-molecule permeants to determine the location and solvent characteristics of the barrier domain within various model biomembranes (Mayer et al., 2000; Xiang and Anderson, 1994b, 1998, 2000; Xiang et al., 1992, 1998). These experiments have shown that the barrier hydrophobicity in liquid crystalline bilayers (e.g., egg lecithin) resembles that in a relatively nonpolar hydrocarbon solvent (e.g., 1,9-decadiene). The incorporation of a transmembrane peptide (e.g., gramicidin A) into the bilayer or a phase transition to a more rigid gel phase increases barrier domain hydrophilicity (Hiang and Anderson, 1998, 2000_). Unfortunately, the molecular mechanisms responsible for these observations are still poorly understood. For example, the resemblance of the egg lecithin barrier domain solvent properties to decadiene may reflect the presence of double bonds in acyl chains of egg lecithin lipids or the proximity of the barrier domain to the polar interface in the head group region. The increased barrier domain hydrophilicity accompanying the liquid crystallinegel phase transition may be due to changes in properties at a fixed location or a shift of the barrier region toward the more hydrophilic head group region. Similarly, the increase in hydrophilicity accompanying the introduction of a transmembrane peptide may again reflect either a shift in the barrier domain location or the existence of a new transport pathway along the lipid-peptide interface. To fully understand these and other membrane phenomena, free energies of solvation in water and in the barrier domain of lipid bilayers are required.

In recent years, numerous calculations of free energies of hydration have been conducted (Chipot et al., 1994; Worth and Richards, 1994; Ding et al., 1995; Erion and Reddy, 1998; Wang et al., 2001). Biological applications of free energy calculations have mainly focused on protein structures and the binding of small molecules to proteins (Fox et al., 1997; Raghavan et al., 1992; Sugita and Kitao, 1998; Wang et al., 1999). However, studies on solute solvation in lipid membranes have recently been conducted in this and other laboratories (Gambu and Roux, 1997; Jedlovszky and Mezei, 2000; Marrink and Berendsen, 1994, 1996; Pohorille and Wilson, 1996; Xiang and Anderson, 1994a, 1999) mostly using particle-insertion methods or potentials of mean-force calculations. The thermodynamic perturbation, thermodynamic integration, and various combinations of these two techniques (Mezei and Beveridge, 1986), which have been widely used for free energy studies in water, proteins, and other thermodynamic systems, are also suited for free energy studies of solutes in biomembranes, especially for carrying out mutations to obtain relative free energy changes such as functional group contributions.

In this study, we have used a combined all-atom molecular dynamics (MD) simulation, thermodynamic integration, and perturbation approach to calculate contributions of small functional groups to free energies of solvation in water and in a dipalmitoylphosphatidylcholine (DPPC) bilayer membrane. Our focus is on obtaining group contributions to the free energies of transfer from water into the ordered acyl chain region of the lipid membrane to evaluate the hydrophobicity and other related properties within the region that we have previously suggested constitutes the barrier domain with respect to membrane transport. The small changes in solute structure that occur during mutation processes make this type of calculation in lipid membranes feasible. Acquiring individual solvation components in the aqueous and membrane phases has the added advantage that the effects of molecular interactions in these two separate phases, especially the membrane phase, on the transfer free energies can be more confidently assessed. The free energies of solvation for six functional groups (X = -Cl, -CH₃, -OCH₃, -CN, -OH, -COOH) covalently attached to the -methylene of p-toluic acid have been calculated. With the exception of the -methyl-substituted analog (p-ethyl benzoic acid), this series of compounds was used in previous transport experiments to generate experimental estimates of functional group contributions to the free energy of transfer from water into the barrier domain of lipid bilayers (Xiang and Anderson, 1994b). The simulated results can therefore be compared directly to experimental values. This comparative study serves to test the fundamental hypothesis that the ordered acyl chain region is the barrier domain for molecular permeation across biomembranes (Marrink and Berendsen, 1996; Xiang, 1999; Xiang and Anderson, 1994b). The hydrophobicity scale established based on such functional group contribution data may also play a critical role in predicting protein stability in biomembranes (White and Wimley, 1998).

	COMPUTATIONAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Relative free energy calculations

Processes of solute transfer from water into a lipid membrane are very difficult to simulate directly because it takes a prohibitively long time to sample a meaningful amount of configurational space within a lipid membrane or a water phase because most hydrophilic (hydrophobic) compounds prefer to stay within the aqueous phase (lipid membrane). It is also difficult to follow a typical experimental strategy of measuring the free energies of solvation for each chemical entity (e.g., G°, G°, G°, and G° in Fig. 1) as making solvent molecules (e.g., water or lipids) appear around a solute is a very large perturbation of the system. However, to determine relative free energies representing functional group contributions to the transfer free energy, nonphysical but computationally more tractable processes can be designed to form thermodynamic cycles (e.g., upper and lower cycles in Fig. 1) which sum to zero regardless of the path taken because free energy is a state function. Thus, the group contributions to the free energy of solvation in water and a lipid membrane can be obtained from the free energy changes for mutation from RH to RX in the gas phase, water, and the membrane:

(2)

(3)

The free energies of mutation (G°, G°, and G°) in each phase (gas, water, and membrane) were calculated by means of finite difference thermodynamic integration (Mezei, 1987). This method combines aspects of both the perturbation method and thermodynamic integration to improve the convergence and accuracy of relative free energy calculations. In this method, the free energy change, G°^RHRX, due to the mutation RHRX is expressed as an integral,

(4)

where  is a coupling parameter that controls the continuous conversion from RH to RX. In the present studies, this represents a complete mutation from unsubstituted p-toluic acid to -substituted p-toluic acid in a given thermodynamic ensemble. The angled brackets   in the above equation denote the ensemble average at a given  value. The derivative of the free energy function G, G°/, is simulated using the perturbation method that is ideally suited to calculate free energy differences between two states separated by a small (Mezei, 1986),

(5)

where E() is the Hamiltonian of the system at a particular  value. It is interesting to note that, in some previous free energy calculations, was chosen as the  increment, _i  i1, which requires that the  increment must be small enough to ensure that the perturbation method is valid regardless of the simulation length. In this study, a very small (= 0.005) is used, which is independent of the  increment and allows the use of a relatively small number of windows for a complete mutation. Various energy and structural parameters for the mutating atoms in the solute were a linear function of  with these parameters at  = 0 representing those at the initial state (RH) and at  = 1 representing those at the final state (RX). Computing G° (_i) at different  values between  = 0 and  = 1, the total free energy change G°^RHRX can be estimated. The  values at which the free energy derivatives are calculated were chosen on the basis of a Gaussian-Legendre quadrature method (Press et al., 1986),

(6)

where {w(_i)} are the Gaussian-Legendre weight factors, and N is the number of quadrature points. As shown in Fig. 2, a mutation was started from the unsubstituted p-toluic acid and ended when the substituent X grew to a full size. More detailed mutation procedures for the different substituents are reported in Table 1.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   Thermodynamic cycles used in the calculation of substituent contributions to the free energies of transfer from vapor into water, from water into a lipid membrane, and from vapor into a lipid membrane.

View larger version (15K): [in this window] [in a new window]   FIGURE 2   A schematic illustration of the mutation of an unsubstituted p-toluic acid into an -position substituted (X) p-toluic acid. The diagram also shows the atom labels used in Table 3 for assignments of partial atom charges in these solutes.

View this table: [in this window] [in a new window]   TABLE 1   Mutation from p-toluic acid to -substituted (X) p-toluic acid*

Computational details

All MD simulations were conducted using the Discover (v.2.98) program (Molecular Simulations Inc., San Diego, CA) with a CVFF force field. The construction of initial structures and graphic analyses of simulation results were mostly performed with the Insight II program (v.97, Molecular Simulations Inc.). For simulations in water, a flexible, three-site SPC-like water model was used (Lau et al., 1994), in which the internal geometry corresponds to that of the TIP3P model while the intermolecular potential is that of the rigid SPC model. This water model has been characterized by Stouch and coworkers (Alper et al., 1993; Lau et al., 1994) and found to give the structure and dynamics of liquid water in good agreement with both experiment and the computed properties for other commonly used water models. A solute molecule to be mutated was soaked in the center of a cubic box of neat water with a cell length of 30 Å on each side corresponding to a total of approximately 887 water molecules surrounding the solute molecule. A lipid bilayer membrane was constructed by arranging 98 all-trans DPPC molecules into various lattice sites in two monolayers and by inserting a single molecule of -hydroxy-p-toluic acid manually into the bilayer using Insight II. The bilayer system was then energy minimized for 2000 iterations using both steepest descent and conjugate gradient algorithms to remove any high-energy contacts that might have formed during the system-construction process. This was followed by a 100-ps (60 ps at 500 K/40 ps at 323 K) dynamics run with the phosphate atoms fixed at their original positions. A slab of water molecules with the same surface area as the lipid bilayer was then placed on each side of the membrane with a total of 2800 water molecules, after which the system was subjected to a series (20) of 2-ps dynamics runs with the length of the system in the direction of the bilayer normal gradually reduced. The system so built was first subjected to a 120-ps test run for mutation from -hydroxy-p-toluic acid to p-toluic acid. The bilayer assembly generated after this test run was used for subsequent mutations from p-toluic acid to all -substituted p-toluic acid analogs. Another 100-ps dynamic run at 323 K and 1 bar was conducted for each mutation starting from different configurations,

Partial point charges for the series of p-toluic acid analogs were fitted with CHELP (Chirlian and Francl, 1987) to the quantum mechanical electrostatic (ES) potential computed on optimized structures using the Gaussian 98 program (Gaussian, Inc., Carnegie, PA) and ab initio 6-31G* wave functions. Partial charges for atoms in the lipid molecules were adopted from those determined by Stouch et al. using an ES potential-derived method with the 6-31G basis set (Stouch et al., 1991; Williams and Stouch, 1993). These partial charges combined with the CVFF force field, in which an atomic type of cg rather than c was used for carbon atoms in acyl chains, were found to generate energy minima near the experimentally determined crystal structures. Comparison of four commonly used force fields (CVFF, CHARMM22, OPLS, and GROMOS87) for their ability to reproduce experimental free energies of hydration by Helms and Wade (1997) indicated that a combination of the CVFF force field and solute partial charges generated by the ES potential fits gives free energies of hydration for small polar and nonpolar molecules in reasonable agreement with experiment.

Newton's equations of motion for all the atoms in the system were solved using the Verlet leapfrog algorithm (Verlet, 1967) with a time step of 1.0 fs and three-dimensional periodic boundary conditions. The dielectric constant was set at 1.0. Constant temperature (323 K) was maintained using direct velocity scaling (equilibration stage) and a temperature-bath coupling method with a time constant of 0.1 ps (data collection stage) (Berendsen et al., 1984). For simulations in water and lipid membranes, constant pressure (1 bar) was maintained using Berendsen's method of pressure control with a relaxation constant of 0.5 ps (Berendsen et al., 1984). The nonbonded interaction energies were calculated using an 11-Å charge-group-based cutoff and a quintic spline function to smoothly turn off the interaction. A neighbor list was created and updated automatically such that no atoms outside a buffer region (between 10 and 11 Å) could move close enough to interact without detection. Both intra- and intermolecular interactions were explicitly taken into account, which minimizes the bond shrink problem commonly encountered when intragroup contributions are neglected (Pearlman and Kollman, 1991). To start the free energy simulation, the system prepared above was first energy minimized using 1000 iterations each of steepest descents and conjugate gradients. For most simulations, after the minimization steps, the system was equilibrated for 10-50 ps followed by 10-50 ps of data collection for each  value. While simulation runs at each  value were longer than some of the previous free energy studies, a smaller number of six  values was used for each complete mutation because a previous study by Chipot et al. (1996) has shown that it may be preferable to use a limited number of  values and longer sampling times, rather than numerous values with fewer samplings. The overall time for each complete mutation run ranged from 120 to 600 ps. Both forward and backward perturbations were made for each  value, and the free energy derivatives were calculated as the average of the two simulated results. All calculations were performed on SGI workstations located in the authors' laboratory and the College of Pharmacy, University of Utah, and an SGI Origin 2000 and Power Challenger located at the Center for High Performance Computing, University of Utah.

A major obstacle for an accurate evaluation of relative free energy is the poor convergence due to insufficient sampling of system configuration space. This depends strongly on the extent of solute structural changes during the mutation process (Reddy and Erion, 1999). Although the extent of solute structural changes in this study (a replacement of a single functional group) is relatively small, it is still necessary to evaluate the convergence behaviors. In this study, three methods were used to examine the convergence behaviors of the free energy simulations. First, the hysteresis effect, or the difference in the free energy change resulting from forward (+) and backward () perturbations at each particular  value, G( + )/  G(  )/, was evaluated for all the simulations. If sampling or equilibration were poor, or the perturbation too large, large differences would be expected for a perturbation depending on the direction in which it was evaluated. Second, free energy simulations were run with different simulation lengths or starting from different configurations because insufficient sampling of configuration space has been shown to be the major source of poor convergence (Mitchell and McCammon, 1991; Pearlman and Kollman, 1991; Reddy and Erion, 1999). Finally, the variation of the Boltzmann factor, exp(E( ± )/k_BT), over time in the data collection phase was recorded to assess convergence during a given simulation. Systematic drift in the Boltzmann factor is usually symptomatic of a lack of convergence.

The errors (err) in reported free energies, G°^RHRX, were calculated from the standard deviations (SD) of exp(E/k_BT) obtained during the sampling of phase space:

(7)

where the standard deviation in G/_i is expressed as

(8)

and n is the number of sampling points, and var stands for variance.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Free energies of solvation in water

Some bulk water properties obtained by averaging all the simulations in water during the data collection phase are presented in Table 2. The simulated water density at 323 K, 0.987 ± 0.002 g/cm³, is very close to the experimental value of 0.988 g/cm³ (Lide, 1995). The mutations RHRX for six different substituents in Fig. 2 were carried out in water and in vacuum to obtain the relative free energies of solvation for these substituents. The charge distributions for each molecule are listed in Table 3. The charges of the -CH₃ hydrogens in p-toluic acid and p-ethyl benzoic acid were made equivalent by averaging over the three hydrogens because of their ability to interconvert by rotation. Table 3 illustrates that a substituent X attached to the -carbon in p-toluic acid can alter the charge distribution in the rest of the molecule, though this effect is generally small and declines with the distance from the substituent.

Three methods described earlier were used to examine the convergence behaviors. Two runs were made for each mutation process with 10,000 and 50,000 equilibrium and data collection steps, respectively, at each  value. A representative plot of the free energy derivatives, G/_i, versus _i for the mutation RHROH in water calculated forward and backward is shown in Fig. 3 along with the results in the lipid membrane. The differences in the free energy derivatives in water or the lipid membrane and vacuum are shown in Fig. 4 versus . Group contributions to the free energy of solvation were calculated from a weighted summation of these values. More detailed convergence results are presented in Table 4. In most cases, the hysteresis effects and the errors in G/_i are the largest at the beginning of the mutation where, except in the case of the chloro substituent, dummy atoms were transformed into real particles. The simulation errors and the hysteresis effects are generally of a similar order of magnitude and small (<5%) in comparison to the magnitude of the free energy change during a mutation. However, as shown by the results in Tables 4 and 5, the cumulative effects of simulation errors on the uncertainties of the group contributions to free energies of solvation (the differences of the free energy changes in water and vacuum) are much larger (~10-100%). Also evident in Table 4, the statistical errors in the free energy calculations decreased substantially (~2-3-fold) when the simulation length for a complete mutation was increased from 120 to 600 ps, while the hysteresis effects were only modestly diminished. The errors appear to vary inversely with the square root of the simulation length, suggesting that the errors are mostly random in nature.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   Representative plots of the calculated G/ values versus  for the forward and backward mutation of RHROH in water (×, ) and a DPPC membrane (+, ).

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Representative plots of the differences between the calculated G/ values in water and vacuum () and in the lipid membrane and vacuum () versus .

View this table: [in this window] [in a new window]   TABLE 4   Dependence of the convergence behavior of the free energy simulation for the mutation of RH  RX in water on simulation length

If the equilibration time at a given _i value is insufficient, one would expect a systematic drift in the Boltzmann factor, exp(E_i/k_BT), over time during the data-collection stage. The profiles of the Boltzmann factor versus time were generated during every mutation process and saved in a .tbl file for subsequent analyses. Figure 5 shows representative profiles for the mutation RHROH in water at selected  values. Very little drift was observed in nearly every case, even at a shorter simulation length of 120 ps, suggesting that the equilibration time (= 10 ps) at each  value was sufficient. Free energy calculations were also performed starting at different initial configurations for the mutation RHRX. The differences in the simulated free energies starting from three different initial configurations and simulated for 600 ps, respectively, were within 0.5 kcal/mol, which is comparable to the statistical errors. Thus, the simulations appear to have converged.

View larger version (11K): [in this window] [in a new window]   FIGURE 5   Representative profiles of the ensemble-averaged Boltzmann factor, exp(E/kBT), versus time for the mutation of RHROH in water at three different values of i: Solid line, i = 0.03377; short-dashed line, i = 0.3807; and long-dashed line, i = 0.9662.

To further examine the accuracy of the methodology used in the free energy calculation, mutations of methane to methanol in water and vacuum were performed for a total of 600 ps each following the same procedure as that used for the p-toluic acid analogs and compared with experimental results (Hine and Mookerjee, 1975; Ben-Naim and Marcus, 1984). The experimentally observed free energies of solvation for methane and methanol in water are 2.0 and 5.1 kcal/mol, respectively, representing an -OH contribution of 7.1 kcal/mol. Similar -OH contributions to the free energy of solvation in water (6.8 to 7.0 kcal/mol) were found in separate experimental studies using longer chain aliphatic alcohols (Butler, 1937; Diamond and Wright, 1969; Rytting et al., 1978). The value obtained from the above simulation was 6.5 ± 0.3 kcal/mol, which is slightly higher than the experimental values, but not significantly beyond the computational uncertainty. A somewhat higher free energy of solvation (4.6 kcal/mol) for methanol was also found from Monte Carlo simulations by Carlson et al. (1993). These results suggest that the present method for calculating free energies is in reasonable agreement with experiment and other simulation methods.

The functional group contributions to the free energies of transfer from the gaseous state (vacuum) into water are listed in Table 5. The addition of a nonpolar methylene group increases the solvation free energy by 0.5 ± 0.5 kcal/mol, which is comparable to the value reported in the literature of 0.15-0.18 kcal/mol (Davis et al., 1972) within the calculation errors, though the -CH₂- contribution may vary by several hundred cal per mole depending on the hydrocarbon chain length (Davis et al., 1974). Analysis of the energetic components indicates that the change in solute-solvent interaction energy upon mutation RHRCH₃ is a negative value of 1.2 ± 0.5 kcal/mol, more favorable for the hydration of the mutated solute. Reported experimental molar enthalpy changes range from 0.68 to 1.5 kcal/mol (Davis et al., 1974). These results suggest that the entropic contribution is the dominant term in the free energy of solvation of the methylene group in water, though the exact molecular mechanism(s) responsible for this entropy-driven process remains a subject of debate (Grant and Higuchi, 1990). A similar observation was made for the substitution of a chlorine atom where the free energy of hydration is 0.2 kcal/mol. Although considerable data exist for the aqueous solubility of halogenated aromatic compounds, the halogenated aliphatic compounds are poorly represented (Prosser and Davis, 1992). The free energy of solvation in water for an aromatic chlorine is 1.0 kcal/mol based on the study by Tsonopoulos and Prausnitz (1971), and a principal components analysis of Henry's law constants by Suzuki et al. (1992) predicted a group contribution for -Cl bonded to aliphatic compounds of 1.0 kcal/mol. The replacement of an -hydrogen by a chlorine atom increases the ES interaction energy with surrounding water molecules by an average of 0.9 kcal/mol.

The substituents, -CN and -OCH₃, are polar and capable of accepting hydrogen bonds from water. Indeed, the hydration free energies for these groups reported in Table 5 are 4.0 ± 0.5 and 2.5 ± 0.6 kcal/mol, respectively. The literature value for -O- contribution to the free energy of hydration is 3.9 kcal/mol (Davis et al., 1974). If the contributions of -O- and -CH₃ in the group -OCH₃ were additive, the combined free energy of hydration for -OCH₃ would be approximately 3.7 kcal/mol. Substituents capable of donating and accepting hydrogen bonds, such as -OH and -COOH, interact very favorably with water, exhibiting incremental free energies of 5.9 and 7.4 kcal/mol, respectively. Additional 0.5-ns simulations of p-toluic acid analogs in water were performed, and the van der Waals (vdW) and ES interaction components between the solutes and surrounding water molecules were analyzed to better understand these results. Figure 6 shows a bar graph of the average nonbonded interaction energy and its vdW and ES components for -hydroxy-p-toluic acid in water and the differences of these energy components between -hydroxy-p-toluic acid and p-toluic acid. As expected, the dominant contribution favoring solvation upon replacement of an -H atom by an -OH group is the electrostatic interaction between the -OH group and water. Interestingly, the vdW interaction energy actually increases with the addition of an -OH. This suggests that the formation of an oriented hydrogen bond between the -OH and water induces reorganization of the solute and neighboring water molecules such that the vdW interaction is somewhat less favorable.

View larger version (22K): [in this window] [in a new window]   FIGURE 6   (A) The average nonbonded interaction energy and its vdW and ES components for -hydroxy-p-toluic acid in water. (B) The differences of these energy components between -hydroxy-p-toluic acid and unsubstituted p-toluic acid.

In general, the simulated free energy changes for relatively hydrophilic substituents (e.g., -OCH₃, -CN, -OH, and -COOH) are ~0.9-1.9 kcal/mol larger than the experimental values, though they are closer to those obtained indirectly from a principal components analysis of Henry's law constants (Suzuki et al., 1992). The larger free energies obtained may result from the closeness of the aromatic ring to the substituent X at the  position, which may affect the charge distribution and thereby the free energy of hydration. Indeed, as shown in Table 3, different substituents have a noticeable effect on the partial charges on the aromatic ring and the rest of the molecule, though these changes in partial charge appear to decrease with distance from the substituent and are barely detectable for the carboxylic group on the other end of the molecule. This hypothesis is further supported by octanol/water partition experiments in which Iwasa et al. (1965) found that the  constants for polar substituents X obtained from the series of C₆H₅CH₂X are generally larger by 0.0-0.30 units (representing 0-0.42 kcal/mol in free energy) than those from the series of C₆H₅CH₂CH₂X.

Free energies of solvation in a DPPC lipid bilayer

Previous studies of solute solvation in lipid membranes reported by this and other laboratories have used particle insertion and potential of mean-force methods (Marrink and Berendsen, 1994, 1996; Pohorille and Wilson, 1996; Gambu and Roux, 1997; Xiang and Anderson, 1999; Jedlovszky and Mezei, 2000). The thermodynamic perturbation, thermodynamic integration, and various combinations of these two techniques (Mezei and Beveridge, 1986), which have been widely used for free energy studies in water, proteins, and other thermodynamic systems, are also suited for free energy simulations in biomembranes, though to the best of our knowledge, there appears to be no reported work in lipid membranes. The need for enormous computational resources may be partly responsible for this.

Figure 7 shows a snapshot of a hydrated DPPC membrane within which a single -hydroxy-p-toluic acid molecule resides. An additional 500-ps run was conducted using this system to evaluate the state of equilibration. The time evolution for the average interaction energy, E_int, of the intercalated solute, -hydroxy-p-toluic acid, in the bilayer interior after the equilibration of the bilayer system is presented in Fig. 8. The interaction energy fluctuated around its initial value and no significant drift was found. This suggests that the solute saw a thermodynamically constant environment during subsequent simulations after the initial equilibration step.

View larger version (115K): [in this window] [in a new window]   FIGURE 7   A snapshot of a hydrated DPPC membrane within the ordered acyl chain region of which an hydroxy- p-toluic acid molecule resides.

View larger version (16K): [in this window] [in a new window]   FIGURE 8   The time evolution for the average interaction energy, Eint, of the intercalated solute, -hydroxy-p-toluic acid, in the bilayer interior after the equilibration of the bilayer system.

Some structural properties of the simulated membrane, which were averaged over all the membrane simulations, are reported in Table 2. The lipid membrane was simulated at a constant pressure of 1 bar along all three Cartesian directions. Although a constant surface tension ensemble has been suggested as a more relevant model (Chiu et al., 1995; Feller et al., 1995), Tieleman and Berendsen (1996) have demonstrated that the use of constant-pressure or surface-tension ensemble makes little practical difference in a simulation. The maintenance of a proper surface density appears to be more important. The surface density of the DPPC membrane simulated in this study varied within a narrow range of 60-64 Ų per lipid with an average of 61.5 Å² per lipid. Experimental studies of similar lipid membranes in a liquid-crystalline phase give values in the range of 57-72 Ų per lipid depending on the technique used in the measurement. A commonly accepted value is somewhere between 62 and 64 Å² per lipid (Nagle et al., 1996). Table 2 also shows that the vector between the phosphorus and nitrogen atoms in the bilayer is, on the average, nearly parallel to the membrane surface with a small tilt toward the adjacent water phase, consistent with experiments and other MD simulations (Smondyrev and Berkowitz, 1999). The number of gauche defects per chain (last 12 C-C bonds) was found to be 2.8, which is significantly lower than the value 3.7 obtained from a statistical mechanical analysis of the experimental deuterium order parameters (Schindler and Seelig, 1975). A similar result (2.6) was obtained in the MD simulation of a liquid-crystalline DPPC bilayer by Tu et al. (1995). Although the exact cause for the discrepancy is still not clear, it may arise partly from the simplified treatment of the isomerization and overall rotation in the statistical mechanical analysis using the Marcelja model (Xiang and Anderson, 1995).

Membrane structure can be further characterized by the order parameter, S, for the orientational angle _i of the plane formed by the ith -CH₂- segment in the lipid molecule with respect to normal to the membrane interface:

(9)

The segmental order parameter S_mol is commonly defined as (Seelig and Niederberger, 1974)

(10)

The present simulation was conducted using the same lipid (DPPC) and the same temperature (323 K) as used in an NMR experiment reported by Seelig and Seelig (1974). It is therefore reasonable to compare the order parameters (S_mol) obtained in the present study with their values. This comparison is shown in Fig. 9. The order parameters obtained here generally agree with the experimental values, especially in the ordered middle chain region where the solute of interest resides, though a somewhat larger order parameter was found at the end of the chains. The presence of a foreign molecule in the membrane interior may disturb the organization of lipids in the membrane, but this can only be observed when solute concentration is much higher than that used in the present study. Indeed, studies by Cascales, Bassolino-Klimas and coworkers (Alper et al., 1993; Cascales et al., 1998) have shown that, with lipid/solute ratios of 9:1 and 24:1, respectively, which are much smaller than the ratio used in this study (98:1), only a slight decrease in the order of the membranes was observed, and it occurred only in the deepest zones of the hydrocarbon region (Cascales et al., 1998).

View larger version (14K): [in this window] [in a new window]   FIGURE 9   The order parameters (2SCD), averaged over both acyl chains, versus the segment number of the lipid acyl chains. , Present study; +, deuterated NMR experiment by Seelig and Seelig (1974).

A lipid membrane is heterogeneous along the direction normal to the bilayer interface. Abundant evidence suggests that, among various distinctive regions (cf., the hydrated head groups, the glycerol linkages, the ordered acyl chains, and the relatively disordered region near the center of the membrane), the ordered acyl chain region is the least favorable for partitioning of polar molecules, and, therefore, this region constitutes the barrier domain for transport of polar solute molecules (Marqusee and Dill, 1986; Marrink and Berendsen, 1994; Xiang and Anderson, 1994a). Various thermodynamic equilibrium methods for measuring partition coefficients of solutes into biomembranes have limited success in probing hydrophobicity properties of the barrier domain because solute molecules partitioned into lipid membranes tend to reside near a region of minimum free energy (e.g., near the head group region for polar solutes [Diamond and Katz, 1974]) rather than in the barrier domain where the free energy of solvation is maximum. Because of these difficulties, organic solvents, in particular octanol, hexadecane, and decadiene, have been used to model various biophases.

To obtain free energies of solvation within this ordered acyl-chain region from the present simulations, the targeted solute was restrained to move within a small region such that the position of the -substituent X was centered ~8-9 Å away from the average position of the glycerol groups in the same monolayer by imposing a harmonic force (50 kcal/Ų) on the aromatic carbon adjacent to the carboxylic acid group in substituted and unsubstituted p-toluic acids. Because the thickness of the entire acyl chain region is 32 ± 3 Å, about 70% of which is highly ordered as demonstrated by the order parameter profile in Fig. 9, the -substituent X in the solute was located in the ordered-chain region during the free energy calculations. Although fixing the position of a targeted solute significantly altered free energy results (~1 kcal/mol), restraining it to move within a small region yielded results almost identical to those obtained in simulations without constraints.

Simulations of adequately hydrated lipid bilayers require much larger system assemblies in comparison to simulations in water alone. As a consequence, the CPU time per time step is much longer than that for simulations in water. Fortunately, the dominant force governing solute motion within the acyl-chain region of a lipid membrane is the vdW interaction (see below), which leads to faster convergence than when both strong ES and vdW interactions are present. Thus, simulation times (120 and 360 ps) equal to or shorter than those for simulations in water were used for each complete mutation. Sampling a sufficient rotational space is a prerequisite for free energy convergence. To this end, we have calculated the rotational correlation times for the solute, -hydroxy-p-toluic acid, using the vector between two carbon atoms at opposite sides of the two substituents covalently bonded to the benzene ring and using the MD trajectories obtained. Figure 10 shows the autocorrelation functions, C(t), in water and the DPPC bilayer at 323 K versus time. The rotational correlation times estimated from fits to a single exponential function are 7.5 ± 0.3 ps (r² = 0.90) and 21.9 ± 0.4 ps (r² = 0.95) in water and the DPPC bilayer, respectively, suggesting that the simulation runs of 120-600 ps are adequate to sample sufficient conformational space. The correlation time obtained in the bilayer interior is on the same order of magnitude as the ones for similar-sized solutes such as benzene within different regions of a DMPC bilayer at 320 K (25 ps) obtained via an MD simulation by Bassolino-Klimas et al. (1993), benzyl alcohol in the ordered acyl-chain region and near the center of a DMPC bilayer at 325 K (3-8 ps) via an MD simulation by Cascales et al. (1998), and di-tert-butyl nitroxide in DPPC bilayers at 323 K (20 ps) obtained experimentally (Dix et al., 1978).

View larger version (27K): [in this window] [in a new window]   FIGURE 10   The autocorrelation functions, C(t), for rotation of solute, -hydroxy-p-toluic acid, in water and the DPPC bilayer at 323 K versus time.

The convergence behaviors were examined by the same techniques as those used for mutations in water. The results for all the substituents are summarized in Table 6. A representative plot of the free energy derivatives, G/_i, versus the _i value for the mutation RHROH in the DPPC bilayer membrane (Fig. 3) demonstrates that, as in water, the hysteresis effects for forward and backward perturbations were generally small in comparison to the free energy derivatives during the mutation. However, as shown in Fig. 4, differences in the free energy derivatives in the membrane and vacuum, from which the group contributions to free energy of solvation were obtained, are much smaller. Consequently, the cumulative effects of simulation errors on the uncertainties of the group contributions to free energies of solvation are much larger (~25-200%). Increasing the simulation length from 120 to 360 ps significantly improved the statistical precision of the simulated free energy changes, though, as in water, the hysteresis effects changed only modestly. Changes of the Boltzmann factor, exp(E_i/k_BT), over time at different _i values, some of which are plotted in Fig. 11 for RHROH in the membrane, were nearly flat with no significant upward or downward drift for both short and long simulation times at each _i value, indicating that the equilibration time is sufficient. Simulations of RHROH lasting for 300 ps each were also performed starting at three different initial configurations. The simulated free energies were within 0.8 kcal/mol, which is somewhat larger than similar simulations in water but still on the same order of magnitude as the statistical errors. Combining all these results, it appears that the simulations had converged.

View this table: [in this window] [in a new window]   TABLE 6   Dependence of the convergence behavior of the free energy simulation for the mutation RH  RX in a DPPC lipid bilayer on simulation length

View larger version (11K): [in this window] [in a new window]   FIGURE 11   Representative profiles of the ensemble-averaged Boltzmann factor, exp(E/kBT), versus time for the mutation of RHROH in a DPPC lipid bilayer at three different values of i. Solid line, i = 0.03377; short-dashed line, i = 0.3807; long-dashed line, i = 0.9662.

Functional group contributions to the free energies of transfer from a gas phase into the ordered acyl-chain region of the DPPC bilayer membrane are listed in Table 7. All of the group contributions are negative, varying within a relatively narrow range from 0.3 to 2.3 kcal/mol for different substituents. The calculated -CH₂- contribution to the solvation free energy, 0.9 ± 0.6 kcal/mol, is close to the values of 0.87 and 0.70 kcal/mol obtained from the solubilities of n-alkane gases in egg lecithin membranes (Miller et al., 1977) and in a nonpolar hydrocarbon solvent, respectively. A similar free energy contribution was found for the -Cl substituent (0.3 kcal/mol).

Because few literature values are available for the solvation of more hydrophilic substituents such as -OH, -COOH, -OCH₃, and -CN in the interior of bilayer membranes, values calculated from the above simulations are also compared in Table 7 to data available for the transfer of solutes from vapor into a nonpolar hydrocarbon solvent. The -OH value obtained in this study (1.7 ± 0.4 kcal/mol) agrees well with that obtained experimentally for solvation in nonpolar hydrocarbons (1.3 to 1.6 kcal/mol) but is significantly larger than the value (4.4 kcal/mol) obtained for the transfer of n-alcohols from vapor to hydrated DMPC membranes (Diamond and Katz, 1974; Diamond et al., 1974). This discrepancy undoubtedly reflects the preference of n-alcohols for the membrane-water interface where the ES interaction for the -OH group is more favorable. This may also explain why the contribution of the -CH₂- group to solvation free energies of n-alcohols (0.41 kcal/mol) is significantly smaller than that obtained for the series of n-alkanes in egg PC/PA bilayers (0.87 kcal/mol) (Miller et al., 1977) and in a nonpolar hydrocarbon solvent (0.70 kcal/mol) (Pescar and Martin, 1966). Indeed, these results may reflect the loss of degrees of freedom on anchoring the alcohols at the more rigidly structured head-group interface and the fact that the bilayer interior may be a better solvent for hydrocarbons than is the head-group region. An analysis of solute-membrane interaction energies (see below) indicated that the vdW interaction energy for solutes in the membrane increases upon the substitution of a functional group by an amount ranging from 0.8 kcal/mol for -Cl to 0.5 kcal/mol for -OCH₃.

Free energies of transfer between water and a DPPC bilayer

The hydrophobicity within the barrier domain region of lipid membranes is an important property governing molecular transport rates across biomembranes. Membrane hydrophobicity may also play a role in regulating the activity of membrane proteins and other molecules of biological interest. Numerous models for membrane hydrophobicity have been proposed, most of which are based on free energies of transfer for amino acids and peptide units from water to the biophase of interest and their correlation with those from water to other better known bulk solvents (Nozaki and Tanford, 1971; White and Wimley, 1998). A similar approach has been taken in this study, in which the use of small functional groups as probes has the added advantage that the effects of probe size on solute solubility are minimal and more rapid convergence in simulations can be expected.

Group contributions to the standard free energies of transfer from water into the membrane interior, (G°_wm)_X = G°  G°, can be obtained from the substituent contributions to free energies of solvation in water (G°) and in the DPPC bilayer membrane (G°) as presented in the preceding sections. The results are compiled in Table 8 along with literature values for transport across eggPC bilayers and for partitioning from water into various organic solvents for comparison. It is evident that group contributions to the free energies of transfer from water to a less polar solvent are very sensitive to the solvent's properties, especially for polar, hydrogen-bonding substituents such as -OH and -COOH. This sensitivity can therefore be utilized to determine hydrophobicity properties in lipid membranes. In particular, the following linear free energy relationship originally developed by Collander (1954),

(11)

can be used to assess the similarity of the two thermodynamic systems, s₁ and s₂, where (G°) is the group contribution to the standard free energy of transfer from water into another thermodynamic system s₁. Although a good correlation is an indication of solvent similarity, a more quantitative measure is the slope s, often referred to as the selectivity constant, or hydrophobicity scale. An s value of >1 implies a system s₁ more hydrophobic, and therefore more selective, than the system s₂, and vice versa for s < 1. An s value of ~1 suggests similar hydrophobicity for the two systems. Because most of the group contribution data (cf., Table 8) in the following correlation analyses were obtained in previous transport and partitioning experiments using the same series of -substituted p-toluic acid analogs (Xiang and Anderson, 1994b; Xiang et al., 1998), any effect of the aromatic ring in these solutes on the group contributions would be expected to cancel.