Molecular Dynamics Simulation Provides a Possible Structure for Substance M i k e Peptides in Aqueous Solution OLLE TELEMAN',* and CLAUS-WILHELM VON DER LlETH*

'Department of Physical Chemistry 2, Chemical Centre, PO6 124, S-221 00 Lund, Sweden; 'Deutsches Krebsforschungszentrum,Zentrale Spektroskopie, Im Neuenheimer Feld, D-6900 Heidelberg, Federal Republic of Germany

SYNOPSIS

A hypothetical conformation of the undecapeptide Substance P in aqueous solution is generated by molecular dynamics simulation for 284 ps. The conformation takes explicit solvent interactions into account as well as entropic effects to the extent that phase space is sampled in simulation. The initial conformation is taken from energy minimization studies and modified. In spite of fluctuations through 180" in some backbone dihedral angles, the peptide settles with all backbone dihedrals within i 60' from the initial values. In 130 ps, the radius of gyration decreases from 6.2 A to 5.5 A, whereas only fluctuation (i.2 A ) is observed during the last 150 ps. The root-mean-square deviation at optimal superposition for a pair of conformations from the last 150 ps is 0.6 8, based on backbone atoms. The final structure is close-knit, nearly globular, and stabilized by several longlived hydrogen bonds. The simulation conformation agrees with the scarce experimental data including a large number of structure-activity relationships. Thus, the simulation conformation is a likely candidate for one of the several conformations, the existence of which has been deduced from nuclear magnetic resonance data. Simulation results and experimental modification studies suggest that Phe 8 and Leu 10 are involved in the primary binding of SP to its receptors.

INTRODUCTION

lateral brain ventricle of rat and increased muscular tone in the respiratory tract with concomitantly inSubstance P (referred to as SP), an undecapeptide creased insufflation pressure. I t is strongly natriwith the sequence H-Arg-Pro-Lys-Pro-Gln-Gln-uretic and diuretic and stimulates salivation and Phe-Phe-Gly-Leu-Met-NH2, has aroused considexocrine pancreatic secretion as well as most smooth erable interest. Pernow' lists 662 references and muscles in the gastrointestinal tract. In view of its hundreds of papers are published every year. SP is wide spectrum of physiological actions, SP is of great present in most parts of the central nervous system pharmacological interest. Although the main thrust of all mammals and, in the periphery, in the primary of Substance P research deals with its physiology, sensory neurons and the neurons of the gastroinattempts t o characterize and classify SP receptors testinal tract. Substance P acts as a neurotranshave been made.2-5 Structure-activity relationships mitter or modulator, is involved in nociception and have been searched for by means of synthetic SP causes membrane depolarization of mesenteric ganwith amino acid substitutions, by means of fragglia, often accompanied by intense neuronal disments and of fragments with substitutions.6-'" It charges. Substance P induces miosis in intact rabbit appears that SP has to be modified extensively, typeye, vasodilation, hypertension when injected in the ically a t three or more residues, in order to either increase its selectivity or obtain active antagonists." T h e C-terminal part is the more important for most 0 1990 J o h n Wiley & Sons, Inc. SP-influenced processes, but the N-terminal part is CCC 0006-3525/90/1-20013-11 $04.00 reported to be involved in some processes, such as Biopolymers, Vol. 30, 13-23 (1990) * T o whom correspondence should be addressed. the so-called antistress effect.13The C-terminal part 13

14

TELEMAN AND VON DER LIETH

is less hydrophilic, and experimental and model studies14-16indicate that the interaction of SP with a membrane leads to a n insertion into the hydrophobic region of 6 to 9 of the C-terminal residues, which then form a n a-helix. IR, CD, and NMR data 17," seem to imply an extended conformation in DMSO, and a partial a-helix in methanol and in sodium dodecanoate sulphate ( S D S ) micelles. In aqueous solution, CD data do not show any pronounced secondary structure and Chassaing et al. infer a rapid equilibrium between different conformers characterized by different hydrogen bonding patterns.17 No crystal structure has been reported, but a synthetic gene encoding SP has been expressed in E. colilg and Nakanishi reports extensive studies of the natural gene.20Attempts have been made to find favorable conformations of SP or fragments by minimization of model potential energy functions.21,22 Here, we present a 284 ps Molecular Dynamics ( M D ) simulation of CH,CO-SP in aqueous solution, where the peptide and all water molecules have been treated explicitly. The main aims were ( 1) to assess the effects on conformation from the explicit description of the solvent and from entropy, which both are taken into account in the MD simulation as opposed to the energy minimization studies; ( 2 ) to see whether the rapid equilibrium proposed by Chassaing et al.17occurs on the simulation time scale and ( 3 ) to relate molecular properties to observed structure-activity relationships.

S P C model was used, but with intramolecular flexibility." The internal bond and bond angle vibrations were treated explicitly, and, assuming harmonicity, the potential function is given by Eq. ( 2 ) , where x, and x , , , ~are the actual bond

length and the equilibrium distance, respectively, with similar notation for angles. The force constants k , and kb and the equilibrium values were taken from the The interactions due to internal rotation were handled by means of periodic dihedral (torsion) potentials, Eq. ( 3 ) , were \ I / k is a dihedral angle and ck,,a n

interaction parameter, k runs over all dihedral angles a t each bond. T h e parameters in Eqs. (1-3) are available as supplementary material to Ahlstrom et al. (cf. ref. 23). No explicit hydrogen-bond potential was used, instead hydrogen bonds are assumed to be described adequately by Coulomb and LennardJones terms.

Simulation Technique

