Rugged Free Energy Landscape Research Articles

Networks of ion-pairs are responsible for the large differences in the thermal stability of two structurally similar aminopeptidases.

Aminopeptidases are an important class of enzymes for protein metabolism. Leucyl aminopeptidase (PepL) preferably removes leucine from the N-terminus of small peptides. PepL of Lacticaseibacillus casei was observed to be thermally unstable, while a structurally similar aminopeptidase T (AmpT) of Thermus thermophilus is highly stable. To understand the molecular interaction responsible for the large difference in their stability, molecular dynamics simulations were carried out to study the thermal stability of PepL and AmpT at 300 K to 450 K temperature range over 100 ns. PepL sampled a larger conformational space with a rugged free-energy landscape, while AmpT navigated a smoother energy landscape to reach the global minimum. The RMSD, RMSF, radius of gyration and principal component analysis suggested large movements in PepL than in AmpT with an increase in temperature. Analysis of residue-interaction network revealed AmpT possessing a greater number of low, medium and high energy contacts in comparison to PepL. AmpT showed a higher abundance of ion-pair clusters and ionic residues per cluster compared to PepL. Moreover, AmpT retained a greater number of high-energy contacts at elevated temperatures. These findings showed that the inherently lower stability of PepL originates from a comparatively smaller number of contacts and can be pivotal in engineering PepL for higher stability.

Monte Carlo simulations of glass-forming liquids beyond Metropolis.

Monte Carlo simulations are widely employed to measure the physical properties of glass-forming liquids in thermal equilibrium. Combined with local Monte Carlo moves, the Metropolis algorithm can also be used to simulate the relaxation dynamics, thus offering an efficient alternative to molecular dynamics. Monte Carlo simulations are, however, more versatile because carefully designed Monte Carlo algorithms can more efficiently sample the rugged free energy landscape characteristic of glassy systems. After a brief overview of Monte Carlo studies of glass-formers, we define and implement a series of Monte Carlo algorithms in a three-dimensional model of polydisperse hard spheres. We show that the standard local Metropolis algorithm is the slowest and that implementing collective moves or breaking detailed balance enhances the efficiency of the Monte Carlo sampling. We use time correlation functions to provide a microscopic interpretation of these observations. Seventy years after its invention, the Monte Carlo method remains the most efficient and versatile tool to compute low-temperature properties in supercooled liquids.

Transient excited states of the metamorphic protein Mad2 and their implications for function.

The spindle checkpoint complex is a key surveillance mechanism in cell division that prevents premature separation of sister chromatids. Mad2 is an integral component of this spindle checkpoint complex that recognizes cognate substrates such as Mad1 and Cdc20 in its closed (C-Mad2) conformation by fastening a "seatbelt" around short peptide regions that bind to the substrate recognition site. Mad2 is also a metamorphic protein that adopts not only the fold found in C-Mad2, but also a structurally distinct open conformation (O-Mad2) which is incapable of binding substrates. Here, we show using chemical exchange saturation transfer (CEST) and relaxation dispersion (CPMG) NMR experiments that Mad2 transiently populates three other higher free energy states with millisecond lifetimes, two in equilibrium with C-Mad2 (E1 and E2) and one with O-Mad2 (E3). E1 is a mimic of substrate-bound C-Mad2 in which the N-terminus of one C-Mad2 molecule inserts into the seatbelt region of a second molecule of C-Mad2, providing a potential pathway for autoinhibition of C-Mad2. E2 is the "unbuckled" conformation of C-Mad2 that facilitates the triage of molecules along competing fold-switching and substrate binding pathways. The E3 conformation that coexists with O-Mad2 shows fluctuations at a hydrophobic lock that is required for stabilizing the O-Mad2 fold and we hypothesize that E3 represents an early intermediate on-pathway towards conversion to C-Mad2. Collectively, the NMR data highlight the rugged free energy landscape of Mad2 with multiple low-lying intermediates that interlink substrate-binding and fold-switching, and also emphasize the role of molecular dynamics in its function.

Opposing roles of organic salts on mini-protein structure.

