Training machine learning potentials for reactive systems: A Colab tutorial on basic models.

In the last several years, there has been a surge in the development of machine learning potential (MLP) models for describing molecular systems. We are interested in a particular area of this field - the training of system-specific MLPs for reactive systems - with the goal of using these MLPs to accelerate free energy simulations of chemical and enzyme reactions. To help new members in our labs become familiar with the basic techniques, we have put together a self-guided Colab tutorial (https://cc-ats.github.io/mlp_tutorial/), which we expect to be also useful to other young researchers in the community. Our tutorial begins with the introduction of simple feedforward neural network (FNN) and kernel-based (using Gaussian process regression, GPR) models by fitting the two-dimensional Müller-Brown potential. Subsequently, two simple descriptors are presented for extracting features of molecular systems: symmetry functions (including the ANI variant) and embedding neural networks (such as DeepPot-SE). Lastly, these features will be fed into FNN and GPR models to reproduce the energies and forces for the molecular configurations in a Claisen rearrangement reaction.

Potent inhibitors targeting cyclin-dependent kinase 9 discovered via virtual high-throughput screening and absolute binding free energy calculations.

Due to the crucial regulatory mechanism of cyclin-dependent kinase 9 (CDK9) in mRNA transcription, the development of kinase inhibitors targeting CDK9 holds promise as a potential treatment strategy for cancer. A structure-based virtual screening approach has been employed for the discovery of potential novel CDK9 inhibitors. First, compounds with kinase inhibitor characteristics were identified from the ZINC15 database via virtual high-throughput screening. Next, the predicted binding modes were optimized by molecular dynamics simulations, followed by precise estimation of binding affinities using absolute binding free energy calculations based on the free energy perturbation scheme. The binding mode of molecule 006 underwent an inward-to-outward flipping, and the new binding mode exhibited binding affinity comparable to the small molecule T6Q in the crystal structure (PDB ID: 4BCF), highlighting the essential role of molecular dynamics simulation in capturing a plausible binding pose bridging docking and absolute binding free energy calculations. Finally, structural modifications based on these findings further enhanced the binding affinity with CDK9. The results revealed that enhancing the molecule's rigidity through ring formation, while maintaining the major interactions, reduced the entropy loss during the binding process and, thus, enhanced binding affinities.

Hydrogen Abstraction by Alkoxyl Radicals: Computational Studies of Thermodynamic and Polarity Effects on Reactivities and Selectivities.

Density functional theory calculations (ωB97X-D) are reported for the reactions of methoxy, tert-butoxy, trichloroethoxy, and trifluoroethoxy radicals with a series of 26 C-H bonds in different environments characteristic of a variety of hydrocarbons and substituted derivatives. The variations in activation barriers are analyzed with modified Evans-Polanyi treatments to account for polarity and unsaturation effects. The treatments by Roberts and Steel and by Mayer have inspired the development of a simple treatment involving the thermodynamics of reactions, the difference between the reactant radical and product radical electronegativities, and the absence or presence of α-unsaturation. The three-parameter equation (ΔH⧧ = 0.52ΔHrxn(1 - d) - 0.35ΔχAB2 + 10.0, where d = 0.44 when there is α-unsaturation to the reacting C-H bond), correlates well with quantum mechanically computed barriers and shows the quantitative importance of the thermodynamics of reactions (dictated by the reactant and the product bond dissociation energies) and polar effects.

Near Infrared Light-Driven Photothermal Effect on Homochiral Au/TiO2 Nanotube Arrays for Enantioselective Desorption.

Chiral enantiomers have different effects on biological processes. Enantiomer separation is significant and necessary. Herein, a photothermal (PT) effect-derived enantioselective desorption strategy based on homochiral Au/TiO2 nanotubes (NTs) is developed. Using 3,4-dihydroxyphenylalanine (DOPA) as the model enantiomer, an obvious selective desorption of L/D-DOPA can be achieved by the NIR light-triggered local temperature enhancement. Molecular docking simulation further verifies that the distinct affinity precipitated by the different hydrogen bonds between homochiral sorbent and target enantiomers is the origin of enantioselective desorption. This desorption strategy provides a green and alternative approach for the selective separation of chiral molecules.

Free energy change estimation: The Divide and Conquer MBAR method.

In the present study, the Divide and Conquer MBAR (DC-MBAR) method is proposed to predict the free energies based on the data sampled by multi-states simulations. For DC-MBAR method, the overlap between any two alchemical states is calculated first and those with sufficient overlap are defined as the adjacent states. Unlike the traditional MBAR method, which calculates the free energy of each state using all the data at once, DC-MBAR focuses on predicting the free energy changes between adjacent states. To estimate the free energy changes accurately, the other states with overlaps with the two adjacent states bigger than the defined threshold are included in the MBAR equation. At a specific threshold, the free energies predicted by DC-MBAR are very close to those calculated by the traditional MABR method. Furthermore, DC-MBAR scheme can reduce both the computation and memory cost. One important characteristic of DC-MBAR method is linear scaling, which means the CPU time with the change of the number of states is a straight-line relation. As the pair-based calculations are mutually independent and parallelizable, all accessible CPU cores on the HPC cluster could be utilized, which makes DC-MBAR strategy more efficient.

Atomic structures and deletion mutant reveal different capsid-binding patterns and functional significance of tegument protein pp150 in murine and human cytomegaloviruses with implications for therapeutic development.

Cytomegalovirus (CMV) infection causes birth defects and life-threatening complications in immunosuppressed patients. Lack of vaccine and need for more effective drugs have driven widespread ongoing therapeutic development efforts against human CMV (HCMV), mostly using murine CMV (MCMV) as the model system for preclinical animal tests. The recent publication (Yu et al., 2017, DOI: 10.1126/science.aam6892) of an atomic model for HCMV capsid with associated tegument protein pp150 has infused impetus for rational design of novel vaccines and drugs, but the absence of high-resolution structural data on MCMV remains a significant knowledge gap in such development efforts. Here, by cryoEM with sub-particle reconstruction method, we have obtained the first atomic structure of MCMV capsid with associated pp150. Surprisingly, the capsid-binding patterns of pp150 differ between HCMV and MCMV despite their highly similar capsid structures. In MCMV, pp150 is absent on triplex Tc and exists as a "Λ"-shaped dimer on other triplexes, leading to only 260 groups of two pp150 subunits per capsid in contrast to 320 groups of three pp150 subunits each in a "Δ"-shaped fortifying configuration. Many more amino acids contribute to pp150-pp150 interactions in MCMV than in HCMV, making MCMV pp150 dimer inflexible thus incompatible to instigate triplex Tc-binding as observed in HCMV. While pp150 is essential in HCMV, our pp150-deletion mutant of MCMV remained viable though with attenuated infectivity and exhibiting defects in retaining viral genome. These results thus invalidate targeting pp150, but lend support to targeting capsid proteins, when using MCMV as a model for HCMV pathogenesis and therapeutic studies.

Inhibitor of DNA-binding family regulates the prognosis of ovarian cancer.

Aim: The mechanistic role of inhibitor of DNA binding or differentiation (ID) family in ovarian cancer (OC) has remained unclear. Materials & methods: We used the Oncomine, GEPIA, Kaplan-Meier Plotter, cBioPortal, SurvExpress, PROGgene V2, TIMER, and FunRich to evaluate the prognostic value of IDs in patients with OC. Results: the mRNA transcripts of all IDs were markedly downregulated in OC compared with normal tissue. The prognostic value of IDs was also explored within the subtypes, pathological stages, clinical stages and TP53 mutational status. The group with low-risk IDs showed relatively good overall survival (OS) compared with the high-risk group. Conclusion: ID1/3/4 may be exploited as promising prognostic biomarkers and therapeutic targets in OC patients.

Affordable Ab Initio Path Integral for Thermodynamic Properties via Molecular Dynamics Simulations Using Semiempirical Reference Potential.

Path integral molecular dynamics (PIMD) is becoming a routinely applied method for incorporating the nuclear quantum effect in computer simulations. However, direct PIMD simulations at an ab initio level of theory are formidably expensive. Using the protonated 1,8-bis(dimethylamino)naphthalene molecule as an example, we show in this work that the computational expense for the intramolecular proton transfer between the two nitrogen atoms can be remarkably reduced by implementing the idea of reference-potential methods. The simulation time can be easily extended to a scale of nanoseconds while maintaining the accuracy on an ab initio level of theory for thermodynamic properties. In addition, postprocessing can be carried out in parallel on massive computer nodes. A 545-fold reduction in the total CPU time can be achieved in this way as compared to a direct PIMD simulation at the same ab initio level of theory.

Phase-error-compensation-based surface recovery algorithm using spectrum selection for white light interferometry.

White light interferometry is a well-established surface recovery technique. In this paper, a white light signal processing algorithm based on phase error compensation using spectrum selection is proposed. The derived nonlinear phase distribution from the correlogram is modeled as the combination of random errors and systemic deviations. By developing a new, to the best of our knowledge, recovery algorithm, the phase noise can be separated from the linear map and significantly attenuated. Based on the proposed algorithm, the spectrum features of white light LEDs and halogen lamps are investigated in detail. The inner products defined by three selected points are employed to generate a coefficient to evaluate the linearity of an unwrapped phase map within a certain spectrum region. The optimal spectrum range corresponding to the best measurement performance can then be located where the coefficient approximates 1 and the spectrum energy stays relatively high. The simulations are carried out under different levels of SNR and scan step noises, which show that the new method can effectively reduce additional disturbance from the recovered topography. In experiments, the system with the proposed method is first calibrated by a step height standard (VLSI, 182.7±2.0nm) with the repeatability of 0.44 nm. A silicon wafer and three roughness standards are also tested to further verify the robustness of the new method.

Reproducing the low-temperature excitation energy transfer dynamics of phycoerythrin 545 light-harvesting complex with a structure-based model Hamiltonian.

To elucidate the energy transfer mechanism of the PE545 light-harvesting complex, an exciton model is constructed with the full Hamiltonian obtained from structure-based calculations. The electronic couplings and spectral densities are evaluated on the basis of the site energies and transition dipole moments obtained from our recent Molecular Dynamics-Quantum Mechanical/Molecular Mechanical (MD-QM/MM) study [Tong et al., J. Phys. Chem. B 123, 2040-2049 (2019)]. The polarized protein-specific charge model is employed both in the MD simulation and in the QM/MM calculations to account for the environmental fluctuation of the protein scaffold. The energy transfer pathways are, thus, derived, which agree well with the phenomenological models based on the spatial organization of the chromophores and the experimental observations. Moreover, the simulated linear absorption spectra using the dissipaton equation of motion approach agree well with the experimental ones, and the resulting population dynamics indicates that an optimal energy transfer efficiency is reproduced.

Accelerated computation of free energy profile at ab initio quantum mechanical/molecular mechanical accuracy via a semi-empirical reference potential. II. Recalibrating semi-empirical parameters with force matching.

An efficient and accurate reference potential simulation protocol is proposed for producing ab initio quantum mechanical/molecular mechanical (AI-QM/MM) quality free energy profiles for chemical reactions in a solvent or macromolecular environment. This protocol involves three stages: (a) using force matching to recalibrate a semi-empirical quantum mechanical (SE-QM) Hamiltonian for the specific reaction under study; (b) employing the recalibrated SE-QM Hamiltonian (in combination with molecular mechanical force fields) as the reference potential to drive umbrella samplings along the reaction pathway; and (c) computing AI-QM/MM energy values for collected configurations from the sampling and performing weighted thermodynamic perturbation to acquire an AI-QM/MM corrected reaction free energy profile. For three model reactions (identity SN2 reaction, Menshutkin reaction, and glycine proton transfer reaction) in aqueous solution and one enzyme reaction (Claisen arrangement in chorismate mutase), our simulations using recalibrated PM3 SE-QM Hamiltonians well reproduced QM/MM free energy profiles at the B3LYP/6-31G* level of theory all within 1 kcal mol-1 with a 20 to 45 fold reduction in the computer time.

Brain temporal complexity in explaining the therapeutic and cognitive effects of seizure therapy.

Over 350 million people worldwide suffer from depression, a third of whom are medication-resistant. Seizure therapy remains the most effective treatment in depression, even when many treatments fail. The utility of seizure therapy is limited due to its cognitive side effects and stigma. The biological targets of seizure therapy remain unknown, hindering design of new treatments with comparable efficacy. Seizures impact the brains temporal dynamicity observed through electroencephalography. This dynamicity reflects richness of information processing across distributed brain networks subserving affective and cognitive processes. We investigated the hypothesis that seizure therapy impacts mood (depressive symptoms) and cognition by modulating brain temporal dynamicity. We obtained resting-state electroencephalography from 34 patients (age = 46.0 ± 14.0, 21 females) receiving two types of seizure treatments-electroconvulsive therapy or magnetic seizure therapy. We used multi-scale entropy to quantify the complexity of the brain's temporal dynamics before and after seizure therapy. We discovered that reduction of complexity in fine timescales underlined successful therapeutic response to both seizure treatments. Greater reduction in complexity of fine timescales in parieto-occipital and central brain regions was significantly linked with greater improvement in depressive symptoms. Greater increase in complexity of coarse timescales was associated with greater decline in cognition including the autobiographical memory. These findings were region and timescale specific. That is, change in complexity in occipital regions (e.g. O2 electrode or right occipital pole) at fine timescales was only associated with change in depressive symptoms, and not change in cognition, and change in complexity in parieto-central regions (e.g. Pz electrode or intra and transparietal sulcus) at coarser timescale was only associated with change in cognition, and not depressive symptoms. Finally, region and timescale specific changes in complexity classified both antidepressant and cognitive response to seizure therapy with good (80%) and excellent (95%) accuracy, respectively. In this study, we discovered a novel biological target of seizure therapy: complexity of the brain resting state dynamics. Region and timescale dependent changes in complexity of the brain resting state dynamics is a novel mechanistic marker of response to seizure therapy that explains both the antidepressant response and cognitive changes associated with this treatment. This marker has tremendous potential to guide design of the new generation of antidepressant treatments.

Efficient Computation of Free Energy Surfaces of Diels⁻Alder Reactions in Explicit Solvent at Ab Initio QM/MM Level.

For Diels⁻Alder (DA) reactions in solution, an accurate and converged free energy (FE) surface at ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) level is imperative for the understanding of reaction mechanism. However, this computation is still far too expensive. In a previous work, we proposed a new method termed MBAR+wTP, with which the computation of the ai FE profile can be accelerated by several orders of magnitude via a three-step procedure: (I) an umbrella sampling (US) using a semi-empirical (SE) QM/MM Hamiltonian is performed; (II) the FE profile is generated using the Multistate Bennett Acceptance Ratio (MBAR) analysis; and (III) a weighted Thermodynamic Perturbation (wTP) from the SE Hamiltonian to the ai Hamiltonian is performed to obtain the ai QM/MM FE profile using weight factors from the MBAR analysis. In this work, this method is extended to the calculations of two-dimensional FE surfaces of two Diels⁻Alder reactions of cyclopentadiene with either acrylonitrile or 1-4-naphthoquinone at ai QM/MM level. The accurate activation free energies at the ai QM/MM level, which are much closer to the experimental measurements than those calculated by other methods, indicate that this MBAR+wTP method can be applied in the studies of complex reactions in condensed phase with much-enhanced efficiency.

