Showing papers by "David L. Mobley published in 2019"

Best Practices for Foundations in Molecular Simulations [Article v1.0].

Abstract: This document provides a starting point for approaching molecular simulations, guiding beginning practitioners to what issues they need to know about before and while starting their first simulations, and why those issues are so critical. This document makes no claims to provide an adequate introduction to the subject on its own. Instead, our goal is to help people know what issues are critical before beginning, and to provide references to good resources on those topics. We also provide a checklist of key issues to consider before and while setting up molecular simulations which may serve as a foundation for other best practices documents.

Best Practices for Alchemical Free Energy Calculations [Article v1.0].

Abstract: Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Toward Learned Chemical Perception of Force Field Typing Rules.

Abstract: Molecular mechanics force fields define how the energy and forces in a molecular system are computed from its atomic positions, thus enabling the study of such systems through computational methods...

Comparison of affinity ranking using AutoDock-GPU and MM-GBSA scores for BACE-1 inhibitors in the D3R Grand Challenge 4

Abstract: Molecular docking has been successfully used in computer-aided molecular design projects for the identification of ligand poses within protein binding sites. However, relying on docking scores to rank different ligands with respect to their experimental affinities might not be sufficient. It is believed that the binding scores calculated using molecular mechanics combined with the Poisson–Boltzman surface area (MM-PBSA) or generalized Born surface area (MM-GBSA) can predict binding affinities more accurately. In this perspective, we decided to take part in Stage 2 of the Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) to compare the performance of a quick scoring function, AutoDock4, to that of MM-GBSA in predicting the binding affinities of a set of $$\beta$$-Amyloid Cleaving Enzyme 1 (BACE-1) ligands. Our results show that re-scoring docking poses using MM-GBSA did not improve the correlation with experimental affinities. We further did a retrospective analysis of the results and found that our MM-GBSA protocol is sensitive to details in the protein-ligand system: (i) neutral ligands are more adapted to MM-GBSA calculations than charged ligands, (ii) predicted binding affinities depend on the initial conformation of the BACE-1 receptor, (iii) protonating the aspartyl dyad of BACE-1 correctly results in more accurate binding affinity predictions.

Biomolecular Solvation Structure Revealed by Molecular Dynamics Simulations.

Abstract: To compare ordered water positions from experiment with those from molecular dynamics (MD) simulations, a number of MD models of water structure in crystalline endoglucanase were calculated. The starting MD model was derived from a joint X-ray and neutron diffraction crystal structure, enabling the use of experimentally assigned protonation states. Simulations were performed in the crystalline state, using a periodic 2 × 2 × 2 supercell with explicit solvent. Water X-ray and neutron scattering density maps were computed from MD trajectories using standard macromolecular crystallography methods. In one set of simulations, harmonic restraints were applied to bias the protein structure toward the crystal structure. For these simulations, the recall of crystallographic waters using strong peaks in the MD water electron density was very good, and there also was substantial visual agreement between the boomerang-like wings of the neutron scattering density and the crystalline water hydrogen positions. An unrestrained simulation also was performed. For this simulation, the recall of crystallographic waters was much lower. For both restrained and unrestrained simulations, the strongest water density peaks were associated with crystallographic waters. The results demonstrate that it is now possible to recover crystallographic water structure using restrained MD simulations but that it is not yet reasonable to expect unrestrained MD simulations to do the same. Further development and generalization of MD water models for force-field development, macromolecular crystallography, and medicinal chemistry applications is now warranted. In particular, the combination of room-temperature crystallography, neutron diffraction, and crystalline MD simulations promises to substantially advance modeling of biomolecular solvation.

Binding Thermodynamics of Host–Guest Systems with SMIRNOFF99Frosst 1.0.5 from the Open Force Field Initiative

Abstract: Designing ligands that bind their target biomolecules with high affinity and specificity is a key step in small-molecule drug discovery, but accurately predicting protein–ligand binding free energi...

Binding Modes and Metabolism of Caffeine

