Machine Learning of Reactive Potentials.

In the past two decades, machine learning potentials (MLPs) have driven significant developments in chemical, biological, and material sciences. The construction and training of MLPs enable fast and accurate simulations and analysis of thermodynamic and kinetic properties. This review focuses on the application of MLPs to reaction systems with consideration of bond breaking and formation. We review the development of MLP models, primarily with neural network and kernel-based algorithms, and recent applications of reactive MLPs (RMLPs) to systems at different scales. We show how RMLPs are constructed, how they speed up the calculation of reactive dynamics, and how they facilitate the study of reaction trajectories, reaction rates, free energy calculations, and many other calculations. Different data sampling strategies applied in building RMLPs are also discussed with a focus on how to collect structures for rare events and how to further improve their performance with active learning.

Assessment of Embedding Schemes in a Hybrid Machine Learning/Classical Potentials (ML/MM) Approach.

Machine learning (ML) methods have reached high accuracy levels for the prediction of in vacuo molecular properties. However, the simulation of large systems solely through ML methods (such as those based on neural network potentials) is still a challenge. In this context, one of the most promising frameworks for integrating ML schemes in the simulation of complex molecular systems are the so-called ML/MM methods. These multiscale approaches combine ML methods with classical force fields (MM), in the same spirit as the successful hybrid quantum mechanics-molecular mechanics methods (QM/MM). The key issue for such ML/MM methods is an adequate description of the coupling between the region of the system described by ML and the region described at the MM level. In the context of QM/MM schemes, the main ingredient of the interaction is electrostatic, and the state of the art is the so-called electrostatic-embedding. In this study, we analyze the quality of simpler mechanical embedding-based approaches, specifically focusing on their application within a ML/MM framework utilizing atomic partial charges derived in vacuo. Taking as reference electrostatic embedding calculations performed at a QM(DFT)/MM level, we explore different atomic charges schemes, as well as a polarization correction computed using atomic polarizabilites. Our benchmark data set comprises a set of about 80k small organic structures from the ANI-1x and ANI-2x databases, solvated in water. The results suggest that the minimal basis iterative stockholder (MBIS) atomic charges yield the best agreement with the reference coupling energy. Remarkable enhancements are achieved by including a simple polarization correction.

Accelerated Organic Crystal Structure Prediction with Genetic Algorithms and Machine Learning.

We present a high-throughput, end-to-end pipeline for organic crystal structure prediction (CSP)─the problem of identifying the stable crystal structures that will form from a given molecule based only on its molecular composition. Our tool uses neural network potentials to allow for efficient screening and structural relaxation of generated crystal candidates. Our pipeline consists of two distinct stages: random search, whereby crystal candidates are randomly generated and screened, and optimization, where a genetic algorithm (GA) optimizes this screened population. We assess the performance of each stage of our pipeline on 21 molecules taken from the Cambridge Crystallographic Data Centre's CSP blind tests. We show that random search alone yields matches for ≈50% of targets. We then validate the potential of our full pipeline, making use of the GA to optimize the root-mean-square deviation between crystal candidates and the experimentally derived structure. With this approach, we are able to find matches for ≈80% of candidates with 10-100 times smaller initial population sizes than when using random search. Lastly, we run our full pipeline with an ANI model that is trained on a small data set of molecules extracted from crystal structures in the Cambridge Structural Database, generating ≈60% of targets. By leveraging machine learning models trained to predict energies at the density functional theory level, our pipeline has the potential to approach the accuracy of ab initio methods and the efficiency of empirical force fields.

AmberTools.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

Axially Chiral Cannabinoids: Design, Synthesis, and Cannabinoid Receptor Affinity.

The resorcinol-terpene phytocannabinoid template is a privileged scaffold for the development of diverse therapeutics targeting the endocannabinoid system. Axially chiral cannabinols (axCBNs) are unnatural cannabinols (CBNs) that bear an additional C10 substituent, which twists the cannabinol biaryl framework out of planarity creating an axis of chirality. This unique structural modification is hypothesized to enhance both the physical and biological properties of cannabinoid ligands, thus ushering in the next generation of endocannabinoid system chemical probes and cannabinoid-inspired leads for drug development. In this full report, we describe the philosophy guiding the design of axCBNs as well as several synthetic strategies for their construction. We also introduce a second class of axially chiral cannabinoids inspired by cannabidiol (CBD), termed axially chiral cannabidiols (axCBDs). Finally, we provide an analysis of axially chiral cannabinoid (axCannabinoid) atropisomerism, which spans two classes (class 1 and 3 atropisomers), and provide first evidence that axCannabinoids retain─and in some cases, strengthen─affinity and functional activity at cannabinoid receptors. Together, these findings present a promising new direction for the design of novel cannabinoid ligands for drug discovery and exploration of the complex endocannabinoid system.

