DOI: 10.1002/cphc.201500169

Articles

Molecular Modeling of Bioorganometallic Compounds: Thermodynamic Properties of Molybdocene–Glutathione Complexes and Mechanism of Peptide Hydrolysis Dimas Su‚rez* and Natalia D†az[a] The computational study of bioinorganic complexes between transition metals and flexible ligands is still challenging, given that, besides requiring extensive conformational searches, the treatment of metal–ligand bonds demands the application of quantum chemical methods. Herein, the adducts formed between molybdocene, which exhibits antitumor activity and reacts with thiol groups to give stable water-soluble complexes, and the tripeptide glutathione, which is a major source of biological thiols, are studied. Conformational searches are

performed using the semiempirical PM6 method followed by geometry optimizations and single-point calculations using density functional theory methods. In addition, molecular dynamics simulations of the molybdocene–glutathione complex involved in the regioselective hydrolysis of the Cys–Gly linkage are performed in explicit solvent. The reactive process is also studied theoretically on cluster models of both the molybdocene-bound and the free peptide.

1. Introduction The high structural diversity of coordination compounds, the broad spectrum of metal–ligand interactions, together with the Lewis acidity, redox activity and spectroscopic properties of transition metals, among other properties, confer unique characteristics to the inorganic complexes that play essential roles in biological systems.[1] These features, which can be greatly expanded by taking advantage of biologically exotic metals and/or ligands, have allowed new possibilities in the design of artificial bioinorganic compounds that interact with specific proteins or nucleic acids.[2, 3] Thus, bioorganometallic compounds with low or no toxicity that are stable in air and in aqueous solutions have been introduced as potential antibiotics and drugs against cancer, HIV and malaria.[4] The tethering of organometallic complexes to protein-affinity ligands can also provide efficient and selective protein inhibitors that might be promising candidates for the design of catalytic drugs that specifically bind and degrade proteins that are implicated in disease.[5–7] Other compounds are of interest for biotechnological applications, for example, organometallic reagents capable of hydrolyzing the amide linkages of proteins under near physiological conditions of temperature and pH,[5, 8, 9] which can be useful for protein sequencing and peptide mapping. Metallocenes of the form Cp2MCl2 (Cp = h5-C5H5, M = Ti, V, Nb, Mo) are neutral complexes with a distorted tetrahedral geometry that constitute a bioactive family of inorganic com[a] Dr. D. Su‚rez, Dr. N. D†az Departamento de Qu†mica F†sica y Anal†tica Universidad de Oviedo Juli‚n Claver†a 8, 33006 Oviedo (Spain) E-mail: [email protected] Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201500169.

ChemPhysChem 2015, 16, 1646 – 1656

pounds. Their role as antitumor compounds, especially the titanocenes (M = Ti), has long been recognized.[10, 11] Furthermore, Cp2MCl2 complexes are considerably diverse in their biological and chemical activities as antitumor compounds, protein inhibitors, homogeneous catalysts, and so forth.[12–20] One of the most studied metallocenes in recent decades is molybdocene dichloride (Cp2MoCl2). Although it has poor solubility in water, this complex has attracted substantial attention as a watersoluble organometallic catalyst/promoter of important organic reactions,[12] such as nitrile hydration,[13] carbon monoxide oxidation,[14] and carboxylic and phosphate ester hydrolysis.[15] Both the hydrolytic chemistry of molybdocene,[15, 16] as well as its coordination chemistry[15, 17] with nucleobases, nucleotides, amino acids and small peptides, have been studied experimentally under physiological conditions with a special aim to determine which molybdocene species are likely to be formed in biological environments. From these studies, it is known that water molecules readily displace the two chloro ligands of Cp2MoCl2 and that the Cp2Mo2 + moiety is remarkably stable, unlike those of titanocenes. More specifically, the hydrolysis of Cp2MoCl2 at neutral pH leads to the formation of a [Cp2Mo(OH)(OH2)] + complex in equilibrium with a dimer complex, [Cp2Mo(m-OH)2–MoCp2]2 + .[15] [Cp2Mo(OH)(OH2)] + is the precursor of other molybdocene complexes formed with amino acids and peptides. Of the possible metal-coordinating groups in biomolecules (e.g. amino, carboxylate, thiol, phosphate (through O) and heterocyclic (N) groups), molybdocenes bind preferentially to deprotonated thiol groups, resulting in highly stable water-soluble complexes[18, 19] that are not cytotoxic and have thus far not found clinical application.[20] In this computational study, we analyzed the molybdocenethiol complexes formed from the tripeptide glutathione (g-Glu–Cys–Gly; GSH),[18, 19] which are of particular interest be-

1646

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles cause: 1) Cp2MoCl2 is expected to be converted to Cp2Mo(GSH)1–2 in vivo given that GSH is a major source of biological thiols[15] and 2) the Cp2Mo2 + moiety attached to the thiolate group of GSH promotes the regioselective hydrolysis of the Cys–Gly peptide bond at 40–60 8C.[19] Our main goal was twofold: first, to provide detailed information about the molecular structure, conformational properties and free energies of the various molybdocene-GSH complexes that have been detected experimentally in aqueous solution, and second, to determine the most likely reaction mechanism(s) and the associated energy barriers for the hydrolysis of the molybdocenebound GSH peptide. Experimental and computational modeling based on quantum chemical methods have been used in fruitful and complementary ways to study important organometallic complexes and catalysts that act upon relatively small and/or rigid organic ligands.[21, 22] Recent theoretical studies have rationalized the mechanisms of two molybdocene-assisted processes—the oxidation of carbon monoxide[23] and the hydration of nitriles.[24] However, from a computational chemistry point of view, the achievement of our goal of characterizing molybdocene–GSH complexes could have methodological interest given that de novo predictions of the structure and molecular properties of bioorganometallic complexes involving transition metals and flexible ligands (e.g. oligopeptides) remain challenging. The reason for this is that whereas extensive conformational search calculations are necessary for identifying the most populated conformers of each complex in aqueous solution, the presence of multi-hapto or metal–ligand bonding largely hinders the direct application of the fast molecular mechanics (MM) methods suitable for extensive sampling; more-expensive quantum chemical methods are required. As a compromise, the equilibrium geometries and free energies of the molybdocene–GSH complexes are determined using a recently developed protocol[25] that combines a molecular dynamics (MD) conformational search at the semiempirical quantum mechanical PM6 level with density functional theory (DFT) geometry optimizations and single-point calculations. Nonspecific solvent effects are included by means of a solvent continuum model coupled with the PM6 or DFT Hamiltonians. However, it was found that, prior to the quantum chemical study of the hydrolysis reaction, classical MD simulations in explicit solvent using MM potentials are also needed to study the hydration structure of the prereactive molybdocene complex. Overall, the results from our high-throughput computational approach, which combines a variety of quantum mechanical (QM) and MM approaches, will complement experimental information to give a complete view of the interactions between molybdocene and GSH, which, in turn, might be useful for assessing the biological activity of other molybdocene–thiol derivatives and/or their possibilities and limitations as synthetic analogues of metalloproteases. In addition, the details of the computational protocol and its performance might foster the intensive application of computational techniques for probing bioorganometallic systems.

ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