We investigated the effects of 1-ethyl-3-methylimidazolium chloride ([EMIM][Cl]) and choline chloride ([Chol][Cl]) on the local environment and conformational landscapes of Trp-cage and Trpzip4 mini-proteins using experimental and computational approaches. Fluorescence experiments and computational simulations revealed distinct behaviors of the mini-proteins in the presence of these organic salts. [EMIM][Cl] showed a strong interaction with Trp-cage, leading to fluorescence quenching and destabilization of its native structural interactions. Conversely, [Chol][Cl] had a negligible impact on Trp-cage fluorescence at low concentrations but increased it at high concentrations, indicating a stabilizing role. Computational simulations elucidated that [EMIM][Cl] disrupted the hydrophobic core packing and decreased proline-aromatic residue contacts in Trp-cage, resulting in a more exposed environment for Trp residues. In contrast, [Chol][Cl] subtly influenced the hydrophobic core packing, creating a hydrophobic environment near the tryptophan residues. Circular dichroism experiments revealed that [Chol][Cl] stabilized the secondary structure of both mini-proteins, although computational simulations did not show significant changes in secondary content at the explored concentrations. The simulations also demonstrated a more rugged free energy landscape for Trp-cage and Trpzip4 in [EMIM][Cl], suggesting destabilization of the tertiary structure for Trp-cage and secondary structure for Trpzip4. Similar fluorescence trends were observed for Trpzip4, with [EMIM][Cl] quenching fluorescence and exhibiting stronger interaction, while [Chol][Cl] increased the fluorescence at high concentrations. These findings highlight the interplay between [EMIM][Cl] and [Chol][Cl] with the mini-proteins and provide a detailed molecular-level understanding of how these organic salts impact the nearby surroundings and structural variations. Understanding such interactions is valuable for diverse applications, from biochemistry to materials science.

Liquid-Like States in Micelle-Forming Diblock Copolymer Melts.

Large cell self-consistent field theory (SCFT) solutions for a neat, micelle-forming diblock copolymer melt, initialized using the structure of a Lennard-Jones fluid, reveal the existence of a vast number of liquid-like states, with free energies of order 10-3 kBT per chain higher than the body-centered cubic (bcc) state near the order-disorder transition (ODT). Computation of the structure factor for these liquids at temperatures below the ODT indicates that their intermicellar distance is slightly swollen compared to bcc. In addition to providing a mean-field picture of the disordered micellar state, the number of liquid-like states and their near-degeneracy with the equilibrium bcc morphology suggest that self-assembly of micelle-forming diblock copolymers navigates a rugged free energy landscape with many local minima. This picture provides a basis for the anomalously slow ordering kinetics of particle-forming diblock copolymer melts observed in experiments.

Exploration, Representation, and Rationalization of the Conformational Phase Space of N-Glycans.

Despite their fundamental biological relevance, structure-property relationships in N-glycans are fundamentally lacking, and their highly multidimensional compositional and conformational phase spaces remain largely unexplored. The torsional flexibility of the glycosidic linkages and the ring dynamics result in wide, rugged free-energy landscapes that are difficult to sample in molecular dynamics simulations. We show that a novel enhanced-sampling scheme combining replica exchange with solute and collective-variable tempering, enabling transitions over all relevant energy barriers, delivers converged distributions of solvated N-glycan conformers. Several dimensionality-reduction algorithms are compared and employed to generate conformational free-energy maps in two dimensions. Together with an originally developed conformation-based nomenclature scheme that uniquely identifies glycan conformers, our modeling procedure is applied to reveal the effect of chemical substitutions on the conformational ensemble of selected high-mannose-type and complex glycans. Moreover, the structure-prediction capabilities of two commonly used glycan force fields are assessed via the theoretical prediction of experimentally available nuclear magnetic resonance J-coupling constants. The results especially confirm the key role of ω and ψ torsion angles in discriminating between different conformational states and suggest an intriguing correlation between the torsional and ring-puckering degrees of freedom that may be biologically relevant.

Overcoming Free-Energy Barriers with a Seamless Combination of a Biasing Force and a Collective Variable-Independent Boost Potential.

