J. Phys. Chem. B 2001, 105, 711-714
711
Simulation of Hydrated BPTI at High Pressure: Changes in Hydrogen Bonding and Its Relation with NMR Experiments Massimo Marchi* Section de Biophysique des Prote´ ines et des Membranes, DBCM, DSV, CEA, Centre d’EÄ tudes, Saclay, 91191 Gif-sur-YVette Cedex, France
Kazuyuki Akasaka Department of Molecular Science, Graduate School of Science and Technology, Kobe UniVersity, Kobe 657-8501, Japan ReceiVed: July 18, 2000; In Final Form: September 8, 2000
This paper is focused on the structural changes occurring in hydrated bovine pancreatic trypsin inhibitor (BPTI) at varying hydrostatic pressure. We have carried out a series of nanosecond constant pressure simulations of BPTI at increasing pressure employing two popular force fields, AMBER and CHARMM. Results concerning the hydrogen bonds lengths and distribution underline differences in the two force fields. Overall, CHARMM predicts a compression of the BPTI backbone hydrogen bond network, while AMBER computes a slight increase in the length of these bonds. Compared to the experimental NMR estimate of hydrogen bond shortening, CHARMM, and to a minor extent AMBER, follow the experimental trend in the β -sheet and helix H2 regions of BPTI. Notwithstanding, both force fields show a rather large overall deviation from NMR findings.
I. Introduction During the past few years, considerable effort has been made to elucidate the mechanism for protein denaturation at high pressure1,2 and to understand the origin of the pressure effects on the structure and volume of proteins.3-7 Compared to temperature denaturation, pressure produces only changes in volume (i.e., in interatomic distances among the atoms of the protein), while temperature changes thermal energy which implies a variation of volumes due to thermal expansion. Thus, variations of pressure at constant temperature allow to distinguish between energy and volume effects.6 Experimentally, the effects of pressure on protein structure are investigated by a variety of techniques. In particular, highresolution NMR spectroscopy at varying hydrostatic pressure has provided information on the structural changes occurring in proteins at high pressure.8-12 Of significant importance are changes in the hydrogen bond network under pressure that can be gathered by two-dimensional 1H NMR measurements. These changes have been interpreted by Li et al.10 in bovine pancreatic trypsin inhibitor (BPTI) in terms of a systematic shortening of the NH‚‚‚OdC hydrogen bonds. Although their estimate of hydrogen bond changes is not very accurate, their experimental results can be used as a test of the behavior of modern protein force fields with changing pressure and indicate how the potential parameters should be improved. In the past few years, we have developed stable and accurate molecular dynamics (MD) techniques based on the r-RESPA (reversible REference System Propagation Algorithm)13 class of algorithms and suitable to simulate hydrated biomolecules.14-16 More recently, we have derived new r-RESPA algorithms for * To whom correspondence should be addressed. E-mail: marchi@ villon.saclay.cea.fr.
MD simulations in the constant pressure and temperature ensembles (NPT).17 Motivated by the paper by Li et al.,10 we have investigated the predictive capabilities of modern force fields toward hydrogen bonding in biomolecules at high pressure. In this study we use these new computational techniques, as implemented in our MD simulation package ORAC,16,18 to simulate the high-pressure NMR experiment on hydrated BPTI. We then compare results obtained using two popular force fields due to Cornell et al.19 and MacKerell et al.20 with experiment. Since an alternative interpretation of the NMR shifts observed by Li et al.10 would imply the increase in the number of hydrogen bonds, we have also investigated changes in the hydrogen-bonding network of the solvent and the protein. II. Methods Simulations in the isobaric-isothermal (NPT) ensemble can be carried out with the extended system approach that involves adding extra dynamical variables to the systems to control temperature and pressure.21-24 Although other approaches are possible,25 only techniques of this class are able to generate the desired equilibrium distribution functions. Throughout this study we have used the atom group scaling algorithm fully described in ref 17 and implemented in our MD program ORAC.16,17 This is a five-time-step integration scheme with a Liouvillean separation in three nonbonded shells, with particle mesh Ewald26,27 to handle electrostatic interactions, and constraints on covalent bonds entailing hydrogens. The computational advantage of using constraints on covalent hydrogen bonds is discussed in refs 17 and 18 In order to investigate the effect of pressure on hydrated BPTI, we ran a series of simulations by preparing and equilibrating the system in an identical manner for both force fields. In
10.1021/jp002539p CCC: $20.00 © 2001 American Chemical Society Published on Web 12/30/2000
712 J. Phys. Chem. B, Vol. 105, No. 3, 2001
Marchi and Akasaka
TABLE 1: Computed Xrms Relative to Cr’s (〈XrmsCA〉) Obtained from the Four Constant Pressure Runsa force field pressure 〈XrmsCA〉
AMBER P ) 0.1 0.75
P ) 500 0.89
CHARMM P ) 0.1 0.88
P ) 500 0.93
a Lengths are in Å and pressures in MPa. The statistical maximum error is typically less than 0.4 %.
particular, our system consisted of one molecule of BPTI surrounded by 2346 TIP3P28 water molecules in a cubic periodic box of side length set, initially, to 43 Å. After a brief structural minimization of the X-ray structure,29 we heated the system gradually from 20 to 300 K in four steps of 10 ps each with microcanonical runs. Following this phase, we equilibrated for 50 ps in the NPT ensemble at P ) 0.1 MPa and T ) 300 K. Two runs at low (P 0.1 MPa) and high (P ) 500 MPa) pressure were then started from the same final configuration of the equilibration run. At low pressure, we continued the run and acquired statistical data for one additional nanosecond. At high pressure instead, the run was started by imposing a new pressure and equilibrating for 50 ps. After this initial stage, an additional simulation 1 ns long was run. For all production runs, coordinates and cell parameters were recorded every 240 fs for statistical analysis. This amounts to a total of 4166 configurations being saved for each 1 ns run. All simulations discussed in this paper were carried out on a CRAY T3E 600 machine with the parallel version of ORAC and required in total 4 days of CPU time on 16 processors including calculations for both force fields.
Figure 1. Backbone hydrogen bond (HB) distances in BPTI. The two overlapping histograms give the NMR maximum and minimum HB lengths. Triangles and diamonds are results from simulations with the AMBER and CHARMM force fields, respectively. For comparison we also report HB estimated from the crystallographic data.
III. Results and Discussion In Table 1 we report the deviation of the simulated BPTI averaged structures from X-ray coordinates29 (or Xrms). These averages were computed for each run from the recorded coordinates which were rotated and translated so as to obtain a maximum overlap with the CR of the reference X-ray structure. Confirming our previous results on BPTI crystal,30 the AMBER force field at atmospheric pressure reproduces the BPTI structure better than the CHARMM force field with differences in Xrms of about 0.1 Å. At higher pressure both force fields predict a deviation from X-ray data at P ) 0.1 MPa of about 0.9 Å. In this work, hydrogen bonds occurring between an acceptor (AA-A) and donor (D-H) are arbitrarily defined according to two parameters: (i) the distance dH between acceptor (A) and donor (D) atoms; and (ii) the angle θH between the acceptor and donor axis, i.e., AA-A and H-D. For each force field and pressure investigated here the calculation of the backbone hydrogen bond lengths was carried out on the corresponding average protein structure. These configurations were obtained by averaging the protein coordinates over the whole length of each run. As for the Xrms calculation, translational drifts and rotations of the protein occurring during the simulations were eliminated in the process. Before discussing the effect of pressure on the protein hydrogen bonds (HB), we first present a comparison between the calculated average BPTI backbone NH‚‚‚OdC HB lengths at P ) 0.1 MPa with the NMR31 and X-ray results29,35 Since hydrogens coordinates are not available from X-ray, we generated the hydrogen position using purely geometrical criteria. These coordinates were used to compute the hydrogen bond length in Figure 1. Only bonds with lengths dH < 3 Å and with the N-H‚‚‚O and H‚‚‚OdC angles θH > 90° are shown in Figure 1. We first notice that the AMBER and CHARMM HB lengths are close to the NMR data within the β-sheet and the last helix
Figure 2. Comparison between computed and experimental estimates of backbone hydrogen bond shortening. Only the hydrogen bonds within the chosen parameters in at least one of the two simulations are shown.
(H2), while differences are larger in the first helix region (H1). For both force fields the β-sheet region fares better, except near B-turns. Here, the HB of ala27 and ala16 were off scale for CHARMM and AMBER, respectively, and that of leu29 for AMBER was 0.4 Å longer than the NMR result. On a rootmean-square (rms) sense CHARMM and AMBER perform equally, with an rms deviation of 0.3 Å from the NMR data for both. We point out that differences between estimated HB obtained from X-ray and NMR are about 3 times smaller. We compare our results at different pressure with backbone HB compression data obtained by 1H NMR chemical shift measurements at P ) 0.1 and 200 MPa. This experimental study relates the HB “formation shift” with the shortening of the bond through an empirical cubic dependence from the H‚‚‚O distance of the NsH‚‚‚OdC10 hydrogen bond. In spite of the uncertainties, this approach is able to clearly relate to well defined structural changes in the protein structure. In Figure 2 we compare the HB shortening attained by CHARMM and AMBER with the experimental estimate. Here, we accounted for differences in experimental (200 MPa) and simulated (500 MPa) pressure by assuming a linear compression regime and uniformly rescaling the NMR HB lengths. The choice of a simulated pressure higher than that of the experiments was made to minimize the statistical error on the average distance changes between structures computed in our investigation.
Simulation of Hydrated BPTI at High Pressure
J. Phys. Chem. B, Vol. 105, No. 3, 2001 713
Figure 3. Radial distribution functions for proteinsprotein (dpp(r), bottom) and protein-water (dpw(r), top) hydrogen bonds obtained from simulations with the AMBER and CHARMM force fields. Distribution functions obtained with AMBER are in parts a and c, while parts b and d display results for CHARMM.
We first note a similar trend in the results between the two force field, although fluctuations up and down the distance change scale are much more pronounced for AMBER than for CHARMM. Averaging over the whole backbone hydrogen bonds included here, we find that CHARMM and AMBER predicts bond changes of -0.027 (compression) and 0.007 (slight expansion) Å, respectively. This compares with an experimental estimate of -0.042 Å at 500 MPa (corresponding to -0.017 Å at 200 MPa). Moreover, both force fields predict an expansion in the H1 helix region (residues cys5 and leu6) and close to strands and B-turns (residues gly28, tyr35 and asn44), while experiment estimates a compression. On the contrary, simulation predicts large compression (larger than experiment) for the glu7 residue. Finally, the experimental shortening trend in the β-sheet and helix H2 regions is followed by the CHARMM and, to a minor extent, by the AMBER result, but the overall rms deviation from experiment is rather high for both, i.e., 0.09 and 0.14 Å for CHARMM and AMBER, respectively. Previous high-pressure simulations32,33 have pointed out that while hydrogen bonding between protein and water increased with pressure, it decreased within the protein itself. In order to investigate this pressure effect on the hydrogen bonding network, we have computed the radial distribution functions (d(r)) of the HB's at P ) 0.1 and 500 MPa. These functions are defined such that d(r) dr gives the number of HB's contained in the radial shell between r and r + dr. As shown by the dpw(r)'s in Figure 3a,b, both CHARMM and AMBER predict a similar increase of the protein-water hydrogen bonding. The picture is different for the corresponding protein-protein HB's. Indeed, AMBER in Figure 3c shows an increase, although little, in hydrogen bonding with increasing pressure. Conversely, CHARMM predicts a decrease of dpp(r) at high pressure (see Figure 3d). In both cases, however, pressure effects on the hydrogen bonding within the protein itself are of minor importance with respect of induced changes in the proteinwater interactions. Since an alternative interpretation of the NMR shifts observed by Li et al. would imply the increase of the number of hydrogen bonds, our finding in Figure 3, a and c, gives support to their analysis in term of HB shortening. We point out here that the protein-protein HB radial distribution functions discussed before include contributon not only from backbone HB’s but also from all intramolecular contributions.
Figure 4. CR difference distance matrix for BPTI computed with the AMBER (top) and CHARMM (bottom) force fields.
Because of the small changes in the intraprotein hydrogen bonds at varying pressure, it is of interest to look at the structural changes in BPTI as computed with the two force fields. In doing so, we single out the behavior of distinct regions within the protein induced by pressure changes by computing the difference distance matrix defined as
dij ) |rPi 1 - rPj 1| - |rPi 2 - rPj 2|
(1)
where rPi 1 is the average position of the ith atom at pressure P1. If a portion of the protein is invariant under pressure, the elements of the difference distance matrix corresponding to distances within the invariant structure will be close to zero. Instead of plotting dij in three dimensions, we prefer to use a two-dimensional representation where different colors correspond to values of dij within certain arbitrary thresholds. Moreover, for clarity and simplicity, in the computation of the difference distance matrix we have considered only distances between CR’s atoms. In Figure 4 we show the matrix dij computed from the BPTI average structures obtained with CHARMM and AMBER at P ) 0.1 and 500 MPa. We first notice that the compression patterns of the two force fields have both similarities and differences. Indeed, AMBER predicts the least compression of BPTI and its difference matrix show much larger regions of yellow and white (both signaling expansion) than that of CHARMM. At the same time, we observe that black and red regions of high compression are larger again in the AMBER difference matrix. Both results are consistent with the HB compression data shown in Figure 2 which shows large fluctuations in the HB compression/expansion along the protein and an overall small expansion for AMBER. Similarities between the two results can also be observed. Both for AMBER and CHARMM, the final two residues belonging to helix H2 compress strongly with respect to the β-sheet regions and the loop region before helix H2. In addition, the distance between the H2 helix and the β-sheets decreases upon compression for both potentials. For CHARMM this effects is evinced by the two green regions on top of the plot, while for AMBER the same two regions are interlaced alse with
714 J. Phys. Chem. B, Vol. 105, No. 3, 2001 red and black. A similar behavior is seen only in the CHARMM calculation also for the H1 helix. IV. Conclusion In this paper we have tackled the problem of how hydrogen bond lengths and distributions change upon hydrostatic compression. Our work has been motivated by recent NMR highpressure experiments on BPTI. Our calculations, not showing any significant increase in the number of intra-protein hydrogen bonds, gives support to Li et al.’s interpretation of the 1H NMR shifts as distance shortening, not increase in the number of hydrogen bonds. While structural deviations of BPTI crystal at increased and decreased pressure were seen by others for the AMBER force field,34 our computation has shown clear differences in the high pressure results obtained with AMBER and CHARMM. Although both potentials show slight changes in the radial distribution function of the intra-protein hydrogen bonds, only CHARMM predicts an overall shortening of these bonds. On the other hand, at room pressure, the backbone hydrogen bond lengths, computed with the two force fields, were for both in rather poor agreement with the lengths measured by NMR and X-ray. From our simulation and previous unrelated studies35 it is clear that the accuracy of modern force fields, among which CHARMM and AMBER, is not sufficient to study quantitatively structural changes of biomolecules. New developments in the potential functions, possibly involving polarizability, are therefore needed. We hope that our study will provide some of the necessary framework for this work to be undertaken. References and Notes (1) Weber, G. J. Phys. Chem. 1993, 97, 7108. (2) Peng, X. D.; Jonas, J.; Silva, J. L. Biochemistry 1994, 33, 8223. (3) Millero, F.; Ward, G.; Chetirkin, P. J. Biol. Chem. 1976, 251, 4001. (4) Gekko, K.; Noguchi, H. J. Phys. Chem. 1979, 83, 2706. (5) Gavish, B.; Gratton, E.; Hardy, C. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 750. (6) Frauenfelder, H.; Alberding, N.; Ansari, A.; Braunstein, D.; Cowen, B.; Hong, M.; Iben, I.; Johnson, J.; Luck, S.; Marden, M.; Mourant, J.; Ormos, P. et al. J. Phys. Chem 1990, 94, 1024. (7) Paci, E.; Marchi, M. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 11609. (8) Samarasinghe, S. D.; Campbell, D. M.; Jonas, A.; Jonas, J. Biochemistry 1992, 31, 7773.
Marchi and Akasaka (9) Akasaka, K.; Tezuka, T.; Yamada, H. J. Mol. Biol. 1997, 271, 671. (10) Li, H.; Yamada, H.; Akasaka, K. Biochemistry 1998, 37, 1167. (11) Li, H.; Yamada, H.; Akasaka, K. Biophys. J. 1999, 77, 2801. (12) Akasaka, K.; Li, H.; Yamada, H.; Thoresen, R. L. T.; Woodward, C. K. Protein Sci. 1999, 8, 1946. (13) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (14) Procacci, P.; Marchi, M. J. Chem. Phys. 1996, 104, 3003. (15) Procacci, P.; Darden, T.; Marchi, M. J. Phys. Chem. 1996, 100, 10464. (16) Procacci, P.; Darden, T.; Paci, E.; Marchi, M. J. Comput. Chem. 1997, 18, 1827. (17) Marchi, M.; Procacci, P. J. Chem. Phys. 1998, 109, 5194. (18) Procacci, P.; Marchi, M. In AdVances in the Computer Simulations of Liquid Crystals; Zannoni, C., Pasini, P., Eds.; Kluwer Academic: Dodrecht, The Netherlands, 1999; Proceedings of the NATO-ASI School, Erice 11-21 June 1998; p 333. (19) Cornell, W. D.; Cieplak, P.; Bavly, C. I.; Gould, I. R.; K. M. M. Jr., Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollmann, P. J. Am. Chem. Soc. 1995, 117, 5179. (20) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Gao, H.; Ha, S.; JosephMcCarty, D.; Kuchnir, L. et al. J. Phys. Chem. 1998, B102, 3586. (21) Andersen, H. J. Chem. Phys. 1980, 72, 2384. (22) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (23) Nose´, S. J. Chem. Phys. 1984, 81, 511. (24) Hoover, W. Phys. ReV. A 1985, 31, 1695. (25) Berendsen, H.; Postma, J.; van Gunsteren, W.; DiNola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684. (26) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (27) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 101, 8577. (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. JCP 1983, 79, 926. (29) Wlodawer, A.; Deisenhofer, J.; Huber, R. J. Mol. Biol. 1987, 193, 145. (30) Ceccarelli, M.; Marchi, M. J. Phys. Chem. 1997, 101, 2105. (31) Berndt, K. D.; Guntert, P.; Orbons, L. P. M.; Wuthrich, K. J. Mol. Biol. 1992, 227, 757. (32) Kitchen, D. B.; Reed, L. H.; Levy, R. M. Biochemistry 1992, 31, 10083. (33) Wroblowski, B.; Diaz, J. F; Abd, K. H.; Engelborghs, Y. Proteins 1996, 25, 446. (34) York, D. M.; Darden, T. A.; Pedersen, L. G. In Modelling of Biomolecular Structures and Mechanisms; Pullman, A., Ed.; Kluwer Academic Press: Dordrecht, The Netherlands, 1995; p 203. (35) Beachy, M. D.; Chasman, D.; Murphy, R. B.; Halgren, T. A.; Friesner, R. A. J. Am. Chem. Soc. 1997, 119, 5908. (36) Since hydrogen coordinates are not available from X-ray, we generated the hydrogen position using purely geometrical criteria. These coordinates were used to compute the hydrogen bond length in Figure 1.
