Weighted Ensemble Simulation: Review of Methodology, Applications, and Software

AI in cellular engineering and reprogramming

During the last decade, artificial intelligence (AI) has increasingly been applied in biophysics and related fields, including cellular engineering and reprogramming, offering novel approaches to understand, manipulate, and control cellular function. The potential of AI lies in its ability to analyze complex datasets and generate predictive models. AI algorithms can process large amounts of data from single-cell genomics and multiomic technologies, allowing researchers to gain mechanistic insights into the control of cell identity and function. By integrating and interpreting these complex datasets, AI can help identify key molecular events and regulatory pathways involved in cellular reprogramming. This knowledge can inform the design of precision engineering strategies, such as the development of new transcription factor and signaling molecule cocktails, to manipulate cell identity and drive authentic cell fate across lineage boundaries. Furthermore, when used in combination with computational methods, AI can accelerate and improve the analysis and understanding of the intricate relationships between genes, proteins, and cellular processes. In this review article, we explore the current state of AI applications in biophysics with a specific focus on cellular engineering and reprogramming. Then, we showcase a couple of recent applications where we combined machine learning with experimental and computational techniques. Finally, we briefly discuss the challenges and prospects of AI in cellular engineering and reprogramming, emphasizing the potential of these technologies to revolutionize our ability to engineer cells for a variety of applications, from disease modeling and drug discovery to regenerative medicine and biomanufacturing.

BAD-NEUS: Rapidly converging trajectory stratification

An issue for molecular dynamics simulations is that events of interest often involve timescales that are much longer than the simulation time step, which is set by the fastest timescales of the model. Because of this timescale separation, direct simulation of many events is prohibitively computationally costly. This issue can be overcome by aggregating information from many relatively short simulations that sample segments of trajectories involving events of interest. This is the strategy of Markov state models (MSMs) and related approaches, but such methods suffer from approximation error because the variables defining the states generally do not capture the dynamics fully. By contrast, once converged, the weighted ensemble (WE) method aggregates information from trajectory segments so as to yield unbiased estimates of both thermodynamic and kinetic statistics. Unfortunately, errors decay no faster than unbiased simulation in WE as originally formulated and commonly deployed. Here, we introduce a theoretical framework for describing WE that shows that the introduction of an approximate stationary distribution on top of the stratification, as in nonequilibrium umbrella sampling (NEUS), accelerates convergence. Building on ideas from MSMs and related methods, we generalize the NEUS approach in such a way that the approximation error can be reduced systematically. We show that the improved algorithm can decrease the simulation time required to achieve the desired precision by orders of magnitude.

Simulation-driven design of stabilized SARS-CoV-2 spike S2 immunogens

The full-length prefusion-stabilized SARS-CoV-2 spike (S) is the principal antigen of COVID-19 vaccines. Vaccine efficacy has been impacted by emerging variants of concern that accumulate most of the sequence modifications in the immunodominant S1 subunit. S2, in contrast, is the most evolutionarily conserved region of the spike and can elicit broadly neutralizing and protective antibodies. Yet, S2's usage as an alternative vaccine strategy is hampered by its general instability. Here, we use a simulation-driven approach to design S2-only immunogens stabilized in a closed prefusion conformation. Molecular simulations provide a mechanistic characterization of the S2 trimer's opening, informing the design of tryptophan substitutions that impart kinetic and thermodynamic stabilization. Structural characterization via cryo-EM shows the molecular basis of S2 stabilization in the closed prefusion conformation. Informed by molecular simulations and corroborated by experiments, we report an engineered S2 immunogen that exhibits increased protein expression, superior thermostability, and preserved immunogenicity against sarbecoviruses.

WEDAP: A Python Package for Streamlined Plotting of Molecular Simulation Data

Given the growing interest in path sampling methods for extending the time scales of molecular dynamics (MD) simulations, there has been great interest in software tools that streamline the generation of plots for monitoring the progress of large-scale simulations. Here, we present the WEDAP Python package for simplifying the analysis of data generated from either conventional MD simulations or the weighted ensemble (WE) path sampling method, as implemented in the widely used WESTPA software package. WEDAP facilitates (i) the parsing of WE simulation data stored in highly compressed, hierarchical HDF5 files and (ii) incorporates trajectory weights from WE simulations into all generated plots. Our Python package consists of multiple user-friendly interfaces: a command-line interface, a graphical user interface, and a Python application programming interface. We demonstrate the plotting features of WEDAP through a series of examples using data from WE and conventional MD simulations that focus on the HIV-1 capsid protein's C-terminal domain dimer as a showcase system. The source code for WEDAP is freely available on GitHub at https://github.com/chonglab-pitt/wedap.

QM/CG-MM: Systematic Embedding of Quantum Mechanical Systems in a Coarse-Grained Environment with Accurate Electrostatics

Quantum Mechanics/Molecular Mechanics (QM/MM) can describe chemical reactions in molecular dynamics (MD) simulations at a much lower cost than ab initio MD. Still, it is prohibitively expensive for many systems of interest because such systems usually require long simulations for sufficient statistical sampling. Additional MM degrees of freedom are often slow and numerous but secondary in interest. Coarse-graining (CG) is well-known to be able to speed up sampling through both reduction in simulation cost and the ability to accelerate the dynamics. Therefore, embedding a QM system in a CG environment can be a promising way of expediting sampling without compromising the information about the QM subsystem. Sinitskiy and Voth first proposed the theory of Quantum Mechanics/Coarse-grained Molecular Mechanics (QM/CG-MM) with a bottom-up CG mapping. Mironenko and Voth subsequently introduced the DFT-QM/CG-MM formalism to couple a Density Functional Theory (DFT) treated QM system and to an apolar environment. Here, we present a more complete theory that addresses MM environments with significant polarity by explicitly accounting for the electrostatic coupling. We demonstrate our QM/CG-MM method with a chloride-methyl chloride SN2 reaction system in acetone, which is sensitive to solvent polarity. The method accurately recapitulates the potential of mean force for the substitution reaction, and the reaction barrier from the best model agrees with the atomistic simulations within sampling error. These models also have generalizability. In two other reactive systems that they have not been trained on, the QM/CG-MM model still achieves the same level of agreement with the atomistic QM/MM models. Finally, we show that in these examples the speed-up in the sampling is proportional to the acceleration of the rotational dynamics of the solvent in the CG system.

Identification of small-molecule ligand-binding sites on and in the ARNT PAS-B domain

