Development of CHARMM-Compatible Force-Field Parameters for Cobalamin and Related Cofactors from Quantum Mechanical Calculations

Corrinoid cofactors such as cobalamin are used by many enzymes and are essential for most living organisms. Therefore, there is broad interest in investigating cobalamin–protein interactions with molecular dynamics simulations. Previously developed parameters for cobalamins are based mainly on crystal structure data. Here, we report CHARMM-compatible force field parameters for several corrinoids developed from quantum mechanical calculations. We provide parameters for corrinoids in three oxidation states, Co³⁺, Co²⁺, and Co¹⁺, and with various axial ligands. Lennard-Jones parameters for the cobalt center in the Co­(II) and Co­(I) states were optimized using a helium atom probe, and partial atomic charges were obtained with a combination of natural population analysis (NPA) and restrained electrostatic potential (RESP) fitting approaches. The Force Field Toolkit was used to optimize all bonded terms. The resulting parameters, determined solely from calculations of cobalamin alone or in water, were then validated by assessing their agreement with density functional theory geometries and by analyzing molecular dynamics simulation trajectories of several corrinoid proteins for which X-ray crystal structures are available. In each case, we obtained excellent agreement with the reference data. In comparison to previous CHARMM-compatible parameters for cobalamin, we observe a better agreement for the fold angle and lower RMSD in the cobalamin binding site. The approach described here is readily adaptable for developing CHARMM-compatible force-field parameters for other corrinoids or large biomolecules

Toward Quantitatively Accurate Calculation of the Redox-Associated Acid–Base and Ligand Binding Equilibria of Aquacobalamin

Redox processes in complex transition metal-containing species are often intimately associated with changes in ligand protonation states and metal coordination number. A major challenge is therefore to develop consistent computational approaches for computing pH-dependent redox and ligand dissociation properties of organometallic species. Reduction of the Co center in the vitamin B12 derivative aquacobalamin can be accompanied by ligand dissociation, protonation, or both, making these properties difficult to compute accurately. We examine this challenge here by using density functional theory and continuum solvation to compute Co–ligand binding equilibrium constants (<i>K</i>_on/off), p<i>K</i>_as, and reduction potentials for models of aquacobalamin in aqueous solution. We consider two models for cobalamin ligand coordination: the first follows the hexa, penta, tetra coordination scheme for Co^III, Co^II, and Co^I species, respectively, and the second model features saturation of each vacant axial coordination site on Co^II and Co^I species with a single, explicit water molecule to maintain six directly interacting ligands or water molecules in each oxidation state. Comparing these two coordination schemes in combination with five dispersion-corrected density functionals, we find that the accuracy of the computed properties is largely independent of the scheme used, but including only a continuum representation of the solvent yields marginally better results than saturating the first solvation shell around Co throughout. PBE performs best, displaying balanced accuracy and superior performance overall, with RMS errors of 80 mV for seven reduction potentials, 2.0 log units for five p<i>K</i>_as and 2.3 log units for two log <i>K</i>_on/off values for the aquacobalamin system. Furthermore, we find that the BP86 functional commonly used in corrinoid studies suffers from erratic behavior and inaccurate descriptions of Co–axial ligand binding, leading to substantial errors in predicted p<i>K</i>_as and <i>K</i>_on/off values. These findings demonstrate the effectiveness of the present approach for computing electrochemical and thermodynamic properties of a complex transition metal-containing cofactor

Quantum Chemical Calculation of p<i>K</i>_as of Environmentally Relevant Functional Groups: Carboxylic Acids, Amines, and Thiols in Aqueous Solution

