Deep reinforcement learning of transition states.

Combining reinforcement learning (RL) and molecular dynamics (MD) simulations, we propose a machine-learning approach, called RL‡, to automatically unravel chemical reaction mechanisms. In RL‡, locating the transition state of a chemical reaction is formulated as a game, and two functions are optimized, one for value estimation and the other for policy making, to iteratively improve our chance of winning this game. Both functions can be approximated by deep neural networks. By virtue of RL‡, one can directly interpret the reaction mechanism according to the value function. Meanwhile, the policy function allows efficient sampling of the transition path ensemble, which can be further used to analyze reaction dynamics and kinetics. Through multiple experiments, we show that RL‡ can be trained tabula rasa hence allowing us to reveal chemical reaction mechanisms with minimal subjective biases.

Comparison of the Microsolvation of CaX2 (X = F, Cl, Br, I) in Water: Size-Selected Anion Photoelectron Spectroscopy and Theoretical Calculations.

To understand the microsolvation of alkaline-earth dihalides in water and provide information about the dependence of solvation processes on different halides, we investigated CaBr2(H2O)n-, CaI2(H2O)n-, and CaF2(H2O)n- (n = 0-6) clusters using size-selected anion photoelectron spectroscopy and conducted theoretical calculations on these clusters and their neutrals. The results are compared with those of CaCl2(H2O)n-/0 clusters reported previously. It is found that the vertical detachment energies (VDEs) of CaCl2(H2O)n-, CaBr2(H2O)n-, and CaI2(H2O)n- show a similar trend with increasing cluster size, while the VDEs of CaF2(H2O)n- show a different trend. The VDEs of CaF2(H2O)n- are much lower than those of CaCl2(H2O)n-, CaBr2(H2O)n-, and CaI2(H2O)n-. A detailed probing of the structures shows that a significant increase of the Ca-X distance (separation of Ca2+-X- ion pair) in CaCl2(H2O)n-/0, CaBr2(H2O)n-/0, and CaI2(H2O)n-/0 clusters occurred at about n = 5. However, for CaF2(H2O)n-/0, no abrupt change of the Ca-F distance with the increasing cluster size has been observed. In CaCl2(H2O)6-/0, CaBr2(H2O)6-/0, and CaI2(H2O)6-/0, the Ca atom coordinates directly with 5 H2O molecules. However, in CaF2(H2O)n-/0, the Ca atom coordinates directly with only 2 or 3 H2O molecules. The similarity or differences in the structures and coordination numbers are consistent with the fact that CaCl2, CaBr2, and CaI2 have similar solubility, while CaF2 has much lower solubility.

A Perspective on Deep Learning for Molecular Modeling and Simulations.

Deep learning is transforming many areas in science, and it has great potential in modeling molecular systems. However, unlike the mature deployment of deep learning in computer vision and natural language processing, its development in molecular modeling and simulations is still at an early stage, largely because the inductive biases of molecules are completely different from those of images or texts. Footed on these differences, we first reviewed the limitations of traditional deep learning models from the perspective of molecular physics and wrapped up some relevant technical advancement at the interface between molecular modeling and deep learning. We do not focus merely on the ever more complex neural network models; instead, we introduce various useful concepts and ideas brought by modern deep learning. We hope that transacting these ideas into molecular modeling will create new opportunities. For this purpose, we summarized several representative applications, ranging from supervised to unsupervised and reinforcement learning, and discussed their connections with the emerging trends in deep learning. Finally, we give an outlook for promising directions which may help address the existing issues in the current framework of deep molecular modeling.

Hydration processes of barium chloride: Size-selected anion photoelectron spectroscopy and theoretical calculations of BaCl2-water clusters.

In order to understand the hydration processes of BaCl2, we investigated BaCl2(H2O)n - (n = 0-5) clusters using size-selected anion photoelectron spectroscopy and theoretical calculations. The structures of neutral BaCl2(H2O)n clusters up to n = 8 were also investigated by theoretical calculations. It is found that in BaCl2(H2O)n -/0, the Ba-Cl distances increase very slowly with the cluster size. The hydration process is not able to induce the breaking of a Ba-Cl bond in the cluster size range (n = 0-8) studied in this work. In small BaCl2(H2O)n clusters with n ≤ 5, the Ba atom has a coordination number of n + 2; however, in BaCl2(H2O)6-8 clusters, the Ba atom coordinates with two Cl atoms and (n - 1) water molecules, and it has a coordination number of n + 1. Unlike the previously studied MgCl2(H2O)n - and CaCl2(H2O)n -, negative charge-transfer-to-solvent behavior has not been observed for BaCl2(H2O)n -, and the excess electron of BaCl2(H2O)n - is mainly localized on the Ba atom rather on the water molecules. No observation of Ba2+-Cl- separation in current work is consistent with the lower solubility of BaCl2 compared to MgCl2 and CaCl2. Considering the BaCl2/H2O mole ratio in the saturated solution, one would expect that about 20-30 H2O molecules are needed to break the first Ba-Cl bond in BaCl2.

