INTRODUCTION

Top Abstract Introduction Methods Results Discussion Conclusion References

Predicting the tertiary structure of membrane proteins remains a challenging problem. Some insight into the structure of helical membrane proteins has been achieved by using two types of analysis based on hydrophobicity. One type of analysis, as exemplified by the Goldman-Engelman-Steitz (GES) scale, attempts to define a relative hydrophobicity for each amino acid that can then be combined within a sliding window to determine the likelihood of a particular stretch of amino acids crossing the membrane bilayer (Engelman et al., 1986). Other hydrophobicity scales have been proposed, and applied with varying success, to membrane proteins (e.g., Segrest and Feldmann, 1974; von Heijne, 1981; Argos et al., 1982; Kyte and Doolittle, 1982; Eisenberg et al., 1984). A second type of analysis extends hydrophobicity arguments to predict the orientation of helices by using hydrophobic moments (Eisenberg, 1984; Eisenberg et al., 1984). The hydrophobic moments are based on the idea of using residue hydrophobicity with a vector sense to estimate the net direction and magnitude of -helix amphiphilicity, and thus predict which portions are most likely to be lipid facing.

The GES scale was explicitly developed for application to membrane proteins (Engelman et al., 1986). Whereas most hydrophobicity scales are based largely on relative partitioning of side-chain analogs or amino acid residues between solvents, the GES scale took into account helical secondary structure as well as hydrophobic and hydrophilic interactions (Engelman et al., 1986). The hydrophobic term is derived by assuming that the effective surface area of a particular side chain within an -helix is proportional to the free energy of transfer. The hydrophilic term is composed of two types of interactions. The first contribution includes the effect of electrostatics through the Born model. The second estimates hydrogen bonding interactions that contribute to the transfer process. The result, frequently cited and used, assigns a free energy of transfer for individual amino acids within an -helical structure. The assigned numbers are used within a sliding window context, where 20 amino acids are taken as a group, and the likelihood of their forming a membrane-spanning region crossing the bilayer is calculated from the sum of their GES values. The conceptual basis behind this scale is the two-stage model of membrane protein folding (Engelman and Steitz, 1981; Popot and Engelman, 1990). Recently this type of concept has been improved by the application of a neural net trained on a set of membrane proteins and with information from multiple sequence alignment (Rost et al., 1995). The neural net achieved 95% accuracy at predicting membrane-spanning regions in a set of -helical membrane proteins (Rost et al., 1995).

The concept of hydrophobic moments was developed to quantify the nonuniformity of hydrophobic residue distributions within the helix, in the hope that this would aid prediction of stable structures (Eisenberg, 1984). The first application of the concept was to differentiate surface-associated -helices from bilayer-crossing helices. Several modeling groups have tried to combine this information into a form that can be used for structure prediction. An example of this is the work of Taylor et al. (1994), where helices were initially positioned based on multiple sequence alignment and hydrophobic moment. The additional assumption was made that a greater mutation rate is allowed in the lipid-exposed helical faces relative to helix-helix interaction faces.

Molecular dynamics simulations of neat lipid and lipid-protein systems remain a challenging area. Recent progress in these simulations has been summarized in a book and two review articles (Merz and Roux, 1996; Pastor, 1994; Merz, 1997). The largest set of work has involved simulations of neat bilayer systems. Simulations of lipid-protein systems require potential functions and simulation conditions able to emulate the neat bilayer environment (see Shen et al., 1997, for a recent simulation). The debate on the most appropriate ensemble for simulating bilayers is still not fully resolved (Chiu et al., 1995; Feller et al., 1995; Tu et al., 1995; Jahnig, 1996; Feller and Pastor, 1996; Tielman and Berendsen, 1996). Despite the evolving nature of this field, it is clear that simulations have progressed to the point where significant insight into the molecular details of bilayer systems has been achieved. Comparisons between experiment and simulations, where possible, have shown that the calculations are in agreement with measured quantities and can provide additional insight into molecular details that are not available experimentally (e.g. Venable et al., 1993; Woolf and Roux, 1996).

