INTRODUCTION

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The binding of peptides to lipid bilayers is a subject of great interest, from both a biomedical and biophysical standpoint. Small, surface-active peptides that bind and affect the fluid order of membranes can lyse and destroy invasive microorganisms (Epand and Vogel, 1999) and inhibit (Damodaran and Merz, 1995) or promote (Colotto et al., 1996; Epand et al., 1999) membrane fusion. More generally, detailed investigations of peptide-bilayer binding can provide insight into vital biological phenomena such as protein-membrane interactions, channel formation, and nonspecific transport of molecules across cellular interfaces. Although experimental techniques such as NMR (Seelig, 1977; Opella, 1997; Klassen and Opella, 1997), x-ray diffraction (Wallace and Janes, 1991), neutron diffraction (King and White, 1986; Jacobs and White, 1989), and calorimetry (Seelig, 1997) have provided a wealth of data on peptide-bilayer interactions, there remains a need for detailed, molecular-level information on the structure and energetics of peptide-bilayer complexes. Moreover, the mechanisms and kinetics of the binding process should be analyzed. The most direct way of obtaining such microscopic information is through computer simulation techniques. In particular, molecular dynamics (MD) has become a popular method for studying peptide-bilayer systems, with many MD studies of bilayer-bound peptides to be found in the literature (Bernèche et al., 1998; Damodaran et al., 1995; Damodaran and Merz, 1995; Kothekar, 1996; Kothekar et al., 1996; Huang and Loew, 1995; Tieleman et al., 1998; Forrest et al., 1998; Tieleman et al., 1999b).

To perform an MD simulation on a peptide-bilayer system, one commonly uses available experimental data and computational manipulation (Bernèche et al., 1998; Woolf and Roux, 1994) to generate a reasonable structure for the peptide already bound to the bilayer interface or inserted into the bilayer interior. Following a suitable equilibration procedure, a production simulation can be carried out from which one obtains information on the structure, energetics, and dynamics of the system and its components. In situations where only low-resolution structural data are available, however, atomic details of the prepared system must necessarily be specified by the investigator. Due to present-day limitations on simulation length, the system is frequently unable to evolve significantly from the starting configuration, and therefore any bias in the starting conditions will affect the results. As noted by Damodaran and Merz (1995) many simulations with different starting configurations should be carried out to avoid this bias and make more meaningful comparisons to experiment. Another shortcoming of standard MD simulations on these systems is the lack of information provided on the binding process itself. Ideally, one would like to carry out simulations in which the peptide is gradually placed into contact with the bilayer and extract thermodynamic information on the binding process using a method like thermodynamic integration (van Gunsteren and Berendsen, 1987a). Unfortunately, the computational expense of such simulations can become prohibitive when full atomic detail is used for the surrounding solvent (Juffer et al., 1996). Also, the effect of ionic strength is commonly ignored.

In this paper we present the results of two separate simulations of a peptide-bilayer system. In the first of these simulations, we have initially placed the peptide outside the interface of the bilayer in a region of bulk water density. Within 2 (ns, the peptide becomes adsorbed (bound) to the surface of the bilayer. The details of the peptide-bilayer complex thus formed are solely a result of the interactions between the various components of the system. As such, this simulation partially circumvents the shortcomings of more conventional simulations discussed above, in that 1) the bound complex is relatively free of bias and 2) information on the binding process is provided. Using a bound configuration of the peptide at the bilayer surface from the first simulation, the N-terminus of the peptide was deprotonated and a second simulation was carried out. A deeper penetration of the deprotonated peptide into the bilayer was observed. Again, this occurred solely as a result of the interactions between the peptide and bilayer. Corresponding to the protonation state of the N-terminus, we will hereafter refer to these simulations as the protonated and deprotonated simulations and the peptide-bilayer complexes observed therein as the protonated and deprotonated complexes, respectively. Two additional simulations of the peptide in water and the bilayer in water enable us to compare features of the complexes with that of the free peptide and bilayer, respectively.

The peptide we have chosen for our study is sandostatin (also known as octreotide), a cyclic analog of the peptide hormone somatostatin (Maurer et al., 1982). The bilayer is composed entirely of 1-palmityl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipids and is therefore electrostatically neutral. Fig. 1 shows schematic diagrams of sandostatin and a POPC lipid. The binding of sandostatin to lipid structures composed of neutral lipids (POPC), negatively charged lipids (1-palmityl-2-oleoyl-sn-glycero-3-phosphoglycerol (POPG)), and mixtures thereof has been studied extensively by Beschiaschvili and Seelig (1990, 1991, 1992) using a variety of experimental techniques. Indeed, the choice of this particular system for the current study was motivated by the breadth and quality of their experimental data, which includes information on the secondary structure, orientation, and protonation state of the bound peptide as well as measurements of lipid headgroup orientation. Another beneficial feature of the system is the disulfide bridge of sandostatin, which locks the peptide into a rigid hairpin structure. The relatively small conformational change occurring in the peptide upon binding simplified the experimental thermodynamic analyses (Beschiaschvili and Seelig, 1992) and in a like manner helps to ensure that the conformational space of the peptide is adequately sampled during the simulations. This would not be the case for a larger, linear peptide in which helix-coil or other conformational changes may occur as a result of binding (Epand and Vogel, 1999; Hwang and Vogel, 1998). As we will show in the Discussion, we have found very good agreement between properties of the peptide-bilayer complexes in the simulations and those obtained from experiment. Moreover, features of the binding process have been analyzed from the first 2 ns of the protonated simulation.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   Schematics of sandostatin and a POPC lipid. Segments of the lipid are named according to Jacobs and White (1989), with the modification that all methylene segments on both chains are included in the CH(2) group.

	SIMULATION DETAILS

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS DISCUSSION CONCLUSIONS REFERENCES

All simulations were carried out using the GROMACS simulation package (Berendsen et al., 1995), both in serial mode on a single DEC  500au workstation and in parallel on a cluster of four identical processors. All analyses were performed using software provided in the GROMACS package or with code written specifically for the project. To investigate the binding of sandostatin to the POPC bilayer we have carried out simulations of the following systems: 1) sandostatin in water, 2) fully charged sandostatin in the presence of the bilayer (with the peptide initially placed outside the bilayer interfacial region), 3) bilayer-bound sandostatin with deprotonated N-terminus, and 4) the POPC bilayer in water. Details of all four simulations are provided below.

Sandostatin/water (SW) simulation