First-in-class multifunctional TYMS nonclassical antifolate inhibitor with potent in vivo activity that prolongs survival.

Although thymidylate synthase (TYMS) inhibitors have served as components of chemotherapy regimens, the currently available inhibitors induce TYMS overexpression or alter folate transport/metabolism feedback pathways that tumor cells exploit for drug resistance, limiting overall benefit. Here we report a small molecule TYMS inhibitor that i) exhibited enhanced antitumor activity as compared with current fluoropyrimidines and antifolates without inducing TYMS overexpression, ii) is structurally distinct from classical antifolates, iii) extended survival in both pancreatic xenograft tumor models and an hTS/Ink4a/Arf null genetically engineered mouse tumor model, and iv) is well tolerated with equal efficacy using either intraperitoneal or oral administration. Mechanistically, we verify the compound is a multifunctional nonclassical antifolate, and using a series of analogs, we identify structural features allowing direct TYMS inhibition while maintaining the ability to inhibit dihydrofolate reductase. Collectively, this work identifies nonclassical antifolate inhibitors that optimize inhibition of thymidylate biosynthesis with a favorable safety profile, highlighting the potential for enhanced cancer therapy.

Novel Lennard-Jones Parameters for Cysteine and Selenocysteine in the AMBER Force Field.

Cysteine is a common amino acid with a thiol group that plays a pivotal role in a variety of scenarios in redox biochemistry. In contrast, selenocysteine, the 21st amino acid, is only present in 25 human proteins. Classical force-field parameters for cysteine and selenocysteine are still scarce. In this context, we present a methodology to obtain Lennard-Jones parameters for cysteine and selenocysteine in different physiologically relevant oxidation and protonation states. The new force field parameters obtained in this work are available at https://github.com/MALBECC/AMBER-parameters-database. The parameters were adjusted to reproduce water radial distribution functions obtained by density functional theory ab initio molecular dynamics. We validated the results by evaluating the impact of the choice of parameters on the structure and dynamics in classical molecular dynamics simulations of representative proteins containing catalytic cysteine/selenocysteine residues. There are significant changes in protein structure and dynamics depending on the parameters choice, specifically affecting the residues close to the catalytic sites.

Elucidating the Protonation State of the γ-Secretase Catalytic Dyad.

γ-Secretase (GS) is an intramembrane aspartyl protease that participates in the sequential cleavage of C99 to generate different isoforms of the amyloid-ß (Aß) peptides that are associated with the development of Alzheimer's disease. Due to its importance in the proteolytic processing of C99 by GS, we performed pH replica exchange molecular dynamics (pH-REMD) simulations of GS in its apo and substrate-bound forms to sample the protonation states of the catalytic dyad. We found that the catalytic dyad is deprotonated at physiological pH in our apo form, but the presence of the substrate at the active site displaces its monoprotonated state toward physiological pH. Our results show that Asp257 acts as the general base and Asp385 as the general acid during the cleavage mechanism. We identified different amino acids such as Lys265, Arg269, and the PAL motif interacting with the catalytic dyad and promoting changes in its acid-base behavior. Finally, we also found a significant pKa shift of Glu280 related to the internalization of TM6-CT in the GS-apo form. Our study provides critical mechanistic insight into the GS mechanism and the basis for future research on the genesis of Aß peptides and the development of Alzheimer's disease.

Auto3D: Automatic Generation of the Low-Energy 3D Structures with ANI Neural Network Potentials.

Computational programs accelerate the chemical discovery processes but often need proper three-dimensional molecular information as part of the input. Getting optimal molecular structures is challenging because it requires enumerating and optimizing a huge space of stereoisomers and conformers. We developed the Python-based Auto3D package for generating the low-energy 3D structures using SMILES as the input. Auto3D is based on state-of-the-art algorithms and can automatize the isomer enumeration and duplicate filtering process, 3D building process, geometry optimization, and ranking process. Tested on 50 molecules with multiple unspecified stereocenters, Auto3D is guaranteed to find the stereoconfiguration that yields the lowest-energy conformer. With Auto3D, we provide an extension of the ANI model. The new model, dubbed ANI-2xt, is trained on a tautomer-rich data set. ANI-2xt is benchmarked with DFT methods on geometry optimization and electronic and Gibbs free energy calculations. Compared with ANI-2x, ANI-2xt provides a 42% error reduction for tautomeric reaction energy calculations when using the gold-standard coupled-cluster calculation as the reference. ANI-2xt can accurately predict the energies and is several orders of magnitude faster than DFT methods.

