Fast,metadynamics‐based method for prediction of the stereochemistry‐dependent relative free energies of ligand–receptor interactions

The computational approach applicable for the molecular dynamics (MD)‐based techniques is proposed to predict the ligand–protein binding affinities dependent on the ligand stereochemistry. All possible stereoconfigurations are expressed in terms of one set of force‐field parameters [stereoconfiguration‐independent potential (SIP)], which allows for calculating all relative free energies by only single simulation. SIP can be used for studying diverse, stereoconfiguration‐dependent phenomena by means of various computational techniques of enhanced sampling. The method has been successfully tested on the β₂‐adrenergic receptor (β₂‐AR) binding the four fenoterol stereoisomers by both metadynamics simulations and replica‐exchange MD. Both the methods gave very similar results, fully confirming the presence of stereoselective effects in the fenoterol‐β₂‐AR interactions. However, the metadynamics‐based approach offered much better efficiency of sampling which allows for significant reduction of the unphysical region in SIP. © 2014 Wiley Periodicals, Inc.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.     

Ab initio prediction of protein-ligand binding structures by replica-exchange umbrella sampling simulations

We have developed a prediction method for the binding structures of ligands with proteins. Our method consists of three steps. First, replica-exchange umbrella sampling simulations are performed along the distance between a putative binding site of a protein and a ligand as the reaction coordinate. Second, we obtain the potential of mean force (PMF) of the unbiased system using the weighted histogram analysis method and determine the distance that corresponds to the global minimum of PMF. Third, structures that have this global-minimum distance and energy values around the average potential energy are collected and analyzed using the principal component analysis. We predict the binding structure as the global-minimum free energy state on the free energy landscapes along the two major principal component axes. As test cases, we applied our method to five protein-ligand complex systems. Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank.     

Hamiltonian replica‐permutation method and its applications to an alanine dipeptide and amyloid‐β(29–42) peptides

We propose the Hamiltonian replica‐permutation method (RPM) (or multidimensional RPM) for molecular dynamics and Monte Carlo simulations, in which parameters in the Hamiltonian are permuted among more than two replicas with the Suwa‐Todo algorithm. We apply the Coulomb RPM, which is one of realization of the Hamiltonian RPM, to an alanine dipeptide and to two amyloid‐β(29–42) molecules. The Hamiltonian RPM realizes more efficient sampling than the Hamiltonian replica‐exchange method. We illustrate the protein misfolding funnel of amyloid‐β(29–42) and reveal its dimerization pathways. © 2013 Wiley Periodicals, Inc.     

Application of replica exchange umbrella sampling to protein structure refinement of nontemplate models

We provide an assessment of a computational strategy for protein structure refinement that combines self‐guided Langevin dynamics with umbrella‐potential biasing replica exchange using the radius of gyration as a coordinate (R_g‐ReX). Eight structurally nonredundant proteins and their decoys were examined by sampling conformational space at room temperature using the CHARMM22/GBMV2 force field to generate the ensemble of structures. Two atomic statistical potentials (RWplus and DFIRE) were analyzed for structure identification and compared to the simulation force‐field potential. The results show that, while the R_g‐ReX simulations were able to sample conformational basins that were more structurally similar to the X‐ray crystallographic structures than the starting first‐order ranked decoys, the potentials failed to detect these basins from refinement. Of the three potential functions, RWplus yielded the highest accuracy for recognition of structures that refined to an average of nearly 20% increase in native contacts relative to the starting decoys. The overall performance of R_g‐ReX is compared to an earlier study of applying temperature‐based replica exchange to refine the same decoy sets and highlights the general challenge of achieving consistently the sampling and detection threshold of 70% fraction of native contacts. © 2013 Wiley Periodicals, Inc.     

GENESIS 1.1: A hybrid‐parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms

GENeralized‐Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all‐atom and coarse‐grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time‐step integration, and hybrid (CPU + GPU) computing. The string method and replica‐exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free‐energy pathway and obtaining free‐energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research. © 2017 Wiley Periodicals, Inc.     

How hot? Systematic convergence of the replica exchange method using multiple reservoirs

