Dependency of ligand free energy landscapes on charge parameters and solvent models

We performed replica-exchange molecular dynamics (REMD) simulations of six ligands to examine the dependency of their free energy landscapes on charge parameters and solvent models. Six different charge parameter sets for each ligand were first generated by RESP and AM1-BCC methods using three different conformations independently. RESP charges showed some conformational dependency. On the other hand, AM1-BCC charges did not show conformational dependency and well reproduced the overall trend of RESP charges. The free energy landscapes obtained from the REMD simulations of ligands in vacuum, Generalized-Born (GB), and TIP3P solutions were then analyzed. We found that even small charge differences can produce qualitatively different landscapes in vacuum condition, but the differences tend to be much smaller under GB and TIP3P conditions. The simulations in the GB model well reproduced the landscapes in the TIP3P model using only a fraction of the computational cost. The protein-bound ligand conformations were rarely the global minimum states, but similar conformations were found to exist in aqueous solution without proteins in regions close to the global minimum, local minimum or intermediate states.     

Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics

Recently, we have introduced a new method, metadynamics, which is able to sample rarely occurring transitions and to reconstruct the free energy as a function of several variables with a controlled accuracy. This method has been successfully applied in many different fields, ranging from chemistry to biophysics and ligand docking and from material science to crystal structure prediction. We present an important development that speeds up metadynamics calculations by orders of magnitude and renders the algorithm much more robust. We use multiple interacting simulations, walkers, for exploring and reconstructing the same free energy surface. Each walker contributes to the history-dependent potential that, in metadynamics, is an estimate of the free energy. We show that the error on the reconstructed free energy does not depend on the number of walkers, leading to a fully linear scaling algorithm even on inexpensive loosely coupled clusters of PCs. In addition, we show that the accuracy and stability of the method are much improved by combining it with a weighted histogram analysis. We check the validity of our new method on a realistic application.     

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     

Adlayer morphologies and free energy landscapes of clusters of bis-fullerenes on model gold surfaces

There have been a few experimental reports of self-assembled adlayers of bis-fullerene molecules on solid substrates. Most of these studies suggest the adsorbate molecules are lying down on the surface, with the fullerene moieties almost close packed. However, very little theoretical work has been carried out on such systems, and little is known about the roles played by different parts of the potential energy in driving the self-assembly. We carry out a Temperature Replica Exchange Monte Carlo study here of two representative bis-fullerene molecules on a metal substrate. We use a coarse-grained model potential energy function, in which certain parameters can be varied within the range of their experimental uncertainty. The molecules investigated consist of two fullerene moieties bonded by a rigid bridging group. In particular, the effect of the strength of the fullerene interaction E(FG) with the substrate (nominally Au(111)) has been investigated in detail. To ensure efficient sampling of the rugged potential energy surfaces encountered in the simulations, we utilize replica exchange techniques. These enable us to construct free energy landscapes for the system. We find that for relatively low values of E(FG) the molecules form standing-up adlayers. By contrast, for higher values of E(FG), lying-down adlayers dominate. For one molecule, two different crystalline adlayer morphologies have been identified. The detailed structure of the lying-down layer is a function of the temperature and of the group used to bridge the fullerene moieties.     

Long-term observation of fluorescence of free single molecules to explore protein-folding energy landscapes

A method was developed to detect fluorescence intensity signals from single molecules diffusing freely in a capillary cell. A unique optical system based on a spherical mirror was designed to enable quantitative detection of the fluorescence intensity. Furthermore, "flow-and-stop" control of the sample can extend the observation time of single molecules to several seconds, which is more than 1000 times longer than the observation time available using a typical confocal method. We used this method to scrutinize the fluorescence time series of the labeled cytochrome c in the unfolded state. Time series analyses of the trajectories based on local equilibrium state analysis revealed dynamically differing substates on a millisecond time scale. This system presents a new avenue for experimental characterization of the protein-folding energy landscape.     

Graph‐theoretical identification of dissociation pathways on free energy landscapes of biomolecular interaction