Abstract: A correct estimate of ligand binding modes and a ratio of their occupancies is crucial for calculations of binding free energies. The newly developed method BLUES combines molecular dynamics with nonequilibrium candidate Monte Carlo. Nonequilibrium candidate Monte Carlo generates a plethora of possible binding modes and molecular dynamics enables the system to relax. We used BLUES to investigate binding modes of caffeine in the active site of its metabolizing enzyme Cytochrome P450 1A2 with the aim of elucidating metabolite-formation profiles at different concentrations. Because the activation energies of all sites of metabolism do not show a clear preference for one metabolite over the others, the orientations in the active site must play a key role. In simulations with caffeine located in a spacious pocket above the I-helix, it points N3 and N1 to the heme iron, whereas in simulations where caffeine is in close proximity to the heme N7 and C8 are preferably oriented toward the heme iron. We propose a mechanism where at low caffeine concentrations caffeine binds to the upper part of the active site, leading to formation of the main metabolite paraxanthine. On the other hand, at high concentrations two molecules are located in the active site, forcing one molecule into close proximity to the heme and yielding metabolites theophylline and trimethyluretic acid. Our results offer an explanation of previously published experimental results.

Enhancing Side Chain Rotamer Sampling Using Nonequilibrium Candidate Monte Carlo.

Abstract: Molecular simulations are a valuable tool for studying biomolecular motions and thermodynamics. However, such motions can be slow compared to simulation time scales, yet critical. Specifically, adequate sampling of side chain motions in protein binding pockets is crucial for obtaining accurate estimates of ligand binding free energies from molecular simulations. The time scale of side chain rotamer flips can range from a few ps to several hundred ns or longer, particularly in crowded environments like the interior of proteins. Here, we apply a mixed nonequilibrium candidate Monte Carlo (NCMC)/molecular dynamics (MD) method to enhance sampling of side chain rotamers. The NCMC portion of our method applies a switching protocol wherein the steric and electrostatic interactions between target side chain atoms and the surrounding environment are cycled off and then back on during the course of a move proposal. Between NCMC move proposals, simulation of the system continues via traditional molecular dynamics. Here, we first validate this approach on a simple, solvated valine-alanine dipeptide system and then apply it to a well-studied model ligand binding site in T4 lysozyme L99A. We compute the rate of rotamer transitions for a valine side chain using our approach and compare it to that of traditional molecular dynamics simulations. Here, we show that our NCMC/MD method substantially enhances side chain sampling, especially in systems where the torsional barrier to rotation is high (≥10 kcal/mol). These barriers can be intrinsic torsional barriers or steric barriers imposed by the environment. Overall, this may provide a promising strategy to selectively improve side chain sampling in molecular simulations.

Assessing the Conformational Equilibrium of Carboxylic Acid via Quantum Mechanical and Molecular Dynamics Studies on Acetic Acid.

Abstract: Accurate hydrogen placement in molecular modeling is crucial for studying the interactions and dynamics of biomolecular systems. The carboxyl functional group is a prototypical example of a functional group that requires protonation during structure preparation. To our knowledge, when in their neutral form, carboxylic acids are typically protonated in the syn conformation by default in classical molecular modeling packages, with no consideration of alternative conformations, though we are not aware of any careful examination of this topic. Here, we investigate the general belief that carboxylic acids should always be protonated in the syn conformation. We calculate and compare the relative energetic stabilities of syn and anti acetic acid using ab initio quantum mechanical calculations and atomistic molecular dynamics simulations. We focus on the carboxyl torsional potential and configurations of microhydrated acetic acid from molecular dynamics simulations, probing the effects of solvent, force field (GAFF vs GAFF2), and partial charge assignment of acetic acid. We show that while the syn conformation is the preferred state, the anti state may in some cases also be present under normal NPT conditions in solution.

Understanding the role of intermolecular interactions between lissoclimides and the eukaryotic ribosome

Abstract: Natural products that target the eukaryotic ribosome are promising therapeutics to treat a variety of cancers. It is therefore essential to determine their molecular mechanism of action to fully understand their mode of interaction with the target and to inform the development of new synthetic compounds with improved potency and reduced cytotoxicity. Toward this goal, we have previously established a short synthesis pathway that grants access to multiple congeners of the lissoclimide family. Here we present the X-ray co-crystal structure at 3.1 A resolution of C45, a potent congener with two A-ring chlorine-bearing stereogenic centers with 'unnatural' configurations, with the yeast 80S ribosome, intermolecular interaction energies of the C45/ribosome complex, and single-molecule FRET data quantifying the impact of C45 on both human and yeast ribosomes. Together, these data provide new insights into the role of unusual non-covalent halogen bonding interactions involved in the binding of this synthetic compound to the 80S ribosome.