METHODS Interaction Potential The interaction potential is that of Ahlstrom et a1.,23 although modified so that all atoms (including nonpolar hydrogens) are treated explicitly. T h e potential is based on atomic interactions. A pair ofatoms, i and j , belonging to different molecules or separated by more than two bonds interact via a LennardJones ( L J ) and a Coulomb potential (Eq. 1) . t,] and uLJare the usual Lennard-Jones parameters and qz is the partial charge of the atom i.

Note that, as the solvent water was included explicitly in the simulations, the relative dielectric permittivity t, was 1. The parameter values were based on data in the for water the

T h e Molecular Dynamics ( M D ) program used was that of Teleman and Jonsson31 and will be described only briefly. T h e trajectory is a solution to the Newtonian equation of motion and was obtained numerically by means of a fourth order predictor-corrector Gear algorithm. The program is based on independent atoms and hence includes covalent interactions in terms of harmonic and periodic potentials. Usually, this requires a rather small time step because of rapid bond and bond-angle vibrations. In order to circumvent this problem, we have devised a double time-step algorithm,31 where different time steps are used for slowly and rapidly varying degrees of freedom. Here, covalent bonds and bond angles were integrated with a time step of 0.2 fs ( 1 fs = s ) , and all other degrees of freedom with one of 1.2 fs. The two-time step algorithm gives radial distribution functions and translational diffusion coefficients identical to those obtained from a one-step algorithm with a time step of 0.2 fs. Of the simulation total of 283.5 ps, 14.7 ps was dis-

MOLECULAR DYNAMICS SIMULATION OF SUBSTANCE P

carded as equilibration and 268.8 ps was used for analysis. Also a large part of the analyzed 268.8 ps has to be considered a s equilibration. As discussed below, times refer to the analysis part of the simulation only, so t h a t 0 ps means the start of the analysis part and not the start of the entire simulation. A total of 115 cpu hours were used on a n IBM 3090150 VF. Initial Configuration

As no crystal structure of SP has been reported, the starting conformation was generated artificially. Sets of conformations of the fragments SP ( 1-4) and SP (7-11 ) were taken from Cotrait3' and Cotrait and Hospital, 33 respectively. The two fragments were joined by addition of Gln 5 and Gln 6 in either ahelix or P-sheet conformations. The potential energy surface around each in the total set of SP conformations was searched for a local minimum and the structure with the lowest minimum was chosen. Computer programs PEPCRE, 34 SIMPEP, 35 and AMBER3' were used in the process. T h e peptide was simulated in the form CH3COS P , which we denote AcSP. Because removal of up to five N-terminal residues preserves nearly all SPrelated effects, modification of the N-terminal group is expected to be uncritical. T o our knowledge, AcSP has not been synthesized and tested, but Folkers et a1.I' report an antagonist D-Trp 7,9-Leu"-AcSP (Table I1 of ref., compound 3 ) . This caused a 30fold increase in SP dose t o elicit 50% of maximal response, a t a n antagonist concentration of M, and demonstrates that blocking the N-terminus does not preclude strong binding t o receptors. AcSP and 633 water molecules were simulated under full periodic boundary conditions in a cubic cell with a side length of 27.9 A. This amounts to a concentration of 76 m M that is far above the selfaggregation concentration in water of SP,17,37 but any such eff'ects would not occur on the simulation time scale. Simulation parameters and averages are given in 'Fable 1.

RESULTS General Structure

The deformation of the peptide was monitored by means of the radius of gyration ( r c ) , moments of inertia, and root mean-square deviations a t optimal superposition. The radius of gyration decreases from 6.21 A for the initial conformation to 5.55 rt 0.004

Table I

15

Simulation Parameters and Averages

Total number of atoms Box size Time steps Equilibration time Analyzed part of trajectory Velocity scaling interval Average scaling factor Truncation distance in potential Neighbor-list update interval Average temperature Average pressure

2099 27.9*27.9*27.9 8, 0.2 fs and 1.2 fs 14.7 ps 268.8 ps 0.096 ps 0.994 10.0 A 4.8 fs 311.7 K -29. 12. MPa

*

A during the last 150 ps. The time evolution of the principal moments of inertia is uneventful. T h e peptide conformation is not elongated, (Fig. 1)and the largest moment of inertia is 50% larger than the smallest. The general structure of the peptide a t different times was compared by means of the root-mean square deviation a t optimal rigid body superposition of two conformations:

where m, is the mass of atom i and r,,a and r,,b are its position in the two conformations. The sum in Eq. (4)ran over all backbone carbon and nitrogen atoms of residues 1-11. The rms deviation was calculated for conformations from the trajectory with respect t o the initial conformation, these deviations will be referred to as Rl,start,and to the final ( R l , e n d ) (Fig. 2 ) . T h e backbone is a t almost all times during the analysis period closer to the final conformation than the initial. Between 100 and 260 ps Rl,endand Rl,sta,tfluctuate rather than change. The fluctuation time scale is 10 ps, which can be seen also from the duration of the decay in R,,endto 0 a t 268.8 ps. The rms deviation with respect to the conformation a t the beginning of the analysis period (R,,(), not shown in Fig. 2 ) is smaller than from 0 to 50 ps. From 100 ps and forthwith it is larger than and Rt,stsrt,so that, in fact, the final conformation is more similar to the initial one than to that a t 0 ps.

-

Backbone

The backbone dihedral angles were calculated to provide Ramachandran diagrams. Average 4 and $ angles are given in Table I1 together with those of

16

TELEMAN AND VON DER LIETH

Figure 1 Stereogram showing SP a t 230.4 ps. The asterisk is the center of mass. The phenylalanine a t the top is Phe 8.

the starting conformation. The full set of coordinates is available from the authors. All peptide bonds remain close to the trans position. T h e backbone conformation a t the end of the simulation can be described approximately as P’ap‘a/3’’~P’/3’ya’a’ following the notation of C ~ t r a i twhere , ~ ~ a a n d a’are