Developing accurate quantum chemical approaches for calculating p<i>K</i>_as is of broad interest. Useful accuracy can be obtained by using density functional theory (DFT) in combination with a polarizable continuum solvent model. However, some classes of molecules present problems for this approach, yielding errors greater than 5 p<i>K</i> units. Various methods have been developed to improve the accuracy of the combined strategy. These methods perform well but either do not generalize or introduce additional degrees of freedom, increasing the computational cost. The Solvation Model based on Density (SMD) has emerged as one of the most commonly used continuum solvent models. Nevertheless, for some classes of organic compounds, e.g., thiols, the p<i>K</i>_as calculated with the original SMD model show errors of 6–10 p<i>K</i> units, and we traced these errors to inaccuracies in the solvation free energies of the anions. To improve the accuracy of p<i>K</i>_as calculated with DFT and the SMD model, we developed a scaled solvent-accessible surface approach for constructing the solute–solvent boundary. By using a “direct” approach, in which all quantities are computed in the presence of the continuum solvent, the use of thermodynamic cycles is avoided. Furthermore, no explicit water molecules are required. Three benchmark data sets of experimentally measured p<i>K</i>_a values, including 28 carboxylic acids, 10 aliphatic amines, and 45 thiols, were used to assess the optimized SMD model, which we call SMD with a scaled solvent-accessible surface (SMD_sSAS). Of the methods tested, the M06-2X density functional approximation, 6-31+G­(d,p) basis set, and SMD_sSAS solvent model provided the most accurate p<i>K</i>_as for each set, yielding mean unsigned errors of 0.9, 0.4, and 0.5 p<i>K</i> units, respectively, for carboxylic acids, aliphatic amines, and thiols. This approach is therefore useful for efficiently calculating the p<i>K</i>_as of environmentally relevant functional groups

Hydrolysis of DFP and the Nerve Agent (<i>S</i>)‑Sarin by DFPase Proceeds along Two Different Reaction Pathways: Implications for Engineering Bioscavengers

Organophosphorus (OP) nerve agents such as (<i>S</i>)-sarin are among the most highly toxic compounds that have been synthesized. Engineering enzymes that catalyze the hydrolysis of nerve agents (“bioscavengers”) is an emerging prophylactic approach to diminish their toxic effects. Although its native function is not known, diisopropyl fluorophosphatase (DFPase) from Loligo vulgaris catalyzes the hydrolysis of OP compounds. Here, we investigate the mechanisms of diisopropylfluorophosphate (DFP) and (<i>S</i>)-sarin hydrolysis by DFPase with quantum mechanical/molecular mechanical umbrella sampling simulations. We find that the mechanism for hydrolysis of DFP involves nucleophilic attack by Asp229 on phosphorus to form a pentavalent intermediate. P–F bond dissociation then yields a phosphoacyl enzyme intermediate in the rate-limiting step. The simulations suggest that a water molecule, coordinated to the catalytic Ca²⁺, donates a proton to Asp121 and then attacks the tetrahedral phosphoacyl intermediate to liberate the diisopropylphosphate product. In contrast, the calculated free energy barrier for hydrolysis of (<i>S</i>)-sarin by the same mechanism is highly unfavorable, primarily because of the instability of the pentavalent phosphoenzyme species. Instead, simulations suggest that hydrolysis of (<i>S</i>)-sarin proceeds by a mechanism in which Asp229 could activate an intervening water molecule for nucleophilic attack on the substrate. These findings may lead to improved strategies for engineering DFPase and related six-bladed β-propeller folds for more efficient degradation of OP compounds

Exploring Covalent Allosteric Inhibition of Antigen 85C from <i>Mycobacterium tuberculosis</i> by Ebselen Derivatives