Amid collective-variable (CV)-based importance-sampling algorithms, a hybrid of the extended adaptive biasing force and the well-tempered metadynamics algorithms (WTM-eABF) has proven particularly cost-effective for exploring the rugged free-energy landscapes that underlie biological processes. However, as an inherently CV-based algorithm, this hybrid scheme does not explicitly accelerate sampling in the space orthogonal to the chosen CVs, thereby limiting its efficiency and accuracy, most notably in those cases where the slow degrees of freedom of the process at hand are not accounted for in the model transition coordinate. Here, inspired by Gaussian-accelerated molecular dynamics (GaMD), we introduce the same CV-independent harmonic boost potential into WTM-eABF, yielding a hybrid algorithm coined GaWTM-eABF. This algorithm leans on WTM-eABF to explore the transition coordinate with a GaMD-mollified potential and recovers the unbiased free-energy landscape through thermodynamic integration followed by proper reweighting. As illustrated in our numerical tests, GaWTM-eABF effectively overcomes the free-energy barriers in orthogonal space and correctly recovers the unbiased potential of mean force (PMF). Furthermore, applying both GaWTM-eABF and WTM-eABF to two biologically relevant processes, namely, the reversible folding of (i) deca-alanine and (ii) chignolin, our results indicate that GaWTM-eABF reduces the uncertainty in the PMF calculation and converges appreciably faster than WTM-eABF. Obviating the need of multiple-copy strategies, GaWTM-eABF is a robust, computationally efficient algorithm to surmount the free-energy barriers in orthogonal space and maps with utmost fidelity the free-energy landscape along selections of CVs. Moreover, our strategy that combines WTM-eABF with GaMD can be easily extended to other biasing-force algorithms.

Volume-scaled common nearest neighbor clustering algorithm with free-energy hierarchy.

The combination of Markov state modeling (MSM) and molecular dynamics (MD) simulations has been shown in recent years to be a valuable approach to unravel the slow processes of molecular systems with increasing complexity. While the algorithms for intermediate steps in the MSM workflow such as featurization and dimensionality reduction have been specifically adapted to MD datasets, conventional clustering methods are generally applied to the discretization step. This work adds to recent efforts to develop specialized density-based clustering algorithms for the Boltzmann-weighted data from MD simulations. We introduce the volume-scaled common nearest neighbor (vs-CNN) clustering that is an adapted version of the common nearest neighbor (CNN) algorithm. A major advantage of the proposed algorithm is that the introduced density-based criterion directly links to a free-energy notion via Boltzmann inversion. Such a free-energy perspective allows a straightforward hierarchical scheme to identify conformational clusters at different levels of a generally rugged free-energy landscape of complex molecular systems.

Simulated tempering with irreversible Gibbs sampling techniques.

We present here two novel algorithms for simulated tempering simulations, which break the detailed balance condition (DBC) but satisfy the skewed detailed balance to ensure invariance of the target distribution. The irreversible methods we present here are based on Gibbs sampling and concern breaking DBC at the update scheme of the temperature swaps. We utilize three systems as a test bed for our methods: a Markov chain Monte Carlo simulation on a simple system described by a one-dimensional double well potential, the Ising model, and molecular dynamics simulations on alanine pentapeptide (ALA5). The relaxation times of inverse temperature, magnetic susceptibility, and energy density for the Ising model indicate clear gains in sampling efficiency over conventional Gibbs sampling techniques with DBC and also over the conventionally used simulated tempering with the Metropolis-Hastings (MH) scheme. Simulations on ALA5 with a large number of temperatures indicate distinct gains in mixing times for inverse temperature and consequently the energy of the system compared to conventional MH. With no additional computational overhead, our methods were found to be more efficient alternatives to the conventionally used simulated tempering methods with DBC. Our algorithms should be particularly advantageous in simulations of large systems with many temperature ladders, as our algorithms showed a more favorable constant scaling in Ising spin systems as compared with both reversible and irreversible MH algorithms. In future applications, our irreversible methods can also be easily tailored to utilize a given dynamical variable other than temperature to flatten rugged free energy landscapes.

