Predicting the Stability Constants of Metal-Ion Complexes from First Principles

The most important experimental quantity describing the thermodynamics of metal-ion binding with various (in)­organic ligands, or biomolecules, is the stability constant of the complex (β). In principle, it can be calculated as the free-energy change associated with the metal-ion complexation, i.e., its uptake from the solution under standard conditions. Because this process is associated with the interactions of charged species, large values of interaction and solvation energies are in general involved. Using the standard thermodynamic cycle (in vacuo complexation and solvation/desolvation of the reference state and of the resulting complexes), one usually subtracts values of several hundreds of kilocalories per mole to obtain final results on the order of units or tens of kilocalories per mole. In this work, we use density functional theory and Møller–Plesset second-order perturbation theory calculations together with the conductor-like screening model for realistic solvation to calculate the stability constants of selected complexes[M­(NH₃)₄]²⁺, [M­(NH₃)₄(H₂O)₂]²⁺, [M­(Imi)­(H₂O)₅]²⁺, [M­(H₂O)₃(His)]⁺, [M­(H₂O)₄(Cys)], [M­(H₂O)₃(Cys)], [M­(CH₃COO)­(H₂O)₃]⁺, [M­(CH₃COO)­(H₂O)₅]⁺, [M­(SCH₂COO)₂]^2–with eight divalent metal ions (Mn²⁺, Fe²⁺, Co²⁺, Ni²⁺, Cu²⁺, Zn²⁺, Cd²⁺, and Hg²⁺). Using the currently available computational protocols, we show that it is possible to achieve a <i>relative</i> accuracy of 2–4 kcal·mol^–1 (1–3 orders of magnitude in β). However, because most of the computed values are affected by metal- and ligand-dependent systematic shifts, the accuracy of the “absolute” (uncorrected) values is generally lower. For metal-dependent systematic shifts, we propose the specific values to be used for the given metal ion and current protocol. At the same time, we argue that ligand-dependent shifts (which cannot be easily removed) do not influence the metal-ion selectivity of the particular site, and therefore it can be computed to within 2 kcal·mol^–1 average accuracy. Finally, a critical discussion is presented that aims at potential caveats that one may encounter in theoretical predictions of the stability constants and highlights the perspective that theoretical calculations may become both competitive and complementary tools to experimental measurements

Predicting the Stability Constants of Metal-Ion Complexes from First Principles

The most important experimental quantity describing the thermodynamics of metal-ion binding with various (in)­organic ligands, or biomolecules, is the stability constant of the complex (β). In principle, it can be calculated as the free-energy change associated with the metal-ion complexation, i.e., its uptake from the solution under standard conditions. Because this process is associated with the interactions of charged species, large values of interaction and solvation energies are in general involved. Using the standard thermodynamic cycle (in vacuo complexation and solvation/desolvation of the reference state and of the resulting complexes), one usually subtracts values of several hundreds of kilocalories per mole to obtain final results on the order of units or tens of kilocalories per mole. In this work, we use density functional theory and Møller–Plesset second-order perturbation theory calculations together with the conductor-like screening model for realistic solvation to calculate the stability constants of selected complexes[M­(NH₃)₄]²⁺, [M­(NH₃)₄(H₂O)₂]²⁺, [M­(Imi)­(H₂O)₅]²⁺, [M­(H₂O)₃(His)]⁺, [M­(H₂O)₄(Cys)], [M­(H₂O)₃(Cys)], [M­(CH₃COO)­(H₂O)₃]⁺, [M­(CH₃COO)­(H₂O)₅]⁺, [M­(SCH₂COO)₂]^2–with eight divalent metal ions (Mn²⁺, Fe²⁺, Co²⁺, Ni²⁺, Cu²⁺, Zn²⁺, Cd²⁺, and Hg²⁺). Using the currently available computational protocols, we show that it is possible to achieve a <i>relative</i> accuracy of 2–4 kcal·mol^–1 (1–3 orders of magnitude in β). However, because most of the computed values are affected by metal- and ligand-dependent systematic shifts, the accuracy of the “absolute” (uncorrected) values is generally lower. For metal-dependent systematic shifts, we propose the specific values to be used for the given metal ion and current protocol. At the same time, we argue that ligand-dependent shifts (which cannot be easily removed) do not influence the metal-ion selectivity of the particular site, and therefore it can be computed to within 2 kcal·mol^–1 average accuracy. Finally, a critical discussion is presented that aims at potential caveats that one may encounter in theoretical predictions of the stability constants and highlights the perspective that theoretical calculations may become both competitive and complementary tools to experimental measurements