Bimodal Evans-Polanyi Relationships in Dioxirane Oxidations of sp3 C-H: Non-perfect Synchronization in Generation of Delocalized Radical Intermediates.

The selectivities in C-H oxidations of a variety of compounds by DMDO have been explored with density functional theory. There is a linear Evans-Polanyi-type correlation for saturated substrates. Activation energies correlate with reaction energies or, equivalently, BDEs (ΔH‡sat = 0.91*BDE - 67.8). Unsaturated compounds, such as alkenes, aromatics, and carbonyls, exhibit a different correlation for allylic and benzylic C-H bonds (ΔH‡unsat = 0.35*BDE - 13.1). Bernasconi's Principle of Non-Perfect Synchronization (NPS) is found to operate here. The origins of this phenomenon were analyzed by a Distortion/Interaction model. Computations indicate early transition states for H-abstractions from allylic and benzylic C-H bonds, but later transition states for the saturated. The reactivities are mainly modulated by the distortion energy and the degree of dissociation of the C-H bond. While the increase in barrier with higher BDE is not unexpected from the Evans-Polanyi relationship, two separate correlations, one for saturated compounds, and one for unsaturated leading to delocalized radicals, were unexpected.

Evaluation of the Coupled Two-Dimensional Main Chain Torsional Potential in Modeling Intrinsically Disordered Proteins.