Biomolecular association and dissociation reactions take place on complicated interaction free energy landscapes that are still very hard to characterize computationally. For large enough distances, though, it often suffices to consider the six relative translational and rotational degrees of freedom of the two particles treated as rigid bodies. Here, we computed the six‐dimensional free energy surface of a dimer of water‐soluble alpha‐helices by scanning these six degrees of freedom in about one million grid points. In each point, the relative free energy difference was computed as the sum of the polar and nonpolar solvation free energies of the helix dimer and of the intermolecular coulombic interaction energy. The Dijkstra graph algorithm was then applied to search for the lowest cost dissociation pathways based on a weighted, directed graph, where the vertices represent the grid points, the edges connect the grid points and their neighbors, and the weights are the reaction costs between adjacent pairs of grid points. As an example, the configuration of the bound state was chosen as the source node, and the eight corners of the translational cube were chosen as the destination nodes. With the strong electrostatic interaction of the two helices giving rise to a clearly funnel‐shaped energy landscape, the eight lowest‐energy cost pathways coming from different orientations converge into a well‐defined pathway for association. We believe that the methodology presented here will prove useful for identifying low‐energy association and dissociation pathways in future studies of complicated free energy landscapes for biomolecular interaction. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Archetypal energy landscapes: dynamical diagnosis

Recent studies have identified several motifs for potential energy surfaces corresponding to distinct dynamic and thermodynamic properties. The corresponding disconnectivity graphs were identified as "palm tree," "willow tree," and "banyan tree" patterns. In the present contribution we present a quantitative analysis of the relation between the topography and dynamics for each of these motifs. For the palm tree and willow tree forms we find that the arrangement of the stationary points in the monotonic sequences with respect to the global minimum is the most important factor in establishing the kinetic properties. However, the results are somewhat different for motifs involving a rough surface with several deep basins (banyan tree motif), with large barriers relative to the energy differences between minima. Here it is the size of the barrier for escape from the region relative to the barriers at the bottom that is most important. The present results may be helpful in distinguishing between the dynamics of "structure seeking" and "glass forming" systems.     

Convergence of folding free energy landscapes via application of enhanced sampling methods in a distributed computing environment

We have implemented the serial replica exchange method (SREM) and simulated tempering (ST) enhanced sampling algorithms in a global distributed computing environment. Here we examine the helix-coil transition of a 21 residue alpha-helical peptide in explicit solvent. For ST, we demonstrate the efficacy of a new method for determining initial weights allowing the system to perform a random walk in temperature space based on short trial simulations. These weights are updated throughout the production simulation by an adaptive weighting method. We give a detailed comparison of SREM, ST, as well as standard MD and find that SREM and ST give equivalent results in reasonable agreement with experimental data. In addition, we find that both enhanced sampling methods are much more efficient than standard MD simulations. The melting temperature of the Fs peptide with the AMBER99phi potential was calculated to be about 310 K, which is in reasonable agreement with the experimental value of 334 K. We also discuss other temperature dependent properties of the helix-coil transition. Although ST has certain advantages over SREM, both SREM and ST are shown to be powerful methods via distributed computing and will be applied extensively in future studies of complex bimolecular systems.     

Exploring protein energy landscapes with hierarchical clustering

In this work we present a new method for investigating local energy minima on a protein energy landscape. The CABS (CAlpha, CBeta and the center of mass of the Side chain) method was employed for generating protein models, but any other method could be used instead. Cα traces from an ensemble of models are hierarchical clustered with the HCPM (Hierarchical Clustering of Protein Models) method. The efficiency of this method for sampling and analyzing energy landscapes is shown.     

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.     

Intrinsic energy landscapes of amino acid side-chains

