Replica exchange and expanded ensemble simulations as Gibbs sampling: Simple improvements for enhanced mixing

Replica Exchange of Expanded Ensembles: A Generalized Ensemble Approach with Enhanced Flexibility and Parallelizability

Generalized ensemble methods such as Hamiltonian replica exchange (HREX) and expanded ensemble (EE) have been shown effective in free energy calculations for various contexts, given their ability to circumvent free energy barriers via nonphysical pathways defined by states with different modified Hamiltonians. However, both HREX and EE methods come with drawbacks, such as limited flexibility in parameter specification or the lack of parallelizability for more complicated applications. To address this challenge, we present the method of replica exchange of expanded ensembles (REXEE), which integrates the principles of HREX and EE methods by periodically exchanging coordinates of EE replicas sampling different yet overlapping sets of alchemical states. With the solvation free energy calculation of anthracene and binding free energy calculation of the CB7-10 binding complex, we show that the REXEE method achieves the same level of accuracy in free energy calculations as the HREX and EE methods, while offering enhanced flexibility and parallelizability. Additionally, we examined REXEE simulations with various setups to understand how different exchange frequencies and replica configurations influence the sampling efficiency in the fixed-weight phase and the weight convergence in the weight-updating phase. The REXEE approach can be further extended to support asynchronous parallelization schemes, allowing looser communications between larger numbers of loosely coupled processors such as cloud computing and therefore promising much more scalable and adaptive executions of alchemical free energy calculations. All algorithms for the REXEE method are available in the Python package ensemble_md, which offers an interface for REXEE simulation management without modifying the source code in GROMACS.

Advanced computational approaches to understand protein aggregation

Protein aggregation is a widespread phenomenon implicated in debilitating diseases like Alzheimer's, Parkinson's, and cataracts, presenting complex hurdles for the field of molecular biology. In this review, we explore the evolving realm of computational methods and bioinformatics tools that have revolutionized our comprehension of protein aggregation. Beginning with a discussion of the multifaceted challenges associated with understanding this process and emphasizing the critical need for precise predictive tools, we highlight how computational techniques have become indispensable for understanding protein aggregation. We focus on molecular simulations, notably molecular dynamics (MD) simulations, spanning from atomistic to coarse-grained levels, which have emerged as pivotal tools in unraveling the complex dynamics governing protein aggregation in diseases such as cataracts, Alzheimer's, and Parkinson's. MD simulations provide microscopic insights into protein interactions and the subtleties of aggregation pathways, with advanced techniques like replica exchange molecular dynamics, Metadynamics (MetaD), and umbrella sampling enhancing our understanding by probing intricate energy landscapes and transition states. We delve into specific applications of MD simulations, elucidating the chaperone mechanism underlying cataract formation using Markov state modeling and the intricate pathways and interactions driving the toxic aggregate formation in Alzheimer's and Parkinson's disease. Transitioning we highlight how computational techniques, including bioinformatics, sequence analysis, structural data, machine learning algorithms, and artificial intelligence have become indispensable for predicting protein aggregation propensity and locating aggregation-prone regions within protein sequences. Throughout our exploration, we underscore the symbiotic relationship between computational approaches and empirical data, which has paved the way for potential therapeutic strategies against protein aggregation-related diseases. In conclusion, this review offers a comprehensive overview of advanced computational methodologies and bioinformatics tools that have catalyzed breakthroughs in unraveling the molecular basis of protein aggregation, with significant implications for clinical interventions, standing at the intersection of computational biology and experimental research.

Enhancing torsional sampling using fully adaptive simulated tempering

Enhanced sampling algorithms are indispensable when working with highly disconnected multimodal distributions. An important application of these is the conformational exploration of particular internal degrees of freedom of molecular systems. However, despite the existence of many commonly used enhanced sampling algorithms to explore these internal motions, they often rely on system-dependent parameters, which negatively impact efficiency and reproducibility. Here, we present fully adaptive simulated tempering (FAST), a variation of the irreversible simulated tempering algorithm, which continuously optimizes the number, parameters, and weights of intermediate distributions to achieve maximally fast traversal over a space defined by the change in a predefined thermodynamic control variable such as temperature or an alchemical smoothing parameter. This work builds on a number of previously published methods, such as sequential Monte Carlo, and introduces a novel parameter optimization procedure that can, in principle, be used in any expanded ensemble algorithms. This method is validated by being applied on a number of different molecular systems with high torsional kinetic barriers. We also consider two different soft-core potentials during the interpolation procedure and compare their performance. We conclude that FAST is a highly efficient algorithm, which improves simulation reproducibility and can be successfully used in a variety of settings with the same initial hyperparameters.

Calculations of Absolute Free Energies, Enthalpies, and Entropies for Drug Binding

Computer simulation methods can aid in the rational design of drugs aimed at a specific target, typically a protein. The affinity of a drug for its target is given by the free energy of binding. Binding can be further characterized by the enthalpy and entropy changes in the process. Methods exist to determine exact free energies, enthalpies, and entropies that are dependent only on the quality of the potential model and adequate sampling of conformational degrees of freedom. Entropy and enthalpy are roughly an order of magnitude more difficult to calculate than the free energy. This project combines a replica exchange method for enhanced sampling, designed to be efficient for protein-sized systems, with free energy calculations. This approach, replica exchange with dynamical scaling (REDS), uses two conventional simulations at different temperatures so that the entropy can be found from the temperature dependence of the free energy. A third replica is placed between them, with a modified Hamiltonian that allows it to span the temperature range of the conventional replicas. REDS provides temperature-dependent data and aids in sampling. It is applied to the bromodomain-containing protein 4 (BRD4) system. We find that for the force fields used, the free energies are accurate but the entropies and enthalpies are not, with the entropic contribution being too positive. Reproducing the entropy and enthalpy of binding appears to be a more stringent test of the force fields than reproducing the free energy.

Decoding the dynamics of BCL9 triazole stapled peptide