The starting structure for sandostatin was taken from Brookhaven Protein Databank (PDB) (Bernstein et al., 1977) entry 1SOC (Melacini et al., 1997). Interaction parameters within the peptide were taken from the Gromos-87 force field (van Gunsteren and Berendsen, 1987b) using the modifications of van Buuren et al. (1993). Explicit hydrogens were included on aromatic residues (van der Spoel et al., 1996). Because aliphatic carbons are treated as extended atoms, the equilibrium values for the improper dihedrals on the C atoms of D-Trp1 and D-Phe3 were explicitly changed to the negative of the value for L-amino acids in the force field. This maintained the chirality of these residues throughout the simulation. The reduced C-terminus was built onto the peptide using standard parameters for methylene, hydroxyl oxygen, and hydroxyl hydrogen atoms from the force field. A molecular topology file for the sandostatin peptide complete with these modifications will be made available at the online GROMACS topology repository (http://rugmd4.chem.rug.nl/TILDEgmx/topologies.html). The peptide was placed in a rectangular box 2.6 × 2.2 × 3.0 nm and solvated using multiple equilibrated copies of 216 SPC (Berendsen et al., 1981) water molecules. Any water molecule whose oxygen atom approached a peptide atom to a distance smaller than the sum of their respective van der Waals radii was then removed, resulting in a total of 449 waters. To neutralize the positive charges on the N-terminus and lysine side chain of sandostatin the GROMACS program genion (van der Spoel et al., 1999) was used to replace two of the remaining waters with chloride counterions. The energy of the final system, which consisted of 1488 atoms, was minimized using a steepest-descents algorithm and a convergence criterion of 0.001 kJ mol¹. Following energy minimization, the system was subjected to 20 ps of position-restrained dynamics, in which an isotropic force constant of 1000 kJ mol¹ nm¹ was applied to all peptide atoms. This allowed equilibration of the solvent. During restrained dynamics the temperature was controlled through weak coupling to a bath of constant temperature (Berendsen et al., 1984) using a coupling time _T and reference temperature T_o of 0.1 ps and 300 K, respectively. The peptide, water, and counterions were all coupled separately to the temperature bath. Analogously, the system was maintained at constant pressure through weak coupling to a pressure bath using P_o = 1.0 bar and _p = 0.5 ps. The center of mass motion of the system was removed every step to maintain the effective simulation temperature at 300 K. The LINCS algorithm (Hess et al., 1997) was used to restrict all covalent bonds to their equilibrium lengths, and the water molecules were kept rigid with the SETTLE algorithm (Miyamoto and Kollman, 1992), allowing for a time step of 2 fs. A twin-range, group-based cutoff was employed for the evaluation of all nonbonded interactions, with cutoffs of 1.0 and 1.8 nm for short-range (Lennard-Jones and Coulomb) and long-range (Coulomb) forces, respectively. Long-range forces were updated during generation of the neighbor list, which was done every 20 fs. The final coordinates from the restrained run served as the initial configuration for a 2.5-ns production run, the first 500 ps of which were used as equilibration and were excluded from analyses.

Protonated simulation

The starting structure for sandostatin in the protonated simulation was taken from the final frame of the SW simulation. The starting structure for the lipid bilayer consisted of 128 POPC lipids (64 per leaflet) and was taken from a 1.6-ns simulation of the bilayer solvated with 3655 water molecules (D. P. Tieleman, personal communication). This equilibrated bilayer has been used as a starting point for the simulation of transbilayer helices and proteins by Tieleman and colleagues (Tieleman et al., 1998; Forrest et al., 1998). As in their studies, we have used the lipid force-field parameters of Berger et al. (1997) to describe the lipid-lipid, lipid-water, and lipid-protein interactions and retained the Gromos-87 force field (van Gunsteren and Berendsen, 1987b) for the intramolecular protein interactions. For our simulation, we retained approximately the same x and y dimensions of the bilayer system mentioned above, but increased the size of the simulation box on either side of the bilayer (z dimension) to allow space for sandostatin and solvent. Thus, the bilayer was centered in a rectangular box of dimensions 6.34 × 6.38 × 9.92 nm, with the bilayer normal oriented parallel to the z axis. This provided a region approximately 2.8 nm thick on either side of the bilayer. The peptide was then placed above one leaflet (hereafter referred to as the upper leaflet) so that it was centered with respect to the x and y dimensions. Because periodic boundary conditions were used in all simulations, we oriented the peptide with its long axis roughly parallel to the membrane surface. Aside from this specification, which helped to keep the system size reasonable and to ensure that the peptide would not interact with the images of the bilayer above the central simulation box, the orientation of the peptide was completely arbitrary. That is, no attempt was made to place the peptide in an energetically favorable position or orientation with respect to the bilayer. An estimate of the peptide's translational diffusion coefficient from the SW simulation (see above) showed that the mean distance traveled by the peptide in 2 ns was ~1.5 nm. Thus, the location of the peptide along the system's z axis was chosen such that the center of mass of the peptide was 1.5 nm away from the average z coordinate of the phosphatidylcholine nitrogen atoms in the upper leaflet. At this location, the minimum distance between the bilayer and peptide was approximately 1 nm. Next, the peptide and bilayer were solvated with water molecules and two chloride counterions in the manner described for the SW simulation. However, due to the amount of free volume in the hydrocarbon region of the bilayer, this procedure allowed some water molecules to be artificially placed in the bilayer interior. Based on the neutron diffraction data of Jacobs and White (1989), we additionally removed any waters that penetrated the bilayer deeper than any phosphate atom. This resulted in a system with 6731 water molecules and a total of 26,954 atoms. Aside from the cutoff for long-range nonbonded interactions, which was reduced to 1.7 nm, all energy minimization parameters were identical to those employed in the SW simulation. This value of the cutoff, which was deemed suitable for bilayer simulations in a systematic study of simulation parameters by Tieleman and Berendsen (1996), was used in both of the peptide-bilayer simulations. However, it should be noted that the use of short cutoffs can result in inaccurate quantitative results for interfacial systems (Feller et al., 1996). Other methods for treating electrostatic interactions, such as the three-dimensional Ewald summation and its variants are becoming steadily more popular. Currently, new methods being developed to cope with the two-dimensional periodicity of interfacial systems (Wong and Pettitt, 2000; Yeh and Berkowitz, 1999) have not yet achieved the level of standards. We feel that the use of a 1.7-nm cutoff, which is quite long by current standards, is justified given the following: 1) following a diffusive motion of the peptide toward the bilayer in the first 0.5 ns of the protonated simulation, the atoms of the peptide all remain well within the cutoff distance of the charged bilayer atoms in both bilayer simulations, and 2) the results presented are primarily qualitative in nature, being analyses of relative changes in various properties of the system as a result of the binding process, rather than an attempt to accurately calculate quantitative properties of the system. Valuable results obtained from similar systems using electrostatic cutoffs are prevalent in the recent literature (Lindahl and Edholm, 2000; Shrivastava and Sansom, 2000; Tieleman et al., 2000; Capener et al., 2000; Shrivastava et al., 2000; Forrest et al., 2000).

