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

Detection of Putative Ligand Dissociation Pathways in Proteins Using Site-Identification by Ligand Competitive Saturation

Drug efficacy often correlates better with dissociation kinetics than binding affinity alone. To study binding kinetics computationally, it is necessary to identify all of the possible ligand dissociation pathways. The site identification by ligand competitive saturation (SILCS) method involves the precomputation of a set of maps (FragMaps), which describe the free energy landscapes of typical chemical functionalities in and around a target protein or RNA. In the current work, we present and implement a method to use SILCS to identify ligand dissociation pathways, termed "SILCS-Pathway." The A* pathfinding algorithm is utilized to enumerate ligand dissociation pathways between the ligand binding site and the surrounding bulk solvent environment defined on evenly spaced points around the protein based on a Fibonacci lattice. The cost function for the A* algorithm is calculated using the SILCS exclusion maps and the SILCS grid free energy scores, thereby identifying paths that account for local protein flexibility and potential favorable interactions with the ligand. By traversing all evenly distributed bulk solvent points around the protein, we located all possible dissociation pathways and clustered them to identify general ligand unbinding pathways. The procedure is verified by using proteins studied previously with enhanced sampling molecular dynamics (MD) techniques and is shown to be capable of capturing important ligand dissociation routes in a highly computationally efficient manner. The identified pathways will act as the foundation for determining ligand dissociation kinetics using SILCS free energy profiles, which will be described in a subsequent article.

Thermodynamic Interpolation: A Generative Approach to Molecular Thermodynamics and Kinetics

Using normalizing flows and reweighting, Boltzmann generators enable equilibrium sampling from a Boltzmann distribution, defined by an energy function and thermodynamic state. In this work, we introduce thermodynamic interpolation (TI), which allows for generating sampling statistics in a temperature-controllable way. We introduce TI flavors that work directly in the ambient configurational space, mapping between different thermodynamic states or through a latent, normally distributed reference state. Our ambient-space approach allows for the specification of arbitrary target temperatures, ensuring generalizability within the temperature range of the training set and demonstrating the potential for extrapolation beyond it. We validate the effectiveness of TI on model systems that exhibit metastability and nontrivial temperature dependencies. Finally, we demonstrate how to combine TI-based sampling to estimate free energy differences through various free energy perturbation methods and provide corresponding approximated kinetic rates, estimated through generator extended dynamic mode decomposition (gEDMD).

Molecular simulations reveal intricate coupling between agonist-bound β-adrenergic receptors and G protein

G protein-coupled receptors (GPCRs) and G proteins transmit signals from hormones and neurotransmitters across cell membranes, initiating downstream signaling and modulating cellular behavior. Using advanced computer modeling and simulation, we identified atomistic-level structural, dynamic, and energetic mechanisms of norepinephrine (NE) and stimulatory G protein (Gs) interactions with β-adrenergic receptors (βARs), crucial GPCRs for heart function regulation and major drug targets. Our analysis revealed distinct binding behaviors of NE within β1AR and β2AR despite identical orthosteric binding pockets. β2AR had an additional binding site, explaining variations in NE binding affinities. Simulations showed significant differences in NE dissociation pathways and receptor interactions with the Gs. β1AR binds Gs more strongly, while β2AR induces greater conformational changes in the α subunit of Gs. Furthermore, GTP and GDP binding to Gs may disrupt coupling between NE and βAR, as well as between βAR and Gs. These findings may aid in designing precise βAR-targeted drugs.

Multiple Chaperone DnaK-FliC Flagellin Interactions are Required for Pseudomonas aeruginosa Flagellum Assembly and Indicate a New Function for DnaK

The DnaK (Hsp70) protein is an essential ATP-dependent chaperone foldase and holdase found in most organisms. In this study, combining multiple experimental approaches we determined FliC as major interaction partner of DnaK in the opportunistic bacterial pathogen Pseudomonas aeruginosa. Implementing immunofluorescence microscopy and electron microscopy techniques DnaK was found extracellularly associated to the assembled filament in a regular pattern. dnaK repression led to intracellular FliC accumulation and motility impairment, highlighting DnaK essentiality for FliC export and flagellum assembly. SPOT-membrane peptide arrays coupled with artificial intelligence analyses suggested a highly dynamic DnaK-FliC interaction landscape involving multiple domains and transient complexes formation. Remarkably, in vitro fast relaxation imaging (FReI) experiments mimicking ATP-deprived extracellular environment conditions exhibited DnaK ATP-independent holdase activity, regardless of its co-chaperone DnaJ and its nucleotide exchange factor GrpE. We present a model for the DnaK-FliC interactions involving dynamic states throughout the flagellum assembly stages. These results expand the classical view of DnaK chaperone functioning and introduce a new participant in the Pseudomonas flagellar system, an important trait for bacterial colonisation and virulence.

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

