J. Phys. Chem. B 2001, 105, 10709-10714
10709
Modeling Zn2+-Cysteinate Complexes in Proteins Todor Dudev† and Carmay Lim*,†,‡ Institute of Biomedical Sciences, Academia Sinica, Taipei 11529, Taiwan, R.O.C., and Department of Chemistry, National Tsing Hua UniVersity, Hsinchu 300, Taiwan, R.O.C. ReceiVed: June 1, 2001
Before computational results can be reliably interpreted, it is critical to calibrate the theoretical calculations with respect to appropriate experimental data. Such calibrations have been carried out for biologically important zinc-cysteinate complexes in this work. Calculations using different quantum mechanical methods were carried out to determine the theory/basis set that can reproduce the X-ray structure of a zinc-cysteinate-like complex and the measured gas-phase deprotonation free energy of H2S. The S-VWN/6-311++G(d,p) method was found to reproduce the X-ray geometry of bis-(ethane-1,2-dithiolato-S,S′)-zinc(II). In particular, it yielded an average Zn-S bond distance of 2.34 Å for bis-(ethane-1,2-dithiolato-S,S′)-zinc(II) and [Zn(CH3S-)4]2-, in excellent agreement with experiment. In contrast, B3-LYP with the same basis set overestimated the average Zn-S bond distance by about 0.1 Å. With the 6-311++G(d,p) basis set, MP2 and post-MP2 methods could reproduce the experimental gas-phase deprotonation free energy of H2S, while DFT methods such as B3LYP and S-VWN yielded less accurate values. Furthermore, a set of effective radii for zinc and atoms of water and HS- consistent with S-VWN/6-311++G(d,p) geometries and NBO charges as well as MP2/6311++G(d,p)//S-VWN/6-311++G(d,p) energies has been obtained. These radii predicted the correct free energy of Zn2+ binding to dianionic 2,3-dimercapto-1-propanol in solution. The results obtained here should help in modeling the structural and thermodynamical properties of zinc-cysteinate binding sites. Moreover, the strategy described in this work could be applied in modeling other metal-binding sites in proteins.
Introduction Due largely to the importance of metals in biological systems, there has been growing interest toward modeling metal-binding sites in proteins. A wide range of problems has been tackled ranging from structure elucidation1-13 to (free) energy evaluations2,7,8,12-21 and reaction profile analyses of various reactions.22,23 Except for a few studies,3,12,13,15,22 most studies do not take the environment (buried or solvent exposed) of the metal-binding site in the protein into account, but employ quantum mechanical calculations to study the metal and its first and sometimes second coordination layer.8,20 Due to the computational expense of quantum mechanical calculations, simple organic molecules have been used to model the amino acid residues most commonly found coordinated to the metal of interest. These are HCOO- or CH3COO- for deprotonated aspartic/glutamic acid side chain; HS- or CH3Sfor ionized cysteine side chain; NH3, imidazole, or methylimidazole for neutral histidine side chain; and formamide for asparagine/glutamine side chain as well as the backbone peptide group. As an example, the ZnCys4 binding site has been modeled by [Zn(HS-)4]2- and/or [Zn(CH3S-)4]2-.3,5,6 The most popular method employed in the quantum mechanical studies of metal complexes6,13,15,18-22 is density functional theory (DFT) in conjunction with the hybrid 3-parameter exchange-correlation functional, B3-LYP,24,25 since it has been found to be cost effective.26 The key limitation in quantum modeling of isolated metal complexes is the neglect of the protein environment, which may * Corresponding author. † Academia Sinica. ‡ National Tsing Hua University.
significantly influence the gas-phase properties derived from ab initio calculations. The influence of the protein on the properties of the metal binding site can be taken into account using combined quantum mechanical and molecular mechanics (QM/MM) calculations. Such an approach has been used to study the coordination number of the catalytic zinc ion in alcohol dehydrogenase.3 However, QM/MM calculations on metalloproteins are currently computationally intensive, thus prohibiting the use of a high theory level with a large basis set to describe the metal binding site. An alternative approach is to treat the protein as a dielectric continuum. This approximation has been used in our previous works12,13,15,27 to compute the reaction free energy in an environment characterized by a dielectric constant ) x, ∆Gx, according to the following thermodynamic cycle: SCHEME 1
Here, ∆G1 is the gas-phase free energy; it can be evaluated using ab initio methods and standard statistical mechanical formulas.28 ∆Gxsolv is the free energy for transferring a molecule in the gas phase to a continuous medium characterized by a dielectric constant, x. Its electrostatic component can be estimated by
10.1021/jp012090f CCC: $20.00 © 2001 American Chemical Society Published on Web 09/29/2001
10710 J. Phys. Chem. B, Vol. 105, No. 43, 2001
Dudev and Lim
solving Poisson’s equation using finite difference methods.29,30 On the other hand, its nonelectrostatic part can be estimated by the solvation free energy of a neutral molecule of similar shape, or it can be implicitly included by adjusting the effective radii of the solute atoms, which define the dielectric boundary in the finite-difference Poisson calculations, to reproduce experimental solvation data. Thus, ∆Gx can be computed from
∆Gx ) ∆G1 + ∆Gxsolv(Products) - ∆Gxsolv(Reactants) (1) There are several advantages in using Scheme 1. First, a high theory level with a large basis set can be used to (a) optimize the geometry of the metal and its ligands and (b) compute the gas-phase free energy, ∆G1. Second, using continuum dielectric methods, it is relatively fast to compute reaction free energies covering the entire spectrum of , which, in turn, represent the reaction occurring in the gas phase ( ) 1), in a buried or partially buried metal-binding site ( ) 2 or 4) and in a fully solvent-exposed metal site ( ) 80). The trends in these free energies can yield useful insights into metal binding and selectivity in proteins, as demonstrated in our previous works.12-15 For example, the free energies for exchanging a metal-bound water with model amino acid ligands in various dielectric media ( ) 1 to 80) showed that a low dielectric medium favors the inner-sphere binding of protein ligands to the metal;12 in the case of magnesium, they delineated its most preferable set of inner-sphere ligands.14 As another example, the free energies of isomerization between hexa- and tetracoordinated Mg2+ and Zn2+ complexes have helped to determine the most preferable ground-state coordination geometry for these complexes in a protein environment.13 Furthermore, the free energies of Mg2+ T Zn2+ exchange in model binding sites have elucidated the factors governing Zn2+ vs Mg2+ selectivity in metalloproteins.15 Before computational results can be reliably interpreted, it is critical to calibrate the theoretical calculations with respect to appropriate experimental data. In this work, we present such calibrations for zinc-cysteinate-like complexes. These complexes have been chosen since cysteine is the ligand most commonly found coordinated to Zn2+ in zinc finger proteins, which include (a) the cellular or transcription factor class, (characterized by a Cys-Cys-His-His metal-binding motif), (b) the retroviral class (Cys-Cys-His-Cys), and (c) the steroid receptor class (Cys-Cys-Cys-Cys).31,32 Furthermore, up to four cysteine residues are also found ligated to Zn2+ in other zinc proteins.33,34 Another reason zinc-cysteinate-like complexes have been chosen is because previous quantum mechanical studies have yielded conflicting results regarding the Zn-S bond distances in model [Zn(CH3S)4]2- and [Zn(HS)4]2systems (see also Discussion).3,5,6 Thus, it is not clear what ab initio theory/basis set can predict the experimental geometry of isolated Zn-cysteinate complexes and to what extent the gas-phase calculations can reproduce crystallographic or solution structural parameters. Furthermore, to the best of our knowledge, there has so far been no attempt to calibrate thermodynamical parameters for reactions between Zn2+ and cysteine-like ligands in the gas phase or in condensed media. Hence, we have carried out a systematic study on the performance of different quantum mechanical methods (theory/ basis set) in reproducing the experimental geometry and gasphase free energy of zinc-cysteinate-like complexes. In computing the solvation free energies in Scheme 1 using continuum dielectric methods, the effective solute radii that determine the dielectric boundary are uncertain. Hence, they were adjusted to reproduce available experimental data for zinc complexes in solution. The calibration data and procedure are described in
Figure 1. Ball and stick diagram of the fully optimized S-VWN/6311++G(d,p) structure of bis-(ethane-1,2-dithiolato-S,S′)-zinc(II), which has C1 symmetry. Bond distances (Å) in italics and bold correspond to the X-ray and S-VWN/6-311++G(d,p) structures, respectively.
the Methods section. The results from the calibration studies should help in modeling the structural and thermodynamical properties of zinc-cysteinate binding sites. For example, they could be employed in (i) mapping the initial steps of zinc-finger folding which is triggered by Zn2+-Cys binding,32 (ii) elucidating the protonation state of zinc-cysteine cores, and (iii) determining the factors that govern metal selectivity in zinccysteinate binding sites. Furthermore, the strategy outlined here is general and could be applied in modeling other metal-binding sites in proteins. Methods Geometry Calibration and Prediction. The Cambridge Structure Database was searched for high-resolution X-ray structures of Zn2+ bound to sulfur-containing ligands that can be used as targets for calibrating the predicted geometries. The compound that most resembled the [Zn(CH3S)4]2- complex of interest was bis-(ethane-1,2-dithiolato-S,S′)-zinc(II) (CSD code DUFYUY, Figure 1), whose structure had been solved to an R factor of 3.6%.35 Furthermore, its average Zn-S bond distance (2.34 Å) is the same as the average Zn-S bond distance evaluated from 22 ZnCys4 binding sites found in the Protein Data Bank.36 Hence, this molecule was fully geometry optimized using various levels of theory and basis sets. The theories tested include Hartree-Fock (HF), second-order Møller-Plesset (MP2), Slater exchange functional with the Vosko-Wilk-Nusiar correlation functional (S-VWN), B3-LYP, B3-P86, and B-VWN. Since addition of polarization functions on both metal and atoms involved in π-bonding to the metal has been shown to significantly affect the geometry of the metal,37 they were included in the basis sets tested, which ranged from 6-31+G(d) to 6-311++G(2d,2p) (see below and Table 1). The theory/basis set that reproduced the X-ray structure of bis-(ethane-1,2dithiolato-S,S′)-zinc(II) was then employed to predict the structure of [Zn(HS)4]2- and [Zn(CH3S)4]2-. Gas-Phase Free Energy Calibration. Since the theory/basis set level that produces accurate geometries may not yield the proper energetics of zinc complexes in the gas phase, calibrations against experimental thermodynamical parameters for gasphase zinc-complex reactions have to be performed. Unfortunately, such experimental information could not be found, but the free energies for deprotonating H2S and CH3SH in the gasphase have been measured.38-40 Compared to the pKa of CH3SH (10.3/14.241,42), the pKa of H2S (6.9/7.043,44) is closer to that of cysteine (8.145/8.446). Hence, H2S may represent cysteine better than CH3SH in modeling reactions that affect the protonation state of cysteine. The gas-phase free energy for deprotonating H2S was computed using eq 5 as described below (Table 2).
Modeling Zn2+-Cysteinate Complexes
J. Phys. Chem. B, Vol. 105, No. 43, 2001 10711
TABLE 1: Calculated and Experimental Structural Parameters for Bis-(ethane-1,2-dithiolato-S,S′)-zinc(II) Complexa method
〈Zn-S〉
〈S-C〉
〈C-C〉
〈S-Zn-S〉b
〈S-Zn-S〉c
〈Zn-S-C〉
HF/6-31+G(d) HF/6-311++G(2d,2p) B-VWN/6-31+G(d) B3-LYP/6-31+G(d) B3-LYP/6-31++G(2d,2p) B3-LYP/6-311++G(2d,2p) B3-P86/6-31+G(d) MP2/6-31+G(d) MP2/6-31++G(d,p) S-VWN/6-31+G(d) S-VWN/6-311++G(d,p) S-VWN/6-311++G(2d,2p) experimentd
2.457 2.460 2.484 2.433 2.433 2.432 2.400 2.356 2.356 2.338 2.342 2.336 2.34(2)
1.831 1.832 1.872 1.848 1.844 1.846 1.833 1.827 1.826 1.821 1.818 1.814 1.81(1)
1.530 1.525 1.550 1.535 1.534 1.529 1.526 1.528 1.527 1.514 1.509 1.508 1.50
91.3 91.1 92.5 92.9 92.8 92.7 93.2 94.4 94.4 94.4 94.2 94.2 94.1(1)
119.3 119.4 118.5 118.3 118.4 118.4 118.1 117.5 117.5 117.4 117.6 117.6 117.6(1)
96.5 96.5 95.8 95.9 95.9 96.0 95.9 95.8 95.7 95.8 95.8 95.8 95.2
a Average values. Bond lengths are in Å and bond angles in degrees. b Bond angle including sulfur atoms from the same ligand. c Bond angle including sulfur atoms from different ligands. d From Rao, 1986.35 Experimental errors, where available, are given in brackets.
TABLE 2: Gas-Phase Free Energies (in kcal/mol) at 298.15 K for H2S f HS- + H+ a 6-31+G(d) S-VWN B3-LYP MP2 MP3 MP4SDTQ CCSD(T) QCISD(T) experiment
6-311++G(d,p)
6-311++G(2d,2p)
6-311++G(3df,3pd)
332.4 338.8 340.1 342.8 343.2 343.6 343.6
333.1 339.6 339.3 341.5 341.4 341.8 341.8
328.7 331.1 335.3 337.7 335.6 344.1 338.0 347.0 338.1 346.9 338.3 347.1 338.2 347.0 345.4 ( 4,38 345.6 ( 2,40 347.1 ( 239
a The free energies in the table were evaluated using S-VWN/6-311++G(d,p) geometries and vibrational frequencies with single-point calculations of the energy at different theory levels (see Methods).
Solution Free Energy Calibrations and Prediction. pKa of H2S. The experimental number in water (6.943 or 7.044) was used to calibrate the effective radii, Reff, for the atoms of HS-. The pKa of H2S in water was computed according to
2.303RT pKa(H2S) ) ∆G1(H2S) - ∆Gsolv80(H2S) + ∆Gsolv80(H+) + ∆Gsolv80(HS-) (2) where R is the gas constant, T is the temperature, ∆G1 is the gas-phase free energy for H2S f HS- + H+, and ∆Gsolv80(X) is the hydration free energy of X. In computing the pKa of H2S, ∆G1 was estimated using eq 5 below with ∆Eelec evaluated at the MP2/6-311++G(d,p)//S-VWN/6-311++G(d,p) level, while the hydration free energies for H2S (-0.7 kcal/mol47) and H+ (-262.5 kcal/mol48) were taken from experiment. The ∆Gsolv80(HS-) was estimated by solving the Poisson equation using finite difference methods (see below) and iteratively adjusting Reff of H and S- until ∆Gsolv80(HS-) yielded a pKa of H2S close to the experimental number. Hydration Free Energy of Zn2+. The experimental value (-484.6 kcal/mol49) was used to calibrate the Reff for zinc and water atoms. The hydration free energy of zinc, ∆Gsolv80(Zn2+), was computed according to
∆Gsolv80(Zn2+) ) ∆G1(ZnW62+) + ∆Gsolv80(ZnW62+) (3) where ∆G1 is the gas-phase free energy for Zn2+ + 6(H2O) f [Zn(H2O)6]2+, and ∆Gsolv80(ZnW62+) is the hydration free energy of [Zn(H2O)6]2+. In eq 3, the first term was estimated using eq 5 with ∆Eelec evaluated at the MP2/6-311++G(d,p)//S-VWN/ 6-311++G(d,p) level, while the second term was estimated by solving the Poisson equation using finite difference methods (see below). The Reff for zinc and water atoms were adjusted to obtain a value of ∆Gsolv80(ZnW62+) that in conjunction with the computed ∆G1 yielded the experimental hydration free energy of Zn2+.
TABLE 3: Effective Solute Radii for Continuum Dielectric Calculations calibration quantity pKa (H2S) ∆Gsolv80(Zn2+) (kcal/mol)
experiment 6.9b 7.0c -484.6d
theory 6.7 -484.7
Reff (Å)a 1.05 (H) 2.29 (S-) 1.40 (Zn) 0.93 (Hwater) 1.67 (Owater)
a
Radius adjusted to reproduce the quantity in the first column for S-VWN/6-311++G(d,p) geometries and NBO charges on the solute atoms. b From Pearson, 1986.43 c From Atkins, 1999.44 d From Burgess, 1978.49
Reaction Free Energy of Zn2+ Binding to DMP2-. The effective radii in Table 3, which had been calibrated to reproduce the pKa of H2S and the hydration free energy of Zn2+, were then tested by predicting the free energy change upon Zn2+ binding to dianionic 2,3-dimercapto-1-propanol (DMP2-). In aqueous solution Zn2+ is experimentally found to be octahedrally coordinated to six water molecules,13,50 and therefore the zincaqua complex was modeled as [Zn(H2O)6]2+. However, upon binding to a negatively charged ligand in a buried or solventexposed site, zinc has been shown to change its coordination number from six to four.13 Hence a tetrahedral geometry was assumed for zinc complexed to a DMP2- ligand, and the binding of Zn2+ to DMP2- was modeled as
[Zn(H2O)6]2+ + DMP2- f [Zn(H2O)2(DMP)]0 + 4H2O (4) The free energy for reaction 4 in water was calculated according to Scheme 1 and eq 1. The gas-phase free energy for reaction 4 was estimated using eq 5 with ∆Eelec evaluated at the MP2/ 6-311++G(d,p)//S-VWN/6-311++G(d,p) level. The ∆Gsolv80(Zn2+) was computed using eq 3 as described above, while the hydration free energies of DMP2- and [Zn(H2O)2(DMP)]0 were
10712 J. Phys. Chem. B, Vol. 105, No. 43, 2001
Dudev and Lim
estimated by solving the Poisson equation using finite difference methods (see below). The solvation free energy of water (-6.3 kcal/mol) was taken from experiment.51 Quantum Mechanical Calculations. Full geometry optimization (without any symmetry constraints) for all complexes, except [Zn(HS)4]2- and [Zn(CH3S)4]2- complexes that have S4 symmetry, was carried out using the Gaussian 98 program.52 The free energy corresponding to a given gas-phase reaction, ∆G1, was computed according to
∆G1 ) ∆Eelec + ∆ZPE + ∆ET - T∆S
(5)
In eq 5, ∆Eelec, ∆ZPE, ∆ET, and ∆S are the differences in the electronic energy, zero-point energy, thermal energy, and entropy between the product(s) and reactant(s), respectively. The geometries of the reactant(s) and product(s) were fully optimized at the S-VWN/6-311++G(d,p) level (see Results), and their vibrational frequencies were evaluated at the same level of theory. No imaginary frequency was found in any of the molecules studied. Using the S-VWN/6-311++G(d,p) geometry of a molecule, its Eelec was estimated by a single-point calculation at a specified theory level (see Results). The zeropoint energy, thermal energy, and entropy corrections were evaluated at 298.15 K using standard statistical mechanical formulas28 with the S-VWN/6-311++G(d,p) frequencies scaled by an empirical factor of 0.9833.53 A counterpoise correction was applied to the gas-phase binding energy for Zn2+ + 6(H2O) f [Zn(H2O)6]2+ to account for basis set superposition error,54 which is likely to be significant in complex formation reactions than in exchange reactions such as eq 4. Continuum Dielectric Calculations. These employed a 71 Å × 71 Å × 71 Å lattice with a grid spacing of 0.25 Å, ab initio geometries and natural bond orbital (NBO) atomic charges.55 The low-dielectric region of the solute was defined as the region inaccessible to contact by a 1.4 Å-radius sphere rolling over the molecular surface. This region was assigned a dielectric constant of two (in ) 2) to account for the electronic polarizability of the solute. The molecular surface was defined by effective solute radii, Reff, which were initially obtained from the CHARMM (version 22)56 van der Waals radii. As mentioned above, the Reff of zinc as well as the atoms of HS- and water were adjusted to reproduce experimental data. Poisson’s equation was solved with out equal to 1 or 80 and in ) 2. The difference between the computed electrostatic potentials in aqueous solution (out ) 80) and in the gas phase (out ) 1) yielded the hydration free energy ∆Gsolv80 of the molecule. Results Geometry Calibration. Table 1 compares selected structural parameters for the bis-(ethane-1,2-dithiolato-S,S′)-zinc(II) complex from full geometry optimization at different levels of theory/basis set with respective X-ray values.35 The results in Table 1 show that the theoretical methods tested can be divided into two groups depending on their ability to reproduce the experimental geometry. Calculations employing HF theory as well as B-VWN, B3-LYP, and B3-P86 functionals significantly overestimate the bond distances in the complex. The effect is especially pronounced for the average Zn-S bond distance, which is overestimated, in some cases (HF and B-VWN calculations) by more than 0.1 Å. The calculated bond angles, especially S-Zn-S, also deviate from the respective experimental values (by 1 to 3°). Increasing the size of the basis set from 6-31+G(d) to 6-311++G(2d,2p) in the HF or B3-LYP calculations does not improve the agreement with experiment.
Figure 2. Ball and stick diagram of fully optimized S-VWN/6311++G(d,p) structures of (A) [Zn(HS-)4]2- and (B) [Zn(CH3S-)4]2-, which have S4 symmetry. The four Zn-S and S-C bond lengths are the same in each structure.
In contrast, calculations employing MP2 and S-VWN57,58 yield structural parameters close to the experimental observables. In particular, the average Zn-S distance differs from the experimental value by only 0.02 Å at the MP2 level, while it matches the experimental value using the S-VWN functional with the 6-311++G(d,p) basis set. The S-VWN/6-311++G(d,p) method also reproduces the individual bond lengths observed in the X-ray structure to within 0.03 Å (Figure 1). For both MP2 and S-VWN calculations, variations in the basis set have little effect on the optimized geometry. Thus for the 6-31+G(d) basis set and larger, including electron correlation is crucial in determining the correct mean Zn-S distance as evidenced by a bond distance of 2.46 Å at the HF/6-31+G(d) level compared to 2.36 Å at the MP2 level using the same basis. Compared to MP2 calculations, the S-VWN calculations with the same basis set are faster and less computationally demanding (especially in evaluating frequencies). Hence, the S-VWN/6311++G(d,p) method was chosen to compute the structures and vibrational frequencies of the other molecules studied in this work. Geometry Prediction. The S-VWN/6-311++G(d,p) method was used to predict the geometries in [Zn(HS)4]2- and [Zn(CH3S)4]2- complexes, which model ZnCys4 cores in steroidreceptor-type zinc fingers and in some enzymes such as alcohol dehydrogenase. The fully optimized structures are shown in Figure 2. The average Zn-S bond distances in [Zn(HS)4]2- and [Zn(CH3S)4]2- are 2.357 and 2.341 Å, respectively. These agree with the average Zn-S bond distance of 2.34 Å evaluated from 22 ZnCys4 binding sites found in the Protein Data Bank.36 Gas-Phase Free Energy Calibration. Comparison between the computed and experimental gas-phase free energies for H2S f HS- + H+ in Table 2 shows that although the S-VWN method can yield accurate geometries of molecules containing sulfur, it fails to produce correct free energies. The calculated ∆G1 values, even with the largest basis set used, significantly underestimate the experimental number (by 12 to 17 kcal/mol). The other DFT method in Table 2, B3-LYP, performs better than S-VWN yielding free energies that are closer to the observed value, but they still fall outside the experimental error bar. The rest of the methods in Table 2 (MP2, third-order Møller-Plesset MP3, fourth-order Møller-Plesset with single, double, triple, and quadruple substitutions MP4SDTQ, coupled cluster with single, double and triple excitations CCSD(T), and quadratic configuration interaction with single, double, and triple substitutions QCISD(T)) produce more accurate free energies than the DFT methods. In particular, these methods in conjunction with the 6-311++G(d,p) basis set yield free energies within
Modeling Zn2+-Cysteinate Complexes
J. Phys. Chem. B, Vol. 105, No. 43, 2001 10713 Discussion
Figure 3. Ball and stick diagram of fully optimized S-VWN/6311++G(d,p) structures of (A) dianionic 2,3-dimercapto-1-propanol (DMP2-) and (B) the 1:1 complex between Zn2+ and DMP2-. Both structures have C1 symmetry.
experimental uncertainty. However, the free energy computed at the MP2/6-311++G(d,p) level (344.1 kcal/mol) is ∼3 kcal/ mol less than that computed with higher-level correlation methods (MP3, MP4SDTQ, CCSD(T), and QCISD(T)) and the same basis, but ∼3 kcal/mol greater than post-MP2/6-311++G(3df,3pd) free energies (Table 2). Since a single-point calculation at the MP3 level on a molecular system with 390 basis functions needs ∼1 GB of computer memory and more than 20 GB of disk space, post-MP2 calculations seem currently impractical for bigger systems. Therefore, MP2/6-311++G(d,p) correction to the reaction energy was employed in computing free energies in condensed media (see below). Solution Free Energy Calibrations. pKa of H2S. Using an effective radius of 1.05 Å for H and 2.29 Å for S-, the hydration free energy of HS- was calculated to be -73.2 kcal/mol. This value was used in eq 2 together with the “MP2/6-311++G(d,p)//S-VWN/6-311++G(d,p)” free energy in Table 2 (344.1 kcal/mol) and experimental hydration free energies for H2S (-0.7 kcal/mol47) and H+ (-262.5 kcal/mol48) to yield a pKa of 6.7 for H2S, which is close to the experimental value (6.9/7.0, Table 3). Hydration Free Energy of Zn2+. The gas-phase formation free energy of [Zn(H2O)6]2+ computed using eq 5 with ∆Eelec evaluated at the MP2/6-311++G(d,p)//S-VWN/6-311++G(d,p) level was found to be -251.9 kcal/mol. Using an effective radius of 1.40 Å for Zn as well as 0.93 and 1.67 Å for the hydrogen and oxygen atoms of water, respectively, the hydration free energy of [Zn(H2O)6]2+ was estimated to be -232.8 kcal/mol. This value together with ∆G1(ZnW62+) in eq 3 yielded a Zn2+ hydration free energy of -484.7 kcal/mol, in accord with the experimental value (-484.6 kcal/mol, Table 3). Solution Free Energy Prediction. The structures of DMP2- and [Zn(H2O)2(DMP)]0 fully optimized at the S-VWN/ 6-311++G(d,p) level are shown in Figures 3A and 3B, respectively. The gas-phase free energy for reaction 4, computed using eq 5 with Eelec estimated by a single-point MP2/6311++G(d,p) calculation, is -387.8 kcal/mol. Using the effective radii in Table 3 to define the dielectric boundary in the continuum dielectric calculations, the hydration free energies of DMP2- and [Zn(H2O)2(DMP)]0 were estimated to be -196.1 and -34.4 kcal/mol, respectively (see Methods). Using the experimental solvation free energy of water (-6.3 kcal/mol) and the computed hydration free energy of [Zn(H2O)6]2+ (see above), the free energy for reaction 4 according to eq 1 is equal to -387.8 + 232.8 + 196.1 - 34.4 - 25.2 ) -18.5 kcal/mol. This number is in agreement with the experimental value of -18.3 kcal/mol, evaluated from the stability constant of the [Zn(DMP)]0 complex (3 × 1013 M-1).59
For zinc hydrates, B3-LYP calculations with various basis sets can reproduce the observed Zn-O distance in solution (2.10 ( 0.07 Å) or in the crystalline state (2.11 ( 0.07 Å) as they yield an average Zn-O distance of 2.12 Å,19 2.12 Å,13 and 2.13 Å18 with the 6-311+G(d), 6-31++G(2d,2p) and 6-311++G(d,p) basis sets, respectively. However, B3-LYP as well as B-VWN, B3-P86, and HF calculations significantly overestimate the average X-ray Zn-S distance (2.34 Å) in bis(ethane-1,2-dithiolato-S,S′)-zinc(II) (Table 1). This is in accord with previous studies, which found an average Zn-S distance in [Zn(HS-)4]2- and [Zn(CH3S-)4]2- complexes of 2.47 to 2.48 Å at the HF/6-31G(d) or HF/6-31+G(d,p) level,3 and 2.54 Å at the B3-LYP/LANL2DZ level.6 The results in Table 1 show that the S-VWN or MP2 method can yield relatively accurate Zn-S bond distances in zinctetrathiolate complexes. In particular, the S-VWN functional together with the 6-311++G(d,p) basis set yields average Zn-S bond distances for bis-(ethane-1,2-dithiolato-S,S′)-zinc(II) and [Zn(CH3S-)4]2- that are in excellent agreement with experiment (see Results). This is in accord with previous work,22 which found an average Zn-S distance of 2.36 to 2.37 Å when [Zn(CH3S-)4]2- was optimized with the VWN58 functional in conjunction with the DZVP2 and DZVP basis sets.60 These results suggest that experimental Zn-S distances can be obtained using the appropriate level of theory (S-VWN) and basis set (6-311++G(d,p)) but without the need to include dielectric effects as for zinc hydrates. This is in contrast to a previous study, which found that “the experimental bond lengths can be reproduced only with methods that account explicitly for the enzyme”.3 Note, however, that the average Zn-S distance in isolated zinc-cysteinate-like complexes in the Cambridge Structure Database, such as bis-(ethane-1,2dithiolato-S,S′)-zinc(II), is the same as that evaluated from 22 ZnCys4 binding sites found in the Protein Data Bank.36 The good agreement between the computed and experimental average Zn-S bond distance in [Zn(HS)4]2- and especially [Zn(CH3S-)4]2- complexes indicates that both HS- and CH3Sare suitable for modeling the geometries of zinc-cysteinate complexes. However, as noted above (see Methods), H2S may be better than CH3SH for modeling reactions that affect the protonation state of cysteine since its pKa (6.9/7.043,44) is closer to that of cysteine (8.145/8.446). There is apparently no unique cysteine model that can reproduce the entire spectrum of its properties. Therefore, for each particular reaction of interest care should be taken in choosing the most appropriate cysteine model. The results in Table 2 show that MP2 and post-MP2 methods perform better than DFT methods such as B3-LYP and S-VWN in reproducing the experimental gas-phase free energy for deprotonating H2S. The quality of the basis set is also important in evaluating correct free energies of the systems studied. Basis sets of triple-ζ quality with diffuse functions on both light and heavy atoms are needed to properly reproduce the experimental free energies of sulfur-containing ligands. The most costeffective method in computing the gas-phase proton affinity of HS- was found to be MP2/6-311++G(d,p)//S-VWN/6311++G(d,p). Effective solute radii used in defining the dielectric boundary in continuum dielectric calculations generally depend on the quantum mechanical method used to (a) optimize the geometry, (b) determine the charges on the solute atoms, and (c) compute the gas-phase free energy. Hence the effective radii listed in Table 3 should be used in conjunction with S-VWN/6-311++G(d,p) geometries and NBO charges as well as MP2/6-311++G-
10714 J. Phys. Chem. B, Vol. 105, No. 43, 2001 (d,p)//S-VWN/6-311++G(d,p) energies. This combination predicted the correct free energy of Zn2+ binding to DMP2- in solution (reaction 4). Note that since the effective radii in Table 3 have been adjusted to reproduce experimental hydration free energies, they implicitly incorporate the nonelectrostatic contributions to the net solvation free energy. We are currently using the combination of methods and parameters found here to study the protonation state of zinc-cysteine cores and the initial steps of zinc-finger folding. Acknowledgment. We thank Dr. J. Wright for technical assistance with the calculations. We are grateful to Drs. D. Bashford, M. Sommer, and M. Karplus for the program to solve the Poisson equation. T.D. is supported by a fellowship from the Institute of Biomedical Sciences. This work is supported by the Institute of Biomedical Sciences at Academia Sinica, the National Center for High Performance Computing, and the National Science Council, Republic of China (NSC-88-2113M-001). References and Notes (1) Ryde, U. Int. J. Quantum Chem. 1994, 52, 1229. (2) Ryde, U. Biophys. J. 1999, 77, 2777-2787. (3) Ryde, U. Eur. Biophys. J. 1996, 24, 213-221. (4) Deerfield, D. W., II; Pedersen, L. G. J. Mol. Struct. (THEOCHEM) 1997, 419, 221-226. (5) Topol, I. A.; Casas-Finet, J. R.; Gussio, R.; Burt, S. K.; Erickson, J. W. J. Mol. Struct. (THEOCHEM) 1998, 423, 13-28. (6) Cini, R. J. Biomol. Struct. Dyn. 1999, 16, 1225-1237. (7) Cini, R.; Musaev, D. G.; Marzilli, L. G.; Morokuma, K. THEOCHEM 1997, 392, 55. (8) Tiraboschi, G.; Roques, B.-P.; Gresh, N. J. Comput. Chem. 1999, 20, 1379-1390. (9) Sabolovic, J.; Liedl, K. R. Inorg. Chem. 1999, 38, 2764-2774. (10) Yliniemela, A.; Uchimaru, T.; Hirose, T.; Baldwin, B. W.; Tanabe, K. THEOCHEM 1996, 369, 9. (11) Shamovsky, I. L.; Ross, G. M.; Riopelle, R. J.; Weaver, D. F. J. Am. Chem. Soc. 1999, 121, 9797. (12) Dudev, T.; Lim, C. J. Phys. Chem. B 2000, 104, 3692-3694. (13) Dudev, T.; Lim, C. J. Am. Chem. Soc. 2000, 122, 11146-11153. (14) Dudev, T.; Cowan, J. A.; Lim, C. J. Am. Chem. Soc. 1999, 121, 7665-7673. (15) Dudev, T.; Lim, C. J. Phys. Chem. B 2001, 105, 4446-4452. (16) Peschke, M.; Blades, A. T.; Kebarle, P. J. Am. Chem. Soc. 2000, 122, 1492-1505. (17) Garmer, D. R.; Gresh, N. J. Am. Chem. Soc. 1994, 116, 35563567. (18) Bock, C. W.; Katz, A. K.; Markham, G. D.; Glusker, J. P. J. Am. Chem. Soc. 1999, 121, 7360-7372. (19) Rulisek, L.; Havlas, Z. J. Phys. Chem. A 1999, 103, 1634-1639. (20) El Yazal, J.; Roe, R. R.; Pang, Y.-P. J. Phys. Chem. B 2000, 104, 6662-6667. (21) Rulisek, L.; Havlas, Z. J. Am. Chem. Soc. 2000, 122, 10428-10439. (22) Topol, I. A.; Nemukhin, A. V.; Chao, M.; Iyer, L. K.; Tawa, G. J.; Burt, S. K. J. Am. Chem. Soc. 2000, 122, 7087-7094. (23) Krauss, M.; Garmer, D. R. J. Am. Chem. Soc. 1991, 113, 642. (24) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (25) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. 1988, B37, 785. (26) Ziegler, T. Can. J. Chem. 1995, 73, 743-761. (27) Dudev, T.; Lim, C. J. Phys. Chem. A 1999, 103, 8093-8100.
Dudev and Lim (28) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (29) Gilson, M. K.; Honig, B. Biopolymers 1986, 25, 2097. (30) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 56105620. (31) Coleman, J. E. Annu. ReV. Biochem. 1992, 61, 897-946. (32) Berg, J. M.; Godwin, H. A. Annu. ReV. Biophys. Biomol. Struct. 1997, 26, 357. (33) Lipscomb, W. N.; Strater, N. Chem. ReV. 1996, 96, 2375-2433. (34) Christianson, D. W. AdV. Protein Chem. 1991, 42, 281-355. (35) Rao, C. P.; Dorfman, J. R.; Holm, R. H. Inorg. Chem. 1986, 25, 428. (36) Lin, P.; Dudev, M.; Dudev, T.; Lim, C., in preparation. (37) Ujaque, G.; Maseras, F.; Lledos, A. Int. J. Quantum Chem. 2000, 77, 544-551. (38) Rosenstock, H. M.; Draxl, K.; Steiner, B. W.; Herron, J. T. J. Phys. Chem. Ref. Data 1977, 6, 736. (39) Bartmess, J. E.; McIver, R. T., Jr. Gas-phase ion chemistry; Academic Press: New York, 1979; Vol. 2. (40) Cumming, J. B.; Kebarle, P. Can. J. Chem. 1978, 56, 1. (41) Howard, P. H.; Meylan, W. M. Handbook of Physical Properties of Organic Chemicals; Lewis Publishers: Boca Raton, 1997. (42) Irving, R. J.; Nelender, L.; Wadso, I. Acta Chem. Scand. 1964, 18, 769. (43) Pearson, R. G. J. Am. Chem. Soc. 1986, 108, 6109-6114. (44) Atkins, P. W. Physical Chemistry; Sixth edition ed.; Oxford University Press: Oxford, 1999. (45) Lide, D. R. Handbook of Chemistry and Physics; CRC Press: Boca Raton, 1998. (46) Dawson, R. M. C.; Elliot, D. C.; Elliot, W. H.; Jones, K. M. Data for Biochemical Research; 3rd ed.; Clarendon Press: Oxford, 1995. (47) Chambers, C. C.; Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 16385-16398. (48) Gilson, M. K.; Honig, B. Proteins: Struct., Func., Genet. 1988, 4, 7-18. (49) Burgess, M. A. Metal ions in solution; Ellis Horwood: Chichester, England, 1978. (50) Marcus, Y. J. Chem. Soc., Faraday Trans. 1991, 87, 2995-2999. (51) Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. J. Solut. Chem. 1981, 10, 563. (52) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, M. C.; Strain, K. N.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, ReVision A.5; Gaussian, Inc.: Pittsburgh, 1998. (53) Wong, M. W. Chem. Phys. Lett. 1996, 256, 391. (54) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (55) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899-926. (56) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (57) Slater, J. C. Quantum theory of molecules and solids. Vol. 4: The self-consistent field for molecules and solids; McGraw-Hill: New York, 1974. (58) Vosko, S. J.; Wilk, L.; Nusair, M. Can. J. Chem. 1980, 58, 1200. (59) Leussing, D. L.; Tischer, T. N. J. Am. Chem. Soc. 1961, 83, 65. (60) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J. Chem. 1992, 70, 560.
