Acidity constants and protonation sites of cyclic dinucleotides determined by capillary electrophoresis, quantum chemical calculations, and NMR spectroscopy.

Cyclic dinucleotides (CDNs) are important second messengers in bacteria and eukaryotes. Detailed characterization of their physicochemical properties is a prerequisite for understanding their biological functions. Herein, we examine acid-base and electromigration properties of selected CDNs employing capillary electrophoresis (CE), density functional theory (DFT), and nuclear magnetic resonance (NMR) spectroscopy to provide benchmark pKa values, as well as to unambiguously determine the protonation sites. Acidity constants (pKa ) of the NH+ moieties of adenine and guanine bases and actual and limiting ionic mobilities of CDNs were determined by nonlinear regression analysis of the pH dependence of their effective electrophoretic mobilities measured by CE in aqueous background electrolytes in a wide pH range (0.98-11.48), at constant temperature (25°C), and constant ionic strength (25 mM). The thermodynamic pKa values were found to be in the range 3.31-4.56 for adenine and 2.28-3.61 for guanine bases, whereas the pKa of enol group of guanine base was in the range 10.21-10.40. Except for systematic shifts of ∼2 pKa , the pKa values calculated by the DFT-D3//COSMO-RS composite protocol that included large-scale conformational sampling and "cross-morphing" were in a relatively good agreement with the pKa s determined by CE and predict N1 atom of adenine and N7 atom of guanine as the protonation sites. The protonation of the N1 atom of adenine and N7 atom of guanine in acidic background electrolytes (BGEs) and the dissociation of the enol group of guanine in alkaline BGEs was confirmed also by NMR spectroscopy.

Predicting Effects of Site-Directed Mutagenesis on Enzyme Kinetics by QM/MM and QM Calculations: A Case of Glutamate Carboxypeptidase II.

Quantum and molecular mechanics (QM/MM) and QM-only (cluster model) modeling techniques represent the two workhorses in mechanistic understanding of enzyme catalysis. One of the stringent tests for QM/MM and/or QM approaches is to provide quantitative answers to real-world biochemical questions, such as the effect of single-point mutations on enzyme kinetics. This translates into predicting the relative activation energies to 1-2 kcal·mol-1 accuracy; such predictions can be used for the rational design of novel enzyme variants with desired/improved characteristics. Herein, we employ glutamate carboxypeptidase II (GCPII), a dizinc metallopeptidase, also known as the prostate specific membrane antigen, as a model system. The structure and activity of this major cancer antigen have been thoroughly studied, both experimentally and computationally, which makes it an ideal model system for method development. Its reaction mechanism is quite well understood: the reaction coordinate comprises a "tetrahedral intermediate" and two transition states and experimental activation Gibbs free energy of ∼17.5 kcal·mol-1 can be inferred for the known kcat ≈ 1 s-1. We correlate experimental kinetic data (including the E424H variant, newly characterized in this work) for various GCPII mutants (kcat = 8.6 × 10-5 s-1 to 2.7 s-1) with the energy profiles calculated by QM/MM and QM-only (cluster model) approaches. We show that the near-quantitative agreement between the experimental values and the calculated activation energies (ΔH⧧) can be obtained and recommend the combination of the two protocols: QM/MM optimized structures and cluster model (QM) energetics. The trend in relative activation energies is mostly independent of the QM method (DFT functional) used. Last but not least, a satisfactory correlation between experimental and theoretical data allows us to provide qualitative and fairly simple explanations of the observed kinetic effects which are thus based on a rigorous footing.

Enzymatic Synthesis of 3'-5', 3'-5' Cyclic Dinucleotides, Their Binding Properties to the Stimulator of Interferon Genes Adaptor Protein, and Structure/Activity Correlations.