Real-time observation of ligand-induced allosteric transitions in a PDZ domain

While allostery is of paramount importance for protein regulation, the underlying dynamical process of ligand (un)binding at one site, resulting time evolution of the protein structure, and change of the binding affinity at a remote site are not well understood. Here the ligand-induced conformational transition in a widely studied model system of allostery, the PDZ2 domain, is investigated by transient infrared spectroscopy accompanied by molecular dynamics simulations. To this end, an azobenzene-derived photoswitch is linked to a peptide ligand in a way that its binding affinity to the PDZ2 domain changes upon switching, thus initiating an allosteric transition in the PDZ2 domain protein. The subsequent response of the protein, covering four decades of time, ranging from ∼1 ns to ∼μs, can be rationalized by a remodeling of its rugged free-energy landscape, with very subtle shifts in the populations of a small number of structurally well-defined states. It is proposed that structurally and dynamically driven allostery, often discussed as limiting scenarios of allosteric communication, actually go hand-in-hand, allowing the protein to adapt its free-energy landscape to incoming signals.

Combined Force-Frequency Sampling for Simulation of Systems Having Rugged Free Energy Landscapes.

An adaptive, machine learning-based sampling method is presented for simulation of systems having rugged, multidimensional free energy landscapes. The method's main strength resides in its ability to learn both from the frequency of visits to distinct states and the generalized force estimates that arise in a system as it evolves in phase space. This is accomplished by introducing a self-integrating artificial neural network, which generates an estimate of the free energy directly from its derivatives. The usefulness of the proposed combined approach is examined in the context of two concrete examples, namely, an alanine dipeptide molecule in water and a polymer diffusing through a narrow pore. This new method is found to be robust, faster, and more accurate than approaches that rely only on frequency-based or generalized force-based estimations. After combining the proposed approach with overfill protection and support for sparse data storage and training, the method is shown to be more effective than comparable, previously available techniques and capable of scaling efficiently to larger numbers of collective variables.

Advances in enhanced sampling molecular dynamics simulations for biomolecules

Molecular dynamics simulation has emerged as a powerful computational tool for studying biomolecules as it can provide atomic insights into the conformational transitions involved in biological functions. However, when applied to complex biological macromolecules, the conformational sampling ability of conventional molecular dynamics is limited by the rugged free energy landscapes, leading to inherent timescale gaps between molecular dynamics simulations and real biological processes. To address this issue, several advanced enhanced sampling methods have been proposed to improve the sampling efficiency in molecular dynamics. In this review, the theoretical basis, practical applications, and recent improvements of both constraint and unconstrained enhanced sampling methods are summarized. Further-more, the combined utilizations of different enhanced sampling methods that take advantage of both approaches are also briefly discussed.

Replica-Exchange Methods for Biomolecular Simulations.

In this study, a replica-exchange method was developed to overcome conformational sampling difficulties in computer simulations of spin glass or other systems with rugged free-energy landscapes. This method was then applied to the protein-folding problem in combination with molecular dynamics (MD) simulation. Owing to its simplicity and sampling efficiency, the replica-exchange method has been applied to many other biological problems and has been continuously improved. The method has often been combined with other sampling techniques, such as umbrella sampling, free-energy perturbation, metadynamics, and Gaussian accelerated MD (GaMD). In this chapter, we first summarize the original replica-exchange molecular dynamics (REMD) method and discuss how new algorithms related to the original method are implemented to add new features. Heterogeneous and flexible structures of an N-glycan in a solution are simulated as an example of applications by REMD, replica exchange with solute tempering, and GaMD. The sampling efficiency of these methods on the N-glycan system and the convergence of the free-energy changes are compared. REMD simulation protocols and trajectory analysis using the GENESIS software are provided to facilitate the practical use of advanced simulation methods.

Numerical algorithms for solving self-consistent field theory reversely for block copolymer systems.