We have devised a systematic approach to converge a replica exchange molecular dynamics simulation by dividing the full temperature range into a series of higher temperature reservoirs and a finite number of lower temperature subreplicas. A defined highest temperature reservoir of equilibrium conformations is used to help converge a lower but still hot temperature subreplica, which in turn serves as the high‐temperature reservoir for the next set of lower temperature subreplicas. The process is continued until an optimal temperature reservoir is reached to converge the simulation at the target temperature. This gradual convergence of subreplicas allows for better and faster convergence at the temperature of interest and all intermediate temperatures for thermodynamic analysis, as well as optimizing the use of multiple processors. We illustrate the overall effectiveness of our multiple reservoir replica exchange strategy by comparing sampling and computational efficiency with respect to replica exchange, as well as comparing methods when converging the structural ensemble of the disordered Aβ_21‐30 peptide simulated with explicit water by comparing calculated Rotating Overhauser Effect Spectroscopy intensities to experimentally measured values. © 2009 Wiley Periodicals, Inc. J Comput Chem 31: 620–627, 2010     

Large‐scale asynchronous and distributed multidimensional replica exchange molecular simulations and efficiency analysis

We describe methods to perform replica exchange molecular dynamics (REMD) simulations asynchronously (ASyncRE). The methods are designed to facilitate large scale REMD simulations on grid computing networks consisting of heterogeneous and distributed computing environments as well as on homogeneous high‐performance clusters. We have implemented these methods on NSF (National Science Foundation) XSEDE (Extreme Science and Engineering Discovery Environment) clusters and BOINC (Berkeley Open Infrastructure for Network Computing) distributed computing networks at Temple University and Brooklyn College at CUNY (the City University of New York). They are also being implemented on the IBM World Community Grid. To illustrate the methods, we have performed extensive (more than 60 ms in aggregate) simulations for the beta‐cyclodextrin‐heptanoate host‐guest system in the context of one‐ and two‐dimensional ASyncRE, and we used the results to estimate absolute binding free energies using the binding energy distribution analysis method. We propose ways to improve the efficiency of REMD simulations: these include increasing the number of exchanges attempted after a specified molecular dynamics (MD) period up to the fast exchange limit and/or adjusting the MD period to allow sufficient internal relaxation within each thermodynamic state. Although ASyncRE simulations generally require long MD periods (>picoseconds) per replica exchange cycle to minimize the overhead imposed by heterogeneous computing networks, we found that it is possible to reach an efficiency similar to conventional synchronous REMD, by optimizing the combination of the MD period and the number of exchanges attempted per cycle. © 2015 Wiley Periodicals, Inc.     

Ergodicity and model quality in template‐restrained canonical and temperature/Hamiltonian replica exchange coarse‐grained molecular dynamics simulations of proteins

Molecular simulations restrained to single or multiple templates are commonly used in protein‐structure modeling. However, the restraints introduce additional barriers, thus impairing the ergodicity of simulations, which can affect the quality of the resulting models. In this work, the effect of restraint types and simulation schemes on ergodicity and model quality was investigated by performing template‐restrained canonical molecular dynamics (MD), multiplexed replica‐exchange molecular dynamics, and Hamiltonian replica exchange molecular dynamics (HREMD) simulations with the coarse‐grained UNRES force field on nine selected proteins, with pseudo‐harmonic log‐Gaussian (unbounded) or Lorentzian (bounded) restraint functions. The best ergodicity was exhibited by HREMD. It has been found that non‐ergodicity does not affect model quality if good templates are used to generate restraints. However, when poor‐quality restraints not covering the entire protein are used, the improved ergodicity of HREMD can lead to significantly improved protein models. © 2017 Wiley Periodicals, Inc.     

Calculation of binding free energies of inhibitors to plasmepsin II

An understanding at the atomic level of the driving forces of inhibitor binding to the protein plasmepsin (PM) II would be of interest to the development of drugs against malaria. To this end, three state of the art computational techniques to compute relative free energies-thermodynamic integration (TI), Hamiltonian replica-exchange (H-RE) TI, and comparison of bound versus unbound ligand energy and entropy-were applied to a protein-ligand system of PM II and several exo-3-amino-7-azabicyclo[2.2.1]heptanes and the resulting relative free energies were compared with values derived from experimental IC(50) values. For this large and flexible protein-ligand system, the simulations could not properly sample the relevant parts of the conformational space of the bound ligand, resulting in failure to reproduce the experimental data. Yet, the use of Hamiltonian replica exchange in conjunction with thermodynamic integration resulted in enhanced convergence and computational efficiency compared to standard thermodynamic integration calculations. The more approximate method of calculating only energetic and entropic contributions of the ligand in its bound and unbound states from conventional molecular dynamics (MD) simulations reproduced the major trends in the experimental binding free energies, which could be rationalized in terms of energetic and entropic characteristics of the different structural and physico-chemical properties of the protein and ligands.     

