Multidimensional Free-energy Landscapes Research Articles

Machine Learning Derived Collective Variables for the Study of Protein Homodimerization in Membrane.

The accurate calculation of equilibrium constants for protein-protein association is of fundamental importance to quantitative biology and remains an outstanding challenge for computational biophysics. Traditionally, equilibrium constants have been computed from one-dimensional free energy surfaces derived from sampling along a single collective variable. Importantly, recent advances in enhanced sampling methodology have facilitated the characterization of multidimensional free energy landscapes, often exposing multiple thermodynamically important minima missed by more restrictive sampling methods. A key to the effectiveness of this multidimensional sampling approach is the identification of collective variables that effectively define the configurational space of dissociated and associated states. Here we present the application of two machine learning methods for the unbiased determination of collective variables for enhanced sampling and analysis of protein-protein association. Our results both validate prior work, based on intuition derived collective variables, and demonstrate the effectiveness of the machine learning methods for the identification of collective variables for association reactions in complex biomolecular systems.

Effective data-driven collective variables for free energy calculations from metadynamics of paths.

A variety of enhanced sampling (ES) methods predict multidimensional free energy landscapes associated with biological and other molecular processes as a function of a few selected collective variables (CVs). The accuracy of these methods is crucially dependent on the ability of the chosen CVs to capture the relevant slow degrees of freedom of the system. For complex processes, finding such CVs is the real challenge. Machine learning (ML) CVs offer, in principle, a solution to handle this problem. However, these methods rely on the availability of high-quality datasets-ideally incorporating information about physical pathways and transition states-which are difficult to access, therefore greatly limiting their domain of application. Here, we demonstrate how these datasets can be generated by means of ES simulations in trajectory space via the metadynamics of paths algorithm. The approach is expected to provide a general and efficient way to generate efficient ML-based CVs for the fast prediction of free energy landscapes in ES simulations. We demonstrate our approach with two numerical examples, a 2D model potential and the isomerization of alanine dipeptide, using deep targeted discriminant analysis as our ML-based CV of choice.

Conformational Changes in Surface-Immobilized Proteins Measured Using Combined Atomic Force and Fluorescence Microscopy.

Biological organisms rely on proteins to perform the majority of their functions. Most protein functions are based on their physical motions (conformational changes), which can be described as transitions between different conformational states in a multidimensional free-energy landscape. A comprehensive understanding of this free-energy landscape is therefore of paramount importance for understanding the biological functions of proteins. Protein dynamics includes both equilibrium and nonequilibrium motions, which typically exhibit a wide range of characteristic length and time scales. The relative probabilities of various conformational states in the energy landscape, the energy barriers between them, their dependence on external parameters such as force and temperature, and their connection to the protein function remain largely unknown for most proteins. In this paper, we present a multimolecule approach in which the proteins are immobilized at well-defined locations on Au substrates using an atomic force microscope (AFM)-based patterning method called nanografting. This method enables precise control over the protein location and orientation on the substrate, as well as the creation of biologically active protein ensembles that self-assemble into well-defined nanoscale regions (protein patches) on the gold substrate. We performed AFM-force compression and fluorescence experiments on these protein patches and measured the fundamental dynamical parameters such as protein stiffness, elastic modulus, and transition energies between distinct conformational states. Our results provide new insights into the processes that govern protein dynamics and its connection to protein function.

How and When Does an Enzyme React? Unraveling α-Amylase Catalytic Activity with Enhanced Sampling Techniques