The analysis of molecular dynamics trajectories frequently involves consideration of deviations from the starting structure and determination of dynamic and average structural properties. Analysis of interaction energies is less commonly performed, although the approach was used in the early development of the molecular dynamics field (Brooks et al., 1988). Recently the approach was used to analyze the energetic details of gramicidin A-dimyristoylphosphatidylcholine interactions (Woolf and Roux, 1996). The method may also be used, by separate consideration of van der Waal (vdW) and electrostatic interactions, to provide an estimate of relative free energies by linear response techniques (Aqvist et al., 1994). Although the calculation of interaction energies can provide insights into the strength of connections between components of a mixed system, there are certain limitations to the method that should be made clear. First, although the interaction energies provide useful information, they are not being used, in this paper, to probe the thermodynamic measurables of a system. They are not equivalent to a relative free energy calculation using thermodynamic perturbation theory (e.g., Kollman, 1993). Second, the calculation of an interaction energy distribution function has the same sampling problems common to all molecular dynamics calculations. If sampling is inadequate, the distribution functions may not be fully converged. Despite these limitations, the interaction energy distributions can give insight into important details of a molecular dynamics system. For example, the distribution of the interaction energy between helix and dimyristoylphosphatidylcholine (DMPC) provides an estimate of the likelihood of observing tightly interacting lipids, and such lipids may be important to structure and function.

The current paper should be viewed as a first step toward a detailed understanding of the preferred environment for -helices as a function of amino acid sequence. The analysis that follows uses the 10 simulations of individual -helices of bacteriorhodopsin in explicit DMPC (Woolf, 1997). The trajectories are used to calculate the interaction energy probability distributions between subsets of atoms in the system. The results provide new insights into the microcopic details of -helices spanning membrane bilayers and comments on both hydrophobicity scales and hydrophobic moments.

	METHODS

Top Abstract Introduction Methods Results Discussion Conclusion References

The average structural and dynamic features of the 10 simulations analyzed in this paper were described previously (Woolf, 1997). The program CHARMm (version 23) was used for all calculations. The previous and related papers should be consulted for a more detailed description of the methods used for starting the trajectory and producing the full simulation (Woolf and Roux, 1994, 1996; Belohorcov et al., 1997; Woolf, 1997). For completeness, the methods are reviewed below.

Initial structures

The Henderson bacteriorhodopsin structure (Henderson et al., 1990; pdb number 1BRD) was used as the starting point for each -helix. The -helices were extracted from the Brookhaven entry, and each helix was aligned with its axis normal to the bilayer. The more recently refined bacteriorhodopsin structure (Grigorieff et al., 1996; 2BRD) was not available at the time the simulations of this paper were completed (Fall, 1995).

System construction and dynamics

The simulations used the NVE (constant number of atoms, volume, energy; microcanonical ensemble). This choice of simulation system required an estimate of the lateral square area for each -helix and DMPC lipid to set the pressure. To begin the current simulations, initial placement of the lipids was determined by a 500-ps simulation of large vdW spheres to simulate the DMPC headgroups. The -helix was kept rigid. The spheres were confined to planes surrounding the phosphate region of the eventual bilayer. This created a set of initial placements for the lipids responsive to the distribution of side chains and shape of the -helix. The set of 12 initial positions was then used for placement of the DMPC phosphate atoms. A set of 2000 DMPC lipid conformations provided by Pastor and co-workers was used to initiate the simulations from conformers representative of the liquid crystalline state (Hardy and Pastor, 1994). This leads to a more rapid relaxation than would a start from an all-trans conformation (e.g., Heller et al., 1993).

The lipids were systematically rotated and translated to reduce the number of bad contacts. A bad contact was heavy atoms within 2.6 Å of one another. After the systematic rotations and translations, a gradual awareness of neighboring atoms was included in the system by using a series of short nonbonded cut-offs. Equilibration dynamics was pursued with a gradually decreasing series of harmonic restraints. The first 25 ps of dynamics used a 1.0 kcal/mol-Ų restraint on all heavy atoms and Langevin dynamics with a collision frequency of 5/ps¹. This low value for the collision frequency has been shown to maximally allow for relaxation of the system while providing a temperature coupling to the heat bath (Loncharich et al., 1992). The next 25 ps used velocity scaling, with no restraints on atom positions. This was followed by 50 ps with no restraints and no temperature coupling. The temperature remained constant throughout the simulation of each system.