BCL9 is a key protein in Wnt signaling pathway. It acts as a transcriptional co-activator to β-catenin, and dysregulation in this pathway leads to tumor growth. Inhibiting such a protein-protein interaction is considered as a therapeutic challenge. The interaction between β-catenin and BCL9 is facilitated by a 23-residue helical domain from BCL9 and a hydrophobic groove of β-catenin. To prevent this interaction, a peptide that mimics the alpha-helical domain of BCL9 can be designed. Stapling is considered a successful strategy in the pursuit of designing such peptides in which amino acids side are stitched together using chemical moieties. Among the various types of cross-linkers, triazole is the most rapid and effective one synthesized via click reaction. However, the underlying interactions behind maintaining the secondary structure of stapled peptides remain less explored. In the current work, we employed the molecular dynamics simulation to study the conformational behavior of the experimentally synthesized single and double triazole stapled BCL9 peptide. Upon the addition of a triazole staple, there is a significant reduction in the conformational space of BCL9. The helical character of the stapled peptide increases with an increase in separation between the triazole cross-linkers. Also, we encompassed the Replica Exchange with Solute Tempering (REST2) simulation to validate the high-temperature response of the stapled peptide. From REST2, the PCA and t-SNE show the reduction in distinct cluster formation on the addition of triazole staple. Our study infers further development of these triazole-stapled BCL9 peptides into effective inhibitors to target the interaction between β-catenin and BCL9.

Kartograf: A Geometrically Accurate Atom Mapper for Hybrid-Topology Relative Free Energy Calculations