The inherent antigen cross-reactivity of the T cell receptor (TCR) is balanced by high specificity. Surprisingly, TCR specificity often manifests in ways not easily interpreted from static structures. Here we show that TCR discrimination between an HLA-A*03:01 (HLA-A3)-restricted public neoantigen and its wild-type (WT) counterpart emerges from distinct motions within the HLA-A3 peptide binding groove that vary with the identity of the peptide's first primary anchor. These motions create a dynamic gate that, in the presence of the WT peptide, impedes a large conformational change required for TCR binding. The neoantigen is insusceptible to this limiting dynamic, and, with the gate open, upon TCR binding the central tryptophan can transit underneath the peptide backbone to the opposing side of the HLA-A3 peptide binding groove. Our findings thus reveal a novel mechanism driving TCR specificity for a cancer neoantigen that is rooted in the dynamic and allosteric nature of peptide/MHC-I binding grooves, with implications for resolving long-standing and often confounding questions about T cell specificity.

Insights into Translocation of Arginine-Rich Cell-Penetrating Peptides across a Model Membrane

It is well-known that membrane deformation and water pores contribute to the spontaneous translocation of arginine-rich cell-penetrating peptides (CPPs). We confirm this through the observation of the spontaneous translocation of single R9 (nona-arginine) and Tat (48-60) peptides across a model membrane using the weighted ensemble (WE) method within all-atom molecular dynamics (MD) simulations. Furthermore, we demonstrate that membrane deformation and the presence of a water pore reduce the effective charge of the CPP and the bending rigidity of the model membrane during translocation. We find that R9 disturbs the model membrane more than Tat (48-60), leading to more efficient translocation of R9 across the model membrane.

System size effects on the free energy landscapes from molecular dynamics of phase-separating bilayers

The "lipid raft" hypothesis proposes that cell membranes contain distinct domains of varying lipid compositions, where "rafts" of ordered lipids and cholesterol coexist with disordered lipid regions. Experimental and theoretical phase diagrams of model membranes have revealed multiple coexisting phases. Molecular dynamics (MD) simulations can also capture spontaneous phase separation of bilayers. However, these methods merely determine the sign of the free energy change upon phase separation-whether or not it is favorable-but not the amplitude. Recently, we developed a workflow to compute the free energy of phase separation from MD simulations using the weighted ensemble method. However, while theoretical treatments generally focus on infinite systems and experimental measurements on mesoscopic to macroscopic systems, MD simulations are comparatively small. Therefore, if we are to put the results of these calculations into the appropriate context, we need to understand the effects the finite size of the simulation has on the computed free energy landscapes. In this study, we investigate this phenomenon by computing free energy profiles for a model phase-separating system as a function of system size, ranging from 324 to 10 110 lipids. The results suggest that, within the limits of statistical uncertainty, bulk-like behavior emerges once the systems contain roughly 4000 lipids.

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.

Investigation on the molecular, electronic and spectroscopic properties of rosmarinic acid: an intuition from an experimental and computational perspective

Various drugs such as corticosteroids, salbutamol, and β2 agonist are available for the treatment of asthma an inflammatory disease and its symptoms, although the ingredient and the mode of action of these drugs are not clearly elucidated. Hence this research aimed at carrying out improved scientific research with respect to the use of natural product rosmarinic acid which poses minima, side effects. Herein, we first carried out extraction, isolation, and spectroscopic (FT-IR, 1H-NMR and 13C-NMR) investigation, followed by molecular modeling analysis on the naturally occurring rosmarinic acid extracted from Rosmarinus officinalis. A detailed comparison of the experimental and theoretical vibrational analysis has been carried out using five DFT functionals: BHANDH, HSEH1PBE, M06-2X, MPW3PBE and THCTHHYB with the basis set 6-311++G (d, p) to investigate into the structural, reactivity, and stability of the isolated compound. Frontier molecular orbital analysis and appropriate quantum descriptors were calculated. Results showed that the compound was more stable at M06-2X and more reactive at HSEH1PBE with an energy gap of 6.43441 eV and 3.8047 eV, respectively, which was later affirmed by the global quantum reactivity parameters. From natural bond orbital analysis, π* → π* is the major contributor to electron transition with the summation perturbation energy of 889.57 kcal/mol, while π → π* had the perturbation energy totaling of 145.3 kcal/mol. Geometry analysis shows BHANDH to have lower bond length values and lesser deviation from 120° in carbon-carbon angle. The potency of the title molecule as an asthma drug was tested via a molecular docking approach and the binding score of -8.2 kcal/mol was observed against -7.0 of salbutamol standard drug, suggesting romarinic acid as a potential natural organic treatment for asthma.Communicated by Ramaswamy H. Sarma.

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.

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.

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.

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