Transcription factors are challenging to target with small-molecule inhibitors due to their structural plasticity and lack of catalytic sites. Notable exceptions include naturally ligand-regulated transcription factors, including our prior work with the hypoxia-inducible factor (HIF)-2 transcription factor, showing that small-molecule binding within an internal pocket of the HIF-2α Per-Aryl hydrocarbon Receptor Nuclear Translocator (ARNT)-Sim (PAS)-B domain can disrupt its interactions with its dimerization partner, ARNT. Here, we explore the feasibility of targeting small molecules to the analogous ARNT PAS-B domain itself, potentially opening a promising route to modulate several ARNT-mediated signaling pathways. Using solution NMR fragment screening, we previously identified several compounds that bind ARNT PAS-B and, in certain cases, antagonize ARNT association with the transforming acidic coiled-coil containing protein 3 transcriptional coactivator. However, these ligands have only modest binding affinities, complicating characterization of their binding sites. We address this challenge by combining NMR, molecular dynamics simulations, and ensemble docking to identify ligand-binding "hotspots" on and within the ARNT PAS-B domain. Our data indicate that the two ARNT/transforming acidic coiled-coil containing protein 3 inhibitors, KG-548 and KG-655, bind to a β-sheet surface implicated in both HIF-2 dimerization and coactivator recruitment. Furthermore, while KG-548 binds exclusively to the β-sheet surface, KG-655 can additionally bind within a water-accessible internal cavity in ARNT PAS-B. Finally, KG-279, while not a coactivator inhibitor, exemplifies ligands that preferentially bind only to the internal cavity. All three ligands promoted ARNT PAS-B homodimerization, albeit to varying degrees. Taken together, our findings provide a comprehensive overview of ARNT PAS-B ligand-binding sites and may guide the development of more potent coactivator inhibitors for cellular and functional studies.

Calculating Protein-Ligand Residence Times through State Predictive Information Bottleneck Based Enhanced Sampling

Understanding drug residence times in target proteins is key to improving drug efficacy and understanding target recognition in biochemistry. While drug residence time is just as important as binding affinity, atomic-level understanding of drug residence times through molecular dynamics (MD) simulations has been difficult primarily due to the extremely long time scales. Recent advances in rare event sampling have allowed us to reach these time scales, yet predicting protein-ligand residence times remains a significant challenge. Here we present a semi-automated protocol to calculate the ligand residence times across 12 orders of magnitude of time scales. In our proposed framework, we integrate a deep learning-based method, the state predictive information bottleneck (SPIB), to learn an approximate reaction coordinate (RC) and use it to guide the enhanced sampling method metadynamics. We demonstrate the performance of our algorithm by applying it to six different protein-ligand complexes with available benchmark residence times, including the dissociation of the widely studied anticancer drug Imatinib (Gleevec) from both wild-type Abl kinase and drug-resistant mutants. We show how our protocol can recover quantitatively accurate residence times, potentially opening avenues for deeper insights into drug development possibilities and ligand recognition mechanisms.

Good Rates From Bad Coordinates: The Exponential Average Time-dependent Rate Approach

Our ability to calculate rate constants of biochemical processes using molecular dynamics simulations is severely limited by the fact that the time scales for reactions, or changes in conformational state, scale exponentially with the relevant free-energy barrier heights. In this work, we improve upon a recently proposed rate estimator that allows us to predict transition times with molecular dynamics simulations biased to rapidly explore one or several collective variables (CVs). This approach relies on the idea that not all bias goes into promoting transitions, and along with the rate, it estimates a concomitant scale factor for the bias termed the "CV biasing efficiency" γ. First, we demonstrate mathematically that our new formulation allows us to derive the commonly used Infrequent Metadynamics (iMetaD) estimator when using a perfect CV, where γ = 1. After testing it on a model potential, we then study the unfolding behavior of a previously well characterized coarse-grained protein, which is sufficiently complex that we can choose many different CVs to bias, but which is sufficiently simple that we are able to compute the unbiased rate directly. For this system, we demonstrate that predictions from our new Exponential Average Time-Dependent Rate (EATR) estimator converge to the true rate constant more rapidly as a function of bias deposition time than does the previous iMetaD approach, even for bias deposition times that are short. We also show that the γ parameter can serve as a good metric for assessing the quality of the biasing coordinate. We demonstrate that these results hold when applying the methods to an atomistic protein folding example. Finally, we demonstrate that our approach works when combining multiple less-than-optimal bias coordinates, and adapt our method to the related "OPES flooding" approach. Overall, our time-dependent rate approach offers a powerful framework for predicting rate constants from biased simulations.

STORMM: Structure and topology replica molecular mechanics for chemical simulations

The Structure and TOpology Replica Molecular Mechanics (STORMM) code is a next-generation molecular simulation engine and associated libraries optimized for performance on fast, vectorized central processor units and graphics processing units (GPUs) with independent memory and tens of thousands of threads. STORMM is built to run thousands of independent molecular mechanical calculations on a single GPU with novel implementations that tune numerical precision, mathematical operations, and scarce on-chip memory resources to optimize throughput. The libraries are built around accessible classes with detailed documentation, supporting fine-grained parallelism and algorithm development as well as copying or swapping groups of systems on and off of the GPU. A primary intention of the STORMM libraries is to provide developers of atomic simulation methods with access to a high-performance molecular mechanics engine with extensive facilities to prototype and develop bespoke tools aimed toward drug discovery applications. In its present state, STORMM delivers molecular dynamics simulations of small molecules and small proteins in implicit solvent with tens to hundreds of times the throughput of conventional codes. The engineering paradigm transforms two of the most memory bandwidth-intensive aspects of condensed-phase dynamics, particle-mesh mapping, and valence interactions, into compute-bound problems for several times the scalability of existing programs. Numerical methods for compressing and streamlining the information present in stored coordinates and lookup tables are also presented, delivering improved accuracy over methods implemented in other molecular dynamics engines. The open-source code is released under the MIT license.

Prebound State Discovered in the Unbinding Pathway of Fluorinated Variants of the Trypsin-BPTI Complex Using Random Acceleration Molecular Dynamics Simulations