Amino acid side-chain conformational properties influence the overall structural and dynamic properties of proteins and, therefore, their biological functions. In this study, quantum mechanical (QM) potential energy surfaces for the rotation of side-chain χ(1) and χ(2) torsions in dipeptides in the alphaR, beta, and alphaL backbone conformations were calculated. The QM energy surfaces provide a broad view of the intrinsic conformational properties of each amino acid side-chain. The extent to which intrinsic energetics dictates side-chain orientation was studied through comparisons of the QM energy surfaces with χ(1) and χ(2) free energy surfaces from probability distributions obtained from a survey of high resolution crystal structures. In general, the survey probability maxima are centered in minima of the QM surfaces as expected for sp(3) (or sp(2) for χ(2) of Asn, Phe, Trp, and Tyr) atom centers with strong variations between amino acids occurring in the energies of the minima indicating intrinsic differences in rotamer preferences. High correlations between the QM and survey data were found for hydrophobic side-chains except Met, suggesting minimal influence of the protein and solution environments on their conformational distributions. Conversely, low correlations for polar or charged side-chains indicate a dominant role of the environment in stabilizing conformations that are not intrinsically favored. Data also link the presence of off-rotamers in His and Trp to favorable interactions with the backbone. Results also suggest that the intrinsic energetics of the side-chains of Phe and Tyr may play important roles in protein folding and stability. Analyses on whether intrinsic side-chain energetics can influence backbone preference identified a strong correlation for residues in the alphaL backbone conformation. It is suggested that this correlation reflects the intrinsic instability of the alphaL backbone such that assumption of this backbone conformation is facilitated by intrinsically favorable side-chain conformations. Together our results offer a broad overview of the conformational properties of amino acid side-chains and the QM data may be used as target data for force field optimization.     

Analyzing energy landscapes for folding model proteins

A new benchmark 20-bead HP model protein sequence (on a square lattice), which has 17 distinct but degenerate global minimum (GM) energy structures, has been studied using a genetic algorithm (GA). The relative probabilities of finding particular GM conformations are determined and related to the theoretical probability of generating these structures using a recoil growth constructor operator. It is found that for longer successful GA runs, the GM probability distribution is generally very different from the constructor probability, as other GA operators have had time to overcome any initial bias in the originally generated population of structures. Structural and metric relationships (e.g., Hamming distances) between the 17 distinct GM are investigated and used, in conjunction with data on the connectivities of the GM and the pathways that link them, to explain the GM probability distributions obtained by the GA. A comparison is made of searches where the sequence is defined in the normal (forward) and reverse directions. The ease of finding mirror image solutions are also compared. Finally, this approach is applied to rationalize the ease or difficulty of finding the GM for a number of standard benchmark HP sequences on the square lattice. It is shown that the relative probabilities of finding particular members of a set of degenerate global minima depend critically on the topography of the energy landscape in the vicinity of the GM, the connections and distances between the GM, and the nature of the operators used in the chosen search method.     

Conformational energy landscapes of liquid crystal molecules

The conformational energy landscape of the prototypical nematic liquid crystal 4-pentyl-4cyanobiphenyl (5CB) is studied using first principles computer modelling. It is found that the most favourable conformation occurs when the two constituent phenyl rings are inclined at an angle of 31 with respect to each other. Also, the orientation of the alkyl chain is found to have an important influence on the ring-ring torsional potential. We fit the energy surface of these coupled torsions to yield an accurate intramolecular potential for use in empirical modelling. To test the strength of the coupling between the alkyl tail and the phenyl rings and the cyano group, we also calculate potentials for the relative orientation of the phenyl rings in biphenyl and cyanobiphenyl (0CB). Our calculations are performed using density functional theory using pseudo-potentials and the generalized gradient approximation to exchange and correlation. The molecular electronic wavefunction is expanded in terms of a plane wave basis set. We compare our results with recent NMR and Gaussian-based quantum chemistry calculations where available.     

Global optimization and the energy landscapes of Dzugutov clusters

The global minima of clusters bound by a Dzugutov potential form non-compact polytetrahedral clusters mainly composed of interpenetrating and face-sharing 13-atom icosahedra. As the size increases, these icosahedral units first form linear arrays, then two-dimensional rings, then three-dimensional networks. Characterization of the energy landscapes of these clusters shows that they are particularly rough and generally exhibit a multiple-funnel topography. These results provide new insights into the structure and dynamics of bulk supercooled Dzugutov liquids and the form of the bulk phase diagram.     