Following energy minimization, a force constant of 1000 kJ mol ¹ nm¹ was applied to the peptide atoms for 500 ps of position-restrained dynamics, allowing equilibration of the bilayer and surrounding solvent. The final frame of the restrained run served as the starting configuration for a standard MD simulation. Temperature coupling, restraint algorithms, neighbor list update, and time step were all implemented as in the SW simulation, with the peptide, lipids, water, and counterions being coupled separately to the temperature bath. In contrast, anisotropic pressure coupling was employed by scaling all directions separately to 1 atm with a time constant of 0.5 ps, corresponding to zero surface tension. The total volume and energy of the system were monitored as criteria for equilibration. Both of these properties relaxed to constant values (data not shown) in the first 250 ps of simulation, which were discarded. The remaining 6 ns of the simulation were used for analysis.

Deprotonated simulation

All simulation parameters for the deprotonated simulation were identical to those used in the protonated run. The starting structure was taken from the protonated trajectory at the 4-ns mark, at which time the peptide was stably adsorbed to the bilayer surface (see Results). The N-terminus of the peptide was deprotonated, and new atomic charges from the force field were assigned to the relevant atoms. One of the counterions was replaced by a water molecule to achieve electrostatic neutrality. Energy minimization and 500 ps of restrained dynamics were carried out as discussed above. Following an equilibration period of 250 ps, 6.0 ns of dynamics were performed on this system. Thus, the combined simulation time for the protonated and deprotonated trajectories is 12 ns.

Bilayer/water (BW) simulation

All simulation parameters in the bilayer water simulation were identical to those used in the protonated simulation. The starting structure for the system was generated by cutting the peptide out of the initial configuration of the protonated simulation and resolvating the system in the manner discussed for the peptide-bilayer simulations. The two chloride atoms were replaced by waters, resulting in a total of 6765 water molecules. Following energy minimization, the system was subjected to 500 ps of standard MD as equilibration. Two additional nanoseconds were used for analysis.

	RESULTS

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Peptide-bilayer distance

Snapshots from the protonated and deprotonated simulations in Fig. 2 show the gradual approach of the peptide to the bilayer surface and a representative structure of the deprotonated complex, respectively. The same information is displayed quantitatively in Fig. 3, which shows the distance between the center of mass of the peptide and that of the upper leaflet phosphate atoms for both simulations. Despite some apparent movements of the peptide away from the bilayer in the first nanosecond of the protonated run, the peptide-bilayer distance decreases in an essentially linear fashion with respect to time up until 2 ns. At this time, the distance between sandostatin and the bilayer levels off. From Fig. 2, it is apparent that this is due to adsorption of the peptide to the bilayer surface. The peptide-bilayer distance then remains essentially constant for the next 3 ns and then decreases further between 5 and 6 ns. The relatively small fluctuations in the distance between 2 and 6 ns would seem to indicate that the adsorbed complex is stable on the time scale of the simulation. However, the slight decrease in the distance during the last nanosecond of the simulation indicates that the peptide would like to immerse itself further into the bilayer phase.

View larger version (63K): [in this window] [in a new window]   FIGURE 2   Snapshots of the binding process from the protonated trajectory (numbered) and a representative structure of the deprotonated complex. For clarity, only the upper leaflet of the bilayer and sandostatin are shown. The nitrogen atoms of the phosphocholine headgroups are shown in the space-filled mode. Pictures were made using the program RASMOL (Sayle and Milner-White, 1995).

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Distance between the center of mass of the peptide and that of the upper leaflet phosphate atoms as a function of simulation time. Solid line, protonated run; dashed line, deprotonated run.

The increase in bilayer penetration by the peptide following deprotonation is clearly seen in the peptide-bilayer distance for the deprotonated run. Interestingly, this decrease (approximately 0.5 nm relative to the protonated complex) was achieved mostly during the energy minimization step of the deprotonated simulation (see Simulation Details). Therefore, although sufficient space for penetration of sandostatin into the bilayer existed during the protonated simulation, this was probably prevented by electrostatic repulsions between positive charges on the peptide and lipid headgroups. The peptide-bilayer distance fluctuates more in the first half of the deprotonated simulation than in the protonated complex; logically, the deeper penetration of sandostatin in the deprotonated complex perturbs the bilayer more than the surface-bound peptide and more time is required for the membrane to adjust to the inclusion. The fluctuations settle down between 4 and 6 ns and become comparable to those of the protonated complex. Thus, on the time scale of the simulations, both complexes enjoy a similar stability.

Peptide orientation

To measure the overall orientation of the peptide during the simulation, we have employed a method similar to the one used by Tieleman et al. (1989) for measuring the orientation of aromatic side chains with respect to a bilayer surface. Given a planar molecule, the method is as follows: two parameters, S_N and S_L, describe the alignment (with respect to the bilayer normal) of a vector normal to the molecular plane and a long axis vector lying in the molecular plane, respectively. Both parameters can assume the values 1 (for a vector lying parallel to the bilayer normal) to 0.5 (vector lying perpendicular to the bilayer normal) but are not mutually independent. Fig. 2 of Tieleman et al. (1998) provides a more detailed explanation, using a tyrosine ring as an example. Fig. 4 shows S_N and S_L parameters calculated for sandostatin as a whole, using specific peptide atoms to define a plane and long axis for the molecule (see figure caption for details). Starting from similar values at the start of the protonated simulation, both S_N and S_L undergo large fluctuations within the first nanosecond. In fact, the large changes in the value of S_N correspond closely with the increases in the peptide-bilayer distance in this same time period (see Fig. 3). These fluctuations in distance can therefore be attributed to overall rotation of sandostatin as a whole, rather than translational motion of the peptide away from the bilayer surface. Between 1 and 2 ns, the peptide becomes oriented with its long axis parallel to the bilayer normal. Because the long axis vector points away from the N- and C-termini of the peptide, the negative value of S_N reveals that it is the terminal ends of the peptide that face the bilayer (see also Fig. 2), with the turn around residues Trp4 and Lys5 facing in the opposite direction. Upon adsorption at 2 ns, a rapid orientational transition is observed: the normal to the plane S_N becomes aligned with the bilayer normal. Thus, the long face of the peptide lies against the bilayer surface. The value of S_L also reflects this change, although it is clear that the long axis of the peptide retains some freedom of orientation with respect to the bilayer normal. In the last nanosecond of the simulation, however, the peptide comes to lie completely flat against the bilayer surface. This change is concomitant with the decrease in the peptide-bilayer distance observed during this time period (see Fig. 3), indicating that binding of the peptide to the bilayer favors the flat orientation. Correspondingly, this type of orientation is also observed in the deeply inserted peptide of the deprotonated complex: the S_N and S_L values indicate a rigid, parallel orientation of the peptide with respect to the bilayer surface and do not undergo large fluctuations. Presumably, the environment surrounding sandostatin is much more rigid than in the protonated simulation and exerts a greater hindrance on peptide rotation.