2. Results and Discussion 2.1. Molybdocene–GSH Complexes Molybdocene–GSH complexes are readily obtained at neutral pH and room temperature from equimolar solutions of Cp2MoCl2 and GSH.[18, 19] In these complexes, the Cp2Mo2 + moiety exhibits a large affinity for the GSH thiolate group, but it can also be coordinated by water/hydroxide ligands or the carboxylate and amide groups of GSH. 1H NMR measurements of the Gly-Ha and Cp chemical shifts suggest the presence of three 1:1 adducts ([Cp2Mo(OH)(GSH-S)], [Cp2Mo(GSH-S,O)] + and [Cp2Mo(GSH-S,N)]) and a 1:2 complex ([Cp2Mo(GSH-S)2]; by convention, the ligating atoms are indicated in italics).[19] To further characterize their structure and relative stability, we used the combined semiempirical and DFT computational protocol as described in the Computational Methods section. We assumed that the [Cp2Mo(GSH)] complexes are formed in aqueous solution upon the reaction between the thiol group of the GSH peptide and the [Cp2Mo(OH)(H2O)] + complex, the latter being the reactive species produced by the fast hydrolysis of the precursor molybdocene complex Cp2MoCl2.[16] Specifically, we focused on the 1:1 and 1:2 molybdocene–GSH adducts in which the Mo-bound thiol groups are deprotonated because the corresponding complexes with neutral thiol groups would be much less stable, as suggested by previous calculations and experimental results on the molybdocene–Cys system.[18, 25] The coordination modes and/or protonation states of the [Cp2Mo(GSH)1–2] species that we examined computationally are shown in Scheme 1, which also shows the relative free energies for every complex (these free energy values take into account the conformational averaging). The molecular geometries of the most stable conformer for each species are displayed in Figure 1. The results of the conformational calculations are summarized in Table S1 in the Supporting Information, which contains the relative values of various free energy components for the top 10 most stable conformers of each complex. In principle, formation of the [Cp2Mo(H2O)(GSH-S¢)] + complex (A in Scheme 1) can readily occur through a water–thiol exchange process at the Mo center followed by a fast proton transfer from the Mo-bound thiol group to the Mo¢OH fragment. The net reaction, [Cp2Mo(OH)(H2O)] + + GSH¢ ! [Cp2Mo(H2O)(GSH-S¢)] + H2O, is clearly exoergic, with a calculated DrG value of ¢11.1 kcal mol¢1 in terms of the Boltzmannaveraged free energies of the peptide and molybdocene species. The GSH peptide in complex A, which is coordinated to the metal center through a single Mo¢S bond, can act as a chelating ligand when its Gly carboxylate displaces the Mo-bound water molecule to give complex B. As shown in Figure 1, the presence of the Mo-OH2···¢OOC-Gly contact in the most likely conformer of A seems particularly sufficient for this process to occur. According to our calculations, the A!B transformation is favorable with a DrG value of ¢6.3 kcal mol¢1. As metal-bound water molecules can behave as weak acids, the A complex would be in equilibrium with its conjugate base A’. According to the calculated DrG value, the acid form of A is

1647

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles mol¢1, respectively. Experimentally, the abundances of the A, B, C, and D complexes in the reaction mixture at room temperature have been assigned values of 10 %, 15 %, 20 %, and 15 %, respectively.[19] The estimated errors in the B3LYP-D3 energies for related ligand-exchange processes[25–27] are at least 1–2 kcal mol¢1 and, most likely, additional errors are introduced by the continuum solvent model and the inclusion of the G(H + ) parameter. Similarly, a recent study assessing COSMO-DFT-based calculations of relative stability constants for metal-ion complexes has reported an average accuracy of 2 kcal mol¢1.[28] Hence, the magnitude of some of the computed DG differences ( … 2 kcal mol¢1) among the molybdocene complexes is similar to their expected uncertainty and, therefore, the theoretical results seem consistent with an equilibrium Scheme 1. Proposed equilibria for molybdocene–GSH complexes in water at neutral pH. Relative free energies mixture. In this respect, it might [kcal mol¢1] with respect to the starting reactants are shown. be interesting to note that the Boltzmann averaging of the free energies of the conformers favored by 5.2 kcal mol¢1 at pH 7.0 and, consequently, the can benefit from partial cancellation of random errors, as remetal-bound water molecule in the monodentate molybdocently shown in error analyses of computational methods.[29] cene–GSH complex would be preferentially protonated at neutral pH. Nevertheless, when we compared the geometries of Because the free energy of the bidentate complex B the lowest free-energy conformers of A and A’ (Figure 1), it [Cp2Mo(GSH-S,O)] + is dominated by a single conformer, wherewas clear that the charge neutrality of the Mo coordination as the other complexes populate several conformational states, sphere in A’ induces an important conformational change in the overestimation of the stability of B (¢6.3 kcal mol¢1 with rethe GSH moiety. Interestingly, the backbone NH group of the spect to A) could be due to a worse cancelation of errors, Cys residue in the most stable A’ conformer takes part in which, in turn, would be related to its sharp conformational a Mo-OH···HN-Cys contact and is well situated to anchor the distribution. peptide to the Mo center. It is then reasonable to hypothesize that the destabilization of the amide group induced by the 2.2. Molybdocene-Promoted GSH Hydrolysis: MoIV cation, together with the basicity of the Mo¢OH fragThermochemistry and Stability of Tetrahedral Intermediates ment, triggers the deprotonation of the metal-bound amide group and the release of a water molecule. Thus, we propose As mentioned in the Introduction, the reaction of Cp2MoCl2 that the conjugate base A’ would be the precursor for the forwith the GSH tripeptide at elevated temperature assists the mation of the chelate complex C, whose free energy in solucleavage of the Cys–Gly amide linkage.[19] Thus, 1H NMR meas¢1 tion is 2.6 kcal mol below that of the basic form A’ and urements on equimolar mixtures of Cp2MoCl2 and GSH kept at 60 8C after 1 h show a Gly-Ha signal that is characteristic of 2.6 kcal mol¢1 above that of the acid form A. Another reactive free Gly, whereas Cys-Ha displacements suggest the presence pathway accessible to the complex A’ consists of the ligand-exof the dipeptide complex [Cp2Mo(g-Glu–Cys-S,O)] + (E in change reaction with a second GSH molecule, leading to the formation of the 1:2 adduct (D in Scheme 1). The shift of the Scheme 1). hydroxo ligand in A’ by a thiolate group is also exoergic, the Before studying the reaction mechanism, we estimated the corresponding DrG value being ¢8.0 kcal mol¢1 for the process free energy of the hydrolysis products in order to characterize their thermodynamic stability. Thus, we applied the PM6-DFT A’ + GSH!D + H2O, similar to that of the initial reaction conformational search protocol to the bidentate complex E [Cp2Mo(OH)(H2O)] + + GSH!A + H2O. and estimated its average free energy in solution on the basis From the DrG values in Scheme 1, the stabilities of the B, C of the COSMO-B3LYP calculations. For zwitterionic Gly, we only and D complexes relative to A are ¢6.3, +2.5 and ¢2.8 kcal ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

1648

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles

Figure 1. Views of the COSMO-B3LYP/aug-cc-pVDZ-optimized geometries for the most stable conformers of the molybdocene–GSH complexes. Selected bond distances [æ] are also indicated.

