One-dimensional Free Energy Profile Research Articles

Free Energy Simulations of Receptor-Binding Domain Opening of the SARS-CoV-2 Spike Indicate a Barrierless Transition with Slow Conformational Motions.

Infection by sarbecoviruses begins with the attachment of the homotrimeric viral "spike" protein to the angiotensin-converting enzyme 2 receptor on the surface of mammalian cells. This requires one or more receptor-binding domains (RBDs) to be in the open (up) position. Here, we present the results of long molecular dynamics simulations with umbrella sampling (US) to compute a one-dimensional free energy profile of RBD opening/closing and the associated transition times. After ≃3.58μs of simulation time per US window (∼229 μs in total), which was required to approach trajectory decorrelation, the computed free energy profile was found to be without large barriers. This suggests that the RBD diffuses between the open and closed positions without significant energetic hindrance. This interpretation appears consistent with experiments but is at odds with some previous simulations. Modeling the RBD motion as diffusive dynamics along the computed free energy profile, we find that the overall time required for the transition is only about 2 μs, which is 5 orders of magnitude shorter than experimentally measured transition times. We speculate that the most likely reason for the transition time mismatch is our use of very short glycans, which was required to make the simulations performed here feasible. Despite the long simulation times, the final free energy profile is not fully converged with statistical errors of ≃1.16 kcal/mol, which were found to be consistent with the slow time decay in the autocorrelation of the conformational motions of the protein. The simulation lengths that would be required to obtain fully converged results remain unknown, but the present calculations would benefit from at least an order-of-magnitude extension.

Exhaustive classification and systematic free-energy profile study of single-stranded DNA inter-overhang migration.

Migration of a short single-stranded DNA (ssDNA) between DNA overhangs is a basic molecular process that is widely used in dynamic DNA nanotechnology. The migration rate is sensitive to migration gaits, and limits the speed of dynamic DNA systems like DNA nanowalkers and other functional devices. Here, we identify and exhaustively classify all possible inter-overhang migration gaits of a ssDNA into only four categories based on their intrinsic symmetry. Using the oxDNA package, we conduct a systematic computational study for a typical migrator-overhang system to identify the lowest-energy pathway for all four migration categories. The one-dimensional free-energy profile along this pathway allows a parameter-free estimation of migration rates for all the four categories by the first passage time theory plus benchmarking from experimental rates available for one migration category. The obtained rates indicate a big room to improve DNA nanowalkers' speed above 1 μm per minute. The free-energy profile for each migration category possesses distinct and robust symmetric patterns, which largely decide local barriers, trapping states, and thereby a migration's rate-limiting processes and capacity for directional bias. This study thus provides a unified symmetry-based framework to analyze and optimize ssDNA migrations in kinetics, bias capacity, and structural design for better dynamic DNA nanotechnology.

Predicting and Simulating Mutational Effects on Protein Folding Kinetics.

Mutational perturbations of protein structures, i.e., phi-value analysis, are commonly employed to probe the extent of involvement of a particular residue in the rate-determining step(s) of folding. This generally involves the measurement of folding thermodynamic parameters and kinetic rate constants for the wild-type and mutant proteins. While computational approaches have been reasonably successful in understanding and predicting the effect of mutations on folding thermodynamics, it has been challenging to explore the same on kinetics due to confounding structural, energetic, and dynamic factors. Accordingly, the frequent observation of fractional phi-values (mean of ~0.3) has resisted a precise and consistent interpretation. Here, we describe how to construct, parameterize, and employ a simple one-dimensional free energy surface model that is grounded in the basic tenets of the energy landscape theory to predict and simulate the effect of mutations on folding kinetics. As a proof of principle, we simulate one-dimensional free energy profiles of 806 mutations from 24 different proteins employing just the experimental destabilization as input, reproduce the relative unfolding activation free energies with a correlation of 0.91, and show that the mean phi-value of 0.3 essentially corresponds to the extent of stabilization energy gained at the barrier top while folding.

Capturing Transition Paths and Transition States for Conformational Rearrangements in the Ribosome