Intrinsically disordered proteins (IDPs) carry out crucial biological functions in essential biological processes of life. Because of the highly dynamic and conformationally heterogeneous nature of the disordered states of IDPs, molecular dynamics simulations are becoming an indispensable tool for the investigation of the conformational ensembles and dynamic properties of IDPs. Nevertheless, there is still no consensus on the most reliable force field in molecular dynamics simulations for IDPs hitherto. In this work, the recently proposed AMBER99SB2D force field is evaluated in modeling some disordered polypeptides and proteins by checking its ability to reproduce experimental NMR data. The results highlight that when the ildn side-chain corrections are included, AMBER99SB2D-ildn exhibits reliable results that agree with experiments compared with its predecessors, the AMBER14SB, AMBER99SB, AMBER99SB-ildn, and AMBER99SB2D force fields, and that decreasing the overall magnitude of protein-protein interactions in favor of protein-water interactions is a key ingredient behind the improvement.

Efficient Strategy for the Calculation of Solvation Free Energies in Water and Chloroform at the Quantum Mechanical/Molecular Mechanical Level.

The partitioning of solute molecules between immiscible solvents with significantly different polarities is of great importance. The polarization between the solute and solvent molecules plays an essential role in determining the solubility of the solute, which makes computational studies utilizing molecular mechanics (MM) rather difficult. In contrast, quantum mechanics (QM) can provide more reliable predictions. In this work, the partition coefficients of the side chain analogs of some amino acids between water and chloroform were computed. The QM solvation free energies were calculated indirectly via a series of MM states using the multistate Bennett acceptance ratio (MBAR) and the MM-to-QM corrections were applied at the two endpoints using thermodynamic perturbation (TP). Previously, it has been shown (Jia et al. J. Chem. Theory Comput. 2016, 12, 499-511) that this method provides the minimal variance in the results without running QM simulations. However, if there is insufficient overlap in phase space between the MM and QM Hamiltonians, this method fails. In this work, we propose, for the first time, a quantity termed the reweighting entropy that serves as a metric for the reliability of the TP calculations. If the reweighting entropy is below a certain threshold (0.65 for the solvation free energy calculations in this work), this MM-to-QM correction should be avoided and two alternative methods can be employed by either introducing a semiempirical state or conducting nonequilibrium simulations. However, the results show that the QM methods are not guaranteed to yield better results than the MM methods. Further improvement of the QM methods are imperative, especially the treatment of the van der Waals and the electrostatic interactions between the QM region and the MM region in the first shell. We also propose a scheme for the calculation of the van der Waals parameters for the solute molecules in nonaqueous solvent, which improves the quality of the computed thermodynamic properties. Furthermore, the force field parameters for the sulfur-containing molecules are also optimized.