Enzymatic catalysis is a complex process that can involve multiple conformations of the enzyme:substrate complex and several competitive reaction pathways, resulting in a multi-dimensional free energy landscape. The study of enzymatic activity often requires deep knowledge of the system to establish the catalytic mechanism and identify the possible reactive conformations of the complex. Here, we present an enhanced sampling and machine learning-based approach to explore the catalytic reaction space and characterize the transformation from reactive to non-reactive conformations with minimal a priori knowledge of the system. We applied this approach to study the rate-determining step of the glycolysis reaction of maltopentose catalyzed by human pancreatic α-amylase, an important enzyme in glucose production as well as a major drug target for the treatment of type-II diabetes. We unravel the complexity of the enzymatic reaction, reveal three binding modes of the substrate within the active site, and highlight the role of water in the catalytic process and in the stepwise conversion of reaction-ready to non-reactive conformations. Overall, these insights offer atomistic details on the catalytic mechanism and dynamics of the active site, allowing one to shed light on two fundamental questions in enzymatic catalysis, that is how and when does an enzyme react?

Molecular insights into the stereospecificity of arginine in RNA tetraloop folding.

One of the hypotheses for the homochirality of amino acids in the context of the origin of life is that only a particular stereoisomer provides preferential stability to RNA folding by acting as a chemical chaperon. However, the effect at the molecular level is not well understood. This study provides a molecular understanding of such preferential stability for a small GAAA RNA tetraloop in the presence of chiral arginine through a multidimensional free energy landscape constructed using a combination of umbrella sampling and parallel bias metadynamics (PBMetaD) simulations. We show that the origin of the chirality difference in RNA folding-unfolding dynamics is due to differences in the configurational diversity of RNA in adopting various non-natural conformations that accompany the diverse binding modes of D-arginine and L-arginine. We show that while D-arginine stabilizes the native folded state of RNA, L-arginine destabilizes it. Furthermore, free energy calculations on the binding of D- and L-arginine reveal a specific geometric constraint that helps D-arginine to stack with the terminal base pairs of RNA and pushes L-arginine for groove binding.

Markovian Weighted Ensemble Milestoning (M-WEM): Long-Time Kinetics from Short Trajectories.

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-time-scale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond time-scale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy but it also allows for a scheme to reconstruct a multidimensional free-energy landscape along additional degrees of freedom, which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to 4 orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential applications in kinetics and free-energy-based computational drug design.

Force-Correction Analysis Method for Derivation of Multidimensional Free-Energy Landscapes from Adaptively Biased Replica Simulations.

A methodology is proposed for the calculation of multidimensional free-energy landscapes of molecular systems, based on analysis of multiple molecular dynamics trajectories wherein adaptive biases have been applied to enhance the sampling of different collective variables. In this approach, which we refer to as the Force-Correction Analysis Method (FCAM), local averages of the total and biasing forces are evaluated post hoc, and the latter are subtracted from the former to obtain unbiased estimates of the mean force across collective-variable space. Multidimensional free-energy surfaces and minimum free-energy pathways are then derived by integrating the mean-force landscape with a kinetic Monte Carlo algorithm. To evaluate the proposed method, a series of numerical tests and comparisons with existing approaches were carried out for small molecules, peptides, and proteins, based on all-atom trajectories generated with standard, concurrent, and replica-exchange metadynamics in collective-variable spaces ranging from one to six dimensional. The tests confirm the correctness of the FCAM formulation and demonstrate that calculated mean forces and free energies converge rapidly and accurately, outperforming other methods used to unbias this kind of simulation data.

Unfolding of Dynamical Events in the Early Stage of Insulin Dimer Dissociation.

The dissociation of an insulin dimer is an important biochemical event that could also serve as a prototype of dissociations in similar biomolecular assemblies. We use a recently developed multidimensional free energy landscape for insulin dimer dissociation to unearth the microscopic and mechanistic aspects of the initial stages of the process that could hold the key to understanding the stability and the rate. The following sequence of events occurs in the initial stages: (i) The backbone hydrogen bonds break partially at the antiparallel β-sheet junction, (ii) the two α-helices (chain B) move away from each other while several residues (chain A) move closer, and (iii) a flow of adjacent water molecules occurs into the junction region. Interestingly, the intermonomeric center-to-center distance does not increase, but the number of native contacts exhibits a sharp decrease. Subsequent steps involve further disengagement of hydrophobic groups. This process is slow because of an entropic bottleneck created by the existence of the large configuration space available in the native state (NS), which is inhabited by low-frequency conformational fluctuations. We carry out a density-of-states analyses in the dimer NS to unearth distinctive features not present in the monomers. These low-frequency modes are also responsible for a large entropic stabilization of the NS. Hydrophobic disengagement in the early stage leads to the formation of a twisted intermediate state which itself is a metastable minimum (IS-1). The subsequent progress leads to another dimeric complex (IS-2), which is on the dissociative pathway and characterized by a further decrease in the native contacts. The dissociation process provides insights into the workings of a biomolecular assembly.