From 1D sequence to 3D chromatin dynamics and cellular functions: a phase separation perspective.

The high-order chromatin structure plays a non-negligible role in gene regulation. However, the mechanism, especially the sequence dependence for the formation of varied chromatin structures in different cells remains to be elucidated. As the nucleotide distributions in human and mouse genomes are highly uneven, we identified CGI (CpG island) forest and prairie genomic domains based on CGI densities of a species, dividing the genome into two sequentially, epigenetically, and transcriptionally distinct regions. These two megabase-sized domains also spatially segregate to different extents in different cell types. Forests and prairies show enhanced segregation from each other in development, differentiation, and senescence, meanwhile the multi-scale forest-prairie spatial intermingling is cell-type specific and increases in differentiation, helping to define cell identity. We propose that the phase separation of the 1D mosaic sequence in space serves as a potential driving force, and together with cell type specific epigenetic marks and transcription factors, shapes the chromatin structure in different cell types. The mosaicity in genome of different species in terms of forests and prairies could relate to observations in their biological processes like development and aging. In this way, we provide a bottoms-up theory to explain the chromatin structural and epigenetic changes in different processes.

The Exploration of a New Stable G-Triplex DNA and Its Novel Function in Electrochemical Biosensing.

A G-triplex, a new kind of DNA structure, has been identified as an intermediate in the folding of G-quadruplexes. However, the studies on G-triplexes are still very limited, and the functions and applications of G-triplexes need to be further developed. In this paper, a new G-triplex sequence (5'-CTGGGAGGGAGGGA-3', G3), obtained by truncating four bases (GGGA) from the 3' end of an 18-base G-quadruplex sequence (G4), was found to significantly decrease the diffusion current of methylene blue (MB). In particular, we proved that (a) MB stabilized the structure of G3 and increased the Tm of G3 considerably based on circular dichroism; and (b) MB formed a 1:1 noncovalent complex with G3 based on electrospray ionization mass spectrometry. Moreover, molecular dynamics simulations established reliable speculation in the folding topology of G3 and interaction sites between G3 and MB. Based on the strong affinity of G3 with MB, we further developed a novel function of G3 as an electrochemical signal read-out and applied it in the fabrication of a sensitive homogeneous electrochemical aptasensor for cocaine. The features we observed in the G3/MB complex will serve as a new inspiring guideline for developing functional short G-rich ligands.

Selective enhanced sampling in dihedral energy facilitates overcoming the dihedral energy increase in protein folding and accelerates the searching for protein native structure.

The dihedral energy function is the most influential parameter in molecular mechanics (MM) force field parameter optimization. A selective enhanced sampling of dihedral energy could effectively reflect the influence of dihedral energy settings on protein secondary structure representation, which in turn testifies the availability of the force field in folding simulation. Here, a Dihedral-based Selective Integrated-Tempering-Sampling Molecular Dynamics (D-SITSMD) simulation method is shown to provide a selective enhanced sampling of dihedral energy without introducing large energetic noise. Its capabilities of searching for protein natively folded structure and providing the underlying folding pathway are evaluated through the folding tests of three peptides (chignolin, TC5b, and HP35) with multiple AMBER force fields (FF14SBonlysc, FF99SBildn, or FF03) and the comparison to presented experimental data and REMD simulations. Both above-mentioned capabilities are improved, displaying the potential of D-SITSMD in the studies of in silico protein folding and structure refinement. Additionally, it is commonly observed among the test simulation systems that their folding processes are thermodynamically favorable for non-bonded vdW and electrostatic energies but unfavorable for dihedral energy, such that the folding barrier height correlated with the dihedral energy increase from the unfolded to folded state whereas the unfolding free energy barrier correlated with the combined increase of vdW and electrostatic energies in the unfolding process. It is speculated that the influence of a force field on the folding barrier of a protein is fulfilled mainly through regulating the contribution of dihedral energy to determine the secondary structure formation in the global folding process.

Enhanced sampling in molecular dynamics.