Adiabatic Excited-State Molecular Dynamics with an Explicit Solvent: NEXMD-SANDER Implementation.

We present a method to link the Nonadiabatic EXcited-state Molecular Dynamics (NEXMD) package to the SANDER package supplied by AMBERTOOLS to provide excited-state adiabatic quantum mechanics/molecular mechanics (QM/MM) simulations. NEXMD is a computational package particularly developed to perform simulations of the photoexcitation and subsequent nonadiabatic electronic and vibrational energy relaxation in large multichromophoric conjugated molecules involving several coupled electronic excited states. The NEXMD-SANDER exchange has been optimized in order to achieve excited-state adiabatic dynamics simulations of large conjugated materials in a QM/MM environment, such as an explicit solvent. Dynamics of a substituted polyphenylene vinylene oligomer (PPV3-NO2) in vacuum and different explicit solvents has been used as a test case by performing comparative analysis of changes in its optical spectrum, state-dependent conformational changes, and quantum bond orderings. The method has been tested and compared with respect to previous implicit solvent implementations. Also, the impact on the expansion of the QM region by including a variable number of solvent molecules has been analyzed. Altogether, these results encourage future implementations of NEXMD simulations using the same combination of methods.

Diastereoselective Indole-Dearomative Cope Rearrangements by Compounding Minor Driving Forces.

Reported herein is the discovery of a diastereoselective indole-dearomative Cope rearrangement. A suite of minor driving forces promote dearomatization: (i) steric congestion in the starting material, (ii) alkylidene malononitrile and stilbene conjugation events in the product, and (iii) an unexpected intramolecular π-π* stack on the product side of the equilibrium. The key substrates are rapidly assembled from simple starting materials, resulting in many successful examples. The products are structurally complex and bear vicinal stereocenters generated by the dearomative Cope rearrangement. They also contain a variety of functional groups for interconversion to complex architectures.

Anticancer Agents Derived from Cyclic Thiosulfonates: Structure-Reactivity and Structure-Activity Relationships.

Reported are structure-property-function relationships associated with a class of cyclic thiosulfonate molecules-disulfide-bond disrupting agents (DDAs)-with the ability to downregulate the Epidermal Growth Factor Receptor (HER) family in parallel and selectively induce apoptosis of EGFR+ or HER2+ breast cancer cells. Recent findings have revealed that the DDA mechanism of action involves covalent binding to the thiol(ate) from the active site cysteine residue of members of the protein disulfide isomerase (PDI) family. Reported is how structural modifications to the pharmacophore can alter the anticancer activity of cyclic thiosulfonates by tuning the dynamics of thiol-thiosulfonate exchange reactions, and the studies reveal a correlation between the biological potency and thiol-reactivity. Specificity of the cyclic thiosulfonate ring-opening reaction by a nucleophilic attack can be modulated by substituent addition to a parent scaffold. Lead compound optimization efforts are also reported, and have resulted in a considerable decrease of the IC50 /IC90 values toward HER-family overexpressing breast cancer cells.

Construction of vicinal 4°/3°-carbons via reductive Cope rearrangement.

Herein reported is a strategy for constructing vicinal 4°/3° carbons via reductive Cope rearrangement. Substrates have been designed which exhibit Cope rearrangement kinetic barriers of ∼23 kcal mol-1 with isoenergetic favorability (ΔG ∼ 0). These fluxional/shape-shifting molecules can be driven forward by chemoselective reduction to useful polyfunctionalized building blocks.

Development and Evaluation of Geometry Optimization Algorithms in Conjunction with ANI Potentials.