considered a single conformer. The calculated DrG value for the hydrolytic cleavage A!E + Gly, amounts to ¢11.8 kcal mol¢1, which is significantly larger in absolute terms than the relative DG values of the molybdocene–tripeptide complexes. The overall free energy for the net reaction, [Cp2Mo(OH)(OH2)] + g-Glu–Cys–Gly![Cp2Mo(g-Glu–Cys-S,O)] + Gly + H2O, amounts to ¢22.9 kcal mol¢1 and, therefore, our calculations unambiguously predict that the hydrolysis of the GSH peptide promoted by molybdocene is thermodynamically favored, which is in agreement with experiment as 80 % of GSH is hydrolyzed after 200–250 h of reaction time. For comparison, we also calculated the free-energy change for the hydrolysis of the free GSH peptide. To this end, the conformational states and free energy of the g-Glu–Cys dipeptide were determined using the semiempirical-DFT protocol. The resulting DrG value for the reaction g-Glu–Cys–Gly + H2O!g-Glu–Cys + Gly is only ¢0.2 kcal mol¢1. Therefore, it is apparent that the metal complex, which remains bound to the dipeptide product, plays a crucial role in promoting the regioselective cleavage of the GSH peptide by increasing the thermodynamic driving force of the hydrolytic process by approximately 22 kcal mol¢1 in absolute terms. The fact that the GSH hydrolysis in the presence of molybdocene does not proceed appreciably at room temperature indicates that the reaction must overcome an important activation energy barrier. The mechanistic proposal assumes that the Mo-bound hydroxide would be the required nucleophile and that the intramolecular attack of Mo-OH towards the carbonylic C atom of the Cys residue would lead first to the formation ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

of a tetrahedral intermediate and the posterior release of Gly. In this scenario, the molybdocene–tripeptide complexes B and C, and the dimer complex D, would be inactive for hydrolysis.[19] Prior knowledge of the structure and stability of the putative tetrahedral intermediates should help us outline the best strategy for studying the reaction pathway. For this reason, we performed preliminary geometry optimizations of two tetrahedral intermediates built from the most stable conformers of A and A’, which differ in the presence of the metal-bound water or hydroxide moieties, respectively. In these calculations, some bond distances were initially constrained [e.g. Mo¢S(Cys) = 2.55 æ, Mo¢O(Nuc) = 2.15 æ, O(Nuc)¢C(Cys) = 1.40 æ; Nuc = nucleophilic] to facilitate optimization of the tetrahedral intermediate, but these constraints were removed during a second optimization stage. We also assumed that the formation of the tetrahedral intermediate from the lowest energy conformer of A is accompanied by proton transfer from the Mo-bound water to the carbonylic O atom in the scissile peptide bond. We found that the intermediate derived from A is associated to a stable energy minimum on the potential energy surface, which exhibits a short H-bond contact between the Gly carboxylate group and the metal-bound O atom. Its relative energy (DE) with respect to A is 19.6 kcal mol¢1 at the COSMOB3LYP/au-cc-pVTZ//COSMO-B3LYP/au-cc-pVDZ level, which is a moderate DE difference that suggests a relevant kinetic role for this intermediate (see Scheme 2). In contrast, the tetrahe-

Scheme 2. Proposed tetrahedral intermediates involved in the hydrolysis reaction of the molybdocene–GSH complex. The relative energy DE [kcal mol¢1] of the tetrahedral intermediate derived from A is also indicated.

dral intermediate from the A’ complex was shown to be highly unstable because 1) the relative energy of the intermediate structure optimized with the constrained bonds is + 45.9 kcal mol¢1, and 2) when the optimization is resumed with no constraints, the intermediate structure reverts to the initial complex A’. Therefore, our calculations strongly suggest that the hydrolytically active species corresponds to the A complex bearing the Mo¢OH2 nucleophile. 2.3. Molecular Dynamics Simulations of Prereactive Complexes The COSMO-B3LYP/aug-cc-pVDZ structure for complex A (Figure 1) does not constitute the optimum starting point for determining the reaction pathway(s) leading to the formation of the tetrahedral intermediate and its cleavage into the hydrolysis products. On the one hand, the relative orientation be-

1649

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles tween the Mo-bound water and the carbonylic O atom seems unsuitable for promoting the Mo¢OH2 !C=O nucleophilic attack. On the other hand, it is not obvious how to consider the kinetic role of explicit water molecules in assisting the required proton-transfer processes because complex A was obtained using a continuum solvent model. Furthermore, peptide cleavage in metal-assisted processes can also occur through the external attack of water molecules,[8, 30] that is, one water molecule from the first solvation layer could act as a nucleophile instead of the Mo-bound water. To investigate the kinetic role of ancillary water molecules, we followed a hybrid approach that consisted of adding one or more explicit water molecules to the reactive solute and placing the resulting complex within the same cavity surrounded by the continuum solvent. Of course, the number and initial location of the explicit waters must be selected carefully and justified according to structural and mechanistic criteria. To this end, we aimed to investigate the solvation of the [Cp2Mo(H2O)(GSH-S)] + complex by means of MD simulations in explicit solvent, which were performed using MM potentials for both the solute and the water molecules in order to allow for sufficient sampling (10 ns) of the solute conformational motions and the accompanying solvent rearrangements. Hence, as described in the Computational Methods section, we derived a transferable MM parameterization of the molybdocene–Cys complex from the results of QM calculations. During the MD simulation, the backbone of the GSH peptide, which is bound to the [Cp2Mo(OH2)2 + ] moiety, adopts preferentially extended conformations that are characterized by an average radius of gyration of 4.05 œ 0.20 æ. The representatives of the top five clusters extracted from the MD trajectory, which account for approximately 80 % of the sampled snapshots, are shown in Figure 2 a. The relative conformation of the Mo complex with respect to the GSH backbone is nearly identical in several of the cluster representatives, having a characteristic intramolecular Mo-OH2···¢OOC-Gly interaction. This contact and the elongated disposition of the GSH backbone were also observed in the most stable conformer of the A complex generated by the PM6-DFT conformational search calculations (Figures 1 and 2 a). Hence, such coincidence in the conformational predictions made using the PM6-DFT approach and the MM–MD calculations increases the confidence level of our computational results. Nonetheless, the MD simulation gives further information concerning the solvent structure and the abundance of intramolecular contacts. For example, the direct Mo-OH2···¢OOC-Gly contact is found in 24 % of the MD snapshots, whereas a 30 % abundance was measured for the Hbond between the Mo-bound water and the NH amide group of Gly. In addition, the Mo¢O···C¢Cys distance is shorter than 5.0 æ in 40 % of the analyzed snapshots, demonstrating that explicit solvent favors the approach of the Cys carbonyl group to the metal site. Furthermore, the MD simulation revealed that the carboxylate group of the Gly residue holds a water molecule in close proximity to the Cys–Gly amide group. This Gly-COO¢···H2O···C=O-Cys contact, which is found in 47 % of the MD structures, was also observed in the most populated cluster representative (Figure 2 b), suggesting that the external ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

Figure 2. a) Superposition of the cluster representatives corresponding to the most populated clusters extracted from the [Cp2Mo(OH2)(GSH)] simulation. The thickness of the models corresponds to the population of each cluster. b) Ball-and-stick view of the Mo site in the representative structure of the most populated cluster. Selected specific distances [æ] are indicated. The dummy atom at the geometric center of each cyclopentadienyl ring mimics the presence of the h5-Cp···Mo bond. c) Superposition of the cluster representatives from the GSH simulation. d) Ball-and-stick view of the selected GSH cluster for studying the hydrolysis of the Cys–Gly amide bond.

