An improved replica-exchange sampling method: temperature intervals with global energy reassignment

Fundamental Redesign of the TIGER2hs Kernel to Address Severe Parameter Sensitivity

Replica exchange molecular dynamics simulations are one of the most popular approaches to enhance conformational sampling of molecular systems. Applications range from protein folding to protein-protein or other host-guest interactions, as well as binding free energy calculations. While these methods are computationally expensive, highly accurate results can be obtained. We recently developed TIGER2hs, an improved version of the temperature intervals with global exchange of replicas (TIGER2) algorithm. This method combines the replica-based enhanced sampling in an explicit solvent with a hybrid solvent energy evaluation. During the exchange attempts, bulk water is replaced by an implicit solvent model, allowing sampling with significantly less replicas than parallel tempering (REMD). This enables accurate enhanced sampling calculations with only a fraction of computational resources compared to REMD. Our latest results highlight several issues with sampling imbalance and parameter sensitivity within the original TIGER2 exchange algorithms that affect the overall state populations. A high sensitivity on replica number and maximum temperature is eliminated by changing to a pairwise exchange kernel (PE) without additional sorting. Simulations are controlled by adjusting the average temperature change per exchange ⟨ΔT/χ⟩ to below 30 K to mimic a controlled temperature mixing of replicas similar to REMD. Thus, this parameter provides an applicable property for selecting combinations of replica number and maximum temperature to adjust simulations for best accuracy, with flexible resource investment. This increases the robustness of the method and ensures results in excellent agreement with REMD, as demonstrated for three different peptides.

An Implementation of Replica Exchange with Dynamical Scaling for Efficient Large-Scale Simulations

An implementation of the replica exchange with dynamical scaling (REDS) method in the commonly used molecular dynamics program GROMACS is presented. REDS is a replica exchange method that requires fewer replicas than conventional replica exchange while still providing data over a range of temperatures and can be used in either constant volume or constant pressure ensembles. Details for running REDS simulations are given, and an application to the human islet amyloid polypeptide (hIAPP) 11-25 fragment shows that the model efficiently samples conformational space.

Ligand and structure based virtual screening of chemical databases to explore potent small molecule inhibitors against breast invasive carcinoma using recent computational technologies

Breast carcinoma is the most common invasive cancer to affect the women in the North America and the world. Cancer of breast is the number one cancer overall with estimated 1.5 lakh new cases during 2016. The success of the current endocrine therapies is often limited due to the development of resistance. Therefore, there is a need to develop new lead compounds for breast cancer treatment. As 70% of breast carcinoma is ER+, and it is well known previously that estrogen receptor alpha (ERα) is overexpressed in ER + cases, so in the current work we attempt to develop some novel potent analogues against ERα. To achieve this, we have adopted an integrative computational approach that involves multiple sequence alignment, virtual screening (ligand and structure based), molecular docking, fingerprint based clustering and molecular dynamics simulation. The approach envisaged vital information about the binding site residues, conserved sequence among different species, ligand and protein conformations, binding energy of compound to bind into the active site of the receptor. Molecular docking analysis revealed that some analogues exhibited significant binding towards ERα. The top docked complexes showing good docking scores, hydrogen bond and hydrophobic interactions were selected for molecular dynamics simulation studies. RMSD revealed that the systems were quite stable with RMSD value below 3 Å. The RMSF analysis calculated residue wise fluctuations and revealed that the residues are flexible enough to interact with the ligand. The residue at C-terminal showed more flexibility as compared to other residues. To confirm binding of these analogues, MMGBSA analysis was performed which revealed binding energy of the ligands. Further, per-residue decomposition energy analysis revealed that Glu353, Leu346, Leu387 and Arg394 contributed towards ligand binding. The results visibly indicated that MMGBSA can act as filter in virtual screening experiments and play a major role in facilitating drug discovery.

Potential target identification for breast cancer and screening of small molecule inhibitors: A bioinformatics approach