Oversampling Free Energy Perturbation Simulation in Determination of the Ligand-Binding Free Energy

Determination of the ligand-binding affinity is an extremely interesting problem. Normally, the free energy perturbation (FEP) method provides an appropriate result. However, it is of great interest to improve the accuracy and precision of this method. In this context, temperature replica exchange molecular dynamics implementation of the FEP computational approach, which we call replica exchange free energy perturbation (REP) was proposed. In particular, during REP simulations, the system can easily escape from being trapped in local minima by exchanging configurations with high temperatures, resulting in significant improvement in the accuracy and precision of protein–ligand binding affinity calculations. The distribution of the decoupling free energy was enlarged, and its mean values were decreased. This results in changes in the magnitude of the calculated binding free energies as well as in alteration in the binding mechanism. Moreover, the REP correlation coefficient with respect to experiment ( R^REP = 0.85 ± 0.15 ) is significantly boosted in comparison with the FEP one ( R^FEP = 0.64 ± 0.30 ). Furthermore, the root-mean-square error (RMSE) of REP is also smaller than FEP, RMSE^REP = 4.28 ± 0.69 versus RMSE^FEP = 5.80 ± 1.11 kcal/mol, respectively. © 2019 Wiley Periodicals, Inc.     

Free‐energy calculations of residue mutations in a tripeptide using various methods to overcome inefficient sampling

Previous free‐energy calculations have shown that the seemingly simple transformation of the tripeptide KXK to KGK in water holds some unobvious challenges concerning the convergence of the forward and backward thermodynamic integration processes (i.e., hysteresis). In the current study, the central residue X was either alanine, serine, glutamic acid, lysine, phenylalanine, or tyrosine. Interestingly, the transformation from alanine to glycine yielded the highest hysteresis in relation to the extent of the chemical change of the side chain. The reason for that could be attributed to poor sampling of φ₂/ψ₂ dihedral angles along the transformation. Altering the nature of alanine's C_β atom drastically improved the sampling and at the same time led to the identification of high energy barriers as cause for it. Consequently, simple strategies to overcome these barriers are to increase simulation time (computationally expensive) or to use enhanced sampling techniques such as Hamiltonian replica exchange molecular dynamics and one‐step perturbation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

The temperature intervals with global exchange of replicas empirical accelerated sampling method: parameter sensitivity and extension to a complex molecular system

The recently developed "temperature intervals with global exchange of replicas" (TIGER2) algorithm is an efficient replica-exchange sampling algorithm that provides the freedom to specify the number of replicas and temperature levels independently of the size of the system and temperature range to be spanned, thus making it particularly well suited for sampling molecular systems that are considered to be too large to be sampled using conventional replica exchange methods. Although the TIGER2 method is empirical in nature, when appropriately applied it is able to provide sampling that satisfies the balance condition and closely approximates a Boltzmann-weighted ensemble of states. In this work, we evaluated the influence of factors such as temperature range, temperature spacing, replica number, and sampling cycle design on the accuracy of a TIGER2 simulation based on molecular dynamics simulations of alanine dipeptide in implicit solvent. The influence of these factors is further examined by calculating the properties of a complex system composed of the B1 immunoglobulin-binding domain of streptococcal protein G (protein G) in aqueous solution. The accuracy of a TIGER2 simulation is particularly sensitive to the maximum temperature level selected for the simulation. A method to determine the appropriate maximum temperature level to be used in a TIGER2 simulation is presented.     

Efficacy of independence sampling in replica exchange simulations of ordered and disordered proteins

Recasting temperature replica exchange (T‐RE) as a special case of Gibbs sampling has led to a simple and efficient scheme for enhanced mixing (Chodera and Shirts, J. Chem. Phys., 2011, 135, 194110). To critically examine if T‐RE with independence sampling (T‐REis) improves conformational sampling, we performed T‐RE and T‐REis simulations of ordered and disordered proteins using coarse‐grained and atomistic models. The results demonstrate that T‐REis effectively increase the replica mobility in temperatures space with minimal computational overhead, especially for folded proteins. However, enhanced mixing does not translate well into improved conformational sampling. The convergences of thermodynamic properties interested are similar, with slight improvements for T‐REis of ordered systems. The study re‐affirms the efficiency of T‐RE does not appear to be limited by temperature diffusion, but by the inherent rates of spontaneous large‐scale conformational re‐arrangements. Due to its simplicity and efficacy of enhanced mixing, T‐REis is expected to be more effective when incorporated with various Hamiltonian‐RE protocols. © 2017 Wiley Periodicals, Inc.     

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.     

