Enhanced conformational sampling method for proteins based on the TaBoo SeArch algorithm: Application to the folding of a mini‐protein,chignolin

The conformational samplings are indispensible for obtaining reliable canonical ensembles, which provide statistical averages of physical quantities such as free energies. However, the samplings of vast conformational space of biomacromolecules by conventional molecular dynamics (MD) simulations might be insufficient, due to their inadequate accessible time‐scales for investigating biological functions. Therefore, the development of methodologies for enhancing the conformational sampling of biomacromolecules still remains as a challenging issue in computational biology. To tackle this problem, we newly propose an efficient conformational search method, which is referred as TaBoo SeArch (TBSA) algorithm. In TBSA, an inverse energy histogram is used to select seeds for the conformational resampling so that states with high frequencies are inhibited, while states with low frequencies are efficiently sampled to explore the unvisited conformational space. As a demonstration, TBSA was applied to the folding of a mini‐protein, chignolin, and automatically sampled the native structure (C_α root mean square deviation < 1.0 Å) with nanosecond order computational costs started from a completely extended structure, although a long‐time 1‐µs normal MD simulation failed to sample the native structure. Furthermore, a multiscale free energy landscape method based on the conformational sampling of TBSA were quantitatively evaluated through free energy calculations with both implicit and explicit solvent models, which enable us to find several metastable states on the folding landscape. © 2015 Wiley Periodicals, Inc.     

Unfolding the conformational behavior of peptide dendrimers: insights from molecular dynamics simulations

We present here the first comprehensive structural characterization of peptide dendrimers using molecular simulation methods. Multiple long molecular dynamics simulations are used to extensively sample the conformational preferences of five third-generation peptide dendrimers, including some known to bind aquacobalamine. We start by analyzing the compactness of the conformations thus sampled using their radius of gyration profiles. A more detailed analysis is then performed using dissimilarity measures, principal coordinate analysis, and free energy landscapes, with the aim of identifying groups of similar conformations. The results point to a high conformational flexibility of these molecules, with no clear "folded state", although two markedly distinct behaviors were found: one of the dendrimers displayed mostly compact conformations clustered into distinct basins (rough landscape), while the remaining dendrimers displayed mainly noncompact conformations with no significant clustering (downhill landscape). This study brings new insight into the conformational behavior of peptide dendrimers and may provide better routes for their functional design. In particular, we propose a yet unsynthesized peptide dendrimer that might exhibit enhanced ability to coordinate aquocobalamin.     

Hybrid Monte Carlo with adaptive temperature in mixed-canonical ensemble: Efficient conformational analysis of RNA

A hybrid Monte Carlo method with adaptive temperature choice is presented that exactly generates the distribution of a mixed-canonical ensemble composed of two canonical ensembles at low and high temperature. The analysis of resulting Markov chains with the reweighting technique shows an efficient sampling of the canonical distribution at low temperature whereas the high temperature component facilitates conformational transitions, which allows shorter simulation times. The algorithm is tested by comparing analytical and numerical results for the small n-butane molecule before simulations are performed for a triribonucleotide. Sampling the complex multiminima energy landscape of this small RNA segment, we observe enforced crossing of energy barriers. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1689–1697, 1998     

Single molecule conformational dynamics of adenylate kinase: energy landscape, structural correlations, and transition state ensembles

We developed a coarse grained two-well model to study the single molecule protein conformational dynamics in microscopic detail at the residue level, overcoming the often encountered computational bottleneck. In particular, we explored the underlying conformational energy landscape of adenylate kinase, a crucial protein for signal transduction in the cell, and identified two major kinetic pathways for the conformational switch between open and closed states through either the intermediate state or the transient state. Based on the parameters fitted to the room-temperature experimental data, we predicted open and closed kinetic rates at the whole temperature ranges from 10 to 50 degrees C, which agree well with the experimental turnover numbers. After uncovering the underlying mechanism for conformational dynamics and exploring the structural correlations, we found the crucial dynamical interplay between the nucleoside monophosphate binding domain (NMP) and the ATP-binding domain (LID) in controlling the conformational switch. The key residues and contacts responsible for the conformational transitions are identified by following the time evolution of the two-dimensional spatial contact maps and characterizing the transition state as well as intermediate structure ensembles through phi value analysis. Our model provides a general framework to study the conformational dynamics of biomolecules and can be applied to many other systems.     