The nonbonded interactions were cut at 12 Å with electrostatics shifted and v dW interactions switched. The vdW interactions used a switching region of 4 Å. The list was generated to 14 Å, and an automatic update of the nonbonded interactions was performed if the deviations exceeded 1 Å. Atom-based cutoffs were used. Previous simulations have shown that these choices produce reasonable results with the CHARMm potential function (e.g., Venable et al., 1993; Woolf and Roux, 1994, 1996). The empirical parameters developed in the Karplus group were used (Schlenkrich et al., 1996).

Conformations were saved every 50 fs throughout the time from 100 ps to 350 ps for eight simulations. Two simulations were extended from 350 ps to 600 ps. In all cases, the SHAKE algorithm was used with a time step of 2 fs. The simulations, on average, required 1.8 CPU h/ps on a SGI R4400 machine.

Because three helices contained charged residues that might be within the hydrophobic part of the bilayer, additional simulations were performed to address the differences with the presence of charged side chains at those locations. In particular, this meant that the lysine residues at positions 3 and 4 in helix B, the aspartate residues at positions 6 and 17 in helix C, and the aspartate residue at position 8 in helix D were simulated in both charged and neutral form.

Analysis

Interaction energies were calculated between the environment and either the full helix, a single residue, or a helical stripe in each of the 10 simulations. The vdW, electrostatic, and total interaction energies were determined. The approach was used in slightly different ways for the full helix and the residue/stripe interactions.

The full helix interaction energies were calculated between individual DMPC molecules and the entire helix. A nonbonded cutoff of 100 Å was used to verify that the results were not biased by missing events due to a particular cutoff distance. For each conformation of the trajectory, the calculated total energy was used for assignment to a bin. The bin width was 1.0 kcal/mol, and the number of events in each bin was normalized by the number of conformations used for the calculation. This results in a distribution function that when integrated returns the total number of lipids (including images) surrounding the helix (84). Thus a particular entry from the distribution provides an estimate of the average number of lipid molecules expected to have that interaction in a typical equilibrium conformation. The full trajectory was divided into five subsets to estimate confidence bars for each point of the distribution function (e.g., Loncharich et al., 1992). The approach was to calculate the mean for each of the five subtrajectories, and then from the set of five means, to calculate a system mean and standard deviation. Confidence bars were then assigned by comparison to Student's t-test; these represent the confidence that values within the mean ± the 80% confidence bar are correct representations of the true distribution.

The helical stripes approach determined the interaction of 45° slices of the helix, thought of as a cylinder with a normal along z, with the environment. The interactions were computed for both lipid and water. A further division into the electrostatic and vdW components was made. For each conformation of the trajectory, the interaction energy was computed for the ~60 atoms in a given slice. The next set was determined with a 5° rotation and so forth until the full 360° was computed. This was continued for the full trajectory. A similar analysis was made of the contributions to the interaction total in each slice from subsets of atoms belonging to each residue. This second calculation was used for the contour plots that describe the contributions from each residue to the full interaction energy within a given angular bin.

The interaction energy between residues or stripes and the environment was calculated with images and the same cutoff options as used in the simulation. The individual residue selections included the backbone as well as side chain. The collection of interaction energies was then binned (2.5 kcal/mol increment) and normalized by the number of conformations to produce the interaction energy probability density. From this density, the mean and standard deviation were calculated. This effectively assumes that the distribution is well described by a Gaussian. Although this choice of mean and standard deviation to represent the probability distribution was adequate for many of the residues, it did not fit all of the residues. Nonetheless, the error in the presentation of the interaction probability density functions was deemed small enough for present purposes to support the analysis.

The contour plots (prepared with Mathematica; Wolfram Research) present the contributions of particular residues to the total interaction energy for a helical stripe. The 10 contour lines are drawn with shading such that a darker filling is a stronger interaction energy. The contribution of a particular residue to a given bin angle included only those parts of the residue that had been selected as being within the pie shaped wedge. Thus, if only one atom of, e.g., a leucine was selected as being within the wedge being calculated, then only the one atom's contribution from the residue would be included in the contour plot.

	RESULTS

