WESTPA 2.0: High-Performance Upgrades for Weighted Ensemble Simulations and Analysis of Longer-Timescale Applications

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.

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.

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.

Amber free energy tools: Interoperable software for free energy simulations using generalized quantum mechanical/molecular mechanical and machine learning potentials

We report the development and testing of new integrated cyberinfrastructure for performing free energy simulations with generalized hybrid quantum mechanical/molecular mechanical (QM/MM) and machine learning potentials (MLPs) in Amber. The Sander molecular dynamics program has been extended to leverage fast, density-functional tight-binding models implemented in the DFTB+ and xTB packages, and an interface to the DeePMD-kit software enables the use of MLPs. The software is integrated through application program interfaces that circumvent the need to perform "system calls" and enable the incorporation of long-range Ewald electrostatics into the external software's self-consistent field procedure. The infrastructure provides access to QM/MM models that may serve as the foundation for QM/MM-ΔMLP potentials, which supplement the semiempirical QM/MM model with a MLP correction trained to reproduce ab initio QM/MM energies and forces. Efficient optimization of minimum free energy pathways is enabled through a new surface-accelerated finite-temperature string method implemented in the FE-ToolKit package. Furthermore, we interfaced Sander with the i-PI software by implementing the socket communication protocol used in the i-PI client-server model. The new interface with i-PI allows for the treatment of nuclear quantum effects with semiempirical QM/MM-ΔMLP models. The modular interoperable software is demonstrated on proton transfer reactions in guanine-thymine mispairs in a B-form deoxyribonucleic acid helix. The current work represents a considerable advance in the development of modular software for performing free energy simulations of chemical reactions that are important in a wide range of applications.

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.

Computational insights into the conformational transition of STING: Mechanistic, energetic considerations, and the influence of crucial mutations

STING (stimulator of interferon genes) is a crucial protein in the innate immune system's response to viral and bacterial infections. In this study, we investigated the mechanistic and energetic mechanism of the conformational transition process of STING activated by cGAMP binding. We found that the STING connector region undergoes an energetically unfavorable rotation during this process, which is compensated by the favorable interaction between cGAMP and the STING ligand binding domain. We further studied several disease-causing mutations and found that the V155 M mutation facilitates a smoother transition in the STING connector region. However, the V147L mutation exhibits unfavorable conformational transition energy, suggesting it may hinder STING activation pathway that relies on connector region rotation. Despite being labeled as hyperactive, the widespread prevalence of V147L/V147I mutations across species implies a neutral character, indicating complexity in its role. Overall, our analysis deepens the understanding of STING activation within the connector region, and targeting this region with compounds may provide an alternative approach to interfering with STING's function.

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 .

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.

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.

Domain coupling in activation of a family C GPCR

The G protein-coupled metabotropic glutamate receptors form homodimers and heterodimers with highly diverse responses to glutamate and varying physiological function. The molecular basis for this diversity remains poorly delineated. We employ molecular dynamics, single-molecule spectroscopy, and hydrogen-deuterium exchange to dissect the pathway of activation triggered by glutamate. We find that activation entails multiple loosely coupled steps and identify a novel pre-active intermediate whose transition to the active state forms dimer interactions that set signaling efficacy. Such subunit interactions generate functional diversity that differs across homodimers and heterodimers. The agonist-bound receptor is remarkably dynamic, with low occupancy of G protein-coupling conformations, providing considerable headroom for modulation of the landscape by allosteric ligands. Sites of sequence diversity within the dimerization interface and diverse coupling between activation rearrangements may contribute to precise decoding of glutamate signals and transients over broad spatial and temporal scales.

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.

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.

LPATH: A semi-automated 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 semi-automated 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 alanine dipeptide, LPATH provides physically reasonable classes of pathways and corresponding probabilities.

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.

Uncovering the Kinetic Characteristics and Degradation Preference of PROTAC Systems with Advanced Theoretical Analyses

Proteolysis-targeting chimeras (PROTACs), which can selectively induce the degradation of target proteins, represent an attractive technology in drug discovery. A large number of PROTACs have been reported, but due to the complicated structural and kinetic characteristics of the target-PROTAC-E3 ligase ternary interaction process, the rational design of PROTACs is still quite challenging. Here, we characterized and analyzed the kinetic mechanism of MZ1, a PROTAC that targets the bromodomain (BD) of the bromodomain and extra terminal (BET) protein (Brd2, Brd3, or Brd4) and von Hippel-Lindau E3 ligase (VHL), from the kinetic and thermodynamic perspectives of view by using enhanced sampling simulations and free energy calculations. The simulations yielded satisfactory predictions on the relative residence time and standard binding free energy (r_p > 0.9) for MZ1 in different Brd^BD-MZ1-VHL ternary complexes. Interestingly, the simulation of the PROTAC ternary complex disintegration illustrates that MZ1 tends to remain on the surface of VHL with the BD proteins dissociating alone without a specific dissociation direction, indicating that the PROTAC prefers more to bind with E3 ligase at the first step in the formation of the target-PROTAC-E3 ligase ternary complex. Further exploration of the degradation difference of MZ1 in different Brd systems shows that the PROTAC with higher degradation efficiency tends to leave more lysine exposed on the target protein, which is guaranteed by the stability (binding affinity) and durability (residence time) of the target-PROTAC-E3 ligase ternary complex. It is quite possible that the underlying binding characteristics of the Brd^BD-MZ1-VHL systems revealed by this study may be shared by different PROTAC systems as a general rule, which may accelerate rational PROTAC design with higher degradation efficiency.

