Discrete path sampling

Energy landscapes of a hairpin peptide including NMR chemical shift restraints

Methods recently introduced to improve the efficiency of protein structure prediction simulations by adding a restraint potential to a molecular mechanics force field introduce additional input parameters that can affect the performance. Here we investigate the changes in the energy landscape as the relative weight of the two contributions, force field and restraint potential, is systematically altered, for restraint functions constructed from calculated nuclear magnetic resonance chemical shifts. Benchmarking calculations were performed on a 12-residue peptide, tryptophan zipper 1, which features both secondary structure (a β-hairpin) and specific packing of tryptophan sidechains. Basin-hopping global optimization was performed to assess the efficiency with which lowest-energy structures are located, and the discrete path sampling approach was employed to survey the energy landscapes between unfolded and folded structures. We find that inclusion of the chemical shift restraints improves the efficiency of structure prediction because the energy landscape becomes more funnelled and the proportion of local minima classified as native increases. However, the funnelling nature of the landscape is reduced as the relative contribution of the chemical shift restraint potential is increased past an optimal value.

Coarse-Grained HiRE-RNA Model for ab Initio RNA Folding beyond Simple Molecules, Including Noncanonical and Multiple Base Pairings

HiRE-RNA is a coarse-grained model for RNA structure prediction and the dynamical study of RNA folding. Using a reduced set of particles and detailed interactions accounting for base-pairing and stacking, we show that noncanonical and multiple base interactions are necessary to capture the full physical behavior of complex RNAs. In this paper, we give a full account of the model and present results on the folding, stability, and free energy surfaces of 16 systems with 12 to 76 nucleotides of increasingly complex architectures, ranging from monomers to dimers, using a total of 850 μs of simulation time.

Intrinsically disordered energy landscapes

Analysis of an intrinsically disordered protein (IDP) reveals an underlying multifunnel structure for the energy landscape. We suggest that such 'intrinsically disordered' landscapes, with a number of very different competing low-energy structures, are likely to characterise IDPs, and provide a useful way to address their properties. In particular, IDPs are present in many cellular protein interaction networks, and several questions arise regarding how they bind to partners. Are conformations resembling the bound structure selected for binding, or does further folding occur on binding the partner in a induced-fit fashion? We focus on the p53 upregulated modulator of apoptosis (PUMA) protein, which adopts an α-helical conformation when bound to its partner, and is involved in the activation of apoptosis. Recent experimental evidence shows that folding is not necessary for binding, and supports an induced-fit mechanism. Using a variety of computational approaches we deduce the molecular mechanism behind the instability of the PUMA peptide as a helix in isolation. We find significant barriers between partially folded states and the helix. Our results show that the favoured conformations are molten-globule like, stabilised by charged and hydrophobic contacts, with structures resembling the bound state relatively unpopulated in equilibrium.

Exploring the potential energy landscape of the Thomson problem via Newton homotopies

Locating the stationary points of a real-valued multivariate potential energy function is an important problem in many areas of science. This task generally amounts to solving simultaneous nonlinear systems of equations. While there are several numerical methods that can find many or all stationary points, they each exhibit characteristic problems. Moreover, traditional methods tend to perform poorly near degenerate stationary points with additional zero Hessian eigenvalues. We propose an efficient and robust implementation of the Newton homotopy method, which is capable of quickly sampling a large number of stationary points of a wide range of indices, as well as degenerate stationary points. We demonstrate our approach by applying it to the Thomson problem. We also briefly discuss a possible connection between the present work and Smale's 7th problem.

Analysis of the Contrasting Pathogenicities Induced by the D222G Mutation in 1918 and 2009 Pandemic Influenza A Viruses

In 2009, the D222G mutation in the hemagglutinin (HA) glycoprotein of pandemic H1N1 influenza A virus was found to correlate with fatal and severe human infections. Previous static structural analysis suggested that, unlike the H1N1 viruses prevalent in 1918, the mutation did not compromise binding to human α2,6-linked glycan receptors, allowing it to transmit efficiently. Here we investigate the interconversion mechanism between two predicted binding modes in both 2009 and 1918 HAs, introducing a highly parallel intermediate network search scheme to construct kinetically relevant pathways efficiently. Accumulated mutations at positions 183 and 224 that alter the size of the binding pocket are identified with the fitness of the 2009 pandemic virus carrying the D222G mutation. This result suggests that the pandemic H1N1 viruses could gain binding affinity to the α2,3-linked glycan receptors in the lungs, usually associated with highly pathogenic avian influenza, without compromising viability.