The weighted ensemble (WE) family of methods is one of several statistical mechanics-based path sampling strategies that can provide estimates of key observables (rate constants and pathways) using a fraction of the time required by direct simulation methods such as molecular dynamics or discrete-state stochastic algorithms. WE methods oversee numerous parallel trajectories using intermittent overhead operations at fixed time intervals, enabling facile interoperability with any dynamics engine. Here, we report on the major upgrades to the WESTPA software package, an open-source, high-performance framework that implements both basic and recently developed WE methods. These upgrades offer substantial improvements over traditional WE methods. The key features of the new WESTPA 2.0 software enhance the efficiency and ease of use: an adaptive binning scheme for more efficient surmounting of large free energy barriers, streamlined handling of large simulation data sets, exponentially improved analysis of kinetics, and developer-friendly tools for creating new WE methods, including a Python API and resampler module for implementing both binned and "binless" WE strategies.

Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling

Gaussian-accelerated molecular dynamics (GaMD) is a well-established enhanced sampling method for molecular dynamics simulations that effectively samples the potential energy landscape of the system by adding a boost potential, which smoothens the surface and lowers the energy barriers between states. GaMD is unable to give time-dependent properties such as kinetics directly. On the other hand, the weighted ensemble (WE) method can efficiently sample transitions between states with its many weighted trajectories, which directly yield rates and pathways. However, convergence to equilibrium conditions remains a challenge for the WE method. Hence, we have developed a hybrid method that combines the two methods, wherein GaMD is first used to sample the potential energy landscape of the system and WE is subsequently used to further sample the potential energy landscape and kinetic properties of interest. We show that the hybrid method can sample both thermodynamic and kinetic properties more accurately and quickly compared to using either method alone.

Free Energy Analyses of Cell-Penetrating Peptides Using the Weighted Ensemble Method

Cell-penetrating peptides (CPPs) have been widely used for drug-delivery agents; however, it has not been fully understood how they translocate across cell membranes. The Weighted Ensemble (WE) method, one of the most powerful and flexible path sampling techniques, can be helpful to reveal translocation paths and free energy barriers along those paths. Within the WE approach we show how Arg9 (nona-arginine) and Tat interact with a DOPC/DOPG(4:1) model membrane, and we present free energy (or potential mean of forces, PMFs) profiles of penetration, although a translocation across the membrane has not been observed in the current simulations. Two different compositions of lipid molecules were also tried and compared. Our approach can be applied to any CPPs interacting with various model membranes, and it will provide useful information regarding the transport mechanisms of CPPs.

A Possible Mechanism of Graphene Oxide to Enhance Thermostability of D-Psicose 3-Epimerase Revealed by Molecular Dynamics Simulations

Thermal stability is a limiting factor for effective application of D-psicose 3-epimerase (DPEase) enzyme. Recently, it was reported that the thermal stability of DPEase was improved by immobilizing enzymes on graphene oxide (GO) nanoparticles. However, the detailed mechanism is not known. In this study, we investigated interaction details between GO and DPEase by performing molecular dynamics (MD) simulations. The results indicated that the domain (K248 to D268) of DPEase was an important anchor for immobilizing DPEase on GO surface. Moreover, the strong interactions between DPEase and GO can prevent loop α1′-α1 and β4-α4 of DPEase from the drastic fluctuation. Since these two loops contained active site residues, the geometry of the active pocket of the enzyme remained stable at high temperature after the DPEase was immobilized by GO, which facilitated efficient catalytic activity of the enzyme. Our research provided a detailed mechanism for the interaction between GO and DPEase at the nano–biology interface.

A glycan gate controls opening of the SARS-CoV-2 spike protein