Connectivity and binding‐site recognition: Applications relevant to drug design

Here, we describe a family of methods based on residue–residue connectivity for characterizing binding sites and apply variants of the method to various types of protein–ligand complexes including proteases, allosteric‐binding sites, correctly and incorrectly docked poses, and inhibitors of protein–protein interactions. Residues within ligand‐binding sites have about 25% more contact neighbors than surface residues in general; high‐connectivity residues are found in contact with the ligand in 84% of all complexes studied. In addition, a k‐means algorithm was developed that may be useful for identifying potential binding sites with no obvious geometric or connectivity features. The analysis was primarily carried out on 61 protein–ligand structures from the MEROPS protease database, 250 protein–ligand structures from the PDBSelect (25%), and 30 protein–protein complexes. Analysis of four proteases with crystal structures for multiple bound ligands has shown that residues with high connectivity tend to have less variable side‐chain conformation. The relevance to drug design is discussed in terms of identifying allosteric‐binding sites, distinguishing between alternative docked poses and designing protein interface inhibitors. Taken together, this data indicate that residue–residue connectivity is highly relevant to medicinal chemistry. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Orthogonal sampling in free‐energy calculations of residue mutations in a tripeptide: TI versus λ‐LEUS

In a recent article (Bieler et al., J. Chem. Theory Comput. 2014, 10, 3006), we introduced a combination of λ‐dynamics and local‐elevation umbrella‐sampling termed λ‐LEUS to calculate free‐energy changes associated with alchemical processes using molecular dynamics simulations. This method was suggested to be more efficient than thermodynamic integration (TI), because the dynamical variation of the alchemical variable λ opens up pathways to circumvent barriers in the orthogonal space (defined by the N – 1 degrees of freedom that are not subjected to the sampling enhancement), a feature λ‐LEUS shares with Hamiltonian replica‐exchange (HR) approaches. However, the mutation considered, hydroquinone to benzene in water, was no real challenge in terms of orthogonal‐space properties, which were restricted to solvent‐relaxation processes. In the present article, we revisit the comparison between TI and λ‐LEUS considering non‐trivial mutations of the central residue X of a KXK tripeptide in water (with X = G, E, K, S, F, or Y). Side‐chain interactions that may include salt bridges, hydrogen bonds or steric clashes lead to slow relaxation in the orthogonal space, mainly in the two‐dimensional subspace spanned by the central and ψ dihedral angles of the peptide. The efficiency enhancement afforded by λ‐LEUS is confirmed in this more complex test system and can be attributed explicitly to the improved sampling of the orthogonal space. The sensitivity of the results to the nontrivial choices of a mass parameter and of a thermostat coupling time for the alchemical variable is also investigated, resulting in recommended ranges of 50 to 100 u nm² and 0.2 to 0.5 ps, respectively. © 2015 Wiley Periodicals, Inc.     

Q‐DockLHM: Low‐resolution refinement for ligand comparative modeling

The success of ligand docking calculations typically depends on the quality of the receptor structure. Given improvements in protein structure prediction approaches, approximate protein models now can be routinely obtained for the majority of gene products in a given proteome. Structure‐based virtual screening of large combinatorial libraries of lead candidates against theoretically modeled receptor structures requires fast and reliable docking techniques capable of dealing with structural inaccuracies in protein models. Here, we present Q‐Dock^LHM, a method for low‐resolution refinement of binding poses provided by FINDSITE^LHM, a ligand homology modeling approach. We compare its performance to that of classical ligand docking approaches in ligand docking against a representative set of experimental (both holo and apo) as well as theoretically modeled receptor structures. Docking benchmarks reveal that unlike all‐atom docking, Q‐Dock^LHM exhibits the desired tolerance to the receptor's structure deformation. Our results suggest that the use of an evolution‐based approach to ligand homology modeling followed by fast low‐resolution refinement is capable of achieving satisfactory performance in ligand‐binding pose prediction with promising applicability to proteome‐scale applications. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010