Deciphering the Self-Cleavage Reaction Mechanism of Hairpin Ribozyme.

Hairpin ribozyme catalyzes the reversible self-cleavage of phosphodiester bonds which plays prominent roles in key biological processes involving RNAs. Despite impressive advances on ribozymatic self-cleavage, critical aspects of its molecular reaction mechanism remain controversially debated. Here, we generate and analyze the multidimensional free energy landscape that underlies the reaction using extensive QM/MM metadynamics simulations to investigate in detail the full self-cleavage mechanism. This allows us to answer several pertinent yet controversial questions concerning activation of the 2'-OH group, the mechanistic role of water molecules present in the active site, and the full reaction pathway including the structures of transition states and intermediates. Importantly, we find that a sufficiently unrestricted reaction subspace must be mapped using accelerated sampling methods in order to compute the underlying free energy landscape. It is shown that lower-dimensional sampling where the bond formation and cleavage steps are coupled does not allow the system to sufficiently explore the landscape. On the basis of a three-dimensional free energy surface spanned by flexible generalized coordinates, we find that 2'-OH is indirectly activated by adjacent G8 nucleobase in conjunction with stabilizing H-bonding involving water. This allows the proton of the 2'-OH group to directly migrate toward the 5'-leaving group via a nonbridging oxygen of the phosphodiester link. At variance with similar enzymatic processes where water wires connected to protonable side chains of the protein matrix act as transient proton shuttles, no such de/reprotonation events of water molecules are found to be involved in this ribozymatic transesterification. Overall, our results support an acid-catalyzed reaction mechanism where A38 nucleobase directly acts as an acid whereas G8, in stark contrast, participates only indirectly via stabilizing the nascent nucleophile for subsequent attack. Moreover, we conclude that self-cleavage of hairpin ribozyme follows an AN + DN two-step associative pathway where the rate-determining step is the cleavage of the phosphodiester bond. These results provide a major advancement in our understanding of the unique catalytic mechanism of hairpin ribozyme which will fruitfully impact on the design of synthetic ribozymes.

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.

Reconstructing free-energy landscapes for cyclic molecular motors using full multidimensional or partial one-dimensional dynamic information.

Diffusion on a free-energy landscape is a fundamental framework for describing molecular motors. In the landscape framework, energy conversion between different forms of energy, e.g., chemical and mechanical, is explicitly described using multidimensional nonseparable potential landscapes. We present a k-space method for reconstructing multidimensional free-energy landscapes from stochastic single-molecule trajectories. For a variety of two-dimensional model potential landscapes, we demonstrate the robustness of the method by reconstructing the landscapes using full dynamic information, i.e., simulated two-dimensional stochastic trajectories. We then consider the case where the stochastic trajectory is known only along one dimension. With this partial dynamic information, the reconstruction of the full two-dimensional landscape is severely limited in the majority of cases. However, we reconstruct effective one-dimensional landscapes for the two-dimensional model potentials. We discuss the interpretation of the one-dimensional landscapes and identify signatures of energy conversion. Finally, we consider the implications of these results for biological molecular motors experiments.

