Leveraging the Information from Markov State Models To Improve the Convergence of Umbrella Sampling Simulations

Kemeny Constant-Based Optimization of Network Clustering Using Graph Neural Networks

The recent trend in using network and graph structures to represent a variety of different data types has renewed interest in the graph partitioning (GP) problem. This interest stems from the need for general methods that can both efficiently identify network communities and reduce the dimensionality of large graphs while satisfying various application-specific criteria. Traditional clustering algorithms often struggle to capture the complex relationships within graphs and generalize to arbitrary clustering criteria. The emergence of graph neural networks (GNNs) as a powerful framework for learning representations of graph data provides new approaches to solving the problem. Previous work has shown GNNs to be capable of proposing partitionings using a variety of criteria. However, these approaches have not yet been extended to Markov chains or kinetic networks. These arise frequently in the study of molecular systems and are of particular interest to the biomolecular modeling community. In this work, we propose several GNN-based architectures to tackle the GP problem for Markov Chains described as kinetic networks. This approach aims to maximize the Kemeny constant, which is a variational quantity and it represents the sum of time scales of the system. We propose using an encoder-decoder architecture and show how simple GraphSAGE-based GNNs with linear layers can outperform much larger and more expressive attention-based models in this context. As a proof of concept, we first demonstrate the method's ability to cluster randomly connected graphs. We also use a linear chain architecture corresponding to a 1D free energy profile as our kinetic network. Subsequently, we demonstrate the effectiveness of our method through experiments on a data set derived from molecular dynamics. We compare the performance of our method to other partitioning techniques, such as PCCA+. We explore the importance of feature and hyperparameter selection and propose a general strategy for large-scale parallel training of GNNs for discovering optimal graph partitionings.

Lithium Pre-Storage Enables High Initial Coulombic Efficiency and Stable Lithium-Enriched Silicon/Graphite Anode

Application of silicon-based anodes is significantly challenged by low initial Coulombic efficiency (ICE) and poor cyclability. Traditional pre-lithiation reagents often pose safety concerns due to their unstable chemical nature. Achieving a balance between water-stability and high ICE in prelithiated silicon is a critical issue. Here, we present a lithium-enriched silicon/graphite material with an ultra-high ICE of ≥110 % through a high-stable lithium pre-storage methodology. Lithium pre-storage prepared a nano-drilled graphite material with surficial lithium functional groups, which can form chemical bonds with adjacent silicon during high-temperature sintering. This results in an unexpected O-Li-Si interaction, leading to in situ pre-lithiation of silicon nanoparticles and providing high stability in air and water. Additionally, the lithium-enriched silicon/graphite materials impart a combination of high ICE, high specific capacity (620 mAh g-1), and long cycling stability (>400 cycles). This study opens up a promising avenue for highly air- and water-stable silicon anode prelithiation methods.

OGRe: Optimal Grid Refinement Protocol for Accurate Free Energy Surfaces and Its Application in Proton Hopping in Zeolites and 2D COF Stacking

While free energy surfaces are the crux of our understanding of many chemical and biological processes, their accuracy is generally unknown. Moreover, many developments to improve their accuracy are often complicated, limiting their general use. Luckily, several tools and guidelines are already in place to identify these shortcomings, but they are typically lacking in flexibility or fail to systematically determine how to improve the accuracy of the free energy calculation. To overcome these limitations, this work introduces OGRe, a Python package for optimal grid refinement in an arbitrary number of dimensions. OGRe is based on three metrics that gauge the confinement, consistency, and overlap of each simulation in a series of umbrella sampling (US) simulations, an enhanced sampling technique ubiquitously adopted to construct free energy surfaces for hindered processes. As these three metrics are fundamentally linked to the accuracy of the weighted histogram analysis method adopted to generate free energy surfaces from US simulations, they facilitate the systematic construction of accurate free energy profiles, where each metric is driven by a specific umbrella parameter. This allows for the derivation of a consistent and optimal collection of umbrellas for each simulation, largely independent of the initial values, thereby dramatically increasing the ease-of-use toward accurate free energy surfaces. As such, OGRe is particularly suited to determine complex free energy surfaces with large activation barriers and shallow minima, which underpin many physical and chemical transformations and hence to further our fundamental understanding of these processes.