Top Abstract Introduction Methods Results Discussion Conclusion References

Ten molecular dynamics simulations of individual -helices from bacteriorhodopsin in explicit DMPC bilayers are analyzed in terms of interaction energy distributions. These distributions are for the interaction between all atoms of the helix and individual DMPC molecules, for helical stripes of helix atoms interacting with DMPC, and for individual amino acids interacting with DMPC or water. The previous calculations (Woolf, 1997) showed intriguing differences in average structural and dynamic properties between the 10 simulations. In particular, differences in RMS deviations from the average structure varied with the sequence of amino acids along the helix. Particularly large deviations for helix D were observed and rationalized by the lack of aromatic residues interacting with the interface region. Helix A had the smallest RMS deviations throughout the trajectory and had the most complete canonical hydrogen bonding pattern for an -helix. Analysis of cylinders fit to two turns of the helix showed that two of the three proline-containing helices may have fluctuations about the proline axis. These differences emerged in the trajectories, despite the use of identical simulation conditions. This paper continues the analysis of the 10 trajectories by examining interaction energies.

Full helix interaction energies

The interaction energy between each helix and individual lipid molecules was calculated over the course of the trajectories. A particularly strong interaction of one helix relative to another would indicate that the strongly interacting helix is more likely to prefer a lipid environment than an environment with helix-helix interactions or water. The interaction energies were binned and the full distribution determined by dividing by the number of conformations. In each case, the full trajectory was divided into five subtrajectories for analysis, so that confidence bars for the mean values could be calculated (e.g. Loncharich et al., 1992). The resulting distribution suggests the average number of lipid molecules interacting with a particular interaction energy at any point in the simulation. These results are shown in Figs. 1 and 2 and in Table 1. There are significant similarities as well as differences between each of the 10 simulations. First, each simulation has a large population of weakly interacting lipids. These are not due to the cutoff used in the simulations. This was checked by calculations of the distributions both with a 12-Å cutoff and with a 100-Å cutoff. The results were similar in overall shape, with only moderate differences (not shown). All of the calculations presented in Figs. 1 and 2 and Table 1 used the longer cutoff.

View larger version (33K): [in this window] [in a new window]   FIGURE 1   Interaction energy distributions between individual DMPC molecules and the full helix for all 10 simulations. The curves are an estimate for the probability of finding a particular number of lipids surrounding the central helix with a particular interaction energy. They have been obtained by computing the interaction energy over the trajectory and binning and normalizing the resulting series of numbers.

View larger version (35K): [in this window] [in a new window]   FIGURE 2   Contributions to the total interaction energy from vdW and electrostatic components. The x axis represents the full interaction energy, and the y axis is the value of the vdW or electrostatic interaction energy within the bins. Thus a total interaction energy of 50 kcal/mol might have 25 kcal/mol of electrostatic energy (shown as darker line) and 25 kcal/mol of vdW energy (shown as a thinner line). Confidence bars (80%) are computed from five subsets of the full trajectory.

View this table: [in this window] [in a new window]   TABLE 1   Average number of DMPC lipids with a particular interaction energy for the indicated -helix

Several differences between the systems are interesting. For example, some systems had significant probabilities of very strong interaction energies. Helix G, for example, had interactions as high as 100 kcal/mol. Others, such as helix B-neutral, did not see interaction energies beyond 40 kcal/mol. Some helices had very clear second maxima within the region from 10 to 20 kcal/mol. For example, helix A and helix G have secondary maxima around 15 kcal/mol. A few helices show relatively smooth decreases in probability from the initial peak, such as helices C-neutral and E.

Fig. 2 illustrates the contribution to each energetic bin of Fig. 1 from vdW and electrostatic components. The vdW contributions are dominated by the alkane-chain-protein interactions, whereas the electrostatic interactions are largely due to interactions with the headgroup. For all 10 simulations, the vdW contributions dominated the interaction energy for small to medium interaction energies. However, for the strongest interactions, the electrostatic contribution was at least as strong, if not the dominant contribution to the total interaction energy. Confidence bars are plotted with each data point and emphasize that the statistics are poor for the strongest interactions.