Although molecular dynamics simulations have become a useful tool in essentially all fields of chemistry, condensed matter physics, materials science, and biology, there is still a large gap between the time scale which can be reached in molecular dynamics simulations and that observed in experiments. To address the problem, many enhanced sampling methods were introduced, which effectively extend the time scale being approached in simulations. In this perspective, we review a variety of enhanced sampling methods. We first discuss collective-variables-based methods including metadynamics and variationally enhanced sampling. Then, collective variable free methods such as parallel tempering and integrated tempering methods are presented. At last, we conclude with a brief introduction of some newly developed combinatory methods. We summarize in this perspective not only the theoretical background and numerical implementation of these methods but also the new challenges and prospects in the field of the enhanced sampling.

Structure of water confined between two parallel graphene plates.

We study, in this paper, the physical properties of water confined between two parallel graphene plates with different slit widths to understand the effects of confinement on the water structure and how bulk properties are reached as the water layer thickens. It was found that the microscopic structures of the interfacial liquid layer close to graphene vary with the slit width. Water tends to locate at the center of the six-membered ring of graphene planes to form triangular patterns, as found by others. The narrower the slit width is, the more pronounced this pattern is, except for the slit width of 9.5 Å, for which a well-defined two-layer structure of water forms. On the other hand, squared structures can be clearly seen in single snapshots at small (6.5 Å and 7.5 Å) but not large slit widths. Even at small slit widths, the square-like geometry is observed only when an average is taken for a short trajectory, and averaging over a long time yields a triangular pattern dictated by the graphene geometry. We estimate the length of time needed to observe two patterns, respectively. We also used the two-phase thermodynamic model to study the variation of entropy of confined water and found that at 8.5 Å, the entropy of confined water is larger than that of bulk water. The rotational entropy of confined water is higher than that of bulk water for all slit widths due to the reduction of the hydrogen bond in the confined space.

Salt- and pH-Triggered Helix-Coil Transition of Ionic Polypeptides under Physiology Conditions.

Controlling the helix-coil transition of polypeptides under physiological conditions is an attractive way toward smart functional materials. Here, we report the synthesis of a series of tertiary amine-functionalized ethylene glycol (EG x)-linked polypeptide electrolytes with their secondary structures tunable under physiological conditions. The resultant polymers, denoted as P(EG xDMA-Glu) ( x = 1, 2, and 3), show excellent aqueous solubility (>20 mg/mL) regardless of their charge states. Unlike poly-l-lysine that can form a helix only at pH above 10, P(EG xDMA-Glu) undergo a pH-dependent helix-coil switch with their transition points within the physiological range (pH ∼5.3-6.5). Meanwhile, P(EG xDMA-Glu) exhibit an unusual salt-induced helical conformation presumably owing to the unique properties of EG x linkers. Together, the current work highlights the importance of fine-tuning the linker chemistry in achieving conformation-switchable polypeptides and represents a facile approach toward stimuli-responsive biopolymers for advanced biological applications.

Molecular dynamics simulation, ab initio calculation, and size-selected anion photoelectron spectroscopy study of initial hydration processes of calcium chloride.

To understand the initial hydration processes of CaCl2, we performed molecular simulations employing the force field based on the theory of electronic continuum correction with rescaling. Integrated tempering sampling molecular dynamics were combined with ab initio calculations to overcome the sampling challenge in cluster structure search and refinement. The calculated vertical detachment energies of CaCl2(H2O)n- (n = 0-8) were compared with the values obtained from photoelectron spectra, and consistency was found between the experiment and computation. Separation of the Cl-Ca ion pair is investigated in CaCl2(H2O)n- anions, where the first Ca-Cl ionic bond required 4 water molecules, and both Ca-Cl bonds are broken when the number of water molecules is larger than 7. For neutral CaCl2(H2O)n clusters, breaking of the first Ca-Cl bond starts at n = 5, and 8 water molecules are not enough to separate the two ion pairs. Comparing with the observations on magnesium chloride, it shows that separating one ion pair in CaCl2(H2O)n requires fewer water molecules than those for MgCl2(H2O)n. Coincidentally, the solubility of calcium chloride is higher than that of magnesium chloride in bulk solutions.

Accurate placement of substrate RNA by Gar1 in H/ACA RNA-guided pseudouridylation.