water could act as the nucleophile. Alternatively, the same external water is also suitably positioned to assist the proton transfers that likely occur upon the nucleophilic attack of the Mo-bound water. The conformational properties of the free GSH peptide in explicit solvent were also investigated using MM–MD methods. Again, the structure of the most populated clusters (Figure 2 c) is similar to that of the most stable conformer identified using the PM6-MD approach. According to the structural analyses, 35 % of the MD snapshots show a water molecule located in the second solvation shell around the carbonylic O atom of the Cys residue, which is properly positioned to act as the nucleophile that attacks the Cys C atom. The nucleophilic water could also relay a proton to the carbonylic O atom with the assistance of a bridging water molecule that interacts directly with that O atom (see the hydrated cluster representative shown in Figure 2 d). 2.4. GSH Hydrolysis: Reaction Mechanisms To determine the most likely reaction mechanism of the molybdocene-promoted hydrolysis of the GSH peptide, we built a cluster model from the MD snapshot shown in Figure 2 b, selecting the molybdocene complex and the ancillary water molecule. Two reaction pathways differing in the identity of the nucleophile were characterized in terms of the most important transition structures (TSs) and tetrahedral intermediates (Scheme 3). The optimized structures at the COSMO-B3LYP/

1650

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles On the COSMO-B3LYP/aug-ccpVDZ potential energy surface, TS2ext is connected with a product complex, Pext, in which a neutral Gly molecule is H-bonded with the neutral carboxylic group of the Cys residue. This intermolecular complex is predicted to be less stable than the initial complex Cext by 15.2 kcal mol¢1. The release of the final products requires the dissociation of Pext, several acid–base reactions and the formation of the bidentate [Cp2Mo(g-Glu–CysS,O)] + structure through an intramolecular ligand-exchange process. We did not study these transformations in detail because they can take place through manifold reaction pathways and are most likely fast processes Scheme 3. Outline of the external (ext) and intramolecular (intra) reaction mechanisms for the hydrolysis of the that overcome smaller energy [Cp2Mo(OH2)(GSH)] + complex. barriers than that of the rate-determining TS for hydrolysis. We computed instead the DrG value associated with the separate products, the resulting value aug-cc-pVDZ level of theory and their relative free energies for being clearly exoergic, ¢24.6 kcal mol¢1 (this free-energy differthe two mechanisms are shown in Figure 3 (the mechanisms are labeled “ext” for nucleophile = external water and “intra” ence does not take into account the conformational Boltzfor nucleophile = Mo-bound water). As described in the Commann averaging). putational Methods section, single-point COSMO-PWPB95-D3/ The second mechanistic route is characterized by the particiaug-cc-pVTZ calculations in combination with COSMO-B3LYP/ pation of the external water molecule in assisting various cc-pDVZ thermal corrections at 60 8C were used to predict the proton-transfer events. To promote the nucleophilic attack of free energy in solution of all the critical structures. the metal-bound water with simultaneous Mo-OH2 !O=C Globally, the COSMO-B3LYP geometry of the prereactive proton transfer, a H-bond bridge (i.e. Mo-OH2···OH2···O=C) complex Cext is close to that of the parent structure taken from seems a prerequisite. Such a contact was not observed during the MD trajectory. The most significant difference arises from the MD simulation of the [Cp2Mo(H2O)(GSH-S)] + complex and, therefore, we manually changed the H-bond pattern of the exthe positioning of the external water molecule, which forms ternal water at Cext and performed a full geometry optimization a new H2O···H2O-Mo contact that, in turn, is probably due to that produced a second prereactive complex, Cintra, which rethe continuum solvent model. Nevertheless, we found that Cext is the prereactive structure leading to a TS structure (TS1ext) for tains the desired H-bond bridge and is 4.9 kcal mol¢1 less stable than Cext. The Cintra structure would be the immediate the nucleophilic attack of the external water toward the carprecursor of the TS for the formation of the tetrahedral interbonylic C atom with simultaneous proton transfer from the mediate (TS1intra in Figure 3), which has a DG barrier of Mo-bound water to the leaving amino group of the Gly residue and, interestingly, with extra assistance of the Gly carbox28.5 kcal mol¢1 with respect to Cext. The equilibrium geometry ylate group that deprotonates the attacking water molecule of TS1intra indicates that the water-assisted proton transfer from (Figure 3 a). The calculated free-energy barrier for TS1ext is Mo-OH2 to the carbonylic O atom is essentially completed while the transition vector is dominated by the formation of 34.4 kcal mol¢1. This transition state is then connected with the O···C bond. In this respect, we note that a TS for the a zwitterionic tetrahedral intermediate (Iext), which has similar proton-transfer event could not be optimized at the COSMOgeometry and stability (Iext is only 2.1 kcal mol¢1 below TS1ext). B3LYP/aug-cc-pVDZ level. We performed an intense computaMoreover, the dissociation of Iext can occur readily passing tional search for such a structure, which did not render a true through a TS for the cleavage of the amide C¢N bond (TS2ext), critical structure, but a shallow region of the potential energy which is only 3.1 kcal mol¢1 less stable than the intermediate surface, at a DG value of approximately 24 kcal mol¢1, which is (Figure 3). Hence, it can be considered that the series of critical structures, TS1ext !Iext !TS2ext, represents a concerted route for characterized by a Mo-OH···H3O + ···O=C configuration. We conclude therefore that the Mo-OH2 !O=C proton transfer medithe hydrolysis of the GSH peptide without accumulation of intermediates. ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