In the current study, we investigated the role of PAK1 (P21 (RAC1) Activated Kinase 1) gene in breast cancer and to this end, we performed differential gene expression analysis of PAK1 in breast cancer tissues compared to the normal adjacent tissue. We also studied its significance in protein-protein interaction (PPI) network, and analysed biological pathways, cellular processes, and role of PAK1 in different diseases. We found PAK1 to have significant role in breast cancer pathways such as integrin signaling, axonal guidance signaling, signaling by Rho family GTPases, ERK5 signaling. Additionally, it has been found as hub gene in PPI network, suggesting its possible regulatory role in breast carcinogenesis. Moreover, PAK1 had role in progression of various diseases as neoplasia, tumorigenesis, lymphatic neoplasia. Thereby, PAK1 can be used as a therapeutic target in breast cancer. Further, we put our efforts in identification of potential small molecules inhibitors against PAK1 by developing a composite virtual screening protocol involving molecular dynamics (MD) and molecular docking. The chemical library of compounds from NCI diversity sets, Pubchem and eMolecules were screened against PAK1 protein and hits which showed good binding affinity were considered for MD simulation study. Moreover, to assess binding of selected hits, MMGBSA (Molecular Mechanics-Generalized Born Surface Area) analysis was performed using AMBER (Assisted Model Building with Energy Refinement) package. MMGBSA calculations exhibited that the identified ligands showed good binding affinity with PAK1. HighlightsThe PAK1 has been found to be upregulated in breast cancer samples and is a potential oncogene playing role in different cellular functions and processes.The molecular docking studies revealed ligands showed good binding affinity towards PAK1 protein.The residues Glu345, Leu347, Thr406, Asp299, Asp393 and Gly350 were found to make H-bond interactions with small molecule inhibitors.The residues Ile276, Val284, Ala297, Tyr346, Leu396 and Asp407 were found to make hydrophobic interactions.The RMSD analysis confirmed stability of complexes throughout 40 ns production period.The MD simulations studies revealed the binding site flexibility, binding free energy of complexes and per-residue contribution in ligand binding.Communicated by Ramaswamy H. Sarma.

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.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Pathways to Structure-Property Relationships of Peptide-Materials Interfaces: Challenges in Predicting Molecular Structures

An in-depth appreciation of how to manipulate the molecular-level recognition between peptides and aqueous materials interfaces, including nanoparticles, will advance technologies based on self-organized metamaterials for photonics and plasmonics, biosensing, catalysis, energy generation and harvesting, and nanomedicine. Exploitation of the materials-selective binding of biomolecules is pivotal to success in these areas and may be particularly key to producing new hierarchically structured biobased materials. These applications could be accomplished by realizing preferential adsorption of a given biomolecule onto one materials composition over another, one surface facet over another, or one crystalline polymorph over another. Deeper knowledge of the aqueous abiotic-biotic interface, to establish clear structure-property relationships in these systems, is needed to meet this goal. In particular, a thorough structural characterization of the surface-adsorbed peptides is essential for establishing these relationships but can often be challenging to accomplish via experimental approaches alone. In addition to myriad existing challenges associated with determining the detailed molecular structure of any molecule adsorbed at an aqueous interface, experimental characterization of materials-binding peptides brings new, complex challenges because many materials-binding peptides are thought to be intrinsically disordered. This means that these peptides are not amenable to experimental techniques that rely on the presence of well-defined secondary structure in the peptide when in the adsorbed state. To address this challenge, and in partnership with experiment, molecular simulations at the atomistic level can bring complementary and critical insights into the origins of this abiotic/biotic recognition and suggest routes for manipulating this phenomenon to realize new types of hybrid materials. For the reasons outlined above, molecular simulation approaches also face challenges in their successful application to model the biotic-abiotic interface, related to several factors. For instance, simulations require a plausible description of the chemistry and the physics of the interface, which comprises two very different states of matter, in the presence of liquid water. Also, it is essential that the conformational ensemble be comprehensively characterized under these conditions; this is especially challenging because intrinsically disordered peptides do not typically admit one single structure or set of structures. Moreover, a plausible structural model of the substrate is required, which may require a high level of detail, even for single-element materials such as Au surfaces or graphene. Developing and applying strategies to make credible predictions of the conformational ensemble of adsorbed peptides and using these to construct structure-property relationships of these interfaces have been the goals of our efforts. We have made substantial progress in developing interatomic potentials for these interfaces and adapting advanced conformational sampling approaches for these purposes. This Account summarizes our progress in the development and deployment of interfacial force fields and molecular simulation techniques for the purpose of elucidating these insights at biomolecule-materials interfaces, using examples from our laboratories ranging from noble-metal interfaces to graphitic substrates (including carbon nanotubes and graphene) and oxide materials (such as titania). In addition to the well-established application areas of plasmonic materials, biosensing, and the production of medical implant materials, we outline new directions for this field that have the potential to bring new advances in areas such as energy materials and regenerative medicine.

Application of advanced sampling and analysis methods to predict the structure of adsorbed protein on a material surface