The serine protease trypsin forms a tightly bound inhibitor complex with the bovine pancreatic trypsin inhibitor (BPTI). The complex is stabilized by the P1 residue Lys15, which interacts with negatively charged amino acids at the bottom of the S1 pocket. Truncating the P1 residue of wildtype BPTI to α-aminobutyric acid (Abu) leaves a complex with moderate inhibitor strength, which is held in place by additional hydrogen bonds at the protein-protein interface. Fluorination of the Abu residue partially restores the inhibitor strength. The mechanism with which fluorination can restore the inhibitor strength is unknown, and accurate computational investigation requires knowledge of the binding and unbinding pathways. The preferred unbinding pathway is likely to be complex, as encounter states have been described before, and unrestrained umbrella sampling simulations of these complexes suggest additional energetic minima. Here, we use random acceleration molecular dynamics to find a new metastable state in the unbinding pathway of Abu-BPTI variants and wildtype BPTI from trypsin, which we call the prebound state. The prebound state and the fully bound state differ by a substantial shift in the position, a slight shift in the orientation of the BPTI variants, and changes in the interaction pattern. Particularly important is the breaking of three hydrogen bonds around Arg17. Fluorination of the P1 residue lowers the energy barrier of the transition between the fully bound state and prebound state and also lowers the energy minimum of the prebound state. While the effect of fluorination is in general difficult to quantify, here, it is in part caused by favorable stabilization of a hydrogen bond between Gln194 and Cys14. The interaction pattern of the prebound state offers insights into the inhibitory mechanism of BPTI and might add valuable information for the design of serine protease inhibitors.

Implementation of Girsanov Reweighting in OpenMM and Deeptime

Classical molecular dynamics (MD) simulations provide invaluable insights into complex molecular systems but face limitations in capturing phenomena occurring on time scales beyond their reach. To bridge this gap, various enhanced sampling techniques have been developed, which are complemented by reweighting techniques to recover the unbiased dynamics. Girsanov reweighting is a reweighting technique that reweights simulation paths, generated by a stochastic MD integrator, without evoking an effective model of the dynamics. Instead, it calculates the relative path probability density at the time resolution of the MD integrator. Efficient implementation of Girsanov reweighting requires that the reweighting factors are calculated on-the-fly during the simulations and thus needs to be implemented within the MD integrator. Here, we present a comprehensive guide for implementing Girsanov reweighting into MD simulations. We demonstrate the implementation in the MD simulation package OpenMM by extending the library openmmtools. Additionally, we implemented a reweighted Markov state model estimator within the time series analysis package Deeptime.

Augmenting Human Expertise in Weighted Ensemble Simulations through Deep Learning based Information Bottleneck

The weighted ensemble (WE) method stands out as a widely used segment-based sampling technique renowned for its rigorous treatment of kinetics. The WE framework typically involves initially mapping the configuration space onto a low-dimensional collective variable (CV) space and then partitioning it into bins. The efficacy of WE simulations heavily depends on the selection of CVs and binning schemes. The recently proposed State Predictive Information Bottleneck (SPIB) method has emerged as a promising tool for automatically constructing CVs from data and guiding enhanced sampling through an iterative manner. In this work, we advance this data-driven pipeline by incorporating prior expert knowledge. Our hybrid approach combines SPIB-learned CVs to enhance sampling in explored regions with expert-based CVs to guide exploration in regions of interest, synergizing the strengths of both methods. Through benchmarking on alanine dipeptide and chignoin systems, we demonstrate that our hybrid approach effectively guides WE simulations to sample states of interest, and reduces run-to-run variances. Moreover, our integration of the SPIB model also enhances the analysis and interpretation of WE simulation data by effectively identifying metastable states and pathways, and offering direct visualization of dynamics.

Thermal isomerization rates in retinal analogues using Ab-Initio molecular dynamics

For a detailed understanding of chemical processes in nature and industry, we need accurate models of chemical reactions in complex environments. While Eyring transition state theory is commonly used for modeling chemical reactions, it is most accurate for small molecules in the gas phase. A wide range of alternative rate theories exist that can better capture reactions involving complex molecules and environmental effects. However, they require that the chemical reaction is sampled by molecular dynamics simulations. This is a formidable challenge since the accessible simulation timescales are many orders of magnitude smaller than typical timescales of chemical reactions. To overcome these limitations, rare event methods involving enhanced molecular dynamics sampling are employed. In this work, thermal isomerization of retinal is studied using tight-binding density functional theory. Results from transition state theory are compared to those obtained from enhanced sampling. Rates obtained from dynamical reweighting using infrequent metadynamics simulations were in close agreement with those from transition state theory. Meanwhile, rates obtained from application of Kramers' rate equation to a sampled free energy profile along a torsional dihedral reaction coordinate were found to be up to three orders of magnitude higher. This discrepancy raises concerns about applying rate methods to one-dimensional reaction coordinates in chemical reactions.

Identification of Small Molecule Ligand Binding Sites On and In the ARNT PAS-B Domain

Transcription factors are generally challenging to target with small molecule inhibitors due to their structural plasticity and lack of catalytic sites. Notable exceptions include several naturally ligand-regulated transcription factors, including our prior work with the heterodimeric HIF-2 transcription factor which showed that small molecule binding within an internal pocket of the HIF-2α PAS-B domain can disrupt its interactions with its dimerization partner, ARNT. Here, we explore the feasibility of similarly targeting small molecules to the analogous ARNT PAS-B domain itself, potentially opening a promising route to simultaneously modulate several ARNT-mediated signaling pathways. Using solution NMR screening of an in-house fragment library, we previously identified several compounds that bind ARNT PAS-B and, in certain cases, antagonize ARNT association with the TACC3 transcriptional coactivator. However, these ligands have only modest binding affinities, complicating characterization of their binding sites. We address this challenge by combining NMR, MD simulations, and ensemble docking to identify ligand-binding 'hotspots' on and within the ARNT PAS-B domain. Our data indicate that the two ARNT/TACC3 inhibitors, KG-548 and KG-655, bind to a β-sheet surface implicated in both HIF-2 dimerization and coactivator recruitment. Furthermore, while KG-548 binds exclusively to the β-sheet surface, KG-655 can additionally bind within a water-accessible internal cavity in ARNT PAS-B. Finally, KG-279, while not a coactivator inhibitor, exemplifies ligands that preferentially bind only to the internal cavity. All three ligands promoted ARNT PAS-B homodimerization, albeit to varying degrees. Taken together, our findings provide a comprehensive overview of ARNT PAS-B ligand-binding sites and may guide the development of more potent coactivator inhibitors for cellular and functional studies.

Enhanced Sampling with Machine Learning

Molecular dynamics (MD) enables the study of physical systems with excellent spatiotemporal resolution but suffers from severe timescale limitations. To address this, enhanced sampling methods have been developed to improve the exploration of configurational space. However, implementing these methods is challenging and requires domain expertise. In recent years, integration of machine learning (ML) techniques into different domains has shown promise, prompting their adoption in enhanced sampling as well. Although ML is often employed in various fields primarily due to its data-driven nature, its integration with enhanced sampling is more natural with many common underlying synergies. This review explores the merging of ML and enhanced MD by presenting different shared viewpoints. It offers a comprehensive overview of this rapidly evolving field, which can be difficult to stay updated on. We highlight successful strategies such as dimensionality reduction, reinforcement learning, and flow-based methods. Finally, we discuss open problems at the exciting ML-enhanced MD interface.