Conformational ensembles and sampled energy landscapes: Analysis and comparison

We present novel algorithms and software addressing four core problems in computational structural biology, namely analyzing a conformational ensemble, comparing two conformational ensembles, analyzing a sampled energy landscape, and comparing two sampled energy landscapes. Using recent developments in computational topology, graph theory, and combinatorial optimization, we make two notable contributions. First, we present a generic algorithm analyzing height fields. We then use this algorithm to perform density‐based clustering of conformations, and to analyze a sampled energy landscape in terms of basins and transitions between them. In both cases, topological persistence is used to manage (geometric) frustration. Second, we introduce two algorithms to compare transition graphs. The first is the classical earth mover distance metric which depends only on local minimum energy configurations along with their statistical weights, while the second incorporates topological constraints inherent to conformational transitions. Illustrations are provided on a simplified protein model (BLN69), whose frustrated potential energy landscape has been thoroughly studied. The software implementing our tools is also made available, and should prove valuable wherever conformational ensembles and energy landscapes are used. © 2015 Wiley Periodicals, Inc.     

Consistent free energy landscapes and thermodynamic properties of small proteins based on a single all-atom force field employing an implicit solvation

We have attempted to improve the PARAM99 force field in conjunction with the generalized Born (GB) solvation model with a surface area correction for more consistent protein folding simulations. For this purpose, using an extended alphabeta training set of five well-studied molecules with various folds (alpha, beta, and betabetaalpha), a previously modified version of PARAM99/GBSA is further refined, such that all native states of the five training species correspond to their lowest free energy minimum states. The resulting modified force field (PARAM99MOD5/GBSA) clearly produces reasonably acceptable conformational free energy surfaces of the training set with correct identifications of their native states in the free energy minimum states. Moreover, due to its well-balanced nature, this new force field is expected to describe secondary structure propensities of diverse folds in a more consistent manner. Remarkably, temperature dependent behaviors simulated with the current force field are in good agreement with the experiment. This agreement is a significant improvement over the existing standard all-atom force fields. In addition, fundamentally important thermodynamic quantities, such as folding enthalpy (DeltaH) and entropy (DeltaS), agree reasonably well with the experimental data.     

Randomized tree construction algorithm to explore energy landscapes

In this work, a new method for exploring conformational energy landscapes is described. The method, called transition-rapidly exploring random tree (T-RRT), combines ideas from statistical physics and robot path planning algorithms. A search tree is constructed on the conformational space starting from a given state. The tree expansion is driven by a double strategy: on the one hand, it is naturally biased toward yet unexplored regions of the space; on the other, a Monte Carlo-like transition test guides the expansion toward energetically favorable regions. The balance between these two strategies is automatically achieved due to a self-tuning mechanism. The method is able to efficiently find both energy minima and transition paths between them. As a proof of concept, the method is applied to two academic benchmarks and the alanine dipeptide.     

Imaging energy landscapes with concentrated diffusing colloidal probes

The ability to locally interrogate interactions between particles and energetically patterned surfaces provides essential information to design, control, and optimize template directed self-assembly processes. Although numerous techniques are capable of characterizing local physicochemical surface properties, no current method resolves interactions between colloids and patterned surfaces on the order of the thermal energy kT, which is the inherent energy scale of equilibrium self-assembly processes. Here, the authors describe video microscopy measurements and an inverse Monte Carlo analysis of diffusing colloidal probes as a means to image three dimensional free energy and potential energy landscapes due to physically patterned surfaces. In addition, they also develop a consistent analysis of self-diffusion in inhomogeneous fluids of concentrated diffusing probes on energy landscapes, which is important to the temporal imaging process and to self-assembly kinetics. Extension of the concepts developed in this work suggests a general strategy to image multidimensional and multiscale physical, chemical, and biological surfaces using a variety of diffusing probes (i.e., molecules, macromolecules, nanoparticles, and colloids).