Rough Free Energy Landscape Research Articles

Resampling schemes in population annealing: Numerical and theoretical results.

The population annealing algorithm is a population-based equilibrium version of simulated annealing. It can sample thermodynamic systems with rough free-energy landscapes more efficiently than standard Markov chain Monte Carlo alone. A number of parameters can be fine-tuned to improve the performance of the population annealing algorithm. While there is some numerical and theoretical work on most of these parameters, there appears to be a gap in the literature concerning the role of resampling in population annealing which this work attempts to close. The two-dimensional Ising model is used as a benchmarking system for this study. At first various resampling methods are implemented and numerically compared. In a second part the exact solution of the Ising model is utilized to create an artificial population annealing setting with effectively infinite Monte Carlo updates at each temperature. This limit is first performed on finite population sizes and subsequently extended to infinite populations. This allows us to look at resampling isolated from other parameters. Many results are expected to generalize to other systems.

Resampling schemes in population annealing — numerical results

Population annealing (PA) is a population-based algorithm that is designed for equilibrium simulations of thermodynamic systems with a rough free energy landscape. It is known to be more efficient in doing so than standard Markov chain Monte Carlo alone. The algorithm has a number of parameters that can be fine-tuned to improve performance. While there is some theoretical and numerical work regarding most of these parameters, there appears to be a gap in the literature concerning the role of resampling in PA. Here, we present a numerical comparison of a number of resampling schemes for PA simulations of the 2D Ising model.

Mapping Ligand Binding Landscapes using Weighted Ensembles of Trajectories

Protein-ligand interactions involve many conformational states. To bind the ligand must search across a rough free energy landscape that contains many metastable minima. A map of this “ligand binding landscape” would be useful during drug design, as it contains a set of possible binding modes that could be stabilized by modifications of the ligand. These maps would also include information about pose connectivity, both with each other and with the unbound state, which can provide insight into ligand binding kinetics, a useful parameter for drug design. Unfortunately, creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The Dickson lab has developed simulation algorithms (e.g. WExplore, REVO) that can broadly and efficiently explore free-energy landscapes and are capable of simulating ligand release pathways that occur on timescales as long as minutes. Importantly, these methods employ only unbiased trajectory segments, allowing for the construction of Markov state models (MSM) and conformation space networks that combine the results of multiple simulations. In this talk I will describe how these methods have revealed new insights into long timescale ligand release pathways, as well as the interconversion pathways between ligand poses. In particular I will focus on the application of WExplore to map the binding landscapes of two bromodomain-inhibitor systems (Brd4-MS436 and Baz2B-ICR7), and the use of Markov state models to predict stable poses. In this work we find that transitions between deeply bound states convert through the unbound state 81% of the time, implying a trial-and-error approach to ligand binding. This finding has consequences for both the way we think about ligand binding processes, and how molecular dynamics should be used in drug design.

Mapping the Ligand Binding Landscape

The interaction between a ligand and a protein involves a multitude of conformational states. To achieve a particular deeply bound pose, the ligand must search across a rough free-energy landscape with many metastable minima. Creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The WExplore enhanced sampling method is well suited to build these maps because it is designed to broadly explore free-energy landscapes and is capable of simulating ligand release pathways that occur on timescales as long as minutes. WExplore also uses only unbiased trajectory segments, allowing for the construction of Markov state models (MSMs) and conformation space networks that combine the results of multiple simulations. Here, we use WExplore to study two bromodomain-inhibitor systems using multiple docked starting poses (Brd4-MS436 and Baz2B-ICR7) and synthesize our results using a series of MSMs using time-lagged independent component analysis. Ranking the starting poses by exit rate agrees with the crystal structure pose in both cases. We also predict the most stable pose using the equilibrium populations from the MSM but find that the prediction is not robust as a function of MSM parameters. The simulated trajectories are synthesized into network models that visualize the entire binding landscape for each system, and we examine transition paths between deeply bound stable states. We find that, on average, transitions between deeply bound states convert through the unbound state 81% of the time, implying a trial-and-error approach to ligand binding. We conclude with a discussion of the implications of this result for both kinetics-based drug discovery and virtual screening pipelines that incorporate molecular dynamics.

Strong disorder leads to scale invariance in complex biological systems

Despite the innate complexity of the cell, emergent scale-invariant behavior is observed in many biological systems. We investigate one example of this phenomenon: the dynamics of large complexes in the bacterial cytoplasm. The observed dynamics of these complexes is scale invariant in three measures of dynamics: mean-squared displacement (MSD), velocity autocorrelation function, and the step-size distribution. To investigate the physical mechanism for this emergent scale invariance, we explore minimal models in which mobility is modeled as diffusion on a rough free-energy landscape in one dimension. We discover that all three scale-invariant characteristics emerge generically in the strong disorder limit. (Strong disorder is defined by the divergence of the ensemble-averaged hop time between lattice sites.) In particular, we demonstrate how the scale invariance of the relative step-size distribution can be understood from the perspective of extreme-value theory in statistics (EVT). We show that the Gumbel scale parameter is simply related to the MSD scaling parameter. The EVT mechanism of scale invariance is expected to be generic to strongly disordered systems and therefore a powerful tool for the analysis of other systems in biology and beyond.