Toward a structural identification of metastable molecular conformations

Interpreting high-dimensional data from molecular dynamics simulations is a persistent challenge. In this paper, we show that for a small peptide, deca-alanine, metastable states can be identified through a neural net based on structural information alone. While processing molecular dynamics data, dimensionality reduction is a necessary step that projects high-dimensional data onto a low-dimensional representation that, ideally, captures the conformational changes in the underlying data. Conventional methods make use of the temporal information contained in trajectories generated through integrating the equations of motion, which forgoes more efficient sampling schemes. We demonstrate that EncoderMap, an autoencoder architecture with an additional distance metric, can find a suitable low-dimensional representation to identify long-lived molecular conformations using exclusively structural information. For deca-alanine, which exhibits several helix-forming pathways, we show that this approach allows us to combine simulations with different biasing forces and yields representations comparable in quality to other established methods. Our results contribute to computational strategies for the rapid automatic exploration of the configuration space of peptides and proteins.

Potential Mean Force from Umbrella Sampling Simulations: What Can We Learn and What Is Missed?

Changes in free energy provide valuable information for molecular recognition, including both ligand-receptor binding thermodynamics and kinetics. Umbrella sampling (US), a widely used free energy calculation method, has long been used to explore the dissociation process of ligand-receptor systems and compute binding free energy. In existing publications, the binding free energy computed from the potential of mean force (PMF) with US simulation mostly yielded "ball park" values with experimental data. However, the computed PMF values are highly influenced by factors such as initial conformations and/or trajectories provided, the reaction coordinate, and sampling of conformational space in each US window. These critical factors have rarely been carefully studied. Here we used US to study the guest aspirin and 1-butanol dissociation processes of β-cyclodextrin (β-CD) and an inhibitor SB2 dissociation from a p38α mitogen-activated protein kinase (MAPK) complex. For β-CD, we used three different β-CD conformations to generate the dissociation path with US windows. For p38α, we generated the dissociation pathway by using accelerated molecular dynamics followed by conformational relaxing with short conventional MD, steered MD, and manual pulling. We found that, even for small β-CD complexes, different β-CD conformations altered the height of the PMF, but the pattern of PMF was not affected if the MD sampling in each US window was well-converged. Because changing the macrocyclic ring conformation needs to rotate dihedral angles in the ring, a bound ligand largely restrains the motion of cyclodextrin. Therefore, once a guest is in the binding site, cyclodextrin cannot freely change its initial conformation, resulting in different absolute heights of the PMF, which cannot be overcome by running excessively long MD simulations for each US window. Moreover, if the US simulations were not converged, the important barrier and minimum were missed. For ligand-protein systems, our studies also suggest that the dissociation trajectories modeled by an enhanced sampling method must maintain a natural molecular movement to avoid biased PMF plots when using US simulations.

Orientational and Folding Thermodynamics via Electric Dipole Moment Restraining