Mechanism of Daunomycin Intercalation into DNA from Enhanced Sampling Simulations

Daunomycin is a widely used anticancer drug, yet the mechanism underlying how it binds to DNA remains contested. 469 all-atom trajectories of daunomycin binding to the DNA oligonucleotide d(GCG CAC GTG CGC) were collected using weighted ensemble (WE)-enhanced sampling. Mechanistic insights were revealed through analysis of the ensemble of trajectories. Initially, the binding process involves a ubiquitous hydrogen bond between the DNA backbone and the NH3+ group on daunomycin. During the binding process, most trajectories exhibited similar structural changes to DNA, including DNA base pair rise, bending, and minor groove width changes. Variability within the ensemble of binding trajectories illuminates differences in the orientation of daunomycin as it initially intercalates; around 10% of trajectories needed minimal rearrangement from intercalation to reaching the fully bound configuration, whereas most needed an additional 1-5 ns to rearrange. The results here emphasize the utility of generating an ensemble of trajectories to discern biomolecular binding mechanisms.

Dynamic allostery in the peptide/MHC complex enables TCR neoantigen selectivity

The inherent cross-reactivity of the T cell receptor (TCR) is balanced by high specificity, which often manifests in confounding ways not easily interpretable from static structures. We show here that TCR discrimination between an HLA-A*03:01 (HLA-A3)-restricted public neoantigen derived from mutant PIK3CA and its wild-type (WT) counterpart emerges from motions within the HLA binding groove that vary with the identity of the peptide's first primary anchor. The motions form a dynamic gate that in the complex with the WT peptide impedes a large conformational change required for TCR binding. The more rigid neoantigen is insusceptible to this limiting dynamic, and with the gate open, is able to transit its central tryptophan residue underneath the peptide backbone to the contralateral side of the HLA-A3 peptide binding groove, facilitating TCR binding. Our findings reveal a novel mechanism driving TCR specificity for a cancer neoantigen that is rooted in the dynamic and allosteric nature of peptide/MHC-I complexes, with implications for resolving long-standing and often confounding questions about the determinants of T cell specificity.

Exploring Biomolecular Conformational Dynamics with Polarizable Force Field AMOEBA and Enhanced Sampling Method Milestoning

Conformational dynamics play a crucial role in determining the behavior of the biomolecules. Polarizable force fields, such as AMOEBA, can accurately capture electrostatic interactions underlying the conformational space. However, applying a polarizable force field in molecular dynamics (MD) simulations can be computationally expensive, especially in studying long-time-scale dynamics. To overcome this challenge, we incorporated the AMOEBA potential with Milestoning, an enhanced sampling method in this work. This integration allows us to efficiently sample the rare and important conformational states of a biomolecule by using many short and independent molecular dynamics trajectories with the AMOEBA force field. We applied this method to investigate the conformational dynamics of alanine dipeptide, DNA, and RNA A-B form conversion. Well-converged thermodynamic and kinetic properties were obtained, including the free energy difference, mean first passage time, and critical transitions between states. Our results demonstrate the power of integrating polarizable force fields with enhanced sampling methods in quantifying the thermodynamic and kinetic properties of biomolecules at the atomic level.

WEDAP: A Python Package for Streamlined Plotting of Molecular Simulation Data

Given the growing interest in path sampling methods for extending the timescales of molecular dynamics (MD) simulations, there has been great interest in software tools that streamline the generation of plots for monitoring the progress of large-scale simulations. Here, we present the WEDAP Python package for simplifying the analysis of data generated from either conventional MD simulations or the weighted ensemble (WE) path sampling method, as implemented in the widely used WESTPA software package. WEDAP facilitates (i) the parsing of WE simulation data stored in highly compressed, hierarchical HDF5 files, and (ii) incorporates trajectory weights from WE simulations into all generated plots. Our Python package consists of multiple user-friendly interfaces: a command-line interface, a graphical user interface, and a Python application programming interface. We demonstrate the plotting features of WEDAP through a series of examples using data from WE and conventional MD simulations that focus on the HIV-1 capsid protein C-terminal domain dimer as a showcase system. The source code for WEDAP is freely available on GitHub at https://github.com/chonglab-pitt/wedap .

Quantitative Prediction of Dissociation Rates of PYK2 Ligands Using Umbrella Sampling and Milestoning

We used umbrella sampling and the milestoning simulation method to study the dissociation of multiple ligands from protein kinase PYK2. The activation barriers obtained from the potential of mean force of the umbrella sampling simulations correlated well with the experimental dissociation rates. Using the zero-temperature string method, we obtained optimized paths along the free-energy surfaces for milestoning simulations of three ligands with a similar structure. The milestoning simulations gave an absolute dissociation rate within 2 orders of magnitude of the experimental value for two ligands but at least 3 orders of magnitude too high for the third. Despite the similarity in their structures, the ligands took different pathways to exit from the binding site of PYK2, making contact with different sets of residues. In addition, the protein experienced different conformational changes for dissociation of the three ligands.

Revisiting Textbook Azide-Clock Reactions: A "Propeller-Crawling" Mechanism Explains Differences in Rates

An ongoing challenge to chemists is the analysis of pathways and kinetics for chemical reactions in solution, including transient structures between the reactants and products that are difficult to resolve using laboratory experiments. Here, we enabled direct molecular dynamics simulations of a textbook series of chemical reactions on the hundreds of ns to μs time scale using the weighted ensemble (WE) path sampling strategy with hybrid quantum mechanical/molecular mechanical (QM/MM) models. We focused on azide-clock reactions involving addition of an azide anion to each of three long-lived trityl cations in an acetonitrile-water solvent mixture. Results reveal a two-step mechanism: (1) diffusional collision of reactants to form an ion-pair intermediate; (2) "activation" or rearrangement of the intermediate to the product. Our simulations yield not only reaction rates that are within error of experiment but also rates for individual steps, indicating the activation step as rate-limiting for all three cations. Further, the trend in reaction rates is due to dynamical effects, i.e., differing extents of the azide anion "crawling" along the cation's phenyl-ring "propellers" during the activation step. Our study demonstrates the power of analyzing pathways and kinetics to gain insights on reaction mechanisms, underscoring the value of including WE and other related path sampling strategies in the modern toolbox for chemists.

One Descriptor to Fold Them All: Harnessing Intuition and Machine Learning to Identify Transferable Lasso Peptide Reaction Coordinates