Previous studies identified ebselen as a potent <i>in vitro</i> and <i>in vivo</i> inhibitor of the <i>Mycobacterium tuberculosis</i> (<i>Mtb</i>) antigen 85 (Ag85) complex, comprising three homologous enzymes required for the biosynthesis of the mycobacterial cell wall. In this study, the <i>Mtb</i> Ag85C enzyme was cocrystallized with azido and adamantyl ebselen derivatives, resulting in two crystallographic structures of 2.01 and 1.30 Å resolution, respectively. Both structures displayed the anticipated covalent modification of the solvent accessible, noncatalytic Cys209 residue forming a selenenylsulfide bond. Continuous difference density for both thiol modifiers allowed for the assessment of interactions that influence ebselen binding and inhibitor orientation that were unobserved in previous Ag85C ebselen structures. The <i>k</i>_inact/<i>K</i>_I values for ebselen, adamantyl ebselen, and azido ebselen support the importance of observed constructive chemical interactions with Arg239 for increased <i>in vitro</i> efficacy toward Ag85C. To better understand the <i>in vitro</i> kinetic properties of these ebselen derivatives, the energetics of specific protein–inhibitor interactions and relative reaction free energies were calculated for ebselen and both derivatives using density functional theory. These studies further support the different <i>in vitro</i> properties of ebselen and two select ebselen derivatives from our previously published ebselen library with respect to kinetics and protein–inhibitor interactions. In both structures, the α9 helix was displaced farther from the enzyme active site than the previous Ag85C ebselen structure, resulting in the restructuring of a connecting loop and imparting a conformational change to residues believed to play a role in substrate binding specific to Ag85C. These notable structural changes directly affect protein stability, reducing the overall melting temperature by up to 14.5 °C, resulting in the unfolding of protein at physiological temperatures. Additionally, this structural rearrangement due to covalent allosteric modification creates a sizable solvent network that encompasses the active site and extends to the modified Cys209 residue. In all, this study outlines factors that influence enzyme inhibition by ebselen and its derivatives while further highlighting the effects of the covalent modification of Cys209 by said inhibitors on the structure and stability of Ag85C. Furthermore, the results suggest a strategy for developing new classes of Ag85 inhibitors with increased specificity and potency

Comparative Informatics Analysis to Evaluate Site-Specific Protein Oxidation in Multidimensional LC–MS/MS Data

Redox proteomics has yielded molecular insight into diseases of protein dysfunction attributable to oxidative stress, underscoring the need for robust detection of protein oxidation products. Additionally, oxidative protein surface mapping techniques utilize hydroxyl radicals to gain structural insight about solvent exposure. Interpretation of tandem mass spectral data is a critical challenge for such investigations, because reactive oxygen species target a wide breadth of amino acids. Additionally, oxidized peptides may be generated in a wide range of abundances since the reactivity of hydroxyl radicals with different amino acids spans 3 orders of magnitude. Taken together, these attributes of oxidative footprinting pose both experimental and computational challenges to detecting oxidized peptides that are naturally less abundant than their unoxidized counterparts. In this study, model proteins were oxidized electrochemically and analyzed at both the intact protein and peptide levels. A multidimensional chromatographic strategy was utilized to expand the dynamic range of oxidized peptide measurements. Peptide mass spectral data were searched by the “hybrid” software packages Inspect and Byonic, which incorporate <i>de novo</i> elements of spectral interpretation into a database search. This dynamic search capacity accommodates the challenge of searching for more than 40 oxidative mass shifts that can occur in a staggering variety of possible combinatorial occurrences. A prevailing set of oxidized residues was identified with this comparative approach, and evaluation of these sites was informed by solvent accessible surface area gleaned through molecular dynamics simulations. Along with increased levels of oxidation around highly reactive “hotspot” sites as expected, the enhanced sensitivity of these measurements uncovered a surprising level of oxidation on less reactive residues

Comparative Informatics Analysis to Evaluate Site-Specific Protein Oxidation in Multidimensional LC–MS/MS Data

