Combined Use of Homo- and Heteronuclear Coupling Constants as Restraints in Molecular Dynamics Simulations DALE F. MIERKE and HORST KESSLER

Organisch Chemisches Institut, Technische Universitat Munchen, Lichtenbergstrasse4, 8046 Garching, Germany

SYNOPSIS

A penalty function for scalar coupling constants has been applied in molecular dynamics simulations as a n experimental constraint. The function is based on the difference between the coupling constant calculated from the dihedral angle and the experimentally measured coupling constant. The method is illustrated on a model cyclic pentapeptide for which J H N . H a and JHN.ce,both about the C#J backbone dihedral angle, have been measured. The function is efficient in producing structures consistent with the scalar couplings, but removed from the conformation observed in solution. This arises from the lack of J restraints for the 4 dihedral angle. Simulations with both nuclear Overhauser effect (NOE) and J-coupling restraints illustrates small but significant differences from simulations using only NOEs. 0 1992 John Wiley & Sons, Inc.

INTRODUCTION The use of nuclear Overhauser effects (NOEs) in molecular dynamic simulations has been employed in the generation of peptide and protein conformations for many years.lP2This method utilizes a penalty function (quadratic or skewed harmonic) for each pair of protons for which an NOE is observed. During the simulation the sum of potential energy function and penalty function is minimized to create a structure or family of structures compatible with the experimental data. To account for the possibility of numerous structures in fast equilibration3 the NOE penalty function has recently been m ~ d i f i e d . ~ It has been known for some time that coupling constants are an additional source of conformational information by use of the Karplus equation and calculation of dihedral angle^.^,^ However, the application of these data as constraints within molecular dynamic simulations has been avoided because of the multiple solutions to the Karplus equation (up to four dihedral angles). One common method is to utilize only extreme coupling constants (i.e., greater than 8.5 or less than 6.0 Hz) and restrain the diRiopolymers, Val. 32, 1277-1282 (1992) 0 1992 John Wiley & Sons, Inc.

CCC 0006-3V25/92/101277-06$04.00

hedral angle to a range of allowed value^.^^^ Although this is easy to implement and has the advantage of only one minimum, assumptions of the appropriate ranges and the lack of possible transitions to the other minima (from the Karplus curve) are drawbacks. A penalty function of the form, &[ 1 cos (@ - S)] (where 0 is the dihedral angle and 6 a phase shift), similar to what is commonly used for the description of dihedral angles, has been applied in molecular dynamics (MD ) simulations? However, to mimic the energy profile of the Karplus equation this method for many coupling constants requires multiple restraints with a careful adjustment of Kd and 6. In addition, the energy barriers between minima will not be exact, which is extremely important, as will be shown below, when more than one coupling constant is available for a dihedral angle. It seems to be more appropriate to apply a penalty function that directly contains the relationship between the coupling constant and dihedral torsion. Such a function has been proposed and applied, but the results were difficult to interpret because of problems with multiple conformations and coupling constants that were not restrictive (i.e., only the homonuclear J H N . H a coupling constants were used). With knowledge of additional scalar couplings about a dihedral angle, including both homonuclear and heteronuclear coupling constants, the allowed tor-

+

1277

1278

MIERKE AND KESSLER