Identifying optimal reaction coordinates for complex conformational changes and protein folding remains an outstanding challenge. This study combines collective variable (CV) discovery based on chemical intuition and machine learning with enhanced sampling to converge the folding free energy landscape of lasso peptides, a unique class of natural products with knot-like tertiary structures. This knotted scaffold imparts remarkable stability, making lasso peptides resistant to proteolytic degradation, thermal denaturation, and extreme pH conditions. Although their direct synthesis would enable therapeutic design, it has not yet been possible due to the improbable occurrence of spontaneous lasso folding. Thus, simulations characterizing the folding propensity are needed to identify strategies for increasing access to the lasso architecture by stabilizing the pre-lasso ensemble before isopeptide bond formation. Herein, harmonic linear discriminant analysis (HLDA) is combined with metadynamics-enhanced sampling to discover CVs capable of distinguishing the pre-lasso fold and converging the folding propensity. Intuitive CVs are compared to iterative rounds of HLDA to identify CVs that not only accomplish these goals for one lasso peptide but also seem to be transferable to others, establishing a protocol for the identification of folding reaction coordinates for lasso peptides.

Calculating Protein-Ligand Residence Times Through State Predictive Information Bottleneck based Enhanced Sampling

Understanding drug residence times in target proteins is key to improving drug efficacy and understanding target recognition in biochemistry. While drug residence time is just as important as binding affinity, atomic-level understanding of drug residence times through molecular dynamics (MD) simulations has been difficult primarily due to the extremely long timescales. Recent advances in rare event sampling have allowed us to reach these timescales, yet predicting protein-ligand residence times remains a significant challenge. Here we present a semi-automated protocol to calculate the ligand residence times across 12 orders of magnitudes of timescales. In our proposed framework, we integrate a deep learning-based method, the state predictive information bottleneck (SPIB), to learn an approximate reaction coordinate (RC) and use it to guide the enhanced sampling method metadynamics. We demonstrate the performance of our algorithm by applying it to six different protein-ligand complexes with available benchmark residence times, including the dissociation of the widely studied anti-cancer drug Imatinib (Gleevec) from both wild-type Abl kinase and drug-resistant mutants. We show how our protocol can recover quantitatively accurate residence times, potentially opening avenues for deeper insights into drug development possibilities and ligand recognition mechanisms.

PaCS-Toolkit: Optimized Software Utilities for Parallel Cascade Selection Molecular Dynamics (PaCS-MD) Simulations and Subsequent Analyses

Parallel cascade selection molecular dynamics (PaCS-MD) is an enhanced conformational sampling method conducted as a "repetition of time leaps in parallel worlds", comprising cycles of multiple molecular dynamics (MD) simulations performed in parallel and selection of the initial structures of MDs for the next cycle. We developed PaCS-Toolkit, an optimized software utility enabling the use of different MD software and trajectory analysis tools to facilitate the execution of the PaCS-MD simulation and analyze the obtained trajectories, including the preparation for the subsequent construction of the Markov state model. PaCS-Toolkit is coded with Python, is compatible with various computing environments, and allows for easy customization by editing the configuration file and specifying the MD software and analysis tools to be used. We present the software design of PaCS-Toolkit and demonstrate applications of PaCS-MD variations: original targeted PaCS-MD to peptide folding; rmsdPaCS-MD to protein domain motion; and dissociation PaCS-MD to ligand dissociation from adenosine A2A receptor.

Distilling dynamical knowledge from stochastic reaction networks

Stochastic reaction networks are widely used in the modeling of stochastic systems across diverse domains such as biology, chemistry, physics, and ecology. However, the comprehension of the dynamic behaviors inherent in stochastic reaction networks is a formidable undertaking, primarily due to the exponential growth in the number of possible states or trajectories as the state space dimension increases. In this study, we introduce a knowledge distillation method based on reinforcement learning principles, aimed at compressing the dynamical knowledge encoded in stochastic reaction networks into a singular neural network construct. The trained neural network possesses the capability to accurately predict the state conditional joint probability distribution that corresponds to the given query contexts, when prompted with rate parameters, initial conditions, and time values. This obviates the need to track the dynamical process, enabling the direct estimation of normalized state and trajectory probabilities, without necessitating the integration over the complete state space. By applying our method to representative examples, we have observed a high degree of accuracy in both multimodal and high-dimensional systems. Additionally, the trained neural network can serve as a foundational model for developing efficient algorithms for parameter inference and trajectory ensemble generation. These results collectively underscore the efficacy of our approach as a universal means of distilling knowledge from stochastic reaction networks. Importantly, our methodology also spotlights the potential utility in harnessing a singular, pretrained, large-scale model to encapsulate the solution space underpinning a wide spectrum of stochastic dynamical systems.

Comparison of Accuracy and Efficiency of Milestoning Variants: Introducing Buffer Milestoning

The Milestoning algorithm is a method for long-time molecular dynamics simulations. It enables the sampling of rare events. The precise calculations of observables depend on accurately determining the first hitting point distribution (FHPD) for each milestone. There is no analytical expression for FHPD, which is estimated numerically. Several variants of Milestoning offer approximations to the FHPD. Here, we examine in detail the FHPD of an exact calculation and Milestoning variants. We also introduce a new version of the Milestoning algorithm, buffer Milestoning, with a comparable cost to conventional Milestoning but higher accuracy. We use the mean first passage time and the free energy to assess the simulation quality, and we compare the accuracy and efficiency of buffer Milestoning to exact calculations, conventional Milestoning, local-passage-time-weighted Milestoning, Markovian Milestoning with Voronoi tessellation, and exact Milestoning. Conventional Milestoning requires milestone decorrelation. If this condition is not satisfied, it is the least accurate approach of all the techniques we examined. We conclude that for a small increase in cost compared to conventional Milestoning, buffer Milestoning provides accurate results for a range of problems, including more correlated milestones and is, therefore, versatile compared to other variants. Local-passage-time-weighted Milestoning provides accuracy similar to that of buffer Milestoning but at an increased simulation cost. Markovian Milestoning with Voronoi tessellation is the most accurate compared with other approximations, but it is less stable for high barriers and more expensive.

Enhanced Sampling of Buried Charges in Free Energy Calculations Using Replica Exchange with Charge Tempering