Relative binding free energy (RBFE) calculations have emerged as a powerful tool that supports ligand optimization in drug discovery. Despite many successes, the use of RBFEs can often be limited by automation problems, in particular, the setup of such calculations. Atom mapping algorithms are an essential component in setting up automatic large-scale hybrid-topology RBFE calculation campaigns. Traditional algorithms typically employ a 2D subgraph isomorphism solver (SIS) in order to estimate the maximum common substructure. SIS-based approaches can be limited by time-intensive operations and issues with capturing geometry-linked chemical properties, potentially leading to suboptimal solutions. To overcome these limitations, we have developed Kartograf, a geometric-graph-based algorithm that uses primarily the 3D coordinates of atoms to find a mapping between two ligands. In free energy approaches, the ligand conformations are usually derived from docking or other previous modeling approaches, giving the coordinates a certain importance. By considering the spatial relationships between atoms related to the molecule coordinates, our algorithm bypasses the computationally complex subgraph matching of SIS-based approaches and reduces the problem to a much simpler bipartite graph matching problem. Moreover, Kartograf effectively circumvents typical mapping issues induced by molecule symmetry and stereoisomerism, making it a more robust approach for atom mapping from a geometric perspective. To validate our method, we calculated mappings with our novel approach using a diverse set of small molecules and used the mappings in relative hydration and binding free energy calculations. The comparison with two SIS-based algorithms showed that Kartograf offers a fast alternative approach. The code for Kartograf is freely available on GitHub (https://github.com/OpenFreeEnergy/kartograf). While developed for the OpenFE ecosystem, Kartograf can also be utilized as a standalone Python package.

Fast free energy estimates from λ-dynamics with bias-updated Gibbs sampling

Relative binding free energy calculations have become an integral computational tool for lead optimization in structure-based drug design. Classical alchemical methods, including free energy perturbation or thermodynamic integration, compute relative free energy differences by transforming one molecule into another. However, these methods have high operational costs due to the need to perform many pairwise perturbations independently. To reduce costs and accelerate molecular design workflows, we present a method called λ-dynamics with bias-updated Gibbs sampling. This method uses dynamic biases to continuously sample between multiple ligand analogues collectively within a single simulation. We show that many relative binding free energies can be determined quickly with this approach without compromising accuracy. For five benchmark systems, agreement to experiment is high, with root mean square errors near or below 1.0 kcal mol-1. Free energy results are consistent with other computational approaches and within statistical noise of both methods (0.4 kcal mol-1 or less). Notably, large efficiency gains over thermodynamic integration of 18-66-fold for small perturbations and 100-200-fold for whole aromatic ring substitutions are observed. The rapid determination of relative binding free energies will enable larger chemical spaces to be more readily explored and structure-based drug design to be accelerated.

Identifying and Overcoming the Sampling Challenges in Relative Binding Free Energy Calculations of a Model Protein:Protein Complex

Relative alchemical binding free energy calculations are routinely used in drug discovery projects to optimize the affinity of small molecules for their drug targets. Alchemical methods can also be used to estimate the impact of amino acid mutations on protein:protein binding affinities, but these calculations can involve sampling challenges due to the complex networks of protein and water interactions frequently present in protein:protein interfaces. We investigate these challenges by extending a graphics processing unit (GPU)-accelerated open-source relative free energy calculation package (Perses) to predict the impact of amino acid mutations on protein:protein binding. Using the well-characterized model system barnase:barstar, we describe analyses for identifying and characterizing sampling problems in protein:protein relative free energy calculations. We find that mutations with sampling problems often involve charge-changes, and inadequate sampling can be attributed to slow degrees of freedom that are mutation-specific. We also explore the accuracy and efficiency of current state-of-the-art approaches─alchemical replica exchange and alchemical replica exchange with solute tempering─for overcoming relevant sampling problems. By employing sufficiently long simulations, we achieve accurate predictions (RMSE 1.61, 95% CI: [1.12, 2.11] kcal/mol), with 86% of estimates within 1 kcal/mol of the experimentally determined relative binding free energies and 100% of predictions correctly classifying the sign of the changes in binding free energies. Ultimately, we provide a model workflow for applying protein mutation free energy calculations to protein:protein complexes, and importantly, catalog the sampling challenges associated with these types of alchemical transformations. Our free open-source package (Perses) is based on OpenMM and is available at https://github.com/choderalab/perses.

Alchemical Free Energy and Hamiltonian Replica Exchange Molecular Dynamics to Compute Hydrofluorocarbon Isotherms in Imidazolium-Based Ionic Liquids

Ionic liquids (ILs) have shown promise for applications that leverage differential gas solubility in an IL solvent, e.g., gas separations. Although most available literature provides Henry's law constants, the ability to efficiently estimate full isotherms is important for engineering design calculations. Molecular simulation can be used as a tool to predict full isotherms of gas in ILs. However, particle insertions or deletions in a charge-dense IL medium and the sluggish conformational dynamics of ILs present two sampling challenges for these systems. We therefore devised a method that uses Hamiltonian replica exchange (HREX) molecular dynamics (MD) combined with alchemical free energy calculations to compute full solubility isotherms of two different hydrofluorocarbons (HFCs) in imidazolium-based IL binary mixtures. This workflow is significantly faster than the Gibbs ensemble Monte Carlo (GEMC) simulations which fail to deal with the slow conformational relaxation caused by the sluggish dynamics of ILs. Multiple free energy estimators, including thermodynamic integration, free energy perturbation, and multistate Bennett acceptance ratio method, provided consistent results. Overall, the simulated Henry's law constant, isotherm curvature, and solubility trends match experimental results reasonably well. We close by calculating the full solubility isotherms of two HFCs in IL mixtures that have not been reported in the literature, demonstrating the potential of this method to be used for solubility prediction and setting the stage for future computational screening studies that search for the "best" IL to separate azeotropic HFC mixtures.

Identifying and overcoming the sampling challenges in relative binding free energy calculations of a model protein:protein complex

Relative alchemical binding free energy calculations are routinely used in drug discovery projects to optimize the affinity of small molecules for their drug targets. Alchemical methods can also be used to estimate the impact of amino acid mutations on protein:protein binding affinities, but these calculations can involve sampling challenges due to the complex networks of protein and water interactions frequently present in protein:protein interfaces. We investigate these challenges by extending a GPU-accelerated opensource relative free energy calculation package (Perses) to predict the impact of amino acid mutations on protein:protein binding. Using the well-characterized model system barnase:barstar, we describe analyses for identifying and characterizing sampling problems in protein:protein relative free energy calculations. We find that mutations with sampling problems often involve charge-changes, and inadequate sampling can be attributed to slow degrees of freedom that are mutation-specific. We also explore the accuracy and efficiency of current state-of-the-art approaches-alchemical replica exchange and alchemical replica exchange with solute tempering-for overcoming relevant sampling problems. By employing sufficiently long simulations, we achieve accurate predictions (RMSE 1.61, 95% CI: [1.12, 2.11] kcal/mol), with 86% of estimates within 1 kcal/mol of the experimentally-determined relative binding free energies and 100% of predictions correctly classifying the sign of the changes in binding free energies. Ultimately, we provide a model workflow for applying protein mutation free energy calculations to protein:protein complexes, and importantly, catalog the sampling challenges associated with these types of alchemical transformations. Our free open-source package (Perses) is based on OpenMM and available at https://github.com/choderalab/perses .

Q-RepEx: A Python pipeline to increase the sampling of empirical valence bond simulations

The exploration of chemical systems occurs on complex energy landscapes. Comprehensively sampling rugged energy landscapes with many local minima is a common problem for molecular dynamics simulations. These multiple local minima trap the dynamic system, preventing efficient sampling. This is a particular challenge for large biochemical systems with many degrees of freedom. Replica exchange molecular dynamics (REMD) is an approach that accelerates the exploration of the conformational space of a system, and thus can be used to enhance the sampling of complex biomolecular processes. In parallel, the empirical valence bond (EVB) approach is a powerful approach for modeling chemical reactivity in biomolecular systems. Here, we present an open-source Python-based tool that interfaces with the Q simulation package, and increases the sampling efficiency of the EVB free energy perturbation/umbrella sampling approach by means of REMD. This approach, Q-RepEx, both decreases the computational cost of the associated REMD-EVB simulations, and opens the door to more efficient studies of biochemical reactivity in systems with significant conformational fluctuations along the chemical reaction coordinate.

The vitamin D receptor as a potential target for the toxic effects of per- and polyfluoroalkyl substances (PFASs): An in-silico study

Due to their persistence and toxicity, perfluoroalkyl and polyfluoroalkyl substances (PFASs) constitute significant hazards to human health and the environment. Their effects include immune suppression, altered hormone levels, and osteoporosis. Recently, the most studied PFAS, perfluorooctanoic acid (PFOA), was shown to competitively binding to the Vitamin D receptor (VDR). VDR plays a crucial role in regulating genes involved in maintaining immune, endocrine, and calcium homeostasis, suggesting it may be a target for at least some of the health effects of PFAS. Hence, this study examined the potential binding of 5206 PFASs to VDR using molecular docking, molecular dynamics, and free energy binding calculations. We identified 14 PFAS that are predicted to interact strongly with VDR, similar to the natural ligands. We further investigated the interactions of VDR with 256 PFASs of established commercial importance. Eighty-three (32%) of these 256 commercially important PFAS were predicted to be stronger binders to VDR than PFOA. At least 16 PFASs of regulatory importance, because they have been identified in water supplies and human blood samples, were also more potent binders to VDR than PFOA. Further, PFASs are usually found together in contaminated drinking water and human blood samples, which raises the concern that multiple PFASs may act together as a mixture on VDR function, potentially producing harmful effects on the immune, endocrine, and bone homeostasis.

Inverse Design of Pore Wall Chemistry To Control Solute Transport and Selectivity

Next-generation membranes for purification and reuse of highly contaminated water require materials with precisely tuned functionality to address key challenges, including the removal of small, charge-neutral solutes. Bioinspired multifunctional membrane surfaces enhance transport properties, but the combinatorically large chemical space is difficult to navigate through trial and error. Here, we demonstrate a computational inverse design approach to efficiently identify promising materials and elucidate design rules. We develop a combined evolutionary optimization, machine learning, and molecular simulation workflow to spatially design chemical functional group patterning in a model nanopore that enhances transport of water relative to solutes. The genetic optimization discovers nonintuitive functionalization strategies that hinder the transport of solutes through the pore, simply by patterning hydrophobic methyl and hydrophilic hydroxyl functional groups. Examining these patterns, we demonstrate that they exploit an unexpected diffusive solute hopping mechanism. This inverse design procedure and the identification of novel molecular mechanisms for pore chemical heterogeneity to impact solute selectivity demonstrate new routes to the design of membrane materials with novel functionalities. More broadly, this work illustrates how chemical design is a powerful strategy to modulate water-mediated surface-solute interactions in complex, soft material systems that are relevant to diverse technologies.

Nonbonded Force Field Parameters from Minimal Basis Iterative Stockholder Partitioning of the Molecular Electron Density Improve CB7 Host-Guest Affinity Predictions

Binding affinity prediction by means of computer simulation has been increasingly incorporated in drug discovery projects. Its wide application, however, is limited by the prediction accuracy of the free energy calculations. The main error sources are force fields used to describe molecular interactions and incomplete sampling of the configurational space. Organic host-guest systems have been used to address force field quality because they share similar interactions found in ligands and receptors, and their rigidity facilitates configurational sampling. Here, we test the binding free energy prediction accuracy for 14 guests with an aromatic or adamantane core and the CB7 host using molecular electron density derived nonbonded force field parameters. We developed a computational workflow written in Python to derive atomic charges and Lennard-Jones parameters with the Minimal Basis Iterative Stockholder method using the polarized electron density of several configurations of each guest in the bound and unbound states. The resulting nonbonded force field parameters improve binding affinity prediction, especially for guests with an adamantane core in which repulsive exchange and dispersion interactions to the host dominate.

From data to noise to data for mixing physics across temperatures with generative artificial intelligence

Using simulations or experiments performed at some set of temperatures to learn about the physics or chemistry at some other arbitrary temperature is a problem of immense practical and theoretical relevance. Here we develop a framework based on statistical mechanics and generative artificial intelligence that allows solving this problem. Specifically, we work with denoising diffusion probabilistic models and show how these models in combination with replica exchange molecular dynamics achieve superior sampling of the biomolecular energy landscape at temperatures that were never simulated without assuming any particular slow degrees of freedom. The key idea is to treat the temperature as a fluctuating random variable and not a control parameter as is usually done. This allows us to directly sample from the joint probability distribution in configuration and temperature space. The results here are demonstrated for a chirally symmetric peptide and single-strand RNA undergoing conformational transitions in all-atom water. We demonstrate how we can discover transition states and metastable states that were previously unseen at the temperature of interest and even bypass the need to perform further simulations for a wide range of temperatures. At the same time, any unphysical states are easily identifiable through very low Boltzmann weights. The procedure while shown here for a class of molecular simulations should be more generally applicable to mixing information across simulations and experiments with varying control parameters.

Comparing the Performance of Different AMBER Protein Forcefields, Partial Charge Assignments, and Water Models for Absolute Binding Free Energy Calculations

Identifying chemical starting points is a vital first step in small molecule drug discovery and can take significant time and money. For this reason, computational approaches to virtual screening are of great interest as they can lower the cost and shorten timeframes. However, simple approaches such as molecular docking and pharmacophore screening are of limited accuracy and provide a low probability of success. Alchemical binding free energies represent a promising approach for virtual screening as they naturally incorporate the key effects of water molecules, protein flexibility, and binding entropy. However, the calculations are technically very challenging, with performance depending on the specific forcefield used. For this reason, it is important that the community has access to benchmark test sets to assess prediction accuracy. In this paper, we present an approach to alchemical binding free energies using OpenMM. We identify effective simulation parameters using an existing BRD4(1) test set and present two new benchmark sets (cMET and PDE2A) that can be used in the community for validation purposes. Our findings also highlight the effectiveness of some AMBER forcefields, in particular, AMBER ff15ipq.

Bayesian-Inference-Driven Model Parametrization and Model Selection for 2CLJQ Fluid Models

A high level of physical detail in a molecular model improves its ability to perform high accuracy simulations but can also significantly affect its complexity and computational cost. In some situations, it is worthwhile to add complexity to a model to capture properties of interest; in others, additional complexity is unnecessary and can make simulations computationally infeasible. In this work, we demonstrate the use of Bayesian inference for molecular model selection, using Monte Carlo sampling techniques accelerated with surrogate modeling to evaluate the Bayes factor evidence for different levels of complexity in the two-centered Lennard-Jones + quadrupole (2CLJQ) fluid model. Examining three nested levels of model complexity, we demonstrate that the use of variable quadrupole and bond length parameters in this model framework is justified only for some chemistries. Through this process, we also get detailed information about the distributions and correlation of parameter values, enabling improved parametrization and parameter analysis. We also show how the choice of parameter priors, which encode previous model knowledge, can have substantial effects on the selection of models, penalizing careless introduction of additional complexity. We detail the computational techniques used in this analysis, providing a roadmap for future applications of molecular model selection via Bayesian inference and surrogate modeling.

Site-Specific Water Dynamics in the First Hydration Layer of an Anti-Freeze Glyco-Protein: A Simulation Study

Antifreeze glycoproteins (AFGPs) inhibit ice recrystallization by a mechanism remaining largely elusive. Dynamics of AFGPs’ hydration water and its involvement in the antifreeze activity, for instance, have not been identified...

Free Energy-Based Computational Methods for the Study of Protein-Peptide Binding Equilibria

This chapter discusses the theory and application of physics-based free energy methods to estimate protein-peptide binding free energies. It presents a statistical mechanics formulation of molecular binding, which is then specialized in three methodologies: (1) alchemical absolute binding free energy estimation with implicit solvation, (2) alchemical relative binding free energy estimation with explicit solvation, and (3) potential of mean force binding free energy estimation. Case studies of protein-peptide binding application taken from the recent literature are discussed for each method.

Using a Coarse-Grained Modeling Framework to Identify Oligomeric Motifs with Tunable Secondary Structure

Coarse-grained modeling can be used to explore general theories that are independent of specific chemical detail. In this paper, we present cg_openmm, a Python-based simulation framework for modeling coarse-grained hetero-oligomers and screening them for structural and thermodynamic characteristics of cooperative secondary structures. cg_openmm facilitates the building of coarse-grained topology and random starting configurations, setup of GPU-accelerated replica exchange molecular dynamics simulations with the OpenMM software package, and features a suite of postprocessing thermodynamic and structural analysis tools. In particular, native contact analysis, heat capacity calculations, and free energy of folding calculations are used to identify and characterize cooperative folding transitions and stable secondary structures. In this work, we demonstrate the capabilities of cg_openmm on a simple 1-1 Lennard-Jones coarse-grained model, in which each residue contains 1 backbone and 1 side-chain bead. By scanning both nonbonded and bonded force-field parameter spaces at the coarse-grained level, we identify and characterize sets of parameters which result in the formation of stable helices through cooperative folding transitions. Moreover, we show that the geometries and stabilities of these helices can be tuned by manipulating the force-field parameters.

Expanded Ensemble Methods Can be Used to Accurately Predict Protein-Ligand Relative Binding Free Energies

Alchemical free energy methods have become indispensable in computational drug discovery for their ability to calculate highly accurate estimates of protein-ligand affinities. Expanded ensemble (EE) methods, which involve single simulations visiting all of the alchemical intermediates, have some key advantages for alchemical free energy calculation. However, there have been relatively few examples published in the literature of using expanded ensemble simulations for free energies of protein-ligand binding. In this paper, as a test of expanded ensemble methods, we compute relative binding free energies using the Open Force Field Initiative force field (codename "Parsley") for 24 pairs of Tyk2 inhibitors derived from a congeneric series of 16 compounds. The EE predictions agree well with the experimental values (root-mean-square error (RMSE) of 0.94 ± 0.13 kcal mol^-1 and mean unsigned error (MUE) of 0.75 ± 0.12 kcal mol^-1). We find that while increasing the number of alchemical intermediates can improve the phase space overlap, faster convergence can be obtained with fewer intermediates, as long as acceptance rates are sufficient. We also find that convergence can be improved using more aggressive updating of biases, and that estimates can be improved by performing multiple independent EE calculations. This work demonstrates that EE is a viable option for alchemical free energy calculation. We discuss the implications of these findings for rational drug design, as well as future directions for improvement.

Molecular Environment-Specific Atomic Charges Improve Binding Affinity Predictions of SAMPL5 Host-Guest Systems

Host-guest systems are widely used in benchmarks as model systems to improve computational methods for absolute binding free energy predictions. Recent advances in sampling algorithms for alchemical free energy calculations and the increase in computational power have made their binding affinity prediction primarily dependent on the quality of the force field. Here, we propose a new methodology to derive the atomic charges of host-guest systems based on quantum mechanics/molecular mechanics calculations and minimal basis iterative stockholder (MBIS) partitioning of the polarized electron density. A newly developed interface between the OpenMM and ORCA software packages provides D-MBIS charges that represent the guest's average electrostatic interactions in the hosts or the solvent. The simulation workflow also calculates the average energy required to polarize the guest in the bound and unbound state. Alchemical free energy calculations using the general Amber force field parameters with D-MBIS charges improve the binding affinity prediction of six guests bound to two octa acid hosts compared to the AM1-BCC charge set after correction with the average energetic polarization cost. This correction originates from the difference in potential energy that is required to polarize the guest in the bound and unbound state and contributes significantly to the binding affinity of anionic guests.

PyAutoFEP: An Automated Free Energy Perturbation Workflow for GROMACS Integrating Enhanced Sampling Methods

Free energy perturbation (FEP) calculations are now routinely used in drug discovery to estimate the relative FEB (RFEB) of small molecules to a biomolecular target of interest. Using enhanced sampling can improve the correlation between predictions and experimental data, especially in systems with conformational changes. Due to the large number of perturbations required in drug discovery campaigns, the manual setup of FEP calculations is no longer viable. Here, we introduce PyAutoFEP, a flexible and open-source tool to aid the setup of RFEB FEP. PyAutoFEP is written in Python3, and automates the generation of perturbation maps, dual topologies, system building and molecular dynamics (MD), and analysis. PyAutoFEP supports multiple force fields, incorporates replica exchange with solute tempering (REST) and replica exchange with solute scaling (REST2) enhanced sampling methods, and allows flexible λ values along perturbation windows. To validate PyAutoFEP, it was applied to a set of 14 Farnesoid X receptor ligands, a system included in the drug design data resource grand challenge 2. An 88% mean correct sign prediction was achieved, and 75% of the predictions had an error below 1.5 kcal/mol. Results using Amber03/GAFF, CHARMM36m/CGenFF, and OPLS-AA/M/LigParGen had Pearson's r values of 0.71 ± 0.13, 0.30 ± 0.27, and 0.66 ± 0.20, respectively. The Amber03/GAFF and OPLS-AA/M/LigParGen results were on par with the top grand challenge 2 submissions. Applying REST2 improved the results using CHARMM36m/CGenFF (Pearson's r = 0.43 ± 0.21) but had little impact on the other force fields. CHARMM36-YF and CHARMM36-WYF modifications did not yield improved predictions compared to CHARMM36m. Finally, we estimated the probability of finding a molecule 1 pK_i better than a lead when using PyAutoFEP to screen 10 or 100 analogues. The probabilities, when compared to random sampling, increased up to sevenfold when 100 molecules were to be screened, suggesting that PyAutoFEP would likely be useful for lead optimization. PyAutoFEP is available on GitHub at https://github.com/lmmpf/PyAutoFEP.