DeepWEST: Deep Learning of Kinetic Models with the Weighted Ensemble Simulation Toolkit for Enhanced Sampling

Recent advances in computational power and algorithms have enabled molecular dynamics (MD) simulations to reach greater time scales. However, for observing conformational transitions associated with biomolecular processes, MD simulations still have limitations. Several enhanced sampling techniques seek to address this challenge, including the weighted ensemble (WE) method, which samples transitions between metastable states using many weighted trajectories to estimate kinetic rate constants. However, initial sampling of the potential energy surface has a significant impact on the performance of WE, i.e., convergence and efficiency. We therefore introduce deep-learned kinetic modeling approaches that extract statistically relevant information from short MD trajectories to provide a well-sampled initial state distribution for WE simulations. This hybrid approach overcomes any statistical bias to the system, as it runs short unbiased MD trajectories and identifies meaningful metastable states of the system. It is shown to provide a more refined free energy landscape closer to the steady state that could efficiently sample kinetic properties such as rate constants.

Translocation of a single Arg[Formula: see text] peptide across a DOPC/DOPG(4:1) model membrane using the weighted ensemble method

It is difficult to observe a spontaneous translocation of cell-penetrating peptides(CPPs) within a short time scale (e.g., a few hundred ns) in all-atom molecular dynamics(MD) simulations because the time required for the translocation of usual CPPs is on the order of minutes or so. In this work, we report a spontaneous translocation of a single Arg[Formula: see text](R9) across a DOPC/DOPG(4:1) model membrane within an order of a few tens ns scale by using the weighted ensemble(WE) method. We identify how water molecules and the orientation of Arg[Formula: see text] play a role in translocation. We also show how lipid molecules are transported along with Arg[Formula: see text]. In addition, we present free energy profiles of the translocation across the membrane using umbrella sampling and show that a single Arg[Formula: see text] translocation is energetically unfavorable. We expect that the WE method can help study interactions of CPPs with various model membranes within MD simulation approaches.

Weighted ensemble: Recent mathematical developments

Weighted ensemble (WE) is an enhanced sampling method based on periodically replicating and pruning trajectories generated in parallel. WE has grown increasingly popular for computational biochemistry problems due, in part, to improved hardware and accessible software implementations. Algorithmic and analytical improvements have played an important role, and progress has accelerated in recent years. Here, we discuss and elaborate on the WE method from a mathematical perspective, highlighting recent results that enhance the computational efficiency. The mathematical theory reveals a new strategy for optimizing trajectory management that approaches the best possible variance while generalizing to systems of arbitrary dimension.

A Suite of Tutorials for the WESTPA 2.0 Rare-Events Sampling Software [Article v2.0]

The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present two sets of tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. The first set of more basic tutorials describes a range of simulation types, from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding. The second set ecompasses six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for larger systems and/or slower processes. The advanced tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of "binless" schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling for systems biology models. Applications of the advanced tutorials include atomistic and non-spatial models, and consist of complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics or systems biology simulations.

#COVIDisAirborne: AI-enabled multiscale computational microscopy of delta SARS-CoV-2 in a respiratory aerosol

We seek to completely revise current models of airborne transmission of respiratory viruses by providing never-before-seen atomic-level views of the SARS-CoV-2 virus within a respiratory aerosol. Our work dramatically extends the capabilities of multiscale computational microscopy to address the significant gaps that exist in current experimental methods, which are limited in their ability to interrogate aerosols at the atomic/molecular level and thus obscure our understanding of airborne transmission. We demonstrate how our integrated data-driven platform provides a new way of exploring the composition, structure, and dynamics of aerosols and aerosolized viruses, while driving simulation method development along several important axes. We present a series of initial scientific discoveries for the SARS-CoV-2 Delta variant, noting that the full scientific impact of this work has yet to be realized.

Iterative Landmark-Based Umbrella Sampling (ILBUS) Protocol for Sampling of Conformational Space of Biomolecules