H/ACA RNA-guided ribonucleoprotein particle (RNP), the most complicated RNA pseudouridylase so far known, uses H/ACA guide RNA for substrate capture and four proteins (Cbf5, Nop10, L7Ae and Gar1) for pseudouridylation. Although it was shown that Gar1 not only facilitates the product release, but also enhances the catalytic activity, the chemical role that Gar1 plays in this complicated machinery is largely unknown. Kinetics measurement on Pyrococcus furiosus RNPs at different temperatures making use of fluorescence anisotropy showed that Gar1 reduces the catalytic barrier through affecting the activation entropy instead of enthalpy. Site-directed mutagenesis combined with molecular dynamics simulations demonstrated that V149 in the thumb loop of Cbf5 is critical in placing the target uridine to the right position toward catalytic D85 of Cbf5. The enzyme elegantly aligns the position of uridine in the catalytic site with the help of Gar1. In addition, conversion of uridine to pseudouridine results in a rigid syn configuration of the target nucleotide in the active site and causes Gar1 to pull out the thumb. Both factors guarantee the efficient release of the product.

Dynamics of spontaneous flipping of a mismatched base in DNA duplex.

DNA base flipping is a fundamental theme in DNA biophysics. The dynamics for a B-DNA base to spontaneously flip out of the double helix has significant implications in various DNA-protein interactions but are still poorly understood. The spontaneous base-flipping rate obtained previously via the imino proton exchange assay is most likely the rate of base wobbling instead of flipping. Using the diffusion-decelerated fluorescence correlation spectroscopy together with molecular dynamics simulations, we show that a base of a single mismatched base pair (T-G, T-T, or T-C) in a double-stranded DNA can spontaneously flip out of the DNA duplex. The extrahelical lifetimes are on the order of 10 ms, whereas the intrahelical lifetimes range from 0.3 to 20 s depending on the stability of the base pairs. These findings provide detailed understanding on the dynamics of DNA base flipping and lay down foundation to fully understand how exactly the repair proteins search and locate the target mismatched base among a vast excess of matched DNA bases.

Conformation switching of AIM2 PYD domain revealed by NMR relaxation and MD simulation.

Protein absent in melanoma 2 (AIM2) is a double-strand DNA (ds DNA) sensor mainly located in cytoplasm of cell. It includes one N terminal PYD domain and one C terminal HIN domain. When the ds DNA such as DNA viruses and bacteria entered cytoplasm, the HIN domain of AIM2 will recognize and bind to DNA, and the PYD domain will bind to ASC protein which will result in the formation of AIM2 inflammasome. Three AIM2 PYD domain structures have been solved, but every structure yields a unique conformation around the α3 helix region. To understand why different AIM2 PYD structures show different conformations in this region, we use NMR relaxation techniques to study the backbone dynamics of mouse AIM2 PYD domain and perform molecular dynamics (MD) simulations on both mouse and human AIM2 PYD structures. Our results indicate that this region is highly flexible in both mouse and human AIM2 PYD domains, and the PYD domain may exist as a conformation ensemble in solution. Different environment makes the population vary among pre-existing conformational substrates of the ensemble, which may be the reason why different AIM2 PYD structures were observed under different conditions. Further docking analysis reveals that the conformation switching may be important for the autoinhibition of the AIM2 protein.

From thermodynamics to kinetics: enhanced sampling of rare events.

Despite great advances in molecular dynamics simulations, there remain large gaps between the simulations and experimental observations in terms of the time and length scales that can be approached. Developing fast and accurate algorithms and methods is of ultimate importance to bridge these gaps. In this Account, we briefly summarize recent efforts in such directions. In particular, we focus on integrated tempering sampling. The efficiency of this sampling method has been demonstrated by applications to a range of chemical and biological problems: protein folding, molecular cluster structure searches, and chemical reactions. The combination of integrated tempering sampling and a trajectory sampling method allows the calculation of rate constants and reaction pathways without predefined collective coordinates.

Efficient sampling over rough energy landscapes with high barriers: A combination of metadynamics with integrated tempering sampling.

In order to efficiently overcome high free energy barriers embedded in a complex energy landscape and calculate overall thermodynamics properties using molecular dynamics simulations, we developed and implemented a sampling strategy by combining the metadynamics with (selective) integrated tempering sampling (ITS/SITS) method. The dominant local minima on the potential energy surface (PES) are partially exalted by accumulating history-dependent potentials as in metadynamics, and the sampling over the entire PES is further enhanced by ITS/SITS. With this hybrid method, the simulated system can be rapidly driven across the dominant barrier along selected collective coordinates. Then, ITS/SITS ensures a fast convergence of the sampling over the entire PES and an efficient calculation of the overall thermodynamic properties of the simulation system. To test the accuracy and efficiency of this method, we first benchmarked this method in the calculation of ϕ - ψ distribution of alanine dipeptide in explicit solvent. We further applied it to examine the design of template molecules for aromatic meta-C-H activation in solutions and investigate solution conformations of the nonapeptide Bradykinin involving slow cis-trans isomerizations of three proline residues.