Communication: optimal parameters for basin-hopping global optimization based on Tsallis statistics

A fundamental problem associated with global optimization is the large free energy barrier for the corresponding solid-solid phase transitions for systems with multi-funnel energy landscapes. To address this issue we consider the Tsallis weight instead of the Boltzmann weight to define the acceptance ratio for basin-hopping global optimization. Benchmarks for atomic clusters show that using the optimal Tsallis weight can improve the efficiency by roughly a factor of two. We present a theory that connects the optimal parameters for the Tsallis weighting, and demonstrate that the predictions are verified for each of the test cases.

Roles of dynamical symmetry breaking in driving oblate-prolate transitions of atomic clusters

This paper explores the driving mechanisms for structural transitions of atomic clusters between oblate and prolate isomers. We employ the hyperspherical coordinates to investigate structural dynamics of a seven-atom cluster at a coarse-grained level in terms of the dynamics of three gyration radii and three principal axes, which characterize overall mass distributions of the cluster. Dynamics of gyration radii is governed by two kinds of forces. One is the potential force originating from the interactions between atoms. The other is the dynamical forces called the internal centrifugal forces, which originate from twisting and shearing motions of the system. The internal centrifugal force arising from twisting motions has an effect of breaking the symmetry between two gyration radii. As a result, in an oblate isomer, activation of the internal centrifugal force that has the effect of breaking the symmetry between the two largest gyration radii is crucial in triggering structural transitions into prolate isomers. In a prolate isomer, on the other hand, activation of the internal centrifugal force that has the effect of breaking the symmetry between the two smallest gyration radii is crucial in triggering structural transitions into oblate isomers. Activation of a twisting motion that switches the movement patterns of three principal axes is also important for the onset of structural transitions between oblate and prolate isomers. Based on these trigger mechanisms, we finally show that selective activations of specific gyration radii and twisting motions, depending on the isomer of the cluster, can effectively induce structural transitions of the cluster. The results presented here could provide further insights into the control of molecular reactions.

Locating landmarks on high-dimensional free energy surfaces

Coarse graining of complex systems possessing many degrees of freedom can often be a useful approach for analyzing and understanding key features of these systems in terms of just a few variables. The relevant energy landscape in a coarse-grained description is the free energy surface as a function of the coarse-grained variables, which, despite the dimensional reduction, can still be an object of high dimension. Consequently, navigating and exploring this high-dimensional free energy surface is a nontrivial task. In this paper, we use techniques from multiscale modeling, stochastic optimization, and machine learning to devise a strategy for locating minima and saddle points (termed "landmarks") on a high-dimensional free energy surface "on the fly" and without requiring prior knowledge of or an explicit form for the surface. In addition, we propose a compact graph representation of the landmarks and connections between them, and we show that the graph nodes can be subsequently analyzed and clustered based on key attributes that elucidate important properties of the system. Finally, we show that knowledge of landmark locations allows for the efficient determination of their relative free energies via enhanced sampling techniques.

The role of binding site on the mechanical unfolding mechanism of ubiquitin

We apply novel atomistic simulations based on potential energy surface exploration to investigate the constant force-induced unfolding of ubiquitin. At the experimentally-studied force clamping level of 100 pN, we find a new unfolding mechanism starting with the detachment between β₅ and β₃ involving the binding site of ubiquitin, the Ile44 residue. This new unfolding pathway leads to the discovery of new intermediate configurations, which correspond to the end-to-end extensions previously seen experimentally. More importantly, it demonstrates the novel finding that the binding site of ubiquitin can be responsible not only for its biological functions, but also its unfolding dynamics. We also report in contrast to previous single molecule constant force experiments that when the clamping force becomes smaller than about 300 pN, the number of intermediate configurations increases dramatically, where almost all unfolding events at 100 pN involve an intermediate configuration. By directly calculating the life times of the intermediate configurations from the height of the barriers that were crossed on the potential energy surface, we demonstrate that these intermediate states were likely not observed experimentally due to their lifetimes typically being about two orders of magnitude smaller than the experimental temporal resolution.

