Open-Source Machine Learning in Computational Chemistry.

The field of computational chemistry has seen a significant increase in the integration of machine learning concepts and algorithms. In this Perspective, we surveyed 179 open-source software projects, with corresponding peer-reviewed papers published within the last 5 years, to better understand the topics within the field being investigated by machine learning approaches. For each project, we provide a short description, the link to the code, the accompanying license type, and whether the training data and resulting models are made publicly available. Based on those deposited in GitHub repositories, the most popular employed Python libraries are identified. We hope that this survey will serve as a resource to learn about machine learning or specific architectures thereof by identifying accessible codes with accompanying papers on a topic basis. To this end, we also include computational chemistry open-source software for generating training data and fundamental Python libraries for machine learning. Based on our observations and considering the three pillars of collaborative machine learning work, open data, open source (code), and open models, we provide some suggestions to the community.

Determining Lennard-Jones Parameters Using Multiscale Target Data through Presampling-Enhanced, Surrogate-Assisted Global Optimization.

Force field-based models are a Newtonian mechanics approximation of reality and are inherently noisy. Coupling models from different molecular scale domains (including single, gas-phase molecules up to multimolecule, condensed phase ensembles) is difficult, which is also the case for finding solutions that transfer well between the scales. In this contribution, we introduce a surrogate-assisted algorithm to optimize Lennard-Jones parameters for target data from different scale domains to overcome the difficulties named above. Specifically, our approach combines a surrogate-assisted global evolutionary optimization method with a presampling phase that takes advantage of one scale domain being less computationally expensive to evaluate. The algorithm's components were evaluated individually, elucidating their individual merits. Our findings show that the process of parametrizing force fields can significantly benefit from both the presampling method, which alleviates the need to have a good initial guess for the parameters, and the surrogate model, which improves efficiency.

Mebendazole's Conformational Space and Its Predicted Binding to Human Heat-Shock Protein 90.

Recent experimental evidence suggests that mebendazole, a popular antiparasitic drug, binds to heat shock protein 90 (Hsp90) and inhibits acute myeloid leukemia cell growth. In this study we use quantum mechanics (QM), molecular similarity, and molecular dynamics (MD) calculations to predict possible binding poses of mebendazole to the adenosine triphosphate (ATP) binding site of Hsp90. Extensive conformational searches and minimization of the five mebendazole tautomers using the MP2/aug-cc-pVTZ theory level resulted in 152 minima. Mebendazole-Hsp90 complex models were subsequently created using the QM optimized conformations and protein coordinates obtained from experimental crystal structures that were chosen through similarity calculations. Nine different poses were identified from a total of 600 ns of explicit solvent, all-atom MD simulations using two different force fields. All simulations support the hypothesis that mebendazole is able to bind to the ATP binding site of Hsp90.

Mebendazole Mediates Proteasomal Degradation of GLI Transcription Factors in Acute Myeloid Leukemia.

The prognosis of elderly AML patients is still poor due to chemotherapy resistance. The Hedgehog (HH) pathway is important for leukemic transformation because of aberrant activation of GLI transcription factors. MBZ is a well-tolerated anthelmintic that exhibits strong antitumor effects. Herein, we show that MBZ induced strong, dose-dependent anti-leukemic effects on AML cells, including the sensitization of AML cells to chemotherapy with cytarabine. MBZ strongly reduced intracellular protein levels of GLI1/GLI2 transcription factors. Consequently, MBZ reduced the GLI promoter activity as observed in luciferase-based reporter assays in AML cell lines. Further analysis revealed that MBZ mediates its anti-leukemic effects by promoting the proteasomal degradation of GLI transcription factors via inhibition of HSP70/90 chaperone activity. Extensive molecular dynamics simulations were performed on the MBZ-HSP90 complex, showing a stable binding interaction at the ATP binding site. Importantly, two patients with refractory AML were treated with MBZ in an off-label setting and MBZ effectively reduced the GLI signaling activity in a modified plasma inhibitory assay, resulting in a decrease in peripheral blood blast counts in one patient. Our data prove that MBZ is an effective GLI inhibitor that should be evaluated in combination to conventional chemotherapy in the clinical setting.