Besides dictating the equilibrium phase diagram, the rugged free-energy landscape of AB block copolymers gives rise to a multitude of non-equilibrium phenomena. Self-consistent field theory (SCFT) can be employed to calculate the mean-field free energy, , of a non-equilibrium unstable state that is characterized by a given spatial density distribution, , in the incompressible system. Such a free-energy functional is the basis of describing the structure formation by dynamic SCFT techniques or the identification of minimum free-energy paths via the string method. The crucial step consists in computing the external potential fields that generate the given density distribution in the corresponding system of non-interacting copolymers, i.e., the potential-to-density relation employed in equilibrium SCFT calculations has to be inverted (reverse SCFT calculation). We describe, generalize, and evaluate the computational efficiency of two different numerical algorithms for this reverse SCFT calculation-the Debye-function algorithm based on the structure factor and the field-theoretic umbrella-potential (FUP) algorithm. In contrast to the Debye-function algorithm, the FUP algorithm only yields the exact mean-field values of the given target densities in the limit of a strong umbrella potential, and we devise a two-step variant of the FUP algorithm that significantly mitigates this issue. For Gaussian copolymers, the Debye-function algorithm is more efficient for highly unstable states that are far away from the equilibrium, whereas the improved FUP algorithm outperforms the Debye-function algorithm closer to metastable states and is easily transferred to more complex molecular architectures.

Layered nested Markov chain Monte Carlo.

A configurational sampling algorithm based on nested layerings of Markov chains (Layered Nested Markov Chain Monte Carlo or L-NMCMC) is presented for simulations of systems characterized by rugged free energy landscapes. The layerings are generated using a set of auxiliary potential energy surfaces. The implementation of the method is demonstrated in the context of a rugged, two-dimensional potential energy surface. The versatility of the algorithm is next demonstrated on a simple, many-body system, namely, a canonical Lennard-Jones fluid in the liquid state. In that example, different layering schemes and auxiliary potentials are used, including variable cutoff distances and excluded-volume tempering. In addition to calculating a variety of properties of the system, it is also shown that L-NMCMC, when combined with a free-energy perturbation formalism, provides a straightforward means to construct approximate free-energy surfaces at no additional computational cost using the sampling distributions of each auxiliary Markov chain. The proposed L-NMCMC scheme is general in that it could be complementary to any number of methods that rely on sampling from a target distribution or methods that exploit a hierarchy of time scales and/or length scales through decomposition of the potential energy.

Kinetics of pattern formation in symmetric diblock copolymer melts.

In equilibrium, copolymers self-assemble into spatially modulated phases with long-range order. When the system is quenched far below the order-disorder transition temperature, however, such an idealized, defect-free structure is difficult to obtain in experiments and simulations, instead a fingerprint-like structure forms. The relaxation toward long-range order is very protracted because it involves numerous thermally activated processes, and the rugged free-energy landscape has been likened to that of glass-forming systems. Using large-scale particle-based simulations of high-aspect-ratio, quasi-two-dimensional systems with periodic boundary condition, we study the kinetics of structure formation in symmetric, lamella-forming diblock copolymers after a quench from the disordered state. We characterize the ordering process by the correlation length of the lamellar structure and its Euler characteristic and observe that the growth of the correlation length and the rate of change of the Euler characteristic significantly slow down in the range of incompatibilities, 15 ≤ χN ≤ 20, studied. The increase of the time scale of ordering is, however, gradual. The density fields of snapshots of the particle-based simulations are used as starting values for self-consistent field theory (SCFT) calculations. The latter converge to the local, metastable minimum of the free-energy basin. This combination of particle-based simulations and SCFT calculations allows us to relate an instantaneous configuration of the particle-based model to a corresponding metastable free-energy minimum of SCFT-the inherent morphology-and we typically observe that a change of a free-energy basin is associated with a change of the Euler characteristic of the particle-based morphology, i.e., changes of free-energy basins are correlated to changes of the domain topology. Subsequently, we employ the string method in conjunction with SCFT to study the minimum free-energy paths (MFEPs) of changes of the domain topology. Since the time scales of relaxing toward the inherent morphology within a free-energy basin and jumps between free-energy basins are not well separated, the MFEP may overestimate the barriers encountered in the course of ordering.

REFOLDING OF HOMOPOLYMER UNDER QUENCHED FORCE

