Density Functional Complexation Study of Metal Ions with Cysteine

Dec 3, 2009 - Jose P. Ruelas-Leyva , Ignacio Contreras-Andrade , Juan I. ... Walkimar de M. Carneiro , Stanislav R. Stoyanov , Leonardo Moreira da Cos...
0 downloads 0 Views 1MB Size
Download PDF
466

J. Phys. Chem. A 2010, 114, 466–473

Density Functional Complexation Study of Metal Ions with Cysteine Henna Pesonen,† Reijo Aksela,‡ and Kari Laasonen†,* Department of Chemistry, PO Box 3000, FIN-90014 UniVersity of Oulu, Finland, and Kemira Oyj, Espoo Research Centre, PO Box 44, FIN-02271 Espoo, Finland ReceiVed: June 18, 2009; ReVised Manuscript ReceiVed: October 6, 2009

Metal complexes of cysteine have been studied using a density functional theory based method together with a continuum solvation model. Complexation geometries for metal complexes of cysteine with Zn2+, Mg2+, Ca2+, Fe3+, and Mn2+ have been determined. Different coordination modes have been considered for MCys2 and MCys3 complexes. Complexation energies have been determined, and they have been corrected with (empirical) metal- and ligand-specific parameters, the latter of which was determined separately for Cys(N,O,S) and Cys(N,S) ligands. The results indicate that the preferred binding mode for Zn2+ with cysteine is bidentate (N,S) type binding in tetrahedral or trigonal bipyramidal geometry while Mg2+ and Ca2+ prefer sulfur-free binding sites in octahedral geometry. Fe3+ prefers binding via sulfur and nitrogen atoms, whereas for Mn2+ several equally stable structures were found. The new correction parameters can be applied for other sulfurcontaining ligands to evaluate the binding strength of a new ligand with metal ions. The observed preferences of metal ions in binding are in agreement with the previous knowledge of the behavior of metal ions. Introduction Cysteine is an amino acid that is found in most proteins. The role of cysteine is significant in many biochemical and biophysical applications, e.g., in medicine where it has been shown to be a potential molecule against cancer as it is able to bind carcinogenic acetaldehyde from tobacco smoke.1 A cysteine molecule contains a thiol group (-SH) which makes it unique among the amino acids. Cysteine is easily oxidized to form cystine in which two cysteine monomers are attached to each other by a disulfide bond (-S-S-). The fully protonated cysteine molecule (HSCH2CH[NH2]COOH, (2-amino-3-mercapto)propanoic acid) contains three dissociable protons. Naturally, the most acidic of these is in the carboxylic group with a log K value of 1.71 (recommended value by IUPAC).2 For two other protonation steps the situation is less clear due to the difficulty of assigning the proton-donating groups that are involved in these reactions. The basic equilibrium of cysteine has been discussed for decades, and nowadays a widely accepted fact is that deprotonation occurs simultaneously from NH3+ and SH groups leading to formation of L2- cysteine. pKa for this step lies around 10-11.2 With three different functional groups, cysteine is an interesting ligand from a metal ion selectivity point of view. According to the traditional HSAB (hard and soft acid and bases) principle of Parr and Pearson,3,4 a thiolate group can be classified as a soft base while two other functional groups (COO- and RNH2) are commonly considered hard bases. Thus, the thiolate group has high affinity for the soft metal ions, e.g., Hg2+ and Pd2+, and its binding to hard metal ions, e.g., Mn2+ and Fe3+, is less probable. However, as the cysteine molecule contains also the two harder groups, it is interesting to compare the most stable complex structures for a series of metal ions of different types. Zn2+ can be classified as a metal ion of intermediate hardness, * Corresponding author: tel, +358 8 553 1640; fax, +358 8 553 1603; e-mail, [email protected]. † University of Oulu. ‡ Kemira Oyj.