The use of standard molecular dynamics simulation methods to predict the interactions of a protein with a material surface have the inherent limitations of lacking the ability to determine the most likely conformations and orientations of the adsorbed protein on the surface and to determine the level of convergence attained by the simulation. In addition, standard mixing rules are typically applied to combine the nonbonded force field parameters of the solution and solid phases the system to represent interfacial behavior without validation. As a means to circumvent these problems, the authors demonstrate the application of an efficient advanced sampling method (TIGER2A) for the simulation of the adsorption of hen egg-white lysozyme on a crystalline (110) high-density polyethylene surface plane. Simulations are conducted to generate a Boltzmann-weighted ensemble of sampled states using force field parameters that were validated to represent interfacial behavior for this system. The resulting ensembles of sampled states were then analyzed using an in-house-developed cluster analysis method to predict the most probable orientations and conformations of the protein on the surface based on the amount of sampling performed, from which free energy differences between the adsorbed states were able to be calculated. In addition, by conducting two independent sets of TIGER2A simulations combined with cluster analyses, the authors demonstrate a method to estimate the degree of convergence achieved for a given amount of sampling. The results from these simulations demonstrate that these methods enable the most probable orientations and conformations of an adsorbed protein to be predicted and that the use of our validated interfacial force field parameter set provides closer agreement to available experimental results compared to using standard CHARMM force field parameterization to represent molecular behavior at the interface.

Peptide interactions with zigzag edges in graphene

Recognition and manipulation of graphene edges enable the control of physical properties of graphene-based devices. Recently, the authors have identified a peptide that preferentially binds to graphene edges from a combinatorial peptide library. In this study, the authors examine the functional basis for the edge binding peptide using experimental and computational methods. The effect of amino acid substitution, sequence context, and solution pH value on the binding of the peptide to graphene has been investigated. The N-terminus glutamic acid residue plays a key role in recognizing and binding to graphene edges. The protonation, substitution, and positional context of the glutamic acid residue impact graphene edge-binding. Our findings provide insights into the binding mechanisms and the design of peptides for recognizing and functionalizing graphene edges.

TIGER2 with solvent energy averaging (TIGER2A): An accelerated sampling method for large molecular systems with explicit representation of solvent

The recently developed "temperature intervals with global exchange of replicas" (TIGER2) accelerated sampling method is found to have inaccuracies when applied to systems with explicit solvation. This inaccuracy is due to the energy fluctuations of the solvent, which cause the sampling method to be less sensitive to the energy fluctuations of the solute. In the present work, the problem of the TIGER2 method is addressed in detail and a modification to the sampling method is introduced to correct this problem. The modified method is called "TIGER2 with solvent energy averaging," or TIGER2A. This new method overcomes the sampling problem with the TIGER2 algorithm and is able to closely approximate Boltzmann-weighted sampling of molecular systems with explicit solvation. The difference in performance between the TIGER2 and TIGER2A methods is demonstrated by comparing them against analytical results for simple one-dimensional models, against replica exchange molecular dynamics (REMD) simulations for sampling the conformation of alanine dipeptide and the folding behavior of (AAQAA)3 peptide in aqueous solution, and by comparing their performance in sampling the behavior of hen egg-white lysozyme in aqueous solution. The new TIGER2A method solves the problem caused by solvent energy fluctuations in TIGER2 while maintaining the two important characteristics of TIGER2, i.e., (1) using multiple replicas sampled at different temperature levels to help systems efficiently escape from local potential energy minima and (2) enabling the number of replicas used for a simulation to be independent of the size of the molecular system, thus providing an accelerated sampling method that can be used to efficiently sample systems considered too large for the application of conventional temperature REMD.

Understanding and Designing the Gold-Bio Interface: Insights from Simulations

Gold nanoparticles (AuNPs) are an integral part of many exciting and novel biomedical applications, sparking the urgent need for a thorough understanding of the physicochemical interactions occurring between these inorganic materials, their functional layers, and the biological species they interact with. Computational approaches are instrumental in providing the necessary molecular insight into the structural and dynamic behavior of the Au-bio interface with spatial and temporal resolutions not yet achievable in the laboratory, and are able to facilitate a rational approach to AuNP design for specific applications. A perspective of the current successes and challenges associated with the multiscale computational treatment of Au-bio interfacial systems, from electronic structure calculations to force field methods, is provided to illustrate the links between different approaches and their relationship to experiment and applications.

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.

Free Energy Landscapes of Alanine Oligopeptides in Rigid-Body and Hybrid Water Models