Recently single molecule force spectroscopy has become an useful tool to study protein, DNA and RNA. However, very little attention was paid to homopolymer which plays an important role in many domains of science. In this paper we make the first attempt to decipher the free energy landscape of homopolymer using the external force as reaction coordinate. The impact of the quenched force on the free energy landscape was studied using simplified coarse-grain Go model. Similar to protein, we have obtained a clear switch from the thermal regime to force-driven regime. The distance between the denatured state and transition state in the temperature-driven regime is smaller than in the force-driven one. Having a rugged free energy landscape without a pronounced funnel the homopolymer folding is much slower than that of protein making study of homopolymer very time consuming.

Slow domain reconfiguration causes power-law kinetics in a two-state enzyme

Protein dynamics are typically captured well by rate equations that predict exponential decays for two-state reactions. Here, we describe a remarkable exception. The electron-transfer enzyme quiescin sulfhydryl oxidase (QSOX), a natural fusion of two functionally distinct domains, switches between open- and closed-domain arrangements with apparent power-law kinetics. Using single-molecule FRET experiments on time scales from nanoseconds to milliseconds, we show that the unusual open-close kinetics results from slow sampling of an ensemble of disordered domain orientations. While substrate accelerates the kinetics, thus suggesting a substrate-induced switch to an alternative free energy landscape of the enzyme, the power-law behavior is also preserved upon electron load. Our results show that the slow sampling of open conformers is caused by a variety of interdomain interactions that imply a rugged free energy landscape, thus providing a generic mechanism for dynamic disorder in multidomain enzymes.

Inference of Calmodulin's Ca2+-Dependent Free Energy Landscapes via Gaussian Mixture Model Validation.

A free energy landscape estimation method based on the well-known Gaussian mixture model (GMM) is used to compare the efficiencies of thermally enhanced sampling methods with respect to regular molecular dynamics. The simulations are carried out on two binding states of calmodulin, and the free energy estimation method is compared with other estimators using a toy model. We show that GMM with cross-validation provides a robust estimate that is not subject to overfitting. The continuous nature of Gaussians provides better estimates on sparse data than canonical histogramming. We find that diffusion properties determine the sampling method effectiveness, such that diffusion-dominated apo calmodulin is most efficiently sampled by regular molecular dynamics, while holo calmodulin, with its rugged free energy landscape, is better sampled by enhanced sampling methods.

Massively parallel multicanonical simulations

Generalized-ensemble Monte Carlo simulations such as the multicanonical method and similar techniques are among the most efficient approaches for simulations of systems undergoing discontinuous phase transitions or with rugged free-energy landscapes. As Markov chain methods, they are inherently serial computationally. It was demonstrated recently, however, that a combination of independent simulations that communicate weight updates at variable intervals allows for the efficient utilization of parallel computational resources for multicanonical simulations. Implementing this approach for the many-thread architecture provided by current generations of graphics processing units (GPUs), we show how it can be efficiently employed with of the order of 104 parallel walkers and beyond, thus constituting a versatile tool for Monte Carlo simulations in the era of massively parallel computing. We provide the fully documented source code for the approach applied to the paradigmatic example of the two-dimensional Ising model as starting point and reference for practitioners in the field. Program summaryProgram Title: cudamucaProgram Files doi:http://dx.doi.org/10.17632/tzhfpdymv9.1Licensing provisions: Creative Commons Attribution license (CC BY 4.0)Programming language: C++, CUDAExternal routines/libraries: NVIDIA CUDA Toolkit 6.5 or newerNature of problem: The program determines weights for a multicanonical simulation of the 2D Ising model to result in a flat energy histogram. A final production run with these weights provides an estimate of the density of states of the model.Solution method: The code uses a parallel variant of the multicanonical method employing many parallel walkers that accumulate a common histogram. The resulting histogram is used to determine the weight function for the next iteration. Once the iteration has converged, simulations visit all possible energies with the same probability.Additional comments including restrictions and unusual features: The system size and size of the population of replicas are limited depending on the memory of the GPU device used. Code repository at https://github.com/CQT-Leipzig/cudamuca.