View larger version (29K): [in this window] [in a new window]   FIGURE 4   Orientation of sandostatin in the protonated (upper) and deprotonated (lower) runs. The vector that determines the value of SN is a normal to the plane defined by the C atoms of Cys2, Lys5, and Cys7. A vector pointing from the S atom of Cys2 to the C atom of Lys5 was used to calculate SL. Values of the parameters are interpreted in the same manner as for aromatic rings (see Tieleman et al., 1998). Solid line, SN; dashed line, SL. Running averages over 50 ps are shown for clarity.

Peptide conformation

Fig. 5 shows the root-mean-square fluctuation (RMSF) values in the peptide atom positions for the SW simulation, the first 2 ns of the protonated run, the last 4 ns of the protonated run, and the entire deprotonated run, respectively. The profiles from all four data sets are similar in that the terminal residues and those located near the top of the turn are relatively mobile, with flexibility decreasing toward the cysteine residues engaged in the disulfide bond. However, the magnitudes of the values in the various profiles differ greatly, depending on whether the peptide is free or bound to the bilayer. The flexibility of all residues is roughly the same in the SW simulation and the first 2 ns of the protonated simulation, with the peptide being only slightly more rigid in the latter. With the exception of the N-terminal phenylalanine, all residues of the peptide become much more restricted in the latter 4 ns of the protonated simulation (protonated complex), and the most significant reductions in flexibility are seen in those residues that penetrate deepest into the bilayer (see below). The data from the deprotonated simulation shows that the apparent translational and rotational rigidity of the deprotonated complex also extends to the intermolecular degrees of freedom: the RMSF values for all residues are reduced relative to the protonated complex. Moreover, the RMSF values for all residues in this complex are roughly similar: even those near the top of the turn have experienced a significant loss of conformational freedom.

View larger version (19K): [in this window] [in a new window]   FIGURE 5   RMSF of sandostatin during the simulations. Solid line, entire SW simulation; dotted line, first 2 ns of the deprotonated run; dashed line, last 4 ns of the deprotonated run; long dashed line, the entire deprotonated run. The RMSF for each residue was averaged over all atoms in the residue.

For both complexes, additional evidence of structural rigidity was seen in the secondary structure and ,  angles of the peptide. While bound, very small fluctuations were seen in the average backbone dihedrals, with each residue sampling small and nearly separate regions of ,  space (data not shown). Consistent with the RMSF plots, the spread in the dihedral angle values was smaller in the deprotonated than in the protonated complex. The secondary structure of sandostatin in both complexes was essentially identical as determined by the DSSP program (Kabsch and Sander, 1983). A single, stable hydrogen bond was formed between the carbonyl oxygen of Cys2 and the nitrogen of Thr8, forming a -bridge structure separated by a bend. No other intramolecular hydrogen bonds were possible, as the main-chain atoms of the remaining residues engaged in extensive hydrogen bonding with the lipids (data not shown). The existence of a single, persistent conformation differs from observations made in the SW simulation and the first 2 ns of the protonated run: during these times, the peptide was observed to sample both of the major conformational families determined from solution NMR studies of sandostatin (Melacini et al., 1997; data not shown).

Buried surface area

Fig. 6 shows the amount of surface area buried at the peptide-bilayer interface for both the protonated and deprotonated trajectories. This was calculated as follows: the peptide and bilayer were taken separately from each recorded configuration of the trajectories and the sum of their surface areas was calculated. The area was then calculated for the peptide and bilayer together and was subtracted from the sum of the two species taken independently, giving the total surface area buried in the complex. The ASC program of Eisenhaber and Argos (1993) was used for all area calculations, using atomic radii from GROMACS. From the data for the protonated run, a linear increase in the buried surface area is observed between 1 and 2 ns, followed by a plateau. The absence of contact between the peptide and bilayer during the first nanosecond of the simulation explains the unrestricted rotation of the peptide during this time period (see Fig. 4). Aside from a slight increase between 5 and 6 ns, the value remains fairly constant for the remainder of the simulation. The interfacial area of the protonated complex is 40% hydrophobic and 60% hydrophilic in character, respectively. The buried surface area increases greatly following deprotonation of the N-terminus, with the value jumping to almost twice that of the protonated complex. Aside from transient changes near the middle and ends of the run, the size of the interface remains essentially constant for the entire simulation. The character of the interface has changed slightly from that observed in the protonated simulation, with hydrophobic and hydrophilic atoms contributing equally to the total.

View larger version (27K): [in this window] [in a new window]   FIGURE 6   Surface area shared by sandostatin and the bilayer as a function of simulation time. Contributions from hydrophilic and hydrophobic surface areas are shown along with the total. An atom was considered hydrophobic if the magnitude of its partial charge was less than 0.2. Running averages over 50 ps are shown for clarity.

Intermolecular forces

Intermolecular forces (Lennard-Jones and Coulomb) on the peptide due to the bilayer, solvent (water and counterions), and their sum is shown in Fig. 7. Only the z components are shown as they act normal to the bilayer surface and thus control the binding process. For the first 750 ps of the protonated simulation, the total force component and its contributions essentially fluctuate around zero, and the peptide is roughly at equilibrium with respect to the normal forces. After this time, however, the bilayer exerts a significantly attractive force on the sandostatin molecule, reaching a minimum when the peptide becomes completely adsorbed at 2 ns. Aside from a small fluctuation at ~2.6 ns, the bilayer force is negative (favors adsorption) and fairly constant between 2 and 5 ns. In the last nanosecond of the simulation, the cause of the decrease in the peptide-bilayer distance (see Fig. 3) becomes clear; the force on the peptide due to the bilayer becomes gradually more negative and draws sandostatin deeper into the lipids.

View larger version (30K): [in this window] [in a new window]   FIGURE 7   The z coordinates of the intermolecular forces acting on the peptide as a function of simulation time. (Upper plot) Protonated trajectory; (Lower plot) Deprotonated trajectory. Thin solid line, force due to the bilayer; dashed line, force due to the solvent (water and counterions); thick solid line, total intermolecular force on the peptide. Running averages over 50 ps are shown for clarity.

The most striking feature of the solvent force component in the protonated simulation is its near mirror-image quality with respect to that of the bilayer between 750 and 2500 ps. For the bulk of this time period the bilayer contribution dominates the total force but is consistently opposed by near-simultaneous solvent repulsions, presumably caused by expulsion of water molecules from the space between the interacting species. Between 2.5 and 5 ns, the solvent force becomes comparatively independent of the bilayer force and oscillates around zero. Concomitant with the increased peptide-bilayer interaction between 5 and 6 ns, this solvent force then disappears entirely.