The Effect of Pulling and Twisting Forces on Chameleon Sequence Peptides**

Chameleon sequences are amino acid sequences found in several distinct configurations in experiment. They challenge our understanding of the link between sequence and structure, and provide insight into structural competition in proteins. Here, we study the energy landscapes for three such sequences, and interrogate how pulling and twisting forces impact the available structural ensembles. Chameleon sequences do not necessarily exhibit multiple structural ensembles on a multifunnel energy landscape when we consider them in isolation. The application of even small forces leads to drastic changes in the energy landscapes. For pulling forces, we observe transitions from helical to extended structures in a very small span of forces. For twisting forces, the picture is much more complex, and highly dependent on the magnitude and handedness of the applied force as well as the reference angle for the twist. Depending on these parameters, more complex and more simplistic energy landscapes are observed alongside more and less diverse structural ensembles. The impact of even small forces is significant, confirming their likely role in folding events. In addition, small forces exerted by the remaining scaffold of a protein may be sufficient to lead to the adoption of a specific structural ensemble by a chameleon sequence.     

Estimating kinetic rates from accelerated molecular dynamics simulations: alanine dipeptide in explicit solvent as a case study

Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.     

Generalized simulated tempering realized on expanded ensembles of non-Boltzmann weights

A generalized version of the simulated tempering operated in the expanded ensembles of non-Boltzmann weights has been proposed to mitigate a quasiergodicity problem occurring in simulations of rough energy landscapes. In contrast to conventional simulated tempering employing the Boltzmann weight, our method utilizes a parametrized, generalized distribution as a workhorse for stochastic exchanges of configurations and subensembles transitions, which allows a considerable enhancement for the rate of convergence of Monte Carlo and molecular dynamics simulations using delocalized weights. A feature of our method is that the exploration of the parameter space encouraging subensembles transitions is greatly accelerated using the dynamic update scheme for the weight via the average guide specific to the energy distribution. The performance and characteristic feature of our method have been validated in the liquid-solid transition of Lennard-Jones clusters and the conformational sampling of alanine dipeptide by taking two types of Tsallis [C. Tsallis, J. Stat. Phys. 52, 479 (1988)] expanded ensembles associated with different parametrization schemes.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

A skewed-momenta method to efficiently generate conformational-transition trajectories

