J. Phys. Chem. 1996, 100, 5123-5133
5123
New Theoretical Methodology for Elucidating the Solution Structure of Peptides from NMR Data. 3. Solvation Effects Hagai Meirovitch* and Eva Meirovitch* Supercomputer Computations Research Institute, Florida State UniVersity, Tallahassee, Florida 32306-4052 ReceiVed: October 9, 1995X
A short linear peptide in solution may populate several stable states (denoted here microstates) in thermodynamic equilibrium. Elucidating its dynamic 3D structure by multidimensional nuclear magnetic resonance (NMR) is complex, since the experimentally measured nuclear Overhauser effect intensities (NOEs) represent averages over the individual contributions. In previous papers (paper 1, Meirovitch, H.; et al. J. Phys. Chem. 1995, 99, 4847; and paper 2, Meirovitch, E.; Meirovitch, H. Biopolymers 1996, 38, 69) we have developed a new theoretical methodology based on statistical mechanical considerations for analyzing NMR data from flexible molecules and applied it to Leu-enkephalin (H-Tyr-Gly-Gly-Phe-Leu-OH) using the potential energy function ECEPP. The corresponding experimental NOEs have been obtained from a cryoprotective solution of this molecule. Our approach is based on a conformational search which identifies a set of significantly different low-energy structures within a certain energy range above the global energy minimun. These structures are taken as “seeds” for Monte Carlo (MC) simulations which span their vicinity; the corresponding samples are called MC microstates. Their free energies, hence populations, are obtained by the local states (LS) method, and the individual contribution of each microstate to the NOEs is calculated. The overall NOE intensity is given by the average of these individual contributions, weighted by the populations of the MC microstates. Here we apply this methodology to the same molecule described by the ECEPP energy and a solvation free energy term for water developed by Wesson and Eisenberg. This term is a summation over products of the solvent-accessible surface area of each atom and its solvation parameter. Since water is the most important solvent in biological systems, investigating the properties of this model is an imperative step in the development of our methodology. Thus, the energy barriers of the solvation model are expected to be lower than those of ECEPP alone; hence it is crucial to verify that the MC microstates of the former model are thermodynamically stable and structurally distinctive (i.e., they do not overlap). Criteria for these purposes, proposed in paper 2, are further developed here and applied to the MC microstates. We verify another crucial point, that converging results for the free energy are obtained with the LS method. In agreement with Leu-enkephalin in water being a random coil, we find that MC microstates of the solvation model have larger entropy and structural diversity than those based on ECEPP (paper 2), and they do not feature intramolecular hydrogen bonds. These physical and computational properties of the solvation model suggest that it would be applicable to other peptide systems as well. We also compare the NOEs predicted by the two models and outline plans for future work.
I. Introduction The dynamic structure of biomolecules in solution can be determined exclusively with multidimensional nuclear magnetic resonance (NMR).1 This technique is well established for globular proteins, which reside in a single microstate, i.e., a limited region of conformational space around their native structure. On the other hand, flexible systems such as peptides may populate seVeral microstates in thermodynamic equilibrium, and the quantitative interpretation of the aVerage nuclear Overhauser enhancement (NOE) effect is complex; it requires identification of the dominant microstates and determination of their populations. The existing methods treat this problem only partially; they can be classified into several categories. A commonly used approach is to acknowledge the prevalence of multiconformational equilibria and interpret the experimental data in qualitative ways (see, for example, ref 2). The 3D structures of several cyclic peptides have been determined in analogy with protein structures, in cases where the NMR data could be interpreted in terms of one predominant conformation (see, for example, ref 3). Recently, deconvolution of NOEs X
Abstract published in AdVance ACS Abstracts, February 15, 1996.
0022-3654/96/20100-5123$12.00/0
enabled the application of this approach to linear peptides.4 In some instances several interconverting structures have been invoked5 and incorporated into the data analysis.5c,d Yet, populations have not been determined with statistical mechanical methods. Peptide conformations have been studied quite extensively with restrained molecular dynamics (MD), using pseudopotentials based on NMR data, and time-averaged restraints (for example, see ref 6). A procedure called MEDUSA has been set forth by Bruschweiler and Ernst.7 The underlying hypothesis is that individual conformations may violate some of the NOE distance restraints, which are fulfilled only the entire dynamic set of substrates. Only pairs of exchanging conformations are considered, and the best combinations in terms of structural similarity are delineated.8 A few additional approaches which treat the problem of conformational multiplicity explicitly have been developed.9,10 In general, a set of reasonable conformations are generated, using any of the currently available force fields. Conformations are retained if their energies do not exceed an arbitrarily set threshold, they differ structurally, and they are consistent with a subset of the NMR restraints. The main disadvantage of these © 1996 American Chemical Society
5124 J. Phys. Chem., Vol. 100, No. 12, 1996 methods is the arbitrariness inherent in the selection of the conformations, which turn the populations into fitting parameters, rather than thermodynamic Variables. Two recent methods of Yang and Havel11 and Cicero et al.12 in several respects also pertain to this category. We have developed a new methodology for treating ensembles of interconverting conformations, which is based on pure statistical mechanical considerations and hence in principle free from the arbitrariness inherent in the “best-fit” methods; in practice one has to use approximations, which can, however, be improved systematically in an efficient way. We have defined the so-called, “localized microstate”, which is the potential energy well in conformational space around an energyminimized structure.13,14 In paper 1 (ref 14) methods (based on the harmonic entropy) were developed, which enable one to determine the relative contribution to the partition function Z of localized microstates with energy minima that pertain to a given energy range above the global energy minimum (GEM). These methods were applied to Leu-enkephalin (H-Tyr-GlyPhe-Leu-OH), for which experimental NOE data from a cryoprotective solution were obtained at 280 K,15 using the empirical potential energy function ECEPP.16 Both the force field and the molecule selected are simple, effective, and practical “tools” suitable for methodology development, which do not impair generality and serve as prototypes for other force fields and molecules. Our study suggests that the partition function can be divided into two parts, Z ) Za + Zb. Za contains the relatively small number of the lowest energy (and free energy) localized microstates and provides the dominant contribution to Z. The probability of each localized microstate is relatively large, and its contribution to medium and long-range NOE intensities is significant. On the other hand, not only is Zb smaller than Za, but the number of localized microstates associated with Zb is extremely large, which means that their individual probability is small, and their contribution to the above NOE intensities is expected to be negligible. This is so because the conformations of higher energy (i.e., those of Zb) display a much larger structural variety than those which belong to Za and can be viewed approximately as an ensemble of random coil structures, which is not expected to contribute to the long-range NOEs. We found that microstates which pertain to the 2 (3) kcal/mol range above the GEM contribute ∼0.60 (∼0.70) of the total partition function. These results provide only a global picture of the stability regions of the conformational space. In order to reproduce the experimental results one has to study in detail the most stable localized microstates, i.e., those that pertain to Za. This we did in paper 2 (ref 17), by adopting for Za an energy range of 2 kcal/mol above the GEM. However, already within this small range, the number of energy-minimized structures is large and many of them are similar. In order to avoid this redundancy and the need to identify all the microstates (an untractable task for larger molecules) we defined the so-called Monte Carlo (MC) microstate, which is a broader region of conformational space that typically includes many localized microstates. More specifically, this is the region spanned by an MC simulation which starts from a minimum energy structure (“seed”) and is characterized by a relatively large dihedral angle variability (see Figure 1). This reflects the intermediate flexibility of the peptide chain, which exceeds that of a stable state (such as the R-helix) but is still much more restricted than that of the random coil. The complete set of MC microstates is obtained from a group of “seeds” within the 2 kcal/mol range that are significantly different. This set of structures has to be determined thoughtfully; i.e., one has to verify that the MC microstates cover the
Meirovitch and Meirovitch
Figure 1. Schematic one-dimensional representation of two large potential energy wells as a function of a coordinate X. These wells are “decorated” by localized microstates, each being the basin of attraction of a corresponding local energy minimum. The localized microstates are denoted intermittently by dashed and solid lines. Each well can be spanned by Monte Carlo (MC) simulations that start from any conformation pertaining to it and not necessarily from the lowest energy seed. These wells are therefore called MC microstates. The important point is that the properties of an MC microstate depend mainly on the general shape (and depth) of the potential well and less on its detailed energy “landscape”. This suggests that MC microstates should be relatively insensitive to the type of force field used (see paper 2).
relevant part of conformational space, are distinct (i.e., do not overlap), and are thermodynamically stable (i.e., the molecule has not moved to a neighboring microstate during the simulation). In paper 2 methods and criteria for these purposes have been proposed; they are further developed in the present paper. In order to obtain populations which are thermodynamically meaningful, it is necessary to calculate the free energies of the MC microstates. This is difficult because of the intermediate flexibility of the chain, which precludes using the harmonic or quasi-harmonic approximations.18 Thermodynamic integration19 is also not useful since it is inefficient in calculating the free energy difference between microstates of a significant structural variance (therefore it has been only used in simple cases19c,d). Therefore we have employed the local states (LS) method20 which, in principle, can handle any chain flexibility. In paper 2 it was demonstrated for the first time that the LS method is, indeed, applicable to microstates with intermediate flexibility. The overall NOEs have been obtained as averages over the MC microstate NOEs, weighted by the LS populations. It has been assumed that dipolar magnetization transfer occurs exclusively between two spins, elongation and rotation of the internuclear vectors are decoupled, and angular fluctuations are ignored.5c Even though ECEPP does not take into account the cryoprotective solvent used in the experiments, we have found in paper 2 reasonable agreement between calculated and experimental NOEs. However, one has to bear in mind that the upper bounds of the experimental restraints, which need to be reproduced by the calculations, have been relaxed considerably, in view of the relatively low integrity of the experimental data. We reiterate that this small peptide and ECEPP were chosen because they constitute convenient tools for developing our new approach. Obviously, the next step requires the inclusion of solvent effects, which is the objective of the present paper. Namely, we reexamine the critical new aspects of our methodology, which have been developed for vacuum conditions in paper 2 and consolidate them for solvent conditions. Ideally solvent effects should be taken into account by using explicit models.19b,21a,b However, to date simulation of ∼20 MC (or MD) microstates of a peptide in a “box” of hundreds of water molecules is unfeasible. Also, the extension of the
Solution Structure of Peptides LS method to such models has not as yet been completed. Therefore one has to resort to more approximate mean-field continuum models, which consider the solvent implicitly. Thus, in order to obtain the electrostatic free energy of solvation the molecule is commonly characterized by a dielectric constant m, the distribution of its charges, and the geometrical boundaries that separate it from the continuum solvent; the latter is defined by s. One then solves numerically the Poisson-Boltzmann equation21 or applies approximations;22 the Langevin model of Warshel and Levitt also belongs to this category.23 The hydrophobic free energy has been taken into account as a product of the solvent-accessible surface area of the molecule and a solvation parameter obtained from data of free energy of transfer of hydrocarbons from the gas phase to water.21a,b,22,24a,b In order to obtain the total free energy, these free energies are added to the intramolecular interaction energy defined by the usual force field. Implementation of these semiempirical24c solvation energies in the framework of various simulation procedures has been carried out, but the calculations are still relatively time consuming.22a,b24a,d,e More approximate fully-empirical24c solvation models have also been developed, in which the above electrostatic and hydrophobic energies are replaced by a summation over products of the solvent-accessible surface area of the individual atoms and the corresponding atomic solvation parameter obtained from experiment;24b,25,26a,e such parameters have been derived thus far only for water. This approach is feasible due to the recent development of efficient methods for calculating the solvent-accessible surface area.26a,b,f,g In particular one should mention the extremely efficient method of Stouten et al.,25e,f which is based on volume rather than surface calculations; thus, MD simulations which consist of their solvation energy and the GROMOS force field required only 50% more computer time than simulations based on the GROMOS energy alone. Obviously, these fully-empirical solvation models are based on approximations. However, in most cases they were found to outperform the usual force fields, leading to relatively strong correlations between the (minimized) energy of a conformation and its root mean square deviation from the X-ray structure.24b,25,26a,b In view of the efficiency and the relative success of the fullyempirical models, we seek to implement them in the framework of our methodology. It would be desirable if we could study a model of Leu-enkephalin based on the ECEPP potential and a solvation energy for the cryoprotective mixture of 10% water and 90% ethylene glycol used in the experiments; unfortunately such a model is not available as yet. Therefore we have decided to study instead the solvation energy function for aqueous solution developed by Wesson and Eisenberg (see next section).25b Leu-enkephalin does not generate in water NOEs indicative of structure, implying that the molecule prevails as a random coil; therefore this model is not expected to reproduce exactly the experimental NOEs. However, since water is the most important solvent in biological systems, such a studysundertaken in this worksconstitutes an imperative step in the development of our methodology. Indeed, this solvation model is expected to demonstrate unique properties. First, its energy barriers are expected to be relatively low, which raises the question whether stable and distinct MC microstates can be defined for this model, as they had been for ECEPP alone. This is a critical aspect of our methodology, which is studied in this paper. In particular, criteria proposed in paper 2 for determining the structural distinctiveness (i.e., nonoverlap) and the thermodynamic stability of MC microstates are further developed here. It should also be pointed out that the solvation terms, which depend on the
J. Phys. Chem., Vol. 100, No. 12, 1996 5125 accessible surface area, are expected to affect the correlations among the dihedral angles; therefore it is not clear whether the LS method, which is based on these correlations, will lead [similar to ECEPP (paper 2)] to converging results for the free energies (hence, populations) of the MC microstates. The investigation of this problem is another important objective of the present paper. Finally, we seek to study various physical properties of the solvation model, such as the entropies, populations, and structural variability of the MC microstates, and the formation of hydrogen bonds and β-turns. Theoretical NOEs are also calculated and compared to those obtained in paper 2 and to the experimental ones (even though the latter are not expected to be reproduced exactly by the solvation model). This provides further information about the properties of the solvation model and the quality of the experimental NOEs. II. Theory and Methods II.1. The Model. As mentioned in the previous section, Leu-enkephalin (H-Tyr-Gly-Gly-Phe-Leu-OH) is modeled by a potential energy function that consists of ECEPP2/ and a solvation free energy term which takes into account solvation effects implicitly. ECEPP assumes rigid geometry (i.e., constant bond lengths and angles) and is based on nonbonded, electrostatic, torsional, and hydrogen-bond potentials.16 The peptide bond angles ω are fixed at 180°; therefore a conformation is defined solely by the 10 backbone dihedral angles φ and ψ and the nine side-chain dihedral angles χ. The solvation free energy Esol is
Esol )
∑ σjAj
(1)
atoms j
where σj is the solvation parameter of atom j and Aj is its conformational dependent solvent-accessible surface area. The atomic parameters σj used in this study were obtained by Wesson and Eisenberg25b from experimental values of free energy of transfer of amino acid side-chain analogue from vapor to water (Wolfenden et al.27). Before fitting these parameters they corrected the transfer free energies by taking into account the entropy of mixing, as suggested by Sharp et al.28 Table 3 of ref 25b shows the solvation parameters used. The model based on the ECEPP energy and Esol will be referred to as the solvation model. This model gives preference to conformations in which hydrophilic (hydrophobic) groups have a large (small) solventexposed surface area. In this study we have used the program FANTOM26a-e of Braun and co-workers, in which the solvation model has been implemented and an efficient procedure for calculating the accessible surface area has been provided. The solvation model assumes that ECEPP provides the correct difference of intramolecular interaction energy between two conformations in vacuum and that the corresponding difference in solvation energy is provided by Esol. In all the previous studies of continuum solvation models Esolv has been added to the ECEPP energy (or to the energy of other force fields) without further parameter optimization. However, different values for the dielectric constant have been used in order to take into account the screening of the electrostatic interaction in vacuum. For example, values of 1 (ref 25d), 4 (ref 25c,e), r (refs 24b, 25f, 26a,b), 4r (ref 22b), and 40 (ref 25b) (r is the distance between two charges) have been used for peptides and proteins, in simulations based on different force fields (GROMOS, AMBER, ECEPP, and CHARMM) and various sets of solvation energy parameters σ. In our calculations we use the experimental value, ) 50.25.
5126 J. Phys. Chem., Vol. 100, No. 12, 1996
Meirovitch and Meirovitch
II.2. The Local States Method. The local states (LS) method has been discussed in detail in previous publications,20 therefore we shall only describe here its main features. This method enables one to calculate approximately the free energy of a macromolecule from a given sample of conformations generated by any simulation technique. For simplicity we shall describe the method as applied to an MC sample of polyglycine (i.e., no side chains) of K backbone dihedral angles, modeled by the solvation model. For convenience, the dihedral angles of the molecules are ordered sequentially from the N- to the C-terminal and are denoted by Rk (k ) 1, K). In the first step one calculates from the sample the variability range ∆Rk of each dihedral angle
∆Rk ) Rk(max) - Rk(min)
(2)
where Rk(max) [Rk(min)] is the maximum [minimum] value of Rk found in the sample. Then the ranges ∆Rk are divided into l equal segments, where l is the discretization parameter. We denote these segments by νk (νk ) 1, l). Thus, an angle Rk is now represented by the segment νk to which it belongs, and a conformation i is expressed by the corresponding vector of segments [ν1(i), ν2(i), ..., νK(i)]. A local state related to νk is a partial conformation which consists of νk and the b angles preceding it along the chain, i.e., the vector (νk, νk-1, ..., νk-b); b is called the correlation parameter. For a given b one can calculate from the sample the number of occurrences n(νk,l) of all the local states from which transition probabilities p(νk|νk-1, ..., νk-b) can be obtained. Now, for each member i of the sample one determines the K local states and the corresponding transition probabilities, whose product defines an approximate probability Pi(b,l) for conformation i K
Pi(b,l) ) ∏p(νk|νk-1, ..., νk-b)
(3)
k)1
Pi(b,l) allows one to define an approximate free energy ˜A functional FA, estimated by F
F hA )
1
n
n
∑Ei(t) + kBT ln Pi(t)(b,l) ) ∑Fi(t)(b,l) ) Eh - TShA
n t)1
t)1
(4)
where n is the sample size, i(t) is the tth conformation in the sample, Ei(t) is the energy of conformation i(t), and SA is the entropy. One can show that FA constitutes a lower bound for the correct free energy F [calculated with the Boltzmann probability, rather than Pi(b,l)]; i.e., it is expected to increase as the approximation improves. One also expects the fluctuation σFA of FA to decrease as the approximation improves; it is estimated by σ j FA,
σ j FA )
[
1
n
∑(Fh A - Fi(t)(b,l))2 n t)1
]
1/2
(5)
These properties of FA and σFA constitute criteria for the statistical reliability of the results. In what follows, for the sake of simplicity, we shall not keep the distinction between the average FA and its fluctuation σFA defined over the ensemble and the corresponding estimations that appear with the bar; thus we shall denote both without the bar. The fact that FA is approximate is not a serious limitation, since in general one is interested in the difference in the free energy between two states, rather than in the absolute values. It has been found20 that in
many cases such differences ∆FAi,j(b,l) between MC microstates i and j (or similar microstates defined by molecular dynamics) A ∆Fi,j (b,l) ) FiA(b,l) - FjA(b,l)
(6)
converge rapidly as the approximation improves, i.e., as b and l increase. Therefore the converged value is expected to remain unchanged for even larger values of these parameters and thus represent faithfully the correct difference, within the statistical error. Finally, notice that the method can also handle peptides with side chains, which requires the definition of two correlation parameters b(back) and b(side), for the backbone and the side chains, respectively.20d,e The population pi of MC microstate i (we use the same notation i for a conformation and a microstate) is calculated from ∆FAi,1, the conVerged value of the difference ∆FAi,1(b,l) (eq 6) between a reference MC microstate (denoted 1) and MC microstate i A A /kBT]/∑ exp[-∆Fi,1 /kBT] pi ) exp[-∆Fi,1
(7)
i
This equation is also used for calculating populations based on the energy, by replacing ∆FAi,1 with the corresponding energy difference for each i. III. Results and Discussion III.1. Conformational Search. Since the present study is based on a potential function which includes the ECEPP terms and solvation energy terms, one should calculate the contribution of localized microstates to the partition function for different energy ranges above the GEM as were carried out for ECEPP in paper 1. However, preliminary calculations have indicated that with the solvation model the process of energy minimization is too time consuming, turning the study unfeasible. The CPU time can be reduced considerable by implementing within FANTOM analytical second derivatives of the energy. This will make possible the use of the Newton Raphson minimizer, instead of the currently used conjugate gradients minimizer. We plan to pursue this in future work; meanwhile, we decided to adopt for the solvation model as well the 2 kcal/mol range above the GEM as the energy range which defines Za (see section I). Our objective is to generate within this range a complete set of energy-minimized conformations which are significantly different (see below). For that we have used the Monte Carlo minimization (MCM) method of Li and Scheraga,29a which is also implemented in FANTOM.26-e With this method a trial conformation (obtained by selecting a random value for one of the dihedral angles of the present conformation) is energy minimized and accepted or rejected with the usual Metropolis criterion.30 While MCM does not lead to a Boltzmann sampling, it has been found to be more efficient than other methods for generating low-energy conformations. FANTOM enables one to carry out two types of MCM simulations, based on the “adaptive” and “shock” protocols, respectively. It should be pointed out that Scheraga’s group was the first to apply MCM to peptides described by ECEPP and an empirical solvation model.29b,c We carried out approximately 150 MCM calculations of 1000 MCM steps each, and several calculations of 2500 steps each, with a roughly equal number of “adaptive” and “shock” schedules. Most of these calculations started from random structures; some started from the lowest energy structures identified in paper 2 using ECEPP. The conformation with the
Solution Structure of Peptides
J. Phys. Chem., Vol. 100, No. 12, 1996 5127
TABLE 1: Backbone Dihedral Angles O and ψ (in deg) of the Five Significantly Different Energy-Minimized Backbone Structures (Motifs) Found within the 2 kcal/mol Range above (and including) -50.25 kcal/mol, Which Is Considered To Be the Global Energy Minimin (GEM) of Leu-enkephalina ψ1
φ1
ψ2
φ2
ψ3
φ3
ψ4
φ4
-83 148 -170 43 74 35 -162 -45 164 -70 -178 60 -85 163 -83 -39 -176 51 -160 149 -177 50 94 25 -99 172 -146 87 152 -72
φ5
ψ5
-94 -19 -85 142 -158 -52 -83 144 -100 -43 -90 137 -158 156 -154 -58 -158 154 -156 117
a
162 side chain conformations are defined for each of these backbone motifs and the energy of the molecule is minimized (for details see the text). The side-chain conformations that lead to minimized energy within the above range appear in Table 2.
TABLE 2: Side-Chain Dihedral Angles (in deg) for the 16 Energy Minimized Structures Found within the 2 kcal/mol Range above the GEMa no.
E
1 -50.52 2 -49.60 3 -48.28 4 5 6 7 8 9 10 11 12
χ11
χ21
χ61
χ14
χ24
χ15
174 -111 178 -53 -48 -54 -65 113 179 -176 -109 -57 -65 -64 -179 179 -93 -177
-50.13 176 -48.75 -177 -48.28 -179 -48.19 -179
76 178 86 177 -93 -179 76 179
-49.68 -49.38 -49.21 -48.92 -48.32
-84 -86 -84 111 -90
-64 -70 67 -65 58
174 -131 -56 168 -104 -172 -60 -70 -55 -66 -81 -52
-179 -178 -179 170 -179 -56 179 -65 -178 -179
13 -49.51 -177 14 -49.51 -178
-92 178 88 -178
15 -48.31 16 -48.27
-81 -178 57 -82 -178 -179
74 71
56 55
χ25 174 170 69 177 72 172 180
χ36 62 58 54
χ45 69 -53 59
-55 -170 -66 -60 -60 69 -54 -50
59 -54 174 61 -171 78 -175 70 54 -60 106 -51 -180 -55 -170 97 -82 61 -71 -180 60 -55 176 63 -171 -97 -177 -97 -177
69 174 68 174
60 60
-89 -98
64 64
57 57
178 178
53 53
a The angles are divided into five groups corresponding to the five backbone motifs presented in Table 1. The backbone dihedral angles for each structure do not differ considerably from the values of Table 1. Excluding φ1 and ψ5, the maximum deviation found is 40°, except for three dihedral angles, where the deviations are 66, 70, and 81°. E is the energy in kcal/mol.
lowest energy of -50.52 kcal/mol, obtained once, was considered to be the GEM. Four additional conformations with significantly different backbone motifs were identified within 2 kcal/mol above GEM; the backbone dihedral angles of these five structures are presented in Table 1. As in paper 2, the criterion for variance of two structures is that at least one angle differs by 60° or more, excluding the intrinsically unrestricted end-chain dihedral angles φ1 and ψ5, the angle χ6 which describes the orientation of the OH group of the Tyr side chain, and the symmetrical dihedral angles χ2 of Phe and χ3 and χ4 of Leu (for further details see related discussion in paper 2). These five structures were submitted to a systematic sidechain grid search based on their ECEPP torsional potentials (accounting for the 2-fold symmetry of χ2(Phe) and the 3-fold symmetries of χ3 and χ4 of Leu). Thus, for each backbone conformation, 162 different side-chain structures were generated as described in paper 2. These structures were subjected to energy minimization with respect to all 19 variables. Within 2 kcal/mol above the GEM, we obtained for the five backbone motifs three, four, four, two, and two different side-chains patterns, respectively. The side-chain dihedral angles of these structures and the corresponding energies are presented in Table 2. Note that for structures 13 and 14 all the corresponding angles, except for χ2 of Tyr, are equal.
TABLE 3: Results in Degrees for ∆rk (Eq 2), the Variability Range of the Dihedral Angles, and for Their Minimum and Maximum Values rk(min) and rk(max), Respectively, Obtained from the MC Sample of the Global Energy Minimum of Leu-enkaphalin (Structure 1 of Table 2) residue
angle
∆Rk
Rk(min)
Rk(max)
Tyr1
φ ψ χ1 χ2 χ6 φ ψ φ ψ φ ψ χ1 χ2 φ ψ χ1 χ2 χ3 χ4
360 146 195 118 360 243 360 185 236 141 360 183 149 125 360 177 162 360 360
-180 -312 -231 -155 -180 -290 -180 -321 -123 -179 -180 -207 -157 -171 -180 -200 -321 -180 -180
180 -166 -35 -37 180 -48 180 -136 113 -39 180 -24 -8 -46 180 -24 -160 180 180
Gly2 Gly3 Phe4
Leu5
III.2. Monte Carlo Microstates. Metropolis Monte Carlo simulations30 have been initiated from the structures of Tables 1 and 2 (the “seeds”) and carried out as described below, enabling us to define the MC microstates mentioned in the Introduction (see also Figure 1). At each MC step a trial conformation is obtained by selecting for each dihedral angle a random value within the range [-D°,D°] around the current value; thus, all the K ) 19 dihedral angles are changed at once, unlike the procedure used in paper 2, where only one angle is treated at each MC step. For Leu-enkephalin the present procedure is more efficient than the previous one by a factor of 6-10 because of the significant reduction in the number of energy calculations (see discussion of this issue in paper 2). As in paper 2, we have used D ) 1° for the glycine angles, which are relatively flexible, and D ) 3° for the remaining dihedral angles. The trial conformation is accepted or rejected, according to the usual Metropolis criterion (using the experimental temperature T ) 280 K). After repeating this procedure ng ) 40 times the current conformation is stored, and the process is continued until ns conformations have been stored (typically ns ) 2200). At this point the simulation is ended, and a new one is started from the same seed. In this way we generated four groups of conformations, each consisting of nine samples (see below) of 2200 structures, i.e., the size of the combined sample is (like in paper 2) n ) 36 × 2200 ) 79 200; this sample defines an MC microstate. Obviously, one seeks to increase ns in order to better simulate the MC microstates. On the other hand, too large a value for ns might lead to an “escape” of the molecule to a neighboring microstate that is already included in the set and should therefore not be taken into account twice. In view of the bias given to the original microstate, it is always necessary to verify that if a neighboring MC microstate is visited, it will be sampled correctly. It should be pointed out that in principle there is no limitation on the size of the MC microstate studied. However, for our purposes, a microstate which is too large has the disadvantage that it does not characterize well enough a specific structure. Thus, an approximate criterion for determining the parameters ng, ns, and D is to require that the range of change, ∆Rk (eq 2), of each angle does not exceed 120°. However, as can be seen from Table 3, most of the ∆Rk values exceed 120°, some of them considerably. Therefore, already a single mi-
5128 J. Phys. Chem., Vol. 100, No. 12, 1996
Meirovitch and Meirovitch
Figure 2. Typical “random coil” population profile obtained for ψ of Leu5 of MC microstate 1 (see Tables 1 and 2). X denotes the seed value from which the simulations are started. It can be seen that the full range ∆χ ) 360° is visited.
Figure 4. Population profile for χ2 of Leu5 of MC microstate 1′ (see text). A second peak, higher than that of the original well, is observed, which means that the molecule has “slipped” into the potential well of MC microstate 1. Microstate 1′ is discarded.
Figure 3. Typical concentrated population profile obtained for χ1 of Tyr1 of MC microstate 1. X is defined in Figure 2. Only a range of ∼40° around X is populated significantly, while the population of the “tails” is low, and their contribution to the approximate free energy FA (eq 4) is thus negligible.
Figure 5. Population profile for χ2 of Leu5 of MC microstate 1. A second small peak is observed around the seed value of this angle in MC microstate 1′.
crostate reflects a chain with intermediate flexibilitysbetween a stable structure and a random coil. We found that these ∆Rk values are larger than those obtained in paper 2, indicating that the entropies of the MC microstates obtained with the solvation model are larger than the entropies of the MC microstates obtained with ECEPP (see also later discussion in III.3). Note, however, that for most angles the ranges that contribute significantly to the free energy FA are smaller than the ranges which appear in Table 3, as discussed later in this section. It should be pointed out that ∆χ3 and ∆χ4 of Leu are 360° in Table 3. This stems from the special way these angles (which are of 3-fold symmetry) are sampled (see paper 2). The above definition of the MC microstate assumes that during our simulation procedure most angles remain in the vicinity of their seed values, while some may cover the full range of 360°. In order to check whether the molecule has moved to a neighboring microstate, one can plot (as in paper 2) for each of the 19 angles the number of occurrences n(νk,l) (see the paragraph preceding eq 3) for the l different segments comprising the range ∆Rk, as illustrated in Figures 2-5 for l ) 35. Figure 2 depicts a “random coil” type plot, consisting of two peaks, with the full range [-180°,180°] visited; it represents the dihedral angle ψ5 of structure 1 of Table 2. Figure 3 shows the plot for χ1 of Tyr of structure 1, illustrating the common case of a population profile concentrated around the seed and having asymmetric “tails”. Because of their low population,
the “tails” do not contribute significantly to the values of FA (eq 4) (however, the relatively large ranges of ∆Rk (eq 2) require increasing the value of the discretization parameter l). Thus, Figure 3 shows that the range of ∼60°, which contributes significantly to FA, is much smaller than the full range of ∆R ∼ 195°. An MC microstate characterized by the two types of profiles described above is considered to be stable. Such profiles enable one to verify that overlap between different MC microstates has not occurred. Pairs of MC microstates are considered and the corresponding profiles are compared; if at least for one angle the significant ranges of the profiles do not overlap, the MC microstates are considered to be different. A special consideration is needed in the case where the molecule moves to a neighbor MC microstate. For example, our side-chain grid search generated a seed (to be denoted structure 1′) of energy -48.90 kcal/mol pertaining to the first backbone motif, which differs from structure 1 of Tables 1 and 2 only by the value of the side-chain angle χ2 of Leu [61° versus 174° (or -299° versus -186° in Figures 4 and 5)]. However, in the course of the MC simulation of structure 1′, the values of this angle stayed only a short time around their seed, moving over rapidly to the region around the seed of structure 1, where they resided a much longer time. This is demonstrated by the profile of χ2 of structure 1′ depicted in Figure 4, where the peak around the seed of 1 is ∼2.8 times higher than the peak around the seed of structure 1′, from which the simulation started (we take the ratio of the peak heights as an approximation for the ratio of the peak areas since the peak widths are comparable).
Solution Structure of Peptides A larger population of region 1 has already been obtained for ns ) 1000, with a ratio of ∼1.4, which was found to increase to 1.9 and 2.4 for ns ) 1500 and 1800, respectively. Since our simulation is biased toward the seed of structure 1′, this “slide” means that the MC microstate is unstable; therefore structure 1′ was discarded. Note, however, that the potential energy well of 1′ is taken into account, at least partially, by the MC microstate of 1, as demonstrated in Figure 5 by the profile of χ2 of structure 1. Here a small “bump” is observed around the seed of structure 1′, indicating that this region is also populated, although significantly less than the region around the seed of structure 1 (-186°). The ratios 11.6, 11.0, and 10.5 between the double peaks (obtained for ns ) 1500, 1800, and 2200, respectively) associated with structure 1 are significantly larger than their counterparts associated with structure 1′. Also, the rates at which the peak ratios of structure 1 decrease as a function of ns are much lower than the rates at which the peak ratios of structure 1′ increase as a function of ns. This suggests that 10.5 is closer than 2.8 to the correct ratio, which is expected to be attained by increasing ns beyond 2200. Occurrence profiles with double peaks have been obtained for six out of the 16 structures of Table 3, mostly for φ4, ψ4 or φ5 (i.e., less than 5% of the variables exhibit this behavior). However, the molecule did not “slip” into one of the existing MC microstates, as has been the case with structures 1 and 1′. We checked whether these extra peaks represent additional low potential energy wells, missed perhaps by the conformational search. Namely, we subjected to energy minimization the structures obtained by the various possible combinations of these peaks; however, no new energy wells were found within the 2 kcal/mol range above the GEM. Also, the fact that the ratios between the peak heights converge with increasing ns is a strong indication that the double peaks are sampled correctly. In eight cases such a convergence has been observed. As an example, for the double peaks of structure 3, the ratios are 2.7, 2.2, 2.2, and 2.1 (for φ4) and 1.33, 1.32, 1.31, and 1.36 (for φ5), for ns ) 1000, 1500, 1800, and 2200, respectively (all the following ratio quartets will be given for these ns values). Similar behavior has been observed for structure 5 for the same angles and structure 6 for φ4 and ψ4. These MC microstates are therefore considered to be stable. On the other hand, convergence has not been obtained for φ4 of microstate 4 (the only variable with a double peak), where the ratios between the original and the extra peak decrease from 3 to 2.7, 2.3, and 2.0. The same occurs for microstate 7 for φ3 (4.7, 4.5, 3.6, and 3.1) and ψ4 (4.8, 3.8, 3.4, and 3), while the ratios of the double peak of φ4 do converge. For φ4 of microstate 9 the ratios decrease from 1.9 to 1.8, 1.7, and 1.5 (they are stable (1.3) for φ5). It should be pointed out that in each case the decrease in the peak height ratios is moderate, with the original peak being always larger than the extra peak. This should be compared with microstate 1′, where the extra peak is larger than the original peak already for ns ) 1000. As a further test we checked whether the free energy differences ∆FAi,1(b,l) (eq 6) between such an MC microstate i and a reference microstate, say 1, are stable as a function of ns; this would lead to stable populations. We calculated these free energy differences for ns ) 1000, 1500, 1800, and 2200 for MC microstates 4, 7, and 9, as well as the other microstates. In all these calculations ∆FAi,1(b,l) was typically found to fluctuate randomly as a function of ns within a 0.03 kcal/(mol‚residue) range; only in two cases did the results for ns ) 1000 deviate beyond this range. Note that 0.03 kcal/(mol‚residue) is the error in the converged values of ∆FAi,1(b,l), as discussed in the next section. It is also necessary to verify that those NOEs which
J. Phys. Chem., Vol. 100, No. 12, 1996 5129 depend on φ4, ψ4, and φ5 are not affected considerably by the unstable double peaks. To this end we calculated virtual distances for ns ) 1000, 1500, 1800, and 2200, to find that the variations do not exceed the error margin of 0.2 Å. We therefore decided to include structures 4, 7, and 9 in our set of microstates. This stability of ∆FAi,1(b,l) is not self-evident. We found that FA always decreases (due to stronger increase in TSA than in the energy) with increasing ns. Because of the relatively small change in the populations of the double peaks, the decrease in FA is approximately the same for all the structures, including 4, 7, and 9. This indicates that even though the definition of an MC microstate depends on various parameters, it is not an arbitrary entity, since within certain parameter ranges, particular stability requirements must be fulfilled. The statistical quality of an MC sample is determined by its size and the extent of correlation in the variables of interest. We calculated the autocorrelation function for the energy
CE,E(t) ) [〈E(s)E(s+t)〉 - 〈E〉2]/[〈E2〉 - 〈E〉2]
(8)
where 〈 〉 denotes an average over the MC microstate. Similar correlations were obtained for the various approximations ln Pi(b,l) (eqs 3 and 4). We found that the correlations for the energy at time t ) 1 are relatively large, between 0.71 and 0.83, decaying slowly to ∼0.37-0.61, respectively, at t ) 10. These correlations are larger than those observed in paper 2 for Leuenkephalin described by ECEPP, where the corresponding results are 0.58-0.73 for t ) 1 and 0.26-0.47 for t ) 10. Comparable correlations are observed for ln Pi(b,l) (replacing E in eq 8), where, as expected, they decrease as the approximation improves (i.e., as the parameters b and l increase). Thus, for l ) 35 and b(back) ) b(side) ) 2, the correlations are within the range 0.66-0.80 for t ) 1, decreasing to 0.35-0.49 for t ) 10, respectively. The correlations obtained with the solvation model exceed, on average, by 0.12 the corresponding correlations obtained with ECEPP. This can stem not only from the difference in the potential energy but also from the difference in the MC procedures used. In paper 2 a single angle was treated at each MC step, and the acceptance rate there was therefore relatively largesabove 0.9 for most angles and above 0.8 for all the angles. This should be compared with an acceptance rate of ∼0.65 obtained in the present study, where all the angles are changed at once. In the next section we show that in spite of these relatively large correlations, reliable results for ∆FAi,j(b,l) (eq 6) can be obtained from our samples. The MC simulations were carried out using the program FANTOM;26a-e the simulation of 79 200 conformations requires ∼80 h CPU time on an IBM SP2 RISC 6000 workstation. In practice the potential energy surface of a macromolecule is unknown. A stable (or metastable) region is commonly defined by a computer simulation sample (or a trajectory) initiated from a particular seed.18d,e,g,19 However, such a definition (which has also been adopted in this work) depends on the sample size and therefore is somewhat arbitrary (see, however, relevant discussion above) and determination of the errors is not straightforward (see detailed discussion in ref 20e). In order to reduce this arbitrariness we generate all the MC microstates with the same set of parameters (ns, D, etc.), check that microstates do not overlap, and verify (see next section) that within a certain range, the results are relatively insensitive to the choice of these parameters. III.3. Calculation of Free Energy and Populations. In order to obtain the relative contribution of the 16 MC microstates to the partition function, we calculated their free energy by the
5130 J. Phys. Chem., Vol. 100, No. 12, 1996
Meirovitch and Meirovitch
A TABLE 4: Differences in Free Energy, ∆F6,1 (b,l) (Eq 6) [in kcal(mol‚residue)] between MC Microstates 6 and 1a
0
1
2
3
b(side) l/b(back)
1
2
3
1
2
3
1
2
3
1
2
3
12 16 20 25 35 45
0.18 0.19 0.17 0.17 0.16 0.16
0.17 0.17 0.15 0.15 0.14 0.14
0.17 0.17 0.16
0.16 0.16 0.14 0.14 0.13 0.13
0.15 0.15 0.13 0.12 0.12 0.12
0.15 0.15 0.14
0.15 0.16 0.14 0.14 0.13 0.13
0.14 0.14 0.12 0.12 0.11 0.12
0.14 0.14 0.13
0.15 0.16 0.14
0.14 0.14 0.13
0.14 0.15 0.14
a b(back) and b(side) are the correlation parameters for the backbone and side chains, respectively. b(side) ) 0 means that the transition probabilities are uncorrelated. l is the discretization parameter. The results in each column show convergence as l increases. For every section of b(side), the results in each row converge as b(back) is increased. The same happens for the values of a given b(back), as a function of b(side). The converged value of ∆FA is 0.12 ( 0.03 kcal/(mol‚residue).
LS method (see also ref 20c-e). The average NOE intensities corresponding to the individual MC microstates were also calculated (see III.4). As in paper 2, for each microstate a sequence of approximations for the free energy FA (eq 4) was obtained, with backbone and side chain correlation parameters 0 e b(back) e 4 and 0 e b(side) e 3 (see II.2) and discretization parameters 6 e l e 45. Each calculation for a given value of l and all the approximations that depend on b(back) and b(side) required ∼5 min CPU on an IBM 375 RISC 6000 workstation. We verified that FA and its fluctuation σFA (eq 5) always increases and decreases, respectively, as the approximation improves, i.e., as b(back), b(side), and l increase. This is an indication that the results are statistically reliable. We find that the difference between the worst approximation FA [b(back) ) b(side) ) 0, l ) 1)] (i.e., no correlations and no discretization, see II.2) and the best approximation [b(side) ) b(back) ) 2, l ) 45] lies within the range 10.09-14.05 kcal/mol, which is much larger than the statistical error of FA. This demonstrates that our approximation is effective. The corresponding decrease in the fluctuations is of the order of 0.59-0.89 kcal/mol. In order to obtain the populations of the MC microstates, one has to calculate the difference in the approximate free energy FA between a reference microstate and all the other MC microstates. We calculated ∆FAi,1(b,l) (eq 6), choosing structure 1 of Table 2 as a reference, for various approximations (b,l), seeking convergence of the results, as the approximation improves. The converged values are considered to be the true differences within the statistical errors. Results for ∆FA6,1 are presented in Table 4; convergence to ∆FA6,1 ) 0.12 ( 0.03 kcal/(mol‚residue) is clearly demonstrated. We generated 15 similar tables for all the other MC microstates and found convergence in all of them and statistical errors that are not larger than (0.03 kcal/(mol‚residue). In order to check the dependence of our results on the sample size ns, we calculated results for FA and ∆FAi,1(b,l) for seven microstates i, based on smaller samples of ns ) 1800, 1500, and 1000. For each microstate the results for ∆FAi,1 were found to be stable, fluctuating within a range that is not larger than 0.03 kcal/(mol‚ residue). Only in two cases did the results for ns ) 1000 deviate slightly from this range. This supports our estimated error of (0.03 kcal/(mol‚residue). In Table 5 we provide the populations of the 16 MC microstates. As in Table 2, the results are grouped according to the backbone motifs of the seeds, i.e., the structures from which the MC simulations were started; within each group the results appear in an increasing order of the corresponding (minimized) seed energies Eseed i . These energies together with and the the average energies of the MC microstates Eave i converged differences of the LS free energies ∆FAi,1 (eq 6) also appear in the table. These quantities lead (using eq 7) to three LS ave sets of populations denoted pseed i , pi , and pi , respectively.
TABLE 5: Results for the Populations of the 16 Microstates, Organized in Five Groups Corresponding to the Different Backbone Motifs and Presented in an Increasing a Order of the Seed Energy Eseed i i
Eiseed
Eiave
A ∆Fi,1
piseed
piave
piLS
1 2 3
-50.52 -49.60 -48.28
-43.08 -43.20 -42.17
0.00 -0.35 -0.20
0.37 0.07 0.01
0.09 0.11 0.02
0.10 0.19 0.14
4 5 6 7
-50.13 -48.75 -48.28 -48.19
-42.39 -41.77 -41.30 -41.22
0.50 1.35 0.60 0.75
0.18 0.02 0.01 0.01
0.03 0.01 0.00 0.00
0.04 0.01 0.03 0.04
8 9 10 11 12
-49.68 -49.38 -49.21 -48.92 -48.32
-43.19 -42.80 -42.98 -42.38 -41.73
-0.35 0.05 0.75 -0.05 1.45
0.08 0.05 0.03 0.02 0.01
0.11 0.06 0.08 0.03 0.01
0.19 0.09 0.03 0.11 0.01
13 14
-49.51 -49.51
-43.60 -43.48
1.30 1.20
0.06 0.06
0.24 0.19
0.01 0.01
15 16
-48.31 -48.27
-42.17 -41.07
2.15 2.25
0.01 0.01
0.02 0.00
0.002 0.002
A a ave is the Ei is the average energy of MC microstate i. ∆Fi,1 A converged value of the free energy differences ∆Fi,l (b,l) (eq 4 and 6) with respect to i ) 1, the microstate with the lowest seed energy. All these results are in kcal/mol. The populations based on these quantities obtained with eq 7 are piseed, piave, and piLS. The statistical errors of Eiave are smaller than (0.02 and those for piave are smaller than the last A are (0.15 kcal/mol, which lead to digit presented. The errors of ∆Fi,1 LS errors of up to (0.03 in pi for all the structures except for 2, 3, and 8, where the errors are not larger than (0.05. These errors are correlated due to the normalization condition (see text for details).
Notice that only pLS i takes into account entropic effects. The statistical error of the average energy is smaller than (0.02 kcal/ mol; this induces errors in the third digital place in the corresponding populations pave i . Obviously, all these errors are correlated, due to the normalization condition of the probabilities. (0.55) is Table 5 reveals that most of the probability of pseed i contributed by seeds 1 and 4, which have the lowest energies. are different and less concentrated, with The results for pave i structures 13 and 14 providing the largest population of 0.43 is only 0.12). The results for pLS (their contribution to pseed i i seed patterns; the discrepandiffer from both the pi and the pave i cies are most significant for structures 13 and 14 (but also for 2, 3, 8, and 11), demonstrating that entropic effects are important and should not be neglected. Note that the first three backbone motifs contribute significantly (0.43, 0.12, and 0.43, respectively); with the ECEPP model (paper 2) we found that the backbone motif of the GEM contributes 0.98 of the population. This indicates that the solvation model provides a better description of the random coil nature of Leu-enkephalin in aqueous solution, a picture also supported by comparing the
Solution Structure of Peptides average entropy TSA (eq 4) of the MC microstates here and in paper 2. Thus, for both the lowest approximation [l ) 1, b(side) ) b(back) ) 0] and the better approximation [l ) 35, b(back) ) b(side) ) 2], at T ) 280 K the average entropy obtained with the solvation model is larger by ∼2 kcal/mol than that obtained with the ECCEP model (see also discussion in III.2 and the results of Table 3). In fact, this discrepancy would have been even larger if the value of ng had been increased in the solvation model simulations. This would have compensated for the lower acceptance rate and higher energy and entropy autocorrelations in the solvation model calculations, as compared with the ECEPP calculations. III.4. Results for NOEs. In this work (as in paper 2) we adopt the model for the motional state of the molecule proposed by Kessler et al.5c (see also refs 8, 31, and 32). This model consists of the “initial rate approximation”, the assumptions that intramicrostate conformational exchange is much faster, whereas intermicrostate exchange is much slower, than the overall rotational reorientation and that the various types of motions are uncorrelated (a detailed discussion of the nature of these motions appears in paper 2). The theory outlined in ref 5c is applicable provided that the rotational correlation time is longer than the inverse Larmor frequency, which is a valid hypothesis for the Leu-enkephalin system under investigation, as indicated by Motta et al.33 Note that angular modulations32b of the internuclear vectors are being ignored in this model, as discussed later on. Thus, the dipolar cross-relaxation rate of MC microstate i is in direct proportion to the average 〈1/r3〉2i , calculated over the conformations within microstate i; an effective is defined as reff ) 〈1/r3〉-1/3 . The interatomic distance reff i i i overall dipolar cross-relaxation rate Γ of the 16 MC microstates eff -6 is the weighted average Γ ) ∑pLS i (ri ) . In the “initial rate approximation” Γ is proportional to the NOE intensity I, where the proportionality constant is determined by scaling intensities with respect to reference peaks. I-6 is defined as the virtual theoretical distance, which will be compared with the corresponding virtual experimental distance. As has been pointed out above, in principle, it is necessary to account for both radial and angular modulations; the former increase, whereas the latter decrease, the NOE. However, a recent theoretical study34 has shown that the generalized angular order parameters of the internuclear vectors of the pentapeptide VPGDV are typically higher than 0.9 and the distance errors introduced by ignoring the angular part are of the order of 15% in the worst case. Moreover, the maximum error introduced in our calculations is expected to be smaller, for two reasons. First, the errors in ref 34 are calculated with respect to rates that are based on 〈r-6〉i, while our rates consist of 〈r-3〉2i , which lead to smaller errors (compare eqs 38 and 42 of ref 34). Second, the 2.18-ns MD trajectory of YPGDV in ref 34 defines a microstate that is much more flexible than the present MC microstates; this can be deduced from the relatively flat distributions of the dihedral angles of Gly3 within the whole range of [-180°,180°] Figure 2c,d of ref 34). Therefore the angular order parameters in our case are expected to be larger than those of YPGDV in ref 34. Since the experimental NOEs of Leu-enkephalin are not very accurate and precise and our model describes this molecule in water rather than in the experimental cryoprotective solvent, we did not consider it warranted to carry out the extensive extra calculations to obtain site-specific angular correlations (individual motional parameters need to be calculated for all the internuclear vectors, of all the MC microstates). We intend to take into account the angular correlations as well as other refinement procedures, such as calculations based on a complete
J. Phys. Chem., Vol. 100, No. 12, 1996 5131 TABLE 6: Sequential Virtual Distances for Cr-NH, NH-NH, and Cβ-NH, Based on Experimental Data (Table 3 of Ref 20) and Theoretical Calculationsa 1A-2N 3A-4N 4A-5N 2N-3N 1B-2N 4B-5N exp. sol. ECEPP
2.8 2.4 2.4
3.8 3.1 2.9
2.8 3.3 3.4
3.5 2.6 3.1
4.6 3.3 3.8
4.6 3.2 3.4
a The atom notations are defined in the text. sol (ECEPP) denotes NOEs calculated with the solvation (ECEPP) model. The statistical error in the theoretical distances is not larger than 0.2 Å.
TABLE 7: Sequential Virtual Distances Not Included in Table 6a exptl sol ECEPP a
3A-4A
4A-5A
3A-4B
4B-5A
4.6 4.4 4.3
3.6 4.6 4.7
5.6 5.3 5.0
4.6 4.8 4.7
Otherwise as in Table 6.
relaxation matrix35 (instead of assuming the “initial rate approximation”), in future work, when better experimental data will be used and the microstates will be simulated by MD (rather than the Metropolis MC method). We have calculated theoretical counterparts to the interresidual experimental connectivities presented in Table 3 of ref 15. Table 6 contains the sequential contacts CR-NH, NHNH, and Cβ-NH, which are likely to have substantial contributions from both folded and unfolded states.4 The remaining connectivities appear in Table 7 (sequential) and in Table 8 (long range). The theoretical data shown in Tables 6-8 are virtual distances, calculated from the 16 microstates of the solvation model as outlined previously; for comparison we also provide the corresponding results from paper 2 (based on 21 MC microstates obtained with ECEPP). The experimental data in these tables are upper-bound values of virtual distances,15 which only defines strong medium, and weak NOEs. The translation of the latter to distances and the bound relaxation adopted are described in detail in paper 2. The distances are denoted by the residue number and one or several letters, which specify either the heavy atoms to which the hydrogen atoms are bonded or the pseudo-atoms defined by the experiment. Thus, A represents the CR carbon, N the amide nitrogen, M the pseudoatom comprising the two methyl groups of Leu, B the β methylene pseudoatoms, R the benzene ring pseudoatoms of Tyr and Phe, and BG the pseudoatom of Leu containing the spectroscopically unresolved β and γ protons. This notation allowed us to designate both sequential and long-range NOEs in a similar fashion. The experimental and ECEPP results of Tables 6 and 7 have been discussed in paper 2, where the reader can get the details. It should be pointed out that the distances obtained with the solvation model are comparable to those obtained with ECEPP, where the maximal deviation between the results is 0.5 Å (1B2N). A theoretical results is considered to violate the corresponding experimental value if it exceeds the latter. The disagreement with experiment occurs for both models for 4A5N (Table 6) and 4A-5A (Table 7) (see paper 2). The longrange distances are shown in Table 8. With both models three out of the 13 distances are violated. They are 4A-5M, 1A5BG, and 1A-5M for the solvation model with violations of 1.2, 0.1, and 0.4 Å, respectively, and 4A-5M, 1R-5M, and 4A-5BG for the ECEPP model, with violations of 1.7, 0.5, and 0.4 Å, respectively. Altogether [ignoring the contact 4A5A for which the experimental result is probably an artifact (see paper 2)], five out of 23 distances are not satisfied by both models, which means that they provide a comparable fit to
5132 J. Phys. Chem., Vol. 100, No. 12, 1996
Meirovitch and Meirovitch
TABLE 8: Long-Range Virtual Distancesa exptl sol ECEPP
exptl sol ECEPP a
4R-5M
4B-5BG
4R-5N
1R-5M
4A-5BG
1R-2N
7.6 6.3 6.8
6.6 4.9 5.6
5.6 5.1 4.5
7.6 5.2 8.1
5.6 5.6 6.0
4.6 4.0 2.9
3A-4R
4R-5A
4R-5BG
4A-5M
1A-5BG
1A-5M
5M-4B
6.6 6.6 6.2
5.6 5.4 5.1
7.6 6.0 5.9
5.6 6.8 7.3
6.8 6.9 4.2
6.8 7.2 5.3
7.8 5.6 6.3
Otherwise as in Table 6.
experiment. Notice, however, that the corresponding distances of these models in Table 8 are in general different where the maximal deviation of 8.1-5.2 ) 2.9 Å occurs for 1R-5M. It is also of interest to discuss the formation of hydrogen bonds and β-turns. A potential donor and acceptor are considered to be hydrogen bonded if they are not further apart than 2.4 Å. In paper 2 (Table 9), we found five NH-CO distances within the 2.4-Å range. Conversely, with the present solvation model, not a single NH-CO distance associated with the 16 structures considered was found to be below 2.4 Å, in accord with the inherent prediction that peptide-water, rather than peptide-peptide, bonds should be formed. The different types of β-turns in proteins have been classified by Wilmot and Thornton.36 Using their definition we find with the solvation model that the backbones of all the seed structures, except for structure 2, feature a 1-4 type IV turn, similar to the motifs 2-4 identified with the ECEPP model in paper 2. On the other hand, backbone motif 1 of the ECEPP model, which has ∼90% of the population, forms a 2-5 type II′ turn. This discussion on hydrogen bonds and turns indicates that the stable structures obtained with the solvation and ECEPP models differ significantly. This is supported by the considerable discrepancy between virtual distances obtained with the two models, where 16 distances in Table 6-8 differ between 0.3 and 2.9 Å. Since both models provide a comparable fit to the experiment, it appears that the experimental results available for Leu-enkephalin are not sufficiently accurate and precise to discriminate between them. IV. Summary In the present study we have carried out, to the best of our knowledge, the first extensiVe conformational search of a peptide based entirely on the total energy, which includes both intramolecular and solvation terms. (In previous studies the energy minimization was carried out with respect to the force field energy alone, while the total energy was calculated at the minimum; e.g., see refs 22b and 24b. Only limited searches have been carried out with the entire energy; e.g., see ref 26f.) We have further developed methods which enable one to examine the thermodynamic stability of the MC microstates and their structural distinctiveness. We found that, similar to ECEPP, most of the MC microstates of the solvation model are stable and that their free energy can be obtained reliably with the local states method. This is demonstrated by the converging (for improving approximations) results obtained for differences in FA, which lead to the populations of the MC microstates. As has been pointed out in the Introduction, these findings are not trivial, and they are crucial for future applications of the solvation model. As expected, the solvation model based on the parameters of Wesson and Eisenberg provides a much better description than ECEPP of Leu-enkephalin in water, which is known experimentally to be a random coil. Thus, the solvation model
generates three significantly populated backbone motifs, whereas with ECEPP only a single backbone motif dominates. Also, the MC microstates obtained with the solvation model feature higher entropy than the ECEPP microstates. The lack of intramolecular hydrogen bonds in the seeds of the solvation model, versus four hydrogen bonds displayed by the ECEPP seeds, is also in accord with the expected random coil characteristics of the solute, as well as the existence of intermolecular peptide-water hydrogen bonds. These findings suggest that implicit solvation models may be adequate for treating structural properties of peptides and that they might be improved by studying the fit between theoretical and experimental NOEs. However, the experimental results available for Leu-enkephalin are not accurate enough for such studies. Thus, although structural features such as interatomic distances, predicted by the solvation model and ECEPP, display considerable variance, the two models provide comparable fit to the experimental upper distance bounds. Unfortunately, experimental NOEs for short linear peptides available in the literature are, in general, comparable in quality to the Leu-enkephalin data. Consultation with experimentalists has made it clear that in many cases the published NOEs have been obtained with relatively low accuracy and precision, not because of experimental limitations, but because adequate methods of analyses were not available. For this reason the extra effort required to improve the integrity of the experimental NOEs was not deemed justified. We hope that by offering a methodology for interpreting the NOE data in terms of the dynamic structure of the molecule investigated, we have provided the necessary motivation for acquiring, whenever possible, higher quality NOEs from solutions of linear peptides. Until such data become available, we will apply our methodology to more complete, accurate, and precise NMR data which have been obtained from solutions of cyclic peptides and also to flexible loops in proteins using the Stouten et al. solvation model.25d Our methodology is based on pure thermodynamic considerations, and it does not rely on experimental results. On the other hand, with the commonly used analysis, several structures that satisfy the experimental data are typically determined by distance geometry or restrained molecular dynamics.3 Since for larger peptides a complete conformational search is unfeasible, one can combine the conventional methods with our approach by applying our methodology to the subspace of conformations that are compatible with the experimental data. Acknowledgment. We are grateful to Dr. Max Va´squez for valuable discussions. We thank Dr. T. C. Oppe for his help in installing the program FANTOM. We acknowledge support from the Florida State University Supercomputer Computations Research Institute, which is partially funded by the US Department of Energy (DOE) under Contract DE-FC05-85ER250000. This work was also partially supported by DOE Grant DEFG05-95ER62070.
Solution Structure of Peptides References and Notes (1) (a) Wu¨thrich, K. NMR of Proteins and Nucleic Acids; John Wiley and Sons: New York, 1986. (b) Clore, G. M.; Gronenborn, A. M. Annu. ReV. Biophys. Chem. 1991, 20, 29-63. (c) Bax, A.; Grzesiek, S. Acc. Chem. Res. 1993, 26, 131-138. (d) Oschkinat, H.; Muller, T.; Diekmann, T. Angew. Chem., Int. Ed. Engl. 1994, 33, 277-293. (2) (a) Dyson, H. J.; Rance, M.; Houghten, R. A.; Lerner, R. A.; Wright, P. E. J. Mol. Biol. 1988, 201, 161-200. (b) Stradley, S. J.; Rizo, J.; Bruch, M. D.; Stroup, A. N.; Gierasch, L. M. Biopolymers 1990. (3) (a) Peishoff, C. E.; Ali, F. E.; Bean, J. W.; Calvo, R.; D’Ambrosio, C. A.; Eggleston, D. S.; Hwang, S. M.; Kline, T. P.; Koster, P. F.; Nichols, A.; Powers, D.; Romoff, T.; Samanen, J. M.; Stadel, J.; Vasko, J. A.; Kopple, K. D. J. Med. Chem. 1992, 35, 3962-3969. (b) Sautilis, J.; Mierke, D. F.; Byk, G.; Gilon, H.; Kessler, H. J. Am. Chem. Soc. 1992, 114, 4818-4818. (4) Yao, J.; Dyson, H. J.; Wright, P. E. J. Mol. Biol. 1994, 243, 754766. (5) (a) Kessler, H.; Matter, H.; Gemmeker, G.; Kottenhahn, M.; Bats, J. W. J. Am. Chem. Soc. 1992, 114, 4805-4818. (b) Bogusky, M. J.; Naylor, A. M.; Pitzenberger, S. M.; Nutt, R. F.; Brady, S. F.; Colton, C. D.; Sisko, J. T.; Anderson, P. S.; Veber, D. F. Int. J. Peptide Res. 1992, 39, 63-76. (c) Kessler, H.; Griesinger, C.; Lautz, J.; Muller, A.; van Gunsteren, W. F.; Berendsen, H. J. C. J. Am. Chem. Soc. 1988, 110, 33933396. (d) Mierke, D. F.; Kurz, M.; Kessler, H. J. Am. Chem. Soc. 1994, 116, 1042-1049. (6) Mierke, D. F.; Scheek, R. M.; Kessler, H. Biopolymers 1994, 34, 559-563. (7) (a) Bru¨schweiler, R.; Blackledge, M.; Ernst, R. R. J. Biomol. NMR 1991, 1, 3-11. (b) Bru¨schweiler, R. Structural Dynamics in Biomolecules Monitored by Nuclear Magnetic Resonance Relaxation. Ph.D. Thesis, ETH, Zurich, 1991. (8) Blackledge, M. J.; Bru¨schweiler, R.; Griesinger, C.; Schmidt, J. M.; Xu, P.; Ernst, R. R. Biochemistry 1993, 32, 10960-10974. (9) (a) Nikiforovich, G. V.; Prakash, O.; Gehrig, C. A.; Hruby, V. J. J. Am. Chem. Soc. 1993, 115, 3399-3406. (10) Landis, C.; Allured, V. S. J. Am. Chem. Soc. 1991, 113, 94939499. (11) Yang, J.; Havel, T. F. J. Biomol. NMR 1993, 3, 355-360. (12) Cicero, D. O.; Barbato, G.; Bazzo, R. J. Am. Chem. Soc. 1995, 117, 1027-1033. (13) Va´squez, M.; Meirovitch, E.; Meirovitch, H. J. Phys. Chem. 1994, 98, 9380-9382. (14) Meirovitch, H.; Meirovitch, E.; Lee, J. J. Phys. Chem. 1995, 99, 4847-4854. (15) Picone, D.; D’Ursi, A.; Motta, A.; Tancredi, T.; Temussi, P. A. Eur. J. Biochem. 1990, 192, 433-439. (16) (a) Momany, F. A.; McGuire, R. F.; Burgess, A. W.; Scheraga, H. A. J. Phys. Chem. 1975, 79, 2361-2381. (b) Ne´methy, G.; Pottle, M. S.; Scheraga, H. A. J. Phys. Chem. 1983, 87, 1833-1887. (c) Sippl, M. J.; Ne´methy, G.; Scheraga, H. A. J. Phys. Chem. 1984, 88, 6231-6233. (17) Meirovitch, E.; Meirovitch, H. Biopolymers 1996, 38, 69-88. (18) (a) Goj, N.; Scheraga, H. A. J. Chem. Phys. 1969, 51, 4751-4767. (b) Gibson, K. D.; Scheraga, H. A. Physiol. Chem., Phys. 1969, 1, 109126. (c) Goj, N.; Scheraga, H. A. Macromolecules 1976, 9, 535-542. (d) Karplus, M.; Kushick, J. N. Macromolecules 1981, 14, 325-332. (e) DiNola, A.; Berendsen, H. J. C.; Edholm, O. Macromolecules 1984, 17, 2044-2050. (f) Hagler, A. T.; Stern, P. S.; Sharon, R.; Becker, J. M.; Naider, F. J. Am. Chem. Soc. 1979, 101, 6842-6852. (g) Case, D. A. Current Opinion Struct. Biol. 1994, 4, 285-290. (19) (a) Beveridge, D. L.; DiCapua, F. M. Annu. ReV. Biophys., Biophys. Chem. 1989, 18, 431-498. (b) Computer Simulation of Biomolecular Systems; van Gunsteren, W. F., Weiner, P. K., Eds.; ESCOM: Leiden, 1989. (c) Tobias, D. J.; Sneddon, S. F.; Brooks, C. L., III. J. Mol. Biol. 1990, 216, 783-796. (d) Tobias, D.; Brooks, C. L., III. Biochemistry 1991, 30, 6059-6070. (20) (a) Meirovitch, H. Chem. Phys. Lett. 1977, 45, 389-392. (b) Meirovitch, H. Phys. ReV. A 1985, 32, 3709-3715. (c) Meirovitch, H.;
J. Phys. Chem., Vol. 100, No. 12, 1996 5133 Va´squez, M.; Scheraga, H. A. Biopolymers 1987, 26, 651-671. (d) Meirovitch, H.; Kitson, D. H.; Hagler, A. T. J. Am. Chem. Soc. 1992, 114, 5386-5399. (e) Meirovitch, H.; Koerber, S. C.; Rivier, J.; Hagler, A. T. Biopolymers 1994, 34, 815-839. (21) (a) Smith, P. E.; Pettitt, B. M. J. Phys. Chem. 1994, 98, 97009711. (b) Van Gunsteren, W. F.; Luque, F. J.; Timms, D.; Torda, A. E. Annu. ReV. Biophys. Biomol. Struct. 1994, 23, 847-863. (c) Sharp, K. A.; Honig, B. Annu. ReV. Biophys. Chem. 1990, 19, 301-332. (d) Davis, M. E.; McCammon, J. A. Chem. ReV. 1990 90, 509-521. (22) (a) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (b) Abagyan, R.; Totrov, M. J. Mol. Biol. 1994, 235, 983-1002. (c) Purisima, E. O.; Nilar, S. H. J. Compt. Chem. 1995, 16, 681-689. (23) (a) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227-249. (b) Warshel, A.; Russell, S. T. Q. ReV. Biophys. 1984, 17, 283-422. (24) (a) Humphreys, D. D.; Friesner, R. A.; Berne, B. J. J. Phys. Chem. 1995, 99, 10674-10685. (b) Smith, K. C.; Honig, B. Proteins 1994, 18, 119-132. (c) Simonson, T.; Bru¨nger, A. T. J. Phys. Chem. 1994, 98, 46834694. (d) Sharp, K. J. Comput. Chem. 1991, 12, 454-468. (e) Gilson, M. K.; Davis, M. E.; Luty, B. A.; McCammon, J. A. J. Phys. Chem. 1993, 97, 3591-3600. (25) (a) Eisenberg, D.; McLaughlin, A. D. Nature 1986, 319, 199203. (b) Wesson, L.; Eisenberg, D. Protein Sci. 1992, 1, 227-235. (c) Ooi, T.; Oobatake, M.; Ne´methy, G.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 3086-3090. (d) Stouten, P. F. W.; Frommel, C.; Nakamura, H.; Sander, C. Mol. Simul. 1993, 10, 97-120. (e) Vila, J.; Williams, R. L.; Va´squez, M.; Scheraga, H. A. Proteins 1991, 10, 199218. (f) Luty, B. A.; Wasserman, Z. R.; Stouten, P. F. W.; Hodge, C. N.; ZaCharias, M.; McCammon, J. A. J. Compt. Chem. 1995, 16, 454-464. (g) Schiffer, C. A.; Caldwell, J. W.; Kollman, P. A.; Stroud, R. M. Mol. Simul. 1993, 10, 121-150. (26) (a) von Freyberg, B.; Richmond, T. J.; Braun, W. J. Mol. Biol. 1993, 233, 275-292. (b) von Freyberg, B.; Braun, W. J. Comput. Chem. 1993, 14, 510-521. (c) von Freyberg, B.; Schaumann, T.; Braun, W. FANTOMsUser’s Manual and Instructions; ETH Zu¨rich: Zu¨rich, 1993. (d) Schaumann, T.; Braun, W.; Wu¨thrich, K. Biopolymers 1990, 29, 679694. (e) von Freyberg, B.; Braun, W. J. Comput. Chem. 1991, 12, 10651076. (f) Prot, G.; Cheng, B.; Gibson, K. D.; Vila, J.; Palmer, K. A.; Nayeem, A.; Maigret, B.; Scheraga, H. A. J. Comput. Chem. 1992, 13, 1-11. (g) Eiaenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. J. Comput. Chem. 1995, 16, 273-284. (27) Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B. Biochemistry 1981, 20, 849-855. (28) Sharp, K. A.; Nicholls, A.; Friedman, R.; Honig, B. Biochemistry 1991, 30, 9686-9697. (29) (a) Li, Z.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 6611-6615. (b) Li, Z.; Scheraga, H. A. J. Mol. Struct. (Theochem) 1988, 179, 333-352. (c) Nayeem, A.; Vila, J.; Scheraga, H. A. J. Compt. Chem. 1991, 12, 594-605. (30) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087-1092. (31) Bru¨schweiler, R.; Roux, B.; Blackledge, M.; Griesinger, C.; Karplus, M. J. Am. Chem. Soc. 1992, 114, 2289-2302. (32) (a) Schmidt, J. M.; Bru¨schweiler, R.; Ernst, R. R.; Dunbrack, R. L., Jr.; Joseph, D.; Karplus, M. J. Am. Chem. Soc. 1993, 115, 8747-8756. (b) Lipari, G.; Szabo, A. J. Am. Chem. Soc. 1982, 104, 4546-4559. (c) Tropp, J. J. Chem. Phys. 1980, 72, 6035-6043. (33) Motta, A.; Picone, D.; Tancredi, T.; Temussi, P. A. J. Magn. Reson. 1987, 75, 364-370. (34) Abseher, R.; Lu¨demann, S.; Schreiber, H.; Steinhauser, O. J. Am. Chem. Soc. 1994, 116, 4006-4018. (35) (a) Keepers, J. W.; James, T. L. J. Magn. Reson. 1984, 57, 404426. (b) Boelens, R.; Koning, T. M. G.; Kaptein, R. J. Mol. Struct. 1988, 173, 299-311. (36) Wilmot, C. M.; Thornton, J. M. J. Mol. Biol. 1988, 203, 221-232.
JP953016Y