Analysis and optimization of population annealing.

Population annealing is an easily parallelizable sequential Monte Carlo algorithm that is well suited for simulating the equilibrium properties of systems with rough free-energy landscapes. In this work we seek to understand and improve the performance of population annealing. We derive several useful relations between quantities that describe the performance of population annealing and use these relations to suggest methods to optimize the algorithm. These optimization methods were tested by performing large-scale simulations of the three-dimensional (3D) Edwards-Anderson (Ising) spin glass and measuring several observables. The optimization methods were found to substantially decrease the amount of computational work necessary as compared to previously used, unoptimized versions of population annealing. We also obtain more accurate values of several important observables for the 3D Edwards-Anderson model.

Population annealing simulations of a binary hard-sphere mixture.

Population annealing is a sequential Monte Carlo scheme well suited to simulating equilibrium states of systems with rough free energy landscapes. Here we use population annealing to study a binary mixture of hard spheres. Population annealing is a parallel version of simulated annealing with an extra resampling step that ensures that a population of replicas of the system represents the equilibrium ensemble at every packing fraction in an annealing schedule. The algorithm and its equilibration properties are described, and results are presented for a glass-forming fluid composed of a 50/50 mixture of hard spheres with diameter ratio of 1.4:1. For this system, we obtain precise results for the equation of state in the glassy regime up to packing fractions φ≈0.60 and study deviations from the Boublik-Mansoori-Carnahan-Starling-Leland equation of state. For higher packing fractions, the algorithm falls out of equilibrium and a free volume fit predicts jamming at packing fraction φ≈0.667. We conclude that population annealing is an effective tool for studying equilibrium glassy fluids and the jamming transition.

GPU-Accelerated Population Annealing Algorithm: Frustrated Ising Antiferromagnet on the Stacked Triangular Lattice

The population annealing algorithm is a novel approach to study systems with rough free-energy landscapes, such as spin glasses. It combines the power of simulated annealing, Boltzmann weighted differential reproduction and sequential Monte Carlo process to bring the population of replicas to the equilibrium even in the low-temperature region. Moreover, it provides a very good estimate of the free energy. The fact that population annealing algorithm is performed over a large number of replicas with many spin updates, makes it a good candidate for massive parallelism. We chose the GPU programming using a CUDA implementation to create a highly optimized simulation. It has been previously shown for the frustrated Ising antiferromagnet on the stacked triangular lattice with a ferromagnetic interlayer coupling, that standard Markov Chain Monte Carlo simulations fail to equilibrate at low temperatures due to the effect of kinetic freezing of the ferromagnetically ordered chains. We applied the population annealing to study the case with the isotropic intra- and interlayer antiferromagnetic coupling (J 2 /|J 1 | = −1). The reached ground states correspond to non-magnetic degenerate states, where chains are antiferromagnetically ordered, but there is no long-range ordering between them, which is analogical with Wannier phase of the 2D triangular Ising antiferromagnet.

Population annealing: Theory and application in spin glasses.

Population annealing is an efficient sequential Monte Carlo algorithm for simulating equilibrium states of systems with rough free-energy landscapes. The theory of population annealing is presented, and systematic and statistical errors are discussed. The behavior of the algorithm is studied in the context of large-scale simulations of the three-dimensional Ising spin glass and the performance of the algorithm is compared to parallel tempering. It is found that the two algorithms are similar in efficiency though with different strengths and weaknesses.

Robust Estimation of Diffusion-Optimized Ensembles for Enhanced Sampling

Molecular dynamics and Monte Carlo simulation can provide detailed structural insights into biological systems by exploring the distribution of accessible conformational states. However, due to the rough free energy landscapes characterizing many systems, obtaining converged statistics in such simulation remains an important challenge. In order to address this issue, over the past decades, a number of enhanced sampling strategies have been proposed, including methods like simulated annealing, replica exchange, metadynamics and generalized ensemble approaches such as multicanonical (flat histogram) sampling. Most of these methods were designed to soften the landscape, making it easier to overcome the free energy barriers. Recently, it was proposed that additional performance gains could be obtained by taking the position-dependent diffusion coefficient along the reaction coordinate into account, thus placing greater emphasis on regions diffusing slowly so that the round trip rate is maximized (Trebst et al. 2004. 10.1103/PhysRevE.70.046701). Although the diffusion-optimized ensemble approach has shown great promise, the practical application of it has been hindered by the serious challenge of estimating the diffusion coefficient accurately, especially at the preliminary stage of a simulation when only a small part of the phase space has been explored. In this study, a simple yet robust solution to this problem is presented. We propose replacing the standard global estimation procedure with a local estimation technique that was previously introduced in the context of trajectory analysis and reaction coordinate optimization (Krivov et al. 2008. 10.1073/pnas.0800228105). We demonstrate that the precision of the proposed estimation method is dramatically increased, requiring considerable fewer data compared to current state-of-the-art procedures. Finally, we show that with our method, diffusion-optimized sampling arises as a natural generalization of the standard multicanonical ensemble, estimating a simple histogram of transitions rather than the usual histogram of states.