The z components of the intermolecular force for the deprotonated simulation are shown in the lower plot of Fig. 7. Like the last nanosecond of the protonated simulation, the contribution of the solvent to the total force is small and essentially disappears after the first nanosecond. The value of the bilayer contribution is negative throughout, indicating a net force attempting to draw sandostatin even deeper into the bilayer.

Density plots

The four plots in Fig. 8 show the density distributions of various groups along the bilayer normal calculated as averages over the first and last 500 ps of the protonated (upper plots) and deprotonated (lower plots) runs, respectively. Before discussing the location of the peptide in the interface, it should be mentioned that the distribution of lipid components observed in our simulations are slightly too wide compared with the neutron diffraction data of Jacobs and White (1989), as exemplified by the phosphate-headgroup densities, which are approximately 0.5 nm too far from the bilayer center. This is due to the small area per lipid observed in our simulations, which is 60 Å², compared with the experimental estimate of 68 Å² (Evans et al., 1987). However, as with most of our results, we are interested mainly in changes occurring as a result of peptide binding and in particular the partitioning of specific peptide groups in the interface. Thus, this discrepancy has little effect on our overall conclusions. Due to the large distance between sandostatin and the bilayer surface, there is no overlap of the peptide and lipid densities during the first 500 ps of the protonated run. These data are included to provide a source of comparison for the other plots and to confirm that the peptide is initially placed outside the interfacial region of the bilayer. In the densities calculated between 5.5 and 6 ns of the protonated run there is noticeable overlap between several side chains and the headgroups of the lipids. During this time period the Phe3 side chain penetrates the bilayer further than all other residues, having significant overlap with the PO^*₄Chol (see Fig. 1) headgroup density and to a lesser degree the C==O*Glyc and CH₂ regions. From the same plot it can be seen that the N-terminus of sandostatin is located near the outer edge of the bilayer at the PO^*₄Chol/water interface. In the first 500 ps of the deprotonated simulation, all side chains insert themselves more deeply into the bilayer than in the protonated complex. It is interesting to note, however, that the position of each side chain relative to the others along the bilayer normal remains very similar. Phe3, which still penetrates the interface the furthest, now has significant overlap with both the C==O*Glyc and hydrocarbon tail region of the lipids. This is also true of the Cys2, Cys7, and Phe1 side chains. The Lys side chain moves closer to the bilayer than Trp4 during the first 500 ps of the deprotonated simulation, despite the positive charges on the outer edge of the PO^*₄Chol density, whereas the tryptophan side chain has now moved to a position similar to that of Phe3 in the protonated simulation. In this position the indole group engages in significant hydrogen bonding with both the water and headgroup phosphates (data not shown), which partially prevents this side chain from penetrating the inner bilayer core.

View larger version (35K): [in this window] [in a new window]   FIGURE 8   Distribution of system components along the bilayer normal. (A and B) First and last 500 ps of the protonated trajectory, respectively. (C and D) First and last 500 ps of the deprotonated trajectory, respectively. The scale on the left and right of the graphs correspond to the lipid and peptide components, respectively. For each peptide residue, only the density of the side-chain atoms is shown.

The situation in the last 500 ps of the deprotonated run remains largely the same, except for a few significant differences. First, there is a sharpening of all the side chain densities, indicating that the location of the peptide along the bilayer normal has stabilized compared with the first 500 ps of the trajectory. Second, although all the features of the first 500 ps remain, there has been an overall movement of the peptide away from the bilayer core. Third, the side chain of Lys5 is once more located in the water phase.

Aromatic side-chain orientation

S(_N) and S(_L) parameters for the three aromatic residues of sandostatin as a function of simulation time are shown in Fig. 9. Concentrating first on the results for the protonated run, it is clear that there is significant orientational freedom for both Phe residues in the first 500 ps of the simulation, whereas the side chain of Trp4 appears to be slightly more rigid during this time. It should be noted, however, that the peptide as a whole rotates rather freely in this portion of the simulation. Although changes in the orientation of the side chains occur on a much faster time scale than overall rotation, changes in peptide orientation could effect the values of S_N and S_L. Directly after the first 500 ps, rapid and significant changes occur in the orientation of all three residues, particularly the phenylalanines. The long axes and planes of both these residues become rigidly perpendicular to the bilayer surface, whereas the parameters for Trp4 spike and then fluctuate over a wide range of values. This situation lasts up until 1.5 ns. Between ~1.5 and 2 ns, the orientations of all three residues change once again. Uniformly, all of the aromatic planes orient themselves more or less parallel to the bilayer surface. Consistent with its deep penetration of the bilayer, the plane of Phe3 appears to have the most rigid orientation of the three residues. In comparison, the planes of Phe1 and Trp4, which are located further out from the bilayer interior toward the water phase, fluctuate rather significantly about their orientations. The long axes of all three residues lie parallel to the plane of the bilayer for the most part but are reasonably free to move, particularly between 2.5 and 5.5 ns. The orientation of both Phe residues indicate changes occurring in the protonated complex in the last nanosecond of the simulation. In particular, the plane of Phe1 undergoes rotation around the long axis at approximately 5.5 ns. For the remainder of the simulation the plane of the ring remains perpendicular to the bilayer surface. Slightly different effects are seen in the orientation of the Phe3 side chain, with an overall movement of both the plane and long axis from a perpendicular to a parallel orientation. Conversely, the orientation of Trp4 remains more or less stable during this time.

View larger version (53K): [in this window] [in a new window]   FIGURE 9   Orientation of the aromatic peptide side chains from the protonated and deprotonated trajectories. Solid line, SN; dashed line, SL. Running averages over 50 ps are shown for clarity.

The overall picture for all three residues during the deprotonated run is similar to that between 2.5 and 5.5 ns of the protonated simulation. Once again, all of the aromatic planes and long axes lay parallel to the bilayer surface. Aside from some very rapid and transient changes occurring in all residues at the beginning and/or end of the simulation, no appreciable changes in orientation are observed. As in the protonated run, these changes occur more or less simultaneously in all three residues. All three residues experience about the same freedom in the orientation of their respective planes, whereas the long axis of Phe1 seems to be much more rigidly oriented than those of the Phe3 and Trp4.

Phospholipid headgroup orientation