and it is known to form stable complexes with sulfur donors, e.g., in metalloenzymes. During the last half century, the metal complexes of cysteine have been widely studied, mostly by experimental methods. The equilibrium studies differ from each other at least by method, ionic medium, and temperature. For critical evaluation of the available equilibrium data, the reader is referred to the report of IUPAC.2 Theoretical studies related to cysteine include, e.g., conformational and optical studies.5,6 In studies concerning binding properties,7,8 the cysteine molecule is often reduced to H2S, CH3SH, or CH3CH2SH. Computational ab initio level complexation studies including complete cysteine molecule(s) are only few; e.g., Spezia et al. have studied binding of Co2+ and Co3+ to cysteine with static calculations at the DFT/B3LYP level and with Car-Parrinello molecular dynamics (CPMD).9,10 Recently we have studied metal complexes of various amino polycarboxylic and polycarboxylic acids, including, e.g., EDTA (ethylenediamine tetraacetic acid), DTPA (diethylenetriamine pentaacetic acid), and EDDHA (ethylenediamine-N,N′-bis(2hydroxyphenylacetic acid), using density functional methods and a continuum solvation model (Conductor-like screening model, COSMO).11-14 The structure of the Ca-ISA complex has been further studied with CPMD.11 These ligands are very effective for the complexation of hard transition metals (Fe3+, Mn2+, etc.), but unfortunately, they are not very selective and, therefore, they form stable complexes also with other metals, e.g., with alkaline earth metals Mg2+ and Ca2+. Sometimes lack of selectivity can be a serious problem, for example, in pulp bleaching process using hydrogen peroxide and peracids. Fe3+ and Mn2+ participate in decomposition of the bleaching chemicals,15 and thus, their removal or inactivation by chelation is an important step of the process. Other metals present in the pulp, Ca2+ and Mg2+, do not participate in decomposition of bleaching chemicals but they do form stable complexes with the chelating agents, thus reducing their availability for the chelation of Fe3+ and Mn2+. In order to further develop the method and to extend the group of ligands, we have now chosen cysteine as a model ligand for

10.1021/jp905733d  2010 American Chemical Society Published on Web 12/03/2009

Metal Complexes of Cysteine sulfur-containing ligands. The aim is to study metal-specific preferences for different binding modes. Furthermore, the complexation energies will be determined and they will be compared with the available experimental results. Our goal is not to describe the complexation behavior of cysteine at physiological conditions but, instead, use this unique and wellknown amino acid to further test and develop our method to be more useful for sulfur-containing ligands also. This first cysteine study has a focus on the static calculations and the local coordination geometries with different metals; structural properties of selected complexes in aqueous environment will be further studied with ab initio molecular dynamics. Methods All the geometry optimizations have been carried out using the program Turbomole, version 5.91.16 We have employed the density functional of Becke and Perdew (BP86)17,18 and the polarized valence triple-ζ (TZVP) basis set for all atoms. RI approximation19,20 has been used to speed up the calculations. We have tested three other density functionals (PBE, B3LYP, and PBE0) for selected MCys2/3 complexes, and the results are presented separately after the BP86 results. All the optimizations have been performed using the COSMO21 (Conductor-like screening model) solvent model as it is implemented in Turbomole. Water has been used as the solvent. The optimized radii for H, O, C, N, and S atoms exist, and the COSMO radii for metal ions have been determined with hexaaquo metal complexes in a way that the metal radii were chosen to reproduce the hexaaquo ion solvation energies without the explicit water molecules (See refs 11–13 for the details.) The optimized radii for Mg2+, Ca2+, Zn2+, Fe3+, and Mn2+ are 1.456, 1.722, 1.342, 1.391, and 1.467 Å, respectively. Note that the radii for Ca2+ and Zn2+ have been reoptimized to correspond to the solvation energy obtained with the COSMO within version 5.91 of Turbomole. Transition metal ions have several different spin configurations. In this study we have employed spin unrestricted formalism for the complexes of Mn2+ and Fe3+ and for the closed shell systems, i.e., for alkaline earth metal and Zn2+ complexes, spin restricted formalism. Both Mn2+ and Fe3+ preferred high-spin configurations in their complexes. Results In this section we will first describe the general procedure for the calculations. Then the results for MCys2/3 complexes obtained with the BP86/TZVP + COSMO approach are presented, first geometries and relative energies of the structures with different coordination modes followed by the energetics. To test the solvent description, we have performed some calculations with explicit water molecules. These results are presented after the complexation energies. At the end of the section, reliability of the BP86 functional is evaluated with respect to three other density functionals which were tested for selected complexes. General Procedure. We have decided to ignore the quite complex ionization stability of free cysteine, and thus, we use throughout the study the L2- form of cysteine in which the carboxylic acid and thiol groups are assumed deprotonated while the amino group remains as -NH2. This cysteinato structure corresponds to very alkaline conditions and does not directly represent the physiological conditions of the amino acids but, on the other hand, is consistent with the conditions used in our previous studies.11-14

J. Phys. Chem. A, Vol. 114, No. 1, 2010 467 In order to study metal-specific preferences on a cysteine binding site, we optimized metal complexes with two and three cysteine molecules with different coordination modes. The optimizations were carried out without any symmetry constraints, and identical starting structures were used for all metals. It is well-known that Zn2+ favors metal complexes with low coordination numbers (4-6) while Mg2+ and Mn2+ prefer sixcoordinated complex structures.22 Ca2+ and Fe3+ can adopt also seven- and eight-coordinated structures, but these are excluded in this study. Cysteine has three potential metal-binding groups, and thus, it can act as a monodentate, bidentate, or tridentate ligand. If bidentante binding is assumed, actual binding sites are (1) -NH2 and -COO-, (2) -COO- and -S-, and (3) -NH2 and -S-. With two ligands, we have considered the following bonding patterns and geometries: (1) six-coordinated structure, all possible groups bound to the metal ion (the isomer of type N2O2S2); (2) four-coordinated structure with two equal groups from both ligands bonding to the metal (N2O2, N2S2, and O2S2); (3) four-coordinated structure with the coordination modes N2O1S1, N1O2S1, and N1O1S2. Furthermore, there are three possible geometric isomers for the six-coordinated octahedral complex (1) depending on the arrangement of the equal donor atoms around the metal ion (i.e., cis and trans isomers). Tetrahedral complexes (2) and (3), on the other hand, do not have geometric isomers. The case of three ligands is a bit more complicated since there are in total nine functional groups that can bind to the metal. The structural isomers (all in octahedral geometry) can be divided in groups as follows: (1) isomers of type N3O3, N3S3, and O3S3; (2) isomers of type N2O2S2, and (3) isomers of type N3S2O1, N3O1S2, N2O3S1, N1O3S2, N2O1S3, and N1O2S3. Moreover, all of these isomers include at least two other geometric isomers (cis-trans isomers) depending on the spatial arrangement of the donor atoms in the first coordination sphere. All of these options have been studied, but only the minimum energy configurations for each metal-ligand combination are presented for both cases (with two and three ligands, respectively). Geometries and Relative Energies. Cysteine is known to be an efficient complexing agent for zinc, and mostly due to the biological relevance, ZnCys complexes have been the object for many equilibrium studies.2 Cysteine complexes of four other metals (Mg2+, Ca2+, Fe3+, and Mn2+), on the other hand, are of minor relevance, at least from the physiological point of view. In this section, we will first present in detail the most stable structures for cysteine complexes of zinc and then the structures for other metals at a more general level. Figure 1 shows the minimum energy geometries for [ZnCys2]2and [ZnCys3]4-. Figure 2 presents a structural comparison for the optimized [ZnCys2]2- complex with a crystal structure.23 Metal-ligand bond lengths are shown next to the bonds in both figures. The most stable structure for [ZnCys2]2- complex is a tetrahedral complex in which sulfur and nitrogen atoms from both cysteine ligands are bound to the metal. In comparison to the crystal structure, metal-ligand bonds are slightly longer in the COSMO optimized structure than in the crystal structure. Orientation of the carboxylate groups differs also in a way that in the COSMO optimized structure charged COO- groups point directly to the solvent and they are further away from each other than in the crystal structure. Also with three cysteine ligands Zn2+ binds preferentially via sulfur and nitrogen atoms. The most stable structure is a five-coordinated structure in trigonal bipyramidal geometry with three nitrogen atoms and two sulfur

468

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Pesonen et al. TABLE 1: Relative Stabilities (values in kJ mol-1) of the [MCys2]2-/- Complexesc

Figure 1. The minimum energy structures for [ZnCys2]2- and [ZnCys3]4-. Zn-S and Zn-N bond lengths (in angstroms) are shown next to the bonds. Blue spheres depict nitrogen atoms, red spheres oxygen, and the yellow ones sulfur. Hydrogens are shown in white and carbons in gray.

Figure 2. The COSMO optimized [ZnCys2]2- and a crystal structure23 in comparison. Metal-ligand bond lengths (in angstroms) are shown next to the bonds. Crystal structure is reduced as it is crystallized as Na salt with explicit H2O molecules.

Figure 3. The most stable structures for Mg2+, Ca2+, Mn2+, and Fe3+ with two and three cysteine ligands. In MCys2 complexes, 2a1 is the preferred structure by all four metals. The other structures, 2a2 and 2a3, are within 10 kJ mol-1 from the most stable one, 2a1, for all the metals. With three ligands the structures 3c1 and 3c2 are equally stable for Mg2+ and Ca2+, Fe3+ prefers 3b2, while for Mn2+ the four structures 3b1, 3b2, 3h3, and 3i2 are energetically very close to each other (within 6 kJ mol-1).

atoms bound to the metal. The initial geometry was an octahedron with (N3O1S2) bonding pattern, and thus, detachment of the COO- group occurred during the optimization. This preference of Zn2+ toward tetrahedral and trigonal bipyramidal geometries over the octahedral geometry has been observed earlier in several metalloenzymes.24 The most stable structures for MCys2 and MCys3 complexes of Mg2+, Ca2+, Fe3+, and Mn2+ are shown in Figure 3. These metals seem to have quite different preferences when binding to two or three cysteine ligands than Zn2+ does. With two ligands, all four metals preferred octahedral complex structure 2a1 with the bonding pattern N2O2S2. However, the structures 2a1, 2a2, and 2a3 are energetically very close to each other,

structure

Zn2+

Mg2+

Ca2+

Fe3+

Mn2+

binding

2a1 2a2 2a3 2b 2c 2d 2d + 2H2O 2e 2f 2g

42 10 49 49 0 71

0 6 1 57 54 12 -26 52 45 55

0 5 2 40 26 9 -35 29 44 33

0 5 10 62 67 123 88 114 142 69

0 5 5 27 3 41

N2O2S2 N2O2S2a N 2 O2 S2 O 2 S2 N 2 S2 N 2 O2 N2O2O2wb N 1 O2 S1 N 2 O1 S1 N 1 O1 S2

57 57 28

44 175 20

Initial binding mode N2O2S2, for Zn2+ final binding mode N2S2. 2d structure optimized with two explicit water molecules in trans positions. c The most stable structures are shown in Figures 1, 3, and 4 and the rest in Supporting Information. The absolute complexation energies of the lowest energy structures are given in Table 4. a

b

for Mg2+, Ca2+, and Mn2+ within 6 kJ mol-1. Also for Fe3+, the difference between 2a1 and 2a3 is quite small, only ∼10 kJ mol-1. With three ligands, alkaline earth metals Mg2+ and Ca2+ show a clear preference toward sulfur-free binding site of the cysteine ligand, and thus, for these ions the most stable structures are 3c1 and 3c2. The transition metal Fe3+ prefers bonding to sulfur and nitrogen atoms in geometry 3b2. For Mn2+, on the other hand, we found several structures (3b1, 3b2, 3h3, and 3i2) that are energetically very close (within 6 kJ mol-1) to each other. In all these complexes three nitrogen atoms are bound to the metal, and thus, the structures differ only with regard to the number of metal-binding S- and COO- groups. The preferential binding sites of metal ions with cysteine can be compared more in detail by reviewing the relative stabilities of the complexes with different metal-ligand coordination modes. Table 1 shows the relative stability of the [MCys2]2-/complexes. The initial geometries for the complexes 2a1, 2a2, and 2a3 are octahedral while 2b-2g are four-coordinated complexes started in tetrahedral geometry. Table 1 shows the clear preference of Zn2+ toward sulfurrich binding sites. Sulfur-free structure 2d is more than 70 kJ mol-1 higher in energy than the minimum energy structure 2c with N2S2 binding site. The structure of the 2a2 complex changed during the optimization in a way that both COOdetached from the metal, and thus, the optimized structure has the same bonding pattern as the most stable structure 2c. The energy difference of 10 kJ mol-1 originates from the bent positions of the detached COO- groups in 2a2 structure while in 2c structure these groups point directly out to the solvent. For Mg2+ and Ca2+, the trend is almost the opposite and these metals prefer maximum coordination to nitrogen and oxygen. For all the other metals except Zn2+, octahedral structures 2a1, 2a2, and 2a3 with two tridentate cysteine ligands were the most favorable ones. For Mn2+, also the structure 2c with N2S2 binding mode is within 5 kJ mol-1 from the most stable structure 2a1. Table 2 presents relative stabilities for the [MCys3]4-/3complexes. As the metals may prefer different coordination geometries with a certain binding mode, relative energies presented in the table refer to the minimum energy structure for each binding mode. In general, the trends observed in Table 2 are similar to those seen in the case of MCys2 complexes. The preference of alkaline earth metals Mg2+ and Ca2+ toward a sulfur-free binding site is clear. These metal ions are known to bind overwhelmingly to oxygen rather than nitrogen and sulfur.22 Also Fe3+ has a clear

Metal Complexes of Cysteine

J. Phys. Chem. A, Vol. 114, No. 1, 2010 469

TABLE 2: Relative Stabilities (in kJ mol-1) for the [MCys3]4-/3- Complexesd structure 3a2/3a1 (Zn2+) 3b2a 3c2b 3d2 3e3 3f2/3f3 (Mn2+) 3g1 (Zn2+)/3g2/3g3 (Fe3+) 3h1 (Zn2+)/3h3 3i2 3j2 (Mg2+)/3j3 (Ca2+, Mn2+)/3j4 (Zn2+, Fe3+)

Zn2+ Mg2+ Ca2+ Fe3+ Mn2+

∆ECO ) ∆ECSM ) ECSM(MLn-m) + 6ECSM(H2O) ECSM(M(H2O)6n+) - ECSM(Lm-) (1)