The accelerated weight histogram method for alchemical free energy calculations

The accelerated weight histogram method is an enhanced sampling technique used to explore free energy landscapes by applying an adaptive bias. The method is general and easy to extend. Herein, we show how it can be used to efficiently sample alchemical transformations, commonly used for, e.g., solvation and binding free energy calculations. We present calculations and convergence of the hydration free energy of testosterone, representing drug-like molecules. We also include methane and ethanol to validate the results. The protocol is easy to use, does not require a careful choice of parameters, and scales well to accessible resources, and the results converge at least as quickly as when using conventional methods. One benefit of the method is that it can easily be combined with other reaction coordinates, such as intermolecular distances.

Generalizing the Discrete Gibbs Sampler-Based λ-Dynamics Approach for Multisite Sampling of Many Ligands

In this work, the discrete λ variant of the Gibbs sampler-based λ-dynamics (d-GSλD) method is developed to enable multiple functional group perturbations to be investigated at one or more sites of substitution off a common ligand core. The theoretical framework and special considerations for constructing discrete λ states for multisite d-GSλD are presented. The precision and accuracy of the d-GSλD method is evaluated with three test cases of increasing complexity. Specifically, methyl → methyl symmetric perturbations in water, 1,4-benzene hydration free energies and protein-ligand binding affinities for an example HIV-1 reverse transcriptase inhibitor series are computed with d-GSλD. Complementary MSλD calculations were also performed to compare with d-GSλD's performance. Excellent agreement between d-GSλD and MSλD is observed, with mean unsigned errors of 0.12 and 0.22 kcal/mol for computed hydration and binding free energy test cases, respectively. Good agreement with experiment is also observed, with errors of 0.5-0.7 kcal/mol. These findings support the applicability of the d-GSλD free energy method for a variety of molecular design problems, including structure-based drug design. Finally, a discussion of d-GSλD versus MSλD approaches is presented to compare and contrast features of both methods.