The 3'-5', 3'-5' cyclic dinucleotides (3'3'CDNs) are bacterial second messengers that can also bind to the stimulator of interferon genes (STING) adaptor protein in vertebrates and activate the host innate immunity. Here, we profiled the substrate specificity of four bacterial dinucleotide synthases from Vibrio cholerae (DncV), Bacillus thuringiensis (btDisA), Escherichia coli (dgcZ), and Thermotoga maritima (tDGC) using a library of 33 nucleoside-5'-triphosphate analogues and then employed these enzymes to synthesize 24 3'3'CDNs. The STING affinity of CDNs was evaluated in cell-based and biochemical assays, and their ability to induce cytokines was determined by employing human peripheral blood mononuclear cells. Interestingly, the prepared heterodimeric 3'3'CDNs bound to the STING much better than their homodimeric counterparts and showed similar or better potency than bacterial 3'3'CDNs. We also rationalized the experimental findings by in-depth STING-CDN structure-activity correlations by dissecting computed interaction free energies into a set of well-defined and intuitive terms. To this aim, we employed state-of-the-art methods of computational chemistry, such as quantum mechanics/molecular mechanics (QM/MM) calculations, and complemented the computed results with the {STING:3'3'c-di-ara-AMP} X-ray crystallographic structure. QM/MM identified three outliers (mostly homodimers) for which we have no clear explanation of their impaired binding with respect to their heterodimeric counterparts, whereas the R2 = 0.7 correlation between the computed ΔG'int_rel and experimental ΔTm's for the remaining ligands has been very encouraging.

Synthesis and Biological Evaluation of Phosphoester and Phosphorothioate Prodrugs of STING Agonist 3',3'-c-Di(2'F,2'dAMP).

Cyclic dinucleotides (CDNs) are second messengers that bind to the stimulator of interferon genes (STING) and trigger the expression of type I interferons and proinflammatory cytokines. Here we evaluate the activity of 3',3'-c-di(2'F,2'dAMP) and its phosphorothioate analogues against five STING allelic forms in reporter-cell-based assays and rationalize our findings with X-ray crystallography and quantum mechanics/molecular mechanics calculations. We show that the presence of fluorine in the 2' position of 3',3'-c-di(2'F,2'dAMP) improves its activity not only against the wild type (WT) but also against REF and Q STING. Additionally, we describe the synthesis of the acyloxymethyl and isopropyloxycarbonyl phosphoester prodrugs of CDNs. Masking the negative charges of the CDNs results in an up to a 1000-fold improvement of the activities of the prodrugs relative to those of their parent CDNs. Finally, the uptake and intracellular cleavage of pivaloyloxymethyl prodrugs to the parent CDN is rapid, reaching a peak intracellular concentration within 2 h.

Conformational energies and equilibria of cyclic dinucleotides in vacuo and in solution: computational chemistry vs. NMR experiments.

Performance of computational methods in modelling cyclic dinucleotides - an important and challenging class of compounds - has been evaluated by two different benchmarks: (1) gas-phase conformational energies and (2) qualitative agreement with NMR observations of the orientation of the χ-dihedral angle in solvent. In gas-phase benchmarks, where CCSD(T) and DLPNO-CCSD(T) methods have been used as the reference, most of the (dispersion corrected) density functional approximations are accurate enough to justify prioritizing computational cost and compatibility with other modelling options as the criterion of choice. NMR experiments of 3'3'-c-di-AMP, 3'3'-c-GAMP, and 3'3'-c-di-GMP show the overall prevalence of the anti-conformation of purine bases, but some population of syn-conformations is observed for guanines. Implicit solvation models combined with quantum-chemical methods struggle to reproduce this behaviour, probably due to a lack of dynamics and explicitly modelled solvent, leading to structures that are too compact. Molecular dynamics simulations overrepresent the syn-conformation of guanine due to the overestimation of an intramolecular hydrogen bond. Our combination of experimental and computational benchmarks provides "error bars" for modelling cyclic dinucleotides in solvent, where such information is generally difficult to obtain, and should help gauge the interpretability of studies dealing with binding of cyclic dinucleotides to important pharmaceutical targets. At the same time, the presented analysis calls for improvement in both implicit solvation models and force-field parameters.

Ligand Strain and Its Conformational Complexity Is a Major Factor in the Binding of Cyclic Dinucleotides to STING Protein.