Table 1 contains the numeric values for the distribution functions of Fig. 1 at selected values of the x axis. For example, the average number of DMPC lipid molecules that interact with 1.0 ± 0.5 kcal/mol with the helix varies from a high of 31.3 to a low of 14.1 for helices E and C-charged. A strong interaction, for example, the average number of lipid molecules with 30.0 kcal/mol interaction energy, varies from a high of 0.22 for helix A to a low of 0.03 for helix G.

These results are important for general considerations of protein-lipid interactions. For example, the mean-field theories that have been constructed for protein-lipid interactions have concentrated on the alkane-chain-protein type of interaction to the exclusion of the electrostatic interactions between the headgroup and the rest of the helix (e.g. Owicki et al., 1978; Zuckermann and Pink, 1980). The results of this analysis suggest that there is a large range of interaction energies with a mix of energetic type. For weak to moderate interaction energies (roughly +5 to 30 kcal/mol), the vdW interaction is the dominant type of interaction energy. For strong interactions, the electrostatic interaction energy is the dominant form.

Individual residue interactions

Most hydrophobicity scales treat the membrane as a slab of bulk hydrocarbon, surrounded by bulk water. Their application to atomic-level questions is limited by the neglect of the interface between lipid and water, which may extend over a 10-15-Å range (White, 1994). Although the present calculations do not reveal the full relative free energy for movement of amino acid residues from bulk water to bilayer, they do address some of the same issues. A strong interaction energy indicates that the particular amino acid is well situated within the bilayer environment, whereas a weak interaction energy indicates a less stable amino acid-environment interaction. Two presentations for amino acid-environment interactions are shown. The first (Fig. 3) shows the interaction energies for each helix along the length of the helix. The second (Figs. 4-7) pools the residues from all 10 simulations into an analysis by residue type.

View larger version (45K): [in this window] [in a new window]   FIGURE 3   Individual residue interaction energies for all 10 simulations. The interaction energy between a particular amino acid and DMPC was calculated, binned, and normalized. The data points represent the mean value of the distribution, and the lines the rms deviations from the mean. Two y axes are used. The left-hand y axis refers to the filled circles for the noncharged interaction energies. The right-hand y axis refers to the open squares for the charged residue interactions. In each case the residue number is along the x axis.

View larger version (34K): [in this window] [in a new window]   FIGURE 4   Individual residue interaction energies for nonpolar amino acids. The range of interaction energies observed for each amino acid is indicated as a function of relative position in the bilayer. The circles are the mean and the lines from the circles an indication of the rms deviation. The number of residues considered and the average over all residues is also shown. The left-hand column is for interactions with DMPC; the right-hand column is for interactions with water.

View larger version (31K): [in this window] [in a new window]   FIGURE 5   A plot similar to that in Fig. 4, with the aromatic amino acids.

View larger version (26K): [in this window] [in a new window]   FIGURE 6   A plot similar to that in Fig. 4, for the charged amino acids.

View larger version (35K): [in this window] [in a new window]   FIGURE 7   A plot similar to that in Fig. 4, for the remaining amino acids.

Fig. 3 shows the mean and rms interaction energies for residues along the length of all 10 helices. A general result across all 10 helices is smaller interaction energies and rms fluctuations for interaction energies within the central regions of the bilayer relative to those regions in the interface. This trend is clearest for helices E, C, and B. The presence of charged residues led to strong coupling to the environment, which dominated the other interactions. These charge interactions are indicated by the open squares in the figure and refer to the axis values on the right-hand side of each panel. The charge state affected other residue interactions beyond the charged residue, as can be seen by comparison between neutral and charged versions of helices B, C, and D in the figure.

Combining the individual residue interactions gives insight into the microscopic interactions underlying the GES and similar hydropathy scales. The results are shown in Figs. 4-7 and Table 2. The figures show that the decreased luctuations at the center of the bilayer are not specific to particular amino acids, but instead are present in all of the amino acid types. That is, the energetic fluctuations tend to be greater near the ends for each type of amino acid. A further comparison is made in Table 2. The columns represent the relative hydrophobicity assigned to an amino acid by the GES scale, and the estimated average interaction energy for the amino acids at the interface and in the interior of the bilayer. The table shows that for the hydrophobic amino acids, there is a consistency in the ordering of amino acid type with the GES scale.

