help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Sankararamakrishnan, R.
Right arrow Articles by Weinstein, H.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Sankararamakrishnan, R.
Right arrow Articles by Weinstein, H.

Biophys J, November 2000, p. 2331-2344, Vol. 79, No. 5

Molecular Dynamics Simulations Predict a Tilted Orientation for the Helical Region of Dynorphin A(1-17) in Dimyristoylphosphatidylcholine Bilayers

Ramasubbu Sankararamakrishnan and Harel Weinstein

Department of Physiology and Biophysics, Mount Sinai School of Medicine, New York University, New York, New York 10029 USA


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

The structural properties of the endogenous opioid peptide dynorphin A(1-17) (DynA), a potential analgesic, were studied with molecular dynamics simulations in dimyristoylphosphatidylcholine bilayers. Starting with the known NMR structure of the peptide in dodecylphosphocholine micelles, the N-terminal helical segment of DynA (encompassing residues 1-10) was initially inserted in the bilayer in a perpendicular orientation with respect to the membrane plane. Parallel simulations were carried out from two starting structures, systems A and B, that differ by 4 Å in the vertical positioning of the peptide helix. The complex consisted of ~26,400 atoms (dynorphin + 86 lipids + ~5300 waters). After >2 ns of simulation, which included >1 ns of equilibration, the orientation of the helical segment of DynA had undergone a transition from parallel to tilted with respect to the bilayer normal in both the A and B systems. When the helix axis achieved a ~50° angle with the bilayer normal, it remained stable for the next 1 ns of simulation. The two simulations with different starting points converged to the same final structure, with the helix inserted in the bilayer throughout the simulations. Analysis shows that the tilted orientation adopted by the N-terminal helix is due to specific interactions of residues in the DynA sequence with phospholipid headgroups, water, and the hydrocarbon chains. Key elements are the "snorkel model"-type interactions of arginine side chains, the stabilization of the N-terminal hydrophobic sequence in the lipid environment, and the specific interactions of the first residue, Tyr. Water penetration within the bilayer is facilitated by the immersed DynA, but it is not uniform around the surface of the helix. Many water molecules surround the arginine side chains, while water penetration near the helical surface formed by hydrophobic residues is negligible. A mechanism of receptor interaction is proposed for DynA, involving the tilted orientation observed from these simulations of the peptide in the lipid bilayer.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

Endogenous opioid compounds are involved in the perception of pain and in the modulation of behavior and neuroendocrine function; they exert their actions through specialized opioid receptors in the class of G protein-coupled receptors (Brownstein, 1993; Williams et al., 1999; Stein, 1999). Among the three major types of opioid receptors (µ, delta , and kappa ), the ligands of the kappa -receptor are believed to have low abuse potential and to cause milder forms of dependence compared to the narcotic µ-opioid ligand morphine (Millan, 1990; Naqvi et al., 1998). The peptide dynorphin A is an endogenous ligand selective for the kappa -opioid receptor (Chavkin et al., 1982; Chavkin and Goldstein, 1981), and its potential as an analgesic has made it an attractive target for research since its discovery more than two decades ago (Cox et al., 1975). The structures of opioid ligands such as dynorphin have been studied with spectroscopic methods in various solvents (Saviano et al., 1999; Segawa et al., 1995; Yan et al., 1999; Tessmer and Kallick, 1997). The naturally occurring Dynorphin A has 17 amino acids with the sequence H-Tyr-Gly-Gly-Phe-Leu-Arg-Arg-Ile-Arg-Pro-Lys-Leu-Lys-Trp-Asp-Asn-Gln-OH (Goldstein et al., 1981). It has been shown that the peptide incorporating the first 13 residues of the natural peptide dynorphin (Dynorphin A-(1-13)) has practically the same pharmacological profile as its parent peptide (Chavkin and Goldstein, 1981). The smaller N-terminal fragments (the first eight and nine residues of dynorphin A) were also shown to be selective ligands for the kappa -binding site (Corbett et al., 1982). A variety of biochemical and pharmacological studies have identified residues that play important roles in the activity and/or potency of the dynorphin peptides (Snyder et al., 1992; Gairin et al., 1988; Nakajima et al., 1988; Chavkin and Goldstein, 1981; Turcotte et al., 1984). However, an understanding of dynorphin-kappa -receptor interactions at the atomic level is still elusive, particularly because opioid receptors, being integral membrane proteins, have resisted structural characterization.