To the fast calculation of the solvation free energy. Combining expanded ensembles with L2MC

A more efficient version of the Expanded Ensembles method for calculation of free energy in molecular-mechanical simulations is proposed. The method is based on the Horowitz L2MC approach to accelerate movement along the "alchemical" coordinate. It is possible to achieve the same efficiency of the algorithm both with the optimal number of "windows" and with a larger number of them compared to the original algorithm. Since the optimal number of windows is unknown a priory, the proposed algorithm is more robust than the traditional one. We can choose the number of windows in excess and do not worry about the loss of efficiency. We illustrate the method's efficiency with the computation of the hydration free energies of pyridine and water.

Conformational Preferences of an Intrinsically Disordered Protein Domain: A Case Study for Modern Force Fields

Molecular simulations of intrinsically disordered proteins (IDPs) are challenging because they require sampling a very large number of relevant conformations, corresponding to a multitude of shallow minima in a flat free energy landscape. However, in the presence of a binding partner, the free energy landscape of an IDP can be dominated by few deep minima. This characteristic imposes high demands on the accuracy of the force field used to describe the molecular interactions. Here, as a model system for an IDP that is unstructured in solution but folds upon binding to a structured interaction partner, the transactivation domain of c-Myb was studied both in the unbound (free) form and when bound to the KIX domain. Six modern biomolecular force fields were systematically tested and compared in terms of their ability to describe the structural ensemble of the IDP. The protein force field/water model combinations included in this study are AMBER ff99SB-disp with its corresponding water model that was derived from TIP4P-D, CHARMM36m with TIP3P, ff15ipq with SPC/E_b, ff99SB*-ILDNP with TIP3P and TIP4P-D, and FB15 with TIP3P-FB water. Comparing the results from REST2-enhanced sampling simulations with experimental CD spectra and secondary chemical shifts reveals that the ff99SB-disp force field can realistically capture the broad and mildly helical structural ensemble of free c-Myb. The structural ensembles yielded by CHARMM36m, ff99SB*-ILDNP together with TIP4P-D water, and FB15 are also mildly helical; however, each of these force fields can be assigned a specific subset of c-Myb residues for which the simulations could not reproduce the experimental secondary chemical shifts. In addition, microsecond-timescale MD simulations of the KIX/c-Myb complex show that most force fields used preserve a stable helix fold of c-Myb in the complex. Still, all force fields predict a KIX/c-Myb complex interface that differs slightly from the structures provided by NMR because several NOE-derived distances between KIX and c-Myb were exceeded in the simulations. Taken together, the ff99SB-disp force field in the first place but also CHARMM36m, ff99SB*-ILDNP together with TIP4P-D water, and FB15 can be suitable choices for future simulation studies of the coupled folding and binding mechanism of the KIX/c-Myb complex and potentially also other IDPs.