Buried ionizable groups in proteins often play important structural and functional roles. However, it is generally challenging to study the detailed molecular mechanisms solely based on experimental measurements. Free energy calculations using atomistic simulations, on the other hand, complement experimental studies and can provide high temporal and spatial resolution information that can lead to mechanistic insights. Nevertheless, it is also well recognized that sufficient sampling of such atomistic simulations can be challenging, considering that structural changes related to the buried charges may be very slow. In the present study, we describe a simple but effective enhanced sampling technique called replica exchange with charge tempering (REChgT) with a novel free energy method, multisite λ dynamics (MSλD), to study two systems containing buried charges, pKa prediction of a small molecule, orotate, in complex with the dihydroorotate dehydrogenase, and relative stability of a Glu-Lys pair buried in the hydrophobic core of two variants of Staphylococcal nuclease. Compared to the original MSλD simulations, the usage of REChgT dramatically increases sampling in both conformational and alchemical spaces, which directly translates into a significant reduction of wall time to converge the free energy calculations. This study highlights the importance of sufficient sampling toward developing improved free energy methods.

Elucidating molecular mechanisms of protoxin-II state-specific binding to the human NaV1.7 channel

Human voltage-gated sodium (hNaV) channels are responsible for initiating and propagating action potentials in excitable cells, and mutations have been associated with numerous cardiac and neurological disorders. hNaV1.7 channels are expressed in peripheral neurons and are promising targets for pain therapy. The tarantula venom peptide protoxin-II (PTx2) has high selectivity for hNaV1.7 and is a valuable scaffold for designing novel therapeutics to treat pain. Here, we used computational modeling to study the molecular mechanisms of the state-dependent binding of PTx2 to hNaV1.7 voltage-sensing domains (VSDs). Using Rosetta structural modeling methods, we constructed atomistic models of the hNaV1.7 VSD II and IV in the activated and deactivated states with docked PTx2. We then performed microsecond-long all-atom molecular dynamics (MD) simulations of the systems in hydrated lipid bilayers. Our simulations revealed that PTx2 binds most favorably to the deactivated VSD II and activated VSD IV. These state-specific interactions are mediated primarily by PTx2's residues R22, K26, K27, K28, and W30 with VSD and the surrounding membrane lipids. Our work revealed important protein-protein and protein-lipid contacts that contribute to high-affinity state-dependent toxin interaction with the channel. The workflow presented will prove useful for designing novel peptides with improved selectivity and potency for more effective and safe treatment of pain.

Integrating cellular electron microscopy with multimodal data to explore biology across space and time

Biology spans a continuum of length and time scales. Individual experimental methods only glimpse discrete pieces of this spectrum but can be combined to construct a more holistic view. In this Review, we detail the latest advancements in volume electron microscopy (vEM) and cryo-electron tomography (cryo-ET), which together can visualize biological complexity across scales from the organization of cells in large tissues to the molecular details inside native cellular environments. In addition, we discuss emerging methodologies for integrating three-dimensional electron microscopy (3DEM) imaging with multimodal data, including fluorescence microscopy, mass spectrometry, single-particle analysis, and AI-based structure prediction. This multifaceted approach fills gaps in the biological continuum, providing functional context, spatial organization, molecular identity, and native interactions. We conclude with a perspective on incorporating diverse data into computational simulations that further bridge and extend length scales while integrating the dimension of time.

Biomolecular dynamics in the 21st century

The relevance of motions in biological macromolecules has been clear since the early structural analyses of proteins by X-ray crystallography. Computer simulations have been applied to provide a deeper understanding of the dynamics of biological macromolecules since 1976, and are now a standard tool in many labs working on the structure and function of biomolecules. In this mini-review we highlight some areas of current interest and active development for simulations, in particular all-atom molecular dynamics simulations.

Combining stochastic resetting with Metadynamics to speed-up molecular dynamics simulations

Metadynamics is a powerful method to accelerate molecular dynamics simulations, but its efficiency critically depends on the identification of collective variables that capture the slow modes of the process. Unfortunately, collective variables are usually not known a priori and finding them can be very challenging. We recently presented a collective variables-free approach to enhanced sampling using stochastic resetting. Here, we combine the two methods, showing that it can lead to greater acceleration than either of them separately. We also demonstrate that resetting Metadynamics simulations performed with suboptimal collective variables can lead to speedups comparable with those obtained with optimal collective variables. Therefore, applying stochastic resetting can be an alternative to the challenging task of improving suboptimal collective variables, at almost no additional computational cost. Finally, we propose a method to extract unbiased mean first-passage times from Metadynamics simulations with resetting, resulting in an improved tradeoff between speedup and accuracy. This work enables combining stochastic resetting with other enhanced sampling methods to accelerate a broad range of molecular simulations.

Effects of altered backbone composition on the folding kinetics and mechanism of an ultrafast-folding protein

Sequence-encoded protein folding is a ubiquitous biological process that has been successfully engineered in a range of oligomeric molecules with artificial backbone chemical connectivity. A remarkable aspect of protein folding is the contrast between the rapid rates at which most sequences in nature fold and the vast number of conformational states possible in an unfolded chain with hundreds of rotatable bonds. Research efforts spanning several decades have sought to elucidate the fundamental chemical principles that dictate the speed and mechanism of natural protein folding. In contrast, little is known about how protein mimetic entities transition between an unfolded and folded state. Here, we report effects of altered backbone connectivity on the folding kinetics and mechanism of the B domain of Staphylococcal protein A (BdpA), an ultrafast-folding sequence. A combination of experimental biophysical analysis and atomistic molecular dynamics simulations performed on the prototype protein and several heterogeneous-backbone variants reveal the interplay among backbone flexibility, folding rates, and structural details of the transition state ensemble. Collectively, these findings suggest a significant degree of plasticity in the mechanisms that can give rise to ultrafast folding in the BdpA sequence and provide atomic level insights into how protein mimetic chains adopt an ordered folded state.

Revolutionizing Peptide-Based Drug Discovery: Advances in the Post-AlphaFold Era

Peptide-based drugs offer high specificity, potency, and selectivity. However, their inherent flexibility and differences in conformational preferences between their free and bound states create unique challenges that have hindered progress in effective drug discovery pipelines. The emergence of AlphaFold (AF) and Artificial Intelligence (AI) presents new opportunities for enhancing peptide-based drug discovery. We explore recent advancements that facilitate a successful peptide drug discovery pipeline, considering peptides' attractive therapeutic properties and strategies to enhance their stability and bioavailability. AF enables efficient and accurate prediction of peptide-protein structures, addressing a critical requirement in computational drug discovery pipelines. In the post-AF era, we are witnessing rapid progress with the potential to revolutionize peptide-based drug discovery such as the ability to rank peptide binders or classify them as binders/non-binders and the ability to design novel peptide sequences. However, AI-based methods are struggling due to the lack of well-curated datasets, for example to accommodate modified amino acids or unconventional cyclization. Thus, physics-based methods, such as docking or molecular dynamics simulations, continue to hold a complementary role in peptide drug discovery pipelines. Moreover, MD-based tools offer valuable insights into binding mechanisms, as well as the thermodynamic and kinetic properties of complexes. As we navigate this evolving landscape, a synergistic integration of AI and physics-based methods holds the promise of reshaping the landscape of peptide-based drug discovery.

