Lifting the Curse of Dimensionality on Enhanced Sampling of Reaction Networks with Parallel Bias Metadynamics

Adaptive accelerated reactive molecular dynamics driven by parallel collective variables overcoming dimensionality explosion

ReaxFF reactive molecular dynamics has significantly advanced the exploration of chemical reaction mechanisms in complex systems. However, it faces several challenges: (1) the prevalent use of excessively high temperatures (>2000 K), (2) a time scale considerably shorter than the experimental timeframes (nanoseconds vs seconds), and (3) the constraining impact of dimensionality growth due to collective variables on the expansiveness of research systems. To overcome these issues, we introduced Parallel Collective Variable-Driven Adaptive Accelerated Reaction Molecular Dynamics (PCVR), which integrates metadynamics with ReaxFF. This method incorporates bond distortion based on each bond type for customized Collective Variable (CV) parameterization, facilitating independent parallel acceleration. Simultaneously, the sampling was confined to fixed cutoff ranges for distinct bond distortions, effectively overcoming the challenge of the CV dimensionality explosion. This extension enhances the applicability of ReaxFF to non-strongly coupled systems with numerous reaction energy barriers and mitigates the system size limitations. Using accelerated reactive molecular dynamics, the oxidation of ester-based oil was simulated with 31 808 atoms at 500 K for 64 s. This achieved 61% efficiency compared to the original ReaxFF and was ∼37 times faster than previous methods. Unlike ReaxFF's high-temperature constraints, PCVR accurately reveals the pivotal role of oxygen in ester oxidation at industrial temperatures, producing polymers consistent with the sludge formation observed in ester degradation experiments. This method promises to advance reactive molecular dynamics by enabling simulations at lower temperatures, extending to second-level timescales, and accommodating systems with millions of atoms.

Thermodynamic Basis for the Stabilization of Helical Peptoids by Chiral Sidechains

Peptoids are a class of highly customizable biomimetic foldamers that retain properties from both proteins and polymers. It has been shown that peptoids can adopt peptide-like secondary structures through the careful selection of sidechain chemistries, but the underlying conformational landscapes that drive these assemblies at the molecular level remain poorly understood. Given the high flexibility of the peptoid backbone, it is essential that methods applied to study peptoid secondary structure formation possess the requisite sensitivity to discriminate between structurally similar yet energetically distinct microstates. In this work, a generalizable simulation scheme is used to robustly sample the complex folding landscape of various 12mer polypeptoids, resulting in a predictive model that links sidechain chemistry with preferential assembly into one of 12 accessible backbone motifs. Using a variant of the metadynamics sampling method, four peptoid dodecamers are simulated in water: sarcosine, N-(1-phenylmethyl)glycine (Npm), (S)-N-(1-phenylethyl)glycine (Nspe), and (R)-N-(1-phenylethyl)glycine (Nrpe)─to determine the underlying entropic and energetic impacts of hydrophobic and chiral peptoid sidechains on secondary structure formation. Our results indicate that the driving forces to assemble Nrpe and Nspe sequences into polyproline type-I helices in water are found to be enthalpically driven, with small benefits from an entropic gain for isomerization and steric strain due to the presence of the chiral center. The minor entropic gains from bulky chiral sidechains in Nrpe- and Nspe-containing peptoids can be explained through increased configurational entropy in the cis state. However, overall assembly into a helix is found to be overall entropically unfavorable. These results highlight the importance of considering the many various competing interactions in the rational design of peptoid secondary structure building blocks.

Origins of Conformational Heterogeneity in Peptoid Helices Formed by Chiral N-1-Phenylethyl Sidechains

N-substituted glycines (polypeptoids) containing chiral hydrophobic sidechains are known to fold into biomimetic alpha helices. These helix formers often produce conformationally heterogeneous structures and are difficult to characterize at a sub-nanometer resolution. Previously, peptoid N-1-phenylethyl (R)-enantiomer sidechains (Nspe) were inferred from various experiments to form right-handed helices and (S)-enantiomers (Nrpe), left-handed helices. Prior computational work for N(s/r)pe oligomers has struggled to reproduce this trend. Herein, quantum mechanics calculations and molecular dynamics simulations are used to understand the origins of this discrepancy. Results from DFT and molecular mechanics calculations on a variety of Nspe and Nrpe oligomers as a function of chain length are in agreement, showing that Nspe and Nrpe prefer left- and right-handed helices, respectively. Additional metadynamics simulations are used to study Nrpe and Nspe oligomers folding in water. These results show that the free-energy driving forces for assembly into a helical backbone configuration are very small (within ∼k_BT). Lastly, we compare DFT calculations for other experimentally characterized peptoid sidechains, N(r/s)sb, N(r/s)tbe, and N(r/s)npe. In this analysis, we show that peptoid sidechains determined to be more robust experimentally (tbe and npe) have helical preferences opposite the trend seen in less robust assemblies formed by N(r/s)pe and N(r/s)sb chemistries. The more robust tbe and nnpe favor the (S)-enantiomer to right-handed and the (R)-enantiomers to left-handed helices.

Hierarchical Materials from High Information Content Macromolecular Building Blocks: Construction, Dynamic Interventions, and Prediction