Redox proteomics has yielded molecular insight into diseases of protein dysfunction attributable to oxidative stress, underscoring the need for robust detection of protein oxidation products. Additionally, oxidative protein surface mapping techniques utilize hydroxyl radicals to gain structural insight about solvent exposure. Interpretation of tandem mass spectral data is a critical challenge for such investigations, because reactive oxygen species target a wide breadth of amino acids. Additionally, oxidized peptides may be generated in a wide range of abundances since the reactivity of hydroxyl radicals with different amino acids spans 3 orders of magnitude. Taken together, these attributes of oxidative footprinting pose both experimental and computational challenges to detecting oxidized peptides that are naturally less abundant than their unoxidized counterparts. In this study, model proteins were oxidized electrochemically and analyzed at both the intact protein and peptide levels. A multidimensional chromatographic strategy was utilized to expand the dynamic range of oxidized peptide measurements. Peptide mass spectral data were searched by the “hybrid” software packages Inspect and Byonic, which incorporate <i>de novo</i> elements of spectral interpretation into a database search. This dynamic search capacity accommodates the challenge of searching for more than 40 oxidative mass shifts that can occur in a staggering variety of possible combinatorial occurrences. A prevailing set of oxidized residues was identified with this comparative approach, and evaluation of these sites was informed by solvent accessible surface area gleaned through molecular dynamics simulations. Along with increased levels of oxidation around highly reactive “hotspot” sites as expected, the enhanced sensitivity of these measurements uncovered a surprising level of oxidation on less reactive residues

Exploring Covalent Allosteric Inhibition of Antigen 85C from <i>Mycobacterium tuberculosis</i> by Ebselen Derivatives

Previous studies identified ebselen as a potent <i>in vitro</i> and <i>in vivo</i> inhibitor of the <i>Mycobacterium tuberculosis</i> (<i>Mtb</i>) antigen 85 (Ag85) complex, comprising three homologous enzymes required for the biosynthesis of the mycobacterial cell wall. In this study, the <i>Mtb</i> Ag85C enzyme was cocrystallized with azido and adamantyl ebselen derivatives, resulting in two crystallographic structures of 2.01 and 1.30 Å resolution, respectively. Both structures displayed the anticipated covalent modification of the solvent accessible, noncatalytic Cys209 residue forming a selenenylsulfide bond. Continuous difference density for both thiol modifiers allowed for the assessment of interactions that influence ebselen binding and inhibitor orientation that were unobserved in previous Ag85C ebselen structures. The <i>k</i>_inact/<i>K</i>_I values for ebselen, adamantyl ebselen, and azido ebselen support the importance of observed constructive chemical interactions with Arg239 for increased <i>in vitro</i> efficacy toward Ag85C. To better understand the <i>in vitro</i> kinetic properties of these ebselen derivatives, the energetics of specific protein–inhibitor interactions and relative reaction free energies were calculated for ebselen and both derivatives using density functional theory. These studies further support the different <i>in vitro</i> properties of ebselen and two select ebselen derivatives from our previously published ebselen library with respect to kinetics and protein–inhibitor interactions. In both structures, the α9 helix was displaced farther from the enzyme active site than the previous Ag85C ebselen structure, resulting in the restructuring of a connecting loop and imparting a conformational change to residues believed to play a role in substrate binding specific to Ag85C. These notable structural changes directly affect protein stability, reducing the overall melting temperature by up to 14.5 °C, resulting in the unfolding of protein at physiological temperatures. Additionally, this structural rearrangement due to covalent allosteric modification creates a sizable solvent network that encompasses the active site and extends to the modified Cys209 residue. In all, this study outlines factors that influence enzyme inhibition by ebselen and its derivatives while further highlighting the effects of the covalent modification of Cys209 by said inhibitors on the structure and stability of Ag85C. Furthermore, the results suggest a strategy for developing new classes of Ag85 inhibitors with increased specificity and potency

Cluster-Continuum Calculations of Hydration Free Energies of Anions and Group 12 Divalent Cations

