A combined semiempirical and DFT computational protocol for studying bioorganometallic complexes: Application to molybdocene–cysteine complexes

Abstract

Computational methods can help in the design of new bioorganometallic compounds. However, the presence of multihapto or σ/π metal-ligand bonding still precludes the direct application of either pure molecular mechanics (MM) or hybrid quantum mechanics-MM methods to study the flexibility of biomolecules in complex with organometallics. Herein, we present a computational protocol aimed to the evaluation of the relative free energies of bioorganometallic compounds, which explores the conformational space by means of Molecular Dynamics simulations using the semiempirical PM6 method coupled with the COnductor-like Screening MOdel solvation model followed by density functional theory (DFT) calculations including the DFT-D3 dispersion energy correction on the most stable conformers. This protocol is applied to investigate the complexes formed between cysteine and the molybdocene entity Cp2Mo(2+) (Cp=η(5)-C5H5), which has received considerable attention due to its potential antitumor activity and chemical reactivity. The calculated structures and free energies are in agreement with those of the experimentally detected molybdocene-cysteine adducts and allow us to investigate the reaction mechanism for the formation of the most stable 1:1 and 1:2 adducts.

Molecular Modeling of Bioorganometallic Compounds: Thermodynamic Properties of Molybdocene–Glutathione Complexes and Mechanism of Peptide Hydrolysis

“…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 [Cp 2 Mo(GSH) 1–2 ] species that we examined computationally are shown in Scheme , 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.…”

“…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 processes25–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.…”

“…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 γ‐Glu–Cys peptides in aqueous solution.…”

“…In previous work, we have examined the ability of the B3LYP functional complemented with the DFT‐D3 dispersion correction for studying mononuclear [Cp 2 Mo(L)(L′)] complexes, where L, L′=H 2 O, OH − , H 2 S, SH − and Cl − 25. 26 For these complexes, the B3LYP+DFT‐D3 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 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 protocol25 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.…”