Through geometrical interpretation of quadrupole splittings from POPC lipids with selectively deuterated headgroups, it has been shown that the binding of positively charged species to bilayers causes the lipid headgroups to protrude out into the water phase (Seelig et al., 1987; Scherer and Seelig, 1989; Beschiaschvili and Seelig, 1991). In short, the quadrupole splitting of the headgroup -methylene decreases in a linear fashion with respect to the mole fraction of bound, positively charged amphiphile, whereas that of the -methylene increases. Negatively charged amphiphiles have the opposite effect in both cases. For positive species the changes in the splittings can be interpreted as a rotation of the P-N⁺ headgroup dipole away from the surface of the bilayer into the aqueous phase. In the absence of charged binding species, it is estimated that these dipoles lie within 30° of the bilayer surface in a more or less parallel orientation (Scherer and Seelig, 1989). To determine the effect of sandostatin on the lipid headgroups, we first divided the simulations into windows of 500 ps and counted the number of times that any P-N⁺ dipole in the bilayer is significantly aligned with the z axis, with the N⁺ end pointing out into the water phase. Due to their interaction with the peptide, the headgroups of the upper leaflet will presumably undergo more significant changes in orientation than those of the lower leaflet. Fig. 10 shows the results of the calculation for the protonated run. Comparing the plot for the upper leaflet with that of the lower leaflet reveals that, in fact, the P-N⁺ dipoles of the latter point out from the bilayer surface more often than those of the former in the first 1.5 ns. It should be noted that a similar calculation for the BW simulation revealed similar differences between the orientation of the headgroups in the two leaflets (data not shown), and it is therefore unlikely that the peptide is responsible for these initial differences in the protonated simulation. Between 1.5 and 2.0 ns, the number of occurrences in the upper leaflet increases drastically until the value for both leaflets is approximately equal. The end of this time frame coincides with the peptide reaching a stable minimum distance from the upper bilayer surface (see Fig. 3). In the next 2.0 ns, the value for the upper leaflet experiences another drastic increase, reaching a maximum between 2.5 and 3.0 ns. In contrast, the value for the lower leaflet remains fairly constant and close to the values observed for the BW simulation (data not shown). Past 3.0 ns in the simulation, the effect of the peptide on the upper leaflet headgroups relaxes to a value similar to that of the lower leaflet. Interestingly, there was very little effect on the orientation of headgroups in either leaflet during the deprotonated simulation; for the entire 6 ns, values for the upper and lower leaflets remained fairly constant and similar to those of the BW simulation (data not shown).

View larger version (17K): [in this window] [in a new window]   FIGURE 10   (A) Orientation of the lipid headgroups as a function of simulation time. For every 500-ps window of the protonated trajectory, the number of occurrences in which any P-N+ vector is significantly oriented along the z axis and pointing out from the bilayer is counted. To be counted, the cosine of the angle between the vector and the z axis had to be greater than or equal to 0.9. (B) Orientation of the lipid headgroups with respect to distance from the peptide (see text for details). Bin 1 starts at a value of 0. The width of each bin is 0.3 nm.

Given the headgroup orientation effect observed in the protonated simulation, two possibilities exist concerning the spatial distribution of the phenomenon: 1) sandostatin has a small, essentially uniform effect on all headgroups, or 2) the effect is dependent on the distance between the headgroup and the peptide. To investigate these possibilities, we first divided the lipids of the upper leaflet into bins based on the minimum distance between the P-N⁺ headgroup vector and any peptide atom. For each bin, we then calculated the average orientation (as the cosine of the angle formed between the P-N⁺ vector and the bilayer normal) between 1.5 and 3.0 ns. The results are shown in the lower plot of Fig. 10 and reveal that the lipids within 0.3 nm of the peptide lay almost perfectly flat against the bilayer surface. Radiating outward to a distance of 0.3-1.2 nm from the peptide, the alignment of headgroups along the bilayer normal increases slightly and is fairly constant within this range. Between 1.2 and 1.5 nm, however, the degree to which the headgroups protrude from the bilayer surface increases significantly. The localization of the effect at this distance proves to be quite sharp, with the headgroups lying more or less flat at greater distances.

That transient and spatially localized nature of the headgroup orientation effect in the protonated simulation differs from the experimental results. When measuring the quadrupole splittings of selectively deuterated carbon atoms in the POPC headgroups, single, well-defined spectra are obtained for a given concentration of bound peptide, indicating that sandostatin undergoes a rapid two-dimensional diffusion in the plane of the bilayer and thus affects all headgroups uniformly within the order of 10⁵ s (Beschiaschvili and Seelig, 1991). Obviously, the time scale of this diffusion lies far outside the reach of our simulations and may serve to explain the transient nature of the effect observed in the simulation. Also, the time scales of lipid motions (such as the overall rotation about the bilayer normal) are such that the sampling over states required for the accurate calculation of some NMR observable data is not possible from a nanosecond-scale simulation. Nevertheless, we have calculated quadrupole splittings for both the - and -methylenes of the headgroups. Even a qualitative agreement between the calculated and experimental values would help to confirm that the headgroup protrusion observed in the upper plot of Fig. 10 occurs as a result of the same conformational changes inferred from the experimentally observed quadrupolar splittings (Scherer and Seelig, 1989). For calculation of the splittings from the simulation, the following relationship was used:
which is valid for lipids in a bilayer oriented with the normal parallel to the external magnetic field, provided the asymmetry parameter of the electrostatic field gradient is small and can be ignored, which is commonly the case (Seelig, 1977). S_zz is the relevant element of the order parameter tensor, which can be calculated for each methylene directly from the simulations, and the collection of constants in brackets is the deuteron quadrupole coupling constant, for which a value of 170 kHz was used (Seelig, 1977). For the protonated simulation, the average splittings of the upper leaflet headgroups between 1.5 and 3 ns (the time period of reorientation) were calculated. These were compared with the values obtained from the same leaflet in the last 1.5 ns of the BW simulation. The calculated quadrupolar splitting constants for the -methylenes in the protonated and BW simulations were 46.4 and 36.8 kHz, respectively. Both of these values lie far outside the experimentally determined range for POPC bilayers with various charged amphiphiles, which is (+ charge) 11.6 kHz    +10 kHz ( charge) and which varies with the charge and amount of amphiphile bound. A similar situation is seen for the - methylenes: values of 46.4 and 34.1 were calculated for the protonated and BW simulations, respectively, with the experimental range being (+ charge) 14.4 kHz     2 kHz ( charge). Despite the discrepancies between the calculated values and the experimentally determined ranges, which is due to the limited time scale of the simulation, the magnitude of the relative uncertainty in each of the calculated values (calculated from the RMS values in a manner analogous to that of the lipid chain order parameters; see below) is less than 5% in each case, making it possible to calculate statistically significant differences between the values for the protonated and BW simulations. These differences are  = 9.6 and  = 12.4, respectively. Thus, the influence of sandostatin on the quadrupolar splitting of each headgroup segment is qualitatively reproduced in the simulations, with decreasing (becoming more negative) and increasing upon binding of the peptide to the bilayer.

Lipid order parameters