Understanding aqueous phase processes involving group 12 metal cations is relevant to both environmental and biological sciences. Here, quantum chemical methods and polarizable continuum models are used to compute the hydration free energies of a series of divalent group 12 metal cations (Zn²⁺, Cd²⁺, and Hg²⁺) together with Cu²⁺ and the anions OH^–, SH^–, Cl^–, and F^–. A cluster-continuum method is employed, in which gas-phase clusters of the ion and explicit solvent molecules are immersed in a dielectric continuum. Two approaches to define the size of the solute–water cluster are compared, in which the number of explicit waters used is either held constant or determined variationally as that of the most favorable hydration free energy. Results obtained with various polarizable continuum models are also presented. Each leg of the relevant thermodynamic cycle is analyzed in detail to determine how different terms contribute to the observed mean signed error (MSE) and the standard deviation of the error (STDEV) between theory and experiment. The use of a constant number of water molecules for each set of ions is found to lead to predicted relative trends that benefit from error cancellation. Overall, the best results are obtained with MP2 and the Solvent Model D polarizable continuum model (SMD), with eight explicit water molecules for anions and 10 for the metal cations, yielding a STDEV of 2.3 kcal mol^–1 and MSE of 0.9 kcal mol^–1 between theoretical and experimental hydration free energies, which range from −72.4 kcal mol^–1 for SH^– to −505.9 kcal mol^–1 for Cu²⁺. Using B3PW91 with DFT-D3 dispersion corrections (B3PW91-D) and SMD yields a STDEV of 3.3 kcal mol^–1 and MSE of 1.6 kcal mol^–1, to which adding MP2 corrections from smaller divalent metal cation water molecule clusters yields very good agreement with the full MP2 results. Using B3PW91-D and SMD, with two explicit water molecules for anions and six for divalent metal cations, also yields reasonable agreement with experimental values, due in part to fortuitous error cancellation associated with the metal cations. Overall, the results indicate that the careful application of quantum chemical cluster-continuum methods provides valuable insight into aqueous ionic processes that depend on both local and long-range electrostatic interactions with the solvent

Cluster-Continuum Calculations of Hydration Free Energies of Anions and Group 12 Divalent Cations

Understanding aqueous phase processes involving group 12 metal cations is relevant to both environmental and biological sciences. Here, quantum chemical methods and polarizable continuum models are used to compute the hydration free energies of a series of divalent group 12 metal cations (Zn²⁺, Cd²⁺, and Hg²⁺) together with Cu²⁺ and the anions OH^–, SH^–, Cl^–, and F^–. A cluster-continuum method is employed, in which gas-phase clusters of the ion and explicit solvent molecules are immersed in a dielectric continuum. Two approaches to define the size of the solute–water cluster are compared, in which the number of explicit waters used is either held constant or determined variationally as that of the most favorable hydration free energy. Results obtained with various polarizable continuum models are also presented. Each leg of the relevant thermodynamic cycle is analyzed in detail to determine how different terms contribute to the observed mean signed error (MSE) and the standard deviation of the error (STDEV) between theory and experiment. The use of a constant number of water molecules for each set of ions is found to lead to predicted relative trends that benefit from error cancellation. Overall, the best results are obtained with MP2 and the Solvent Model D polarizable continuum model (SMD), with eight explicit water molecules for anions and 10 for the metal cations, yielding a STDEV of 2.3 kcal mol^–1 and MSE of 0.9 kcal mol^–1 between theoretical and experimental hydration free energies, which range from −72.4 kcal mol^–1 for SH^– to −505.9 kcal mol^–1 for Cu²⁺. Using B3PW91 with DFT-D3 dispersion corrections (B3PW91-D) and SMD yields a STDEV of 3.3 kcal mol^–1 and MSE of 1.6 kcal mol^–1, to which adding MP2 corrections from smaller divalent metal cation water molecule clusters yields very good agreement with the full MP2 results. Using B3PW91-D and SMD, with two explicit water molecules for anions and six for divalent metal cations, also yields reasonable agreement with experimental values, due in part to fortuitous error cancellation associated with the metal cations. Overall, the results indicate that the careful application of quantum chemical cluster-continuum methods provides valuable insight into aqueous ionic processes that depend on both local and long-range electrostatic interactions with the solvent