1.51

residues with a-helix dihedral angles, L means a lefthanded a-helix, 0’and 0” are residues with /3sheet dihedrals and y is a reverse y-turn. The changes in $I and J, angles are within 60”, in most residues much less, with the exception of CP in Arg 1 and of J, in Met 11, which, being end dihedrals, enjoy more flexibility. At the start of the equilibration, the entire peptide group between Phe 8 and Gly 9 rotates by 180” ( t h e Phe 8 4 and Gly 9 q5 angles change concertedly, whereas the peptide dihedral is always trans) and then back again after 18 ps. An = 60” rotation of the peptide group between Phe 7 and Phe 8 occurred after 90 ps. The $ angle of the end residue, Met 11, fluctuates extensively.

Hydrogen Bonds

I

1

I

0

100

2 00

I

I

300

Figure 2 Root-mean square deviation in A for optimal superposition of pairs of conformations. The deviation is plotted as a function of the time a t which the trajectory was sampled for one of the conformations in the pair, the other was either the initial conformation (dashed line ) or the final (full line). Time unit is picoseconds.

T h e existence of hydrogen bonds was investigated for all atom pairs where one atom ( t h e acceptor) is a n oxygen, a nitrogen, or the sulphur of Met 11 and the other a hydrogen bonded to a n oxygen or a nitrogen ( t h e donor). T h e geometry of a hydrogen bond is described by the hydrogen-acceptor distance and the donor-hydrogen-acceptor angle. We used an asymmetric criterion such that a bond is considered to be formed when the distance between hydrogen and acceptor becomes less than some distance d,

MOLECULAR DYNAMICS SIMULATION OF SUBSTANCE P

-200 I 0

1

I

20

I

I

40

I

I

60 t(ps)

Figure 3 Logarithms of the time-correlation functions for the first and second order Legendre polynomials for the vector from N,, to Ns2 in Arg 1. Characteristic times are about twice as long when fitting between 20 and 65 ps than when fitting between 5 and 15 ps. Inspection of the logarithm shows, however, that the existence of two separate reorientation processes on the time scale 5-65 ps cannot be deduced from the time correlation function, and that the variation in characteristic time reflects its statistical uncertainty.

and the angle a t the same time larger than some angle 8,. The bond is considered to cease when the distance becomes larger than d,, with d, > d,, or when the angle becomes smaller than Be, with 8, < 8,. The reason for this is that unfavorable encounters may not result in a hydrogen bond, whereas, on the other hand, a bond already formed may vibrate and occasionally be very long without actually breaking. The values for d, and d, were 3.2 and 3.6 A, and for 0, and 8,120" and 75". T h e angle dependency rules out relations where the charge dipole interaction is unfavorable. The definition of a hydrogen bclnd is not trivial, for a discussion see Ahlstrom and aLZ3 and references therein. Intrapeptide hydrogen bonds are given in Table 111. The peptide structure in the last 150 ps is stabilized by several hydrogen bonds. Most backboneto-backbone hydrogen bonds connect residues that are either neighbors or next neighbors in the amino acid sequence. A hydrogen bond exists between 0 in Phe 8 and H in Met 11, and the C-terminal end of the backbone is a mixture of a &bend and the start of a n a-helix. The bond is modified by hydrogen bonds between 0 in Phe 8 and HI and H2 in the Cterminal amide (this group will be referred to as Ami 1 2 ) and between H in Phe 8 and N in Ami 12. Apart from the bend no regular secondary structure

17

is observed. An important hydrogen bond between 0 in Leu 10 and H in Lys 3 exists almost throughout the simulation. The ring closed by this bond has 26 backbone atoms. The side-chain of Arg 1is stabilized by hydrogen bonds between H, and the four H,;s and 0 in Pro 4. The NH; of Lys 3 is hydrogenbonded very well to 0 in Gln 5 , Otl and 0 in Gln 6, and 0 in Met 11.The H, atoms of Gln 6 are hydrogen-bonded frequently to 0 in Met 11. A summation of hydrogen-bond populations shows the following. If we denote SP( 1-5) by I and SP(6-11 ) by 11, the number of 1-1, 1-11, and 11-11 hydrogen bonds is 2.0, 2.8, and 7.9, respectively, so that the C-terminal half is stabilized by six more hydrogen bonds than is the N-terminal half. This may be important in connection with the fact that SP ( 6-1 1) produces most of the effects of SP itself. The behavior of peptides to water hydrogen bonds was similar to that observed in simulations of Parvalbumin3' and Calbindin Dgk.3g Most are shortlived and a few last up to 20 ps. As expected for an undecapeptide, there are no trapped water molecules. Several carbonyl oxygens are exposed and solvated. During the last 125 ps, 0 in Arg I, 0, in Gln 5 , and 0 in Met 11, engage on the average in 1, 2, and 1 hydrogen bonds to water molecules, respectively. Also exposed polar hydrogens are solvated, H in Met 11 is, on the average, part of one hydrogen bond. Dynamics

The overall reorientation of the peptide was monitored by time-correlation functions of molecule-fixed Table I1 Dihedral Angles in Degrees Based on the Last 125 ps of the Trajectory Residue

@

\L

Arg 1 Pro 2 Lys 3 Pro 4 Gln 5 Gln 6 Phe 7 Phe 8 Gly 9 Leu 10 Met 11

-132 f 12 (0) -42 2 12 (-72) -127 2 12 (-138) -60 f 13 (-72) -88 f 14 (-124) 80 -t 16 (53) -101 f 12 (-128) -100 f 13 (-140) 87.k 9 (77) -129 f 13 (-122) -138 5 16 (-135)