STING (stimulator of interferon genes) is a key regulator of innate immunity that has recently been recognized as a promising drug target. STING is activated by cyclic dinucleotides (CDNs) which eventually leads to expression of type I interferons and other cytokines. Factors underlying the affinity of various CDN analogues are poorly understood. Herein, we correlate structural biology, isothermal calorimetry (ITC) and computational modeling to elucidate factors contributing to binding of six CDNs-three pairs of natural (ribo) and fluorinated (2'-fluororibo) 3',3'-CDNs. X-ray structural analyses of six {STING:CDN} complexes did not offer any explanation for the different affinities of the studied ligands. ITC showed entropy/enthalpy compensation up to 25 kcal mol-1 for this set of similar ligands. The higher affinities of fluorinated analogues are explained with help of computational methods by smaller loss of entropy upon binding and by smaller strain (free) energy.

Protein-Ligand Interactions in the STING Binding Site Probed by Rationally Designed Single-Point Mutations: Experiment and Theory.

STING protein (stimulator of interferon genes) plays an important role in the innate immune system. A number of potent compounds regulating its activity have been reported, mostly derivatives of cyclic dinucleotides (CDNs), natural STING agonists. Here, we aim to provide complementary information to large-scale "ligand-profiling" studies by probing the importance of STING-CDN protein-ligand interactions on the protein side. We examined in detail six typical CDNs each in complex with 13 rationally devised mutations in STING: S162A, S162T, Y167F, G230A, R232K, R232H, A233L, A233I, R238K, T263A, T263S, R293Q, and G230A/R293Q. The mutations switch on and off various types of protein-ligand interactions: π-π stacking, hydrogen bonding, ionic pairing, and nonpolar contacts. We correlated experimental data obtained by differential scanning fluorimetry, X-ray crystallography, and isothermal titration calorimetry with theoretical calculations. This enabled us to provide a mechanistic interpretation of the differences in the binding of representative CDNs to STING. We observed that the G230A mutation increased the thermal stability of the protein-ligand complex, indicating an increased level of ligand binding, whereas R238K and Y167F led to a complete loss of stabilization (ligand binding). The effects of the other mutations depended on the type of ligand (CDN) and varied, to some extent. A very good correlation (R2 = 0.6) between the experimental binding affinities and interaction energies computed by quantum chemical methods enabled us to explain the effect of the studied mutations in detail and evaluate specific interactions quantitatively. Our work may inspire development of high-affinity ligands against the common STING haplotypes by targeting the key (sometimes non-intuitive) protein-ligand interactions.

Mapping Conformational Space of All 8000 Tripeptides by Quantum Chemical Methods: What Strain Is Affordable within Folded Protein Chains?

To gain more insight into the physicochemical aspects of a protein structure from the first principles, conformational space of all 8000 "capped" tripeptides (i.e., N-Ac-X1X2X3-NH-CH3, where Xi is one of the 20 natural amino acids) was investigated computationally. An enormous dataset (denoted P-CONF_1.6M and containing close to 1 600 000 conformers in total) has been obtained by employing a composite protocol combining density functional theory, semiempirical quantum mechanics (SQM), and state-of-the-art solvation methods with 1000 K molecular dynamics (MD) used to generate initial structures (200 snapshots for each tripeptide). This allowed us to present the first rigorous QM-based glimpse at the vast conformational space spanned by small protein fragments. The same computational procedure was repeated for tripeptide fragments taken from the SCOPe database of three-dimensional protein folds, by restraining them to their geometry in a protein. Such complementary data allowed us to compare the distribution of conformational strain energies of unrestrained tripeptidic fragments "in solvent" with those in existing protein chains. Besides providing a rigorous (ab initio) proof of a few well-known concepts and hypotheses concerning protein structures, such as the distribution of (φ, ψ) angles in Ramachandran plots, we have made several observations that came as a certain surprise: (1) distribution of conformational energies does not significantly differ between the "unbiased/unrestrained" conformers obtained from MD sampling in solvent and the biased conformers, i.e., those of a given tripeptide obtained from protein structures; (2) conformational (strain) energy window up to ∼20 to 25 kcal·mol-1 is readily available to tripeptide fragments within the context of a protein chain; (3) overpopulation in certain regions of Ramachandran plot was observed for the unbiased conformers. Last but not least, the massive dataset of accurate (DFT-D3//COSMO-RS) conformational (free) energies of ∼1.6 M peptide conformers, P-CONF_1.6M, obtained throughout this work may serve as excellent dataset for calibrating and benchmarking of popular force fields.

Enzymatic Preparation of 2'-5',3'-5'-Cyclic Dinucleotides, Their Binding Properties to Stimulator of Interferon Genes Adaptor Protein, and Structure/Activity Correlations.

Cyclic dinucleotides are second messengers in the cyclic GMP-AMP synthase (cGAS)-stimulator of interferon genes (STING) pathway, which plays an important role in recognizing tumor cells and viral or bacterial infections. They bind to the STING adaptor protein and trigger expression of cytokines via TANK binding kinase 1 (TBK1)/interferon regulatory factor 3 (IRF3) and inhibitor of nuclear factor-κB (IκB) kinase (IKK)/nuclear factor-κB (NFκB) signaling cascades. In this work, we describe an enzymatic preparation of 2'-5',3'-5'-cyclic dinucleotides (2'3'CDNs) with use of cyclic GMP-AMP synthases (cGAS) from human, mouse, and chicken. We profile substrate specificity of these enzymes by employing a small library of nucleotide-5'-triphosphate (NTP) analogues and use them to prepare 33 2'3'CDNs. We also determine affinity of these CDNs to five different STING haplotypes in cell-based and biochemical assays and describe properties needed for their optimal activity toward all STING haplotypes. Next, we study their effect on cytokine and chemokine induction by human peripheral blood mononuclear cells (PBMCs) and evaluate their cytotoxic effect on monocytes. Additionally, we report X-ray crystal structures of two new CDNs bound to STING protein and discuss structure-activity relationship by using quantum and molecular mechanical (QM/MM) computational modeling.

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).

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 macrocyclic 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 (GDFT/COSMO-RS) for all of the conformers of each compound were obtained by a composite protocol based on in vacuo quantum mechanics (DFT-D3 method in a large basis set), standard gas-phase thermodynamics, and the COSMO-RS solvation model. The GDFT/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 GDFT/COSMO-RS) than the latter (i.e., plain SQM MD). The CPU time needed to fully evaluate one medium-sized compound (∼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.

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), -CH2PO2HCH2PO2H2 (H2L(1)), or phosphinic-phosphonic acid, -CH2PO2HCH2PO3H2 (H3L(2)), 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 KCuL = 25.8 and 27.7 for H2L(1) and H3L(2), respectively). Kinetic study revealed an unusual three-step complex formation mechanism. The initial equilibrium step leads to out-of-cage complexes with Cu(2+) bound by the phosphorus-containing pendant arm. These species quickly rearrange to an in-cage complex with cyclam conformation II, which isomerizes to another in-cage complex with cyclam conformation I. The first in-cage complex is quantitatively formed in seconds (pH ≈5, 25 °C, Cu:L = 1:1, cM ≈ 1 mM). At pH >12, I isomers undergo nitrogen atom inversion, leading to III isomers; the structure of the III-[Cu(HL(2))] complex in the solid state was confirmed by X-ray diffraction analysis. In an alkaline solution, interconversion of the I and III 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 H3L(2) 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 H2L(1) and H3L(2) with (64)Cu was fast and efficient, even at room temperature, giving specific activities of around 70 GBq of (64)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 H3nota and H4dota 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.