binding

78 40 36 60 69 42 32

116 49 0 74 40 98 83

44 29 0 36 21 47 46

85 0 94 74 92 89 34

87 0 11 56 36 87 43

O 3 S3 N 3 S3 N 3 O3 O 3 S 2N 1 O 3 N2 S1 S 3O 2N 1 S 3N 2O 1

0 31 28

32 17 52

17 14 15

32 60 59

1 6 23

N3S2O1/N3S2c N 3 O2 S1 N 2 O2 S2

a For Zn2+ and Ca2+ the structure 3b1 is within 2 kJ mol-1. b For all five metals, structure 3c1 is within 9 kJ mol-1. c For Zn2+ complex, detachment of COO- occurred during the optimization, and thus, the final binding mode is N3S2. d The energy refers to the lowest energy structure of each binding mode.

preference, and it favors the octahedral structure with three equal cysteine ligands with (N,S) binding each. For Zn2+, on the other hand, the most stable structure is five-coordinated in trigonal bipyramidal geometry with binding via three nitrogens and two sulfur atoms. The three least stable structures for [ZnCys3]2are obtained with the maximum number of metal-binding oxygen atoms. As was seen in the case of MCys2 complexes, the “binding profile” of Mn2+ seems to be somewhere between alkaline earth metals and other transition metals as it has several equally stable structures. This observation is well in agreement, e.g., with the results of Bock et al.25 who, based on the CSD data search and calculations, concluded that Mn2+ has coordination preferences similar to those of Mg2+, but Mn2+ seems to be somewhat softer than Mg2+. For the octahedral complexes presented in Table 2, each binding mode includes at least two geometric isomers. These isomers were often found to have energies that were quite close to the most stable one. For example, the meridional and facial isomers of MCys3 complexes, 3b1 and 3b2, are for all the metals within 18 kJ mol-1. The structures 3c1 and 3c2 are even closer, within 9 kJ mol-1. A cysteine molecule is a relatively small ligand and so there is no steric hindrance between the ligands, at least in MCys2 and MCys3 complexes. Obviously, besides the type and number of metal-coordinating atoms, there are other factors that have an effect on the stability of a certain metal complex. One of these is the size of chelate ring. When cysteine acts as a bidentate ligand and binds via O and S atoms, a six-membered chelate ring is formed. The two other binding modes ((N,O) and (N,S)) lead to formation of five-membered chelate rings. It is worth noting that for all five metals, the MCys3 structure with (O3S3) binding and three sixmembered chelate rings was visibly higher in energy than the most stable structures. For the two smallest cations, Mg2+ and Zn2+, this binding mode was the least stable one while the largest cation, Ca2+, showed relatively flexible binding behavior with cysteine and had all the binding modes within 50 kJ mol-1. The relative energies can be used as a tool for describing metal-specific preferences in binding with cysteine ligand. However, if one wants to study absolute binding strength of the ligand, complexation energies need to be determined. The following section has the focus on complexation energies. Complexation Energies. The complexation energy ∆ECO is defined as