We present a novel computational method, the skewed-momenta method (Skew'M), which applies a bias to the Maxwell distribution of initial momenta used to generate ensembles of trajectories. As a result, conformational transitions are accentuated and kinetic properties are calculated more effectively. The connection to the related puddle jumping method is discussed. A reweighting scheme permits the exact calculation of kinetic properties. Applications are presented for the rapid calculation of rate constants for molecular isomerization, and for the efficient reconstruction of free-energy profiles using a straightforward modification of the Jarzynski identity.     

RNACluster: An integrated tool for RNA secondary structure comparison and clustering

RNA structure comparison is a fundamental problem in structural biology, structural chemistry, and bioinformatics. It can be used for analysis of RNA energy landscapes, conformational switches, and facilitating RNA structure prediction. The purpose of our integrated tool RNACluster is twofold: to provide a platform for computing and comparison of different distances between RNA secondary structures, and to perform cluster identification to derive useful information of RNA structure ensembles, using a minimum spanning tree (MST) based clustering algorithm. RNACluster employs a cluster identification approach based on a MST representation of the RNA ensemble data and currently supports six distance measures between RNA secondary structures. RNACluster provides a user-friendly graphical interface to allow a user to compare different structural distances, analyze the structure ensembles, and visualize predicted structural clusters.     

Visualizing energy landscapes with metric disconnectivity graphs

The visualization of multidimensional energy landscapes is important, providing insight into the kinetics and thermodynamics of a system, as well the range of structures a system can adopt. It is, however, highly nontrivial, with the number of dimensions required for a faithful reproduction of the landscape far higher than can be represented in two or three dimensions. Metric disconnectivity graphs provide a possible solution, incorporating the landscape connectivity information present in disconnectivity graphs with structural information in the form of a metric. In this study, we present a new software package, PyConnect, which is capable of producing both disconnectivity graphs and metric disconnectivity graphs in two or three dimensions. We present as a test case the analysis of the 69‐bead BLN coarse‐grained model protein and show that, by choosing appropriate order parameters, metric disconnectivity graphs can resolve correlations between structural features on the energy landscape with the landscapes energetic and kinetic properties. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Energy landscapes of dynamic ensembles of rolling triplet repeat bulge loops: implications for DNA expansion associated with disease states

DNA repeat domains can form ensembles of canonical and noncanonical states, including stable and metastable DNA secondary structures. Such sequence-induced structural diversity creates complex conformational landscapes for DNA processing pathways, including those triplet expansion events that accompany replication, recombination, and/or repair. Here we demonstrate further levels of conformational complexity within repeat domains. Specifically, we show that bulge loop structures within an extended repeat domain can form dynamic ensembles containing a distribution of loop positions, thereby yielding families of positional loop isomers, which we designate as "rollamers". Our fluorescence, absorbance, and calorimetric data are consistent with loop migration/translocation between sites within the repeat domain ("rollamerization"). We demonstrate that such "rollameric" migration of bulge loops within repeat sequences can invade and disrupt previously formed base-paired domains via an isoenthalpic, entropy-driven process. We further demonstrate that destabilizing abasic lesions alter the loop distributions so as to favor "rollamers" with the lesion positioned at the duplex/loop junction, sites where the flexibility of the abasic "universal hinge" relaxes unfavorable interactions and/or facilitates topological accommodation. Another strategic siting of an abasic site induces directed loop migration toward denaturing domains, a phenomenon that merges destabilizing domains. In the aggregate, our data reveal that dynamic ensembles within repeat domains profoundly impact the overall energetics of such DNA constructs as well as the distribution of states by which they denature/renature. These static and dynamic influences within triplet repeat domains expand the conformational space available for selection and targeting by the DNA processing machinery. We propose that such dynamic ensembles and their associated impact on DNA properties influence pathways that lead to DNA expansion.     

Sparsity‐weighted outlier FLOODing (OFLOOD) method: Efficient rare event sampling method using sparsity of distribution

As an extension of the Outlier FLOODing (OFLOOD) method [Harada et al., J. Comput. Chem. 2015, 36, 763], the sparsity of the outliers defined by a hierarchical clustering algorithm, FlexDice, was considered to achieve an efficient conformational search as sparsity‐weighted “OFLOOD.” In OFLOOD, FlexDice detects areas of sparse distribution as outliers. The outliers are regarded as candidates that have high potential to promote conformational transitions and are employed as initial structures for conformational resampling by restarting molecular dynamics simulations. When detecting outliers, FlexDice defines a rank in the hierarchy for each outlier, which relates to sparsity in the distribution. In this study, we define a lower rank (first ranked), a medium rank (second ranked), and the highest rank (third ranked) outliers, respectively. For instance, the first‐ranked outliers are located in a given conformational space away from the clusters (highly sparse distribution), whereas those with the third‐ranked outliers are nearby the clusters (a moderately sparse distribution). To achieve the conformational search efficiently, resampling from the outliers with a given rank is performed. As demonstrations, this method was applied to several model systems: Alanine dipeptide, Met‐enkephalin, Trp‐cage, T4 lysozyme, and glutamine binding protein. In each demonstration, the present method successfully reproduced transitions among metastable states. In particular, the first‐ranked OFLOOD highly accelerated the exploration of conformational space by expanding the edges. In contrast, the third‐ranked OFLOOD reproduced local transitions among neighboring metastable states intensively. For quantitatively evaluations of sampled snapshots, free energy calculations were performed with a combination of umbrella samplings, providing rigorous landscapes of the biomolecules. © 2015 Wiley Periodicals, Inc.     

Hybridizing rapidly exploring random trees and basin hopping yields an improved exploration of energy landscapes

The number of local minima of the potential energy landscape (PEL) of molecular systems generally grows exponentially with the number of degrees of freedom, so that a crucial property of PEL exploration algorithms is their ability to identify local minima, which are low lying and diverse. In this work, we present a new exploration algorithm, retaining the ability of basin hopping (BH) to identify local minima, and that of transition based rapidly exploring random trees (T‐RRT) to foster the exploration of yet unexplored regions. This ability is obtained by interleaving calls to the extension procedures of BH and T‐RRT, and we show tuning the balance between these two types of calls allows the algorithm to focus on low lying regions. Computational efficiency is obtained using state‐of‐the art data structures, in particular for searching approximate nearest neighbors in metric spaces. We present results for the BLN69, a protein model whose conformational space has dimension 207 and whose PEL has been studied exhaustively. On this system, we show that the propensity of our algorithm to explore low lying regions of the landscape significantly outperforms those of BH and T‐RRT. © 2015 Wiley Periodicals, Inc.     

Calculation of free energy landscapes: A histogram reweighted metadynamics approach

We present an efficient method for the calculation of free energy landscapes. Our approach involves a history‐dependent bias potential, which is evaluated on a grid. The corresponding free energy landscape is constructed via a histogram reweighting procedure a posteriori. Because of the presence of the bias potential, it can be also used to accelerate rare events. In addition, the calculated free energy landscape is not restricted to the actual choice of collective variables and can in principle be extended to auxiliary variables of interest without further numerical effort. The applicability is shown for several examples. We present numerical results for the alanine dipeptide and the Met‐Enkephalin in explicit solution to illustrate our approach. Furthermore, we derive an empirical formula that allows the prediction of the computational cost for the ordinary metadynamics variant in comparison with our approach, which is validated by a dimensionless representation. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Comparison of free energy methods for molecular systems

We present a detailed comparison of computational efficiency and precision for several free energy difference (DeltaF) methods. The analysis includes both equilibrium and nonequilibrium approaches, and distinguishes between unidirectional and bidirectional methodologies. We are primarily interested in comparing two recently proposed approaches, adaptive integration, and single-ensemble path sampling to more established methodologies. As test cases, we study relative solvation free energies of large changes to the size or charge of a Lennard-Jones particle in explicit water. The results show that, for the systems used in this study, both adaptive integration and path sampling offer unique advantages over the more traditional approaches. Specifically, adaptive integration is found to provide very precise long-simulation DeltaF estimates as compared to other methods used in this report, while also offering rapid estimation of DeltaF. The results demonstrate that the adaptive integration approach is the best overall method for the systems studied here. The single-ensemble path sampling approach is found to be superior to ordinary Jarzynski averaging for the unidirectional, "fast-growth" nonequilibrium case. Closer examination of the path sampling approach on a two-dimensional system suggests it may be the overall method of choice when conformational sampling barriers are high. However, it appears that the free energy landscapes for the systems used in this study have rather modest configurational sampling barriers.     

Molecular properties from conformational ensembles. 1. Dipole moments of molecules with multiple internal rotations

We present a novel method for constructing the stable conformational space of small molecules with many rotatable bonds that uses our iterative stochastic elimination (ISE) algorithm, a robust stochastic search method capable of finding ensembles of best solutions for large combinatorial problems. To validate the method, we show that ISE reproduces the best conformers found in a fully exhaustive search, as well as compare computed dipole moments to experimental values, based on molecular ensembles and their Boltzmann distributions. Results were also compared to the alternative molecular dynamics and simulated annealing methods. Our results clarify that many low energy conformations may be required to reproduce molecular properties, while single low energy conformers or ensembles of low energy conformers cannot account for the experimental properties of flexible molecules. Whereas ISE well reproduces conformations that are not separated by very large energy barriers, it has not been successful in reproducing conformations of strained molecules.     

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.     

Role of topology, nonadditivity, and water-mediated interactions in predicting the structures of alpha/beta proteins

The folding of alpha/beta proteins involves most of the commonly known structural and dynamic complexities of the protein energy landscapes. Thus, the interplay among different structural components, taking into account the cooperative interactions, is important in determining the success of protein structure prediction. In this work we present further developments of our knowledge-based force field for alpha/beta proteins, introducing more realistic modeling of many-body interactions governing the folding of beta-sheets. The model's innovations highlight both specific topological characteristics of secondary structures and the generic nonadditive interactions that are mediated by water. We also investigate how a coarse biasing of the protein morphology can be used to understand the role of heterogeneity in protein collapse. Analysis of the simulation results for three test alpha/beta proteins indicates that the addition of the topological and many-body ingredients to the model helps to greatly reduce the roughness in the energy landscape. Consequently, high quality candidate structures for alpha/beta proteins can be generated from simulated annealing runs, using very modest amounts of computer time.