Residue-Specialized Membrane Poration Kinetics of Melittin and Its Variants: Insight from Mechanistic Landscapes**Supported by the National Natural Science Foundation of China under Grant Nos. 21422404, 21774092, U1532108, and 21728502; the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions; and the Natural Science Foundation of Jiangsu Province of China under Grant Nos. BK20171207 and BK20171210

Pore-forming peptides have promising potentials for biomedical uses due to their ability to permeabilize cell membranes. However, to molecularly engineer them for practical applications is still blocked by the poor understanding of the specific roles of individual residues in peptides’ activity. Herein, using an advanced computational approach that combines Coarse-Grained molecular dynamics and well-tempered metadynamics, the membrane activities of melittin, a representative pore-forming peptide, and its gain-of-function variants, are characterized from the kinetics and thermodynamics perspectives. Unbiased simulations elucidate the molecular details of peptide-induced membrane poration; during which, some vital intermediate states, including the aggregation and U-shape configuration formation of peptides in the membrane, are observed and further applied as collective variables to construct the multi-dimensional free energy landscapes of the peptide-membrane interactions. Such a combination of kinetic and thermodynamic descriptions of the interaction process provides crucial information of residue-specialized contribution in chain conformation and consequently membrane perforation ability of the peptide. It is found that residues at the kink part (e.g. Thr) determine the chain flexibility and U-shape bending of the peptide, while residues near the C-terminus (e.g. Arg and Lys) are responsible for recruiting neighboring peptides for inter-molecular cooperation; the probable reaction pathway and the poration efficiency are consequently regulated. These results are helpful for a comprehensive understanding of the complicated molecular mechanism of pore-forming peptides and pave the way to rationally design and/or engineer the peptides for practical applications.

Cleaning Up Mechanistic Debris Generated by Twister Ribozymes Using Computational RNA Enzymology.

The catalytic properties of RNA have been a subject of fascination and intense research since their discovery over 30 years ago. Very recently, several classes of nucleolytic ribozymes have emerged and been characterized structurally. Among these, the twister ribozyme has been center-stage, and a topic of debate about its architecture and mechanism owing to conflicting interpretations of different crystal structures, and in some cases conflicting interpretations of the same functional data. In the present work, we attempt to clean up the mechanistic "debris" generated by twister ribozymes using a comprehensive computational RNA enzymology approach aimed to provide a unified interpretation of existing structural and functional data. Simulations in the crystalline environment and in solution provide insight into the origins of observed differences in crystal structures, and coalesce on a common active site architecture, and dynamical ensemble in solution. We use GPU-accelerated free energy methods with enhanced sampling to ascertain microscopic nucleobase pK a values of the implicated general acid and base, from which predicted activity-pH profiles can be compared directly with experiments. Next, ab initio quantum mechanical/molecular mechanical (QM/MM) simulations with full dynamic solvation under periodic boundary conditions are used to determine mechanistic pathways through multi-dimensional free energy landscapes for the reaction. We then characterize the rate-controlling transition state, and make predictions about kinetic isotope effects and linear free energy relations. Computational mutagenesis is performed to explain the origin of rate effects caused by chemical modifications and make experimentally testable predictions. Finally, we provide evidence that helps to resolve conflicting issues related to the role of metal ions in catalysis. Throughout each stage, we highlight how a conserved L-platform structural motif, to- gether with a key L-anchor residue, forms the characteristic active site scaffold enabling each of the catalytic strategies to come together not only for the twister ribozyme, but the majority of the known small nucleolytic ribozyme classes.

Characterization of heterogeneous intermediate ensembles on the guanidinium chloride-induced unfolding pathway of β-lactoglobulin