1651

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles Again this process can be favored by the assistance of the external water. However, the positioning of the external water at Iintra is inadequate for the COH!NH proton transfer pathway and, consequently, we optimized another conformer of the tetrahedral intermediate, Iintra*, in which the external water is placed in the required prereactive position. Because Iintra* and Iintra have a similar free energy (Iintra is only 1.3 kcal mol¢1 less stable than Iintra*), it seems that the solvent rearrangement necessary for the cleavage of the tetrahedral intermediate could occur readily. The Iintra* intermediate is connected with a concerted TS for the water-assisted C-OH!NH proton transfer and the simultaneous cleavage of the C¢N bond (TS2intra in Figure 3 b), which has a DG barrier of 33.4 kcal mol¢1 with respect to Cext (17.5 kcal mol¢1 with reference to Iintra*). These figures indicate that the neutral intermediates Iintra and Iintra* are located in a deep free-energy well delimited by the TSs for their formation/cleavage and that the cleavage of Iintra* is the rate-determining step of the stepwise mechanism associated with the nucleoFigure 3. COSMO-B3LYP/aug-cc-pVDZ structures of the transition structures and intermediates involved in the mophilic attack of the Mo-bound lybdocene-promoted hydrolysis of the [Cp2Mo(OH2)(GSH)] complex. Selected interatomic distances are expressed water. Finally, the critical point ¢1 in æ. Relative free energies [kcal mol ] of the structures with respect to the Cext structure are shown. Pintra derived from TS2intra is a ternary complex formed by zwitterated by the external water is the event that triggers the nucleionic Gly, water and the molybdocene–dipeptide complex ophilic attack described by TS1intra. showing a bidentate mode of binding, its relative DG value According to our calculations, the tetrahedral intermediate being ¢3.2 kcal mol¢1. Therefore, the dissociation of Pintra into ¢1 Iintra connected with TS1intra is 9.6 kcal mol above Cintra. The the separate fragments is the unique step that is formally needed to yield the final hydrolysis products. structure and stability of Iintra, which is a neutral tetrahedral inOur DFT calculations indicate that the stepwise intramolecutermediate characterized by the presence of a 1,2 gem-diol lar mechanism for the cleavage of the Cys–Gly linkage, which functionality, are similar to those of the tetrahedral intermedihas a rate-determining DG barrier of 33.4 kcal mol¢1, is 2.0 kcal ate derived from the A complex (Scheme 2), but including now the external water molecule bridging the two alcohol groups. mol¢1 more favorable than the concerted route associated with We point out again that a remarkable feature in the structure the attack of the external water. These theoretical results can of Iintra is the close H-bond contact between the Gly carboxylbe compared with experiment. Thus, in the presence of molybdocene, 25 % of the GSH molecules are hydrolyzed within 20 h ate and the Mo-OH moiety (the Gly-COO¢···HO-Mo distance is at neutral pH and 60 8C.[19] Using the thermodynamic formulaonly 1.48 æ; 2.51 æ between the heavy atoms) that probably ¢1 contributes to its stabilization. In fact, Iintra is 13.9 kcal mol tion of transition-state theory and assuming a pseudo-firstorder rate law r = k*[A], where k* = k[H2O], the experimentally more stable than TS1intra. determined kinetic data yields an activation free energy of The productive cleavage of Iintra requires a 1,3-transposition 32.8 kcal mol¢1, which is close to our theoretical estimate. of H between one alcohol group and the leaving amino group. ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

1652

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles To better assess the actual kinetic effect of the molybdocene complex on the peptide hydrolysis reaction, we also studied the water-assisted hydrolysis of the free GSH peptide. On the basis of the structural data produced by the GSH MD simulation in explicit solvent, we built a cluster model comprising the GSH peptide and two water molecules (see Figure 2 d) that solvate the carbonyl group of the Cys residue. From this structure, the DFT geometry optimizations and energy calculations outline a water-assisted stepwise mechanism via a neutral 1,2gem-diol tetrahedral intermediate (Figure 4). This reaction

Figure 4. COSMO-B3LYP/aug-cc-pVDZ structures of the transition structures and intermediates involved in the water-assisted hydrolysis of GSH. Selected interatomic distances are expressed in æ. Relative free energies [kcal mol¢1] of the structures with respect to the CGSH structure are shown.

pathway, which resembles the intramolecular mechanism of molybdocene-promoted hydrolysis, and other processes theoretically determined for the neutral hydrolysis of amide bonds,[31] has a rate-determining energy barrier of 37.6 kcal mol¢1 that corresponds to the TS for the formation of the tetrahedral intermediate (I1GSH in Figure 4), which is 23.9 kcal mol¢1 less stable than the initial complex. During the second step, the TS2GSH structure for the water-assisted Cys-C¢OH! HN-Gly proton transfer (DG = 31.6 kcal mol¢1) precedes a weakly stable tetrahedral intermediate that is only 0.9 kcal mol¢1 below the TS for the cleavage of the C¢N bond. The neutral tetrahedral intermediates I1GSH–I1*GSH are 8–10 kcal ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

mol¢1 less stable than the molybdocene-bound intermediates, and the free energy barrier of the water-assisted process is 4.2 kcal mol¢1 higher than that of the molybdocene-assisted reaction. Hence, the molybdocene moiety has a moderate kinetic effect on the hydrolysis of GSH.

3. Conclusions Our computational results yield structural and energetic information that allow us to fill in some of the gaps from the experimental characterization of the mixture of 1:1 and 1:2 adducts formed between Cp2MoCl2 and GSH. Thus, we propose the following sequence of equilibrium processes for the interconversion of the 1:1 complexes: B$A$A’$C, in which the A complex, [Cp2Mo(OH2)(GSH-S)], plays a pivotal role (Scheme 1). Such a proposal is supported by the structural features of the most stable conformers shown in Figure 1: the Gly carboxylate group in A is appropriately oriented to displace the water ligand, whereas the Cys–Gly amide group in the A’ complex interacts with the basic Mo¢OH group, which is sufficient for the formation of the bidentate complex C [Cp2Mo(GSH-S,N)]. Another piece of information provided by the calculations is that the aquo complex A would be the dominant monodentate 1:1 adduct at neutral pH given that the equivalent hydroxo complex A’ is predicted to be less stable by 5.2 kcal mol¢1. Other free energy differences among the 1:1 (A, C) and the 1:2 (D) adducts are within the accuracy limit of the DFT calculations ( … 2–3 kcal mol¢1), except that for the chelate complex B (¢6.3 kcal mol¢1 with respect to A), the stability of which is probably overestimated. Taking into account both the similarity and uncertainty of the free-energy values, the energetic predictions seem in reasonable agreement with experiment. From the semiempirical-DFT calculations, the most stable conformer of the A adduct can be identified as the prereactive species that undergoes regioselective cleavage. This identification, which refines previous experimental proposals,[19] is supported by the fact that the intramolecular Mo¢OH2 !C=O nucleophilic attack would lead to relatively stable tetrahedral intermediates, in contrast with those presumably formed from the hydroxo complex A’. Furthermore, assigning A as the prereactive species instead of A’ is useful for interpreting the kinetic pH effect on the hydrolysis rate, which nearly doubles if the pH is reduced from 7.4 to 5.2,[19] as a direct consequence of an equilibrium shift C + H + !A, which raises the molar concentration of A. To investigate the mechanism for the molybdocene-assisted hydrolysis of the Cys–Gly amide bond, one or more explicit water molecules should be incorporated into the molecular models because external water can facilitate proton transfer events. To build such a model in a non-arbitrary manner, the hydration structure around the A complex was first studied by means of extended MD simulations in explicit solvent using MM potentials. After having derived a transferable and robust MM parameterization of the molybdocene–thiol complexes, the 10 ns MD simulation of the [Cp2Mo(H2O)(GSH-S)] + complex revealed that the terminal carboxylate group of the Gly residue binds an external water molecule that is ideally located to play