Characterization of Viscous and Excluded Volume effects on the Folding Kinetics of the Tetraloop-Receptor Motif: Single Molecule Tests of Kramer's Theory

It is widely known that there exists a strong relationship between the solvent, including co-solutes (e.g. cations, proteins, nucleic acids), and RNA that cooperatively defines the kinetics of tertiary structure formation. Viscosity, a solvent property that is coupled to dynamics, is critical in RNA folding due to the propensity for populating kinetically trapped species on a rough free energy landscape. Also of importance is the role of macromolecular co-solutes that limit the volume available for molecules to sample (excluded volume effects). The effects of these properties are highlighted by observed in vivo diffusion coefficients that are 4-200x smaller than in aqueous solvents. In this work, Kramers' rate theory is used to describe the viscosity dependence of tetraloop-receptor docking kinetics at the single molecule level. Both rate constants, kfold and kunfold, decrease with increasing viscosity (increasing glycerol %), trends that are predicted by Kramers' theory in the over-damped limit. However, the same measurements made in high molecular weight PEG solutions showed folding rate constants are accelerated by 1-2 orders of magnitude. Scaled particle theory, describing a hard spheres PEG-RNA interaction, quantitatively predicts the stabilizing effect of excluded volume. Temperature dependent measurements show that the thermodynamics of docking are not perturbed (ΔH°= −23(2) kcal/mol and ΔS°= −76(6) cal/mol∗K) even in up to 50% glycerol (∼6 cP) and that solvent activation accounts for 4-6 kcal/mol of the folding enthalpy. In PEG solutions thermodynamics reveal that the folding is stabilized by a reduction in the entropy of the unfolded RNA ensemble (ΔΔS°>0). Thus, these studies validate physical models describing characteristic kinetic and thermodynamic trends for intramolecular structure formation attributed to viscosity and excluded volume.

Unraveling the Beautiful Complexity of Simple Lattice Model Polymers and Proteins Using Wang-Landau Sampling

We describe a class of “bare bones” models of homopolymers which undergo coil-globule collapse and proteins which fold into their native states in free space or into denatured states when captured by an attractive substrate as the temperature is lowered. We then show how, with the use of a properly chosen trial move set, Wang-Landau Monte Carlo sampling can be used to study the rough free energy landscape and ground (native) states of these intriguingly simple systems and thus elucidate their thermodynamic complexity.

Monte Carlo Methods for Rough Free Energy Landscapes: Population Annealing and Parallel Tempering

Parallel tempering and population annealing are both effective methods for simulating equilibrium systems with rough free energy landscapes. Parallel tempering, also known as replica exchange Monte Carlo, is a Markov chain Monte Carlo method while population annealing is a sequential Monte Carlo method. Both methods overcome the exponential slowing associated with high free energy barriers. The convergence properties and efficiency of the two methods are compared. For large systems, population annealing initially converges to equilibrium more rapidly than parallel tempering for the same amount of computational work. However, parallel tempering converges exponentially and population annealing inversely in the computational work so that ultimately parallel tempering approaches equilibrium more rapidly than population annealing.

Population annealing with weighted averages: A Monte Carlo method for rough free-energy landscapes

The population annealing algorithm introduced by Hukushima and Iba is described. Population annealing combines simulated annealing and Boltzmann weighted differential reproduction within a population of replicas to sample equilibrium states. Population annealing gives direct access to the free energy. It is shown that unbiased measurements of observables can be obtained by weighted averages over many runs with weight factors related to the free-energy estimate from the run. Population annealing is well suited to parallelization and may be a useful alternative to parallel tempering for systems with rough free-energy landscapes such as spin glasses. The method is demonstrated for spin glasses.

The glassy wormlike chain

We introduce a new model for the dynamics of a wormlike chain (WLC) in an environment that gives rise to a rough free energy landscape, which we name the glassy WLC. It is obtained from the common WLC by an exponential stretching of the relaxation spectrum of its long-wavelength eigenmodes, controlled by a single parameter . Predictions for pertinent observables such as the dynamic structure factor and the microrheological susceptibility exhibit the characteristics of soft glassy rheology and compare favourably with experimental data for reconstituted cytoskeletal networks and live cells. We speculate about the possible microscopic origin of the stretching, implications for the nonlinear rheology, and the potential physiological significance of our results.