Folding pathway of β-LgA (β-lactoglobulin) evolves through the conformational α→β transition. The α→β transition is a molecular hallmark of various neurodegenerative diseases. Thus, β-LgA may serve as a good model for understanding molecular mechanism of protein aggregation involved in neurodegenerative diseases. Here, we studied the conformational dynamics of β-LgA in 6 M GdmCl at different temperatures using MD simulations. Structural order parameters such as RMSD, Rg, SASA, native contacts (Q), hydrophobic distal-matrix and free-energy landscape (FEL) were used to investigate the conformational transitions. Our results show that GdmCl destabilizes secondary and tertiary structure of β-LgA by weakening the hydrophobic interactions and hydrogen bond network. Multidimensional FEL shows the presence of different unfolding intermediates at 400 K. I1 is long-lived intermediate which has mostly intact native secondary structure, but loose tertiary structure. I2 is structurally compact intermediate formed after the partial loss of secondary structure. The transiently and infrequently buried evolution of W19 shows that intermediate conformational ensembles are structurally heterogeneous. We observed that the intermediate conformations are largely stabilized by non-native H-bonds. The outcome of this work provides the molecular details of intermediates trapped due to non-native interactions that may be regarded as pathogenic conformations involved in neurodegenerative diseases.Communicated by Ramaswamy H. Sarma

Unfolding Dynamics of Ubiquitin from Constant Force MD Simulation: Entropy-Enthalpy Interplay Shapes the Free-Energy Landscape.

Force probe methods are routinely used to study conformational transitions of biomolecules at single-molecule level. In contrast to simple kinetics, some proteins show complex response to mechanical perturbations that is manifested in terms of unusual force-dependent kinetics. Here, we study, via fully atomistic molecular dynamics simulations, constant force-induced unfolding of ubiquitin protein. Our simulations reveal a crossover at an intermediate force (about 400 pN) in the unfolding rate versus force curve. We find by calculation of multidimensional free-energy landscape (FEL) of the protein that the complex unfolding kinetics is intimately related to the force-dependent modifications in the FEL. Pearson correlation coefficient analysis allowed us to identify two appropriate order parameters describing the unfolding transition. The crossover in the rate can be explained in terms of an interplay between entropy and enthalpy with relative importance changing from low force to high force. We rationalize the results by using multidimensional transition-state theory.

Experimental evidence of symmetry breaking of transition-path times

While thermal rates of state transitions in classical systems have been studied for almost a century, associated transition-path times have only recently received attention. Uphill and downhill transition paths between states at different free energies should be statistically indistinguishable. Here, we systematically investigate transition-path-time symmetry and report evidence of its breakdown on the molecular- and meso-scale out of equilibrium. In automated Brownian dynamics experiments, we establish first-passage-time symmetries of colloids driven by femtoNewton forces in holographically-created optical landscapes confined within microchannels. Conversely, we show that transitions which couple in a path-dependent manner to fluctuating forces exhibit asymmetry. We reproduce this asymmetry in folding transitions of DNA-hairpins driven out of equilibrium and suggest a topological mechanism of symmetry breakdown. Our results are relevant to measurements that capture a single coordinate in a multidimensional free energy landscape, as encountered in electrophysiology and single-molecule fluorescence experiments.

Nonnative contact effects in protein folding.

A comprehensive understanding of protein folding includes the knowledge of the formation of individual secondary structures, tertiary structure, and the effects of non-native contacts on these folding events. The measurement of these microscopic events has been posing challenges for experiment and molecular simulation. In this work, we performed enhanced sampling MD simulations for three proteins (NTL9, NuG2b, and CspA) and analyzed minimum free energy paths on multi-dimensional free energy landscapes to explore the underlying folding mechanisms. Consistencies can be seen between the present simulations and the existing experiments as well as other MD simulations. Quantitative analysis reveals the nucleation-condensation folding mechanism indicating the concurrent build-up of secondary and tertiary structures for the three proteins and gives the detailed formation sequence of individual native secondary structure elements. More importantly, nonnative contacts are generally observed among the proteins, creating a nonnative environment to affect the folding of individual secondary structure elements. A general tendency is that the secondary structure element(s) where the maximal nonnative contacts are observed have the largest formation free-energy barrier(s), corresponding to the rate-limiting step(s) of the folding for proteins that follow the nucleation-condensation mechanism. In summary, while native contacts determine the folding mechanism and pathway, non-native contacts play an important role in determining the protein folding thermodynamics by influencing the free energies of individual secondary structure element formation.