1653

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles an active kinetic role. This was confirmed by subsequent DFT calculations of the reaction profiles in which the external water molecule can act as a nucleophile (mechanism ext) or as a bifunctional catalyst in assisting the required proton-transfer steps with the Mo-bound water acting as the nucleophile (mechanism intra). The magnitude of the rate-determining DG barrier, 33.4 kcal mol¢1 for the intra pathway, is in agreement with the experimental estimation of the activation free energy, 32.8 kcal mol¢1. We also note that, although the external attack of the water molecule has a larger barrier, the actual DDG difference is not large (2.0 kcal mol¢1) and we cannot rule out that the ext and intra mechanisms are competitive. In fact the nucleophilic attack of the external water with simultaneous protonation of the leaving amino group by Mo-OH2 could be favored by the conformational properties and hydration structure of the prereactive complex A. The two mechanisms differ dramatically in the stability of the tetrahedral intermediates: whereas the ext mechanism is essentially a concerted process, the tetrahedral intermediate in the intra route would have an extended lifetime. Nevertheless, they also share three essential features: 1) the participation of the Cp2Mo(OH2)2 + fragment attached to the thiolate group of the Cys residue, 2) the active kinetic role of the external water and, equally important, 3) the presence of the C-terminal carboxylate group, which stabilizes the prereactive conformation of A, captures the reactive external water, deprotonates the external water in the ext mechanism and accepts a proton from the Mo-OH fragment in the tetrahedral intermediate of the stepwise intra route. Therefore, the theoretical results suggest that the molybdocene-assisted hydrolysis of peptide bonds can only occur if a solvent-exposed basic group is positioned nearby the thiol ligand. As well as illustrating the molecular details of the peptide hydrolysis reaction, the computational analyses indicate that the molybdocene complex has a moderate kinetic effect. This is clearly seen from a comparison of the critical structures and the free-energy profile of the intra route (Figure 3) with those in the reference hydrolysis pathway involving the free GSH peptide, nucleophilic water and one ancillary water molecule (Figure 4). Thus, coordination of the GSH thiolate group to the [Cp2Mo(OH2)]2 + complex reduces the global free-energy barrier by approximately 4 kcal mol¢1. Although this DDG value results from DFT calculations on the cluster models extracted from individual MD snapshots, and therefore does not take into account conformational averaging effects or the entropic cost for the transfer of water molecules from the bulk to the TSs, it is in agreement with the relatively slow kinetics measured experimentally. Conversely, the calculated DrG value for the hydrolysis reaction [Cp2Mo(OH)(OH2)] + g-Glu–Cys–Gly! [Cp2Mo(g-Glu–Cys–S,O)] + Gly + H2O is ¢22.9 kcal mol¢1, which is somewhat exoergic. As the corresponding DrG in the absence of molybdocene has a value of ¢0.2 kcal mol¢1, molybdocene, which remains anchored to the thiol group of the g-Glu–Cys dipeptide, acts preferentially as a thermodynamic promoter, reinforcing the spontaneity of the GSH hydrolysis. Despite the remaining challenge of the quality of the predicted relative energies among structurally different bioorganometallic complexes, we believe that the theoretical apChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

proach used in this work, which combines the sampling capability of the fast semiempirical QM methods with the structural and energetic predictions made using modern DFT methods, plus the effective description of solvent effects provided by continuum/hybrid models, could be further developed into new high-throughput computational protocols specifically designed to determine the molecular properties of bioinorganic complexes with flexible ligands, as in the recently proposed protocol for estimating the thermodynamics of metabolic reactions.[32] Indeed, such approaches might assist in the search for variants of bioorganometallic compounds through rational ligand design.

Computational Methods Conformational Analyses To find the most likely conformers of the different protonation states and coordination modes of the molybdocene–GSH complexes examined in this work, we used the semiempirical DFT protocol that has been presented in full detail and validated in previous work,[25] and therefore, we only outline here its most important features. The same method was applied to study the conformers of the unbound GSH and g-Glu–Cys peptides in aqueous solution. For each coordination complex or isolated peptide, we built its initial molecular geometry using the GabEdit graphical interface.[33] A ground-state singlet configuration was assigned to all the MoIV d2 complexes. GabEdit was used again to drive a 500 ps (1000 ps for [Cp2Mo(GSH)2] and GSH/g-Glu–Cys) MD simulation, in which the potential energies and forces were calculated externally by the MOPAC2009 program[34] using the semiempirical PM6 Hamiltonian[35] in combination with the COSMO solvation model.[36] A time step of 1 fs was used and the temperature (400 K) was controlled by Andersen’s algorithm, selecting a frequency of stochastic collisions of 20 ps¢1. The GabEdit conformational tool selected a total of 50 snapshots (100 for [Cp2Mo(GSH)2]) that were minimized and scored in terms of their relative COSMO-PM6 energies. After filtering the redundant conformers, all the structures were minimized at the COSMO-B3LYP/cc-pVDZ level of theory[37] and further characterized by numerical frequency calculations using the TURBOMOLE suite of programs (version 5.9).[38] The cc-pVDZ basis set[39, 40] was used for the nonmetal atoms plus cc-pVDZ-PP[41] for Mo, in which the valence electrons are described explicitly by the cc-pVDZ basis set and the core electrons are modeled by the double-z Stuttgart–Koln energy-consistent relativistic pseudopotential. The COSMO method[42] was used to take into account continuum solvent effects on molecular geometries and energies. All the conformers were reoptimized at the COSMO-B3LYP/aug-ccpVDZ level of theory, which improves the level of theory by adding diffuse basis-set functions for the nonmetal atoms. To further refine the quality of the electronic energies, we also performed single-point COSMO-B3LYP/aug-cc-pVTZ calculations on the COSMO-B3LYP/aug-cc-pVDZ geometries. The former basis set combines aug-cc-pVTZ[39, 40] for the nonmetal atoms with cc-pVTZ-PP[41] for Mo. To take into account the important dispersion interactions, we computed the atom-pairwise DFT-D3 dispersion correction (Edisp).[43] The single-point DFT-D3 calculations were done on the B3LYP/aug-cc-pVDZ geometries using the Becke–Johnson damping function[44] to avoid near singularities for small interatomic distances. Due to the large size of the [Cp2Mo(GSH)2] complex, only the

1654

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles top 10 most stable conformers at the COSMO-B3LYP/cc-pVDZ level were considered for the subsequent calculations. In previous work, we have examined the ability of the B3LYP functional complemented with the DFT-D3 dispersion correction for studying mononuclear [Cp2Mo(L)(L’)] complexes, where L, L’ = H2O, OH¢ , H2S, SH¢ and Cl¢ .[25, 26] For these complexes, the B3LYP + DFTD3 calculations provide reaction energies for the relevant protonation and ligand-exchange processes in the gas phase that have a mean unsigned difference of 1.9 kcal mol¢1 with respect to reference values predicted by high-level CCSD(T) calculations. The free energy in aqueous solution of the ith conformer (Gi) for a given complex was estimated by means of the following composite approximation [Eq. (1)]: disp Gi ¼ GCOSMO-B3LYP=aug-cc-pVTZ þ Gtherm COSMO-B3LYP=cc-pVDZ þ E DFT-D3

ð1Þ

where GCOSMO-B3LYP/aug-cc-pVTZ is the sum of the electronic energy plus the solvation free energy in aqueous solution obtained at the COSMO-B3LYP/aug-cc-pVTZ level, Gtherm COSMO-B3LYP=cc-pVDZ is the thermal free-energy correction computed at the COSMO-B3LYP/cc-pVDZ level within the rigid rotor and harmonic oscillator approximations, and E disp DFT-D3 is the dispersion correction. The translational entropy ¢1 ¢1 included in Gtherm to COSMO-B3LYP=cc-pVDZ was reduced by 6.4 kcal mol K take into account the corresponding decrease in the translational entropy ongoing from the gas-phase standard state (1 bar, 298 K, 0.041 m; as typically considered by quantum chemical software) to the required 1 m standard state in solution. For water, a pure liquid standard concentration of 55 m was used. The Gi energies were averaged over the various conformers of each configuration according to the Boltzmann distribution, and the result was augmented with the ¢TSconform term [Eq. (2)]: ‡¼ G

X