Structure-activity relationships of thiostrepton derivatives: implications for rational drug design.

The bacterial ribosome is a major target of naturally occurring thiopeptides antibiotics. Studying thiopeptide (e.g. thiostrepton) binding to the GAR's 23S·L11 ribosomal subunit using docking methods is challenging. Regarding the target, the binding site is composed of a flexible protein-RNA nonbonded interface whose available crystal structure is of medium resolution. Regarding the ligands, the thiopeptides are chemically complex, flexible, and contain macrocycles. In this study we developed a combined MD-docking-MD workflow that allows us to study thiopeptide-23S·L11 binding. It is shown that docking thiostrepton-like ligands to an MD-refined receptor structure instead of the medium resolution crystal leads to better convergence to the native-like docking pose and a better reproduction of experimental binding affinities. By applying an energy decomposition analysis, we identify key structural binding elements within GAR's rRNA-protein binding site and within the ligand structures.

Wolf2Pack--portal based atomistic force-field development.

In this contribution we introduce the technical concept and implementation details concerning the front end of our force-field optimization workflow package for intramolecular degrees of freedom, called Wolf2Pack. The package's design follows our belief that parameter optimization should be a user-driven, but program guided, workflow with specific modular tasks that reduce human errors and save time. Through this design, parameter optimization becomes more reliable and reproducible. Wolf2Pack can integrate common force fields from different research areas, allowing the user to optimize balanced parameters; alternatively users can develop highly specialized force fields that suite their chemical systems. Included in the package's front end is a force-field and molecular database whose contents facilitate parameter optimization. Wolf2Pack can be accessed at www.wolf2pack.com.

Influence of thiostrepton binding on the ribosomal GTPase associated region characterized by molecular dynamics simulation.

The thiostrepton antibiotic inhibits bacterial protein synthesis by binding to a cleft formed by the ribosomal protein L11 and 23S's rRNA helices 43-44 on the 70S ribosome. It was proposed from crystal structures that the ligand restricts L11's N-terminal movement and thus prevents proper translation factor binding. An exact understanding of thiostrepton's impact on the binding site's dynamics at atomistic resolution is still missing. Here we report an all-atom molecular dynamics simulations of the binary L11·rRNA and the ternary L11·rRNA·thiostrepton complex (rRNA = helices 43-44). We demonstrate that thiostrepton directly impacts the binding site's atomic and biomacromolecular dynamics.

Quantum mechanical study of sulfuric acid hydration: atmospheric implications.

The role of the binary nucleation of sulfuric acid in aerosol formation and its implications for global warming is one of the fundamental unsettled questions in atmospheric chemistry. We have investigated the thermodynamics of sulfuric acid hydration using ab initio quantum mechanical methods. For H(2)SO(4)(H(2)O)(n) where n = 1-6, we used a scheme combining molecular dynamics configurational sampling with high-level ab initio calculations to locate the global and many low lying local minima for each cluster size. For each isomer, we extrapolated the Møller-Plesset perturbation theory (MP2) energies to their complete basis set (CBS) limit and added finite temperature corrections within the rigid-rotor-harmonic-oscillator (RRHO) model using scaled harmonic vibrational frequencies. We found that ionic pair (HSO(4)(-)·H(3)O(+))(H(2)O)(n-1) clusters are competitive with the neutral (H(2)SO(4))(H(2)O)(n) clusters for n ≥ 3 and are more stable than neutral clusters for n ≥ 4 depending on the temperature. The Boltzmann averaged Gibbs free energies for the formation of H(2)SO(4)(H(2)O)(n) clusters are favorable in colder regions of the troposphere (T = 216.65-273.15 K) for n = 1-6, but the formation of clusters with n ≥ 5 is not favorable at higher (T > 273.15 K) temperatures. Our results suggest the critical cluster of a binary H(2)SO(4)-H(2)O system must contain more than one H(2)SO(4) and are in concert with recent findings (1) that the role of binary nucleation is small at ambient conditions, but significant at colder regions of the troposphere. Overall, the results support the idea that binary nucleation of sulfuric acid and water cannot account for nucleation of sulfuric acid in the lower troposphere.

Ramachandran-type plots for glycosidic linkages: Examples from molecular dynamic simulations using the Glycam06 force field.

The goals of this article are to (1) provide further validation of the Glycam06 force field, specifically for its use in implicit solvent molecular dynamic (MD) simulations, and (2) to present the extension of G.N. Ramachandran's idea of plotting amino acid phi and psi angles to the glycosidic phi, psi, and omega angles formed between carbohydrates. As in traditional Ramachandran plots, these carbohydrate Ramachandran-type (carb-Rama) plots reveal the coupling between the glycosidic angles by displaying the allowed and disallowed conformational space. Considering two-bond glycosidic linkages, there are 18 possible conformational regions that can be defined by (alpha, phi, psi) and (beta, phi, psi), whereas for three-bond linkages, there are 54 possible regions that can be defined by (alpha, phi, psi, omega) and (beta, phi, psi, omega). Illustrating these ideas are molecular dynamic simulations on an implicitly hydrated oligosaccharide (700 ns) and its eight constituent disaccharides (50 ns/disaccharide). For each linkage, we compare and contrast the oligosaccharide and respective disaccharide carb-Rama plots, validate the simulations and the Glycam06 force field through comparison to experimental data, and discuss the general trends observed in the plots.

Interactions of Water and Alkanes: Modifying Additive Force Fields to Account for Polarization Effects.

Atomistic biomolecular simulations predominantly utilize additive force fields (FF), where the electrostatic potential is modeled by fixed point charges. Among other consequences, the lack of polarizability in these models undermines the balance of hydrophilic/hydrophobic nonbonded interactions. Simulations of water/alkane systems using the TIP3P water model and CHARMM36 parameters reveal a 1 kcal/mol overestimate of the experimental transfer free energy of water to hexadecane; more recent optimized water models (SPC/E, TIP4P/2005, TIP4P-Ew, TIP3P-FB, TIP4P-FB, OPC, TIP4P-D) overestimate this transfer free energy by approximately 2 kcal/mol. In contrast, the polarizable SWM4-NDP and SWM6 water models reproduce experimental values to within statistical error. As an alternative to explicitly modeling polarizability, this paper develops an efficient automated workflow to optimize pair-specific Lennard-Jones parameters within an additive FF. Water/hexadecane is used as a prototype and the free energy of water transfer to hexadecane as a target. The optimized model yields quantitative agreement with the experimental transfer free energy and improves the water/hexadecane interfacial tension by 6%. Simulations of five different lipid bilayers show a strong increase of water permeabilities compared to the unmodified CHARMM36 lipid FF which consistently improves match with experiment: the order-of-magnitude underestimate for monounsaturated bilayers is rectified and the factor of 2.8-4 underestimate for saturated bilayers is turned into a factor of 1.5-3 overestimate. While agreement with experiment is decreased for the diffusion constant of water in hexadecane, alkane transfer free energies, and the bilayers' area per lipid, the method provides a permeant-specific route to achieve a wide range of heterogeneous observables via rapidly optimized pairwise parameters.

GLYCAM06: a generalizable biomolecular force field. Carbohydrates.

