8050
J. Phys. Chem. B 2001, 105, 8050-8055
Dissecting the Vibrational Entropy Change on Protein/Ligand Binding: Burial of a Water Molecule in Bovine Pancreatic Trypsin Inhibitor Stefan Fischer and Jeremy C. Smith Biocomputing, IWR, Im Neuenheimer Feld 368, UniVersita¨ t Heidelberg, D-69120 Heidelberg, Germany
Chandra S. Verma* Structural Biology Laboratory, Department of Chemistry, UniVersity of York, York YO10 5DD, U.K. ReceiVed: June 1, 2001
Using normal-mode analysis, the vibrational entropy change on the burial of a crystallographically wellordered water molecule in bovine pancreatic trypsin inhibitor (BPTI) is dissected. The vibrational entropy content of the complex is 13.4 cal mol-1 K-1 higher than that of the unbound protein. A detailed analysis is performed of how the translational and rotational degrees of freedom of the isolated water molecule are transformed into vibrational modes in the complex. This process is shown to be well described by a model of the complex in which the water molecule librates in a rigid protein cage. These librational modes contribute 9.4 cal mol-1 K-1 to the entropy change. The remaining 4 cal mol-1 K-1 arises from increased protein flexibility due to softening of the delocalized modes, mostly in the frequency range below 50 cm-1. The dominant librational entropy effect suggests a method by which an estimation of the vibrational contribution to ligand binding can be efficiently computed.
Introduction Understanding the binding of small molecules to proteins is of both fundamental and practical interest. The ubiquitous water-protein interaction can be considered an important special case of this and has been the focus of several experimental and theoretical studies. 1-7 In principle, the small size of the water molecule should simplify computational analysis of the underlying thermodynamics. Of particular interest are buried water molecules that fill some of the voids in protein structures.7-11 Characterization of the binding structure, dynamics, and thermodynamics of well ordered, strongly interacting, buried water molecules is likely to provide much information that will also be of use in understanding ligand binding in general. Central to this is an understanding of the entropy change on binding, due to changes in the nature of the degrees of freedom of the ligand. An example of a well-ordered water molecule is found in the bovine pancreatic trypsin inhibitor (BPTI). This deeply buried water molecule (identified as W122 in the PDB entry 5pti12) makes four strong hydrogen bonds with the protein and has served as a probe for investigating the dynamics and thermodynamics of structural water in proteins.13-18 In previous theoretical work,17 the binding of W122 to BPTI was shown to be accompanied by an increase in the total vibrational entropy of the complex relative to the unbound protein. Here the origin of the vibrational entropy change accompanying the proteinwater interaction is dissected. When a water molecule binds, its nine degrees of freedom (thrice the number N of atoms) are transformed from six translational and rotational degrees of freedom and three internal vibrational modes into nine vibrational modes in the complex.19,20 To fully understand the entropy * Corresponding author. E-mail: 44-1904-410519.
[email protected]. FAX:
change, we need to distinguish between the entropy content of the newly introduced vibration modes and that due to changes in the vibrational frequencies of the protein which reflect an intrinsic change in the protein flexibility. To do this, the forms and frequencies of the nine new modes must be determined. It is shown here that the libration modes of the bound water contribute two-thirds of the vibrational entropy increase. The remaining increase is due to an increase in the flexibility of the protein, which is shown to be mostly a softening of the protein vibration modes in the low-frequency range below 50 cm-1. Methods All calculations were performed with the program CHARMM22 using the parameter set 1923 with explicit hydrogens on all aromatic groups.16 The structures of the unbound protein (P) and the complex (PW) with water W122 (W) were prepared by energy minimization, as previously described.17 The present study simulates the behavior in the aqueous phase by applying a distance-dependent dielectric constant ( ) rij) to approximate the solvent-induced screening of the electrostatic interactions and by using standard solution protonation states for the ionizable side-chains. Electrostatic interactions were shifted to zero at 12 Å, and the Lennard-Jones interactions were switched to zero between 8 and 12 Å. The calculated structural changes of BPTI upon binding W122 are small (the RMS deviation from the unbound structure is 0.13 Å over the non-hydrogen atoms) and are localized in the region around W122. The change in the potential energy is ∆E ) -19.8 kcal/mol, reflecting the formation of the four good hydrogen bonds, tetrahedrally arranged around W122, as found in the crystal structures. The vibrational frequencies, νi of each of the three species involved in the reaction P + W f PW were obtained from a normal-mode analysis24 as outlined earlier.17 They correspond to the normal modes of the energy-
10.1021/jp0120920 CCC: $20.00 © 2001 American Chemical Society Published on Web 08/01/2001
Bovine Pancreatic Trypsin Inhibitor
J. Phys. Chem. B, Vol. 105, No. 33, 2001 8051
Figure 1. Cumulative change of the vibrational entropy ∆S(n)vib upon the binding of water W122 to BPTI (as defined in Results section). Frequencies up to 1000 cm-1 are shown, while the insert shows the whole frequency range. The arrows indicate the location of the frequencies of the nine modes given in Table 1.
minimized structures, obtained by diagonalizing the massweighted matrix of the second derivatives of the energy.25 Two normal-mode analyses of the complex were performed. One was with both the protein and the ligand flexible. In the second, the protein was kept fixed and the ligand flexible, for the purpose of characterizing the libration modes of W122. Results Figure 1 shows the cumulative change in the vibrational entropy of the complex (PW) relative to that of the unbound protein (P): n
∆S(n)vib )
∆Si ∑ i)1
(1)
where n is the number of modes in the complex below a given frequency, ν. ∆Si ) Si(PW) - Si(P) is the contribution of the ith vibrational mode, where the entropy content of vibrational mode i is given by19
Si ) (hνi/T)/(ehνi/kT -1) - k ln(1 - e-hνi/kT)
(2)
As higher frequencies are included gradually in the sum, ∆S(n)vib converges to the total vibrational entropy change ∆Svib ) 13.36 cal mol-1 K-1. Although ∆S(n)vib reaches this final plateau value at a frequency of about 1500 cm-1, most of the increase originates from frequencies below 1000 cm-1. ∆Svib does not depend strongly on the conformation of the protein, as shown in our previous study,17 where the standard deviation of ∆Svib over 100 conformers obtained from molecular dynamics was only 2 cal mol-1 K-1. The normal-mode analysis of W122 in the force-field of the fixed protein results in the nine frequencies ωk listed in Table 1, separated into librational modes (k ) 1-6) and internal modes (k ) 7-9). The six libration modes can be viewed as a mixture of rotational and translational motions of the water molecule in its protein “cage” that replace the six lost translational and rotational degrees of freedom of the unbound water. The three internal modes involve deformations of the water molecule. This
TABLE 1: Vibration Modes of Water W122 in Fixed BPTIa modeb k
frequency ωk, cm-1
entropy Svib,c cal mol-1 K-1
1 2 3 4 5 6 7 8 9
102 204 261 366 497 631 1776 3357 3434
3.42 2.11 1.66 1.11 0.67 0.41 0.004 3.4 10-6 2.5 10-6 totald: 9.38
a For the fixed protein analysis, the diagonalization of H, the massweighted Hessian matrix computed for the complex, is reduced to diagonalizing the 9 × 9 block matrix composed of elements of Hij, where i and j are coordinate indices of the atoms of W122. b Modes 1-6 are librational modes, and modes 7-9 are internal modes (angle bending and symmetric and antisymmetric bond stretch); see text. c Vibrational entropy content of each mode. See Results section for the formula of Svib. d Sum of vibrational entropies over all nine modes, SW122 ) Slib + Sinternal.
distinction between libration and internal modes is only possible for a ligand which has very stiff (i.e., high-frequency) internal modes compared to the libration modes. For more flexible ligands, the libration and internal modes are mixed. It can be seen in Table 1 that the internal modes indeed have much higher frequencies than the libration modes, which justifies the distinction between libration and internal modes in the present study. The internal mode frequencies are close to those of an isolated water molecule (1737, 3323, and 3370 cm-1 for the TIP3P model of water used here). To evaluate the contributions to ∆Svib from the new modes in the complex, we also list in Table 1 the vibrational entropy content of each of the nine modes. The sum over all nine modes contributes SW122 ) Slib + Sinternal ) 9.38 cal mol-1 K-1 to ∆Svib. Most of this comes from the six low-frequency libration modes, the contribution from the internal modes being negligible due to their high frequency. We note here that the vibrational entropy of W122 in the field of the fixed protein is rather insensitive to the treatment of the electrostatic interactions or truncations. A
8052 J. Phys. Chem. B, Vol. 105, No. 33, 2001
Fischer et al.
Figure 2. Contributions from modes in the complex to each of the six libration modes of W122 listed in Table 1. The projections pi ) wk‚vi/ (|wk|‚|vi|) are plotted, where wk is one of the six libration modes (whose frequency ωk is indicated by an arrow on the νi scale), and vi is the sub-vector with components of W122 in the ith normal-mode vector of the complex (of frequency νi). (a) w1, ω1 ) 102 cm-1, (b) w2, ω2 ) 204 cm-1, (c) w3, ω3 ) 261 cm-1, (d) w4, ω4 ) 366 cm-1, (e) w5, ω5 ) 497 cm-1, (f) w6, ω6 ) 631 cm-1. Note that because the 3N-6 vectors vi are of dimension 9 (three times the number of atoms in W122), they cannot be all orthogonal so that Σipi2 is not unity.
dielectric constant ) 1 yields SW122 ) 9.8 cal mol-1 K-1, close to the above value obtained with ) rij. In Figure 1, it is possible to distinguish two frequency regions: (1) a region below 50 cm-1 in which ∆S(n)vib increases by about 4 cal mol-1 K-1 and (2) a region between 50 cm-1 and 800 cm-1 in which ∆S(n)vib increases by most of the remaining 9.4 cal mol-1 K-1. The second region covers a frequency range similar to that containing the six water libration modes in the fixed protein analysis (102 to 631 cm-1; see Table 1), and the rise in ∆S(n)vib over this region is about 9 cal mol-1 K-1, corresponding closely to the entropy content of these modes (Slib ) 9.38 cal mol-1 K-1; see Table 1). This suggests that the increase of ∆S(n)vib in this frequency region is due to the introduction of the six new librational modes in the complex, which will be confirmed below. No increase in ∆S(n)vib is seen for the frequency range of the internal water modes, 1700 to 3500 cm-1, because their entropy content is very small. The increase in ∆S(n)vib does not occur in discrete jumps at the libration frequencies listed in Table 1 because the libration modes are coupled with protein modes of similar frequency, so their entropy contribution is smeared out. This is demonstrated in Figure 2, which shows which modes in the complex couple to each of the six libration modes of W122. Each panel in Figure 2 was obtained by projecting the vector of one of the six libration
modes from Table 1 onto every normal-mode vector of the complex (after removing all vector components other than those of W122). It can be seen that a range of modes from the complex project significantly onto each libration mode of W122. The frequency distribution of these modes is centered approximately at the associated librational frequency (indicated by an arrow). For example, the 2nd libration mode (k ) 2), whose frequency is 204 cm-1, is excited mostly in modes of the complex with frequencies ranging from ∼100 to 300 cm-1 (see Figure 2b). This means that, within that frequency range, there are modes in the fully flexible complex which involve a motion of W122 that is similar to the motion of the 2nd libration mode of W122 in the fixed protein. Panels 2a-f show that the frequencies of the complex which participate in the six librations of W122 cover the whole range from 50 to 800 cm-1, over which they are distributed evenly. This confirms that the 9 cal mol-1 K-1 increase of ∆S(n)vib over the 50 to 800 cm-1 range (see Figure 1) arises from the six libration modes and explains why this increase occurs smoothly over that frequency range. Now that the entropy contributions from the nine new modes in the complex have been accounted for, the remaining 4.0 cal mol-1 K-1 increase of ∆S(n)vib in the frequency region below 50 cm-1 can be fully attributed to a softening of these lowfrequency protein modes. This softening can be seen in Figure
Bovine Pancreatic Trypsin Inhibitor
J. Phys. Chem. B, Vol. 105, No. 33, 2001 8053
Figure 3. Relative change in the protein vibration frequencies upon binding W122. Each normal-mode vector ui of the uncomplexed protein is paired with the normal-mode vector cj of the complex, cimax, having the largest projection pj ) ui‚cj/(|ui|‚|cj|). The frequency change is then taken as ∆νi ) ν(cimax) - ν(ui), whose relative value ∆νi/νi is plotted as a function of νi ) ν(ui) (i.e., the frequency of ui).
Figure 4. Cumulative change in the vibrational protein entropy ∆S(n)vib resulting from the frequency shifts of the protein modes shown in Figure 3. ∆S(n)vib ) Σni ∆Si is computed as in Figure 1, but with ∆Si ) Si(νi+∆νi) - Si(νi) (see caption of Figure 3 for definitions of νi and ∆νi), and n is the number of modes in the unbound protein with a frequency below ν.
3, which shows the relative change in the frequencies of the protein upon binding W122. Figure 3 was calculated by assigning by projection each mode in the uncomplexed protein to its corresponding mode in the complex and taking the difference in the frequencies of these two modes. It can be seen that most protein modes below 50 cm-1 undergo a decrease in frequency ranging between -36% and -4%. The relative change in frequency for modes above 200 cm-1 is negligible. The overlap between each mode of the unbound protein and its corresponding mode in the complex (as measured by their projection pimax, see caption of Figure 3) is high for most modes, as indicated by an average over all modes 〈pimax〉i ) 0.86 close
to unity (the standard deviation is 0.16, and the smallest value of pimax is 0.34). A value of pimax close to 1 means that the nature of the motion in a protein mode is not significantly altered by the binding of W122. The change in the vibrational entropy of the protein caused by the decrease in the protein frequencies, ∆Sprot, is expected to be similar to the entropy obtained after subtracting the entropy of bound W122 (SW122 from Table 1) from the total change in the vibrational entropy (∆Svib, see above). This is verified in Figure 4, which shows the cumulative increase in the vibrational entropy of the protein resulting from the frequency shifts of the protein modes, as calculated in Figure 3. It can be seen that
8054 J. Phys. Chem. B, Vol. 105, No. 33, 2001 the entropy increase comes mostly from modes below 100 cm-1 and reaches a total value of ∆Sprot ) 4.6 cal mol-1 K-1, close to the value of ∆Svib - SW122 ) 4.0 cal mol-1 K-1. Discussion Mass spectrometric experiments have been performed on the binding of a single water molecule to fully dehydrated BPTI.15 The value of the total entropy change thus obtained is ∆S ) -62 ( 5 cal mol-1 K-1. In a previous gas-phase computational study,17 the total change in entropy was calculated as the change in vibrational entropy, ∆Svib ) 11.7 cal mol-1 K-1, minus the translational and rotational entropy of unbound water, Strans ) 34.7 cal mol-1 K-1 and Srot ) 10.6 cal mol-1 K-1 (which was erroneously given as 33 cal mol-1 K-117), using the standard formulas.19,20 The resulting ∆S is -33.6 cal mol-1 K-1, significantly smaller than the above experimental value. Further calculations were performed by Mao and co-workers using a similar approach, also yielding values for the enthalpy and entropy changes that do not simultaneously match those determined experimentally.18 Notwithstanding, ref 18 does show that the results for the gas phase are very sensitive to the conformation of the protein, to the protonation state of the ionizable side-chains, and to the location of the binding-site. Unfortunately, little experimental information is available on these aspects to guide theoretical studies in the gas phase. The conformation of the fully dehydrated protein may significantly differ from the highly hydrated crystal state, particularly near the surface, and the protonation states used normally for simulations in the aqueous phase (Glu-, Asp-, Lys+, Arg+) may not be appropriate in the gas phase. Moreover, while W122 can be localized unambiguously in the solvated protein from several crystal structures,21 the location of the preferred binding site of the first water molecule binding to gas-phase dehydrated BPTI is unknown. These uncertainties will have to be resolved before theoretical studies allow a detailed interpretation of the experimental data in the gas phase. The present calculations for the aqueous phase confirm that the binding of W122 with formation of four near-optimal H bonds to BPTI does induce an increase in the protein flexibility for the delocalized vibrations below 50 cm-1. The value of SW122 obtained here (9.4 cal mol-1 K-1) is comparable to the standard entropy of ice, Sice ) 9.9 cal mol-1 K-1,26 indicating that the mobility of the tetrahedrally coordinated W122 in its protein “cage” is similar to that of ice. This is consistent with the experimental order parameters of W122 in BPTI measured by NMR dispersion, from which libration amplitudes have been derived which are essentially the same as those of water molecules in ice Ih.21 That the mobilities of W122 and water in ice are similar was also shown in reaction path calculations of the exchange of the two water hydrogens, which was found to occur in the same manner for W122 in BPTI as for water in ice Ih, i.e., a two-step process rather than a simple C2-flip.16 However, for structural water that is coordinated differently than W122, for example, making fewer H-bonds or interacting with ionized groups, Slib could vary significantly. The present calculations were performed in the harmonic approximation, in which anharmonic features of the potential surface are neglected. One of these anharmonic effects might be conformational locking of the groups directly binding W122, which would give a negative entropy contribution (i.e., favoring dissociation). Assuming, however, that the changes in the conformational entropy of BPTI upon binding can be neglected, then ∆Svib gives an estimate of the standard entropy of W122. The value calculated here of 13.4 cal mol-1 K-1 lies between
Fischer et al. the standard entropies of ice (see above) and of liquid water, Sliq ) 16.7 cal mol-1 K-1,27 as predicted by Dunitz.28 However, the present study shows that the higher entropy compared to that of water in ice is due not to an increased librational mobility as proposed in ref 28, but rather to a change in protein flexibility. This “excess” protein entropy, ∆Sprot ) ∆Svib - SW122 ) 4 cal mol-1 K-1, lowers the entropic cost of transferring a water molecule from the liquid to the BPTI cavity. The resulting entropy of binding, ∆Sprot + (SW122 - Sliq) ) -3.3 cal mol-1 K-1, raises the free energy of binding at 300 K by only 1 kcal/ mol. It remains to be investigated whether this value is common to water molecules tetrahedrally coordinated by neutral groups in different proteins, since the vibrational response of a protein may depend on its particular structure. The agreement between two ways of computing the increase in the vibrational entropy of the protein, one based on the frequency shifts in the fully flexible protein (4.6 cal mol-1 K-1, see Figure 4) and the other based on taking ∆Svib - SW122 (4.0 cal mol-1 K-1), shows that the entropy content of the nine degrees of freedom of W122 in the complex is well approximated by the vibrational entropy obtained by placing W122 in the force-field of fixed BPTI (i.e., SW122 in Table 1). This approach of placing a ligand in the force-field of the fixed protein significantly reduces the computational expense in determining the vibrational entropy of the complex. It should therefore be useful as a first approximation in calculations of the vibrational contribution to binding free energies in cases where the vibrational change is dominated by the ligand (i.e., ∆Sprot is small relative to ∆Svib, as seen here for W122/BPTI) or where ∆Sprot can be expected to be similar for different ligands binding to a same protein site. This approximation has been applied successfully in accurate calculations of the relative binding affinities for trypsin of a series of benzamidine analogues.29 Note Added in Proof Instead of resulting from an increase in protein flexibility, the decrease in protein vibration frequencies observed in Figure 3 could have resulted from the higher mass of the complex due to the presence of water W122. This is not the case, as can be shown by repeating the normal mode analysis of the flexible complex with vanishing masses for the W122 atoms and using the resulting frequencies to calculate ∆Svib. The effect of this is to remove the influence of the ligand mass on the protein frequencies. Moreover, it removes from ∆Svib the entropy pertaining to motions of the ligand atoms in the complex, since the vibrational entropy of a massless particle is 0 (see eq 2 with infinite frequency ν). This yields the intrinsic change in protein vibrational entropy, ∆Sprot ) +3.79 cal/mol/K. This value is positive, thus proving that the flexibility of the protein does indeed increase upon binding of W122. It is close to the estimates of ∆Sprot computed above in different ways (4.0 and 4.6 cal/mol/K), thus showing that the mass effect of W122 on protein frequencies is small. The difference between this value and ∆Svib calculated with the normal mass of W122 (+13.36 cal/mol/K, see Results) yields the intrinsic contribution of the ligand to the vibrational entropy of the complex, SW122 ) 9.58 cal/mol/K, which is close to the value obtained by fixing the protein (9.38 cal/mol/K, see Table 1). Computing ∆Svib with vanishing ligand mass is a simple approach to obtain the intrinsic ∆Sprot for any ligand, even flexible ones, without the approximations involved when the protein is kept fixed. Therefore, it may be proposed as a standard way to dissect the contributions from protein and ligand to the vibrational entropy of the complex.
Bovine Pancreatic Trypsin Inhibitor Acknowledgment. We thank the BBSRC, U.K., for support. References and Notes (1) Brooks, C. L.; Karplus. M. J. Mol. Biol. 1989, 208, 159-181. (2) Levitt, M.; Park, B. H. Structure 1993, 1, 223-226. (3) Daggett, V.; Levitt, M. Realistic Simulations of Native Protein Dynamics in Solution and Beyond. Annu. ReV. Biophys. Biomol. Struct. 1993, 22, 353-380. (4) Zhang, L.; Herman, J. Proteins: Struct., Funct., Genet. 1996, 24, 433-438. (5) Roux, B.; Nina, M.; Pomes, R.; Smith, J. C. Biophys. J. 1996, 71, 670-681. (6) Denisov, P.; Halle, B. Faraday Discuss. 1996, 103, 227-244. (7) Pettitt, B. M.; Makarov, V. A.; Andrews, B. K. Curr. Opin. Struct. Biol. 1998, 8, 218-221. (8) Prevost, M. Folding Design 1998, 3, 345-351. (9) Helms, V.; Wade, R. C. Proteins: Stuct., Funct., Genet. 1998, 32, 381-396. (10) Garcia, A. E.; Hummer, G. Proteins: Stuct., Funct., Genet. 2000, 38, 261-272. (11) Likic, V. A.; Juranic, N.; Macura, S.; Pendergast, F. G Protein Sci. 2000, 9, 497-504. (12) Wlodaver, A.; Walter, J.; Huber, R.; Sjolin, L. J. Mol. Biol. 1984, 180, 301-329. (13) Otting, G.; Wuthrich, K. J. Am. Chem. Soc. 1989, 111, 18711875. (14) Denisov, V. P.; Peters, J.; Horlein, H. D.; Halle, B. Nat. Struct. Biol. 1996, 3, 505-509.
J. Phys. Chem. B, Vol. 105, No. 33, 2001 8055 (15) Woenckhaus, J.; Hudgins, P. R.; Jarrold, M. F. J. Am. Chem. Soc. 1997, 119, 9586-9587. (16) Fischer, S.; Verma, C. S.; Hubbard, R. E. J. Phys. Chem. B 1998, 102, 1797-1805. (17) Fischer, S, Verma, C. S. Proc. Nat. Acad. Sci. U.S.A. 1999, 96, 9613-9615. (18) Mao, Y.; Ratner, M. A.; Jarrold, M. F. J. Am. Chem. Soc. 2000, 122, 2950-2951 (19) McQuarrie, D. A. Statistical Mechanics; Harper Row: New York 1976. (20) Tidor, B.; Karplus, M. J. Mol. Biol. 1994, 238, 405-414 (21) Denisov, V. P.; Venu, K.; Peters, J.; Horlein, H. D.; Halle, B. J. Phys. Chem. B 1997, 101, 9380-9389. (22) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (23) Neria, E.; Fischer, S.; Karplus, M. J. Chem. Phys. 1996, 105, 19021921. (24) Wilson, E. B.. Jr.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: New York, 1980. (25) Brooks, B. R.; Janezic, D.; Karplus, M. J. Comput. Chem. 1995, 16, 1522-1542. (26) Karplus, M.; Janin, J. Protein Eng. 1999, 12, 185-186. (27) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: Oxford, 1969. (28) Dunitz, J. Science 1994, 264, 670. (29) Schwarzl, S. M.; Tschopp, T.; Smith, J. C.; Fischer, S. J. Med. Chem., submitted.
