© Copyright 1998 by the American Chemical Society
VOLUME 102, NUMBER 17, APRIL 23, 1998
LETTERS Computation of the Structure-Dependent pKa Shifts in a Polypentapeptide of the Poly[fv(IPGVG), fE(IPGEG)] Family Jorge A. Vila,†,‡ Daniel R. Ripoll,§ Yury N. Vorobjev,| and Harold A. Scheraga*,‡ UniVersidad Nacional de San Luis, Facultad de Ciencias Fı´sico Matema´ ticas y Naturales, Instituto de Matema´ tica Aplicada San Luis, CONICET, Eje´ rcito de los Andes 950-5700 San Luis, Argentina, Baker Laboratory of Chemistry, Cornell UniVersity, Ithaca, New York 14853-1301, Cornell Theory Center, Cornell UniVersity, Ithaca, New York 14853-3801, and Department of Biochemistry and Biophysics, FLOB, The UniVersity of North Carolina, Chapel Hill, North Carolina 27599 ReceiVed: December 18, 1997; In Final Form: February 18, 1998
We report the results of a theoretical study intended to elucidate the molecular basis of the effect of hydrophobicity on the pKa’s of ionizable groups in polypeptides and proteins. In particular, we have focused on a polypentapeptide of a family associated with the process of free energy transduction. Our theoretical calculations, carried out by using a conformational search method in combination with a fast multigrid boundary element method to solve the Poisson equation, led to results that are in good agreement with the experimental observations.
For many years, calculations of the pKa’s of ionizable groups in proteins and their pH titration behavior have been one of the challenges of biophysical chemistry.1-6 This is still one of the most difficult problems in computational chemistry. As a tradeoff between efficiency and accuracy, the majority of the approaches used to tackle this problem consider the solvent and the solute as continuum media with different dielectric constants. The electrostatic potential is then obtained as a solution of the Poisson equation by numerical methods. Two types of numerical methods are commonly used to accomplish this task, namely, the finite-difference method7-9 and the boundary element method.10-13 In particular, a method for solving the Poisson equation using a multigrid boundary element (MBE) method has been developed in our laboratory.12,14-17 This * Corresponding author. † Universidad Nacional de San Luis. ‡ Baker Laboratory of Chemistry, Cornell University. § Cornell Theory Center, Cornell University. | The University of North Carolina.
method has been used, among other applications, to study the helix-coil transition in polylysine.15,16 Also, in combination with conformational sampling techniques based on Monte Carlo methods, the MBE method has been used to study the interdependence of the pKa’s of titratable groups and the conformational preferences of a short oligopeptide.18 Recently, Urry et al.19 reported large pKa shifts of charged glutamic acid groups associated with changes of composition in polypentapeptides (PPP) of the family poly[fv(IPGVG), fE(IPGEG)] (where the mole fractions fv and fE satisfied the relation fv + fE ) 1). The experimentally determined titration curve clearly shows two regimes with respect to fE, viz., fE > 0.65 and fE < 0.65. The first regime is referred to by the authors as the electrostatic regime and the second one as the hydrophobic regime. In particular, they found that the electrostaticallyinduced pKa shift is not as large as the hydrophobically-induced pKa shift in both pure water and aqueous 0.15 N NaCl. This effect is more dramatic if we compare the shifts at very high fE concentration (∼1) versus the shifts at very low fE concentration
S1089-5647(98)00563-X CCC: $15.00 © 1998 American Chemical Society Published on Web 04/04/1998
3066 J. Phys. Chem. B, Vol. 102, No. 17, 1998
Letters
Figure 1. Stereo image of the lowest-energy conformation (E ) -482.6 kcal/mol) of the set displayed as a blue ribbon. Only the side chains of residues GLU4i (red), Ile1i (white), and Pro2i-3 (white) and their solvent-accessible Connolly surface29 are shown. The N-terminus of the polypeptide chain is indicated by the letter N.
(∼0.05) at 37 °C. While the electrostatic effects leading to a pKa shift are very well-documented and appear as a common feature of any theoretical simulation of this type of process, a theoretical treatment of the effects of hydrophobicity on the pKa’s, as observed for this particular peptide, has not received much attention. A recent review by Urry20 presents an analysis of the large amount of experimental data related to the subject and describes a proposed mechanism to account for the phenomenon. It has been proposed that this family of PPPs adopts a folded β-spiral conformation with approximately three pentamers per turn, with β-turns acting as spacers between turns of the spiral.21-23 According to the authors,21,23 this folded structure leads to an arrangement of hydrophobic residues in close proximity to the side chains of Glu4i residues (where the superscript denotes the position of the residue in the pentapeptide and the subscript denotes the ith pentamer in the sequence). These dominant hydrophobic arrangements near the Glu residues have been ascribed24 to the isoleucine side chains of the residues in position 1 on pentamers i and i +1. These structuredependent intramolecular interactions were assumed by Urry et al.20,24 to be hydration-mediated and mainly responsible for the experimentally observed pKa shift of Glu in PPP and related peptides. To test their assumption, we carried out a series of calculations in which the PPP was modeled as a collection of nine pentamers with fE ) 0.11, i.e., with the sequence (IPGVG)4 (IPGEG) (IPGVG)4. This low value of fE enabled us (a) to exclude the often invoked charge-charge interaction effect on the pKa shift and (b) to determine the structural basis of the largest experimentally observed shift for the Glu residue. The evaluation of the conformational energy follows the procedure recently published;17 i.e., the total free energy E(rp, pH) associated with the conformation rp of the molecule in aqueous solution at a given pH can be defined by considering a three-step thermodynamic process (cavity creation, polarization of the solvent, and alteration of the state of proton binding)
involved in transferring the neutral PPP from the gas phase to aqueous solution as
E(rp, pH) ) Eint(rp) + Fcav(rp) + Fsol(rp) + Finz(rp, pH) (1) where Eint(rp) is the internal conformation energy of the molecule in the absence of solvent, assumed15-17 to correspond to the ECEPP/325 energy of the neutral molecule; Fcav(rp) is the free energy associated with the process of cavity creation when transferring the molecule from the gas phase into the aqueous solution; Fsol(rp) is the free energy associated with the polarization of the aqueous solution; and Finz(rp, pH) is the free energy associated with the change in the state of ionization of the ionizable groups due to the transfer of the molecule from the gas phase to the solvent at a fixed pH. The procedure described previously18 was followed for the calculation of each component of the total energy, E(rp, pH). The pKa shift of a titratable group k, due to the transfer from the solvent to the protein environment, is given by
∆pKak ) [E(PS+k) - E(PS°k)] - [E(S+k) E(S°k)]/[ln 10] kBT (2) where E(S+k) and E(S°k ) are the total electrostatic free energies of the kth single isolated residue in the ionized and neutral states, respectively, in the solvent; E(PS+k) and E(PS°k) are the total electrostatic free energies of the kth residue in the ionized and neutral states, respectively, in the protein environment in the solvent; kB is the Boltzmann constant; and T is the absolute temperature. The total electrostatic free energies are given as the sum of Eint(rp), Fsol(rp), and Finz(rp, pH). It should be noted that the effect of ionic strength was not included in the present study because the largest observed pKa shift19 occurs in the absence of salt. As in our previous work,15,18 a set of 256 low-energy structures was used to represent the ensemble of possible conformations for PPP. Eighteen conformations of this set
Letters correspond to β-spirals, while the remaining 238 were coil structures. These conformations were generated as follows: (a) The dihedral angles for the backbone of the β-spiral conformations were obtained from previously proposed models.21,22 An initial conformation, built by using these dihedral angles, was energy-optimized. This energy-minimized conformation was subsequently used as the starting conformation of a global optimization procedure, namely, the EDMC method.26,27 By keeping the set of dihedral angles of the backbone of the initial conformation fixed at values compatible with a β-spiral, a search was carried out to optimize the dihedral angles associated with the side chains of the molecule. From approximately 1000 generated conformations in a single EDMC run, only the 18 lowest-energy conformations were selected. (b) The 238 coil conformations of the set were selected from a series of four independent runs with the EDMC method, in which all the dihedral angles were considered as variables. The four EDMC runs started from initial conformations in which all the dihedral angles, with the exception of the dihedral angles ω, were randomly generated. The dihedral angles ω were initially set to 180° but were taken as variable during the conformational search. During these runs, 4500 conformations were generated and energy-minimized. All the EDMC runs were carried out with the ECEPP/3 force field plus a model of surface hydration28 (SRFOPT). Additionally, it was assumed that the molecule was neutral. The computation of the total free energy and the pKa shift for the Glu residue, for each of the 256 conformations of the set, was accomplished using eqs 1 and 2. It should be mentioned that the number of possible conformations for this sequence is extremely large. Thus, our exploration of the conformational space was limited because the computations are very time-consuming; as an example, a single evaluation of the free energy (eq 1) requires about 15 CPU minutes on a single node of an IBM/SP2 computer. As a result of this analysis, we found that the Boltzmann average of the pKa’s over all 256 conformations of the set led to a value of 5.65, i.e., a pKa shift of 1.35 from the value of 4.3 for glutamic acid in pure water. This result agrees reasonably well with an experimentally observed19 pKa of about 6.0 for this fE concentration. Figure 1 shows the conformation of (IPGVG)4 (IPGEG) (IPGVG)4 with the lowest total free energy of the set. This conformation corresponds to a β-spiral (E ) -482.6 kcal/mol and pKa ) 5.88). The figure shows that the Glu4i residue interacts mainly with Ile1i and Pro2i-3. These results largely support the assumption of Urry et al.,20,24 that the hydrationmediated interaction of residues Glu4i and Ile1i seems to be responsible for the large pKa shift associated with this conformation. The conformation, however, does not show an interaction of Glu4i with residue Ile1i+1, as was proposed by Urry et al. The latter interaction is replaced by an interaction of Glu4i with the Pro2i-3 residue in this computed conformation. In conclusion, our theoretical result for the pKa shift is in good agreement with the experimental value and provides
J. Phys. Chem. B, Vol. 102, No. 17, 1998 3067 support for the proposed molecular basis for the observed pKa shift in this type of polypeptide. Acknowledgment. This research was supported by grants from the National Science Foundation (MCB95-13167), the National Institutes of Health (1 R03 TW00857-01, GM-14312), the NIH National Center for Research Resources (RR-08102), and the “Agencia Nacional de Promocio´n Cientı´fica y Tecnolo´gica-Argentina” (PICT0030). Support was also received from the National Foundation for Cancer Research. The simulations in this work were carried out on an IBM SP2 supercomputer at the Cornell Theory Center (CTC), which is funded in part by NSF, New York State, the IBM Corporation, the NIH National Center for Research Resources (P41RR-04293), and the CTC Corporate Partnership Program. References and Notes (1) Tanford, C.; Kirkwood, J. G. J. Am. Chem. Soc. 1957, 79, 53335339. (2) Zauhar, R. J.; Morgan, R. S. J. Mol. Biol. 1985, 186, 815-820. (3) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219-10225. (4) Yang, A.-S.; Gunner, M. R.; Sampogna, R.; Sharp, K.; Honig, B. Proteins: Struct., Funct., Genet. 1993, 15, 252-265. (5) Honig, B.; Sharp, K.; Yang, A.-S. J. Phys. Chem. 1993, 97, 11011109. (6) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415-436. (7) Warwicker, J.; Watson, H. C. J. Mol. Biol. 1982, 157, 671-679. (8) Sharp, K. A.; Honig, B. Annu. ReV. Biophys. Biophys. Chem. 1982, 19, 301-332. (9) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435-445. (10) Zauhar, R. J.; Morgan, R. S. J. Comput. Chem. 1988, 9, 171-187. (11) Rashin, A. A. J. Phys. Chem. 1990, 94, 1725-1733. (12) Vorobjev, Y. N.; Grant, J. A.; Scheraga, H. A. J. Am. Chem. Soc. 1992, 114, 3189-3196. (13) Rashin, A. A. Prog. Biophys. Molec. Biol. 1993, 60, 73-200. (14) Vorobjev, Y. N.; Scheraga H. A. J. Phys. Chem. 1993, 97, 48554864. (15) Vorobjev, Y. N.; Scheraga, H. A.; Hitz, B.; Honig, B. J. Phys. Chem. 1994, 98, 10940-10948. (16) Vorobjev, Y. N.; Scheraga H. A.; Honig, B. J. Phys. Chem. 1995, 99, 7180-7187. (17) Vorobjev, Y. N.; Scheraga H. A. J. Comput. Chem. 1997, 18, 569583. (18) Ripoll, D. R.; Vorobjev, Y. N.; Liwo, A.; Vila, J. A.; Scheraga, H. A. J. Mol. Biol. 1996, 264, 770-783. (19) Urry, D. W.; Peng, S.; Parker, T. J. Am. Chem. Soc. 1993, 115, 7509-7510. (20) Urry, D. W. J. Phys. Chem. B 1997, 101, 11007-11028. (21) Chang, D. K.; Venkatachalam, C. M.; Prasad K. U.; Urry, D. W. J. Biomol. Struct. Dyn. 1989, 6, 851-858. (22) Wasserman, Z. R.; Salemme F. R. Biolopymers 1990, 29, 16131631. (23) Urry, D. W. Prog. Biophys. Mol. Biol. 1992, 57, 23-57. (24) Urry, D. W.; Peng, S. Q.; Parker, T. M. Biopolymers 1992, 32, 373-379. (25) Ne´methy, G.; Gibson, K. D.; Palmer, K. A.; Yoon, C. N.; Paterlini, G.; Zagari, A.; Rumsey, S.; Scheraga, H. A. J. Phys. Chem. 1992, 96, 64726484. (26) Ripoll, D. R.; Scheraga, H. A. Biopolymers 1988, 27, 1283-1303. (27) Ripoll, D. R.; Scheraga, H. A. J. Protein Chem. 1989, 8, 263287. (28) Vila, J.; Williams, R. L.; Va´squez, M.; Scheraga, H. A. Proteins: Struct., Funct., Genet. 1991, 10, 199-218. (29) Connolly, M. L. Science 1983, 221, 709-713.