For both simulations, deuterium order parameters were initially calculated by averaging over all lipids in successive time windows of 2 ns each (data not shown). Although the shape of the order parameter profiles thus obtained were in excellent agreement with those obtained from experiment (Lafleur et al., 1990; Seelig and Seelig, 1977; Seelig, 1977) and other simulation studies (Heller et al., 1993; Damodaran et al., 1992; Egberts and Berendsen, 1988), the values themselves were uniformly too high. As with the lipid density distributions, this can be attributed to the small area per lipid we observe in the simulations. This effect of the area per lipid on the order parameters has been discussed by Feller et al. (1996). When compared with profiles obtained from the BW simulation, the order parameters calculated for the protonated and deprotonated complexes using all lipids showed small differences in the palmityl chains; in the protonated complex, the presence of the peptide raised the order parameters near the top of the chain (carbons 2-5) and lowered those in the rest of the chain. In the deprotonated complex, the order parameters near the top of the chain remained close to those of the BW simulation, although the rest of the chain seemed to be disordered by the peptide. Effects observed for the oleoyl side chain were even less noticeable. To examine the possibility of more significant localized changes in lipid order, we carried out an analysis similar to that of Tieleman et al. (1998) in which the carbons were divided into bins based on their distance from any peptide atom. Order parameters for each carbon were then calculated for each bin. As with their results for single transmembrane ()-helices, we found that a disordering effect was observed only for lipids within the first layer of lipids (those within 0.8 nm) surrounding the peptide (data not shown). However, due to the fact that the peptide remains relatively close to the bilayer surface in our simulations, carbons near the terminal end of the chains were rarely placed in the first bin, resulting in poor statistics. To circumvent this, we selected a subset of lipids that made significant van der Waals contacts with the peptide (minimum distance between any peptide and lipid atom was less than 0.6 nm), and calculated order parameters for all atoms in both tails. Once again, successive time windows of 2 ns were used for both the protonated and deprotonated simulations. Fig. 11 shows the order parameters calculated for these nearby lipids as well as the RMS values in the first and last 2 ns of each simulation, which were calculated according to
where S_iL is the instantaneous value of the order parameter for atom (i) in lipid L, and S_iL is the average order parameter for a given atom. N_L and N_F are the number of lipids and coordinate frames used in the averages, respectively (Damodaran et al., 1995). Order parameters calculated from the BW simulation using all lipids are shown in the graphs for comparison. A cursory look at the data seems to show a general trend: in the beginning of both simulations, the order parameters of both chains are significantly and uniformly greater than the BW values. However, as the simulation progresses, order parameters calculated from subsequent time windows become progressively smaller until they are similar or lower than the BW values. In the case of the protonated run, significant lowering of the order parameters occurs near the end of the palmityl chain, with a somewhat smaller effect in the oleoyl chain.

View larger version (37K): [in this window] [in a new window]   FIGURE 11   Deuterium order parameters for lipid chains within van der Waals contact (0.6 nm) of the peptide. Six lipids were included in the calculation for the protonated trajectory (upper plots), and seven were included for the deprotonated trajectory (lower plots). Order parameters were calculated in successive time windows of 2 ns each. Carbon atoms are numbered relative to the ester carbon, which is considered to be carbon 1 of the chain. RMS values for the order parameters calculated for the first and last 2 ns of each simulation are shown as error bars.

For the palmityl chains in the last 2 ns of the deprotonated run, there is an overall lowering of the order parameters along the bulk of the chain as compared with the first 2 ns. Although the values calculated between 2 and 4 ns for this chain cannot be considered significantly lower than those of the first 2 ns, they are roughly intermediate between those from the beginning and end of the simulations and thus seem to indicate an overall trend as the peptide exerts its influences on the lipid chains. In contrast to the results for the protonated simulation, the order parameters near the headgroup rather than terminal region of the oleoyl chains are most affected by sandostatin in the deprotonated run. As with the values calculated using all lipids (see above) the order parameters for all atoms in the nearby lipids are too high at the beginning of the simulations. Once again, this is likely due to the small area per lipid in our simulations. However, if we restrict our attention to relative changes in lipid order as a result of peptide binding, it is clear that the peptide is having a disordering effect.

Interfacial water ordering

Fig. 12 shows the ordering of water molecules along the z axis for the BW, protonated, and deprotonated simulations. The water orientation profiles in the plots differs from those of other simulation analyses in that the orientation of the water molecules with respect to normals pointing out from each bilayer leaflet (as opposed to the z axis of the simulation box) is shown. This merely represents a different frame of reference for the ordering of the water molecules and does not influence the conclusions. Although the use of short electrostatic cutoffs has been shown to induce excessive ordering of water at interfaces (Feller et al., 1996), we do not observe such an effect, and our cutoff value of 1.7 nm seems sufficient to produce good agreement with data obtained in other simulation studies (Zhou and Schulten, 1995; Marrink et al., 1993). As with the order parameters, we are concerned here only with relative changes in the water ordering caused by the binding of the peptide. The characteristic features of the plot include a large positive peak near the top of the carbon chains and a negative peak of approximately the same magnitude further out among the headgroups, representing parallel and antiparallel alignments of the water molecules with the bilayer normal, respectively. Despite the presence of the peptide, the profile remains largely undisturbed during the first 2 ns of the protonated simulation. However, the positive peaks decrease significantly following peptide adsorption, indicating a loss or water ordering near the interface. In the deprotonated simulation, the height and width of these peaks increase again and are only slightly smaller than that of the BW simulation and the first 2 ns of the protonated simulation.

View larger version (24K): [in this window] [in a new window]   FIGURE 12   Orientation of the water molecules with respect to normals on the bilayer leaflets. The simulation box was divided into 50 slices perpendicular to the z dimension and the average angle between the water dipoles and the leaflet normals were calculated (solid line). Where applicable, the dashed line shows the location of the peptide density. (A) BW simulation; (B) First 2 ns of the protonated trajectory; (C) Last 4 ns of the protonated trajectory; (D) Deprotonated trajectory.

	DISCUSSION

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The binding process

As mentioned in the Introduction, the initial location of sandostatin in the protonated run was chosen so as to obtain a spontaneous peptide-bilayer complex that was relatively independent of the starting conditions. The peptide-bilayer distance, orientation, and RMSF (Figs. 3, 4, and 5, respectively) data from the first 2 ns of the simulation attest to the translational, rotational, and conformational freedom of the peptide during this time.

Intermolecular forces

The fluctuation of the intermolecular forces around zero at the beginning of the protonated simulation show that the peptide feels little influence from the bilayer initially. This confirms that the starting conditions of the simulation allow the peptide to adopt a natural structure and orientation at the interface. Despite the positive charge and neutrality of sandostatin and the lipids, respectively, the attractive force existing between the two species at distances less than 1 nm is a key determinant of the binding process. In contrast, the force on the peptide due to the solvent directly opposes the binding and is probably a frictional effect due to the displacement of water molecules from between the two species. We feel that this is a direct observation of the much discussed dehydration force (Israelachvili and Wennerström, 1992, 1996; Marrink et al., 1993). As a frictional effect experienced during the binding process, this force is not to be confused with the hydrophobic effect, which will thermodynamically favor the binding of peptides to bilayers. There is an effective disappearance of the solvent force in the last nanosecond of the protonated run and in the deprotonated run. This, together with the consistently attractive bilayer force, demonstrates the preference of the peptide for the lipid phase over that of the aqueous.