The projection of molecular processes onto a small set of relevant descriptors, the so-called reaction coordinates or collective variables (CVs), is a technique nowadays routinely employed by the biomolecular simulation community. In this work, we implemented two CVs to manipulate the orientation (i.e., angle) (μ⃗_a) and magnitude (|μ⃗|) of the electric dipole moment. In doing so, we studied the thermodynamics of water orientation under the application of external voltages and the folding of two polypeptides at zero-field conditions. The projection of the free-energy [potential of mean force (PMF)] along water orientation defined an upper limit of around 0.3 V for irrelevant thermodynamic effects. On the other hand, sufficiently strong μ⃗_a restraints applied on 12-alanine (Ala₁₂) triggered structural effects because of the alignment of local dipoles; for lower restraints, a full-body rotation is achieved. The manipulation of |μ⃗| produced strong perturbations on the secondary structure of Ala₁₂, promoting an enhanced sampling to its configurational space. Rigorous free-energy calculations in the form of 2-D PMFs for deca-alanine showed the utility of |μ⃗| as a reaction coordinate to study folding in small α helices. As a whole, we propose that the manipulation of both components of the dipole moment, μ⃗_a and |μ⃗|, provides thermodynamics insights into the structural conformation and stability of biomolecules. These new CVs are implemented in the Colvars module, available for NAMD and LAMMPS.

Adaptive ensemble simulations of biomolecules

Recent advances in both theory and computational power have created opportunities to simulate biomolecular processes more efficiently using adaptive ensemble simulations. Ensemble simulations are now widely used to compute a number of individual simulation trajectories and analyze statistics across them. Adaptive ensemble simulations offer a further level of sophistication and flexibility by enabling high-level algorithms to control simulations-based on intermediate results. We review some of the adaptive ensemble algorithms and software infrastructure currently in use and outline where the complexities of implementing adaptive simulation have limited algorithmic innovation to date. We describe an adaptive ensemble API to overcome some of these barriers and more flexibly and simply express adaptive simulation algorithms to help realize the power of this type of simulation.

Folding a viral peptide in different membrane environments: pathway and sampling analyses

Flock House virus (FHV) is a well-characterized model system to study infection mechanisms in non-enveloped viruses. A key stage of the infection cycle is the disruption of the endosomal membrane by a component of the FHV capsid, the membrane active γ peptide. In this study, we perform all-atom molecular dynamics simulations of the 21 N-terminal residues of the γ peptide interacting with membranes of differing compositions. We carry out umbrella sampling calculations to study the folding of the peptide to a helical state in homogenous and heterogeneous membranes consisting of neutral and anionic lipids. From the trajectory data, we evaluate folding energetics and dissect the mechanism of folding in the different membrane environments. We conclude the study by analyzing the extent of configurational sampling by performing time-lagged independent component analysis.

Systematic Dissociation Pathway Searches Guided by Principal Component Modes

We introduce a novel method, Pathway Search guided by Internal Motions (PSIM), that efficiently finds molecular dissociation pathways of a ligand-receptor system with guidance from principal component (PC) modes obtained from molecular dynamics (MD) simulations. Modeling ligand-receptor dissociation pathways can provide insights into molecular recognition and has practical applications, including understanding kinetic mechanisms and barriers to binding/unbinding as well as design of drugs with desired kinetic properties. PSIM uses PC modes in multilayer internal coordinates to identify natural molecular motions that guide the search for conformational switches and unbinding pathways. The new multilayer internal coordinates overcome problems with Cartesian and classical internal coordinates that fail to smoothly present dihedral rotation or generate nonphysical distortions. We used HIV-1 protease, which has large-scale flap motions, as an example protein to demonstrate use of the multilayer internal coordinates. We provide examples of algorithms and implementation of PSIM with alanine dipeptide and chemical host-guest systems, 2-naphthyl ethanol-β-cyclodextrin and tetramethylammonium-cryptophane complexes. Tetramethylammonium-cryptophane has slow binding/unbinding kinetics. Its residence time, the length to dissociate tetramethylammonium from the host, is ∼14 s from experiments, and PSIM revealed 4 dissociation pathways in approximately 150 CPU h. We also searched the releasing pathways for the product glyceraldehyde-3-phosphate from tryptophan synthase, and one complete dissociation pathway was constructed after running multiple search iterations in approximately 300 CPU h. With guidance by internal PC modes from MD simulations, the PSIM method has advantages over simulation-based methods to search for dissociation pathways of molecular systems with slow noncovalent kinetic behavior.