Reaction sampling and reactivity prediction using the stochastic surface walking method

The new theoretical method demonstrates the ability of automated reaction sampling and activity prediction for complex organic reactions.

Metastability, spectrum, and eigencurrents of the Lennard-Jones-38 network

We develop computational tools for spectral analysis of stochastic networks representing energy landscapes of atomic and molecular clusters. Physical meaning and some properties of eigenvalues, left and right eigenvectors, and eigencurrents are discussed. We propose an approach to compute a collection of eigenpairs and corresponding eigencurrents describing the most important relaxation processes taking place in the system on its way to the equilibrium. It is suitable for large and complex stochastic networks where pairwise transition rates, given by the Arrhenius law, vary by orders of magnitude. The proposed methodology is applied to the network representing the Lennard-Jones-38 cluster created by Wales's group. Its energy landscape has a double funnel structure with a deep and narrow face-centered cubic funnel and a shallower and wider icosahedral funnel. However, the complete spectrum of the generator matrix of the Lennard-Jones-38 network has no appreciable spectral gap separating the eigenvalue corresponding to the escape from the icosahedral funnel. We provide a detailed description of the escape process from the icosahedral funnel using the eigencurrent and demonstrate a superexponential growth of the corresponding eigenvalue. The proposed spectral approach is compared to the methodology of the Transition Path Theory. Finally, we discuss whether the Lennard-Jones-38 cluster is metastable from the points of view of a mathematician and a chemical physicist, and make a connection with experimental works.

Exploring energy landscapes: from molecular to mesoscopic systems

We review a comprehensive computational framework to survey the potential energy landscape for systems composed of rigid or partially rigid molecules. Illustrative case studies relevant to a wide range of molecular clusters and soft and condensed matter systems are discussed.

Energy landscapes of planar colloidal clusters

A short-ranged pairwise Morse potential is used to model colloidal clusters with planar morphologies. Potential and free energy global minima as well as rearrangement paths, obtained by basin-hopping global optimisation and discrete path sampling, are characterised. The potential and free energy landscapes are visualised using disconnectivity graphs. The short-ranged potential is found to favour close-packed structures, with the potential energy primarily controlled by the number of nearest neighbour contacts. In the case of quasi-degeneracy the free energy global minimum may differ from the potential energy global minimum. This difference is due to symmetry effects, which result in a higher entropy for structures with lower symmetry.

Design principles for Bernal spirals and helices with tunable pitch

Using the framework of potential energy landscape theory, we describe two in silico designs for self-assembling helical colloidal superstructures based upon dipolar dumbbells and Janus-type building blocks, respectively. Helical superstructures with controllable pitch length are obtained using external magnetic field driven assembly of asymmetric dumbbells involving screened electrostatic as well as magnetic dipolar interactions. The pitch of the helix is tuned by modulating the Debye screening length over an experimentally accessible range. The second design is based on building blocks composed of rigidly linked spheres with short-range anisotropic interactions, which are predicted to self-assemble into Bernal spirals. These spirals are quite flexible, and longer helices undergo rearrangements via cooperative, hinge-like moves, in agreement with experiment.

Structure and properties of DNA in apolar solvents

The study of nucleic acids in low-polarity environments paves the way for novel biotechnological applications of DNA. Here, we use a repertoire of atomistic molecular simulation tools to study the nature of DNA when placed in a highly apolar environment and when transferred from aqueous to apolar solvent. Our results show that DNA becomes stiffer in apolar solvents and suggest that highly negatively charged states, which are the most prevalent in water, are strongly disfavored in apolar solvents and neutral states with conformations not far from the aqueous ones are the dominant forms. Transfer from water to an apolar solvent such as CCl₄ is unlikely to occur, but our results suggest that if forced, the DNA would migrate surrounded by a small shell of water (the higher the DNA charge, the larger the number of water molecules in this shell). Even the neutral form (predicted to be the dominant one in apolar solvents) would surround itself by a small number of highly stable water molecules when moved from water to a highly apolar environment. Neutralization of DNA charges seems a crucial requirement for transfer of DNA to apolar media, and the most likely mechanism to achieve good transfer properties.

Visualization of protein folding funnels in lattice models