Replica exchange molecular dynamics is used to study the effect of different rigid-body (mTIP3P, TIP4P, SPC/E) and hybrid (H1.56, H3.00) water models on the conformational free energy landscape of the alanine oligopeptides (acAnme and acA5nme), in conjunction with the CHARMM22 force field. The free energy landscape is mapped out as a function of the Ramachandran angles. In addition, various secondary structure metrics, solvation shell properties, and the number of peptide-solvent hydrogen bonds are monitored. Alanine dipeptide is found to have similar free energy landscapes in different solvent models, an insensitivity which may be due to the absence of possibilities for forming i-(i + 4) or i-(i + 3) intrapeptide hydrogen bonds. The pentapeptide, acA5nme, where there are three intrapeptide backbone hydrogen bonds, shows a conformational free energy landscape with a much greater degree of sensitivity to the choice of solvent model, though the three rigid-body water models differ only quantitatively. The pentapeptide prefers nonhelical, non-native PPII and β-sheet populations as the solvent is changed from SPC/E to the less tetrahedral liquid (H1.56) to an LJ-like liquid (H3.00). The pentapeptide conformational order metrics indicate a preference for open, solvent-exposed, non-native structures in hybrid solvent models at all temperatures of study. The possible correlations between the properties of solvent models and secondary structure preferences of alanine oligopeptides are discussed, and the competition between intrapeptide, peptide-solvent, and solvent-solvent hydrogen bonding is shown to be crucial in the relative free energies of different conformers.

Replica exchange with guided annealing for accelerated sampling of disordered protein conformations

We critically examine a recently proposed convective replica exchange (cRE) method for enhanced sampling of protein conformation based on theoretical and numerical analysis. The results demonstrate that cRE and related replica exchange with guided annealing (RE-GA) schemes lead to unbalanced exchange attempt probabilities and break detailed balance whenever the system undergoes slow conformational transitions (relative to the temperature diffusion timescale). Nonetheless, numerical simulations suggest that approximate canonical ensembles can be generated for systems with small conformational transition barriers. This suggests that RE-GA maybe suitable for simulating intrinsically disordered proteins, an important class of newly recognized functional proteins. The efficacy of RE-GA is demonstrated by calculating the conformational ensembles of intrinsically disordered kinase inducible domain protein. The results show that RE-GA helps the protein to escape nonspecific compact states more efficiently and provides several fold speedups in generating converged and largely correct ensembles compared to the standard temperature RE.

Optimization of Umbrella Sampling Replica Exchange Molecular Dynamics by Replica Positioning

The positioning of sampling windows in an umbrella sampling simulation has an effect on the rate of convergence and computational efficiency. When such simulation is coupled with a Hamiltonian replica exchange setup, we show that such positioning can be optimized for maximal convergence of the results. We present a method for estimating the exchange acceptance ratio (EAR) between two arbitrary positions on a reaction coordinate in umbrella sampling replica exchange (USRE) molecular dynamics (MD). We designed a scoring function to optimize the position of the set of replicas (windows). By maximizing the scoring function, we make EAR the same for all neighbor replica pairs, increasing the efficiency of the method. We tested our algorithm by sampling a torsion for butane in implicit solvent and by studying a salt bridge in explicit solvent. We found that the optimized set of replicas recovers the correct free energy profile much faster than for equally spaced umbrellas.

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

Development of molecular simulation methods to accurately represent protein-surface interactions: The effect of pressure and its determination for a system with constrained atoms

When performing molecular dynamics simulations for a system with constrained (fixed) atoms, traditional isobaric algorithms (e.g., NPT simulation) often cannot be used. In addition, the calculation of the internal pressure of a system with fixed atoms may be highly inaccurate due to the nonphysical nature of the atomic constraints and difficulties in accurately defining the volume occupied by the unconstrained atoms in the system. The inability to properly set and control pressure can result in substantial problems for the accurate simulation of condensed-phase systems if the behavior of the system (e.g., peptide/protein adsorption) is sensitive to pressure. To address this issue, the authors have developed an approach to accurately determine the internal pressure for a system with constrained atoms. As the first step in this method, a periodically extendable portion of the mobile phase of the constrained system (e.g., the solvent atoms) is used to create a separate unconstrained system for which the pressure can be accurately calculated. This model system is then used to create a pressure calibration plot for an intensive local effective virial parameter for a small volume cross section or "slab" of the system. Using this calibration plot, the pressure of the constrained system can then be determined by calculating the virial parameter for a similarly sized slab of mobile atoms. In this article, the authors present the development of this method and demonstrate its application using the CHARMM molecular simulation program to characterize the adsorption behavior of a peptide in explicit water on a hydrophobic surface whose lattice spacing is maintained with atomic constraints. The free energy of adsorption for this system is shown to be dramatically influenced by pressure, thus emphasizing the importance of properly maintaining the pressure of the system for the accurate simulation of protein-surface interactions.

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.