It has been proposed that the first step in the mechanism of action of peptide hormones would be their accumulation in the lipid bilayer, and that the subsequent interactions between the membrane and membrane-induced ligand conformation of the peptide may determine ligand-receptor interactions (Schwyzer, 1991, 1995). NMR studies of the lipoderivative of the cholecystokinin peptide hormone have shown evidence for the preadsorption of peptide hormones at the cell membrane bilayer before their interaction with their receptors (Moroder et al., 1993). More recently, it has been shown that sodium dodecyl sulfate micelles induce different secondary structure types for agonists and antagonists of mammalian tachykinin NK1 receptor (Whitehead et al., 1998). The conformational differences have been correlated with the binding potencies and biological activity of the peptides. Dynorphin A-(1-13) was found to assume a completely random or extended conformation in water (Renugopalakrishnan et al., 1988) and methanolic solution (Lancaster et al., 1991), respectively. However, dynorphin A-(1-13) was shown by vesicle-mediated hydrophobic labeling to interact with anionic liposomes (Gysin and Schwyzer, 1983). Infrared attenuated total reflection spectroscopy and capacitance minimization experiments on this peptide in neutral lipid membranes suggested a model in which dynorphin A-(1-13) assumed a helical structure from residues 1 to 10 (Erne et al., 1985). Based on these studies and amphiphilic moment calculations (Schwyzer, 1986), Schwyzer suggested that the more hydrophobic N-terminal helical segment is oriented perpendicular to the membrane surface, contacting the hydrophobic membrane layers, whereas the extended C-terminal segment would be in contact with the aqueous phase (Schwyzer, 1995). Recent NMR studies (Tessmer and Kallick, 1997) have shown the existence of an alpha -helical region from residues 3 through 9 when dynorphin A-(1-17) is bound to dodecylphosphocholine (DPC) micelles. However, the structural and dynamic details of the interactions were not elucidated.

Studies of peptides in explicit bilayers with molecular dynamics (MD) simulations have produced details related to various biological mechanisms. Examples include the characterization of gramicidin channels (Woolf and Roux, 1996; Chiu et al., 1999b,c), the two-stage model of membrane protein folding in bacteriorhodopsin helices (Woolf, 1997, 1998a,b), the mechanism of membrane lysis caused by melittin (Berneche et al., 1998; Bachar and Becker, 1999), the action of small fusion-inhibiting peptides (Damodaran and Merz, 1996), and the stability of the channel formed by alamethicin helix bundles (Tieleman et al., 1999). The importance of residues bridging the transmembrane and surface-bound helical segments in bacteriophage Pf1 coat protein (Roux and Woolf, 1996) and the role of the Asn-Pro motif in the seventh transmembrane segment of 5HT2A receptor (Ha Duong et al., 1999) have also been described from such detailed MD simulations. In this study, we have carried out a multinanosecond MD simulation of dynorphin A-(1-17) in dimyristoylphosphatidylcholine (DMPC) bilayers. The NMR structure obtained by Tessmer and Kallick (1997) was used in the simulations. In the initial structure, the dynorphin N-terminal helix was placed inside the DMPC bilayers oriented perpendicular to the membrane, as suggested by Schwyzer (Erne et al., 1985; Schwyzer, 1995). The C-terminal region lay approximately parallel to the membrane surface. Two simulations were carried out in which the initial positioning of the dynorphin helical segment differed by 4 Å from the center of the bilayer. The protocol includes a long equilibration (1.4 to 1.6 ns) followed by a ~2-3-ns production run. Notably, the structures from the two simulations converged, even though the starting structures were different, and the helical segment remained within the bilayer throughout the simulations and in the final 1 ns of the simulations maintained an angle of 50° with the bilayer normal. Detailed analysis of the trajectories reveals the modes of interaction of individual dynorphin residues with the lipid and water, the preferential direction of water penetration in the lipid bilayer, and the role of primary and secondary amphiphilicity of dynorphin in the complex bilayer environment. These findings suggest a role for membrane insertion in the mechanism of interaction of DynA with the opioid receptor.


    MATERIALS AND METHODS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

The simulations were carried out with the all-atom PARAM 22 force field (Mackerell et al., 1998) of CHARMM (Brooks et al., 1983), which includes phospholipids (Schlenkrich et al., 1996) and TIP3P water potentials (Jorgensen et al., 1983). The nonbonded list was generated using a group-based cutoff of 12 Å. Both the electrostatic and van der Waals interactions were smoothly switched from 8 to 11 Å. The SHAKE (Ryckaert et al., 1977) algorithm was used to fix the length of all bonds involving hydrogen atoms. The temperature of the system was set to 330 K, above the gel-liquid phase transition of DMPC (Gennis, 1989). The trajectories were calculated in the microcanonical ensemble with a constant number of particles, volume, and energy (NVE). Although constant pressure simulations are desirable in bilayer simulations, a recent 10-ns constant volume simulation on DPPC bilayers has shown that the system remained stable throughout the simulation (Essmann and Berkowitz, 1999). Many of the peptide bilayer simulations have been carried out successfully in the NVE ensemble (Berneche et al., 1998; Woolf, 1997; Damodaran and Merz, 1996; Ha Duong et al., 1999; Huang and Loew, 1995).