sional space is greatly Here we describe the implementation of a penalty function derived from the Karplus equation and its application to a model cyclic pentapeptide, cyclo [ -Pro ‘-Pro ‘-Ala3Ala4-Ala5-1. This peptide was chosen since both J H N - H a and 3 J ~ ~coupling . ~ p constants have been measured for the 4 dihedral angle and the conformation is well determined by NOES. Another reason for our choice results from the lack of additional problems caused by different populations of side chain rotamers about xl. Hence, we restrict ourselves here to the simplest case assuming one predominant conformation.

ferent atoms.14-16The atomic force F , , from this penalty function required for molecular dynamics simulations can be obtained following standard procedures:

THEORY

SIMULAT I 0N

The penalty function chosen for the coupling constant restraint is similar to that commonly used for NOE restraints lo:

The coupling constant restraint function was tested on a model peptide, cyclo [ -Pro1-Pro2-Ala3-Ala4Ala5-1, which has been extensively examined by nmr.17 The peptide in DMSO exists as two conformations: the major isomer ( 53%) contains cis peptide bonds between Ala5-Pro and Pro ‘-Pro2,while the minor isomer (47% ) contains a cis peptide bond between Pro’-Pro’. For both isomers the homonuclear J H N . H a and heteronuclear JHN.CBscalar coupling constants containing information about the 4 dihedral (see Figure 1) have been measured.” “Random” starting structures were created by application of dihedral angles restraints to the #, of the three alanines to 180”: the cyclic compound cannot obtain these constraints and therefore is of high energy. Although the structures are not random, they are well removed from the structures found in solution, which suits our purposes. The peptide was placed in a periodic truncated

where ri is the Cartesian position vector of the atoms involved in the dihedral angle.



where EJ is the potential energy of the penalty function, K j is the force constant, J i s the coupling constant calculated from the dihedral angle, and Jelp is the experimentally determined scalar coupling constant. The coupling constant is calculated from the dihedral subtended by the coupled atoms 0,using the equations developed by Karplus many years ago5:

J

=A

cos’(0)

+ B cos(0) + C

(2)

where the coefficients A , B , and C have been empirically adjusted for dihedral torsions involving dif-

H

+

0

Figure 1. Schematic representation of the backbone dihedral angles of amino acid residues. Both of the three bond couplings, the homonuclear ’JHN.~* and heteronuclear JHN. cg coupling constants, supply

information about the

I$

dihedral angle.

COMBINED USE OF HOMO- AND HETERONUCLEAR COUPLING CONSTANTS

octahedron of 33.4 nm3 containing 153 DMSO solvent molecules.18 Molecular dynamics were carried out using a step size of 2 fs employing the SHAKE method" with the nonbonded interactions updated every 50 steps with a cutoff distance of 1.0 nm. The simulations were run for 5 ps each a t 1000,750, and then 300 K with a tight coupling to a temperature bath, a relaxation time of 20 fs2" This temperature coupling was relaxed (200 fs) and then the system allowed to equilibrate for 40 ps. The next 100 ps of the simulation were used for the analysis. The force constant was set to 2500 kJ mol-' nmP2and 1.0 kJ mol-' H z - ~for the NOE and coupling constant restraints, respectively. For the major and minor isomers, 18 and 17 NOE distance restraints were utilized. All of the energy minimizations and molecular dynamics simulations were carried out using the GROMOS program."

RESULTS AND DISCUSSION The scalar coupling constants measured for the model compound and the coefficients used in Eq. ( 2 ) are listed in Table 1. In Figure 2 the energy profile of the coupling constants for each of the alanines as a function of the q5 dihedral torsion are plotted. The homonuclear coupling J H N . H n (dashed line) is more than three times larger than the heteronuclear Coupling Constants Table I 3 & ~ . ~ ,and 3JHN.Ca Measured for the Major and Minor Isomers of cyclo[-Pro1-Pro2-Ala3-Ala4-AlaS-] in DMSO

3 J ~ ~ . Torsion" ~ , (Hz) (")

Residue

3 J ~ ~ . c p Torsionb

(H4

(")

2.0

-90, -30 63,177 -60 73,167 -81, -39 68,172

Major ~

1

~

38.0

~

1

~

4

~

1

~

5

6.7

6.6

-153, -87 44,78 -160, -80 32,88 -161, -79 31,89

3.0 2.5

Minor ~

1

~

3

-141, -99

0.2

-160, -120

~

1

~

45.6

2.6

~

1

~

5 10.0

-167, -73 2496 -137, -103

-79, -41 69, 171 -163, -117 -3,43

9.6

-0,40

a

0.4

Values of 9.4, -1.1, and 0.4 were used for A, B, and C in Eq.

(2).

Values of 4.7, -1.2, and -0.2 were used for A, B, and C in Eq. (2).

1279

coupling 3 J H N (dotted .~P line). This is from the smaller range of expected coupling constants and the coefficients ( A , B , C ) used for the latter. Of course, the force constant in Eq. (1) could be adjusted to equalize the energy from the two coupling constants. Despite the lower energy of the heteronuclear coupling, the reduction of the allowed torsional space can be seen from the sum of the two functions (solid line). This is especially clear in Figures 2B and C where the homonuclear 3 J ~ ~ . ~ , produces four well-separated minima, illustrating why such coupling constants (6.7 and 6.6 Hz) are not commonly used as restraints. However, with the inclusion of the 'JHN-CP, three of the minima are excluded. The average values of the dihedral torsions from three different MD simulations are listed in Table 11. During the simulations, ( 1)only the NOEs, ( 2 ) only the coupling constants, or ( 3 ) both of the experimental measurements were used as restraints. The structures from the simulations of the major and minor isomers with the application of the NOEs ( simulation 1) have a small average NOE restraint violation, 8 and 5 pm, for the major and minor isomers, respectively, indicating that the NOEs are consistent with one conformation. The major isomer contains a @II'-turnabout Ala3-Ala4while the minor isomer contains a PVIa turn about Pro'-Pro2.17 The simulations of the major isomer with both NOE and coupling constant restraints (3) illustrate small but significant differences with the results from the simulation with only NOEs ( 1). The q5 torsions of the alanine residues change to values that are minima of the coupling constant restraint term: the largest change is with Ala4,changing from -115" to -80". The average distance restraint violations are identical ( 8 pm) . However, three of the NOEs, all involving the NH of Ala4,show a significant deviation (from simulation 1to 3): Ala4HN-Ala5 HN (208-245 pm) ,Ala4HN-Ala4 H P(295-276 pm) , and Ala4 HN-Pro' H a(357-402 pm) . The first two are in better agreement with the experimentally determined distances while the latter is now 45 pm greater than the upper distance restraint. In addition, the average dihedral angles from simulation 3 are in closer agreement with the standard values for a 011'turn (i.e., the standard $, 4 for Ala3 and Ala4 are 60", -120" and 80", 0", respectively). The rms difference of all of the heavy atoms for the two structures is 37.8 pm. Simulations 1 and 3 of the minor isomer show smaller differences than observed for the major isomer. Again all of the q5 dihedral angles are in minima as determined from the coupling constants and the

75

50

4

B

50

$?

z w

aJ

c 25

25

0

-180

-120

-60

0

60

120

180

-1EO

-120

-60

0

60

120

180

-120

-60

0

60

120

180

125

D 100

75

50

w 25

-180

-120

-60

0

60

120

180

-lEO

4, 125

25

I

-

1 I 1

z

! 1 I 1

I

1

1 I

1

I II

w

I I

;,

:\

0

I

I

\’

-180

I

-120

-60

0

4,

60

I

I

COMBINED USE OF HOMO- AND HETERONUCLEAR COUPLING CONSTANTS

average distance restraint violations from the two simulations are identical. The rms difference of all of the heavy atoms for the two structures is 10.7 pm. The simulations with only coupling constant restraints ( 2 ) produces structures with dihedral angles very similar to those observed from simulations 1 and 3. The exception is Ala5 of the minor isomer where a value of 52" is observed, a different minimum from the coupling constants (see Figure 2F) than found in simulations 1and 3. Despite the similarities of the 6 dihedral angles of the major isomer, the structures are quite different: the rl/ dihedral angles show large deviations. This result is not unexpected since experimental constraints are only applied to the 4 dihedral and MD simulations without restraints (starting from conformations removed from the observed structure) do not converge to structures consistent with the experimental data. Starting with the conformation of the minor isomer from simulation 2 and applying the NOE restraints produces a structure similar to that observed from simulation 3: the 6 dihedral of Ala5 changes from 52" to -108". This is a clear indication that transitions between the different minima of the J restraints are possible (see Figure 2F), which is important for multiple minima functions. Additional simulations were carried out with both NOEs and coupling constant restraints but with the force constant of the heteronuclear coupling constant increased so that the height of the maxima are identical for the two coupling constants (e.g., the force constant of the J H ~ . c awas increased by 3.64, 3.61, and 3.17 for Ala3, Ala4, and Ala5, respectively, for the major isomer). The results from these simulations are identical to those from simulation 3 listed in Table 11. This is not too surprising since the increased force constant will only reinforce the minima shown in Figure 2, observed with the use of the smaller force constant. The importance of using more than one scalar coupling was illustrated with simulations using only the 3 J ~ ~ as .restraints. ~ , For the minor isomer, the Ala4 adopted a 6 value of -157", satisfying the J restraint. Starting from this conformation and application of all of the couplings, the Ala4 adopts a 6

1281

Table I1 Averages of Dihedral Angles from MoIecular Dynamics Simulations of the Major and Minor Isomers of cycle[-Pro1-ProZ-Ala3-Ala4-Ala6-] in DMSO Simulationa Torsion Major Pro'

1

1

*44 * 4

~

1

# 4~

~

1

Pro2 ~

Minor Pro'

~

* * * * *4 * * 4~

4

Pro2 4

~

1 4~

~

1

~

Ala5 4

-72.8 155 -65.0 148 74.3 3 -83.4 4 -115 -19.8 5 -63.3 -42.3 -41.9 130 -79.3 -1.0 3 -96.5 -130 -67.0 4 -39.0 -101 120

2

3

Starting Structure

-36.3 119 -75.8 163 63.0 -91.4 -78.5 -82.6 -84.5 109

-70.7 155 -61.9 142 64.2 -104 -80.7 -21.1 -74.2 -45.7

-126 152 -149 141 180 -68 -168 -104 -177 -120

-43.5 131 -88.3 31.5 -119 -122 -73.4 -39.2 -108 125

-137 141 -121 133 177 -102 -170 94 170 -152

-39.5 131 -72.2 -34.6 -102 -112 -72.0 107 52.0

132

a Values are in degrees. The simulations were run with the application of different experimental constraints: (1)only NOEs, (2) only scalar coupling constants, and (3) both NOEs and coupling constants. For simulation 2, it is important to stress that there are restraints only for the 4 dihedral angles of the alanines.

value of -72", identical with the results from simulation 2 shown in Table 11. The results are in complete accord with the energy profile shown in Figure 2E.

CONCLUSIONS The use of a penalty function for coupling constants has been applied in molecular dynamics simulations.

Figure 2. The potential energy of the coupling constant restraint term as a function of the 4 dihedral angle of ( A ) Ala3, ( B ) Ala4, and ( C ) Ala5 of the major and ( D ) Ala3, ( E ) Ala4,and ( F ) Ala5 of the minor isomers of the cyclic pentapeptide, cyclo[ -Pro'-Pro2-Ala3Ala4-Ala5-1.The energy [using Eqs. ( 1 ) and ( 2 ) with a force constant of 1 kJ mol-' Hz?] for the homonuclear coupling 3JHN.H, (dashed), heteronuclear coupling JHNcs (dotted), and the sum of the two functions (solid) are shown.

1282

MIERKE AND KESSLER

The utility of the function, using the coupling constants directly, has been shown to be quite efficient in the development of structures with minima for the dihedral angles. Such a procedure is a more realistic description of the relationship between the dihedral and coupling constant. This method will become more important as the availability of additional heteronuclear coupling constants increases. An advantage of this method is the simplicity of independently adjusting the energetic contribution (i.e., varying the force constants) of the different J couplings according to the accuracy of the coupling constant and the reliability of the Karplus-type relationship. The Deutsche Forschungsgemeinschaft and Fonds der Chemischen Industrie are gratefully acknowledged for financial support.

REFERENCES 1. Kaptein, R., Zuidenveg, E. R. P., Scheek, R., Boelens, R. & van Gunsteren, W. F. (1985) J . Mol. Biol. 1 8 2 , 179-182. 2. Clore, G. M., Gronenborn, A. M., Brunger, A. T. & Karplus, M. (1985) J.Mol. Biol. 186,435-438. 3. Kessler, H., Griesinger, C., Muller, A., Lautz, J., van Gunsteren, W. F. & Berendsen, H. J. C. (1988) J. Am. Chem. SOC.110, 3393-3398. 4. Torda, A. E., Scheek, R. M. & van Gunsteren, W. F. (1989) Chem. Phys. Lett. 157,289-294. 5. Karplus, M. (1959) J. Chem. Phys. 3 0 , l l - 1 5 . 6. Karplus, M. (1963) J . Am. Chem. SOC.85,2870-2871. 7 . Kline, A. D., Braun, W. & Wuthrich, K. (1988) J . Mol. Biol. 204, 675-694.

8. Moore, J. M., Case, D. A., Chazin, W. J., Gippert, G. P., Havel, T. F., Powls, R. & Wright, P. E. (1988) Science 240, 314-317. 9. de Vlieg, D., Boelens, R., Scheek, R. M., Kaptein, R. &van Gunsteren, W. F. (1986) Isr. J . Chem. 27,181188. 10. Kim, Y. & Prestegard, J. H. (1990) Proteins Struct. Funct. Genet. 8,377-382. 11. Cowburn, D., Live, D. H., Fischman, A. J. & Agosta, W. C. (1983) J . Am. Chem. SOC.1 0 5 , 7435-7442. 12. Kurz, M., Schmieder, P. & Kessler, H. (1991) Angew. Chem. Int. Ed. Engl. 30,1329-1330. 13. Schmieder, P. & Kessler, H. (1992) Biopolymers 3 2 , 435-440. 14. Bystrov, V. F. (1976) Prog. N M R Spectrosc. 1 0 , 4181. 15. Cung, M. T., Marraud, M. & Neel, J. (1974) Macromolecules 7, 606-614. 16. Haasnoot, C. A. G., de Leeuw, F. A. A. M., de Leeuw, H. P. M. & Altona, C. (1980) Tetrahedron 36,27832791. 17. Kurz, M. (1991) Ph.D. thesis, Technical University

Munich. 18. Mierke, D. F. & Kessler, H. (1991) J. Am. Chem. SOC. 113,7466-7470. 19. van Gunsteren, W. F. & Berendsen, H. J. C. (1977) Mol. Phys. 34,1311-1321. 20. Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A. & Haak, J. R. (1984) J . Chem. Phys. 8 1 , 3684-3690. 21. van Gunsteren, W. F. & Berendsen, H. J. C. (1987) Groningen Molecular Simulation Library Manual (GROMOS), Biomos B.V., Groningen.

Received May 22, 1992 Accepted June 4, 1992

Combined use of homo- and heteronuclear coupling constants as restraints in molecular dynamics simulations.

A penalty function for scalar coupling constants has been applied in molecular dynamics simulations as an experimental constraint. The function is bas...
443KB Sizes 0 Downloads 0 Views