TIGER2: an improved algorithm for temperature intervals with global exchange of replicas

An empirical sampling method for molecular simulation based on "temperature intervals with global exchange of replicas" (TIGER2) has been developed to reduce the high demand for computational resources and the low computational efficiency of the conventional replica-exchange molecular dynamics (REMD) method. This new method overcomes the limitation of its previous version, called TIGER, which requires the assumption of constant heat capacity during quenching of replicas from elevated temperatures to the baseline temperature. The robustness of the TIGER2 method is examined by comparing it against a Metropolis Monte Carlo simulation for sampling the conformational distribution of a single butane molecule in vacuum, a REMD simulation for sampling the behavior of alanine dipeptide in explicit solvent, and REMD simulations for sampling the folding behavior of two peptides, (AAQAA)(3) and chignolin, in implicit solvent. The agreement between the results from these conventional sampling methods and the TIGER2 simulations indicates that the TIGER2 algorithm is able to closely approximate a Boltzmann-weighted ensemble of states for these systems but without the limiting assumptions that were required for the original TIGER algorithm. TIGER2 is an efficient replica-exchange sampling method that enables the number of replicas that are used for a replica-exchange simulation to be substantially reduced compared to the conventional REMD method.

An optimized replica exchange molecular dynamics method

We introduce a new way to perform swaps between replicas in replica exchange molecular dynamics simulations. The method is based on a generalized canonical probability distribution function and flattens the potential of mean force along the temperature coordinate, such that a random walk in temperature space is achieved. Application to a Go model of protein A showed that the method is more efficient than conventional replica exchange. The method results in a constant probability distribution of the replicas over the thermostats, yields a minimum round-trip time between extremum temperatures, and leads to faster ergodic convergence.

Calculation of protein heat capacity from replica-exchange molecular dynamics simulations with different implicit solvent models

The heat capacity has played a major role in relating microscopic and macroscopic properties of proteins and their disorder-order phase transition of folding. Its calculation by atomistic simulation methods remains a significant challenge due to the complex and dynamic nature of protein structures, their solvent environment, and configurational averaging. To better understand these factors on calculating a protein heat capacity, we provide a comparative analysis of simulation models that differ in their implicit solvent description and force-field resolution. Our model protein system is the src Homology 3 (SH3) domain of alpha-spectrin, and we report a series of 10 ns replica-exchange molecular dynamics simulations performed at temperatures ranging from 298 to 550 K, starting from the SH3 native structure. We apply the all-atom CHARMM22 force field with different modified analytical generalized Born solvent models (GBSW and GBMV2) and compare these simulation models with the distance-dependent dielectric screening of charge-charge interactions. A further comparison is provided with the united-atom CHARMM19 plus a pairwise GB model. Unfolding-folding transition temperatures of SH3 were estimated from the temperature-dependent profiles of the heat capacity, root-mean-square distance from the native structure, and the fraction of native contacts, each calculated from the density of states by using the weighted histogram analysis method. We observed that, for CHARMM22, the unfolding transition and energy probability density were quite sensitive to the implicit solvent description, in particular, the treatment of the protein-solvent dielectric boundary in GB models and their surface-area-based hydrophobic term. Among the solvent models tested, the calculated melting temperature varied in the range 353-438 K and was higher than the experimental value near 340 K. A reformulated GBMV2 model of employing a smoother molecular-volume dielectric interface was the most accurate in reproducing the native conformation and a two-state folding landscape, although the melting transition temperature did not show the smallest deviation from experiment. For the lower-resolution CHARMM19/GB model, the simulations failed to yield a bimodal energy distribution, yet the melting temperature was observed to be a good estimate of higher-resolution simulation models. We also demonstrate that a careful analysis of a relatively long simulation is necessary to avoid trapping in local minima and to find a true thermodynamic transition temperature.

Molecular simulation of protein-surface interactions: benefits, problems, solutions, and future directions

While the importance of protein adsorption to materials surfaces is widely recognized, little is understood at this time regarding how to design surfaces to control protein adsorption behavior. All-atom empirical force field molecular simulation methods have enormous potential to address this problem by providing an approach to directly investigate the adsorption behavior of peptides and proteins at the atomic level. As with any type of technology, however, these methods must be appropriately developed and applied if they are to provide realistic and useful results. Three issues that are particularly important for the accurate simulation of protein adsorption behavior are the selection of a valid force field to represent the atomic-level interactions involved, the accurate representation of solvation effects, and system sampling. In this article, each of these areas is addressed and future directions for continued development are presented.