pi Gi ¢ TSconform

ð2Þ

i

where pi is the Boltzmann statistical weight of the i-conformer and Sconform accounts for the conformational entropy of the set of conformers,[45] which was calculated according to the Shannon equation [Eq. (3)]: Sconform ¼ ¢R

X

pi lnpi

i

A set of 450 MD snapshots equally spaced during the last 9 ns were grouped into different clusters according to the criterion of a fixed radius of displacement of the backbone with respect to the values of the root-mean-square deviation (RMSD) for the F and Y backbone angles.[51] The structure in each cluster with the lowest deviation is taken as the cluster representative. Other structural analyses were performed using the PTRAJ module of AMBER 10.

ð4Þ

From a representative cluster generated by the MD simulation, we built a cluster model that included the [Cp2Mo(OH2)(GSH-S)] solute and one ancillary water located in the vicinity of the Cys–Gly amide linkage. For this model, we carried out full geometry optimizations of the stable species and transition states that are located along the reaction pathways for the hydrolysis of the Cys–Gly bond. For comparative purposes, we also characterized a water-assisted pathway for the hydrolysis of the free GSH peptide in aqueous solution.

Some free-energy differences reported in this work involve the standard free energy of a proton in aqueous solution. To estimate this quantity, we combined the gas-phase free energy of a proton [(5/2)RT¢TSgas = 1.48–7.76 = ¢6.28 kcal mol¢1 at 298 K, 1 bar] with its solvation free energy [DGsolv(H + )], which, in turn, was considered as a parameter to reproduce the experimental pKa value (5.5) for the first dissociation of the mononuclear [Cp2Mo(H2O)2]2 + (aq) complex.[16] By combining the COSMO-B3LYP free energies of the [Cp2Mo(H2O)2]2 + (aq) and [Cp2Mo(OH)(H2O)] + (aq) species with the experimental pKa value, we derived a DGsolv(H + ) value of ¢253.3 kcal mol¢1, which is close to typical literature values.[46] As similarly done in previous work,[26] we added the term ¢2.303 RT pH to G(H + ) to correct the standard state of H + at pH 7.0. www.chemphyschem.org

The initial structures of the GSH peptide and the [Cp2Mo(OH2)(GSH-S)] complex were taken from the semiempirical DFT conformational calculations and immersed in a octahedral box containing approximately 5000 TIP3P[47] water molecules. Energy minimizations and MD simulations were carried out using the PMEMD program included in the AMBER 10 suite.[48] The full systems were initially relaxed by means of 2500 energy-minimization steps and heated gradually to 300 K during 50 ps of MD. The SHAKE algorithm was used to constrain all R¢H bonds, and periodic boundary conditions were applied to simulate a continuous system. A cutoff of 10.0 æ was defined for the nonbonded interactions and the particle-mesh Ewald (PME) method was used to include the long-range contributions.[49] During the MD simulations, the Berendsen and Langevin methods were used to control the pressure (1 bar) and the temperature (300 K) of the system, respectively.[50] Ten-nanosecond trajectories were computed with a time step of 2 fs and coordinates were saved for analysis every 1000 steps.

Reaction Mechanisms

i

ChemPhysChem 2015, 16, 1646 – 1656

Classical MD simulations in the presence of explicit solvent were used to investigate the solvation of the most likely prereactive structures of the free GSH peptide and the [Cp2Mo(OH2)(GSH-S)] complex. A laborious parameterization work was needed to properly simulate the coordination environment of the MoIV ion using MM methods. The details of the MM parameterization and the resulting parameters, which might be of interest for further computational studies aimed at ascertaining the biological impact and mechanism of action of molybdocene compounds in their interactions with protein targets (e.g. protein kinase C, Cys residues in tubulin and human serum albumin),[15, 17] are reported in the Supporting Information.

ð3Þ

where [Eq. (4)]: ¨ G¦ exp ¢ RTi pi ¼ P ¨ G¦ exp ¢ RTi

Molecular Dynamics Simulations in Explicit Solvent

The critical structures were first optimized at the COSMO-B3LYP/ccpVDZ level, followed by numerical computations of harmonic vibrational frequencies to characterize them as minima or transition states on the potential energy surface, as well as to compute thermal corrections to the free energy (a value of T = 333.15 K was chosen to compare with experimental kinetic data). The critical structures were reoptimized at the COSMO-B3LYP/aug-cc-pVDZ level. Finally, we performed single-point COSMO-PWPB95-D3/augcc-pVTZ calculations on the COSMO-B3LYP/aug-cc-pVDZ geometries. The free energy of the reactive and product complexes, intermediates and transition states, is estimated by adding the COSMO-

1655

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Articles B3LYP/cc-pVDZ thermal contribution to the COSMO-PWPB95-D3 free energy term (i.e. G = GCOSMO-PWPB95-D3/aug-cc-pVTZ + ). Gtherm COSMO-B3LYP=cc-pVDZ The PWPB95-D3 method[52] combines re-parameterized versions of the Perdew–Wang GGA-exchange functional,[53] the Becke95 metaGGA correlation functional,[54] and the DFT-D3 dispersion correction. It is a double-hybrid density functional that includes a second-order perturbation treatment (PT2) based on Kohn–Sham orbitals and adopts the spin-opposite scaled approximation. According to previous benchmark studies, the PWPB95-D3 method is an accurate and robust double-hybrid DFT method for computing thermochemical and activation energies involving molecular systems with main group elements,[55] the estimated mean error of the PWPB95-D3 energies being below 2.0 kcal mol¢1 when tripleor quadruple-z basis sets are used. The PWPB95-D3 method exhibits almost the same performance for the calculation of energy barriers of various bond-activation reactions catalyzed by Ni and Pd cations.[27] Therefore, the PWPB95-D3 DFT method is likely to provide a free-energy profile for the reaction mechanisms studied in this work with similar accuracy. The COSMO-B3LYP geometry optimizations and numerical frequency calculations were made using the TURBOMOLE package. The PWPB95-D3 calculations were performed with the ORCA 3.0.3 program system.[56] The resolution-of-the-identity approximation[57] was activated for the calculation of the PT2 term using the appropriate auxiliary basis set.[58] Solvation effects were accounted for by the COSMO method as implemented in the ORCA program, which produces a solvated set of Kohn–Sham orbitals that are later used to perform a normal PT2 step.