An efficient yet accurate method for producing a large amount of energy data for molecular mechanical force field (MMFF) parameterization is on demand, especially for torsional angle parameters which are typically derived to reproduce ab initio rotational profiles or torsional potential energy surfaces (PESs). Recently, an active learning potential (ANI-1x) for organic molecules which can produce smooth and physically meaningful PESs has been developed. The high efficiency and accuracy make ANI-1x especially attractive for geometry optimization at low cost. To apply the ANI-1x potential in MMFF parameterization, one needs to perform constrained geometry optimization. In this work, we first developed a computational protocol to constrain rotatable torsional angles and other geometric parameters for a molecule whose geometry is described by Cartesian coordinates. The constraint is successfully achieved by force projection for the two conjugated gradient (CG) algorithms. We then conducted large-scale assessments on ANI-1x along with four different optimization algorithms in reproducing DFT energies and geometries for two CG algorithms, CG backtracking line search (CG-BS) and CG Wolfe line search (CG-WS), and two quasi-Newton algorithms, Broyden-Fletcher-Goldfarb-Shanno (BFGS) and low-memory BFGS (L-BFGS). Note that CG-BS is a new algorithm we developed in this work. All four algorithms take the ANI energies and forces to optimize a molecule geometry. Last, we conducted a large-scale assessment of applying ANI-1x in MMFF development in three aspects. First, we performed full optimizations for 100 drug molecules, each consisting of five distinct conformations. The average root-mean-square error (RMSE) between ANI-1x and DFT is about 1.3 kcal/mol, and the root-mean-square displacement (RMSD) of heavy atoms is about 0.35 Å. Second, we generated torsional PESs for 160 organic molecules, and constrained optimizations were performed for up to 18 conformations for each PES. We found that the RMSE of all the conformers is 1.23 kcal/mol. Last, we carried out constrained optimizations for alanine dipeptide with both ϕ and φ angles being frozen. The Ramachandran plots indicate that the two CG algorithms in conjunction with the ANI-1x potential could well reproduce the DFT-optimized geometries and torsional PESs. We concluded that CG-BS and CG-WS are good choices for generating PESs, while CG-WS or BFGS is ideal for performing full geometry optimization. With the continuously increased quality of ANI, it is expected that the computational algorithms and protocols presented in this work will have great applications in improving the quality of an existing small-molecule MMFF.

Exploring the concerted mechanistic pathway for HIV-1 PR-substrate revealed by umbrella sampling simulation.

HIV-1 protease (HIV-1 PR) is an essential enzyme for the replication process of its virus, and therefore considered an important target for the development of drugs against the acquired immunodeficiency syndrome (AIDS). Our previous study shows that the catalytic mechanism of subtype B/C-SA HIV-1 PR follows a one-step concerted acyclic hydrolysis reaction process using a two-layered ONIOM B3LYP/6-31++G(d,p) method. This present work is aimed at exploring the proposed mechanism of the proteolysis catalyzed by HIV-1 PR and to ensure our proposed mechanism is not an artefact of a single theoretical technique. Hence, we present umbrella sampling method that is suitable for calculating potential mean force (PMF) for non-covalent ligand/substrate-enzyme association/dissociation interactions which provide thermodynamic details for molecular recognition. The free activation energy results were computed in terms of PMF analysis within the hybrid QM(DFTB)/MM approach. The theoretical findings suggest that the proposed mechanism corresponds in principle with experimental data. Given our observations, we suggest that the QM/MM MD method can be used as a reliable computational technique to rationalize lead compounds against specific targets such as the HIV-1 protease.

pH Effects and Cooperativity among Key Titratable Residues for Escherichia coli Glycinamide Ribonucleotide Transformylase.

Human glycinamide ribonucleotide transformylase (GAR Tfase) is a regulatory enzyme in the de novo purine biosynthesis pathway that has been extensively studied as an anticancer target. To some extent, inhibition of GAR Tfase selectively targets cancer cells over normal cells and inhibits purine formation and DNA replication. In this study, we investigated E. coli GAR Tfase, which shares high sequence similarity with the human GAR Tfase, and most functional residues are conserved. Herein, we aim to predict the pH-activity curve through a computational approach. We carried out pH-replica exchange molecular dynamics (pH-REMD) simulations to investigate pH-dependent functions such as structural changes, ligand binding, and catalytic activity. To compute the pH-activity curve, we identified the catalytic residues in specific protonation states, referred to as the catalytic competent protonation states (CCPS), which maintain the structure, keep ligands bound, and facilitate catalysis. Our computed population of CCPS with respect to pH matches well with the experimental pH-activity curve. To compute the microscopic pKa values in the catalytically active conformation, we devised a thermodynamic model that considers the coupling between protonation states of CCPS residues and conformational states. These results allow us to correctly identify the general acid and base catalysts and interpret the pH-activity curve at an atomistic level.

pH-Dependent Conformational Changes Lead to a Highly Shifted pKa for a Buried Glutamic Acid Mutant of SNase.