Combine umbrella sampling with integrated tempering method for efficient and accurate calculation of free energy changes of complex energy surface.

Umbrella sampling is an efficient method for the calculation of free energy changes of a system along well-defined reaction coordinates. However, when there exist multiple parallel channels along the reaction coordinate or hidden barriers in directions perpendicular to the reaction coordinate, it is difficult for conventional umbrella sampling to reach convergent sampling within limited simulation time. Here, we propose an approach to combine umbrella sampling with the integrated tempering sampling method. The umbrella sampling method is applied to chemically more relevant degrees of freedom that possess significant barriers. The integrated tempering sampling method is used to facilitate the sampling of other degrees of freedom which may possess statistically non-negligible barriers. The combined method is applied to two model systems, butane and ACE-NME molecules, and shows significantly improved sampling efficiencies as compared to standalone conventional umbrella sampling or integrated tempering sampling approaches. Further analyses suggest that the enhanced performance of the new method come from the complemented advantages of umbrella sampling with a well-defined reaction coordinate and integrated tempering sampling in orthogonal space. Therefore, the combined approach could be useful in the simulation of biomolecular processes, which often involves sampling of complex rugged energy landscapes.

Discriminating trpzip2 and trpzip4 peptides' folding landscape using the two-dimensional infrared spectroscopy: a simulation study.

We analyzed, based on the theoretical spectroscopic modeling, how the differences in the folding landscapes of two ß-hairpin peptides trpzip2 and trpzip4 are reflected in their thermal unfolding infrared measurements. The isotope-edited equilibrium FTIR and two dimensional infrared spectra of the two peptides were calculated, using the nonlinear exciton propagation method, at a series of temperatures. The spectra calculations were based on the configuration distributions generated using the GB(OBC) implicit solvent MD simulation and the integrated tempering sampling technique. Conformational analysis revealed the different local thermal stabilities for these two peptides, which suggested the different folding landscapes. Our study further suggested that the ellipticities of the isotope peaks in the coherent IR signals are more sensitive to these local stability differences compared with other spectral features such as the peak intensities. Our technique can thus be combined with the relevant experimental measurements to achieve a better understanding of the peptide folding behaviors.

PMC-IZ: A Simple Algorithm for the Electrostatics Calculation in Slab Geometric Molecular Dynamics Simulations.

Slab geometric systems are widely utilized in molecular simulations. However, an efficient, straightforward, and accurate method for calculating electrostatic interactions in these systems for molecular dynamics (MD) simulations is still needed. This review introduces a PME-like approach called PMC-IZ, specifically designed for slab geometric systems. Traditional approaches for long-range electrostatic interaction calculations in slab geometry typically involve Ewald summation, where the Gaussian charge density is summed within 3D unit cells and then integrated in the 2D periodic space. In the proposed approach here, the Poisson equation was solved for a single Gaussian charge density within 2Dl periodic space, followed by convolution within 3D unit cells using an effective potential as the convolution kernel for summation. The effective potential ensures that the solution within the region of interest adheres strictly to 2D periodic boundary conditions while inherently possessing 3D periodic boundary condition properties. The PMC-IZ method provides for such systems accurate treatment of electrostatic interactions, overcomes limitations associated with finite vacuum layers, and offers improved computational efficiency. We thus postulate that this method provides a valuable tool for studying electrostatic interactions in slab geometric system MD simulations. It has promising applications in various areas such as surface science, catalysis, and materials research, where accurate modeling of slab geometric electrostatic interactions is essential.

Single-molecule detection reveals knot sliding in TrmD denaturation.

An increasing number of proteins are found to contain a knot in their polypeptide chain. Although some studies have looked into the folding mechanism of knotted proteins, why and how these complex topologies form are still far from being fully answered. Moreover, no experimental information about how the knot moves during the protein-folding process is available. Herein, by combining single-molecule fluorescence resonance energy transfer (smFRET) experiments with molecular dynamics (MD) simulations, we performed a detailed study to characterize the knot in the denatured state of TrmD, a knotted tRNA (guanosine-1) methyltransferase from Escherichia coli, as a model system. We found that the knot still existed in the unfolded state of TrmD, consistent with the results for two other knotted proteins, YibK and YbeA. More interestingly, both smFRET experiments and MD simulations revealed that the knot slid towards the C-terminal during the unfolding process, which could be explained by the relatively strong interactions between the ß-sheet core at the N terminal of the native knot region. The size of the knot in the unfolded state is not larger than that in the native state. In addition, the knot slid in a "downhill" mode with simultaneous chain collapse in the denatured state.