Initial structures

The NMR structure of dynorphin obtained in DPC micelles consisted of an alpha -helical segment (residues 3-9) in the N-terminal region and a type I or type IV beta -turn (residues 14-17) in the C-terminal region (Tessmer and Kallick, 1997). The residues connecting these two segments (residues 10-13) were in random conformation. The initial structure of dynorphin for the present simulations was generated from the internal parameters obtained from these NMR studies. Although residues 1 and 2 were not part of the alpha -helical segment in the NMR structure, the N-terminal helical segment was extended in the initial structure to include these two residues. A type I beta -turn was constructed for the residues 14-17 in the C-terminal end. Two patch residues were created to block the N and C termini with NH2 and COOH groups, respectively.

The initial peptide and hydrated lipid bilayer system was constructed using a protocol that was developed by B. Roux and his collaborators and was described recently (Woolf and Roux, 1996; Berneche et al., 1998; Roux and Woolf, 1996); a brief description is given below in relation to features specific to our system. The nine-residue dynorphin N-terminal helical segment is too short to traverse the membrane fully and will occupy only one-half of the lipid bilayer. This introduces asymmetry into the system and results in a different number of lipids in the top and bottom layers (such an asymmetrical situation occurred for melittin and Pf1 coat protein; Berneche et al., 1998; Roux and Woolf, 1996). Construction of the NVE ensemble requires a reasonable estimate of the cross-sectional area of the system. The value considered for the DMPC cross-sectional area ranged from 61 to 66 Å2 in previous simulation studies (Chiu et al., 1999b; Woolf, 1997; Shen et al., 1997; Damodaran and Merz, 1996; Essmann and Berkowitz, 1999). In the present simulations, we have taken the cross-sectional area of DMPC to be 63.1 Å2. Recently, x-ray diffraction studies have shown that the cross-sectional area for DMPC is 59.7 ± 0.2 Å2 at 30°C (Petrache et al., 1998). Dipalmitoylphosphatidylcholine (DPPC), the headgroup of which is the same as that of DMPC, is found to have a cross-sectional area of 62.9 ± 1.3 Å2 at 50°C (Nagle et al., 1996). The area of DMPC is expected to be larger at higher temperatures, as discussed by Nagle (Petrache et al., 1998; Nagle et al., 1996), and so at a simulation temperature of 330 K (57°C) our value of 63.1 Å2 for DMPC is reasonable. The cross section of dynorphin varied from 120 Å2 to a maximum of 250 Å2. The presence of larger Arg side chains contributed to the cross-sectional area, which is significantly larger than that of poly-Ala (133 Å2) (Shen et al., 1997).

The periodic system consists of a rectangular box with the dimensions X = 53.3 Å, Y = 53.3 Å, and Z = 90 Å. The membrane normal is oriented along the z axis, and the center of the bilayer is at Z = 0 Å. The maximum dynorphin cross-sectional area (250 Å2) is equivalent to approximately four DMPC lipids. In our system, the dynorphin was placed in the upper layer (Z = +ve) with 41 lipids, while the bottom layer consisted of 45 lipids. Individual DMPC molecules were randomly chosen from a library of 2000 preequilibrated (Venable et al., 1993) and prehydrated (Woolf and Roux, 1994) DMPC lipids. The prehydrated lipid molecules included ~20 water molecules around the phosphate and the choline groups of each DMPC. Initial positioning of DMPC headgroups within Z = ±17 Å and removal of bad contacts between individual molecules (lipids, peptide, and water) were achieved as described in previous studies (Woolf and Roux, 1996; Roux and Woolf, 1996; Berneche et al., 1998). The remaining bulk solvent was constructed from a slab of 1474 water molecules equilibrated using the same XY periodic boundary conditions. The water box was translated along the z axis, and its position was adjusted to give the correct total number of waters for the system. Waters that projected into the hydrocarbon interior (Z = ±14 Å) were deleted. The number of water molecules (~62 waters per lipid) is significantly larger to make sure that the C-terminal segment of dynorphin remains fully solvated throughout the simulations within the primary box. With ~5300 water molecules, the final system consisted of a total of ~26,400 atoms.

The N-terminal helical segment of dynorphin was initially oriented perpendicular to the membrane plane within the bilayer. Two simulations were carried out with the center of the Calpha atoms of residues 1-10 of dynorphin initially placed at Z = 14 Å (system A) and Z = 10 Å (system B). In both the cases, the C-terminal residues 11-17 were approximately parallel to the bilayer plane. The systems were refined by energy minimization before the MD simulations started. Periodic boundary conditions were applied in all directions during the MD simulations.

Simulation details