Macrocycle Conformational Sampling by DFT-D3/COSMO-RS Methodology

To find and calibrate a robust and reliable computational protocol for mapping conformational space of medium-sized molecules, exhaustive conformational sampling has been carried out for a series of seven <i>macrocyclic</i> compounds of varying ring size and one acyclic analogue. While five of them were taken from the MD/LLMOD/force field study by Shelley and co-workers (Watts, K. S.; Dalal, P.; Tebben, A. J.; Cheney, D. L.; Shelley, J. C. Macrocycle Conformational Sampling with MacroModel. J. Chem. Inf. Model. 2014, 54, 2680−2696), three represent potential macrocyclic inhibitors of human cyclophilin A. The free energy values (<i>G</i>_{DFT/COSMO‑RS}) for all of the conformers of each compound were obtained by a composite protocol based on <i>in vacuo</i> quantum mechanics (DFT-D3 method in a large basis set), standard gas-phase thermodynamics, and the COSMO-RS solvation model. The <i>G</i>_{DFT/COSMO‑RS} values were used as the reference for evaluating the performance of conformational sampling algorithms: standard and extended MD/LLMOD search (simulated-annealing molecular dynamics with low-lying eigenvector following algorithms, employing the OPLS2005 force field plus GBSA solvation) available in MacroModel and plain molecular dynamics (MD) sampling at high temperature (1000 K) using the semiempirical quantum mechanics (SQM) potential SQM­(PM6-D3H4/COSMO) followed by energy minimization of the snapshots. It has been shown that the former protocol (MD/LLMOD) may provide a more complete set of initial structures that ultimately leads to the identification of a greater number of low-energy conformers (as assessed by <i>G</i>_{DFT/COSMO‑RS}) than the latter (i.e., plain SQM MD). The CPU time needed to fully evaluate one medium-sized <i>compound</i> (∼100 atoms, typically resulting in several hundred or a few thousand conformers generated and treated quantum-mechanically) is approximately 1,000–100,000 CPU hours on today’s computers, which transforms to 1–7 days on a small-sized computer cluster with a few hundred CPUs. Finally, our data sets based on the rigorous quantum-chemical approach allow us to formulate a composite conformational sampling protocol with multiple checkpoints and truncation of redundant structural data that offers superior insights at affordable computational cost

Toward Accurate Conformational Energies of Smaller Peptides and Medium-Sized Macrocycles: MPCONF196 Benchmark Energy Data Set

A carefully selected set of acyclic and cyclic model peptides and several other macrocycles, comprising 13 compounds in total, has been used to calibrate the accuracy of the DFT­(-D3) method for conformational energies, employing BP86, PBE0, PBE, B3LYP, BLYP, TPSS, TPSSh, M06-2X, B97-D, OLYP, revPBE, M06-L, SCAN, revTPSS, BH-LYP, and ωB97X-D3 functionals. Both high- and low-energy conformers, 15 or 16 for each compound adding to 196 in total, denoted as the MPCONF196 data set, were included, and the reference values were obtained by the composite protocol, yielding the CCSD­(T)/​CBS extrapolated energies or their DLPNO-CCSD­(T)/​CBS equivalents in the case of larger systems. The latter was shown to be in near-quantitative (∼0.10–0.15 kcal·mol^–1) agreement with the canonical CCSD­(T), provided the TightPNO setting is used, and, therefore, can be used as the reference for larger systems (likely up to 150–200 atoms) for the problem studied here. At the same time, it was found that many D3-corrected DFT functionals provide results of ∼1 kcal·mol^–1 accuracy, which we consider as quite encouraging. This result implies that DFT-D3 methods can be, for example, safely used in efficient conformational sampling algorithms. Specifically, the DFT-D3/​DZVP-DFT level of calculation seems to be the best trade-off between computational cost and accuracy. Based on the calculated data, we have not found any cheaper variant for the treatment of conformational energies, since the semiempirical methods (including DFTB) provide results of inferior accuracy (errors of 3–5 kcal·mol^–1)

Toward Accurate Conformational Energies of Smaller Peptides and Medium-Sized Macrocycles: MPCONF196 Benchmark Energy Data Set