Protein folding occurs in a very high dimensional phase space with an exponentially large number of states, and according to the energy landscape theory it exhibits a topology resembling a funnel. In this statistical approach, the folding mechanism is unveiled by describing the local minima in an effective one-dimensional representation. Other approaches based on potential energy landscapes address the hierarchical structure of local energy minima through disconnectivity graphs. In this paper, we introduce a metric to describe the distance between any two conformations, which also allows us to go beyond the one-dimensional representation and visualize the folding funnel in 2D and 3D. In this way it is possible to assess the folding process in detail, e.g., by identifying the connectivity between conformations and establishing the paths to reach the native state, in addition to regions where trapping may occur. Unlike the disconnectivity maps method, which is based on the kinetic connections between states, our methodology is based on structural similarities inferred from the new metric. The method was developed in a 27-mer protein lattice model, folded into a 3×3×3 cube. Five sequences were studied and distinct funnels were generated in an analysis restricted to conformations from the transition-state to the native configuration. Consistent with the expected results from the energy landscape theory, folding routes can be visualized to probe different regions of the phase space, as well as determine the difficulty in folding of the distinct sequences. Changes in the landscape due to mutations were visualized, with the comparison between wild and mutated local minima in a single map, which serves to identify different trapping regions. The extension of this approach to more realistic models and its use in combination with other approaches are discussed.

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.

An inversion-relaxation approach for sampling stationary points of spin model Hamiltonians

Sampling the stationary points of a complicated potential energy landscape is a challenging problem. Here, we introduce a sampling method based on relaxation from stationary points of the highest index of the Hessian matrix. We illustrate how this approach can find all the stationary points for potentials or Hamiltonians bounded from above, which includes a large class of important spin models, and we show that it is far more efficient than previous methods. For potentials unbounded from above, the relaxation part of the method is still efficient in finding minima and transition states, which are usually the primary focus of attention for atomistic systems.

Proton transfer pathways, energy landscape, and kinetics in creatine-water systems

We study the exchange processes of the metabolite creatine, which is present in both tumorous and normal tissues and has NH2 and NH groups that can transfer protons to water. Creatine produces chemical exchange saturation transfer (CEST) contrast in magnetic resonance imaging (MRI). The proton transfer pathway from zwitterionic creatine to water is examined using a kinetic transition network constructed from the discrete path sampling approach and an approximate quantum-chemical energy function, employing the self-consistent-charge density-functional tight-binding (SCC-DFTB) method. The resulting potential energy surface is visualized by constructing disconnectivity graphs. The energy landscape consists of two distinct regions corresponding to the zwitterionic creatine structures and deprotonated creatine. The activation energy that characterizes the proton transfer from the creatine NH2 group to water was determined from an Arrhenius fit of rate constants as a function of temperature, obtained from harmonic transition state theory. The result is in reasonable agreement with values obtained in water exchange spectroscopy (WEX) experiments.

Stochastic surface walking method for crystal structure and phase transition pathway prediction

SSW-crystal method for automated structure search and phase transition pathway sampling of crystals.

Observation time scale, free-energy landscapes, and molecular symmetry

When structures that interconvert on a given time scale are lumped together, the corresponding free-energy surface becomes a function of the observation time. This view is equivalent to grouping structures that are connected by free-energy barriers below a certain threshold. We illustrate this time dependence for some benchmark systems, namely atomic clusters and alanine dipeptide, highlighting the connections to broken ergodicity, local equilibrium, and "feasible" symmetry operations of the molecular Hamiltonian.

Exploring Energy Landscapes: Metrics, Pathways, and Normal-Mode Analysis for Rigid-Body Molecules

We present new methodology for exploring the energy landscapes of molecular systems, using angle-axis variables for the rigid-body rotational coordinates. The key ingredient is a distance measure or metric tensor, which is invariant to global translation and rotation. The metric is used to formulate a generalized nudged elastic band method for calculating pathways, and a full prescription for normal-mode analysis is described. The methodology is tested by mapping the potential energy and free energy landscape of the water octamer, described by the TIP4P potential.

Temporal disconnectivity of the energy landscape in glassy systems