To ensure a smooth relaxation of the system toward an equilibrated configuration, harmonic and positional restraints were imposed initially on the selected lipid and peptide atoms. These restraints were gradually relaxed during the equilibration. A time step of 0.002 ps was used. To converge to an equilibrium state, the systems were first coupled to a heat bath at 330 K, and Langevin dynamics simulations of 0.1 ns were carried out. A planar harmonic restraint was applied at the center of mass of the lipid atoms to maintain the planarity of the membrane. Positional harmonic restraints were applied on all Calpha atoms of dynorphin in system A. A cylindrical harmonic restraint was applied on the center of mass of alpha -helical residues 1-9 in system B. Although the restraints on dynorphin were different in the two simulations, they were applied temporarily, and mainly to maintain the alpha -helical conformation of the dynorphin N-terminal segment. The force constants on all of the restraints were gradually reduced during the equilibration, which took 0.9 ns for system B and 1.1 ns for system A. Subsequently, an additional 0.5 ns of equilibration was carried out for each system in the absence of any restraints. During this period, the total energy was in equilibrium and velocity rescaling was observed only occasionally. The structures present at the end of equilibration are shown in Fig. 1.



View larger version (108K):
[in this window]
[in a new window]
 
FIGURE 1   Dynorphin-lipid complex at the end of equilibration runs of 1.6 ns (system A, left) and 1.4 ns (system B, right), respectively. Only nonhydrogen atoms are shown. Lipid headgroup atoms are color-coded yellow (P), red (O), and cyan (N). Lipid acyl chains are green, and dynorphin is in pink. Water oxygen atoms are blue. The helical region of dynorphin is within the lipid bilayers and is oriented approximately parallel to the bilayer normal.

For system A, a further 3-ns MD simulation constituted the production trajectory, whereas a further 2-ns MD simulation was carried out as the production trajectory for system B. The coordinates were saved every 0.1 ps for analysis. The systems were completely free of any restraints during the production runs. The simulations were run on SGI Origin 200 machines with four processors using the parallel version of CHARMM. The computer time required was ~40 min/ps. The total energy, its decomposition into kinetic and potential energies, and the temperature of the systems were monitored throughout the simulations.

Interaction energies

Interaction energies between bacteriorhodopsin helices and DMPC lipids have been used to determine the preferences of individual helical residues for different regions of the lipid bilayers (Woolf, 1998b). They have also been used to study the water penetration within the bilayer when melittin was present (Bachar and Becker, 1999). The interaction energy (Eint (XY)) between any two groups X and Y is calculated as
E<SUB><UP>int</UP>(<UP>XY</UP>)</SUB>=E<SUB><UP>total</UP>(<UP>XY</UP>)</SUB>−(E<SUB><UP>X</UP></SUB>+E<SUB><UP>Y</UP></SUB>) (1)
Where Etotal (XY) is the total potential energy consisting of groups X and Y, and EX and EY are self-energies of X and Y, respectively. For the peptide-bilayer system, the potential energy can be decomposed further into self-energies of peptide, lipid, and water and interaction energies of peptide-lipid, peptide-water, and lipid-water. For such calculations, the nonbonded options were the same as those used in the simulations, and the images were included.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

The structures at the end of the production runs are given in Fig. 2 for systems A and B. The striking change resulting from both simulations is that the helical segment is no longer perpendicular, but is significantly tilted with respect to the membrane plane. However, it should be noted that even after the 3.0-ns (system A) and 2.0-ns (system B) production runs after the long equilibration, the N-terminal helical segment remains imbedded within the bilayer.



View larger version (66K):
[in this window]
[in a new window]
 
FIGURE 2   Orientation of the helical segment of dynorphin within the bilayers at the end of production runs of 3 ns (system A, top) and 2 ns (system B, bottom). Water molecules are omitted for clarity. For color codes, see the legend of Fig. 1. Note that the helical structure of the N-terminal segment imbedded within the bilayers is no longer perpendicular but is significantly tilted with respect to the bilayer normal. The position of the helical portion has converged to the same height in the two simulations (systems A and B).

MD trajectories of potential energies and temperatures are shown for systems A and B in Fig. 3. For system A, the potential energy remains constant for ~1 ns, then it gradually decreases and is stabilized after 2 ns (Fig. 3 A). The temperature trajectory for system A shows that the average temperature remains close to ~ 331 K during the first 1-ns simulation and then rises to remain close to ~334 K after 2 ns (Fig. 3 B). For system B the transition time period is different, with the potential energy decreasing after 0.6 ns of production run and leveling off at 1 ns (Fig. 3 A). This is mimicked by the behavior of the average temperature of system B (Fig. 3 B).



View larger version (48K):
[in this window]
[in a new window]
 
FIGURE 3   Molecular dynamics trajectories of (A) potential energy and (B) temperature for systems A (black) and B (gray).

Dynorphin helix structure and orientation