Acknowledgements We are grateful to the Ministry of Science and Innovation (Spain) for financial support (Project CTQ2010-18231). Keywords: density functional calculations · molecular modeling · organometallic compounds · peptides · semiempirical calculations [1] Biological Inorganic Chemistry: Structure and Reactivity (Eds.: I. Bertini, H. B. Gray, E. I. Stiefel, J. S. Valentine), University Science Books, Sausalito, CA, 2007. [2] K. L. Haas, K. J. Franz, Chem. Rev. 2009, 109, 4921 – 4960. [3] C.-M. Che, F.-M. Siu, Curr. Opin. Chem. Biol. 2010, 14, 255 – 261. [4] C. G. Hartinger, P. J. Dyson, Chem. Soc. Rev. 2009, 38, 391 – 401. [5] J. Suh in Comprehensive Inorganic Chemistry II, Vol. 6, 2nd ed. (Ed.: J. Reedijk, K. Poeppelmeier), Elsevier, Amsterdam, 2013, pp. 779 – 803. [6] J. Suh, W. S. Chei, Curr. Opin. Chem. Biol. 2008, 12, 207 – 213. [7] J. Prakash, J. J. Kodanko, Curr. Opin. Chem. Biol. 2013, 17, 197 – 203. [8] N. M. Milovic´, N. M. Kostic´, J. Am. Chem. Soc. 2003, 125, 781 – 788. [9] N. M. Milovic´, N. M. Kostic´, J. Am. Chem. Soc. 2002, 124, 4759 – 4769. [10] P. Koepf-Maier, H. Koepf, Chem. Rev. 1987, 87, 1137 – 1152. [11] G. Gasser, I. Ott, N. Metzler-Nolte, J. Med. Chem. 2011, 54, 3 – 25. [12] K. L. Breno, T. J. Ahmed, M. D. Pluth, C. Balzarek, D. R. Tyler, Coord. Chem. Rev. 2006, 250, 1141 – 1151. [13] T. J. Ahmed, S. M. M. Knapp, D. R. Tyler, Coord. Chem. Rev. 2011, 255, 949 – 974. [14] K. L. Breno, M. D. Pluth, C. W. Landorf, D. R. Tyler, Organometallics 2004, 23, 1738 – 1746. [15] J. B. Waern, M. M. Harding, J. Organomet. Chem. 2004, 689, 4655 – 4668. [16] C. Balzarek, T. J. R. Weakley, L. Y. Kuo, D. R. Tyler, Organometallics 2000, 19, 2927 – 2931. [17] E. Mel¦ndez, J. Organomet. Chem. 2012, 706 – 707, 4 – 12.

ChemPhysChem 2015, 16, 1646 – 1656

www.chemphyschem.org

[18] J. B. Waern, M. M. Harding, Inorg. Chem. 2004, 43, 206 – 213. [19] A. Erxleben, Inorg. Chem. 2005, 44, 1082 – 1094. [20] J. B. Waern, C. T. Dillon, M. M. Harding, J. Med. Chem. 2005, 48, 2093 – 2099. [21] M. Torrent, M. Sol—, G. Frenking, Chem. Rev. 2000, 100, 439 – 494. [22] D. Balcells, E. Clot, O. Eisenstein, Chem. Rev. 2010, 110, 749 – 823. [23] E. Tilvez, M. I. Menendez, R. Lopez, Eur. J. Inorg. Chem. 2012, 4445 – 4453. [24] E. Tilvez, M. I. Menendez, R. Lopez, Organometallics 2012, 31, 1618 – 1626. [25] D. Su‚rez, N. D†az, R. Lûpez, J. Comput. Chem. 2014, 35, 324 – 334. [26] D. Su‚rez, M. Zakarianezhad, R. Lûpez, Theor. Chem. Acc. 2013, 132, 1409. [27] M. Steinmetz, S. Grimme, ChemistryOpen 2013, 2, 115 – 124. [28] O. Gutten, L. Rul†sˇek, Inorg. Chem. 2013, 52, 10347 – 10355. [29] J. C. Faver, W. Yang, K. M. J. Merz, J. Chem. Theory Comput. 2012, 8, 3769 – 3776. [30] V. Yeguas, P. Campomanes, R. n. Lûpez, N. D†az, D. Su‚rez, J. Phys. Chem. B 2010, 114, 8525 – 8535. [31] V. Navr‚til, V. Klusak, L. Rulisek, Chemistry 2013, 19, 16634 – 16645. [32] A. Jinich, D. Rappoport, I. Dunn, B. Sanchez-Lengeling, R. OlivaresAmaya, E. Noor, A. B. Even, A. Aspuru-Guzik, Sci. Rep. 2014, 4, 7022. [33] A.-R. Allouche, J. Comput. Chem. 2011, 32, 174 – 182. [34] MOPAC2009, version 11.038 L, J. J. P. Stewart, Stewart Computational Chemistry, Colorado Springs, CO, 2009. [35] J. J. P. Stewart, J. Mol. Model. 2007, 13, 1173 – 1213. [36] A. Klamt, G. Schììmann, J. Chem. Soc. Perkin Trans. 1993, 2, 799 – 805. [37] A. D. Becke, J. Chem. Phys. 1993, 98, 5648 – 5652. [38] R. Ahlrichs, M. B•r, M. H•ser, H. Horn, C. Kçlmel, Chem. Phys. Lett. 1989, 162, 165 – 169. [39] T. H. Dunning Jr., J. Chem. Phys. 1989, 90, 1007 – 1013. [40] R. A. Kendall, T. H. Dunning Jr., R. J. Harrison, J. Chem. Phys. 1992, 96, 6796 – 6802. [41] K. A. Peterson, D. Figgen, M. Dolg, H. Stoll, J. Chem. Phys. 2007, 126, 124101. [42] A. Sch•fer, A. Klamt, D. Sattel, J. C. W. Lohrenzc, F. Eckert, Phys. Chem. Chem. Phys. 2000, 2, 2187 – 2193. [43] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104. [44] E. R. Johnson, A. D. Becke, J. Chem. Phys. 2006, 124, 174104. [45] D. Su‚rez, N. D†az, WIREs Comput. Mol. Sci. 2015, 5, 1 – 26. [46] G. J. Tawa, I. A. Topol, S. K. Burt, R. A. Caldwell, A. A. Rashin, J. Chem. Phys. 1998, 109, 4852 – 4863. [47] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, J. Chem. Phys. 1983, 79, 926 – 935. [48] AMBER 10, D. A. Case, T. A. Darden, T. E. Cheatham, III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, M. Crowley, R. C. Walker, W. Zhang, K. M. Merz, B. Wang, S. Hayik, A. Roitberg, G. Seabra, K. F. Kolossvary, K. F. Wong, F. Paesani, J. Vanicek, X. Wu, S. Brozell, T. Steinbrecher, H. Gohlke, L. Yang, C. Tan, J. Mongan, V. Hornak, G. Cui, D. H. Mathews, M. G. Seetin, C. Sagui, V. Babin, P. A. Kollman, University of California, San Francisco, CA, 2008. [49] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G. Pedersen, J. Chem. Phys. 1995, 103, 8577 – 8593. [50] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem. Phys. 1984, 81, 3684 – 3690. [51] M. Feig, J. Karanicolas, C. L. I. Brooks, J. Mol. Graphics Modell. 2004, 22, 377 – 395. [52] L. Goerigk, S. Grimme, J. Chem. Theory Comput. 2011, 7, 291 – 309. [53] J. P. Perdew, K. Burke, Y. Wang, Phys Rev B Condens. Matter Mater. Phys. 1996, 54, 16533 – 16539. [54] A. D. Becke, J. Chem. Phys. 1996, 104, 1040 – 1046. [55] L. Goerigk, S. Grimme, Phys. Chem. Chem. Phys. 2011, 13, 6670 – 6688. [56] F. Neese, WIREs Comput. Mol. Sci. 2012, 2, 73 – 78. [57] O. Vahtras, J. Almlçf, M. W. Feyereisen, Chem. Phys. Lett. 1993, 213, 514 – 518. [58] F. Weigend, A. Kçhn, C. H•ttig, J. Chem. Phys. 2002, 116, 3175 – 3183. Received: February 27, 2015 Published online on April 14, 2015

1656

Ó 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim