Sire: An interoperability engine for prototyping algorithms and exchanging information between molecular simulation programs.

Sire is a Python/C++ library that is used both to prototype new algorithms and as an interoperability engine for exchanging information between molecular simulation programs. It provides a collection of file parsers and information converters that together make it easier to combine and leverage the functionality of many other programs and libraries. This empowers researchers to use sire to write a single script that can, for example, load a molecule from a PDBx/mmCIF file via Gemmi, perform SMARTS searches via RDKit, parameterize molecules using BioSimSpace, run GPU-accelerated molecular dynamics via OpenMM, and then display the resulting dynamics trajectory in a NGLView Jupyter notebook 3D molecular viewer. This functionality is built on by BioSimSpace, which uses sire's molecular information engine to interconvert with programs such as GROMACS, NAMD, Amber, and AmberTools for automated molecular parameterization and the running of molecular dynamics, metadynamics, and alchemical free energy workflows. Sire comes complete with a powerful molecular information search engine, plus trajectory loading and editing, analysis, and energy evaluation engines. This, when combined with an in-built computer algebra system, gives substantial flexibility to researchers to load, search for, edit, and combine molecular information from multiple sources and use that to drive novel algorithms by combining functionality from other programs. Sire is open source (GPL3) and is available via conda and at a free Jupyter notebook server at https://try.openbiosim.org. Sire is supported by the not-for-profit OpenBioSim community interest company.

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

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.

Blinded predictions of standard binding free energies: lessons learned from the SAMPL6 challenge.

In the context of the SAMPL6 challenges, series of blinded predictions of standard binding free energies were made with the SOMD software for a dataset of 27 host-guest systems featuring two octa-acids hosts (OA and TEMOA) and a cucurbituril ring (CB8) host. Three different models were used, ModelA computes the free energy of binding based on a double annihilation technique; ModelB additionally takes into account long-range dispersion and standard state corrections; ModelC additionally introduces an empirical correction term derived from a regression analysis of SAMPL5 predictions previously made with SOMD. The performance of each model was evaluated with two different setups; buffer explicitly matches the ionic strength from the binding assays, whereas no-buffer merely neutralizes the host-guest net charge with counter-ions. ModelC/no-buffer shows the lowest mean-unsigned error for the overall dataset (MUE 1.29 < 1.39 < 1.50 kcal mol-1, 95% CI), while explicit modelling of the buffer improves significantly results for the CB8 host only. Correlation with experimental data ranges from excellent for the host TEMOA (R2 0.91 < 0.94 < 0.96), to poor for CB8 (R2 0.04 < 0.12 < 0.23). Further investigations indicate a pronounced dependence of the binding free energies on the modelled ionic strength, and variable reproducibility of the binding free energies between different simulation packages.

Blinded predictions of host-guest standard free energies of binding in the SAMPL5 challenge.

In the context of the SAMPL5 blinded challenge standard free energies of binding were predicted for a dataset of 22 small guest molecules and three different host molecules octa-acids (OAH and OAMe) and a cucurbituril (CBC). Three sets of predictions were submitted, each based on different variations of classical molecular dynamics alchemical free energy calculation protocols based on the double annihilation method. The first model (model A) yields a free energy of binding based on computed free energy changes in solvated and host-guest complex phases; the second (model B) adds long range dispersion corrections to the previous result; the third (model C) uses an additional standard state correction term to account for the use of distance restraints during the molecular dynamics simulations. Model C performs the best in terms of mean unsigned error for all guests (MUE [Formula: see text]-95 % confidence interval) for the whole data set and in particular for the octa-acid systems (MUE [Formula: see text]). The overall correlation with experimental data for all models is encouraging ([Formula: see text]). The correlation between experimental and computational free energy of binding ranks as one of the highest with respect to other entries in the challenge. Nonetheless the large MUE for the best performing model highlights systematic errors, and submissions from other groups fared better with respect to this metric.

Blinded predictions of distribution coefficients in the SAMPL5 challenge.

In the context of the SAMPL5 challenge water-cyclohexane distribution coefficients for 53 drug-like molecules were predicted. Four different models based on molecular dynamics free energy calculations were tested. All models initially assumed only one chemical state present in aqueous or organic phases. Model A is based on results from an alchemical annihilation scheme; model B adds a long range correction for the Lennard Jones potentials to model A; model C adds charging free energy corrections; model D applies the charging correction from model C to ionizable species only. Model A and B perform better in terms of mean-unsigned error ([Formula: see text] D units - 95 % confidence interval) and determination coefficient [Formula: see text], while charging corrections lead to poorer results with model D ([Formula: see text] and [Formula: see text]). Because overall errors were large, a retrospective analysis that allowed co-existence of ionisable and neutral species of a molecule in aqueous phase was investigated. This considerably reduced systematic errors ([Formula: see text] and [Formula: see text]). Overall accurate [Formula: see text] predictions for drug-like molecules that may adopt multiple tautomers and charge states proved difficult, indicating a need for methodological advances to enable satisfactory treatment by explicit-solvent molecular simulations.

Effect of set up protocols on the accuracy of alchemical free energy calculation over a set of ACK1 inhibitors.

Hit-to-lead virtual screening frequently relies on a cascade of computational methods that starts with rapid calculations applied to a large number of compounds and ends with more expensive computations restricted to a subset of compounds that passed initial filters. This work focuses on set up protocols for alchemical free energy (AFE) scoring in the context of a Docking-MM/PBSA-AFE cascade. A dataset of 15 congeneric inhibitors of the ACK1 protein was used to evaluate the performance of AFE set up protocols that varied in the steps taken to prepare input files (using previously docked and best scored poses, manual selection of poses, manual placement of binding site water molecules). The main finding is that use of knowledge derived from X-ray structures to model binding modes, together with the manual placement of a bridging water molecule, improves the R2 from 0.45 ± 0.06 to 0.76 ± 0.02 and decreases the mean unsigned error from 2.11 ± 0.08 to 1.24 ± 0.04 kcal mol-1. By contrast a brute force automated protocol that increased the sampling time ten-fold lead to little improvements in accuracy. Besides, it is shown that for the present dataset hysteresis can be used to flag poses that need further attention even without prior knowledge of experimental binding affinities.

Reproducibility of Free Energy Calculations across Different Molecular Simulation Software Packages.

Alchemical free energy calculations are an increasingly important modern simulation technique to calculate free energy changes on binding or solvation. Contemporary molecular simulation software such as AMBER, CHARMM, GROMACS, and SOMD include support for the method. Implementation details vary among those codes, but users expect reliability and reproducibility, i.e., for a given molecular model and set of force field parameters, comparable free energy differences should be obtained within statistical bounds regardless of the code used. Relative alchemical free energy (RAFE) simulation is increasingly used to support molecule discovery projects, yet the reproducibility of the methodology has been less well tested than its absolute counterpart. Here we present RAFE calculations of hydration free energies for a set of small organic molecules and demonstrate that free energies can be reproduced to within about 0.2 kcal/mol with the aforementioned codes. Absolute alchemical free energy simulations have been carried out as a reference. Achieving this level of reproducibility requires considerable attention to detail and package-specific simulation protocols, and no universally applicable protocol emerges. The benchmarks and protocols reported here should be useful for the community to validate new and future versions of software for free energy calculations.

Metabolism and hydrophilicity of the polarised 'Janus face' all-cis tetrafluorocyclohexyl ring, a candidate motif for drug discovery.

The metabolism and polarity of the all-cis tetra-fluorocyclohexane motif is explored in the context of its potential as a motif for inclusion in drug discovery programmes. Biotransformations of phenyl all-cis tetra-, tri- and di- fluoro cyclohexanes with the human metabolism model organism Cunninghamella elegans illustrates various hydroxylated products, but limited to benzylic hydroxylation for the phenyl all-cis tetrafluorocyclohexyl ring system. Evaluation of the lipophilicities (log P) indicates a significant and progressive increase in polarity with increasing fluorination on the cyclohexane ring system. Molecular dynamics simulations indicate that water associates much more closely with the hydrogen face of these Janus face cyclohexyl rings than the fluorine face owing to enhanced hydrogen bonding interactions with the polarised hydrogens and water.

Arsenic toxicity and HSP70 expression in Xenopus laevis embryos.

The evaluation of the effect of trace metals on health can be difficult, because of their presence in the environment in various chemical forms. Exposure to arsenic compounds is an example of this complexity, as it can be present in the environment in inorganic and organic forms. The effects of arsenic in vertebrates are complicated by several variables, such as speciation of the element, the exposure route, and the susceptibility of the particular animal species. The embryotoxicity and teratogenicity of three arsenic species - sodium arsenite (NaAsO(2)), disodium hydrogen arsenate (Na(2)HAsO(4)) and dimethylarsinic acid [(CH3)2AsOOH] - were evaluated by the modified frog embryo teratogenic assay on Xenopus (FETAX). We also show how the classical FETAX endpoints, based on morphological and morphometrical analysis, can conveniently be integrated with the study of molecular markers. Possible changes in the expression of the mRNA for the heat-shock protein HSP70, following exposure to NaAsO(2), were examined by using the reverse transcriptase polymerase chain reaction. HSP70 mRNA is strongly induced by arsenic.

Developmental toxicity, uptake and distribution of sodium chromate assayed by frog embryo teratogenesis assay-Xenopus(FETAX).

The embryotoxicity and teratogenicity of Cr(VI) on the survival and morphology of the anuran Xenopus laevis have been assessed by frog embryo teratogenesis assay-Xenopus (FETAX). The lethal median (LC(50)) and teratogenic median (TC(50)) concentration values of Cr(VI) were 890 microM and 260 microM, respectively. The calculated teratogenic index (TI) value was 3.42, suggesting that hexavalent chromium has a teratogenic potential. Malformations of embryos included lifting of the body, coiling of the tail and body oedema. Furthermore, the chromium salt caused significant growth retardation at 25 microM exposure concentrations. The use of radiolabelled (51)Cr(VI) allowed the determination of the time course uptake of Cr in Xenopus exposed to concentrations ranging from 0.025 to 500 microM. The evaluation of its distribution into the body (head-abdomen-tail) was evaluated at different exposure times. Chromium is taken up at 24 h by Xenopus embryos for all concentrations tested. At 48 h post fertilization (stage of larva) the amount of Cr accumulated by the two-day-old larva ranged from 0.42 to 580 pg mg(-1) wet weight at 0.025 and 500 microM respectively. These amounts were lower than those at 24 h (2.77 to 11016 pg mg(-1) wet weight embryo) reaching values of the same order of magnitude at 120 h (five-days-old larva). Since at 48 h Xenopus development leads to a swimming embryo, the observed uptake at 24 h could be the result of the binding of Cr to jelly coat compounds surrounding the embryo body as confirmed by gel filtration experiments on (51)Cr-jelly coat. The interaction of Cr with jelly coat is in agreement with the role of jelly coat in protecting the embryo against pathogen and chemical toxins to ensure fertilization. This work further supports the hypothesis that Cr contamination of surface waters could contribute to explain the reported worldwide depletion of frog population.