The change in orientation of the dynorphin helical segment observed at the end of the production runs (Fig. 2) prompted the analysis of dynorphin helix orientation as a function of time. The average values of backbone dihedral angles phi  and psi  for residues 1-10 in systems A and B are shown in Fig. 4. Residues 4-9 are shown to be in alpha -helical conformation throughout the simulations in both systems A and B, while residues 1-3 are alpha -helical in system B but not in system A. The orientation of the helix in the trajectories is represented by the angle between the helix axis and the bilayer normal that is along the z axis (Fig. 5 A). The vertical position of the peptide is described by the trajectories of the average Z-coordinate of the helical segments (Fig. 5 B).



View larger version (15K):
[in this window]
[in a new window]
 
FIGURE 4   Average backbone dihedral angles phi  () and psi  (black-triangle) calculated for the entire production runs for the dynorphin helical region, shown for system A (black) and system B (gray).



View larger version (34K):
[in this window]
[in a new window]
 
FIGURE 5   MD trajectories of (A) the angle between the dynorphin helix axis and the bilayer normal and (B) the average Z coordinate of the helical segment. The data for systems A and B are shown in black and gray, respectively.

For system A, the helix is parallel to the bilayer normal for ~1 ns, and then a gradual change in orientation is observed from parallel to tilted. After ~2 ns of simulation, the helix had adopted a tilted orientation that made an angle of ~50° with the bilayer normal. In system B, the helix axis changed its parallel orientation and assumed an angle of 30° with respect to the bilayer normal during the 1.4-ns equilibration. During the subsequent run, the average helix orientation changed from ~30° to ~50° and was then very similar in orientation to that observed in system A. The average Z-coordinate of the helical segment indicates that the helix of system A has moved ~2 Å toward the center of the bilayer (Z-coordinate change from ~14 Å to ~12 Å). No such movement was observed for system B, and the average Z remained around 12 Å, similar to the final position of the A system.

The above analysis shows that the A and B systems have converged at the end of the simulations, although they differed with respect to positioning of the helices at the beginning of the simulations. Noting that the final reorientation of the helix was achieved after a total of more than 3 ns of simulation in system A and 2 ns of simulation in system B emphasizes the need for long simulations to observe conformational changes in such complex systems.

The results have revealed the correlation between dynorphin helix orientation, potential energy, and temperature. To determine how the change in helix orientation depends on the interaction energy, we analyzed the contributions from interactions between peptide and lipid, peptide and water, and lipid and water. The specific role of individual residues of dynorphin in determining the interaction geometries was analyzed from their energies of interaction with lipids and water. To track the evolution of these interactions, the time intervals of the MD trajectories are divided into three regions, reflecting the transition from the parallel to the fully tilted orientations of dynorphin helix segment. Regions I and II are the time periods before and after the transition, respectively (for system A, region I: 0-1 ns; region II: 2-3 ns; for system B, region I: 0-0.6 ns; region II: 1-2 ns). The third region represents the period in which the transition occurs. Most of the analysis is reported for regions I and II.

Interaction energy analysis

The average and the rms fluctuations of the interaction energy values are given in Table 1 for regions I and II for both calculated systems. It is clear that in both systems dynorphin interacts more favorably with both the lipids and water in region II. Although the self-energies of the peptide are more favorable in time region I, this is more than compensated for by the interaction energies in region II. Thus the intramolecular interactions are broken in the tilted orientation to form more favorable intermolecular interactions. This is evident from the comparisons of peptide-DMPC and peptide-water interaction energies in regions I and II.


                              
View this table:
[in this window]
[in a new window]
 
TABLE 1   Interaction energy analysis

The interaction energies between individual residues of dynorphin with lipids and water were calculated for both systems A and B (data not shown). The arginines and lysines as well as Asp15 undergo significant interactions with both the lipids and water. While the N-terminal hydrophobic residues interact favorably with the lipids, the interactions with water are more significant for the C-terminal residues. While this pattern of interactions is anticipated from the geometry, quantification reveals that the major interactions of the peptide with both lipids and water are due to the basic residues, arginines and lysines. Thus the three arginines (Arg6, Arg7, and Arg9) and the two lysines (Lys11 and Lys13) contribute ~85% and 75% of peptide-lipid interactions, respectively, in systems A and B, and they contribute ~60% of the peptide-water interactions in both systems. Arginines are present at the carbonyl end of the N-terminal helix, while lysines connect the N-terminal helix and C-terminal turn region. The interaction energy of these residues with the headgroups of DMPC lipids (Table 2) receives positive energy contributions from choline groups, offset by the phosphate groups and to some extent by the carbonyl groups in the headgroup environment. The strong interactions of the basic residues with the phosphate group may be essential in maintaining the helical segment of dynorphin fully inside the bilayer throughout the simulations.


                              
View this table:
[in this window]
[in a new window]
 
TABLE 2   Average interaction energies (in kcal/mol) between charged residues and different components of lipid headgroups, calculated for region II

Density profiles

Average density profiles of various components of the dynorphin-DMPC system along the bilayer normal (z axis) were calculated from the last 1-ns simulations of the production run (region II); these are shown in Fig. 6 A for system B. The properties described in density profile for system B also describe system A. The density profile is in good agreement with experimental results (Lewis and Engelman, 1983; Buldt et al., 1979b; Wiener and White, 1992) and other bilayer simulation studies (Woolf and Roux, 1996; Damodaran and Merz, 1996; Shen et al., 1997; Chiu et al., 1999a; Ha Duong et al., 1999). The distribution of peptide extends from the bulk water to the center of the bilayer in the positive z axis (top layer). The higher protein density observed in the headgroup region (Z = 10-20 Å) is due to the large polar residues in the vicinity of the helix-random coil interface. This is substantiated by the density profiles of individual arginine, lysine, and aspartic acid residues that contribute to the bulk of peptide-lipid and peptide-water interactions (Fig. 6 B). The distribution of these residues falls in the same region as the phosphate and nitrogen profiles and relates to the strong interaction energies of these residues with the DMPC headgroup region identified above.



View larger version (25K):
[in this window]
[in a new window]
 
FIGURE 6   Average density profiles of (A) various components of the dynorphin-DMPC system and (B) selected residues of dynorphin in system B. The profiles were computed over structures saved at 5-ps intervals in region II. Similar density profiles were observed for system A.

An important feature of the water density profile is its asymmetry with respect to the center of the bilayer. In the side containing dynorphin, the water density is shifted toward the center of the bilayer. Previous simulation studies have observed that the presence of a peptide within the bilayer facilitates water penetration into the bilayer (Ha Duong et al., 1999; Bachar and Becker, 1999). This effect is even more evident here because the peptide is present in only one side of the bilayer. In DMPC-gramicidin simulations, a small plateau was observed near 15 Å, indicating the interactions of water with the headgroups (Woolf and Roux, 1996). In our studies, this plateau is slightly broader because of the additional interactions of waters with the basic residues present in this region as described below.

Solvation of specific side chains

The average density profiles of the dynorphin-DMPC system indicated the regions of overlap between different components of hydrated lipids with the peptide (Fig. 6). The microenvironment of the polar and nonpolar sites in the peptide is better understood from the solvation pattern of individual residues. The average numbers of water molecules, acyl chain carbon atoms, and lipid headgroups surrounding each side chain were calculated for the last 1-ns simulation (region II); these are shown for both systems in Fig. 7. It is clear that N-terminal hydrophobic residues are exposed predominantly to the lipid hydrocarbon in both systems. On the other hand, the C-terminal residues are in contact with water, choline, and the phosphate groups. Notably the basic residues that make important contributions to the interaction energy (see above) are exposed to water as well as to all components of the phospholipids. The nonpolar part of the long arginine side chains is surrounded by lipid hydrocarbon, and the positively charged guanidinium group is exposed to water and lipid headgroups (see especially Arg6 and Arg7). This type of interaction is suggested in the "snorkel model" proposed for lysine side chains (Segrest et al., 1990), in which the long and flexible side chains, extending from the hydrocarbon core to the membrane surface, are favorably located near the membrane-water interface.



View larger version (26K):
[in this window]
[in a new window]
 
FIGURE 7   Solvation of dynorphin side chains observed for (A) system A and (B) system B. Contributions from the main components of the membrane system (water oxygens, acyl chain carbons, and phosphate and choline groups) are identified by counting the number of nearest neighbors within a distance of 4 Å around each side chain. The number of contact neighbors was normalized for each side chain with respect to its total number of atoms. Solvation numbers are averaged for all structures saved at 5-ps intervals in region II.

Water penetration

Both density profiles and the details of side-chain interactions indicate that water penetrates deeply inside the top membrane layer where dynorphin is inserted. This is more evident from the average number of water molecules calculated along the z axis (Fig. 8). That the presence of peptide within the bilayers facilitates water penetration has been described from the simulations of melittin (Bachar and Becker, 1999) and a transmembrane segment of the 5HT2A receptor (Ha Duong et al., 1999). Experiments also demonstrate that water molecules more readily penetrate the hydrocarbon core of the membranes in the presence of membrane-bound proteins (Ho and Stubbs, 1992) or small peptides (Jacobs and White, 1989).



View larger version (15K):
[in this window]
[in a new window]
 
FIGURE 8   Average number of water molecules calculated in region II along the bilayer normal for systems A and B. Note that water penetrates deeply inside the top membrane layer, where dynorphin is inserted.