How simple is too simple? Computational perspective on importance of second-shell environment for metal-ion selectivity.

The metal-ion selectivity in biomolecules represents one of the most important phenomena in bioinorganic chemistry. The open question to what extent is the selectivity in the complex bioinorganic structures such as metallopeptides determined by the first-shell ligands of the metal ion is answered herein using six model peptides complexed with the set of divalent metal ions (Mn(2+), Fe(2+), Co(2+), Ni(2+), Cu(2+), Zn(2+), Cd(2+), and Hg(2+)) and their various first-shell representations. By calculating the differences among the free energies of complexation of metal ions in these peptides and their model (truncated) systems it is quantitatively shown that the definition of the first shell is paramount to this discussion and revolves around the chemical nature of the binding site. Despite the vast conceivable diversity of peptidic structures, that suggest certain fluidity of this definition, major contributing factors are identified and assessed based on their importance for capturing metal-ion selectivity. These factors include soft/hard character of ligands and various non-covalent interactions in the vicinity of the binding site. The relative importance of these factors is considered and specific suggestions for effective construction of the models are made. The relationship of first-shell models and their corresponding parent peptides is discussed thoroughly, both with respect to their chemical similarity and potential disparity introduced by generally "non-alignable" conformational flexibility of the two systems. It is concluded that, in special cases, this disparity can be negligible and that heeding the chemical factors contributing to selectivity during construction of the model can successfully result in models that retain the affinity profile for various metal ions with high fidelity.

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(NH3)4](2+), [M(NH3)4(H2O)2](2+), [M(Imi)(H2O)5](2+), [M(H2O)3(His)](+), [M(H2O)4(Cys)], [M(H2O)3(Cys)], [M(CH3COO)(H2O)3](+), [M(CH3COO)(H2O)5](+), [M(SCH2COO)2](2-)--with eight divalent metal ions (Mn(2+), Fe(2+), Co(2+), Ni(2+), Cu(2+), Zn(2+), Cd(2+), and Hg(2+)). Using the currently available computational protocols, we show that it is possible to achieve a relative 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.

