Energy landscapes and persistent minima.

We consider a coarse-graining of high-dimensional potential energy landscapes based upon persistences, which correspond to lowest barrier heights to lower-energy minima. Persistences can be calculated efficiently for local minima in kinetic transition networks that are based on stationary points of the prevailing energy landscape. The networks studied here represent peptides, proteins, nucleic acids, an atomic cluster, and a glassy system. Minima with high persistence values are likely to represent some form of alternative structural morphology, which, if appreciably populated at the prevailing temperature, could compete with the global minimum (defined as infinitely persistent). Threshold values on persistences (and in some cases equilibrium occupation probabilities) have therefore been used in this work to select subsets of minima, which were then analysed to see how well they can represent features of the full network. Simplified disconnectivity graphs showing only the selected minima can convey the funnelling (including any multiple-funnel) characteristics of the corresponding full graphs. The effect of the choice of persistence threshold on the reduced disconnectivity graphs was considered for a system with a hierarchical, glassy landscape. Sets of persistent minima were also found to be useful in comparing networks for the same system sampled under different conditions, using minimum oriented spanning forests.

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.

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.

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.

DL-FIND: an open-source geometry optimizer for atomistic simulations.

Geometry optimization, including searching for transition states, accounts for most of the CPU time spent in quantum chemistry, computational surface science, and solid-state physics, and also plays an important role in simulations employing classical force fields. We have implemented a geometry optimizer, called DL-FIND, to be included in atomistic simulation codes. It can optimize structures in Cartesian coordinates, redundant internal coordinates, hybrid-delocalized internal coordinates, and also functions of more variables independent of atomic structures. The implementation of the optimization algorithms is independent of the coordinate transformation used. Steepest descent, conjugate gradient, quasi-Newton, and L-BFGS algorithms as well as damped molecular dynamics are available as minimization methods. The partitioned rational function optimization algorithm, a modified version of the dimer method and the nudged elastic band approach provide capabilities for transition-state search. Penalty function, gradient projection, and Lagrange-Newton methods are implemented for conical intersection optimizations. Various stochastic search methods, including a genetic algorithm, are available for global or local minimization and can be run as parallel algorithms. The code is released under the open-source GNU LGPL license. Some selected applications of DL-FIND are surveyed.

Refined kinetic transition networks for the GB1 hairpin peptide.

Refinement of databases of connected stationary points to describe global kinetics is discussed for the GB1 hairpin peptide modelled by an empirical potential and an implicit solvent model. Two approaches to the removal of artificial kinetic frustration caused by undersampling are separately applied to an initial database of stationary points. We consider both additional sampling between minima close in energy but separated by high barriers, and the removal of stationary points that do not contribute significantly to the calculated rate constants for the initial database. Results from these two approaches are found to be consistent: the transition networks produced in both cases exhibit structure-seeking properties because most of the initial frustration is removed. Excluding stationary points from the initial database that do not appear on kinetically relevant paths proves to be much less computationally expensive than subsequently finding better connections for them. After application of a coarse-graining scheme that groups together sets of minima separated by low barriers, the calculated folding time is consistent with expectations for beta-hairpins modelled using implicit solvent. The folding mechanism corresponding to the most significant kinetic paths involves early compaction, followed by formation of the turn and then completion of the hydrophobic core.

Folding pathways and rates for the three-stranded beta-sheet peptide Beta3s using discrete path sampling.

The discrete path sampling method was used to investigate the folding of a three-stranded antiparallel beta-sheet peptide, Beta3s, described by an empirical potential and implicit solvent model. After application of a coarse-graining scheme that groups together sets of minima in local equilibrium, the calculated folding time was in reasonable agreement with other simulations and consistent with the experimental upper bound. The folding mechanism exhibited by the most significant discrete paths involves early formation of the C-terminal hairpin followed by docking of the N-terminal strand.

Finding pathways between distant local minima.

We report a new algorithm for constructing pathways between local minima that involve a large number of intervening transition states on the potential energy surface. A significant improvement in efficiency has been achieved by changing the strategy for choosing successive pairs of local minima that serve as endpoints for the next search. We employ Dijkstra's algorithm [E. W. Dijkstra, Numer. Math. 1, 269 (1959)] to identify the "shortest" path corresponding to missing connections within an evolving database of local minima and the transition states that connect them. The metric employed to determine the shortest missing connection is a function of the minimized Euclidean distance. We present applications to the formation of buckminsterfullerene and to the folding of various biomolecules: the B1 domain of protein G, tryptophan zippers, and the villin headpiece subdomain. The corresponding pathways contain up to 163 transition states and will be used in future discrete path sampling calculations.

Global optimization and folding pathways of selected alpha-helical proteins.

The results of basin-hopping global optimization simulations are presented for four small, alpha-helical proteins described by a coarse-grained potential. A step-taking scheme that incorporates the local conformational preferences extracted from a large number of high-resolution protein structures is compared with an unbiased scheme. In addition, the discrete path sampling method is used to investigate the folding of one of the proteins, namely, the villin headpiece subdomain. Folding times from kinetic Monte Carlo simulations and iterative calculations based on a Markovian first-step analysis for the resulting stationary-point database are in good mutual agreement, but differ significantly from the experimental values, probably because the native state is not the global free energy minimum for the potential employed.