Hierarchical materials that exhibit order over multiple length scales are ubiquitous in nature. Because hierarchy gives rise to unique properties and functions, many have sought inspiration from nature when designing and fabricating hierarchical matter. More and more, however, nature's own high-information content building blocks, proteins, peptides, and peptidomimetics, are being coopted to build hierarchy because the information that determines structure, function, and interfacial interactions can be readily encoded in these versatile macromolecules. Here, we take stock of recent progress in the rational design and characterization of hierarchical materials produced from high-information content blocks with a focus on stimuli-responsive and "smart" architectures. We also review advances in the use of computational simulations and data-driven predictions to shed light on how the side chain chemistry and conformational flexibility of macromolecular blocks drive the emergence of order and the acquisition of hierarchy and also on how ionic, solvent, and surface effects influence the outcomes of assembly. Continued progress in the above areas will ultimately usher in an era where an understanding of designed interactions, surface effects, and solution conditions can be harnessed to achieve predictive materials synthesis across scale and drive emergent phenomena in the self-assembly and reconfiguration of high-information content building blocks.

Transient exposure of a buried phosphorylation site in an autoinhibited protein

Autoinhibition is a mechanism used to regulate protein function, often by making functional sites inaccessible through the interaction with a cis-acting inhibitory domain. Such autoinhibitory domains often display a substantial degree of structural disorder when unbound, and only become structured in the inhibited state. These conformational dynamics make it difficult to study the structural origin of regulation, including effects of regulatory post-translational modifications. Here, we study the autoinhibition of the Dbl Homology domain in the protein Vav1 by the so-called acidic inhibitory domain. We use molecular simulations to study the process by which a mostly unstructured inhibitory domain folds upon binding and how transient exposure of a key buried tyrosine residue makes it accessible for phosphorylation. We show that the inhibitory domain, which forms a helix in the bound and inhibited stated, samples helical structures already before binding and that binding occurs via a molten-globule-like intermediate state. Together, our results shed light on key interactions that enable the inhibitory domain to sample a finely tuned equilibrium between an inhibited and a kinase-accessible state.

Reweighted Jarzynski Sampling: Acceleration of Rare Events and Free Energy Calculation with a Bias Potential Learned from Nonequilibrium Work

We introduce a simple enhanced sampling approach for the calculation of free energy differences and barriers along a one-dimensional reaction coordinate. First, a small number of short nonequilibrium simulations are carried out along the reaction coordinate, and the Jarzynski equality is used to learn an approximate free energy surface from the nonequilibrium work distribution. This free energy estimate is represented in a compact form as an artificial neural network and used as an external bias potential to accelerate rare events in a subsequent molecular dynamics simulation. The final free energy estimate is then obtained by reweighting the equilibrium probability distribution of the reaction coordinate sampled under the influence of the external bias. We apply our reweighted Jarzynski sampling recipe to four processes of varying scales and complexities─spanning chemical reaction in the gas phase, pair association in solution, and droplet nucleation in supersaturated vapor. In all cases, we find reweighted Jarzynski sampling to be a very efficient strategy, resulting in rapid convergence of the free energy to high precision.

Closer Look at Inverse Electron Demand Diels–Alder and Nucleophilic Addition Reactions on s-Tetrazines Using Enhanced Sampling Methods

Inverse electron demand [4+2] Diels–Alder (iEDDA) reactions as well as unprecedented nucleophilic (azaphilic) additions of R-substituted silyl-enol ethers (where R is Phenyl, Methyl, and Hydrogen) to 1,2,4,5-tetrazine (s-tetrazine) catalyzed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3 have recently been discovered (Simon et al. in Org Lett 23(7):2426–2430, 2021), where static calculations were employed for calculation of activation energies. In order to have a more realistic dynamic description of these reactions in explicit solution at ambient conditions, in this work we use a semiempirical tight-binding method combined with enhanced sampling techniques to calculate free energy surfaces (FESs) of the iEDDA and azaphilic addition reactions. Relevant products of not only s-tetrazine but also its derivatives such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3-mediated s-tetrazine adducts are investigated. We reconstruct the FESs of the iEDDA and azaphilic addition reactions using metadynamics and blue moon ensemble, and compare the ability of different collective variables (CVs) including bond distances, Social PeRmutation INvarianT (SPRINT) coordinates, and path-CV to describe the reaction pathway. We find that when a bulky Phenyl is used as a substituent at the dienophile the azaphilic addition is preferred over the iEDDA reaction. In addition, we also investigate the effect of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3 in the diene and steric hindrance in the dienophile on the competition between the iEDDA and azaphilic addition reactions, providing chemical insight for reaction design.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

Solvent oligomerization pathways facilitated by electrolyte additives during solid-electrolyte interphase formation

The solid-electrolyte interphase (SEI) layer formation is known to play an important role in determining the lifetime of lithium-ion batteries. A thin, stable SEI layer is linked to overall improved battery performance and longevity, however, the factors and mechanisms that lead to optimal SEI morphology and composition are not well understood. Inclusion of electrolyte additives (fluoroethylene carbonate, FEC; and vinylene carbonate, VC) is often necessary for improving SEI characteristics. To understand how these electrolyte additives impact SEI formation, molecular dynamics (MD) and density functional theory (DFT) simulations were employed to study the reaction networks and oligomerization pathways, respectively, for three systems containing ethylene carbonate (EC), a lithium ion, and FEC or VC. MD simulations suggest radical oligomerization pathways analogous to traditional oligomerization with nucleophilic alkoxide species via SN1 reaction mechanisms. Both SN1 and SN2 mechanisms were studied for all three systems using DFT. Oligomerization reactions were studied with both a standard alkoxide species and a ring-opened EC radical as the nucleophiles and EC, FEC, and VC as the electrophiles. For all cases, FEC and VC exhibited lower free energy barriers and more stable adducts when compared with EC. We conclude that one of the role of additives is to modify the oligomerization process of EC by introducing branching points (FEC) or termination points (VC).