Skewed thermodynamic geometry and optimal free energy estimation

Free energy differences are a central quantity of interest in physics, chemistry, and biology. We develop design principles that improve the precision and accuracy of free energy estimators, which have potential applications to screening for targeted drug discovery. Specifically, by exploiting the connection between the work statistics of time-reversed protocol pairs, we develop near-equilibrium approximations for moments of the excess work and analyze the dominant contributions to the precision and accuracy of standard nonequilibrium free-energy estimators. Within linear response, minimum-dissipation protocols follow the geodesics of the Riemannian metric induced by the Stokes friction tensor. We find that the next-order contribution arises from the rank-3 supra-Stokes tensor that skews the geometric structure such that minimum-dissipation protocols follow the geodesics of a generalized cubic Finsler metric. Thus, near equilibrium, the supra-Stokes tensor determines the leading-order contribution to the bias of bidirectional free-energy estimators.

Structural evidence for pheromone discrimination by the pheromone binding protein 3 from Plutella xylostella

Insect pheromone binding proteins (PBPs) are believed to have a high degree of pheromone selectivity, acting as the first filter to discriminate specific pheromones from other volatile compounds. Herein, we provide evidence using homology-based model for the pheromone discrimination of Plutella xylostella pheromone binding protein 3 (PxPBP3). Combining molecular dynamics simulations and in vitro binding assays, two dominant sites are determined to be essential for the PxPBP3 to discriminate (Z)-11-hexadecenyl acetate (Hexadecenyl) from (Z)-11-hexadecenal (Hexadecenal). As the first key site for pheromone discrimination, Arg111 is indispensable to the PxPBP3-Hexadecenyl interaction. However, its importance in the binding of Hexadecenal to PxPBP3 is greatly reduced. A second site where pheromone discrimination occurs is a small loop (residues 34-38) in PxPBP3. It is shown that the hydrophobic strength provided by three hydrophobic residues (Phe34, Tyr37, and Trp38) in the small loop is significantly biased in the two complexes formed by PxPBP3 and the two pheromones. The discrimination capacity of PxPBP3 indicates that the P. xylostella pheromones may not share the same peri-receptor pathway, although they both show high affinity to PxPBP3.