Interaction of metal ions with biomolecular ligands: how accurate are calculated free energies associated with metal ion complexation?

To address fundamental questions in bioinorganic chemistry, such as metal ion selectivity, accurate computational protocols for both the gas-phase association of metal-ligand complexes and solvation/desolvation energies of the species involved are needed. In this work, we attempt to critically evaluate the performance of the ab initio and DFT electronic structure methods available and recent solvation models in calculations of the energetics associated with metal ion complexation. On the example of five model complexes ([M(II)(CH(3)S)(H(2)O)](+), [M(II)(H(2)O)(2)(H(2)S)(NH(3))](2+), [M(II)(CH(3)S)(NH(3))(H(2)O)(CH(3)COO)], [M(II)(H(2)O)(3)(SH)(CH(3)COO)(Im)], [M(II)(H(2)S)(H(2)O)(CH(3)COO)(PhOH)(Im)](+) in typical coordination geometries) and four metal ions (Fe(2+), Cu(2+), Zn(2+), and Cd(2+); representing open- and closed-shell and the first- and second-row transition metal elements), we provide reference values for the gas-phase complexation energies, as presumably obtained using the CCSD(T)/aug-cc-pVTZ method, and compare them with cheaper methods, such as DFT and RI-MP2, that can be used for large-scale calculations. We also discuss two possible definitions of interaction energies underlying the theoretically predicted metal-ion selectivity and the effect of geometry optimization on these values. Finally, popular solvation models, such as COSMO-RS and SMD, are used to demonstrate whether quantum chemical calculations can provide the overall free enthalpy (ΔG) changes in the range of the expected experimental values for the model complexes or match the experimental stability constants in the case of three complexes for which the experimental data exist. The data presented highlight several intricacies in the theoretical predictions of the experimental stability constants: the covalent character of some metal-ligand bonds (e.g., Cu(II)-thiolate) causing larger errors in the gas-phase complexation energies, inaccuracies in the treatment of solvation of the charged species, and difficulties in the definition of the reference state for Jahn-Teller unstable systems (e.g., [Cu(H(2)O)(6)](2+)). Although the agreement between the experimental (as derived from the stability constants) and calculated values is often within 5 kcal·mol(-1), in more complicated cases, it may exceed 15 kcal·mol(-1). Therefore, extreme caution must be exercised in assessing the subtle issues of metal ion selectivity quantitatively.