Copper(I)-alkyl sulfide and -cysteine tri-nuclear clusters as models for metallo proteins: a structural density functional analysis

After having set up the computational methodology for Cu(I)-sulfur systems as models for copper proteins, namely using the simple ligands H2S, HS−, CH3SH, and CH3S−, the Cu(I)-Cysteine systems have been investigated: [CuI( S -H2Cys) n ]+ (H2Cys, cysteine, NH2,SH,COOH) [CuI( S -HCys) n ]1− n (NH2,S−,COOH). Finally, the structures for bi-nuclear (Et, CH2CH3), and tri-nuclear [CuI( S -SH)]3, , , , (NH2,SH,COOH), (NH2,S−,COOH, and NH2,SH,COO−), as well as (NH2,S−,COO−), were also optimized to mimic the active center for a metallo-chaperone copper transport protein (CopZ). The X-ray structures for the biomolecules were matched fairly well as regards the Cu–S bond distances and Cu…Cu contact distances in the case the model cysteine S atom is deprotonated. Upon protonation of ligand S atoms, the conformation of clusters is altered and might bring about the di- and tri-nuclear core breakage. These findings suggest that subtle protonation/deprotonation steps, i.e. small and/or local pH changes play a significant role for copper transport processes.

Besides the recent discovery of the role of Cu(I) thioether chemistry in copper transport (Arnesano et al., 2003;Davis & O'Halloran, 2008;Peariso, Huffman, Penner-Hahn, & O'Halloran, 2003), the most common motif able to recognize and bind copper in copper-chaperones requires cysteine residues in CXXC arrangement both in prokaryotic and eukaryotic organisms (Arnesano et al., 2002;Banci, Bertini, & Ciofi-Baffoni, 2009). As a consequence, the chemistry of copper transfer in copperchaperones ultimately depends on the properties of Cu (I)-thiolate bonds that allow the thermodynamic and kinetic control of copper-dependent cellular processes. In addition to the interest related to the role of the metal sulfur (thiolate) bonds in the activity of many metalloenzymes, much effort has been dedicated to explore metallo-chaperone copper transport proteins in solution and at the solid state (Banci, Bertini, Cantini, & Ciofi-Baffoni, 2010;Calderone et al., 2005;Davis & O′Halloran, 2008;Horn & Barrientos, 2008;Huffman & O′Halloran, 2001;Puig & Thiele, 2002;Robinson & Winge, 2010;Rosenzweig, 2001;Rosenzweig & O′Halloran, 2000;Tottey, Harvie, & Robinson, 2005). These studies resulted in mechanistic proposal(s) for copper transfer that appears to happen through affinity gradients and low activation barriers achieved by a network of metal-donor covalent, electrostatic, hydrogen bond interactions in the vicinity of the metal binding site (Banci, Bertini, Cantini, & Ciofi-Baffoni, 2010;Boal & Rosenzweig, 2009;Calderone et al., 2005;Cobine, Pierrel, & Winge, 2006;Culotta et al., 1999;Lutsenko, Barnes, Bartee, & Dmitriev, 2007;Huffman & O'Halloran, 2001). Despite the large amount of experimental data available on Cu(I)metallo-proteins, very few or no theoretical studies are available at present in the literature about the description of Cu(I)-thiolate bonds that can be used to deeply understand the copper trafficking pathways in proteins at our knowledge. Density functional methods, especially the hybrid B3LYP exchange-correlation functional, were proven to be accurate and fast for investigating adequate models for the active sites of metallo-enzymes (Siegbahn, 2006;Siegbahn & Himo, 2009).
We wish to report here on selected results from a density functional structural study relevant to Cu(I) complexes of di-and mono-hydrogen-sulfide, methyl sulfide, and cysteine ligands in order to shed light on the mechanism of copper transfer processes in proteins. We do not intend to analyze the mechanism of copper trafficking itself, but rather to get some understanding on the processes via structural optimizations.

Molecular modeling details
All the density functional theory (DFT) computations were performed by using the Gaussian09 (Frisch et al., 2010) packages implemented on IBM SP5 and SP6 machines at CINECA (Inter-University Computing Center, Casalecchio di Reno, Bologna, Italy). All the structures were fully optimized at the gas phase by using the B3LYP method (Frisch & Frisch, 1998) and the 6-31G ⁄⁄ (Frisch & Frisch, 1998) basis set for all the atoms (BS1). The ligands and most of the complex molecules were also calculated at higher level of theory, B3LYP/6-311+ +G ⁄⁄ (BS2) (Frisch & Frisch, 1998) for comparative purposes. The solvent (water) was taken into account by using Polarizable Continuum Model (PCM) method (Tomasi, Mennucci, & Cammi, 2005) and the structures were fully optimized at BS1 W (BS1 + PCM water ) and often even at BS2 W (BS2 + PCM water ) levels of theory.
All the Cu(I) complex structures were computed at spin multiplicity 1, whereas the Cu(II) complex structures were computed at spin multiplicity 2 (an unpaired electron).
The optimization was carried out up to reaching the convergence criteria implemented on Gaussian03/09 (maximum force, 0.000450 Hartrees/Bohr; RMS Force 0.000300 Hartrees/Bohr; maximum displacement 0.001800 Å; and RMS displacement 0.001200 Å). The hessian analysis was also performed and the structure was accepted as 'fully optimized' at no imaginary frequencies. Molecular drawings were obtained by using GaussView03 software (Dennington et al., 2003) implemented on Pentium IV machines.

H 2 S and HS À
The free ligands H 2 S and HS À were optimized at BS1 and BS2 levels of theory (see Supporting material Table S1 and Figure S1) and showed a very good agreement between the two basis sets. The two free ligands were also optimized taking into account the solvent (water) effect through PCM (Tomasi, Mennucci, & Cammi, 2005) The calculated bond distances and angles at BS1 W level agree very well for protonated and not protonated ligands. In conclusion, the changes of basis set do not alter significantly bond distances and angles.
[Cu I (S-SH 2 ) n ] + and [Cu I (S-SH) n ] 1Àn (n = 1-3) The mono-, bi-, and tri-coordinate models ( Figure  S2(a-c) and Table S2) showed elongations by 0.090 Å (or smaller, BS2) on adding an extra ligand. The bicoordinate species, fully optimized to an almost C2 structure (S-Cu-S, 175.7°, BS2), whereas the structure optimization of a tri-coordinate molecule has a quasi-C 3h symmetry (S-Cu-S, 119.9°(av), BS2). The data reported above for the di-hydrogen-sulfide derivatives are in agreement with those previously published by others (Pavelka, Šimánek, Šponer, & Burda, 2006). The computed Cu-S bond distances for [Cu I (S-SH) n ] 1Àn (n = 1-3) models are reported in Figure S2(d)-(f) and Table S2 (Supporting). By starting from different input structures, the same final quasi-C 2 symmetry was obtained for the bi-(mono-hydrogen-sulfide). The tri-coordinate derivative, by starting from a C s structure, converged to a C 3h symmetry. The computations at BS2 give Cu-S vectors significantly longer than the corresponding values at BS1. These structures were also confirmed when the water solvent effect was considered (both BS1 W and BS2 W).
[Cu I-II (S-SH 2 ) 4 ] +/2+ and [Cu I-II (S-SH) 4 ] 3À/2À The complex models for the tetra-coordinate species were computed both as Cu(I) and Cu(II) metal center. Considering a Cu(I) complex, on starting from a planar arrangement, the computations converge to a H 2 S… [Cu I (S-SH 2 ) 2 ] + …H 2 S system (BS2), whereas starting from a pseudo-tetrahedral arrangement, it fully converges to a pseudo-tetrahedron with Cu-S by 2.327 Å (av, Figure S3; 2.428 Å, av, BS2 W). For comparison purposes, the Copper(II) species was also computed and a pseudo-tetrahedral arrangement was optimized (Figure 1(a)). At BS2 W, the convergence is rather oscillating having Cu-S distance by 2.417 Å (av). Several attempts aimed at refining [Cu I (S-SH) 4 ] 3À (by starting from both pseudo-tetrahedral and planar arrangements) do not converge and the structures collapse, always dissociating two HS groups and originating a [Cu I (S-SH) 2 ]-(C 2h ) system. On the contrary, when starting from a pseudo-tetrahedral [Cu II (S-SH) 4 ] 2À (C 1 ) structure, fully convergence is obtained, with the optimized structures in a compressed tetrahedral arrangement of the donors (pseudo-D 2 , Figure 1 The results show that the usage of BS2 and BS2 W basis sets is recommended for low coordination number complexes and small overall number of atoms. See Table S3 (Supporting) for more geometrical parameters.

CH 3 SH and CH 3 S À
The computed S-C, S-H, and C-S-H bond parameters for CH 3 SH free molecule are well reproduced and the change of basis set and inclusion of solvent effects do not alter significantly the results. Instead significant effects are recorded for CH 3 S À free anion ( Figure S4 and Table S4).
[Cu I (S-S(H)CH 3 ) n ] + and [Cu I (S-SCH 3 ) n ] 1Àn and (n = 1-4) Computed Cu-S distances for protonated and not protonated derivatives (n = 1-3) are reported in Tables 1 and S5  and 2 2À . It has to be noticed that computations for [Cu I (S-SCH 2 CH 3 ) 2 ] À , [Cu I (S-S(H)CH 2 CH 3 )(S-SCH 2 CH 3 )], and [Cu I (S-S(H)CH 2 CH 3 ) 2 ] + used as models for the copper efflux regulator (CueR) in bacteria (Rao, Cui, & Xu, 2010), performed via density functional methods with effective core potentials (ECP) for metal and 6-31 + G ⁄ basis set for main group elements (sometimes improved up to 6-311++G ⁄⁄ ), produced structural arrangements similar to those from the present work and even bond distances agree well. A certain degree of elongation by 0.03-0.06 Å is recorded at ECP level.

Results and discussion
After the computational methodology was constructed on the basis of what is reported in molecular modeling details section, the aminoacid cysteine was selected for further calculations in order to study the effect of ligand size expansion and with the aim to include biologically significant ligands (see just below). Selected free ligand and copper(I) and copper(II) (for a few cases) complex Scheme 1. Structural formula for (a) L-H 2 Cys(NH 2 ,SH, COOH); (b) L-HCys À (NH 2 ,S À ,COOH); (c) L-HCys À (NH 2 ,SH, COO À ); and (d) L-Cys 2À (NH 2 ,S À ,COO À ). Structures L-H 2 Cys, ÀHCys À (S À or COO À ), and -Cys 2 -(S À , COO À ) Cysteine is a very flexible molecule and therefore, several possible orientations can be assumed by the backbone. Some computational studies on the selected possible conformations have been reported by others (Dobrowolski, Rode, & Sadlej, 2007;Fernández-Ramos et al., 2000). Here, we wish to present and discuss geometrical parameters for a particular conformation of H 2 Cys (NH 2 ,SH,COOH) computed by starting from a zwitterionic structure (NH þ 3 ,SH,COO À ) (Figure 3(a)) and fully optimized at gas phase, with a S-H…O(=C) intra-molecular interaction (Figure 3(b); O…S, 3.278 Å and S-H…O, 12.63°(BS2), Tables 3 and S7) and those for the corresponding conformation for HCys À (NH 2 ,S À ,COOH) ( Figure 3(c), BS2), where a S…H-O interaction takes place (O…S, 2.936 Å and O-H…S, 159.7°, BS2). Both models presented also an intra-molecular N…O hydrogen bond interaction. The N atom acts as H-acceptor and as H-donor in cysteine and cysteinate, respectively (N…O, 2.610 and 2.691 Å and Ĥ, 123.5 and 110.4°, BS2). A Cys 2À (NH 2 ,S À , COO À ) model similar to that found via X-ray diffraction for the zwitterionic form (input, Figure 3(d)) after removal of SH proton and a proton from the amine function was optimized (see Figure 3(e) for BS2). The geometrical parameters for the ligand are in good agreement with the corresponding ones found for the protonated species. Optimized S-C and (S)C-C bond distances were 1.877 and 1.531 Å (BS2) and 1.861 and 1.539 Å (BS2 W). Calculated geometrical parameters for H 2 Cys(NH 2 ,SH,COOH) model are also in good agreement with the X-ray structure of the zwitterionic form (Tables 3 and S7) (Moggach, Clark, & Parsons, 2005). It is interesting to note that the minimum energy optimized structure (gas phase) for H 2 Cys is obtained by starting from the solid state zwitterionic form (Moggach, Clark, & Parsons, 2005; Figure 3(f)) and that the output is not zwitterionic. On the contrary, upon optimizing the zwitterionic structure from Figure 3(a) and on adding the treatment of solvent effect, the fully optimized structures maintain zwitterionic form (Figure 3(g)). Computed bond parameters are in excellent agreement with the experimental values. Other findings indicate that the final structures were significantly dependent on the starting structures as regards the orientation of thiol, amine, and carboxyl groups and were much dependent on the presence of solvent and on protonation modes. The Cu(I) complex species went to full convergence (at BS1/W, BS2/W levels), for mono-, bi-, and tri-coordinate models with both H 2 Cys(NH 2 ,SH,COOH) and HCys À (NH 2 ,S À ,COOH) ligands. Optimized structure of [Cu I (S-H 2 Cys)] + at BS2 level has Cu-S and S-C bond distances by 2.231 and 1.856 Å (BS2 W; see Tables 4 and S8). Comments on other Cu(I)-H2Cys complexes are reported below in this section ( Figure 4). The results for HCys (NH 2 ,S À ,COOH) derivatives are in agreement with the geometrical parameters reported in section 'Molecular modeling details' for small ligands, in the case of monodentate model, and show a linear and a trigonal arrangement for [Cu I (S-HCys) 2 ] À and [Cu I (S-HCys) 3 ] 2À (Figure 5 (a) and (b)), where S-Cu-S bond angles are 179.9 and 120.0°(av, BS2), respectively (see Tables 5 and S9). The calculated Cu-S bond distances are 2.160, 2.202, and 2.331 Å (av, BS2) for the mono-, bi-, and tri-coordinate complexes, respectively. Even for the metal species, O-H…S hydrogen bonds play significant roles in stabilizing the structures. It is important to note that the mono-coordinate species has been fully optimized both as mono-dentate  Scheme 2. Molecular drawing for copper(I) dimer (a) trimer (b) models from the metallo-chaperone proteins (Singleton et al., 2009;Hearnshaw et al., 2009). Sulfur atoms might be protonated or not protonated. model is stabilized by a strong N…O intermolecular interaction (N…O, 2.528, and 2.498 Å at BS2 and BS2 W, respectively). Going back to the Cu I -H 2 Cys derivatives, on starting from the structure shown in Figure 4(a) the fully optimized model at BS2 level converges to the molecule depicted in Figure 4(b) (Cu-S, 2.241 Å; S-Cu-S, 176.4°, av). The structure is stabilized by O…S and O…N intermolecular interactions (O…S, 3.171 Å, and N…O, 2.665 Å). The solvent effects do not have significant importance at least in the coordination sphere region. Noteworthy, the optimized structure (BS2) for tri-cysteine derivative shows a trigonal arrangement (Figure 4(c)) similar to the other tri-coordinate computed species with smaller ligands (see above). Computed (BS2) Cu-S, S-C, and (S)C-C bond distances are 2.344, 1.860, and 1.539 Å.
The tetra-coordinate H 2 Cys derivative was calculated up to BS2/BS2 W levels and the structures converged to pseudo-tetrahedral arrangements (Figure 6(a)). The Cu-S bond length is 2.426 Å (av, BS2 W). The structure is stabilized by several hydrogen bonds interactions between different cysteine residues.
The computed [Cu I (S-HCys) 4 ] 3À model level shows a certain elongation of a HCys À residue: Cu-S bond dis- (c) computed structure for anionic HCys À (NH 2 ,S À ,COOH) BS2; (d) starting structure for Cys 2À (NH 2 ,S À ,COO À ) from X-ray data (Moggach, Clark, & Parsons, 2005) of the zwitterionic H 2 Cys(NH þ 3 ,SH,COO À ) after removal of SH proton and a NH þ 3 proton; (e) computed structure for anionic Cys 2À (NH 2 ,S À ,COO À ) BS2 by starting from (d); (f) structure for zwitterionic H 2 Cys(NH þ 3 ,SH,COO À ) from X-ray diffraction (Moggach, Clark, & Parsons, 2005); and (g) computed structure for zwitterionic H 2 Cys(NH þ 3 ,SH,COO À ) BS1 W by starting from (a).    Figure 6(b) and 2.464 Å av, BS2 W). Conclusion: Cu(I)-cysteine complexes were investigated in order to simulate the behavior of the aminoacid residue as ligand in enzymes and proteins as regards the coordination properties of SH and S À functions. Of course, the models fail to reproduce exactly the conformations around the amides functions through NH 2 and COOH groups. As regards the bi-, tri-, and tetra-H 2 Cys and HCys À derivatives, the easiest models to manage were those for the tri-derivatives. In all cases, full convergence is quickly reached and the CuS 3 core has almost C 3h symmetry. The equivalence of the coordination bonds does not occur for the tetra-coordinate complexes. This finding suggests that the tri-cysteine models are the most reliable ones and that Cu(I) centers encountered in biological molecules tend to reach that type of arrangement. CopZ Dimer and Trimer Models, ½Cu I (S À S(H) Et) 2 2þ 2 , ½Cu I (S À SEt) 2 2À 2 , [Cu I (S-SH)] 3 , ½Cu I (s ÀSH) 2 3À 3 , ½Cu I (S À S(H)Et) 2 3þ 3 , ½Cu I (S À SEt) 2 3À 3 , ½Cu I (S À H 2 Cys) 2 3þ 3 , ½Cu I (S À HCys) 2 3À 3 . Recently, different crystal structures of the copperchaperone CopZ (Protein Data Bank codes, 2QIF and 3I9Z) were reported Singleton et al., 2009). In the solid state, the protein shows a trinuclear Cu(I) 3 (S-Cys) 6 cluster arising from the interaction of three CopZ molecules where the Cys residues come from the CXXC motif (Scheme 2(b); Singleton et al., 2009;3I9Z) The three Cu(I) ions are bound to three bridging (b) S-Cys atoms and to three S-Cys terminal (peripheral, p) sulfur atoms in a trigonal geometry. In a different crystal form of the same protein a tetranuclear Cu(I) cluster is observed (Protein Data Bank code, 2QIF . In this latter case, the tetra-nuclear entity consists of an inner bi-nuclear Cu (I) core bound to two S-Cys (b) from each subunit (CXXC motif) and of two additional Cu(I) atoms in a 2 S-Cys, N-His O-water four-coordination. Here, two Cu(I) ions (b) display distorted diagonal coordination and two Cu(I) ions (p) display tetrahedral coordination. Earlier structure determinations on Cu(I)-CopZ were performed by X-ray absorption spectroscopy (XAS) in solution . In diluted solution Cu(I)-CopZ dimeric forms of the protein were observed only in high ionic strength where the dinuclear Cu(I) 2 (S-Cys) 4 clusters with trigonal coordination are formed. In the trimeric form the overall number of cysteine ligands is six, while it is four in the CopZ dimer (see Scheme 2). Several models for the dimeric and trimeric CopZ were optimized (Figures 7(a) and (b), 8(a)-(d), 9(a) and (b), 10(a) and (b)) in order to shed light on the structure of the cores. The experimental Cu…Cu contact distances are 3.12 Å (estimated standard deviation, esd, c. 0.14 Å) for the trimer and 2.74 Å (esd c.0.1) for the dimer. The optimized structure (BS2 W, 22nd cycle of refinement) for the di-nuclear model ½Cu I (S À S(H)Et) 2 2þ 2 , (Et, CH 2 CH 3 ) has Cu…Cu contact distance by 3.591 Å, and Cu-S(b) and Cu-S(p) bond distances by 2.348-2.547 and 2.268 Å (Tables 6 and S10, see also Figure 7(a)). The fully optimized (BS2 W) structure for the di-nuclear model ½Cu I (S À SEt) 2 2À 2 , (Figure 7(b)) has Cu…Cu contact distance by 2.717 Å and Cu-S(b) and Cu-S(p) bond distances by 2.377-2.405 and 2.254 Å (av). The computed intermetallic distance is in excellent agreement with the ones found by XAS  and even the computed Cu-S distances agree well with the found ones. Even when computations are performed in the gas phase (BS2) the cluster arrangement goes to full convergence. The calculated Cu…Cu, Cu-S(b), and Cu- S(p) bond distances being 2.912, 2.415, and 2.281 Å (av), respectively. Conclusion: the real cluster is well simulated by the thiolate model (BS2 W). The protonation of the ligand at sulfur reasonably brings about instability and easy dissociation of the cluster.
The structure for the basic tri-nuclear model [Cu I (S-SH)] 3 was fully optimized at BS2 and BS2 W (Figure 8 (a) and Tables 7and S11) and Cu…Cu and Cu-S distances are 2.700 and 2.248 Å and 2.693 and 2.252 Å, respectively, in fairly good agreement with experiment (Singleton et al., 2009). The basic six-membered cycle with three peripheral ligands model ½Cu I (S-SH) 2 3À 3 was fully optimized via BS2 level, to an arrangement having long Cu…Cu distances (4.073 Å) (Figure 8(b) and Tables 7 and S11). The S-Cu-S bond angles (BS2) were 109.8°( endo-cycle) and 123.1 and 126.8°(exo-cycle) in good agreement with the experimental ones (108.3°, 124.6°, and 127.0°, respectively) (Singleton et al., 2009). The computed cyclic structure is much less puckered than the experimental one. In case the water solvent treatment was taken into account (BS2 W), the convergence of the model is full and the computed Cu…Cu bond distance is 3.493 Å, in acceptable agreement with experiment. Also Cu-S-Cu-S torsion angles (72°(av)) are close to the experimental ones.
The ½Cu I (S-SH 2 ) 2 3þ 3 model was fully optimized at BS1 level (Figure 8(c)), but when computed at BS2 level undergoes complete dissociation of the cluster after 34 cycles. A [Cu I (S-SH) b (S-SH 2 ) p ] 3 model was also calculated at BS2 level giving a partially refined structure (after c. 90 cycles), where the Cu…Cu, Cu-S(b), and Cu-S(p) distances are 2.756, 2.265, and 3.063 Å (av), respectively (Figure 8(d)). Interestingly, this analysis shows that on increasing the protonation at sulfur starting from ½Cu I (S À SH) 2 3À 3 , the core arrangement relaxes.
Then the models were upgraded to larger ones, but in some cases the computations were performed at the BS1 W level instead of BS2 W, in order to reduce time of convergence (computational costs).
The computed structure for ½Cu I (S À H 2 Cys) 2 3þ 3 , (NH 2 ,SH,COOH) (BS1 W) converges to a relaxed trimer (Figure 10(a)) in which the Cu…Cu, and Cu-S(b) and Cu-S(p) distances are 3.573, and 2.221 and 2.234 Å (av, see Tables 8 and S12), respectively, in good agreement with experimental data for CopZ protein (it has to be noted that optimization was halted after 490 cycles, when the energy trend was flat for the last 160 cycles and structural parameters did not change appreciably). The simulation for ½Cu I (S-HCys) 2 3À 3 (NH 2 ,S À ,COOH) (BS1 W) successfully converges (Figure 10(b)) and the trimer has Cu…Cu and Cu-S(b) distances of 2.429 and 2.289 Å, respectively.
In order to reach a better simulation for the trinuclear core in the protein, the cysteine molecules were replaced by di-amidate forms (HS-CH 2 -C(H)-(C(=O)-N (H)-CH 3 )-N(H)-C(=O)-CH 3 ) and structurally optimized at BS1 level (preliminary computations). Interestingly, the main characteristics of the structure agree very well with the not amidate H 2 Cys derivative. For example, selected distances range: Cu…Cu, 3.806-4.030 Å; Cu-S (b), 2.233-2.338 Å; and Cu-S(p), 2.284-2.317 Å. The results suggest that the not amidate models represent fairly well the structural properties of the metal cluster of the protein active site. Conclusion: this section shows that deprotonation at all sulfur atoms causes a significant shrinking of the core of trimer. This fact should correspond to a significant stabilization of the core; on the contrary, stepwise protonation at sulfur brings about a swelling at nucleus and a trend toward dissociation.
Gross energy aspects. The values of total electronic energies for selected computed models are reported in Table S13 (Supporting) and selected electronic formation and protonation energies are reported below under the respective sections.
[Cu I (S-SH 2 ) n ] + (n = 1-4) The computed electronic energies of formation are listed in Table 9. On adding a second H 2 S ligand to Cu(I) the bond formation energy is still high (À26.38 kcal/mol, BS2 W). The effect of the BSSE (basis set superposition error) correction is small at BS2 level (c. 2%). The values for the mono-coordinate derivative compare well with corresponding value computed previously from this laboratory for Zn(II) derivatives (c. À88 kcal/mol, B3LYP/Lanl2dz level, gas phase, no BSSE correction) (Cini, 1999). The larger bond energy for the Zn(II) derivative is reasonably due to the larger charge at metal. Going back to Cu(I), the formation energies for the fourth Cu-S bond approach zero (4.08 kcal/mol, BS2 W). [Cu I (S-SH) n ] 1Àn (n = 1-4) and [Cu II (S-SH) n ] 1Àn (n = 3,4) The electronic energies of formation are listed in Table 10. Once the BSSE correction is performed, the overall formation energy for the tri-coordinate species (À200.13 kcal/mol, BS2/BSSE) changes by a mere 0.1%. The addition of the third hydrogen sulfide anion to the bi-coordinate derivative is slightly exo-energetic (À7.42 kcal/mol) at BS1 W. It should be recalled that the structure for the tetra-coordinate species (overall charge 3À) did not converge (BS1, BS2). The addition of the fourth HS À ligand is endo-energetic by + 8.66 kcal/mol. Therefore, the formation of tetra-sulfide is not favored from the overall electronic energy stand-point both in the gas phase and in aqueous sys-   tems. In contrast to Cu(I)-sulfide derivatives, Cu(II)sulfide tetra-coordinated species have high overall formation energies: À232.53 kcal/mol (BS1 W), À666.82 (BS2), and À230.47 (BS2 W). Notwithstanding, these data do not mean that the step formation for tetracoordinate Cu(II) species is exo-energetic. In fact, the  overall formation energy for the trigonal [Cu II (S-SH) 3 ] À is À222.75 kcal/mol (BS2 W). Hence, the formation of the tetra-coordinate species from tri-coordinate is exo-energetic by just À7.72 kcal/mol (BS2 W). Conclusion: from this section, it may be said that formation of four-coordinate M(S-L) 4 species is not favorite with small sized mono-valent (and even di-valent) metal ions. Tri-coordination seems to be a sort of upper-limit.
[Cu I (S-S(H)CH 3 ) n ] + and [Cu I (S-SCH 3 ) n ] 1Àn (n = 1-4) The overall electronic formation energies ( show that the formation of up to tri-coordinate species both for the neutral H 2 S and CH 3 SH ligands, as well as for the anionic forms is predictable because significantly stabilized in terms of electronic energy (once the solvent effect is considered). The formation of the tetra-coordinate species is uncertain because the electronic energy contribution approaches to zero or is even positive.
[Cu I (S-H 2 Cys) n ] + and [Cu I (S-HCys) n ] 1Àn (n = 1-4) The overall electronic formation energies are listed in  mol (no BSSE). A similar pattern is found also when the mono-and bi-coordinated models are investigated for the HCys À (NH 2 ,S À ,COOH) ligand.
The overall formation energy for [Cu I (S-SH)] 3 (see Table 13) computed at BS2 and BS2 W are À701.35 and À246.68 kcal/mol, whereas for computations at BS2/ BSSE the value is À697.44 kcal/mol. When a HS À ligand (peripheral) is added to each Cu(I) center the overall formation energy for ½Cu I (S À SH) 2 3À 3 is À627.92 and À622.30 kcal/mol at BS2 and BS2/BSSE, respectively. One can state that: (i) the solvent effect is very effective, (ii) computations at BS2 and BS2/BSSE levels reveal that the addition of the three exo-cyclic HS À functions is endo-energetic by c. + 75 kcal/mol, and (iii) the basis set superposition effect is not large for such cluster molecules.
The models have been then expanded to EtS À and EtSH ligands, and even to cysteine/cysteinate ligands. For these latter systems, the computations were carried out up to BS1 W only for reducing computational times.
Taking into account EtS À ligand the overall formation energy for ½Cu I (S À SEt) 2 3À 3 is À594.88 kcal/mol (BS1/BSSE), whereas the corresponding value at BS1 W is À482.50 kcal/mol. The effect of solvent reduces the overall formation energy by 70%, in agreement with the corresponding solvent effect for the HS À derivative (76%).As regards EtSH, the overall formation energy for ½Cu I (S À S(H)Et) 2 3þ 3 is À133.38 kcal/mol (BS1/BSSE) and À338.19 kcal/mol (BS1 W, partially optimized). The overall formation energy per Cu(I) center for the EtSH tri-nuclear derivative is therefore À112.73 kcal/mol (BS1 W), whereas the corresponding value for the EtS À species is À160.83 kcal/mol (BS1 W), these latter values being some 6-9 kcal/mol more exo-energetic than the corresponding value for the di-nuclear species. These data are in agreement with lower strain tensions for sixmembered cycles with respect to four-membered cycles.
The overall formation energy for ½Cu I (SÀ H 2 Cys) 2 3þ 3 is À163.09 (BS1/BSSE) and À338.76 kcal/ mol (BS1 W), with relative small effect by the solvent (c. 6%) but large effects from BSSE. As regards ½Cu I (S À HCys) 2 3À 3 , the overall formation energy is À571.50 (BS1/BSSE) and À466.04 kcal/mol (BS1 W), much higher values than for neutral cysteine. It is interesting to note that the protonations for the six ligated cysteinate(À1) correspond to an endo energetic process by c. + 127 kcal/mol (BS1 W, mean c. + 21 kcal/mol per sulfur atom), in agreement with corresponding values found before for mono-nuclear species (Section 3.2.4). Conclusion: we propose that the proton transfer to XS À (X = H,C) linked to Cu(I) can be easily reached as a consequence of low energy transfer processes like those involved in hydrogen bond formations. Once the protein contains di-or tri-nuclear clusters with anionic thio-ligands, the clusters themselves are compact from a structural point of view (see above, structural aspects) and stable enough to guarantee the transport of copper. In case pH decreases or conformational changes at protein chains facilitate the formation of hydrogen bonds, one or more S-cysteine involved in metal clusters can be protonated. This latter events swell the clusters (the protonation status can be mixed, cluster dimensions can vary in a wide range) and make dissociation easier, thus contributing to the copper trafficking mechanism.

General conclusion
The work shows that the existence of bi-nuclear ½Cu I (S À SX) 2 2À 2 and tri-nuclear ½Cu I (S À SX) 2 3À 3 species for important biological species related to copper trafficking can be predicted and rationalized through structure optimization DFT methods. Furthermore, the work shows the importance of protonation equilibria on XS À to modulate the stability of Cu(I)-clusters in such a way to guarantee the accumulation-storage-transport-release mechanisms for copper. At this regard, it has to be emphasized that a recent experimental (thermodynamic) work (Badarau & Dennison, 2011a) showed that the Cu(I) transfer in protein systems 'are dependent from pH and become less favorable as the pH is decreased below,' confirming the findings from the present theoretical study. The concept that relates the release of Copper(I) trough cluster breakage via protonation at sulfur atoms is pictured in Scheme 3. The protonation status at sulfur donors for the dimers and trimers found in solid state for CopZ were probably mixed ones, being Cu…Cu distances c. 3.1 Å from experimental data, c. 2.4 Å from computed not protonated models, and c. 3.6 Å totally protonated computed models (BS1 W) for cysteine derivatives. The increase of protonation at anionic S-ligands can be obtained by protein conformational changes and can cause the partial breakage of Cu 2 S 4 or Cu 3 S 6 clusters followed by copper release. Hence, the preferred forms of S-ligands are predictably not protonated at Cu(I)-transport steps and are protonated at Cu(I)-release steps. Extension on further studies on nature of chemical bond at clusters and on effect of local functionals like BLYP, as well as of salvation effect at attenuated ɛ values (with respect to 78.5, water), are in progress in this laboratory.

ΔEe
ΔEe  bridge the metals, thus relativistic effects should be reduced with respect to compact solid materials. The goal of the present work is not that of reaching the highest degree of accuracy for bond distances and formation energies; rather, we look for trends upon protonation at cysteine sulfur atom. On the other side, it is worth mentioning that the computed structures from the present work compare very well with copper trafficking sites in metal-protein systems reported by Badarau and Dennison (2011b) and Badarau et al. (2010). For example, in a metallo-chaperone Atx1 (Badarau et al., 2010) the H 2 Cys{Cu I (μS-H 2 Cys) 2 }H 2 Cys unit has bridging Cu-S bond distances by 2.25-2.39 Å: this means a good reliability of the present computational strategy for simulating copper chaperons and copper trafficking sites.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2012.689703.