Thrombin - A Molecular Dynamics Perspective

Thrombin is a crucial enzyme involved in blood coagulation, essential for maintaining circulatory system integrity and preventing excessive bleeding. However, thrombin is also implicated in pathological conditions such as thrombosis and cancer. Despite the application of various experimental techniques, including X-ray crystallography, NMR spectroscopy, and HDXMS, none of these methods can precisely detect thrombin's dynamics and conformational ensembles at high spatial and temporal resolution. Fortunately, molecular dynamics (MD) simulation, a computational technique that allows the investigation of molecular functions and dynamics in atomic detail, can be used to explore thrombin behavior. This review summarizes recent MD simulation studies on thrombin and its interactions with other biomolecules. Specifically, the 17 studies discussed here provide insights into thrombin's switch between 'slow' and 'fast' forms, active and inactive forms, the role of Na+ binding, the effects of light chain mutation, and thrombin's interactions with other biomolecules. The findings of these studies have significant implications for developing new therapies for thrombosis and cancer. By understanding thrombin's complex behavior, researchers can design more effective drugs and treatments that target thrombin.

Tracking conformational transitions of the gonadotropin hormone receptors in a bilayer of (SDPC) poly-unsaturated lipids from all-atom molecular dynamics simulations

Glycoprotein hormone receptors [thyrotropin (TSHR), luteinizing hormone/chorionic gonadotropin (LHCGR), and follicle stimulating hormone (FSHR) receptors] are rhodopsin-like G protein-coupled receptors. These receptors display common structural features including a prominent extracellular domain with leucine-rich repeats (LRR) stabilized by β-sheets and a long and flexible loop known as the hinge region (HR), and a transmembrane (TM) domain with seven α-helices interconnected by intra- and extracellular loops. Binding of the ligand to the LRR resembles a hand coupling transversally to the α- and β-subunits of the hormone, with the thumb being the HR. The structure of the FSH-FSHR complex suggests an activation mechanism in which Y335 at the HR binds into a pocket between the α- and β-chains of the hormone, leading to an adjustment of the extracellular loops. In this study, we performed molecular dynamics (MD) simulations to identify the conformational changes of the FSHR and LHCGR. We set up a FSHR structure as predicted by AlphaFold (AF-P23945); for the LHCGR structure we took the cryo-electron microscopy structure for the active state (PDB:7FII) as initial coordinates. Specifically, the flexibility of the HR domain and the correlated motions of the LRR and TM domain were analyzed. From the conformational changes of the LRR, TM domain, and HR we explored the conformational landscape by means of MD trajectories in all-atom approximation, including a membrane of polyunsaturated phospholipids. The distances and procedures here defined may be useful to propose reaction coordinates to describe diverse processes, such as the active-to-inactive transition, and to identify intermediaries suited for allosteric regulation and biased binding to cellular transducers in a selective activation strategy.

Molecular Free Energies, Rates, and Mechanisms from Data-Efficient Path Sampling Simulations

Molecular dynamics is a powerful tool for studying the thermodynamics and kinetics of complex molecular events. However, these simulations can rarely sample the required time scales in practice. Transition path sampling overcomes this limitation by collecting unbiased trajectories and capturing the relevant events. Moreover, the integration of machine learning can boost the sampling while simultaneously learning a quantitative representation of the mechanism. Still, the resulting trajectories are by construction non-Boltzmann-distributed, preventing the calculation of free energies and rates. We developed an algorithm to approximate the equilibrium path ensemble from machine-learning-guided path sampling data. At the same time, our algorithm provides efficient sampling, mechanism, free energy, and rates of rare molecular events at a very moderate computational cost. We tested the method on the folding of the mini-protein chignolin. Our algorithm is straightforward and data-efficient, opening the door to applications in many challenging molecular systems.

LPATH: A Semiautomated Python Tool for Clustering Molecular Pathways

The pathways by which a molecular process transitions to a target state are highly sought-after as direct views of a transition mechanism. While great strides have been made in the physics-based simulation of such pathways, the analysis of these pathways can be a major challenge due to their diversity and variable lengths. Here, we present the LPATH Python tool, which implements a semiautomated method for linguistics-assisted clustering of pathways into distinct classes (or routes). This method involves three steps: 1) discretizing the configurational space into key states, 2) extracting a text-string sequence of key visited states for each pathway, and 3) pairwise matching of pathways based on a text-string similarity score. To circumvent the prohibitive memory requirements of the first step, we have implemented a general two-stage method for clustering conformational states that exploits machine learning. LPATH is primarily designed for use with the WESTPA software for weighted ensemble simulations; however, the tool can also be applied to conventional simulations. As demonstrated for the C7eq to C7ax conformational transition of the alanine dipeptide, LPATH provides physically reasonable classes of pathways and corresponding probabilities.

Adaptive Sampling Methods for Molecular Dynamics in the Era of Machine Learning

Molecular dynamics (MD) simulations are fundamental computational tools for the study of proteins and their free energy landscapes. However, sampling protein conformational changes through MD simulations is challenging due to the relatively long time scales of these processes. Many enhanced sampling approaches have emerged to tackle this problem, including biased sampling and path-sampling methods. In this Perspective, we focus on adaptive sampling algorithms. These techniques differ from other approaches because the thermodynamic ensemble is preserved and the sampling is enhanced solely by restarting MD trajectories at particularly chosen seeds rather than introducing biasing forces. We begin our treatment with an overview of theoretically transparent methods, where we discuss principles and guidelines for adaptive sampling. Then, we present a brief summary of select methods that have been applied to realistic systems in the past. Finally, we discuss recent advances in adaptive sampling methodology powered by deep learning techniques, as well as their shortcomings.

From Byte to Bench to Bedside: Molecular Dynamics Simulations and Drug Discovery

Molecular dynamics (MD) simulations and computer-aided drug design (CADD) have advanced substantially over the past two decades, thanks to continuous computer hardware and software improvements. Given these advancements, MD simulations are poised to become even more powerful tools for investigating the dynamic interactions between potential small-molecule drugs and their target proteins, with significant implications for pharmacological research.

Bispidine as a promising scaffold for designing molecular machines