A new derivation of the GLYCAM06 force field, which removes its previous specificity for carbohydrates, and its dependency on the AMBER force field and parameters, is presented. All pertinent force field terms have been explicitly specified and so no default or generic parameters are employed. The new GLYCAM is no longer limited to any particular class of biomolecules, but is extendible to all molecular classes in the spirit of a small-molecule force field. The torsion terms in the present work were all derived from quantum mechanical data from a collection of minimal molecular fragments and related small molecules. For carbohydrates, there is now a single parameter set applicable to both alpha- and beta-anomers and to all monosaccharide ring sizes and conformations. We demonstrate that deriving dihedral parameters by fitting to QM data for internal rotational energy curves for representative small molecules generally leads to correct rotamer populations in molecular dynamics simulations, and that this approach removes the need for phase corrections in the dihedral terms. However, we note that there are cases where this approach is inadequate. Reported here are the basic components of the new force field as well as an illustration of its extension to carbohydrates. In addition to reproducing the gas-phase properties of an array of small test molecules, condensed-phase simulations employing GLYCAM06 are shown to reproduce rotamer populations for key small molecules and representative biopolymer building blocks in explicit water, as well as crystalline lattice properties, such as unit cell dimensions, and vibrational frequencies.

Efficient and accurate characterization of the Bergman cyclization for several enediynes including an expanded substructure of esperamicin A1.

Incorporation of enediynes into anticancer drugs remains an intriguing yet elusive strategy for the design of therapeutically active agents. Density functional theory was used to locate reactants, products, and transition states along the Bergman cyclization pathways connecting enediynes to reactive para-biradicals. Sum method correction to low-level calculations confirmed B3LYP/6-31G(d,p) as the method of choice in investigating enediynes. Herein described as MI:Sum, calculated reaction enthalpies differed from experiment by an average of 2.1 kcal x mol(-1) (mean unsigned error). A combination of strain energy released across the reaction coordinate and the critical intramolecular distance between reacting diynes explains reactivity differences. Where experimental and calculated barrier heights are in disagreement, higher level multireference treatment of the enediynes confirms lower level estimates. Previous work concerning the chemically reactive fragment of esperamcin, MTC, is expanded to our model system MTC2.

Binding free energy calculation for duocarmycin/DNA complex based on the QPLD-derived partial charge model.

The 3ns unrestrained MD simulations were carried out on the DNA/duocarmycin complex based on (1) the classic RESP charge model, and (2) the QM-polarized ligand docking (QPLD)-based charge model. The RMSDs of the trajectories and the DeltaG(bind) of the QPLD model perform much better than the RESP model, with the DeltaG(bind) estimation for QPLD model (-16.11 kcal/mol) versus DeltaG(bind) estimation for RESP model (-10.05 kcal/mol).

Hydration of OCS with one to four water molecules in atmospheric and laboratory conditions.

Carbonyl sulfide is the most abundant sulfur gas in the atmosphere. We have used MP2 and CCSD(T) theory to study the structures and thermochemistries of carbonyl sulfide interacting with one to four water molecules. We have completed an extensive search for clusters of OCS(H2O)n , where n = 1-4. We located three dimers, two trimers, five tetramers, and four pentamers with the MP2/aug-cc-pVDZ method. In each of the complexes with two or more waters, OCS preferentially interacts with low-energy water clusters. Our results match current theoretical and experimental literature, showing correlation with available geometries and frequencies for the OCS(H2O) species. The CCSD(T)/aug-cc-pVTZ thermochemical values combined with the average amount of OCS and the saturated concentration of H2O in the troposphere, lead to the prediction of 10(6) OCS(H2O) clusters x cm(-3) and 10(2) OCS(H2O)2 clusters x cm(-3) at 298 K. We predict the structures of OCS(H2O)n , n = 1-4 that should predominate in a low-temperature molecular beam and identify specific infrared vibrations that can be used to identify these different clusters.

Thermodynamics of the hydroxyl radical addition to isoprene.