Some practical approaches to treating electrostatic polarization of proteins.

Conspectus Electrostatic interaction plays a significant role in determining many properties of biomolecules, which exist and function in aqueous solution, a highly polar environment. For example, proteins are composed of amino acids with charged, polar, and nonpolar side chains and their specific electrostatic properties are fundamental to the structure and function of proteins. An important issue that arises in computational study of biomolecular interaction and dynamics based on classical force field is lack of polarization. Polarization is a phenomenon in which the charge distribution of an isolated molecule will be distorted when interacting with another molecule or presented in an external electric field. The distortion of charge distribution is intended to lower the overall energy of the molecular system, which is counter balanced by the increased internal energy of individual molecules due to the distorted charge distributions. The amount of the charge redistribution, which characterizes the polarizability of a molecule, is determined by the level of the charge distortion. Polarization is inherently quantum mechanical, and therefore classical force fields with fixed atomic charges are incapable of capturing this important effect. As a result, simulation studies based on popular force fields, AMBER, CHARMM, etc., lack the polarization effect, which is a widely known deficiency in most computational studies of biomolecules today. Many efforts have been devoted to remedy this deficiency, such as adding additional movable charge on the atom, allowing atomic charges to fluctuate, or including induced multipoles. Although various successes have been achieved and progress at various levels has been reported over the past decades, the issue of lacking polarization in force field based simulations is far from over. For example, some of these methods do not always give converged results, and other methods require huge computational cost. This Account reviews recent work on developing polarized and polarizable force fields based on fragment quantum mechanical calculations for proteins. The methods described here are based on quantum mechanical calculations of proteins in solution, but with a different level of rigor and different computational efficiency for the molecular dynamics applications. In the general approach, a fragment quantum mechanical calculation for protein with implicit solvation is carried out to derive a polarized protein-specific charge (PPC) for any given protein structure. The PPC correctly reflects the polarization state of the protein in a given conformation, and it can also be dynamically changed as the protein changes conformation in dynamics simulations. Another approach that is computationally more efficient is the effective polarizable bond method in which only polar bonds or groups can be polarized and their polarizabilities are predetermined from quantum mechanical calculations of these groups in external electric fields. Both methods can be employed for applications in various situations by taking advantage of their unique features.