The development of artificial molecular machines is a challenging endeavor. Herein, we have synthesized a series of bispidine diamides D1-D6 that exhibit rotation reminiscent of a motor motion. Dynamic NMR, X-ray diffraction, quantum mechanical calculations, and molecular dynamics simulations provided insights into their rotational dynamics. All the diamides D1-D6 exhibited mutually independent rotation around the two bispidine arms. However, the rate of rotation and the presence or absence of directionality in amide bond rotation were found to depend on the solvent, temperature, and nature of substitution on the amide carbonyl. These engineered systems may aid in the development of biologically relevant synthetic molecular motors. Studies on homochiral and heterochiral bispidine-peptides revealed that the direction of rotation can be controlled by chirality and the nature of the amino acid.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Error-Controlled Coarse-Graining Dynamics with Mean-Field Randomization

In order to comprehend the stochastic behavior of biological systems, it is essential to accurately infer the dynamics of chemical reaction networks. However, computation of the likelihood remains a bottleneck. In this study, we propose the mean-field randomization procedure as a means of efficiently generating error-controlled coarse-graining dynamics. The error is measured by mutual information between the generated trajectories and the coarse-graining procedure. We demonstrate that the exact dynamics can be recovered by resampling, which eliminates the correlation between the dynamics and the procedure. We developed three algorithms to efficiently generate exact or coarse-graining trajectories within a specified error range. By subjecting our algorithms to testing on chemical reaction systems of varying complexities and scales, we observe that they outperform existing state-of-the-art algorithms, and the efficiency of coarse-graining trajectory generation is only weakly dependent on system scales.

Mutexa: A Computational Ecosystem for Intelligent Protein Engineering

Protein engineering holds immense promise in shaping the future of biomedicine and biotechnology. This Review focuses on our ongoing development of Mutexa, a computational ecosystem designed to enable "intelligent protein engineering". In this vision, researchers will seamlessly acquire sequences of protein variants with desired functions as biocatalysts, therapeutic peptides, and diagnostic proteins through a finely-tuned computational machine, akin to Amazon Alexa's role as a versatile virtual assistant. The technical foundation of Mutexa has been established through the development of a database that combines and relates enzyme structures and their respective functions (e.g., IntEnzyDB), workflow software packages that enable high-throughput protein modeling (e.g., EnzyHTP and LassoHTP), and scoring functions that map the sequence-structure-function relationship of proteins (e.g., EnzyKR and DeepLasso). We will showcase the applications of these tools in benchmarking the convergence conditions of enzyme functional descriptors across mutants, investigating protein electrostatics and cavity distributions in SAM-dependent methyltransferases, and understanding the role of nonelectrostatic dynamic effects in enzyme catalysis. Finally, we will conclude by addressing the future steps and fundamental challenges in our endeavor to develop new Mutexa applications that assist the identification of beneficial mutants in protein engineering.

Understanding the free-energy landscape of phase separation in lipid bilayers using molecular dynamics

Liquid-liquid phase separation inside the cell often results in biological condensates that can critically affect cell homeostasis. Such phase separation events occur in multiple parts of cells, including the cell membranes, where the "lipid raft" hypothesis posits the formation of ordered domains floating in a sea of disordered lipids. The resulting lipid domains often have functional roles. However, the thermodynamics of lipid phase separation and their resulting mechanistic effects on cell function and dysfunction are poorly understood. Understanding such complex phenomena in cell membranes, with their diverse lipid compositions, is exceptionally difficult. For these reasons, simple model systems that can recapitulate similar behavior are widely used to study this phenomenon. Despite these simplifications, the timescale and length scales of domain formation pose a challenge for molecular dynamics (MD) simulations. Thus, most MD studies focus on spontaneous lipid phase separation-essentially measuring the sign (but not the amplitude) of the free-energy change upon separation-rather than directly interrogating the thermodynamics. Here, we propose a proof-of-concept pipeline that can directly measure this free energy by combining coarse-grained MD with enhanced sampling protocols using a novel collective variable. This approach will be a useful tool to help connect the thermodynamics of phase separation with the mechanistic insights already available from MD simulations.

Elucidating protein-ligand binding kinetics based on returning probability theory

The returning probability (RP) theory, a rigorous diffusion-influenced reaction theory, enables us to analyze the binding process systematically in terms of thermodynamics and kinetics using molecular dynamics (MD) simulations. Recently, the theory was extended to atomistically describe binding processes by adopting the host-guest interaction energy as the reaction coordinate. The binding rate constants can be estimated by computing the thermodynamic and kinetic properties of the reactive state existing in the binding processes. Here, we propose a methodology based on the RP theory in conjunction with the energy representation theory of solution, applicable to complex binding phenomena, such as protein-ligand binding. The derived scheme of calculating the equilibrium constant between the reactive and dissociate states, required in the RP theory, can be used for arbitrary types of reactive states. We apply the present method to the bindings of small fragment molecules [4-hydroxy-2-butanone (BUT) and methyl methylthiomethyl sulphoxide (DSS)] to FK506 binding protein (FKBP) in an aqueous solution. Estimated binding rate constants are consistent with those obtained from long-timescale MD simulations. Furthermore, by decomposing the rate constants to the thermodynamic and kinetic contributions, we clarify that the higher thermodynamic stability of the reactive state for DSS causes the faster binding kinetics compared with BUT.

Direct observation of negative cooperativity in a detoxification enzyme at the atomic level by Electron Paramagnetic Resonance spectroscopy and simulation

The catalytic activity of human glutathione S-transferase A1-1 (hGSTA1-1), a homodimeric detoxification enzyme, is dependent on the conformational dynamics of a key C-terminal helix α9 in each monomer. However, the structural details of how the two monomers interact upon binding of substrates is not well understood and the structure of the ligand-free state of the hGSTA1-1 homodimer has not been resolved. Here, we used a combination of electron paramagnetic resonance (EPR) distance measurements and weighted ensemble (WE) simulations to characterize the conformational ensemble of the ligand-free state at the atomic level. EPR measurements reveal a broad distance distribution between a pair of Cu(II) labels in the ligand-free state that gradually shifts and narrows as a function of increasing ligand concentration. These shifts suggest changes in the relative positioning of the two α9 helices upon ligand binding. WE simulations generated unbiased pathways for the seconds-timescale transition between alternate states of the enzyme, leading to the generation of atomically detailed structures of the ligand-free state. Notably, the simulations provide direct observations of negative cooperativity between the monomers of hGSTA1-1, which involve the mutually exclusive docking of α9 in each monomer as a lid over the active site. We identify key interactions between residues that lead to this negative cooperativity. Negative cooperativity may be essential for interaction of hGSTA1-1 with a wide variety of toxic substrates and their subsequent neutralization. More broadly, this work demonstrates the power of integrating EPR distances with WE rare-events sampling strategy to gain mechanistic information on protein function at the atomic level.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.