SARS-CoV-2 infection is controlled by the opening of the spike protein receptor binding domain (RBD), which transitions from a glycan-shielded 'down' to an exposed 'up' state to bind the human angiotensin-converting enzyme 2 receptor and infect cells. While snapshots of the 'up' and 'down' states have been obtained by cryo-electron microscopy and cryo-electron tomagraphy, details of the RBD-opening transition evade experimental characterization. Here over 130 µs of weighted ensemble simulations of the fully glycosylated spike ectodomain allow us to characterize more than 300 continuous, kinetically unbiased RBD-opening pathways. Together with ManifoldEM analysis of cryo-electron microscopy data and biolayer interferometry experiments, we reveal a gating role for the N-glycan at position N343, which facilitates RBD opening. Residues D405, R408 and D427 also participate. The atomic-level characterization of the glycosylated spike activation mechanism provided herein represents a landmark study for ensemble pathway simulations and offers a foundation for understanding the fundamental mechanisms of SARS-CoV-2 viral entry and infection.

A Minimal, Adaptive Binning Scheme for Weighted Ensemble Simulations

A promising approach for simulating rare events with rigorous kinetics is the weighted ensemble path sampling strategy. One challenge of this strategy is the division of configurational space into bins for sampling. Here we present a minimal adaptive binning (MAB) scheme for the automated, adaptive placement of bins along a progress coordinate within the framework of the weighted ensemble strategy. Results reveal that the MAB binning scheme, despite its simplicity, is more efficient than a manual, fixed binning scheme in generating transitions over large free energy barriers, generating a diversity of pathways, estimating rate constants, and sampling conformations. The scheme is general and extensible to any rare-events sampling strategy that employs progress coordinates.

Markov State Models to Elucidate Ligand Binding Mechanism

Molecular dynamics simulations can now routinely access the microsecond timescale, making feasible direct sampling of ligand association events. While Markov State Model (MSM) approaches offer a useful framework for analyzing such trajectory data to gain insight into binding mechanisms, accurate modeling of ligand association pathways and kinetics must be done carefully. We describe methods and good practices for constructing MSMs of ligand binding from unbiased trajectory data and discuss how to use time-lagged independent component analysis (tICA) to build informative models, using as an example recent simulation work to model the binding of phenylalanine to the regulatory ACT domain dimer of phenylalanine hydroxylase. We describe a variety of methods for estimating association rates from MSMs and discuss how to distinguish between conformational selection and induced-fit mechanisms using MSMs. In addition, we review some examples of MSMs constructed to elucidate the mechanisms by which p53 transactivation domain (TAD) and related peptides bind the oncoprotein MDM2.

Ranking of Ligand Binding Kinetics Using a Weighted Ensemble Approach and Comparison with a Multiscale Milestoning Approach

To improve lead optimization efforts in finding the right ligand, pharmaceutical industries need to know the ligand's binding kinetics, such as binding and unbinding rate constants, which often correlate with the ligand's efficacy in vivo. To predict binding kinetics efficiently, enhanced sampling methods, such as milestoning and the weighted ensemble (WE) method, have been used in molecular dynamics (MD) simulations of these systems. However, a comparison of these enhanced sampling methods in ranking ligands has not been done. Hence, a WE approach called the concurrent adaptive sampling (CAS) algorithm that uses MD simulations was used to rank seven ligands for β-cyclodextrin, a system in which a multiscale milestoning approach called simulation enabled estimation of kinetic rates (SEEKR) was also used, which uses both MD and Brownian dynamics simulations. Overall, the CAS algorithm can successfully rank ligands using the unbinding rate constant k_off values and binding free energy ΔG values, as SEEKR did, with reduced computational cost that is about the same as SEEKR. We compare the CAS algorithm simulations with different parameters and discuss the impact of parameters in ranking ligands and obtaining rate constant and binding free energy estimates. We also discuss similarities and differences and advantages and disadvantages of SEEKR and the CAS algorithm for future use.

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na⁺/Cl^- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na⁺/Cl^- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (k_on) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Recent progress in molecular simulation methods for drug binding kinetics

Due to the contribution of drug-target binding kinetics to drug efficacy, there is a high level of interest in developing methods to predict drug-target binding kinetic parameters. During the review period, a wide range of enhanced sampling molecular dynamics simulation-based methods has been developed for computing drug-target binding kinetics and studying binding and unbinding mechanisms. Here, we assess the performance of these methods considering two benchmark systems in detail: mutant T4 lysozyme-ligand complexes and a large set of N-HSP90-inhibitor complexes. The results indicate that some of the simulation methods can already be usefully applied in drug discovery or lead optimization programs but that further studies on more high-quality experimental benchmark datasets are necessary to improve and validate computational methods.