83 F 15 (127) -50 f 15 (-38) 85 k 13 (117) -42 f 1 3 (8) 110 & 13 (136) 66 f 14 (85) 77 f 13 (121) 90 f 11 (106) -77 f 12 (-73) -46 f 13 (-47) 77 f 88 (-88)

Note: Numbers in parentheses are the corresponding angles for the initial configuration. The uncertainties are the root-meansquare deviation in the distribution of angles. All peptide bonds remain close to 180'.

18

TELEMAN AND VON DER LIETH

Table I11 Hydrogen Bonds Based on the Last 154 ps Atom 2

Number of Bonds

Relative Population

H Lys 3 H Gln 5 H. Arg 1 Hnzi Arg 1 Hn22Arg 1 H Gln 6 H Phe 7 H, Lys 3 Hp LYS3 H, Lys 3 HCz2Gln 5 HezzGln 5 Ha LYS3 H, LYS3 Hp LYS3 H, Ami 12 H2 Ami 12 HI Ami 12 H2 Ami 12 Ha LYS3 Hp LYS3 H, LYS3 H Phe 8 HI Ami 12 H2 Ami 12 H Gly 9 H Leu 10 H Met 11 H, Ami 12 H2 Ami 12 H Met 11 H Lys 3 H Phe 8 H, LYS3 Hp LYS3 Hp LYS3 HtZ2Gln 6 H Phe 8 HrT2Gln 6 H Phe 8

1 1 13 15 15 11 19 7 4 2 3 6 8 3 7 9 7 9 13 7 2 4 4 10 8 13 1 21 4 8 1 11 18 6 1 4 6 11 6 8

0.37 0.17 0.37 0.22 0.21 0.38 0.38 0.22 0.29 0.08 0.11 0.11 0.18 0.12 0.19 0.23 0.15 0.22 0.21 0.22 0.09 0.11 0.99 0.16 0.14 0.86 1.00 0.77 0.15 0.27 0.61 0.85 0.66 0.13 0.03 0.15 0.24 0.39 0.10 0.76

Atom 1

N Pro 2 N Pro 4 0 Pro 4

0 Gln 5

N Gln 6 O,, Gln 6

Nr2Gln 6

0 Gln 6

0 Phe 7 0 Phe 8

0 Gly 9 0 Leu 10 N Met 11 0 Gln 11

N Ami 12

Note: The two atoms forming the bond are given together with the number of bonds that last a t least 1ps and the relative population for those bonds, that is, a relative population of 1 indicates that the two atoms are hydrogen-bonded throughout the analysis period.

vectors. The time-correlation function of a vector is defined as

cj(t)= ( P , ( C O S e(s, s + t ) ) ) ,

+

(5)

where O( s, s t ) is the angle between the vector at time s and itself a time t later, P, is the j : th order

Legendre polynomial and the brackets indicate a time average. T h e j : th correlation time, r,, is defined by rou

For stochastic processes, C,( t ) decays exponentially and if several processes contribute to the reorientation, then C,( t ) will contain several exponential components. The characteristic time for such components was obtained by fitting an exponential to C;( t ) in the appropriate time domain. The definition of a molecule-fixed vector is nontrivial for flexible molecules and the following procedure was used. For each conformation, the molecule was moved so that its center of mass coincided with the origin of the laboratory coordinate frame. Three orthonormal vectors were defined arbitrarily at time 0. For each conformation the total angular displacement was used to obtain the rotation that takes the present conformation into the subsequent. In turn, these rotations were applied to the vectors, which in consequence rotate with the molecule as a whole. A detailed formulation of the procedure is given in ref. 40. Characteristic times for the overall motion was obtained from the average C j ( t ) for the three vectors, and are given in Table IV. The overall rotation is faster in the latter part of the trajectory that reflects the increasing smoothness of the peptide surface. Characteristic times for the reorientation of C,H, vectors of the individual residues are given in Table V. In a globular peptide as small as SP, the individual backbone residues are closely coupled to the overall reorientation and no residue has backbone processes more than twice as fast as the global motion. The determination of a characteristic time around 200 ps for individual residues suffers from Table IV

Overall Reorientation

Part of Peptide on which tcfrs were Based All atoms All atoms All atoms N, C,, C in res. 1-11

Part of Trajectory on which tcf:s were Based 0-269 0-134 144-269 0-269

PS

PS PS PS

71

72