Infinite Dilution Activity Coefficients as Constraints for Force Field Parametrization and Method Development.

Abstract: Molecular simulations begin with an underlying energy model or force field and from this can predict diverse physical properties. However, force fields were often developed with relatively limited data sets, yet accuracy for diverse properties across a broad chemical space is desirable; therefore, tests of such accuracy are particularly important. Here, to this end, we calculated 237 infinite dilution activity coefficients (IDACs), comparing with experimental values from NIST's ThermoML database. We found that calculated IDAC values correlate strongly with experiment (Pearson R of 0.92 ± 0.01) and allow us to identify specific functional groups that appear to present challenges to the force field employed. One potentially valuable aspect of IDACs, as compared to solvation free energies, which have been frequently employed as force field tests, is that the same molecules serve both as solutes and solvents in different cases, allowing us to ensure that force fields are not overly tuned to one particular environment or solvent.

Octanol-water partition coefficient measurements for the SAMPL6 Blind Prediction Challenge

Abstract: Partition coefficients describe the equilibrium partitioning of a neutral solute between two immiscible phases. Octanol-water partition coefficients (Kow), or their logarithms (log P), are frequently used as a measure of lipophilicity in drug discovery. The partition coefficient is a physicochemical property that captures the thermodynamics of relative solvation between aqueous and nonpolar phases, and therefore provides an excellent test for physics-based computational models that predict properties of pharmaceutical relevance such as protein-ligand binding affinities or hydration/solvation free energies. The SAMPL6 Part II Octanol-Water Partition Coefficient Prediction Challenge used a subset of kinase inhibitor fragment-like compounds from the SAMPL6 pKa Prediction Challenge in a blind experimental benchmark. Following experimental data collection, the partition coefficient dataset was kept blinded until all predictions were collected from participating computational chemistry groups. A total of 91 submissions were received from 27 participating research groups. This paper presents the octanol-water log P dataset for this SAMPL6 Part II Partition Coefficient Challenge, which consisted of 11 compounds (six 4-aminoquinazolines, two benzimidazole, one pyrazolo[3,4-d]pyrimidine, one pyridine, one 2-oxoquinoline substructure containing compounds) with log P values in the range of 1.95-4.09. We describe the potentiometric log P measurement protocol used to collect this dataset using a Sirius T3, discuss the limitations of this experimental approach, and share suggestions for future log P data collection efforts for the evaluation of computational methods.

The SAMPL6 SAMPLing challenge: Assessing the reliability and efficiency of binding free energy calculations

Abstract: Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange - while displaying very small variance - can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

Liquid-like and rigid-body motions in molecular-dynamics simulations of a crystalline protein

Abstract: To gain insight into crystalline protein dynamics, we performed molecular-dynamics (MD) simulations of a periodic 2 × 2 × 2 supercell of staphylococcal nuclease. We used the resulting MD trajectories to simulate X-ray diffraction and to study collective motions. The agreement of simulated X-ray diffraction with the data is comparable to previous MD simulation studies. We studied collective motions by analyzing statistically the covariance of alpha-carbon position displacements. The covariance decreases exponentially with the distance between atoms, which is consistent with a liquidlike motions (LLM) model, in which the protein behaves like a soft material. To gain finer insight into the collective motions, we examined the covariance behavior within a protein molecule (intraprotein) and between different protein molecules (interprotein). The interprotein atom pairs, which dominate the overall statistics, exhibit LLM behavior; however, the intraprotein pairs exhibit behavior that is consistent with a superposition of LLM and rigid-body motions (RBM). Our results indicate that LLM behavior of global dynamics is present in MD simulations of a protein crystal. They also show that RBM behavior is detectable in the simulations but that it is subsumed by the LLM behavior. Finally, the results provide clues about how correlated motions of atom pairs both within and across proteins might manifest in diffraction data. Overall, our findings increase our understanding of the connection between molecular motions and diffraction data and therefore advance efforts to extract information about functionally important motions from crystallography experiments.

Force Field Partial Charges with Restrained Electrostatic Potential 2 (RESP2)

Abstract: Many molecular simulation force fields represent the charge distributions of molecules with atom-centered partial charges, so simulations with these force fields require that partial charges be assigned to the molecules of interest. The restrained electrostatic potential (RESP) approach is a highly regarded and widely used method of assigning partial charges to varied organic compounds. RESP uses gas-phase HF/6-31G* as the underlying quantum chemical method, intending the resulting overpolarization of molecules to approximate the self-polarization that occurs in the condensed phase setting. However, it is far from clear that this fortuitous overpolarization is optimal or consistent across all compounds. In order to reach a higher level of accuracy, we propose a next generation of this approach, termed RESP2. In RESP2, the charges are derived from higher-level quantum chemical calculations carried out for both gas and aqueous phase, the latter using a continuum solvent model. The polarity of the final charges is tuned by a mixing parameter, δ, which scales the relative contributions of the gas- and aqueous-phase charges. We find that simply substituting RESP2 charges for RESP charges in the context of regular LJ parameters does not lead to clear improvement in liquid-state densities and heats of vaporization but does improve the accuracy of observables expected to depend most strongly on the accuracy of the charge model, i.e., dielectric constants and molecular dipole moments. However, when Lennard-Jones (LJ) parameters are optimized in the context of RESP charges, based on liquid properties, significant improvement in accuracy can be achieved, even with a sharply reduced set of LJ types. We argue that RESP2 with δ≈0.6 (60% aqueous and 40% gas-phase charges) is an accurate and robust method of generating atom-centered partial charges. The present study also highlights the value of optimizing LJ parameters along with the electrostatic model and suggests that a small set of LJ types can be a good starting point for a systematic re-optimization of this important nonbonded term.

Structure of a Mycobacterium tuberculosis Heme-Degrading Protein, MhuD, Variant in Complex with Its Product.

Abstract: Mycobacterium tuberculosis (Mtb), the causative agent of tuberculosis, requires iron for survival. In Mtb, MhuD is the cytosolic protein that degrades imported heme. MhuD is distinct, in both sequence and structure, from canonical heme oxygenases (HOs) but homologous with IsdG-type proteins. Canonical HO is found mainly in eukaryotes, while IsdG-type proteins are predominantly found in prokaryotes, including pathogens. While there are several published structures of MhuD and other IsdG-type proteins in complex with the heme substrate, no structures of IsdG-type proteins in complex with a product have been reported, unlike the case for HOs. We recently showed that the Mtb variant MhuD-R26S produces biliverdin IXα (αBV) rather than the wild-type mycobilin isomers. Given that mycobilin and other IsdG-type protein products like staphylobilin are difficult to isolate in quantities sufficient for structure determination, here we use the MhuD-R26S variant and its product αBV as a proxy to study the IsdG-type protein-product complex. First, we show that αBV has a nanomolar affinity for MhuD and the R26S variant. Second, we determined the MhuD-R26S-αBV complex structure to 2.5 A, which reveals two notable features: (1) two αBV molecules bound per active site and (2) a novel α-helix (α3) that was not observed in previous MhuD-heme structures. Finally, through molecular dynamics simulations, we show that α3 is stable with the proximal αBV alone. MhuD's high affinity for the product and the observed structural and electrostatic changes that accompany substrate turnover suggest that there may be an unidentified class of proteins that are responsible for the extraction of products from MhuD and other IsdG-type proteins.

openforcefield/smirnoff99Frosst: Version 1.0.9

Abstract: Addresses #89: Fixes torsion t56 , where SMIRKS [!1:1]-[#7X4,#7X3:2]-[#6X4;r3:3]-[*:4] should be [!#1:1]-[#7X4,#7X3:2]-[#6X4;r3:3]-[*:4] . The first means "not an isotope with mass 1" ( !1 ), but we intend for it to apply to all hydrogens, so it has been changed to !#1 .

Radiofrequency ablation for lung carcinomas: A retrospective review of a high-risk patient population at a community hospital