Buried surface area

Simple calculations of a molecule or complex often assume a linear relationship between the free energy and the solvent-exposed surface area (Juffer et al., 1995). In their calculations of peptide-bilayer binding free energies, Ben-Tal et al. (1997) further assumed that the change in exposed surface area was linear with respect to the distance between the two species. Based on the linear increase in the area of the peptide-bilayer interface between 1 and 2 ns of the protonated simulation (see Fig. 6) and the linear decrease in the peptide-bilayer distance during this time, we have observed just such a relationship. This result contributes to the validity of the approach and indicates that it could be of general utility in continuum calculations.

Peptide-bilayer complexes

Despite the fact that no procedure was employed to position sandostatin favorably on the bilayer surface, details of both simulated complexes are in excellent agreement with experimental findings of Beschiaschvili and Seelig (1990, 1991, 1992). Taken together, the protonated and deprotonated simulations may closely resemble the process by which sandostatin immerses itself into PC lipids. From calorimetric measurements, Beschiaschvili and Seelig (1992) determined that 0.24 protons were released per peptide during the binding process. Based on standard ±K_a values, this deprotonation most likely occurs at the N-terminus. Thus, as the peptide approaches the bilayer surface, the change in local environment (due to the lipid headgroups and ionic surroundings, the latter of which is absent in our simulations) causes a loss of protons, allowing the peptide to partition more efficiently into the lipid phase. It was not necessary to remove lipid molecules from the bilayer to make room for the deprotonated peptide, indicating that the positive charge on the N-terminus presented the primary barrier to deeper penetration. The location of the N-terminus peak in the density plots supports this conclusion. Similar relationships between peptide charge and the degree of peptide-bilayer interaction have also been observed in other simulations (Bernèche et al., 1998; Damodaran and Merz, 1995) and may indicate a very general effect.

Using a simple binding isotherm that related the degree of binding to surface areas of the peptide and lipids, Beschiaschvili and Seelig (1990) also concluded that the bound peptide presents one of its long faces to the bilayer surface. This is exactly what we observe in both simulations. Moreover, a large rotation was necessary for the peptide to adopt this orientation upon adsorption (See Fig. 6), indicating that it is favored. Assuming that the buried surface area in the complexes results from a perfectly matched interface between the peptide and lipids, the percentage of sandostatin's buried surface area can be estimated, which results in values of 30% and 57% for the protonated and deprotonated complexes, respectively. The latter value is close to the experimental estimate of 60-70% (Beschiaschvili and Seelig, 1990).

Although Beschiaschvili and Seelig (1991) measured the secondary structure of sandostatin bound to charged POPG lipids only, the structural changes in the simulated complexes (relative to the free peptide in the SW simulation) corresponds nicely with the increase in -turn observed from their circular dichroism measurements. As mentioned earlier, they assumed that the conformational entropy loss in the peptide could be effectively ignored in their thermodynamic analyses (Beschiaschvili and Seelig, 1992). Whereas we observe a relative decrease in the conformational freedom of sandostatin following binding and upon deprotonation, it is likely that this is a reasonable assumption for this rigid cyclic peptide. To put the data of Fig. 5 in perspective, Tieleman et al. (1999a,c) made a similar calculation for a 20-residue alamethicin peptide inserted into a bilayer. The RMSF values they calculated for C atoms alone were on the same order of magnitude as those we have obtained for all atoms of sandostatin in both complexes. Nonetheless, given that conformational effects contribute directly to the free energy of peptide-bilayer binding and insertion (Jähnig, 1983; Ben-Shaul et al., 1996), direct calculation of conformational entropy from simulations (Tidor and Karplus, 1994) of rigid and flexible peptides would prove useful in deciding the extent to which they must be considered quantitatively.

Correspondence between the simulated complexes and experimental measurements are not limited to the peptide alone. The effect of the peptide on the orientation of the lipid headgroups was also reproduced in the protonated simulation. As observed in NMR experiments on POPC lipids with selectively perdeuterated methylene groups in the choline moiety (Beschiaschvili and Seelig, 1991), the binding of the peptide in our simulation causes the headgroups to stick out from the bilayer surface into the water phase. Interestingly, this effect was seen only in the protonated simulation and not in the deprotonated run. Moreover, the headgroups relaxed back into a flat orientation shortly after adsorption. This latter phenomenon could be attributed to the fact that the time scale of the simulation is insufficient to observe the rapid, two-dimensional diffusion of the peptide along the bilayer surface. Our result suggests that this rapid two-dimensional diffusion is necessary to sustain the effect, which is caused by peptides localized to the very surface of the bilayer. Also, there is an apparent localization of the effect to lipids within 1.2-1.5 nm away from the peptide (see Fig. 10). It is interesting to note that the protrusion of headgroups from bilayer surfaces may have biological significance: as two bilayer surfaces approach each other, the conformational restriction of these protrusions contribute to a repulsive force between the two membranes, which is entropic in origin (Israelachvili and Wennerström, 1992, 1996; Marrink et al., 1993). Because the binding of peptides to the bilayer surface increases the number of these protrusions, this contribution to the repulsive force would also be expected to increase. From their simulation studies of fusion-inhibition peptides, Damodaran and Merz (1995) suggested that effects of the peptide on the lipid chain order may provide the inhibition mechanism. Based on the experimental observation of the headgroup orientation effect (Seelig et al., 1987; Scherer and Seelig, 1989; Beschiaschvili and Seelig, 1991), the results of our simulations and the theory of the short-range repulsive force between bilayers (Israelachvili and Wennerström, 1992, 1996), we suggest that headgroup orientation could also play a significant role.

Given the paucity of experimental information on the orientation of aromatic residues in membrane-associated peptides, it is difficult to interpret our own results in a general context. In their measurements of side-chain orientation for a number of protein-bilayer systems, Tieleman et al. (1998) found that phenylalanine side chains enjoyed a relatively high degree of flexibility in a wide range of orientations. Tyrosine and tryptophan side chains were more rigid, with individual orientations being determined largely by packing requirements and hydrogen bonding interactions with the lipids and solvent (Tieleman et al., 1998). Although sandostatin contains only three aromatic residues, all of these adopt a similar, parallel orientation with respect to the bilayer surface. Although the apparent stability of these orientations could be an artifact of limited simulation length, Tieleman et al. (1998) observed significant flexibility in a number of residues (including those located in the hydrocarbon core) within shorter or equivalent time scales. Also, the few chang