Multiple environment single system quantum mechanical/molecular mechanical (MESS-QM/MM) calculations. 1. Estimation of polarization energies.

In combined quantum mechanical/molecular mechanical (QM/MM) free energy calculations, it is often advantageous to have a frozen geometry for the quantum mechanical (QM) region. For such multiple-environment single-system (MESS) cases, two schemes are proposed here for estimating the polarization energy: the first scheme, termed MESS-E, involves a Roothaan step extrapolation of the self-consistent field (SCF) energy; whereas the other scheme, termed MESS-H, employs a Newton-Raphson correction using an approximate inverse electronic Hessian of the QM region (which is constructed only once). Both schemes are extremely efficient, because the expensive Fock updates and SCF iterations in standard QM/MM calculations are completely avoided at each configuration. They produce reasonably accurate QM/MM polarization energies: MESS-E can predict the polarization energy within 0.25 kcal/mol in terms of the mean signed error for two of our test cases, solvated methanol and solvated ß-alanine, using the M06-2X or ωB97X-D functionals; MESS-H can reproduce the polarization energy within 0.2 kcal/mol for these two cases and for the oxyluciferin-luciferase complex, if the approximate inverse electronic Hessians are constructed with sufficient accuracy.

Numerical study on the partitioning of the molecular polarizability into fluctuating charge and induced atomic dipole contributions.