An alternative graphical representation of the potential energy landscape (PEL) has been developed and applied to a binary Lennard-Jones glassy system, providing insight into the unique topology of the system's potential energy hypersurface. With the help of this representation one is able to monitor the different explored basins of the PEL, as well as how--and mainly when--subsets of basins communicate with each other via transitions in such a way that details of the prior temporal history have been erased, i.e., local equilibration between the basins in each subset has been achieved. In this way, apart from detailed information about the structure of the PEL, the system's temporal evolution on the PEL is described. In order to gather all necessary information about the identities of two or more basins that are connected with each other, we consider two different approaches. The first one is based on consideration of the time needed for two basins to mutually equilibrate their populations according to the transition rate between them, in the absence of any effect induced by the rest of the landscape. The second approach is based on an analytical solution of the master equation that explicitly takes into account the entire explored landscape. It is shown that both approaches lead to the same result concerning the topology of the PEL and dynamical evolution on it. Moreover, a "temporal disconnectivity graph" is introduced to represent a lumped system stemming from the initial one. The lumped system is obtained via a specially designed algorithm [N. Lempesis, D. G. Tsalikis, G. C. Boulougouris, and D. N. Theodorou, J. Chem. Phys. 135, 204507 (2011)]. The temporal disconnectivity graph provides useful information about both the lumped and the initial systems, including the definition of "metabasins" as collections of basins that communicate with each other via transitions that are fast relative to the observation time. Finally, the two examined approaches are compared to an "on the fly" molecular dynamics-based algorithm [D. G. Tsalikis, N. Lempesis, G. C. Boulougouris, and D. N. Theodorou, J. Chem. Theory Comput. 6, 1307 (2010)].

Exploring the Energy Landscapes of Cyclic Tetrapeptides with Discrete Path Sampling

Cyclic tetrapeptides are an important class of biologically active molecules that exhibit interesting conformational dynamics, with slow interconversion of several different structures. We present calculations on their energy landscapes using discrete path sampling. In acyclic peptides and large cyclic peptides, isomers containing cis-peptide groups are much less stable than the all-trans isomers and separated from them by large barriers. Strain in small cyclic peptides causes the cis and trans isomers to be closer in energy and separated by much lower barriers. If d-amino acids or proline residues are introduced, isomers containing cis-peptides become more stable than the all-trans structures. We also show that changing the polarity of the solvent has a significant effect on the energy landscapes of cyclic tetrapeptides, causing changes in the orientations of the peptide groups and in the degree of intramolecular hydrogen bonding.

Quasi-Continuous Interpolation Scheme for Pathways between Distant Configurations

A quasi-continuous interpolation (QCI) scheme is introduced for characterizing physically realistic initial pathways from which to initiate transition state searches and construct kinetic transition networks. Applications are presented for peptides, proteins, and a morphological transformation in an atomic cluster. The first step in each case involves end point alignment, and we describe the use of a shortest augmenting path algorithm for optimizing permutational isomers. The QCI procedure then employs an interpolating potential, which preserves the covalent bonding framework for the biomolecules and includes repulsive terms between unconstrained atoms. This potential is used to identify an interpolating path by minimizing contributions from a connected set of images, including terms corresponding to minima in the interatomic distances between them. This procedure detects unphysical geometries in the line segments between images. The most difficult cases, where linear interpolation would involve chain crossings, are treated by growing the structure an atom at a time using the interpolating potential. To test the QCI procedure, we carry through a series of benchmark calculations where the initial interpolation is coupled to explicit transition state searches to produce complete pathways between specified local minima.

Force fluctuations impact kinetics of biomolecular systems

A wide array of biological processes occur at rates that vary significantly with force. Instantaneous molecular forces fluctuate due to thermal noise and active processes, leading to concomitant fluctuations in biomolecular rate constants. We demonstrate that such fluctuations have a dramatic effect on the transition kinetics of force-dependent processes. As an illustrative, biologically relevant example, we model the pausing of eukaryotic RNA polymerase as it transcribes nucleosomal DNA. Incorporating force fluctuations in the model yields qualitatively different predictions for the pausing time scales when compared to behavior under the average force alone. We use our model to illustrate the broad range of behaviors that can arise in biomolecular processes that are susceptible to force fluctuations. The fluctuation time scale, which varies significantly for in vivo biomolecular processes, yields very different results for overall rates and dramatically alters the force regime of relevance to the transition. Our results emphasize the importance of transient high-force behavior for determining kinetics in the fluctuating environment of a living cell.