Exploring Protocols to Build Reservoirs to Accelerate Temperature Replica Exchange MD Simulations

Temperature replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method for accelerating biomolecular simulations. During the past 2 decades, several variants of REMD have been developed to further improve the rate of conformational sampling of REMD. One such variant, reservoir REMD (RREMD), was shown to improve the rate of conformational sampling by around 5-20×. Despite the significant increase in the sampling speed, RREMD methods have not been widely used because of the difficulties in building the reservoir and also because of the code not being available on the graphics processing units (GPUs). In this work, we ported the Amber RREMD code onto GPUs making it 20× faster than the central processing unit code. Then, we explored protocols for building Boltzmann-weighted reservoirs as well as non-Boltzmann reservoirs and tested how each choice affects the accuracy of the resulting RREMD simulations. We show that, using the recommended protocols outlined here, RREMD simulations can accurately reproduce Boltzmann-weighted ensembles obtained by much more expensive conventional temperature-based REMD simulations, with at least 15× faster convergence rates even for larger proteins (>50 amino acids) compared to conventional REMD.

Simulated tempering with irreversible Gibbs sampling techniques

We present here two novel algorithms for simulated tempering simulations, which break the detailed balance condition (DBC) but satisfy the skewed detailed balance to ensure invariance of the target distribution. The irreversible methods we present here are based on Gibbs sampling and concern breaking DBC at the update scheme of the temperature swaps. We utilize three systems as a test bed for our methods: a Markov chain Monte Carlo simulation on a simple system described by a one-dimensional double well potential, the Ising model, and molecular dynamics simulations on alanine pentapeptide (ALA5). The relaxation times of inverse temperature, magnetic susceptibility, and energy density for the Ising model indicate clear gains in sampling efficiency over conventional Gibbs sampling techniques with DBC and also over the conventionally used simulated tempering with the Metropolis-Hastings (MH) scheme. Simulations on ALA5 with a large number of temperatures indicate distinct gains in mixing times for inverse temperature and consequently the energy of the system compared to conventional MH. With no additional computational overhead, our methods were found to be more efficient alternatives to the conventionally used simulated tempering methods with DBC. Our algorithms should be particularly advantageous in simulations of large systems with many temperature ladders, as our algorithms showed a more favorable constant scaling in Ising spin systems as compared with both reversible and irreversible MH algorithms. In future applications, our irreversible methods can also be easily tailored to utilize a given dynamical variable other than temperature to flatten rugged free energy landscapes.

Identification of Optimal Ligand Growth Vectors Using an Alchemical Free-Energy Method

In this work, a novel method to rationally design inhibitors with improved steric contacts and enhanced binding free energies is presented. This new method uses alchemical single step perturbation calculations to rapidly optimize the van der Waals interactions of a small molecule in a protein-ligand complex in order to maximize its binding affinity. The results of the optimizer are used to predict beneficial growth vectors on the ligand, and good agreement is found between the predictions from the optimizer and a more rigorous free energy calculation, with a Spearman's rank order correlation of 0.59. The advantage of the method presented here is the significant speed up of over 10-fold compared to traditional free energy calculations and sublinear scaling with the number of growth vectors assessed. Where experimental data were available, mutations from hydrogen to a methyl group at sites highlighted by the optimizer were calculated with MBAR, and the mean unsigned error between experimental and calculated values of the binding free energy was 0.83 kcal/mol.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Capturing the Flexibility of a Protein-Ligand Complex: Binding Free Energies from Different Enhanced Sampling Techniques

Enhanced sampling techniques are a promising approach to obtain reliable binding free-energy profiles for flexible protein-ligand complexes from molecular dynamics (MD) simulations. To put four popular enhanced sampling techniques to a biologically relevant and challenging test, we studied the partial dissociation of an antigenic peptide from the Major Histocompatibility Complex I (MHC I) HLA-B*35:01 to systematically investigate the performance of umbrella sampling (US), replica exchange with solute tempering 2 (REST2), bias exchange umbrella sampling (BEUS, or replica-exchange umbrella sampling), and well-tempered metadynamics (MTD). With regard to the speed of sampling and convergence, the peptide-MHC I complex (pMHC I) under study showcases intrinsic strengths and weaknesses of the four enhanced sampling techniques used. We found that BEUS can best handle the sampling challenges that arise from the coexistence of an enthalpically and an entropically stabilized free-energy minimum in the pMHC I under study. These findings might also be relevant for other flexible biomolecular systems with competing enthalpically and entropically stabilized minima.

Implementation of adaptive integration method for free energy calculations in molecular systems

Estimating free energy differences by computer simulation is useful for a wide variety of applications such as virtual screening for drug design and for understanding how amino acid mutations modify protein interactions. However, calculating free energy differences remains challenging and often requires extensive trial and error and very long simulation times in order to achieve converged results. Here, we present an implementation of the adaptive integration method (AIM). We tested our implementation on two molecular systems and compared results from AIM to those from a suite of other methods. The model systems tested here include calculating the solvation free energy of methane, and the free energy of mutating the peptide GAG to GVG. We show that AIM is more efficient than other tested methods for these systems, that is, AIM results converge to a higher level of accuracy and precision for a given simulation time.

Alchemical Grid Dock (AlGDock): Binding Free Energy Calculations between Flexible Ligands and Rigid Receptors