(PSI

(PSI

196 198 145 201

66 65 49 69

Note: Characteristic times ( T , , T ~ were ) obtained as a fit of an exponential to appropriate parts of the first and second order time-correlation functions ( t c j x ) for vectors artificially constructed to rotate with the molecule as a whole.

19

MOLECULAR DYNAMICS SIMULATION OF SUBSTANCE P

Table V

Backbone Dynamics

t-aD

Residue Arg 1 Pro 2 Lys 3 Pro 4 Gln 5 Gln 6 Phe 7 Phe8 Gly9 Leu 10 Met 11

i2

97 155 136 112 121 161 182 161 133 173 132

32 52 49 36 40 53 60 62 43 56 42

(2-12) (2-10) (5-45) (2-25) (2-15) (2-12) (2-12) (5-40) (2-12) (2-20) (2-20)

71

72

157 223

59 72

(20-70) (12-25)

184 225 248

64 74 81

(15-50) (15-35) (15-90)

DISCUSSION Validation and Relation to Experimental Data

Any molecular modelling result requires validation. The present simulation includes a realistic environment, by means of explicit solvent molecules, as well as entropic effects that are inherent in the simulation approach. Even so, physical relevance rests on the accuracy of the peptide potential model, the quality of which can be judged only from other biomolecule simulations using similar models. Another requirement is that the trajectory samples phase-space adequately, that is, that it covers seconds, minutes, or even months (for, say, the amide proton exchange in B P T I ) . The usefulness of a 250ps simulation is based on the assumption that only the local environment of some minimum need be taken into account. This may be the case for proteins that have only one, well-defined conformation. Peptides the size of SP are sometimes quite flexible. A simulation then has to be long enough to sample all conformations if we want information on equilibrium constants or transition dynamics. Normally,

poor sampling. Thus, Phe 7 seems to undergo a reorientation with a characteristic time 50 ps larger than that of the global. Dynamics for solvent-exposed side-chains is much faster than that of buried side-chains (cf. Table VI and Fig. 3 ) . Lys 3, Gln 5, Gln 6, and Met 11 reorient with first order characteristic times between 20 and 40 ps. Phe 7, which is nearly buried, rotates as the entire peptide, whereas Phe 8, for which there is no room inside the peptide, rotates 10 times faster. The translational diffusion coefficient was obtained from the Einstein relation (Eq. 7 ) Side-Chain Dynamics

Atom Pair

Residue Arg 1 Pro 2 Lys 3 Pro 4 Gln 5 Gln 6 Phe 7 Phe 8 Leu 10 Met 11

(7)

where t is the time and r c ( t ) the position of the peptide center of mass a t time t. The translational diffusion coefficient D was found to be 4.8 f 0.6 X 10-l' m 2 / s , where the uncertainty is the standard deviation in the distribution of diffusion coefficients obtained from different parts of the trajectory.

Note: Characteristic times ( T ~ i2) , were obtained as a fit of an exponential to appropriate parts of the first and second order time-correlation functions (tcfts) for C,-H, vectors. The tcfts were based on the entire trajectory. Numbers in parentheses indicate the fitting region in the t c f x All times in ps. That different characteristic times were obtained when fitting to different time domains in the correlation function reflects the uncertainty in the characteristic times and does not indicate the presence of two separate reorientation processes.

Table VI

1 6t

D = lim-(lrc(0) -rc(t)12)

N,1

N,Z

C,

Nf CY

Ca C, C,1

Nt2

c,

c,

c,

Ne2 C,,

Cd

Ct2

CY Sa

C*l C,

71

72

78 152 38 136 32 30 206 22 74 23

27 62 17 45 13 24 67 22 30

12

(5-15) (5-40) (2-10) (5-40) (1-5) (5-35) (5-40) (5-40) (2-6) (2-10)

71

72

160

47

(20-65)

59

62

(10-90)

122

42

(10-35)

Note: Characteristic times ( T ~ T, J were obtained as a fit of a n exponential to appropriate parts of the first and second order timecorrelation functions (tcfts) for side-chain vectors as indicated. The tcf:s were based on the entire trajectory. Numbers in parentheses indicate the fitting region in the tcfts. All times in ps. That different characteristic times were obtained when fitting to different time domains in the correlation function reflects, in most cases, the uncertainty in the characteristic times and does not indicate the presence of two separate reorientation processes. The logarithm of the correlation function for the side-chain of Arg 1 is shown in Fig. 3.

20

TELEMAN AND VON DER LIETH

such a simulation is not feasible with present comfully to the solvent, which is an excellent reason for puters. However, a shorter simulation does contain low solubility in water. These residues also constiinformation based on the actually sampled part of tute a likely point of attack for the insertion of SP phase-space, that is, some domain around a regional into membranes and other amphiphilic aggregates. minimum. This means that we can model one of the conformations as we model the equilibrium conforDynamics and Comparison to Simulations mation of a protein. In fact, we can do it better beof H142 cause the peptide is smaller, quicker to simulate, The globularity of simulated SP makes for fast roand because its dynamics is faster. tational and translational diffusion. The StokesSpecific validation is difficult in the case of SP. Einstein relation provides an estimate46of r 2 : Experimental data are scarce and not very conclusive. Otter and Kotovych report more detailed rec = 0.42 ps/dalton 7 2 = cM,, sults for the SP(1-4) fragment,41but these do not seem applicable to the entire peptide. Otter and Kowhich gives a value of r2 = 0.58 ns, about 10 times tovych detected a cis-trans equilibrium for the terslower than the simulated 50-70 ps. This discrepminal Pro 4. The transition is slow on the NMR ancy is due partly to the dynamics of the water time scale and the cis form is stabilized by a salt model, as discussed by Teleman and A h l ~ t r O m . ~ ~ bridge between NH; in Lys 3 and the terminal Pro Also note that with a characteristic time of r1 = 200 4 carboxyl group, a mechanism not available to the ps, full reorientational sampling is complete within intact SP. The simulated Pro 4 is trans, and a cis0.5-1 ns and thus more or less within reach. trans rotation would seemingly affect the backbone A study of H142, a decapeptide with the sequence extensively. Recent NMR work on Calbindin Dgk Pro-His-Pro-Phe-His-Leu-Val-Ile-His-Lys that was has revealed a cis-trans equilibrium for Pro 43, and simulated in water48gave an overall rotation about this is also slow on the NMR time scale.42Chassaing 2.5 times slower than that of SP. The initial conet al.17 infer an equilibrium between several conforformation was taken from a crystal structure of the mations. In the trajectory, the SP conformation is complex of H142 with endothiapepsin, in which stable for some 150 ps, and we conclude that, as H142 has an extended conformation. It remains so expected, the transitions of this equilibrium are slow for a long time, which is the reason for the slower on the simulation time scale. Because of this flexioverall rotation. bility, at present, we do not see any way to obtain At or, = .5 5"-sN 1.12, all nuclear Overhauser experimental information on individual SP conforenhancement (NOE) signals disappear from the mations in aqueous solution. NMR ~pectrum.~' For Y = w / 2 ~= 250 MHz and The structure of SP in micelles or membranes is 500 MHz this corresponds to 7 2 = 712 and 356 ps, better characterized but is not amenable to simurespectively. Thus, the scarcity of NOE resonances lation yet. Amphiphile models have to reproduce a may be due not only to large spin separations and delicate equilibrium between entropic and electroconformational transitions, but also to motion on a static effects, and their use is not unproblematic as critical timescale. In solvents where NOE data are found for micelles in water43*44 and for bilayers in fully conclusive, that is, dimethylsulfoxide or pyriwater.45 dine, SP assumes an extended conformation. As for In the simulation, SP lacks pronounced secondary H142, this leads to substantially longer correlation structure that agrees with CD m e a s u r e m e n t ~ . ~ ~ . ~times ~,~~ and thus to increased NOE intensities. Nevertheless, SP is highly folded (again in agreement with the findings of Chassaing & al.17). They Relation to Energy Minimization Studies also found the temperature coefficient of the NMR chemical shift compatible with hydrogen bonds for Attempts have been made to map the conformaH in Gly 9, and the H, : s in Gln 5 and Gln 6. In the tional space of SP by energy m i n i m i ~ a t i o n .For ~~*~~ simulation, H in Gly 9 is hydrogen-bonded freSP, the intrapeptide potential energy is a function quently to 0 in Phe 7; and the H, : s of Gln 6 are of 594 degrees of freedom, of which bond lengths, hydrogen-bonded to 0 in Met 11when the O,] is not angles, and some dihedral angles can be estimated bonded to the Lys 3 NH;. The H, atoms in Gln 5 reasonably accurately. Multidimensional potential energy surfaces possess considerable numbers of lomainly form hydrogen bonds with water molecules. cal minima, and finding the global or even regional At neutral pH, SP aggregates in water above a few mM.17,37 Now, an undecapeptide is far too small minima is nontrivial. Systematically testing each and \k angle with 5 values leads to performing to be able to bury two phenylalanines and one leu2.5 1015minimizations for SP. In order to decrease cine. In the simulation, Phe 8 and Leu 10 are exposed

-

-

MOLECULAR DYNAMICS SIMULATION OF SUBSTANCE P

the amount of computation, it is customary to (1) describe any solvent by some simple continuum model; and ( 2 ) use a build-up procedure, beginning with dipeptides and combining favorable fragments to larger and larger pieces. Furthermore, the physical relevance of any conformation based on energy minimization rests on the assumption that entropy is negligible. The simulated structure is not close to any of the conformations reported by Manavalan and Momany22or by Nikiforovitch et al.50T h e same applies t o the fragment conformations reported by Cotrait & and C ~ t r a i t , but ~ ' here end effects may be important. In most SP conformations, backbone dihedral angles agree for some residues, in a few cases in several, but never in a systematic way, and the number of dihedrals that agree within k 15" never reaches 50%. The simulated structure is irregular and thus hard to find by build-up procedures. It may be that none of the energy-minimized conformations lie in the same domain as that traversed by the simulation, but even if one or more did, one cannot expect exact agreement as this would amount to complete cancellation of entropic effects. Agreement with the conformations of Manavalan and Momany is worse in the region where their fragments A and B were joined together, which is also the region where our two fragments were joined in generating the initial conformation. T h e more realistic environment of the simulated peptide leads to considerable fluctuations and rearrangements in structure, and several changes in backbone dihedrals are far larger than obtained by energy minimization. The importance of a detailed solvent description has been pointed out by Ahlstrom et al.23 The entropy contribution to SP binding can not be assessed until a molecular picture of the interaction of SP with its receptors has been obtained. Discussion of Mechanism

A main topic in SP research is the characterization of receptors and the interaction mechanism with these receptors. Related to this is the finding of selective agonists and antagonists that may find use in therapy or diagnosis. Quantitative structure-activity relationships (QSAR : s ) have been reported for a large number of SP homologs and a search in the literature of July-December 1987 found references to 94 such compounds. For the homologs A1a'-SP with 1 I iI 11 Couture et al. found a change in affinity only for Ala7SP (drastic decrease),Alas-SP, Ala'O-SP, and A1a"SP (moderate decrease) in the guinea pig ileum and the rat mesenteric vein.51In these preparations, the

21

N-terminal half of SP appears uncritical. In the simulation, Phe 7 is buried nearly between the inner part of the Arg 1 side-chain and the backbone of residues 8-10. One would expect replacement by the much smaller alanine side-chain to affect the backbone of residues 7-10 substantially. T h e side-chains of Phe 8 and Leu 10 are prominent and exposed, and, as discussed for aggregation, may be involved in the primary stage of receptor interaction if this has at least some hydrophobic character. Substitution of Phe 8 by tyrosine caused no reduction in affinity, whereas p -methyltyrosine and tryptophan led to moderate decrease^.^' In consequence, the binding mechanism is not a purely hydrophobic interaction. Replacement by D-Trp reduced affinity severely. D-TW can be more or less accommodated by the simulated backbone structure, but the outer half of the six-membered ring will protrude and it will constrain fluctuation in quite a different way from any L-amino acid. Reductions in affinity on substitution of Phe 8 by L-amino acids are moderate, whereas D-amino acids are devastating (ibid.) . As the simulated Phe 8 side-chain is on the outer surface, any substitution need not affect the backbone except in a minor way. D-Phe7-SP is a weak SP antagonist, but suffers from poor affinity. For a n SP homologous SP antagonist to have competitive affinity seems to require several substitutions, such as D-Arg'-D-Trp739-Leu''SP ("Spantide", see ref. 10) with pA2 = 7.1 or DPr0~-D-Trp~,~-Nle"-SP(4-11) (ref. 53) with pA2 = 7.1 on the guinea pig ileum. Dutta and al.7 report antagonists with as many as 9 substituted or modified residues. It seems that, whatever the sequence of a homolog, substitution of Leu 10 leads to loss of affinity or a t most to marginal gain, and that substitution of Phe 8 leads to loss or a t best status quo. On the other hand, substitution of Phe 7 by preferably D-Trp seems necessary to obtain a good antagonist. From these findings and the simulated structure, we suggest that Phe 7 is involved in the function of S P , whereas Phe 8 and Leu 10 are instrumental to the binding of SP to its receptor.

CONCLUSION In about 100 ps, SP settles in a conformation that is stable for the remainder of the simulation, 150 ps. This conformation has not been predicted by energy minimization studies, but agrees with the scarce experimental data on SP in aqueous solution and is a likely candidate for one of the conformations involved in the rapid equilibria put forward on the basis of NMR data.

22

TELEMAN AND VON DER LIETH

Reported structure-activity relationships are consistent with the simulated conformation that qualitatively explains several of them. The exposed position of Phe 8 and Leu 10 side-chains in simulation is consistent with the results of systematic substitution studies in the literature and with the low water solubility of SP. Thus these residues are singularly well suited to take part in the primary interaction with receptors. Several lines of further investigation suggest themselves. First, one may pursue the present simulation further with a view to observing a transition as expected from NMR data. This hope is slight within the limits set by available computers. Second, one may test the present suggested conformation by heating and recooling, a so-called annealingprocess; the peptide is simulated at a higher temperature. From an energy minimum at the higher temperature, the peptide is cooled to see whether the peptide, when back at body temperature, assumes the original conformation. This is perhaps the most promising procedure. Third, one may simulate one or more homologs to investigate differences in folding and dynamics. Last, and perhaps most interesting, one may attempt to test the a-helix formation on insertion into membranes. This study will have to await a better understanding of amphiphile models. Although Substance P is unwilling to crystallize and partake in rapid equilibria in aqueous solution, de facto preventing both X-ray crystallography and NMR spectroscopy from providing detailed molecular information, its moderate size renders it an ideal subject for statistical mechanical simulations. Again, we stress that the simulation is too short to deal with more than one of the physical conformations of SP. The conformation obtained, though, is based on the most complete theoretical treatment available sofar. Thus, statistical mechanics offers a profitable means of probing molecular properties of Substance P. We would like to thank Bo Jonsson, Sture Forsbn, Peter Ahlstrom, Anders Wallquist, and Rolf Ekman for inspired and valuable discussion, and Mats Lundell for drawing the figures. Financial support from the Swedisch Natural Sciences Research Council is acknowledged gratefully.

REFERENCES 1. Pernow, B. (1983) Phurm. Rev. 35,85-141. 2. Regoli, D., Drapeau, G., Dion, S. & D’Orleans-Juste, P. (1987) Life Sci. 40, 109-117. 3. Iversen, L. L. (1985) in Tachykinin Antagonists HBkanson, R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam, pp. 291-304.

4. Brown, J. R., Jordan, C. C., Ward, P. & Whittington, A. R. ( 1985) in Tachykinin Antagonists Hlkanson, R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam, pp. 305-312. 5. Mathison, R., Beny, J.-L., Gulati, N., Huggel, H., Mastrangelo, M., Stofer, P. & von der Weid, P.-E. (1985) in Tachykinin Antagonists HBkanson, R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam, pp. 313-322. 6. Regoli, D., Escher, E. & Mizrahi, J. ( 1984) Pharmacology 28, 301-320; Regoli, D., D’Orlbans-Juste, P., Escher, E. & Mizrahi, J. (1984) Eur. J. Pharm. 97, 161-170; ibid., 171-177; ibid., 179-189. Regoli, D., D’Orlbans-Juste, P., Drapeau, G., Dion, S. & Escher, E. ( 1985 ) in Tachykinin Antagonists HBkanson, R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam. pp. 277-287. 7. Dutta, A. S., Gormley, J. J., Graham, A. S., Briggs, I., Growcott, J. W. & Jamieson, A. (1986) J. Med. Chem. 2 9 , 1163-1171. Dutta, A. S., Gormley, J . J., Graham, A. S., Briggs, I., Growcott, J. W. & Jamieson, A. ibid., 1171-1 178. 8. Cascieri, M. A., Goldenberg, M. M. & Liang, T. (1981) Mol. Pharm. 20,457-459. 9. Cascieri, M. A., Chicchi, G. G., Freidinger, R. M., Dylian Colton, C., Perlow, D. S., Williams, B., Curtis, N. R., McKnight, A. T., Maguire, J. J., Veber, D. F. & Liang, T. (1986) Mol. Pharm. 29, 34-38. 10. Folkers, K., Rosell, S., Hdkanson, R., Chu, J.-Y., Lu, L.-A,, Leander, S., Tang, P.-F. L. & Ljungqvist, A. (1985) in Tachykinin Antagonists HBkanson, R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam, pp. 259-266. 11. Wormser, U., Laufer, R., Hart, Y., Chorev, M., Gilon, C. & Selinger, Z. (1986) E M B O J . 5, 2805-2508. 12. Hall, M. E. & Stewart, J. M. ( 1983 ) Peptides (Fayetteville) 4, 763. 13. Roske, I., Oehme, P., Hecht, K., Nieber, K., Hilse, H. & Wachtel, E. ( 1986) Pharmazie 41, 799-805. 14. Schwyzer, R., Erne, D. & Rolka, K. (1986) Helu. Chim. Acta 69,1789-1797. 15. Rolka, K., Erne, D. & Schwyzer, R. (1986) Helu. Chim. Acta 69, 1798-1806. 16. Erne, D., Rolka, K. & Schwyzer, R. (1986) Helu. Chim. A C ~69, Q 1807-1816. 17. Chassaing, G., Convert, 0. & Lavielle, S. (1986) Eur. J. Biochem. 154, 77-85. 18. Kawaki, H., Otter, A., Beierbeck, H. & Kotovych, G. ( 1986) J. Biomol. Struct. Dyn. 3, 795-803. 19. Kempe, T., Kent, S. B. H., Chow, F., Peterson, S. M., Sundquist, W. I., L’Italien, J . J., Harbrecht, D., Plunkett, D. & DeLorbe, W. J. (1985) Gene 39, 239-245. 20. Nakanishi, S. (1987) Physiol. Reu. 67,1117-1142. 21. Cotrait, M. & Hospital, M. (1986) Znt. J . Pept. Prot. Res. 28, 450-455. 22. Manavalan, P. & Momany, F. A. (1982) Int. J. Pept. Prot. Res. 20, 351-365. 23. Ahlstrom, P., Teleman, O., Jonsson, Bo, & Forsbn, S. (1987) J. A m . Chem. SOC.109, 1541-1551. 24. Hermans, J., Berendsen, H. J. C., van Gunsteren,

MOLECULAR DYNAMICS SIMULATION OF SUBSTANCE P

25. 26. 27.

28. 29. 30.

31. 32. 33. 34.

35.

36. 37.

38.

39. 40. 41.

W. F. & Postma, J. P. M. (1984) Biopolymers 23, 1513-1518. van Gunsteren, W. F. & Karplus, M. ( 1982) Macromolecules 15, 1528-1544. Margenau, M. & Kestner, N. R. ( 1969) Theory of Intermolecular Forces, Pergamon, New York. Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F. & Hermans, J. ( 1981) in Intermolecular Forces Pullman, B., Ed., D. Reidel, Dordrecht, pp. 331-342. Teleman, O., Jonsson, Bo, & Engstrom, S. ( 1987) Mol. Phys. 60,193-203. Dolphin, D. & Wick, A. E. (1977) Tabulation of Infrared Data, Wiley Interscience, New York. Herzberg, G. (1945) Molecular Spectra and Molecular Structure: Infrared and Raman Spectra of Polyatomic Molecules, van Nostrand, Princeton. Teleman, 0. & Jonsson, Bo (1986) J . Comp. Chem. 7, 58-66. Cotrait, M. (1983) Znt. J. Pept. Prot. Res. 22, 110118. Cotrait, M. & Hospital, M. (1982) Biochem. Biophys. Res. Comm. 109, 1123-1128. von der Lieth, C.-W., Palm, J., Sundin, A., Carter, R. E. & Liljefors, T. (1987) J . Mol. Graphics 5, 119125. von der Lieth, C.-W., Carter, R. E. & Liljefors, T. (1988) in Molecular Modelling Bussian, B. M. & Weber, E., EMS., Bundesministerium fur Forschung und Technologie, Bonn, pp. 103-110. Weiner, S. J., Kollman, P. A., Nguyen, D. T. & Case, D. A. (1986) J. Comput. Chem. 7,230-252. Mehlis, B., Rueger, M., Becker, M., Bienert, H., Niedrich, H. & Oehme, P. (1980) Int. J . Pept. Prot. Res. 15, 20. Ahlstrom, P., Teleman, 0. & Jonsson, Bo (1988) J. A m . Chem. SOC.110, 4198-4203. Ahlstrom, P., Teleman, O., Kordel, J., ForsBn, S. & Jonsson, B. (1989) Biochemistry 28, 3205-3211. Teleman, 0. (1990) J. Comp. Chem., 11, Otter, A. & Kotovych, G. (1987) J. Magn. Res. 74, 392-307.

23

42. Chazin, W. J., Kordel, J., Drakenberg, T., Thulin, E., Brodin, P., Grundstrom, T. & ForsBn, S . ( 1989) Proc. Natl. Acad. Sci. U S A 86,2195-2199. 43. Jonsson, Bo, Edholm, 0. & Teleman, 0. (1986) J. Chem. Phys. 85,2259-2271. 44. Watanabe, K., Ferrario, M. & Klein, M. (1988) J . Phys. Chem. 92,819-821. 45. Berendsen, H. J. C., van Gunsteren, W. F., Egberts, E. & De Vlieg, J. (1987) AC S Symp. Ser. 353 106122. 46. Cantor, R. C. & Schimmel, P. R. (1980) in Biophysical Chemistry, part 11, W. H. Freeman and Company, San Francisco, p. 461. 47. Teleman, 0. & Ahlstrom, P. (1986) J. A m . Chem. SOC.108,4333-4341. 48. Lindberg, M., Teleman, O., & Engstrom, S. (1989) Computer Modelling of the decapeptide H142, in preparation. 49. Ernst, R. R., Bodenhausen, G. & Wokaun, A. (1987) Principles of Nuclear Magnetic Resonance in One and Two Dimensions, Clarendon Press, Oxford, pp. 516522. 50. Nikiforovitch, G. V., Balodis, Y. & Chipens, G. I. (1981) Proceedings of the 16th European Symposium Brunfeldt, K., Ed., Scriptor Inc., Copenhagen, pp. 631636. 51. Couture, R., Fournier, A., Magnan, J., St-Pierre, S. & Regoli, D. (1979) Can. J . Physiol. Pharmacol. 57, 1427-1436. 52. Fournier, A., Couture, R., Magnan, J., Gendreau, M., Regoli, D. & St-Pierre, S. (1980) Can. J . Biochem. 58, 272-280. 53. Escher, E., Parent, P., Mizrahi, J., D’OrlBans-Juste, P., Dion, S., Drapeau, G. & Regoli, D. (1985) in Tachykinin Antagonists HQkanson,R. & Sundler, F., Eds., Elsevier Science Publishers, Amsterdam, pp. 267-276.

Received March 22, 1989 Accepted January 31, 1990