ECSM(X) is the COSMO-corrected total energy of species (X) and M and L refer to metal and ligand, respectively. However, it has been shown previously11 that the absolute complexation energies deviate quite a lot from the experimental free energies of complexation, and thus, further correction to energies is needed. An obvious choice would be to estimate the values of some of the omitted contributions, e.g., thermal effects, and include them in the complexation energies. This was done for selected cases in the study of Sillanpa¨a¨ et al.11 and it did not lead to a significant improvement in accuracy of calculated complexation energies. For this reason, all the neglected contributions have been included in one or two correction parameter(s) and these ligand/metal-specific parameters have been used to correct the raw complexation energies ∆ECO. The corrected complexation energies ∆GC1 and ∆GC2 are defined as

˜ ML ∆GC1(Mi, Lj) ) ∆ECO(Mi, Lj) + G i

(2)

˜M + G ˜L ∆GC2(Mi, Lj) ) ∆ECO(Mi, Lj) + G i j

(3)

Note that the naming has been reversed with respect to the previous studies. Equation 3 includes two correction parameters: ˜ L ) while in ˜ M ) and one for each ligand (G one for each metal (G i j ˜ ML . eq 2 there is only one, metal-specific correction parameter G i These parameters are determined by minimizing the difference between experimental and calculated complexation energies as follows

min

∑ (∆ECO(Mi, Lj) - ∆Gexptl(Mi, Lj) + G˜ML )2 i

i,j

(4) min

∑ (∆ECO(Mi, Lj) - ∆Gexptl(Mi, Lj) + G˜M + G˜L )2 i,j

i

j

(5) In these equations i runs over the metals and j over the ligands. All the parameters have been determined previously for the same method applied to amino polycarboxylic acid ligands, including, e.g., EDTA and DTPA.11 Afterward, the method has been successfully applied to metal complexes of oligomeric ligands as well as to other amino carboxylic acids including, for example, different isomers of EDDHA.12-14 The use of previously determined parameters presumes that new ligands are very similar to those that were used in parametrization, i.e., EDTA like ligands having at least COOH and NH2 groups in their binding site. However, as a ligand, cysteine is clearly smaller than previously studied EDTA type ligands and contains a thiol group which is a completely new donor atom to this kind of study. Therefore, the determination of a new ligand parameter for cysteine is necessary. The choice of a correction parameter for a new ligand should be made carefully. Existing parameters may be used if the new ligand is similar to the previously studied ones from the point of view of its binding mode with the metal of interest. To be able to evaluate the different coordination modes of each metal, one has to know metal-specific preferences in the ligand binding site as well as the preferred geometries. A cysteine ligand can act as monodentate, bidentate, or tridentate ligand and in different combinations of binding atoms. Therefore, to be accurate, correction parameters should be calculated for each

470

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Pesonen et al.

˜ L Parameters (in kJ mol-1) ˜ M, G ˜ ML , and G TABLE 3: G j i i metal (i) 2+

Zn Mg2+ Ca2+ Fe3+ Mn2+ a

˜ Ma G i 19 -4 -18 189 19

˜ ML a G i 116 91 77 284 115

˜ Lj G

ligand (j) b

Cys (N,O,S) Cys (N,S)b EDTA (N,O)a ISA (N,O)a DTPA (N,O)a

80 63 99 89 107

Reference 11. b This work.

type of cysteine ligand separately. On the basis of the relative stabilities presented in the previous section (Tables 1 and 2), preferred binding sites are limited to (N,O,S), (N,O), and (N,S) ˜ L need to be determined. and thus, three values for the G j ˜ L for Cys(N,O,S) can be determined with the experimental G j data of FeCys2 and MnCys2 complexes. It is also possible to ˜ Cys(N,S) but (N,O) type binding in MCys3 lacks determine G experimental equilibrium data. However, (N,O) type binding is found in many previously studied ligands, and thus, the ˜ Lj may be at least suitability of a previously developed G ˜ ˜ ˜ Lj are shown considered. Correction parameters GMi, GMLi and G in Table 3. Note that the binding sites indicated in parentheses for EDTA, ISA, and DTPA refer to types of metal-bonding atoms and do not make any reference to the total coordination number which, at least for DTPA, varies from one metal to another. As seen in Table 3, values for the determined correction parameters are 80 kJ mol-1 for tridentate cysteine ligand and 63 kJ mol-1 for the bidentate (N,S) cysteine ligand. The latter can be used for MCys2 and MCys3 complexes as it has been fitted with the log K values for ZnCys2 and FeCys3.27,29 With ˜ L values of EDTA, ISA, and DTPA, cysteine respect to G j corrections are small. This is expected as the cysteine is clearly a smaller ligand than the three amino carboxylic acid ligands. The table shows that correction depends not only on the absolute size of the ligand but more importantly on the denticity of the same ligand. EDTA is a hexadentante ligand while ISA has only five coordination sites. DTPA is an octadentate ligand but often it adopts six- and seven-coordinate structures, especially with smaller metal ions Mg2+ and Zn2+. If these three ligands are considered, only DTPA is able to form a complex with (N3O3) binding that is favored by alkaline earth metals and, thus, ˜ L of DTPA the most reliable one to be applied we consider the G j for the Mg/CaCys3 complexes. Table 4 presents calculated complexation energies ∆ECO, ∆GC1, and ∆GC2 for the most stable MCys2/3 complexes of Mg2+, Ca2+, Zn2+, Fe3+, and Mn2+. ∆Gexptl is shown where available. The binding site is indicated as well as the explicit value of the used correction. All values are expressed in kJ mol-1.

For MCys2 complexes of Mg2+ and Ca2+, the correction of 80 kJ mol-1 for each cysteine is clearly overestimated as the final ∆GC values are positive (∼40 kJ mol-1). For these metals, a more reasonable choice seems to be the use of ∆GC1 instead of ∆GC2. Unfortunately, also for MnCys2 and FeCys2 the ∆GC2 ˜ Cys parameter was fitted with values are inaccurate even if the G the ∆Gexptl for these two complexes. This indicates that in order to obtain further accuracy in complexation energies, the parameters should be fitted separately for M2+ and M3+ metal ions. However, this would require more available experimental data points and more trivalent metal ions to be included in the set of studied ligands. This would be frustrating, as the main objective was to create and use a simple and straightforward method for estimating the complexation energies with different ligands. ∆GC2 for Mg/CaCys3 complexes were calculated with the previously determined DTPA correction parameter as the DTPA binding site is similar to MCys3 complex with (N,O) binding. Now obtained ∆GC2 values of -80 and -70 kJ mol-1 are slightly higher (less negative) than the corresponding ∆GC1 values. However, both ∆GC values predict relatively efficient binding which was not foreseen based on the literature. Zn2+ prefers (N,S) binding in both complexes, and thus, the calculated (N,S) correction can be applied for both cases. However, as ZnCys3 is a five-coordinated structure, the use of a correction for three ligands can be misleading. As seen in Table 4, the use of a correction for three ligands produces a ∆GC2 value (-81 kJ mol-1) which is clearly too high with respect to the experimental value (-121 kJ mol-1). If the correction is estimated for 2.5 ligands instead, a much more accurate value is obtained (-113 kJ mol-1). Naturally, this approximation is quite rough but justified as the correction depends more on actual denticity than size of the ligand. This means that the new ligand correction parameter can be successfully applied, but it presumes that the complexation geometry is predicted well and taken into consideration beforehand. Also the previously developed parameters can be applied where necessary. For MnCys3 we found several equally stable structures with different coordination modes. The most stable one has (N,S) binding, but other binding modes cannot be excluded. The four most stable structures (Figure 3.) have (N3S3), (N3S2O), and (N3O2S) binding, and thus, also the (N,O) binding site is involved in MnCys3 complexes. As seen in Table 4, the use of a cysteine (N,S) type correction for all the ligands in MnCys3 gives a relatively high ∆GC2 value (-55 kJ mol-1), and thus, ˜ Lj parameter(s) could be more application of some other G appropriate for the case of Mn2+.

TABLE 4: Calculated and Experimental Complexation Energies for MCys2 and MCys3 Complexesa metal 2+

Zn Zn2+ Fe3+ Fe3+ Mn2+ Mn2+ Mg2+ Mg2+ Ca2+ Ca2+

structure

binding

∆ECO

˜L G j

∆GC1

∆GC2

∆Gexptl

2c 3h1 2a1 3b2 2a1 3b 2a1 3c2 2a1 3c2

N 2S 2 N3S2b N 2O 2 S 2 N 3S 3 N 2O 2 S 2 various (see Figure 3) N 2O 2 S 2 N 3O 3 N 2O 2 S 2 N 3O 3

-263 -289 -463 -547 -195 -233 -117 -183 -100 -159

63 63 80 63 80 63 80 107 (DTPA) 80 107 (DTPA)

-147 -173 -179 -263 -80 -118 -26 -92 -23 -82

-118 -81/-113c -114 -169 -16 -55 39 -80 42 -70

-107 -121 -81 -183 -49 N/A N/A N/A N/A N/A

a All values are expressed in kJ mol-1. ∆Gexptl calculated using equation ∆Gexptl ) -RT ln10 log β1 for complexation at 0/0.1/0.15 ionic strength and temperature of 293/298 K, log β1 values taken from refs 26-29. b N3S2 binding. c Calculated for 2.5 ligands as one of the three cysteines acts as a monodentate ligand.

Metal Complexes of Cysteine

J. Phys. Chem. A, Vol. 114, No. 1, 2010 471 TABLE 5: Calculated Complexation Energies ∆ECO (in kJ mol-1) for Selected MCys2 and MCys3 Complexes with BP86, PBE, B3LYP, and PBE0a metal 2+

Zn Zn2+ Fe3+ Fe3+ Mn2+ Mn2+ Mg2+ Mg2+ Ca2+ Ca2+ Figure 4. Optimized structures for ZnCys2 and FeCys2 complexes 2c and 2d with explicit water molecules.

If we look at the ∆GC2 values at a more general level, we can see that they are exclusively higher than ∆GC1 values for the same complexes. Nevertheless, ∆GC1 predicts way too strong binding for MCys2/3 complexes, and thus, the use of a single correction parameter is not recommended. Even though, the observed trend in metal binding (Fe3+ > Zn2+ > Mn2+) is produced correctly also with the ∆GC1 values. Explicit Water Molecules. Most of the studied metals prefer a coordination number of six (or even higher) in their complexes, and thus, the optimization of initially tetrahedral MCys2 complexes 2b-2g often produced a distorted octahedral structure with two vacant coordination sites occupied by the implicit solvent. To further test the solvent description and the stability of structures with low coordination numbers, we added explicit water molecules to MCys2 complexes 2c and 2d. In 2c, water molecules were added in trans positions with respect to the two cysteine ligands. The structure was first optimized for Mg2+ and then for other metals using MgCys2 structure as the starting geometry. For 2d complex, H2O molecules were added in both cis and trans positions. The optimized geometries for ZnCys2 and FeCys2 complexes 2c and 2d with explicit water molecules are shown in Figure 4. For Zn2+ the most important finding was the detachment of at least one of the two water molecules from the metal ion. Initial octahedral structure changed during the optimization into a tetrahedral (2c) or trigonal bipyramidal (2d) geometry. Surprisingly, also MnCys2 with water molecules in trans positions (2d) ended up in a (slightly distorted) trigonal bipyramid. For the complexes of Mg2+, Ca2+, and Fe3+ (more or less) octahedral geometries were conserved. As seen in Figure 4, explicit water molecules have often tilted positions and they form internal hydrogen bonds to the ligand when possible. This preference toward internal hydrogen bonds has been observed in our previous studies, e.g., in [Ca-ISA]2complexes.11 In general, the COSMO solvent model has difficulties in describing the first solvation shell when it contains water (solvent) molecules, i.e., when the ligand does not occupy all coordination sites of the metal ion. Calculated complexation energies ∆ECO are at the same level with the values of 2c and 2d. As the water molecules (one or two) do not seem to directly participate in complexation of Zn2+ with two cysteine ligands, we do not report ∆ECO values for these complexes. This is also the case for trans-[MnCys2d(H2O)2], and for two other Mn2+ complexes, calculated ∆ECO are clearly higher (more than 30 kJ mol-1) than the corresponding complexes 2c and 2d. For Mg2+, Ca2+, and Fe3+, on the other hand, the MCys2d(H2O)2 complex with water molecules in trans positions has more than 30 kJ mol-1 lower ∆ECO than the corresponding 2d complex

structure

BP86

PBE

B3LYP

PBE0

2c 3h1 2a1 3b2 2a1 3b2 2a1 3c2 2a1 3c2

-263 -289 -463 -547 -195 -233 -117 -183 -100 -159

-244 -252 -550 -619 -177 -193 -98 -145 -80 -112

-110 -112 -300 -351 -46 -58 19 -36 31 4

-413 -504 -638 -766 -307 -427 -243 -364 -226 -330

a All values have been calculated in BP86/TZVP/COSMO optimized geometries.

without H2O. Even though, these values are still quite far away from the energies of the most stable structures 3c2 (Mg2+ and Ca2+) and 3b2 (Fe3+). For Mg2+ and Ca2+, this MCys2d(H2O)2 structure is the most stable one when all MCys2 complexes are considered, ∆ECO for those complexes are -144 kJ mol-1 (Mg2+) and -136 kJ mol-1 (Ca2+). Tests on Different Density Functionals. BP86 is a wellknown and widely used but relatively old GGA functional. The method we are using was developed with the BP86 functional, and thus, as long as we are applying the previously developed correction parameters, the use of BP86 is necessary. We have tested three more recent density functionals for selected MCys2 and MCys3 complexes and calculated ∆ECO values both with single-point calculations in BP86/TZVP/COSMO optimized geometry and in geometries reoptimized with the functional in question. The studied functionals are GGA functional PBE30 and hybrid functionals B3LYP31 and PBE0.32 The calculated ∆ECO are shown in Table 5. In general, complexation energies are remarkably functional dependent. Two pure GGA functionals, BP86 and PBE, produce energies that are well at the same level, at least for divalent metal ions. The difference between B3LYP and PBE0, on the other hand, is enormous (>250 kJ mol-1) for all five metals. ∆ECO obtained with B3LYP for ZnCys2, ZnCys3, and MnCys2 complexes are very close to the experimental values (Table 4). For Fe3+, deviation is still notable, and thus, no general conclusions can be made about the performance of B3LYP. The ∆ECO values for the reoptimized complexes are slightly (∼12 kJ mol-1 or less) lower than the ones obtained in BP86/TZVP/ COSMO optimized geometry. From the structural point of view, the shortest metal-ligand bond lengths are obtained with PBE0 while B3LYP produces the longest metal-ligand bonds. It would be quite straightforward to reoptimize the correction parameters for another functional but, unfortunately, it would not eliminate other possible source of inaccuracy which can be at least as important as the chosen functional. These contributions are discussed in the next section. Discussion The results indicate that for all the studied metal ions, bidentate binding with cysteine molecule via sulfur and oxygen atoms is excluded. This is probably due to the repulsion between the charged donor atoms and, on the other hand, due to the larger chelate ring size (six-ring) originating from the (O,S) binding mode. In general, when a cysteine ligand acts as a bidental ligand, the structures with five-membered chelate rings are preferred by all the studied metals. Although in MCys2

472

J. Phys. Chem. A, Vol. 114, No. 1, 2010

complexes all the other metals expect Zn2+ preferred tridental binding with cysteine ligands, the most stable structures are found with three cysteine ligands acting as bidental ligands. An exception for this is ZnCys3 for which one of the three ligands in optimized MCys3 structure was bound to the metal in monodentate fashion via N atom even though the initial geometry was octahedral. In general, calculated complexation energies are at a reasonable level and they produce the observed binding strengths correctly. However, the raw complexation energies ∆ECO do not correspond to the actual experimental free energies or enthalpies. This is because the complexation energy calculations with this approach lack several contributions. Namely, these contributions include at least the inaccuracies of the solvation model,33 nonelectrostatic solvation effects, the ZPE, thermal corrections, corrections that derive from different standard states, internal and translational entropies, and, finally, the BSSE. In the study of Sillanpa¨a¨ et al.,11 several of the neglected contributions were estimated for selected cases, but unfortunately, it did not lead to a significant improvement of ∆GC values. Furthermore, accurate calculation of these contributions is very demanding, and they have (some or all of them) been successfully neglected in other similar studies earlier.34-36 The magnitude of these single contributions can be some dozens of kJ mol-1 but, on the other hand, they are often of different signs and, thus, cancel each other. Furthermore, in order to keep the procedure simple and straightforward, we have successfully used empirical corrections to calculated complexation energies in previous studies and in this present study as well. This is the first time when the method is applied to sulfur-containing ligands, and we are confident about the applicability of the method for completely new ligands for evaluations of their metal binding strength but, obviously, more experimental data are needed to further calibrate the method. New correction parameters were determined for tridentate (N,O,S) and bidentate (N,S) type ligands. As in the ZnCys3 complex the three cysteine ligands were not bound in an equal manner, both correction parameters were fitted with only two experimental data points. The ∆GC2 values calculated with Cys(N,O,S) corrections are not satisfactory and they indicate ˜ L should be determined separately for M2+ and M3+ that G j complexes of a certain ligand. The ∆GC2 values calculated with Cys(N,S) correction, on the other hand, are much more accurate, and thus, the correction can be used to other cysteine-like ligands ˜ ML with (N,S) binding site. Also the metal-specific parameters G i can be easily reoptimized to cover a larger group of ligands if needed. However, treatment of very heterogeneous groups of ligands with only one parameter is probably impossible, and despite the promising results in the previous studies, we are now more skeptical about more extensive applicability of single ˜ ML parameters. G i Ligand selectivity plays an important role in many applications, and thus, fast and uncomplicated methods for evaluating the binding strengths are being sought continuously. It is obvious that actual selectivity in the application, e.g., in the pulp bleaching process with several metal ions and chemicals present in the same system, cannot be studied with this approach or with other relatively simple computational procedures, but naturally, the basic knowledge of the binding behavior of relevant metals is a prerequisite for understanding their behavior in the application. As heavy metals can be harmful in various environments, it would be interesting to extend the study to include more metals, such as lead.

Pesonen et al. Conclusions Metal complexes of a well-known and unique amino acid, cysteine, have been studied using a method based on density functional theory (BP86/TZVP) together with an implicit continuum solvent model (Conductor-like screening model, COSMO). Several different binding modes have been considered, and the most stable structures for MCys2/3 complexes of Zn2+, Mg2+, Ca2+, Mn2+, and Fe3+ are reported. The results show that Zn2+ prefers binding via S- and NH2 groups of cysteine, with two ligands in tetrahedral geometry and with three ligands in trigonal bipyramidal geometry. Mg2+ and Ca2+ favor a sulfur-free binding site with COO- and NH2 as metal-binding groups in octahedral geometry. For Mn2+, we found various equally stable structures with (N3O3), (N3O2S1), and (N3S2O1) binding sites. Fe3+, on the other hand, showed a clear preference toward (N3S3) binding site in MCys3 complexes. In addition to complexation geometries, also the complexation energies have been determined. The previously developed parametrization has now been extended and the ligand-specific correction parameters have been determined for cysteine ligands of different types (denticity). Together with the previously developed metal-specific parameters, they have been used to calculate corrected complexation energies ∆GC2. The new parameters are useful especially for the evaluation of binding strength of transition metal ions with sulfur-containing ligands. We find this method appropriate for evaluating the metal-specific preferences on the cysteine binding site. Effective chelation of heavy metals plays an important role in several applications, and this present approach can be very useful for predicting binding strengths and binding geometries when only a limited amount of experimental data is available. Acknowledgment. CSC (the Finnish IT Centre for Science) is acknowledged for the computational resources. Kemira Oyj is gratefully acknowledged for financial support. Special acknowledgment is given to Dr. Rosa Carceller and Ari Juppo at Kemira Oyj for helpful and stimulating discussions. Supporting Information Available: The optimized structures for all ZnCys2 and ZnCy3 complexes. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Salaspuro, V. J.; Hietala, J. M.; Marvola, M. L.; Salaspuro, M. P. Cancer Epidemiol., Biomarkers PreV. 2006, 15, 146. (2) Berthon, G. Pure Appl. Chem. 1995, 67, 1117. (3) Pearson, R. G. J. Am. Chem. Soc. 1963, 85, 3533. (4) Parr, R. G.; Pearson, R. G. J. Am. Chem. Soc. 1983, 105, 7512. (5) Gronert, S.; O’Hair, R. A. J. J. Am. Chem. Soc. 1995, 117, 2071. (6) Pecul, M. Chem. Phys. Lett. 2006, 418, 1. (7) Rulı`sek, L.; Havias, Z. J. Am. Chem. Soc. 2000, 122, 10428. (8) Rulı`sek, L.; Havias, Z. J. Phys. Chem. A 2002, 106, 3855. (9) Spezia, R.; Tournois, G.; Cartailler, T.; Tortajada, J.; Jeanvoine, Y. J. Phys. Chem. A 2006, 110, 9727. (10) Spezia, R.; Bresson, C.; Den Auwer, C.; Gaigeot, M.-P. J. Phys. Chem. B 2008, 112, 6490. (11) Sillanpa¨a¨, A. J.; Aksela, R.; Laasonen, K. Phys. Chem. Chem. Phys. 2003, 5, 3382. (12) Pesonen, H.; Aksela, R.; Laasonen, K. J. Mol. Struct.: THEOCHEM 2007, 804, 101. (13) Pesonen, H.; Sillanpa¨a¨, A.; Aksela, R.; Laasonen, K. Polymer 2005, 46, 12641. (14) Pesonen, H.; Sillanpa¨a¨, A.; Aksela, R.; Laasonen, K. Polymer 2005, 46, 12653. (15) Yuan, Z.; d’Entremont, M.; Ni, Y.; van Heiningen, A. R. P. Pulp Pap. Can. 1997, 98, T408-T413 (24-29) . (16) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165. (17) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (18) Perdew, J. P. Phys. ReV. B 1986, 33, 8822.

Metal Complexes of Cysteine (19) Eichkorn, K.; Treutler, O.; o¨hm, H.; Ha¨ser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283. (20) Eichkorn, K.; Treutler, O.; o¨hm, H.; Ha¨ser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 242, 652. (21) Klamt, A.; Schu¨u¨rmann, G. J. Chem. Soc. Perkin Trans. 2 1993, 799. (22) Kaufman Katz, A.; Glusker, J. P.; Beebe, S. A.; Bock, C. W. J. Am. Chem. Soc. 1996, 118, 5752. (23) Bell, P.; Sheldrick, W. S. Z. Naturforsch., B: Chem. Sci. 1984, 39, 1732. (24) Progress in Inorganic Biochemistry and Biophysics, Vol. 1 Zinc Enzymes; Bertini, I., Luchinat, C., Maret, W., Zeppezauer, M., Eds.; Birkhauser Verlag: Basel, 1986. (25) Bock, C. W.; Kaufman Katz, A.; Markham, G. D.; Glusker, J. P. J. Am. Chem. Soc. 1999, 121, 7360. (26) Li, N. C.; Manning, R. A. J. Am. Chem. Soc. 1955, 77, 5225. (27) Doornbos, D. A., Thesis, University of Rijks, Groningen, Netherlands, 1965. (28) Zucconi, T. D.; Janauer, G. E.; Donahe, S.; Lewkowicz, C. J. Pharm. Sci. 1979, 68, 426.

J. Phys. Chem. A, Vol. 114, No. 1, 2010 473 (29) Tanaka, N.; Kolthoff, I. M.; Striks, W. J. Am. Chem. Soc. 1955, 77, 1996. (30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77 (18), 3865–3868. (31) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (32) Perdew, J. P.; Ernzerhof, M.; Burke, K. J. Chem. Phys. 1996, 105 (22), 9982–9985. (33) By this we mean that none of continuum solvation models can predict accurately the solvation of very different highly charged species as 3+ Fe ion and the ligand. We see this inherent inaccuracy of the continuum solvation model as the main source of error. (34) Martin, R. L.; Hay, P. J.; Lawrence, R. P. J. Phys. Chem. A 1998, 102, 3565. (35) Pliego, J. R.; Riveros, J. M. J. Phys. Chem. A 2001, 105, 7241. (36) Cosentino, U.; Villa, A.; Pitea, D.; Moro, G.; Barone, V. J. Phys. Chem. B 2000, 104, 8001.

JP905733D