Riemann metric approach to optimal sampling of multidimensional free-energy landscapes.

Exploring the free-energy landscape along reaction coordinates or system parameters λ is central to many studies of high-dimensional model systems in physics, e.g., large molecules or spin glasses. In simulations this usually requires sampling conformational transitions or phase transitions, but efficient sampling is often difficult to attain due to the roughness of the energy landscape. For Boltzmann distributions, crossing rates decrease exponentially with free-energy barrier heights. Thus, exponential acceleration can be achieved in simulations by applying an artificial bias along λ tuned such that a flat target distribution is obtained. A flat distribution is, however, an ambiguous concept unless a proper metric is used and is generally suboptimal. Here we propose a multidimensional Riemann metric, which takes the local diffusion into account, and redefine uniform sampling such that it is invariant under nonlinear coordinate transformations. We use the metric in combination with the accelerated weight histogram method, a free-energy calculation and sampling method, to adaptively optimize sampling toward the target distribution prescribed by the metric. We demonstrate that for complex problems, such as molecular dynamics simulations of DNA base-pair opening, sampling uniformly according to the metric, which can be calculated without significant computational overhead, improves sampling efficiency by 50%-70%.

Mapping the conformational free energy of aspartic acid in the gas phase and in aqueous solution

The conformational free energy landscape of aspartic acid, a proteogenic amino acid involved in a wide variety of biological functions, was investigated as an example of the complexity that multiple rotatable bonds produce even in relatively simple molecules. To efficiently explore such a landscape, this molecule was studied in the neutral and zwitterionic forms, in the gas phase and in water solution, by means of molecular dynamics and the enhanced sampling method metadynamics with classical force-fields. Multi-dimensional free energy landscapes were reduced to bi-dimensional maps through the non-linear dimensionality reduction algorithm sketch-map to identify the energetically stable conformers and their interconnection paths. Quantum chemical calculations were then performed on the minimum free energy structures. Our procedure returned the low energy conformations observed experimentally in the gas phase with rotational spectroscopy [M. E. Sanz et al., Phys. Chem. Chem. Phys. 12, 3573 (2010)]. Moreover, it provided information on higher energy conformers not accessible to experiments and on the conformers in water. The comparison between different force-fields and quantum chemical data highlighted the importance of the underlying potential energy surface to accurately capture energy rankings. The combination of force-field based metadynamics, sketch-map analysis, and quantum chemical calculations was able to produce an exhaustive conformational exploration in a range of significant free energies that complements the experimental data. Similar protocols can be applied to larger peptides with complex conformational landscapes and would greatly benefit from the next generation of accurate force-fields.

Determining Protein Folding Pathway and Associated Energetics through Partitioned Integrated-Tempering-Sampling Simulation.

Replica exchange molecular dynamics (REMD) and integrated-tempering-sampling (ITS) are two representative enhanced sampling methods which utilize parallel and integrated tempering approaches, respectively. In this work, a partitioned integrated-tempering-sampling (P-ITS) method is proposed which takes advantage of the benefits of both parallel and integrated tempering approaches. Using P-ITS, the folding pathways of a series of proteins with diverse native structures are explored on multidimensional free-energy landscapes, and the associated thermodynamics are evaluated. In comparison to the original form of ITS, P-ITS improves the sampling efficiency and measures the folding/unfolding thermodynamic quantities more consistently with experimental data. In comparison to REMD, P-ITS significantly reduces the requirement of computational resources and meanwhile achieves similar simulation results. The observed structural characterizations of transition and intermediate states of the proteins under study are in good agreement with previous experimental and simulation studies on the same proteins and homologues. Therefore, the P-ITS method has great potential in simulating the structural dynamics of complex biomolecular systems.