A carefully selected set of acyclic and cyclic model peptides and several other macrocycles, comprising 13 compounds in total, has been used to calibrate the accuracy of the DFT­(-D3) method for conformational energies, employing BP86, PBE0, PBE, B3LYP, BLYP, TPSS, TPSSh, M06-2X, B97-D, OLYP, revPBE, M06-L, SCAN, revTPSS, BH-LYP, and ωB97X-D3 functionals. Both high- and low-energy conformers, 15 or 16 for each compound adding to 196 in total, denoted as the MPCONF196 data set, were included, and the reference values were obtained by the composite protocol, yielding the CCSD­(T)/​CBS extrapolated energies or their DLPNO-CCSD­(T)/​CBS equivalents in the case of larger systems. The latter was shown to be in near-quantitative (∼0.10–0.15 kcal·mol^–1) agreement with the canonical CCSD­(T), provided the TightPNO setting is used, and, therefore, can be used as the reference for larger systems (likely up to 150–200 atoms) for the problem studied here. At the same time, it was found that many D3-corrected DFT functionals provide results of ∼1 kcal·mol^–1 accuracy, which we consider as quite encouraging. This result implies that DFT-D3 methods can be, for example, safely used in efficient conformational sampling algorithms. Specifically, the DFT-D3/​DZVP-DFT level of calculation seems to be the best trade-off between computational cost and accuracy. Based on the calculated data, we have not found any cheaper variant for the treatment of conformational energies, since the semiempirical methods (including DFTB) provide results of inferior accuracy (errors of 3–5 kcal·mol^–1)

Cyclam Derivatives with a Bis(phosphinate) or a Phosphinato–Phosphonate Pendant Arm: Ligands for Fast and Efficient Copper(II) Complexation for Nuclear Medical Applications

Cyclam derivatives bearing one geminal bis­(phosphinic acid), −CH₂PO₂HCH₂PO₂H₂ (H₂<b>L</b>^<b>1</b>), or phosphinic–phosphonic acid, −CH₂PO₂HCH₂PO₃H₂ (H₃<b>L</b>^<b>2</b>), pendant arm were synthesized and studied as potential copper­(II) chelators for nuclear medical applications. The ligands showed good selectivity for copper­(II) over zinc­(II) and nickel­(II) ions (log <i>K</i>_CuL = 25.8 and 27.7 for H₂<b>L</b>^<b>1</b> and H₃<b>L</b>^<b>2</b>, respectively). Kinetic study revealed an unusual three-step complex formation mechanism. The initial equilibrium step leads to <i>out-of-cage</i> complexes with Cu²⁺ bound by the phosphorus-containing pendant arm. These species quickly rearrange to an <i>in-cage</i> complex with cyclam conformation <b>II</b>, which isomerizes to another <i>in-cage</i> complex with cyclam conformation <b>I</b>. The first <i>in-cage</i> complex is quantitatively formed in seconds (pH ≈5, 25 °C, Cu:L = 1:1, <i>c</i>_M ≈ 1 mM). At pH >12, <b>I</b> isomers undergo nitrogen atom inversion, leading to <b>III</b> isomers; the structure of the <b>III</b>-[Cu­(H<b>L</b>^<b>2</b>)] complex in the solid state was confirmed by X-ray diffraction analysis. In an alkaline solution, interconversion of the <b>I</b> and <b>III</b> isomers is mutual, leading to the same equilibrium isomeric mixture; such behavior has been observed here for the first time for copper­(II) complexes of cyclam derivatives. Quantum-chemical calculations showed small energetic differences between the isomeric complexes of H₃<b>L</b>^<b>2</b> compared with analogous data for isomeric complexes of cyclam derivatives with one or two methylphosphonic acid pendant arm(s). Acid-assisted dissociation proved the kinetic inertness of the complexes. Preliminary radiolabeling of H₂<b>L</b>^<b>1</b> and H₃<b>L</b>^<b>2</b> with ⁶⁴Cu was fast and efficient, even at room temperature, giving specific activities of around 70 GBq of ⁶⁴Cu per 1 μmol of the ligand (pH 6.2, 10 min, ca. 90 equiv of the ligand). These specific activities were much higher than those of H₃<b>nota</b> and H₄<b>dota</b> complexes prepared under identical conditions. The rare combination of simple ligand synthesis, very fast copper­(II) complex formation, high thermodynamic stability, kinetic inertness, efficient radiolabeling, and expected low bone tissue affinity makes such ligands suitably predisposed to serve as chelators of copper radioisotopes in nuclear medicine