The aromatic residues do not appear similar to what might be expected from the GES scale. These side chains have a strong interaction energy with the lipid headgroups and water at the interface, as seen in Fig. 5. With the CHARMm empirical potential function, the larger interaction at the interface for Trp and Tyr residues is partly due to the possibility of hydrogen bonds from ring nitrogen/oxygen to water or lipid. The GES scale also predicts a relatively large difference between the Asn, Gln, Ser, Thr residues that is not observed in the simulations. This could be due to a lack of hydrogen-bonding partners for these side chains in the simulations.

Helical stripes

The hydrophobic moment analysis of Eisenberg and co-workers (Eisenberg, 1984; Eisenberg et al., 1984) attempted to distinguish those faces of a membrane -helix that are most likely to interact strongly with lipid from those more likely involved with helix-helix interactions. This concept was tested in the explicit helix-lipid simulations by selecting atoms in 45° wedges of helix interacting with the lipid. Fig. 8 shows the approach. A clear preference for certain wedges of the helix over others would be indicative of the veracity of the concept in this particular case.

View larger version (40K): [in this window] [in a new window]   FIGURE 8   Illustration of the wedge-shaped selection of atoms for evaluating interactions between a cylindrical stripe and the environment.

The results show evidence for helical stripes. Fig. 9 shows the vdW interaction energy between the DMPC and all 10 helices. The presence of minima at different locations for all 10 helices is clearly evident in the bottom curves of each panel. This interaction is greatest for favorable vdW contacts between lipid alkane chain and protein. Equally strongly shown is the presence of a face that is disfavored in the interaction energy. The contour plots show a clear periodicity in the turn length along the y axis. This periodicity is clearly related to the -helix conformation. The darker contour lines are the stronger contributions to the total, and the left-hand axis is the residue number for each helix. It is interesting to note that some helices have a strong interaction stripe extending over the full length of the helix (e.g., helix A), whereas others have a stripe that is dominated by a subset of the full helix (e.g., helix C-neutral).

View larger version (79K): [in this window] [in a new window]   FIGURE 9   Helical stripes of vdW interaction energy determined by selecting atoms within a pie-shaped vertical stripe of each -helix. The interaction energy between protein and DMPC was determined for all atoms within a region of 45°. The regions were rotated by 5° increments around the full helix. The contour plots were generated by considering the contribution of individual residues to the stripe's interaction.

Fig. 10 illustrates the situation for electrostatic interactions of each helix with the lipid environment. Minima are less obviously present, suggesting that the set of bacteriorhodopsin helices is not adopted by evolution to a preferred orientation by electrostatic interactions. A strongly favorable stripe would indicate a set of strong interactions with the glycerol and headgroup regions of the lipid with the protein. Helices A and G have the strongest clear minima. Helix C-charged and helix F have smaller, clearly defined minima. In all cases, the stripe with the strongest electrostatic interaction does not correspond to the strongest vdW interaction.

View larger version (102K): [in this window] [in a new window]   FIGURE 10   Similar to Fig. 9, for electostatic interactions.

Tables 3 and 4 compare the simulation results to the recent (2BRD) bacteriorhodopsin structure. Table 3 compares the lipid-exposed residues of each helix (from visual inspection) to those in the stripes with the strongest calculated lipid interactions. As expected, there is correspondence between the lists. However, because hydrophobic moment analysis neglects helix-helix interactions, the mapping between the lists is not one to one. That is, some of the predicted residues within stripes are involved with helix-helix interactions, and not all possible stripes are seen in the final structure.

Table 4 continues this analysis by presenting data from the 2BRD structure for the 10 lipids that were resolved. It is interesting to note that some lipids had a large amount of contact with helices (e.g., 266, 267), whereas others had only a small number of contacts (e.g., 262, 264). A comparison of the optimal wedges from the simulations with the contacts formed by the lipids in the structure suggests further similarities. For example, helix A is contacted by five lipids (263, 264, 265, 269, and 270). Included in this set of near-contacts is the simulation optimal residues 17 (lipid 270), 3 (lipid 263), 10 (lipid 265), and 6 (lipid 263).

	DISCUSSION

Top Abstract Introduction Methods Results Discussion Conclusion References

The scales for predicting -helical membrane protein secondary structure are based on the assignment of a relative hydrophobicity for amino acids (e.g. Engelman et al., 1986). Ideally, the transfer free energy for a particular amino acid within defined locations of an -helix could be computed from a fully detailed microscopic model. Currently such a calculation is not feasible. Recently the Honig group began the first steps toward detailed calculations of the transfer energy by considering the effects of the -helix orientation on transfer free energies in a continuum electrostatic model for the bilayer environment (Ben-Tal et al., 1996). The current results complement both those calculations and hydrophobicity scales by calculating the molecular details of interaction for individual amino acids at different locations along a transmembrane -helix within the lipid bilayer environment. The interaction energies, although not the full transfer energy, show trends to be expected from the GES scale for hydrophobic residues.

The strongest exceptions to the expectation of the GES scale are the aromatic residues. This is due to the favorable interactions that can be formed by the aromatics at the interface (Fig. 5). This favorable location for Trp and Tyr residues is due to the possibility of hydrogen bonds from ring nitrogen/oxygen as well as hydrophobic interactions with the lipid. This observation has been made before on the basis of structural analysis (Schiffer et al., 1992; Landolt-Marticorena et al., 1993) and from simulations of gramicidin in DMPC (Woolf and Roux, 1994, 1996). Recently, an interfacial hydrophobicity scale for proteins was designed based on these ideas and further experimental evidence (Wimley and White, 1996). This all provides support for an extension to the GES and related scales for the full membrane bilayer that responds to the amino acid location along the helix.

Another exception to the expectations of the GES scale was observed with the polar amino acids. The simulation results suggest a strong similarity among Ser, Pro, Thr, and Met in the center of the bilayer. This contrasts with an assumption used in the GES scale that Ser/Thr/Cys will satisfy internal hydrogen bonds with the carbonyl oxygens of the peptide backbone (Gray and Matthews, 1984; Engelman et al., 1986). From analysis of the side-chain dynamics as well as the interaction energy analysis, there does not appear to be a strong hydrogen-bonding component to the interaction energy of these particular side chains.

Hydrophobic moments have been used for the prediction of membrane protein tertiary structure (e.g., Taylor et al., 1994; Suwa et al., 1995). The current results (Figs. 9 and 10) suggest that the interaction differences between angular stripes of -helix with the environment can be very large. This reinforces the utility of such an analysis for tertiary structure prediction. The comparison in Figs. 9 and 10 between the vdW and the electrostatic components of the stripe interactions suggests that hydrophobic moments could be improved by consideration of the interface region. That is, as in an amended hydrophobicity scale that includes residue location in prediction of membrane spanning regions, it may be possible to improve hydrophobic moment predictions by considering the location of the amino acid within the helix in the calculation of the preferred lipid-exposed face.

Two calculations of hydrophobic moments, with consideration of the solvent, have been reported in the literature. The first of these used a continuum model for the solvent and calculated the preferred face of the -helix in contact with the solvent for structure prediction (Suwa et al., 1995). This was used in the context of evaluating the rotational positions of helices in bacteriorhodopsin. That is, the approach was similar to the Taylor paper (Taylor et al., 1994), but some account was taken of the solvent. The second paper used explicit representations of the solvent as propane (nonpolar) or water for calculating hydrophilic and hydrophobic faces (Efremov and Vergoten, 1995). Four different -helices were used, including helix A from bacteriorhodopsin. Interestingly, the preferred hydrophobic face from the Monte Carlo calculations consisted of residues 2, 6, 10, 13, 17, and 18, which are similar to those predicted from the current calculations (3, 6, 10, 17). Neither the paper by Taylor et al. (1994) nor that by Suwa et al. (1995) provides a listing of the predicted residues in the lipid-exposed faces with their implementations of hydrophobic moments. It appears, from the figures in the papers, that for helix A of bacteriorhodopsin, a similar orientation was calculated in the two cases.