Alchemical Grid Dock (AlGDock) is open-source software designed to compute the binding potential of mean force-the binding free energy between a flexible ligand and a rigid receptor-for a small organic ligand and a biological macromolecule. Multiple BPMFs can be used to rigorously compute binding affinities between flexible partners. AlGDock uses replica exchange between thermodynamic states at different temperatures and receptor-ligand interaction strengths. Receptor-ligand interaction energies are represented by interpolating precomputed grids. Thermodynamic states are adaptively initialized and adjusted on-the-fly to maintain adequate replica exchange rates. In demonstrative calculations, when the bound ligand is treated as fully solvated, AlGDock estimates BPMFs with a precision within 4 kT in 65% and within 8 kT for 91% of systems. It correctly identifies the native binding pose in 83% of simulations. Performance is sometimes limited by subtle differences in the important configuration space of sampled and targeted thermodynamic states. © 2019 Wiley Periodicals, Inc.

The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange-while displaying very small variance-can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations

Riboswitches can regulate gene expression by direct and specific interactions with ligands and have recently attracted interest as potential drug targets for antibacterial. In this work, molecular dynamics (MD) simulations, free energy perturbation (FEP) and molecular mechanics generalized Born surface area (MM-GBSA) methods were integrated to probe the effect of mutations on the binding of ligands to guanine riboswitch (GR). The results not only show that binding free energies predicted by FEP and MM-GBSA obtain an excellent correlation, but also indicate that mutations involved in the current study can strengthen the binding affinity of ligands GR. Residue-based free energy decomposition was applied to compute ligand-nucleotide interactions and the results suggest that mutations highly affect interactions of ligands with key nucleotides U22, U51 and C74. Dynamics analyses based on MD trajectories indicate that mutations not only regulate the structural flexibility but also change the internal motion modes of GR, especially for the structures J12, J23 and J31, which implies that the aptamer domain activity of GR is extremely plastic and thus readily tunable by nucleotide mutations. This study is expected to provide useful molecular basis and dynamics information for the understanding of the function of GR and possibility as potential drug targets for antibacterial.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Ligand Binding Thermodynamic Cycles: Hysteresis, the Locally Weighted Histogram Analysis Method, and the Overlapping States Matrix

Free energy perturbation (FEP) simulations have been widely applied to obtain predictions of the relative binding free energy for a series of congeneric ligands binding to the same receptor, which is an essential component for the lead optimization process in computer-aided drug discovery. In the case of several congeneric ligands forming a perturbation map involving a closed thermodynamic cycle, the summation of the estimated free energy change along each edge in the cycle using Bennett acceptance ratio (BAR) usually will deviate from zero due to systematic and random errors, which is the hysteresis of cycle closure. In this work, the advanced reweighting techniques binless weighted histogram analysis method (UWHAM) and locally weighted histogram analysis method (LWHAM) are applied to provide statistical estimators of the free energy change along each edge in order to eliminate the hysteresis effect. As an example, we analyze a closed thermodynamic cycle involving four congeneric ligands which bind to HIV-1 integrase, a promising target which has emerged for antiviral therapy. We demonstrate that, compared with FEP and BAR, more accurate and hysteresis-free estimates of free energy differences can be achieved by using UWHAM to find a single estimate of the density of states based on all of the data in the cycle. Furthermore, by comparison of LWHAM results obtained from the inclusion of different numbers of neighboring states with UWHAM estimation involving all the states, we show how to determine the optimal neighborhood size in the LWHAM analysis to balance the trade-offs between computational cost and accuracy of the free energy prediction. Even with the smallest neighborhood, LWHAM can improve the BAR free energy estimates using the same input data as BAR. We introduce an overlapping states matrix that is constructed by using the global jump formula of LWHAM and plot its heat map. The heat map provides a quantitative measure of the overlap between pairs of alchemical/thermodynamic states. We explain how to identify and improve the FEP calculations along the edges that most likely cause large systematic errors by using the heat map of the overlapping states matrix and by comparing the BAR and UWHAM estimates of the free energy change.

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding-unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.

Small-molecule targeting of MUSASHI RNA-binding activity in acute myeloid leukemia

The MUSASHI (MSI) family of RNA binding proteins (MSI1 and MSI2) contribute to a wide spectrum of cancers including acute myeloid leukemia. We find that the small molecule Ro 08-2750 (Ro) binds directly and selectively to MSI2 and competes for its RNA binding in biochemical assays. Ro treatment in mouse and human myeloid leukemia cells results in an increase in differentiation and apoptosis, inhibition of known MSI-targets, and a shared global gene expression signature similar to shRNA depletion of MSI2. Ro demonstrates in vivo inhibition of c-MYC and reduces disease burden in a murine AML leukemia model. Thus, we identify a small molecule that targets MSI's oncogenic activity. Our study provides a framework for targeting RNA binding proteins in cancer.

Molecular simulations by generalized-ensemble algorithms in isothermal-isobaric ensemble

Generalized-ensemble algorithms are powerful techniques for investigating biomolecules such as protein, DNA, lipid membrane, and glycan. The generalized-ensemble algorithms were originally developed in the canonical ensemble. On the other hand, not only temperature but also pressure is controlled in experiments. Additionally, pressure is used as perturbation to study relationship between function and structure of biomolecules. For this reason, it is important to perform efficient conformation sampling based on the isothermal-isobaric ensemble. In this article, we review a series of the generalized-ensemble algorithms in the isothermal-isobaric ensemble: multibaric-multithermal, pressure- and temperature-simulated tempering, replica-exchange, and replica-permutation methods. These methods achieve more efficient simulation than the conventional isothermal-isobaric simulation. Furthermore, the isothermal-isobaric generalized-ensemble simulation samples conformations of biomolecules from wider range of temperature and pressure. Thus, we can estimate physical quantities more accurately at any temperature and pressure values. The applications to the biomolecular system are also presented.

Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network

To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 μs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.