Diffusive dynamics on multidimensional rough free energy surfaces

The dynamics of processes relevant to chemistry and biophysics on rough free energy landscapes is investigated using a recently developed algorithm to solve the Smoluchowski equation. Two different processes are considered: ligand rebinding in MbCO and protein folding. For the rebinding dynamics of carbon monoxide (CO) to native myoglobin (Mb) from locations around the active site, the two-dimensional free energy surface (FES) is constructed using extensive molecular dynamics simulations. The surface describes the minima in the A state (bound MbCO), CO in the distal pocket and in the Xe4 pocket, and the transitions between these states and allows to study the diffusion of CO in detail. For the folding dynamics of protein G, a previously determined two-dimensional FES was available. To follow the diffusive dynamics on these rough free energy surfaces, the Smoluchowski equation is solved using the recently developed hierarchical discrete approximation method. From the relaxation of the initial nonequilibrium distribution, experimentally accessible quantities such as the rebinding time for CO or the folding time for protein G can be calculated. It is found that the free energy barrier for CO in the Xe4 pocket and in the distal pocket (B state) closer to the heme iron is approximately 6 kcal/mol which is considerably larger than the inner barrier which separates the bound state and the B state. For the folding of protein G, a barrier of approximately 10 kcal/mol between the unfolded and the folded state is consistent with folding times of the order of milliseconds.

Protein Folding Kinetics and Thermodynamics from Atomistic Simulations

Determining protein folding kinetics and thermodynamics from all-atom molecular dynamics (MD) simulations without using experimental data represents a formidable scientific challenge because simulations can easily get trapped in local minima on rough free energy landscapes. This necessitates the computation of multiple simulation trajectories, which can be independent from each other or coupled in some manner, as, for example, in the replica exchange MD method. Here we present results obtained with a new analysis tool that allows the deduction of faithful kinetics data from a heterogeneous ensemble of simulation trajectories. The method is demonstrated on the decapeptide Chignolin for which we predict folding and unfolding time constants of 1.0 +/- 0.3 and 2.6 +/- 0.4 micros, respectively. We also derive the energetics of folding, and calculate a realistic melting curve for Chignolin.

Heterogeneous Folding of the trpzip Hairpin: Full Atom Simulation and Experiment

The β-hairpin trpzip2 can be tuned continuously from a two-state folder to folding on a rough energy landscape without a dominant refolding barrier. At high denaturant concentration, this extremely stable peptide exhibits a single apparent “two-state” transition temperature when monitored by different spectroscopic techniques. However, under optimal folding conditions the hairpin undergoes an unusual folding process with three clusters of melting transitions ranging from 15 °C to 160 °C, as monitored by 12 different experimental and computational observables. We explain this behavior in terms of a rough free energy landscape of the unfolded peptide caused by multiple tryptophan interactions and alternative backbone conformations. The landscape is mapped out by potentials of mean force derived from replica-exchange molecular dynamics simulations. Implications for deducing cooperativity from denaturant titrations, for the origin of folding cooperativity, and for the folding of thermophilic proteins are pointed out. trpzip is an excellent small tunable model system for the glass-like folding transitions predicted by landscape theory.

Interface pinning and slow ordering kinetics on infinitely ramified fractal structures

We investigate the time-dependent Ginzburg-Landau (TDGL) equation for a nonconserved order parameter on an infinitely ramified (deterministic) fractal lattice employing two alternative methods: the auxiliary field approach and a numerical method of integration of the equations of evolution. In the first case the domain size evolves with time as $L(t)\ensuremath{\sim}{t}^{{1/d}_{w}}$, where ${d}_{w}$ is the anomalous random-walk exponent associated with the fractal and differs from the normal value $2$, which characterizes all Euclidean lattices. Such a power-law growth is identical to the one observed in the study of the spherical model on the same lattice, but fails to describe the asymptotic behavior of the numerical solutions of the TDGL equation for a scalar order parameter. In fact, the simulations performed on a two dimensional Sierpinski carpet indicate that, after an initial stage dominated by a curvature reduction mechanism in the manner of Allen and Cahn [Acta. Metall. 27, 1085 (1979)], the system enters in a regime where the domain walls between competing phases are pinned by lattice defects. The lack of translational invariance determines a rough free-energy landscape, the existence of many metastable minima, and the suppression of the marginally stable modes, which in translationally invariant systems lead to power-law growth and self-similar patterns. On fractal structures, as the temperature vanishes the evolution is frozen since only thermally activated processes can sustain the growth of pinned domains.