The prediction of membrane protein tertiary structure is a difficult and important problem. Several papers have already been mentioned that use hydrophobic moments as part of the analysis (Taylor et al., 1994; Suwa et al., 1995). Other prediction methods have been used (e.g., Herzyk and Hubbard, 1995; Sansom et al., 1995; Tuffery et al., 1994). In all of these cases bacteriorhodopsin was an important test case. The approach of Tuffery et al. (1994) concentrates on helix-helix interactions as the predominant effect. This is similar to an approach used for glycophorin and phospholamban (Treutlein et al., 1992; Adams et al., 1995). Herzyk and Hubbard (1995) using all available experimental information to derive restraints within a reduced representation. In their approach, hydrophobic moments did not seem to aid the approach to a minimum. Sansom et al. (1995) used templates and simulated annealing methods to explore the properties of seven-helix bundles. Another approach (e.g., Cronet et al., 1993; Zhang and Weinstein, 1993) is largely based on bacteriorhodopsin as a template for G-protein-coupled receptors with interactive graphical modeling of changes. Baldwin (1993) used multiple sequence alignment to suggest the lipid-exposed faces as those regions with the greatest variability. Reviews of the forces that contribute most to membrane protein structure and function have frequently emphasized helix-helix interactions as the predominant organizer (e.g., Bormann and Engelman, 1992; Cramer et al., 1992; Lemmon and Engelman, 1994). The results of the current calculations suggest that learning more about the properties of the membrane bilayer interacting with the protein will help to improve our understanding of membrane protein structure and function. Such work will further act to improve the computer prediction of membrane protein tertiary structure. The importance of understanding the lipid bilayer as the medium for membrane protein folding is emphasized by the result that misfolded globular proteins could only be discriminated in the presence of solvent (Novotny et al., 1984).

The concept of a population of "boundary" lipids surrounding a membrane protein was initiated in the early 1970s and advanced through a series of experiments up to the 1980s (e.g., Jost et al., 1973; Kang et al., 1979; Lavialle et al., 1980; Thomas et al., 1982; Post and Dijkema, 1983; Nishiya et al., 1987; Van Gorkom et al., 1990). The current results suggest that a small number of lipids surrounding a central -helix tend to interact considerably more strongly than other nearby lipids or than lipids in a second, "outer" zone. This is seen in Fig. 1, where there is a small but significant population of lipids interacting with energies from 15 to 100 kcal/mol. Nonetheless, it is important to caution that the calculations used image lipids that may not correctly reflect the conformations that a true set of lipids in a second zone around the central helix would adopt. Furthermore, the preferred definition of a "boundary" lipid implies a slow rate of exchange between the lipid closely associated with membrane protein and the bulk lipid. The current molecular dynamics time scale is not sufficient to explore the exchange process. Even with the above cautions, it is intriguing that the different -helices show different degrees of strongly interacting lipids. The energy breakdown in Fig. 2 emphasizes that these strongly interacting lipids are a combination of vdW interactions (mediated by the chains) and electrostatic interactions (mediated by the headgroup).

Further simulations are ongoing to expand the current molecular dynamics "data set" of helices in lipid settings. Important questions regarding the importance of bilayer thickness and the effective ability of a micelle to mimic a bilayer can be addressed with these additional simulations. By simulating the same -helices within a water environment, an estimate will be made for the full thermodynamics of transfer. A further comparison is in progress between the current NVE ensemble calculations and the results from the constant surface tension/pressure ensemble (Feller and Pastor, 1996).

	CONCLUSION

Top Abstract Introduction Methods Results Discussion Conclusion References

Hydrophobicity scales, as currently applied to membrane proteins, assume that a single transfer value can be assigned to each amino acid without regard to where it might be found within a structure spanning the bilayer. This effectively assumes that the -helix is either in water or is a uniform "bilayer-like" solvent for all amino acid residues. In reality, as multiple investigators have suggested (e.g., White, 1994; Woolf and Roux, 1996), the bilayer environment is not a simple uniform medium, but a complex mix of effective environments. In particular, the interfacial region is quite different from the alkane-dominated center of the bilayer. Thus a better picture of the partitioning of amino acids might be achieved with a model for the bilayer that explicitly considers the existence of an interfacial region. Membrane-spanning regions could then be predicted with greater confidence by scoring for location, as well as hydrophobicity in the bilayer.