To reveal the molecular determinants of biological function, one seeks to characterize the interactions that are formed in conformational and chemical transition states. In other words, what interactions govern the molecule’s energy landscape? To accomplish this, it is necessary to determine which degrees of freedom can unambiguously identify each transition state. Here, we perform simulations of large-scale aminoacyl-transfer RNA (aa-tRNA) rearrangements during accommodation on the ribosome and project the dynamics along experimentally accessible atomic distances. From this analysis, we obtain evidence for which coordinates capture the correct number of barrier-crossing events and accurately indicate when the aa-tRNA is on a transition path. Although a commonly used coordinate in single-molecule experiments performs poorly, this study implicates alternative coordinates along which rearrangements are accurately described as diffusive movements across a one-dimensional free-energy profile. From this, we provide the theoretical foundation required for single-molecule techniques to uncover the energy landscape governing aa-tRNA selection by the ribosome.

Energetics of Hydrophilic Protein-Protein Association and the Role of Water.

Hydrophilic protein-protein interfaces constitute a major part of all protein-protein interfaces and are thus of great importance. However, the quantitative characterization of their association is still an ongoing challenge and the driving force behind their association remains poorly characterized. Here, we have addressed the association of hydrophilic proteins and the role of water by means of extensive molecular dynamics simulations in explicit water using three well studied protein complexes; Barnase-Barstar, cytochrome c-cytochrome c peroxidase, and the N-terminal domain of enzyme I-histidine-containing phosphocarrier. The one-dimensional free energy profiles obtained from umbrella sampling simulations are downhill or, in other words, barrierless. Using these one-dimensional free energy profiles, the computed standard free energies of binding are -12.7 ± 1.1 kcal/mol, -9.4 ± 0.7 kcal/mol, and -8.4 ± 1.9 kcal/mol that are in reasonable to very good agreement with the experimental values of -19.6 kcal/mol, -8.8 kcal/mol, and -7.8 kcal/mol. As expected, analysis of the confined water between the hydrophilic complex partners shows that the density and the orientational order parameter deviate noticeably from the bulk values, especially at close separations of the confining proteins.

Multiple barriers in forced rupture of protein complexes

Curvatures in the most probable rupture force (f*) versus log-loading rate (log r(f)) observed in dynamic force spectroscopy (DFS) on biomolecular complexes are interpreted using a one-dimensional free energy profile with multiple barriers or a single barrier with force-dependent transition state. Here, we provide a criterion to select one scenario over another. If the rupture dynamics occurs by crossing a single barrier in a physical free energy profile describing unbinding, the exponent ν, from (1 - f*/f(c))(1/ν) ~ (log r(f)) with f(c) being a critical force in the absence of force, is restricted to 0.5 ≤ ν ≤ 1. For biotin-ligand complexes and leukocyte-associated antigen-1 bound to intercellular adhesion molecules, which display large curvature in the DFS data, fits to experimental data yield ν < 0.5, suggesting that if ligand unbinding is assumed to proceed along one-dimensional pulling coordinate, the dynamics should occur in a energy landscape with multiple-barriers.

Homo-polymers with balanced hydrophobicity translocate through lipid bilayers and enhance local solvent permeability

Recent experimental studies indicate that polymeric structures with a well-adjusted balance of amphiphilic parts may translocate through self-assembled phospholipid bilayers and enhance the passive trans-membrane transport of smaller molecules. Using a coarse grained lattice Monte Carlo model with explicit solvent we investigate self-assembled lipid bilayers interacting with a linear polymer chain under variation of the hydrophobicity of the chain. Here, we focus on the relationship between the chain's hydrophobicity and its translocation behavior through the membrane as well as induced membrane perturbations. We show, that there is an adsorption transition of the polymer at the bilayer interface, where effectively the solvent phase and the tail phase of the bilayer are equally repulsive for the polymer. Close to this adsorption threshold of the polymer both the translocation probability of the polymer as well as the permeability of the membrane with respect to solvent are enhanced significantly. The frequency of polymer translocation events can be understood quantitatively assuming a simple diffusion along a one-dimensional free energy profile, which is controlled by the effective lipophilicity of the chain and the tail-packing in the bilayer's core.

Convergence and error estimation in free energy calculations using the weighted histogram analysis method.

The weighted histogram analysis method (WHAM) has become the standard technique for the analysis of umbrella sampling simulations. In this article, we address the challenges (1) of obtaining fast and accurate solutions of the coupled nonlinear WHAM equations, (2) of quantifying the statistical errors of the resulting free energies, (3) of diagnosing possible systematic errors, and (4) of optimally allocating of the computational resources. Traditionally, the WHAM equations are solved by a fixed-point direct iteration method, despite poor convergence and possible numerical inaccuracies in the solutions. Here, we instead solve the mathematically equivalent problem of maximizing a target likelihood function, by using superlinear numerical optimization algorithms with a significantly faster convergence rate. To estimate the statistical errors in one-dimensional free energy profiles obtained from WHAM, we note that for densely spaced umbrella windows with harmonic biasing potentials, the WHAM free energy profile can be approximated by a coarse-grained free energy obtained by integrating the mean restraining forces. The statistical errors of the coarse-grained free energies can be estimated straightforwardly and then used for the WHAM results. A generalization to multidimensional WHAM is described. We also propose two simple statistical criteria to test the consistency between the histograms of adjacent umbrella windows, which help identify inadequate sampling and hysteresis in the degrees of freedom orthogonal to the reaction coordinate. Together, the estimates of the statistical errors and the diagnostics of inconsistencies in the potentials of mean force provide a basis for the efficient allocation of computational resources in free energy simulations.

Protein dynamics investigated by inherent structure analysis.

Molecular dynamics (MD) simulations provide essential information about the thermodynamics and dynamics of proteins. To construct the free-energy surface from equilibrium trajectories, it is necessary to group the individual snapshots in a meaningful way. The inherent structures (IS) are shown to provide an appropriate discretization of the trajectory and to avoid problems that can arise in clustering algorithms that have been employed previously. The IS-based approach is illustrated with a 30-ns room temperature "native" state MD simulation of a 10-residue peptide in a beta-hairpin conformation. The transitions between the IS are used to construct a configuration space network from which a one-dimensional free-energy profile is obtained with the mincut method. The results demonstrate that the IS approach is useful and that even for this simple system, there exists a nontrivial organization of the native state into several valleys separated by barriers as high as 3 kcal/mol. Further, by introducing a coarse-grained network, it is demonstrated that there are multiple pathways connecting the valleys. This scenario is hidden when the snapshots of the trajectory are used directly with rmsd clustering to compute the free-energy profile. Application of the IS approach to the native state of the PDZ2 signaling domain indicates its utility for the study of biologically relevant systems.

Conformational change path between closed and open forms of C2 domain of coagulation factor V on a two-dimensional free-energy surface

We test a hypothesis that the closed form of the C2 domain of coagulation factor V is more stable than the open form in an aqueous environment using a two-dimensional free-energy calculation with a simple dielectric solvent model. Our result shows that while the free-energy difference between two forms is small, favoring the closed form, a two-dimensional free-energy surface (FES) reveals that a transition state (1.53 kcal/mol) exists between the two conformations. By mapping the one-dimensional order parameter DeltaQ onto the two-dimensional FES, we search the conformational change path with the highest Boltzmann weighting factor between the closed and open form of the factor V C2 domain. The predicted transition path from the closed to open form is not that of simple side chain movements, but instead concerted movements of several loops. We also present a one-dimensional free-energy profile using a collective order parameter, which in a coarse manner locates the energy barriers found on the two-dimensional FES.

Folding Pathway And Free Energy Landscape Of Villin Headpiece Subdomain HP35 Studied by String Method

The string method with swarm of trajectories is applied here to find the most probable folding pathways of the 35-residue villin headpiece subdomain (HP35). The converged pathway, represented by 61 discrete images, fully characterizes the mechanism of HP35 folding. The three alpha helices in HP35 exhibit distinct patterns of formation, and each forms at a different stage in the folding process. The one-dimensional free energy profile is computed by integrating the mean force along the folding pathway and roughly two major energy barriers are observed. The biggest folding energy barrier is estimated to be 4.1 kcal/mol and the second one about 3.3 kcal/mol. The two major free energy barriers divide the whole folding process into three metastable states, namely the unfolded, native and a partially folded state in between, which is characterized by an aromatic core including residue Phe6, Phe10 and Phe17. This three-state picture is consistent with the biphasic kinetics inferred from previous laser temperature jump experiments and a recent computational study using temperature replica exchange MD simulations. Markov states model (MSM) is then built to estimate the rate of folding. We find that folding time scale implied from MSM is much faster than experimental folding time, although good Markovian behavior is observed in present model. Possbile reasons for this difference and alternative ways to improve the MSM are discussed.

Efficient and Direct Generation of Multidimensional Free Energy Surfaces via Adiabatic Dynamics without Coordinate Transformations

Adiabatic free energy dynamics (AFED) was introduced by Rosso et al. [J. Chem. Phys. 2002, 116, 4389] for computing free energy profiles quickly and accurately using a dynamical adiabatic separation between a set of collective variables or reaction coordinates and the remaining degrees of freedom of a system. This approach has been shown to lead to a significant gain in efficiency versus traditional methods such as umbrella sampling, thermodynamic integration, and free energy perturbation for generating one-dimensional free energy profiles. More importantly, AFED is able to generate multidimensional free energy surfaces efficiently via full sweeps of the surface that rapidly map out the locations of the free energy minima. The most significant drawback to the AFED approach is the need to transform the coordinates into a generalized coordinate system that explicitly contains the collective variables of interest. Recently, Maragliano and Vanden-Eijnden built upon the AFED approach by introducing a set of extended phase-space variables, to which the adiabatic decoupling and high temperature are applied [Chem. Phys. Lett. 2006, 426, 168]. In this scheme, which the authors termed "temperature accelerated molecular dynamics" or TAMD, the need for explicit coordinate transformations is circumvented. The ability of AFED and TAMD to generate free energy surfaces efficiently depends on the thermostatting mechanism employed, since both approaches are inherently nonequilibrium due to the adiabatic decoupling. Indeed, Maragliano and Vanden-Eijnden did not report any direct generation of free energy surfaces within the overdamped Langevin dynamics employed by these authors. Here, we show that by formulating TAMD in a manner that is closer to the original AFED approach, including the generalized Gaussian moment thermostat (GGMT) and multiple time-scale integration, multidimensional free energy surfaces for complex systems can be generated directly from the probability distribution function of the extended phase-space variables. The new TAMD formulation, which we term driven AFED or d-AFED, is applied to compare the conformational preferences of small peptides both in gas phase and in solution for three force fields. The results show that d-AFED/TAMD accurately and efficiently generates free energy surfaces in two collective variables useful for characterizing the conformations, namely, the radius of gyration, R(G), and number of hydrogen bonds, N(H).

One-dimensional barrier-preserving free-energy projections of a beta-sheet miniprotein: new insights into the folding process.

The conformational space of a 20-residue three-stranded antiparallel beta-sheet peptide (double hairpin) was sampled by equilibrium folding/unfolding molecular dynamics simulations for a total of 20 micros. The resulting one-dimensional free-energy profiles (FEPs) provide a detailed description of the free-energy basins and barriers for the folding reaction. The similarity of the FEPs obtained using the probability of folding before unfolding (pfold) or the mean first passage time supports the robustness of the procedure. The folded state and the most populated free-energy basins in the denatured state are described by the one-dimensional FEPs, which avoid the overlap of states present in the usual one- or two-dimensional projections. Within the denatured state, a basin with fluctuating helical conformations and a heterogeneous entropic state are populated near the melting temperature at about 11% and 33%, respectively. Folding pathways from the helical basin or enthalpic traps (with only one of the two hairpins formed) reach the native state through the entropic state, which is on-pathway and is separated by a low barrier from the folded state. A simplified equilibrium kinetic network based on the FEPs shows the complexity of the folding reaction and indicates, as augmented by additional analyses, that the basins in the denatured state are connected primarily by the native state. The overall folding kinetics shows single-exponential behavior because barriers between the non-native basins and the folded state have similar heights.

Efficacy of external tetraethylammonium block of the KcsA potassium channel: Molecular and Brownian dynamics studies

Blockade of the KcsA potassium channel by externally applied tetraethylammonium is investigated using molecular dynamics calculations and Brownian dynamics simulations. In KcsA, the aromatic rings of four tyrosine residues located just external to the selectivity filter create an attractive energy well or a binding cage for a tetraethylammonium molecule. We first investigate the effects of re-orienting the four tyrosine residues such that the centers of the aromatic rings face the tetraethylammonium molecule directly. Then, we systematically move the residues inward in both orientations so that the radius of the binding cage formed by them becomes smaller. For each configuration, we construct a one-dimensional free energy profile by bringing in a tetraethylammonium molecule from the external reservoir toward the selectivity filter. The free energy profile is then converted to a one-dimensional potential energy profile, taking the available space between the tyrosine residues and the tetraethylammonium molecule into account. Incorporating this potential energy profile into the Brownian dynamics algorithm, we determine the conductance properties of the channel under various conditions, construct the current-tetraethylammonium-concentration curve and compare it with the experimentally determined inhibitory constant k i for externally applied tetraethylammonium. We show that the experimentally determined binding affinity for externally applied tetraethylammonium can be replicated when each of the four tyrosine residues is moved inward by about 0.7 Å, irrespective of orientation of their aromatic rings.

Brownian Dynamics Theory for Predicting Internal and External Blockages of Tetraethylammonium in the KcsA Potassium Channel

The theory of Brownian dynamics is used to model permeation and the blocking of KcsA potassium channels by tetraethylammonium (TEA). A novel Brownian dynamics simulation algorithm is implemented that comprises two free energy profiles; one profile is seen by the potassium ions and the other by the TEA molecules whose shape is approximated by a sphere. Our simulations reveal that internally applied TEA blocks the passage of K + ions by physically occluding the pore. A TEA molecule in the external reservoir encounters an attractive energy-well created by four tyrosine residues at position 82, in addition to all other attractive and repulsive forces impinging on it. Using Brownian dynamics, we investigate how deep the energy-well needs to be to reproduce the experimentally determined inhibitory constant k i for the TEA blockade of KcsA or the mutant Shaker T449Y. The one-dimensional free energy profile obtained from molecular dynamics is first converted into a one-dimensional potential energy profile, and is then transformed into a three-dimensional free energy profile in Brownian dynamics by adding the short-range potential from the channel walls. When converted, the free energy profile calculated from molecular dynamics gives a well-depth of ∼10 kT. We systematically alter the depths of the profiles, and then use Brownian dynamics simulations to numerically determine the current versus TEA-concentration curves. We show that the sequence of binding and unbinding events of the TEA molecule to the binding pocket can be modeled by a first-order Markov process. The Brownian dynamics simulations also reveal that the probability of a TEA molecule binding to the binding pocket in KcsA potassium channels increases exponentially with TEA concentration and depends also on the applied potential and the K + concentration in the simulation assembly.

One-Dimensional Free-Energy Profiles of Complex Systems: Progress Variables that Preserve the Barriers

We show that the balanced minimum-cut procedure introduced in PNAS 2004, 101, 14766 can be reinterpreted as a method for solving the constrained optimization problem of finding the minimum cut among the cuts with a particular value of an additive function of the nodes on either side of the cut. Such an additive function (e.g., the partition function of the reactant region) can be used as a progress coordinate to determine a one-dimensional profile (FEP) of the free-energy surface of the protein-folding reaction as well as other complex reactions. The algorithm is based on the network (obtained from an equilibrium molecular dynamics simulation) that represents the calculated reaction behavior. The resulting FEP gives the exact values of the free energy as a function of the progress coordinate; i.e., at each value of the progress coordinate, the profile is obtained from the surface with the minimal partition function among the surfaces that divide the full free-energy surface between two chosen end points. In many cases, the balanced minimum-cut procedure gives results for only a limited set of points. An approximate method based on p(fold) is shown to provide the profile for a more complete set of values of the progress coordinate. Applications of the approach to model problems and to realistic systems (beta-hairpin of protein G, LJ38 cluster) are presented.

Fluxes, first passage times, and the reduction of Hill diagrams

The steady-state flux resulting from the coupling of two multistate systems is considered. The dynamics of these systems are described (a) as diffusion along a continuous one-dimensional free-energy profile specified by a conformational coordinate or (b) in terms of transitions between a discrete but arbitrary number of substates. If these multistate systems are connected in a simple way, it is shown that the steady-state flux can be obtained analytically. For both the continuous and discrete cases, the exact flux is shown to be identical to that calculated from a simple kinetic scheme involving only four states, if the effective rate constants of this reduced scheme are appropriately defined in terms of the mean first passage times for moving between various points along the multistate cycles. These results clarify and quantify the manner in which the internal conformational dynamics of two multistate systems influences the steady-state flux.