Abstract: Purpose: The purpose of this study is to retrospectively evaluate the technical efficacy, safety, and treatment outcomes of percutaneous radiofrequency ablation (RFA) of lung tumors in patients not amenable to surgery at an urban community hospital. Materials and Methods: Informed consent and IRB approval was obtained. Eligible tumors were defined as those in patients deemed poor surgical candidates by multidisciplinary consensus or those refusing surgery. Response to treatment was assessed by computed tomography (CT) performed immediately postprocedure and regular intervals up to 36 months later. Complete response was measured as a 30% decrease in mean tumor diameter without evidence of contrast enhancement or tumor growth within the ablation zone as defined by the response evaluation in solid tumors. Patient demographics, technical success, postprocedure complications, and survival were assessed and compared with data available in literature. Results: Twenty-four patients with a total of 29 tumors underwent percutaneous CT guided RFA for biopsy-proven lung malignancies between 2010 and 2016. Complete response was achieved in 82% (14/17) of treated tumors in patients who complied with postprocedure imaging recommendations. Immediate postprocedure complications occurred following 27.6% (8/29) ablations with pneumothorax being the most common, 17.2% (6/29). Mean survival is 28.5 months (95% confidence interval: 19.7–37.3). Progressive disease was seen in 18% (3/17) patients. No immediate treatment mortality was found. No significant difference was found in survival in patients with multiple comorbidities as measured by the Charlson Comorbidity Index. Conclusions: RFA of lung tumors is a well-tolerated procedure with low incidence of minor complications, a good tumor response and survival benefit in selected patients in the community setting. This is a positive endorsement of the potential success of tumor RFA programs outside of the academic setting. In addition, patients with multiple comorbidities should still be considered candidates for RFA as no difference was seen in survival in patients with multiple medical comorbidities.

Liquid-like and rigid-body motions in molecular-dynamics simulations of a crystalline protein

Abstract: To gain insight into crystalline protein dynamics, we performed molecular-dynamics (MD) simulations of a periodic 2×2×2 supercell of staphylococcal nuclease. We used the resulting MD trajectories to simulate X-ray diffraction and to study collective motions. The agreement of simulated X-ray diffraction with the data is comparable to previous MD simulation studies. We studied collective motions by analyzing statistically the covariance of alpha-carbon position displacements. The covariance decreases exponentially with the distance between atoms, which is consistent with a liquid-like motions (LLM) model, in which the protein behaves like a soft material. To gain finer insight into the collective motions, we examined the covariance behavior within a protein molecule (intra-protein) and between different protein molecules (inter-protein). The inter-protein atom pairs, which dominate the overall statistics, exhibit LLM behavior; however, the intra-protein pairs exhibit behavior that is consistent with a superposition of LLM and rigid-body motions (RBM). Our results indicate that LLM behavior of global dynamics is present in MD simulations of a protein crystal. They also show that RBM behavior is detectable in the simulations but that it is subsumed by the LLM behavior. Finally the results provide clues about how correlated motions of atom pairs both within and across proteins might manifest in diffraction data. Overall our findings increase our understanding of the connection between molecular motions and diffraction data, and therefore advance efforts to extract information about functionally important motions from crystallography experiments.

The structure of Mycobacterium tuberculosis heme-degrading protein, MhuD, in complex with product

Abstract: Mycobacterium tuberculosis (Mtb), the causative agent of tuberculosis, requires iron for survival. In Mtb, MhuD is the cytosolic protein that degrades imported heme. MhuD is distinct, both in sequence and structure, from canonical heme oxygenases (HOs) but homologous with IsdG-type proteins. Canonical HO is found mainly in eukaryotes, while IsdG-type proteins are predominantly found in prokaryotes including pathogens. While there are several published structures of MhuD and other IsdG-type proteins in complex with heme substrate, no structures have been reported of IsdG-type proteins in complex with product, unlike HOs. We recently showed that the Mtb variant MhuD-R26S produces biliverdin IXα (αBV) rather than the wild-type (WT) mycobilin isomers as product. Given that mycobilin and other IsdG-type protein products like staphylobilin are difficult to isolate in quantities sufficient for structure determination, here we use the MhuD-R26S variant and its product αBV as a proxy to study the IsdG-type protein/product complex. First we show that αBV has nanomolar affinity for MhuD and the R26S variant. Second we determined the MhuD-R26S-αBV complex structure to 2.5 A, which reveals two notable features (1) two αBV molecules bound per active site and (2) a new α-helix (α3) as compared with the MhuD-heme structure. Finally, by molecular dynamics simulations we show that α3 is stable with the proximal αBV alone. MhuD’s high affinity for its product and structural and electrostatic changes that accompany substrate turnover suggest that there is an unidentified protein that is responsible for product extraction from MhuD and other IsdG-type proteins.