Ionizable residues are rarely present in the hydrophobic interior of proteins, but when they are, they play important roles in biological processes such as energy transduction and enzyme catalysis. Internal ionizable residues have anomalous experimental pKa values with respect to their pKa in bulk water. This work investigates the atomistic cause of the highly shifted pKa of the internal Glu23 in the artificially mutated variant V23E of Staphylococcal Nuclease (SNase) using pH replica exchange molecular dynamics (pH-REMD) simulations. The pKa of Glu23 obtained from our calculations is 6.55, which is elevated with respect to the glutamate pKa of 4.40 in bulk water. The calculated value is close to the experimental pKa of 7.10. Our simulations show that the highly shifted pKa of Glu23 is the product of a pH-dependent conformational change, which has been observed experimentally and also seen in our simulations. We carry out an analysis of this pH-dependent conformational change in response to the protonation state change of Glu23. Using a four-state thermodynamic model, we estimate the two conformation-specific pKa values of Glu23 and describe the coupling between the conformational and ionization equilibria.

Steered molecular dynamic simulations reveal Marfan syndrome mutations disrupt fibrillin-1 cbEGF domain mechanosensitive calcium binding.

Marfan syndrome (MFS) is a highly variable genetic connective tissue disorder caused by mutations in the calcium binding extracellular matrix glycoprotein fibrillin-1. Patients with the most severe form of MFS (neonatal MFS; nMFS) tend to have mutations that cluster in an internal region of fibrillin-1 called the neonatal region. This region is predominantly composed of eight calcium-binding epidermal growth factor-like (cbEGF) domains, each of which binds one calcium ion and is stabilized by three highly conserved disulfide bonds. Crucially, calcium plays a fundamental role in stabilizing cbEGF domains. Perturbed calcium binding caused by cbEGF domain mutations is thus thought to be a central driver of MFS pathophysiology. Using steered molecular dynamics (SMD) simulations, we demonstrate that cbEGF domain calcium binding decreases under mechanical stress (i.e. cbEGF domains are mechanosensitive). We further demonstrate the disulfide bonds in cbEGF domains uniquely orchestrate protein unfolding by showing that MFS disulfide bond mutations markedly disrupt normal mechanosensitive calcium binding dynamics. These results point to a potential mechanosensitive mechanism for fibrillin-1 in regulating extracellular transforming growth factor beta (TGFB) bioavailability and microfibril integrity. Such mechanosensitive "smart" features may represent novel mechanisms for mechanical hemostasis regulation in extracellular matrix that are pathologically activated in MFS.

Selective ring-rearrangement or ring-closing metathesis of bicyclo[3.2.1]octenes.

Explored was the competitive ring-closing metathesis vs. ring-rearrangement metathesis of bicyclo[3.2.1]octenes prepared by a simple and convergent synthesis from bicyclic alkylidenemalono-nitriles and allylic electrophiles. It was uncovered that ring-closing metathesis occurs exclusively on the tetraene-variant, yielding unique, stereochemically and functionally rich polycyclic bridged frameworks, whereas the reduced version (a triene) undergoes ring-rearrangement metathesis to 5-6-5 fused ring systems resembling the isoryanodane core.

NEXMD Software Package for Nonadiabatic Excited State Molecular Dynamics Simulations.

We present a versatile new code released for open community use, the nonadiabatic excited state molecular dynamics (NEXMD) package. This software aims to simulate nonadiabatic excited state molecular dynamics using several semiempirical Hamiltonian models. To model such dynamics of a molecular system, the NEXMD uses the fewest-switches surface hopping algorithm, where the probability of transition from one state to another depends on the strength of the derivative nonadiabatic coupling. In addition, there are a number of algorithmic improvements such as empirical decoherence corrections and tracking trivial crossings of electronic states. While the primary intent behind the NEXMD was to simulate nonadiabatic molecular dynamics, the code can also perform geometry optimizations, adiabatic excited state dynamics, and single-point calculations all in vacuum or in a simulated solvent. In this report, first, we lay out the basic theoretical framework underlying the code. Then we present the code's structure and workflow. To demonstrate the functionality of NEXMD in detail, we analyze the photoexcited dynamics of a polyphenylene ethynylene dendrimer (PPE, C30H18) in vacuum and in a continuum solvent. Furthermore, the PPE molecule example serves to highlight the utility of the getexcited.py helper script to form a streamlined workflow. This script, provided with the package, can both set up NEXMD calculations and analyze the results, including, but not limited to, collecting populations, generating an average optical spectrum, and restarting unfinished calculations.