The distribution of water molecules within 5 Å of the peptide is illustrated in Fig. 9. Although deep penetration of water molecules is observed near the N terminus, it is interesting to note that the penetration is not uniform around the helical surface. More waters are present near polar and hydrophilic residues, and few, if any, waters are observed around one-third of the helical surface formed by the hydrophobic residues. Residues 1-5 are in the hydrophobic side of the primary amphipathic dynorphin sequence. Residue 4, 5, and 8 form the hydrophobic face of the amphipathic helix (residues 4-9). Hydration of Tyr1 seems to play an important role in positioning the peptide in the lipid bilayer. Water access is facilitated near the first three residues by the interaction with the side chain of Tyr or with the backbone at the positions of the subsequent two glycines. The backbone interaction is due either to the absence of a large side chain or to the loss of the helical conformation. The tyrosine side chain can reorient toward the aqueous environment and add to the stabilization of the waters surrounding the backbone of these glycines because of the flexibility in this region. Consequently, the hydrophilic tyrosine side chain is facing the aqueous side of the helix.



View larger version (49K):
[in this window]
[in a new window]
 
FIGURE 9   Water molecules within 5 Å of the peptide's helical region are plotted along with dynorphin at the end of production runs for both systems A (left) and B (right). Two different projections of the structures are shown, perpendicular (top) and parallel (bottom) to the membrane plane. The arginine residues are shown in white. More waters are present near polar and hydrophilic residues. There is less water penetration around one-third of the helical surface formed by hydrophobic residues.

Lipid dynamics

The effect of dynorphin on the lipid bilayer was evaluated from comparisons of various parameters calculated from the trajectories for all lipids and separately for the contact lipids that were identified as those having at least one atom within 5 Å of an element of the dynorphin helix region. Order parameter profiles were calculated for the lipid acyl chains for all lipids and for contact lipids present from the top layer. The results for system B (Fig. 10) are found to be similar for system A (not shown). The main finding is that upper-layer contact lipids are significantly more flexible than the lipids in the bottom layer, which does not contain the peptide. The flexibility seems to be due both to the fact that dynorphin has only partially traversed the lipid membrane and to the tilted conformation, which has made some more space near the center of the bilayer. Supporting these findings are the observations from fluorescence intensity measurements of leakage induced by dynorphin A(1-17) in large unilamellar vesicles composed of pure phosphatidylcholine or phosphatidylserine (Alford et al., 1996). Similarly, a dynorphin analog, E2078, has been found experimentally to increase membrane fluidity in acidic membrane lipids (Asai and Watanabe, 1999).



View larger version (16K):
[in this window]
[in a new window]
 
FIGURE 10   Order parameter profiles for all lipids and contact lipids calculated for system B. Contact lipids are those that have at least one atom within 5 Å of the dynorphin helix region. Order parameters are averaged for region II. Order parameter profiles for system A (not shown) were found to be similar.

Other average lipid parameters calculated for the A and B systems are summarized in Table 3. Whenever differences are observed, values are tabulated separately for contact lipids. Experimentally observable parameters, including the P-N vector orientation, bilayer thickness, and hydrophobic thickness, are similar to those measured in experiments and calculated in previous simulation studies (Buldt et al., 1979a; Akutsu and Nagamori, 1991; Shen et al., 1997; Ha Duong et al., 1999; Chiu et al., 1999b). The increased gauche versus trans ratio found for the contact lipids can be correlated with the increased flexibility observed in the order parameter profiles. Notably, the average Z coordinates of N and P atoms for the contact lipids in the bottom layer show that these lipids have moved ~1-2 Å closer to the center of the bilayer. This rearrangement may occur to fill the space created by the partially traversing, tilted dynorphin helix. Fig. 11 shows dynorphin along with the contact lipids at the end of the production run for system B. 


                              
View this table:
[in this window]
[in a new window]
 
TABLE 3   Average lipid properties calculated for the last nanosecond (region II) of the production run



View larger version (106K):
[in this window]
[in a new window]
 
FIGURE 11   Dynorphin in lipid bilayers at the end of the production run from the simulations of system B, with contact lipids highlighted. Note that contact lipids from the bottom layer have moved ~1 Å toward the center of the bilayer.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

In the present study the MD production runs for both systems A and B were preceded by long equilibrations of 1.6 ns (system A) and 1.4 ns (system B). The orientation of dynorphin remained nearly perpendicular to the membrane plane at the end of equilibration (Fig. 1). In contrast, in the recent simulations of the 11-residue neuropeptide substance P and its analog in sodium dodecyl sulfate micelles (Wymore and Wong, 1999b) and in a TIP3P water/CCl4 biphasic solvent system (Wymore and Wong, 1999a), the peptides, which were originally perpendicular, rearranged during the 500-ps equilibration process to a parallel orientation with respect to the interface. On the basis of these results, the authors stressed that a longer equilibration is needed to overcome the bias caused by the choice of initial configuration. Similarly, Ha Duong et al. also found that short equilibration is insufficient because it resulted in increases in the temperature in simulations of lipid bilayers (Ha Duong et al., 1999).