In order to carry out a detailed analysis of the molecular static polarizability, which is the response of the molecule to a uniform external electric field, the molecular polarizability was computed using the finite-difference method for 21 small molecules, using density functional theory. Within nine charge population schemes (Löwdin, Mulliken, Becke, Hirshfeld, CM5, Hirshfeld-I, NPA, CHELPG, MK-ESP) in common use, the charge fluctuation contribution is found to dominate the molecular polarizability, with its ratio ranging from 59.9% with the Hirshfeld or CM5 scheme to 96.2% with the Mulliken scheme. The Hirshfeld-I scheme is also used to compute the other contribution to the molecular polarizability coming from the induced atomic dipoles, and the atomic polarizabilities in eight small molecules and water pentamer are found to be highly anisotropic for most atoms. Overall, the results suggest that (a) more emphasis probably should be placed on the charge fluctuation terms in future polarizable force field development and (b) an anisotropic polarizability might be more suitable than an isotropic one in polarizable force fields based entirely or partially on the induced atomic dipoles.

Accessing the applicability of polarized protein-specific charge in linear interaction energy analysis.

The reliability of the linear interaction energy (LIE) depends on the atomic charge model used to delineate the Coulomb interaction between the ligand and its environment. In this work, the polarized protein-specific charge (PPC) implementing a recently proposed fitting scheme has been examined in the LIE calculations of the binding affinities for avidin and ß-secretase binding complexes. This charge fitting scheme, termed delta restrained electrostatic potential, bypasses the prevalent numerical difficulty of rank deficiency in electrostatic-potential-based charge fitting methods via a dual-step fitting strategy. A remarkable consistency between the predicted binding affinities and the experimental measurement has been observed. This work serves as a direct evidence of PPC's applicability in rational drug design.