Oxidation of isoprene by the hydroxyl radical leads to tropospheric ozone formation. Consequently, a more complete understanding of this reaction could lead to better models of regional air quality, a better understanding of aerosol formation, and a better understanding of reaction kinetics and dynamics. The most common first step in the oxidation of isoprene is the formation of an adduct, with the hydroxyl radical adding to one of four unsaturated carbon atoms in isoprene. In this paper, we discuss how the initial conformations of isoprene, s-trans and s-gauche, influences the pathways to adduct formation. We explore the formation of pre-reactive complexes at low and high temperatures, which are often invoked to explain the negative temperature dependence of this reaction's kinetics. We show that at higher temperatures the free energy surface indicates that a pre-reactive complex is unlikely, while at low temperatures the complex exists on two reaction pathways. The theoretical results show that at low temperatures all eight pathways possess negative reaction barriers, and reaction energies that range from -36.7 to -23.0 kcal x mol(-1). At temperatures in the lower atmosphere, all eight pathways possess positive reaction barriers that range from 3.8 to 6.0 kcal x mol(-1) and reaction energies that range from -28.8 to -14.4 kcal x mol(-1).

Small Alcohols Revisited: CCSD(T) Relative Potential Energies for the Minima, First- and Second-Order Saddle Points, and Torsion-Coupled Surfaces.

The elucidation of conformations and relative potential energies (rPEs) of small molecules has a long history across a diverse range of fields. Periodically, it is helpful to revisit what conformations have been investigated and to provide a consistent theoretical framework for which clear comparisons can be made. In this paper, we compute the minima, first- and second-order saddle points, and torsion-coupled surfaces for methanol, ethanol, propan-2-ol, and propanol using consistent high-level MP2 and CCSD(T) methods. While for certain molecules more rigorous methods were employed, the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pV5Z theory level was used throughout to provide relative energies of all minima and first-order saddle points. The rPE surfaces were uniformly computed at the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVTZ level. To the best of our knowledge, this represents the most extensive study for alcohols of this kind, revealing some new aspects. Especially for propanol, we report several new conformations that were previously not investigated. Moreover, two metrics are included in our analysis that quantify how the selected surfaces are similar to one another and hence improve our understanding of the relationship between these alcohols.

Structural analysis of human glycoprotein butyrylcholinesterase using atomistic molecular dynamics: The importance of glycosylation site ASN241.

Human butyrylcholinesterase (BChE) is a glycoprotein capable of bioscavenging toxic compounds such as organophosphorus (OP) nerve agents. For commercial production of BChE, it is practical to synthesize BChE in non-human expression systems, such as plants or animals. However, the glycosylation profile in these systems is significantly different from the human glycosylation profile, which could result in changes in BChE's structure and function. From our investigation, we found that the glycan attached to ASN241 is both structurally and functionally important due to its close proximity to the BChE tetramerization domain and the active site gorge. To investigate the effects of populating glycosylation site ASN241, monomeric human BChE glycoforms were simulated with and without site ASN241 glycosylated. Our simulations indicate that the structure and function of human BChE are significantly affected by the absence of glycan 241.

Visualizing potential energy curves and conformations on ultra high-resolution display walls.

In this contribution, we examine how visualization on an ultra high-resolution display wall can augment force-field research in the field of molecular modeling. Accurate force fields are essential for producing reliable simulations, and subsequently important for several fields of applications (e.g. rational drug design and biomolecular modeling). We discuss how using HORNET, a recently constructed specific ultra high-resolution tiled display wall, enhances the visual analytics that are necessary for conformational-based interpretation of the raw data from molecular calculations. Simultaneously viewing multiple potential energy graphs and conformation overlays leads to an enhanced way of evaluating force fields and in their optimization. Consequently, we have integrated visual analytics into our existing Wolf2Pack workflow. We applied this workflow component to analyze how major AMBER force fields (Parm14SB, Gaff, Lipid14, Glycam06j) perform at reproducing the quantum mechanics relative energies and geometries of saturated hydrocarbons. Included in this comparison are the 1996 OPLS force field and our newly developed ExTrM force field. While we focus on atomistic force fields the ideas presented herein are generalizable to other research areas, particularly those that involve numerous representations of large data amounts and whose simultaneous visualization enhances the analysis.