Computer simulations of biomolecules such as molecular dynamics often suffer from insufficient sampling. Due to limited computational resources, insufficient sampling prevents obtaining proper equilibrium distributions of observed properties. To deal with this problem, we proposed a simulation protocol for efficient resampling of collected off-equilibrium trajectories. These trajectories are utilized for the initial mapping of the conformational space, which is later properly resampled by the introduced Iterative Landmark-Based Umbrella Sampling (ILBUS) method. Reconstruction of static equilibrium properties is achieved by the multistate Bennett acceptance ratio (MBAR) method, which enables efficient use of simulated data. The ILBUS protocol is geometry-based and does not demand any additional collective variable or a dimensional-reduction technique. The only requirement is a set of suitably spaced reference conformations, which serve as landmarks in the mapped conformational space. Additionally, the ILBUS protocol encompasses an iterative process that optimizes the force constant used in the umbrella sampling simulation. Such tuning is an inherent feature of the protocol and does not need to be performed by the user in advance. Furthermore, even the simulations with suboptimal force constants can be used in estimates by MBAR. We demonstrate the feasibility and the performance of this approach in the study of the conformational landscape of the alanine dipeptide, met-enkephalin, and adenylate kinase.

Multiagent Reinforcement Learning-Based Adaptive Sampling for Conformational Dynamics of Proteins

Machine learning is increasingly applied to improve the efficiency and accuracy of molecular dynamics (MD) simulations. Although the growth of distributed computer clusters has allowed researchers to obtain higher amounts of data, unbiased MD simulations have difficulty sampling rare states, even under massively parallel adaptive sampling schemes. To address this issue, several algorithms inspired by reinforcement learning (RL) have arisen to promote exploration of the slow collective variables (CVs) of complex systems. Nonetheless, most of these algorithms are not well-suited to leverage the information gained by simultaneously sampling a system from different initial states (e.g., a protein in different conformations associated with distinct functional states). To fill this gap, we propose two algorithms inspired by multiagent RL that extend the functionality of closely related techniques (REAP and TSLC) to situations where the sampling can be accelerated by learning from different regions of the energy landscape through coordinated agents. Essentially, the algorithms work by remembering which agent discovered each conformation and sharing this information with others at the action-space discretization step. A stakes function is introduced to modulate how different agents sense rewards from discovered states of the system. The consequences are three-fold: (i) agents learn to prioritize CVs using only relevant data, (ii) redundant exploration is reduced, and (iii) agents that obtain higher stakes are assigned more actions. We compare our algorithm with other adaptive sampling techniques (least counts, REAP, TSLC, and AdaptiveBandit) to show and rationalize the gain in performance.

Mechanistic Insights into Passive Membrane Permeability of Drug-like Molecules from a Weighted Ensemble of Trajectories

Passive permeability of a drug-like molecule is a critical property assayed early in a drug discovery campaign that informs a medicinal chemist how well a compound can traverse biological membranes, such as gastrointestinal epithelial or restrictive organ barriers, so it can perform a specific therapeutic function. However, the challenge that remains is the development of a method, experimental or computational, which can both determine the permeation rate and provide mechanistic insights into the transport process to help with the rational design of any given molecule. Typically, one of the following three methods are used to measure the membrane permeability: (1) experimental permeation assays acting on either artificial or natural membranes; (2) quantitative structure–permeability relationship models that rely on experimental values of permeability or related pharmacokinetic properties of a range of molecules to infer those for new molecules; and (3) estimation of permeability from the Smoluchowski equation, where free energy and diffusion profiles along the membrane normal are taken as input from large-scale molecular dynamics simulations. While all these methods provide estimates of permeation coefficients, they provide very little information for guiding rational drug design. In this study, we employ a highly parallelizable weighted ensemble (WE) path sampling strategy, empowered by cloud computing techniques, to generate unbiased permeation pathways and permeability coefficients for a set of drug-like molecules across a neat 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine membrane bilayer. Our WE method predicts permeability coefficients that compare well to experimental values from an MDCK-LE cell line and PAMPA assays for a set of drug-like amines of varying size, shape, and flexibility. Our method also yields a series of continuous permeation pathways weighted and ranked by their associated probabilities. Taken together, the ensemble of reactive permeation pathways, along with the estimate of the permeability coefficient, provides a clearer picture of the microscopic underpinnings of small-molecule membrane permeation.

Implications of the unfolded state in the folding energetics of heterogeneous-backbone protein mimetics

Sequence-encoded folding is the foundation of protein structure and is also possible in synthetic chains of artificial chemical composition. In natural proteins, the characteristics of the unfolded state are as important as those of the folded state in determining folding energetics. While much is known about folded structures adopted by artificial protein-like chains, corresponding information about the unfolded states of these molecules is lacking. Here, we report the consequences of altered backbone composition on the structure, stability, and dynamics of the folded and unfolded states of a compact helix-rich protein. Characterization through a combination of biophysical experiments and atomistic simulation reveals effects of backbone modification that depend on both the type of artificial monomers employed and where they are applied in sequence. In general, introducing artificial connectivity in a way that reinforces characteristics of the unfolded state ensemble of the prototype natural protein minimizes the impact of chemical changes on folded stability. These findings have implications in the design of protein mimetics and provide an atomically detailed picture of the unfolded state of a natural protein and artificial analogues under non-denaturing conditions.

Biophysical experiments and atomistic simulation reveal impacts of protein backbone alteration on the ensemble that defines the unfolded state. These effects have implications on folded stability of protein mimetics.