Many peptide simulations in bilayers have been carried out for a few hundred picoseconds (Bachar and Becker, 1999; Woolf, 1997; Chiu et al., 1995; Damodaran and Merz, 1996; Huang and Loew, 1995), leading to various degrees of agreement with experimentally determined parameters. With increased computational power, recent bilayer simulations are being carried out for longer periods ranging from one to a few nanoseconds (Belohorcova et al., 1997; Ha Duong et al., 1999; Shrivastava and Sansom, 2000; Schweighofer and Pohorille, 2000; Chiu et al., 1999b), and stronger claims of agreement with experiment are being made. In the present simulations, a change in orientation of the dynorphin helix within the bilayer started to occur after more than one nanosecond equilibration and one nanosecond production run. This emphasizes the need for long simulation times to observe any significant conformational changes in environments as complex as the present systems. While this raises the question of what would have happened if the present simulations had been continued, say for a period of ~50 ns, the tilted dynorphin helix orientation is likely to represent a stable arrangement like the one that resulted from both the A and B systems. The characteristic primary amphipathic pattern in the dynorphin sequence and the "snorkel model"-type interactions we found for the basic residues in the interface are likely to be responsible for keeping the dynorphin helix within the bilayer. The tilted orientation seems to reflect the stabilizing effects of Tyr1 and Phe4 in the peptide. This is evident from the solvation profiles in which the solvation of Tyr1 by lipid carbonyl groups and Phe4 by acyl chains has increased from region I to region II. Thus a recent study showed that the aromatic residues Phe and Trp have different affinities for different lipid components (Braun and von Heijne, 1999). While Trp has more affinity for the lipid headgroup region, Phe prefers the hydrocarbon core of the bilayer. Because the more hydrophilic Tyr is at position 1 and the Phe at position 4, the preferences for their interactions are better satisfied when the dynorphin helix adopts a tilted orientation with respect to the bilayer normal.

Based on amphiphilic moment calculations, it had been suggested that the orientation of the dynorphin A (1-13) helix is ~15° with respect to the bilayer normal (Schwyzer, 1986). However, the present MD simulation studies, using explicit atomic models, show that the helix orientation deviates even more from the perpendicular orientation to the membrane plane. As solid-state NMR studies of oriented samples of peptides in lipid bilayers have been used successfully to find the orientation of peptides in bilayers (Cross and Opella, 1994; Kovacs et al., 2000), these predictions for dynorphin can be verified directly from such studies with selectively labeled residues.

The preferential orientation of the helical segment of dynorphin at different times may be an important element in the mechanism of interaction of dynorphin with opioid receptors. Thus, in studies of lipo-derivatives of cholecystokinin peptide, the experimental data suggested a two-dimensional migration of membrane-bound ligand to reach the receptor (Moroder et al., 1993). The contemplated lateral penetration of receptor structures by the ligands is preceded by the preadsorption of peptide hormones at the cell membrane bilayer. These inferences and our present results suggest a role for the positioning of the hydrophobic face of the dynorphin helix in the bilayer. If the lateral penetration of the receptor structure through a two-dimensional migration is to occur, then it is possible that the hydrophobic face of the dynorphin helix may first contact the more hydrophobic outer surface of the transmembrane domain of the receptor. This initial contact can then lead to more stable ligand-receptor interactions. This mechanistic speculation suggests specific modes for experimental exploration, based on the present results and the available models for opioid receptors (Filizola et al., 1999; Metzger et al., 1996; Poda and Maigret, 1998; Strahs and Weinstein, 1997). It is also noteworthy that the turn region formed by residues 14-17 does not seem to influence significantly the dynorphin helix orientation within the bilayer. This may explain the fact that dynorphin A(1-17) and dynorphin A(1-13) have similar pharmacological profiles, as their mode of interactions would be the same. However, a direct comparison requires explicit simulation of the dynorphin A(1-13) within the bilayers and an understanding of the effects of specific sequence modifications and of shorter peptide analogs.

    ACKNOWLEDGMENTS

We thank Dr. Ernest Mehler for help and discussions, and Benjamin Goldsteen for the skillful administration of the computer system.

This work was supported by National Institutes of Health grants P01 DA-11470, DA-09083, and K05 DA-00060. Computational support was provided by the Cornell Supercomputer Facility and the Advanced Scientific Computing Laboratory at the Frederick Cancer Research Facility of the National Cancer Institute (Laboratory for Mathematical Biology).

    FOOTNOTES

Received for publication 13 March 2000 and in final form 19 July 2000.

Address reprint requests to Dr. Harel Weinstein, Department of Physiology and Biophysics, Box 1218, Mount Sinai School of Medicine, One Gustave Levy Place, New York, NY 10029. Tel.: 212-241-7018; Fax: 212-860-3369; E-mail: hweinstein{at}inka.mssm.edu.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES