Combining optimal control theory and molecular dynamics for protein folding

A new method to develop low-energy folding routes for proteins is presented. The novel aspect of the proposed approach is the synergistic use of optimal control theory with Molecular Dynamics (MD). In the first step of the method, optimal control theory is employed to compute the force field and the optimal folding trajectory for the Cα atoms of a Coarse-Grained (CG) protein model. The solution of this CG optimization provides an harmonic approximation of the true potential energy surface around the native state. In the next step CG optimization guides the MD simulation by specifying the optimal target positions for the Cα atoms. In turn, MD simulation provides an all-atom conformation whose Cα positions match closely the reference target positions determined by CG optimization. This is accomplished by Targeted Molecular Dynamics (TMD) which uses a bias potential or harmonic restraint in addition to the usual MD potential. Folding is a dynamical process and as such residues make different contacts during the course of folding. Therefore CG optimization has to be reinitialized and repeated over time to accomodate these important changes. At each sampled folding time, the active contacts among the residues are recalculated based on the all-atom conformation obtained from MD. Using the new set of contacts, the CG potential is updated and the CG optimal trajectory for the Cα atoms is recomputed. This is followed by MD. Implementation of this repetitive CG optimization-MD simulation cycle generates the folding trajectory. Simulations on a model protein Villin demonstrate the utility of the method. Since the method is founded on the general tools of optimal control theory and MD without any restrictions, it is widely applicable to other systems. It can be easily implemented with available MD software packages.     

A molecular dynamics approach for the generation of complete protein structures from limited coordinate data

Generation of full protein coordinates from limited information, e.g., the Cα coordinates, is an important step in protein homology modeling and structure determination, and molecular dynamics (MD) simulations may prove to be important in this task. We describe a new method, in which the protein backbone is built quickly in a rather crude way and then refined by minimization techniques. Subsequently, the side chains are positioned using extensive MD calculations. The method is tested on two proteins, and results compared to proteins constructed using two other MD-based methods. In the first method, we supplemented an existing backbone building method with a new procedure to add side chains. The second one largely consists of available methodology. The constructed proteins are compared to the corresponding X-ray structures, which became available during this study, and they are in good agreement (backbone RMS values of 0.5–0.7 Å, and all-atom RMS values of 1.5–1.9 Å). This comparative study indicates that extensive MD simulations are able, to some extent, to generate details of the native protein structure, and may contribute to the development of a standardized methodology to predict reliably (parts of) protein structures when only partial coordinate data are available. © 1994 John Wiley & Sons, Inc.     

Low-mass molecular dynamics simulation: A simple and generic technique to enhance configurational sampling

CLN025 is one of the smallest fast-folding proteins. Until now it has not been reported that CLN025 can autonomously fold to its native conformation in a classical, all-atom, and isothermal–isobaric molecular dynamics (MD) simulation. This article reports the autonomous and repeated folding of CLN025 from a fully extended backbone conformation to its native conformation in explicit solvent in multiple 500-ns MD simulations at 277 K and 1 atm with the first folding event occurring as early as 66.1 ns. These simulations were accomplished by using AMBER forcefield derivatives with atomic masses reduced by 10-fold on Apple Mac Pros. By contrast, no folding event was observed when the simulations were repeated using the original AMBER forcefields of FF12SB and FF14SB. The results demonstrate that low-mass MD simulation is a simple and generic technique to enhance configurational sampling. This technique may propel autonomous folding of a wide range of miniature proteins in classical, all-atom, and isothermal–isobaric MD simulations performed on commodity computers—an important step forward in quantitative biology.     

Effects of surface water on protein dynamics studied by a novel coarse-grained normal mode approach

Normal mode analysis (NMA) has received much attention as a direct approach to extract the collective motions of macromolecules. However, the stringent requirement of computational resources by classical all-atom NMA limits the size of the macromolecules to which the method is normally applied. We implemented a novel coarse-grained normal mode approach based on partitioning the all-atom Hessian matrix into relevant and nonrelevant parts. It is interesting to note that, using classical all-atom NMA results as a reference, we found that this method generates more accurate results than do other coarse-grained approaches, including elastic network model and block normal mode approaches. Moreover, this new method is effective in incorporating the energetic contributions from the nonrelevant atoms, including surface water molecules, into the coarse-grained protein motions. The importance of such improvements is demonstrated by the effect of surface water to shift vibrational modes to higher frequencies and by an increase in overlap of the coarse-grained eigenvector space (the motion directions) with that obtained from molecular dynamics simulations of solvated protein in a water box. These results not only confirm the quality of our method but also point out the importance of incorporating surface structural water in studying protein dynamics.     

Generalized correlation for biomolecular dynamics

Correlated motions in biomolecules are often essential for their function, e.g., allosteric signal transduction or mechanical/thermodynamic energy transport. Because correlated motions in biomolecules remain difficult to access experimentally, molecular dynamics (MD) simulations are particular useful for their analysis. The established method to quantify correlations from MD simulations via calculation of the covariance matrix, however, is restricted to linear correlations and therefore misses part of the correlations in the atomic fluctuations. Herein, we propose a general statistical mechanics approach to detect and quantify any correlated motion from MD trajectories. This generalized correlation measure is contrasted with correlations obtained from covariance matrices for the B1 domain of protein G and T4 lysozyme. The new method successfully quantifies correlations and provides a valuable global overview over the functionally relevant collective motions of lysozyme. In particular, correlated motions of helix 1 together with the two main lobes of lysozyme are detected, which are not seen by the conventional covariance matrix. Overall, the established method misses more than 50% of the correlation. This failure is attributed to both, an interfering and unnecessary dependence on mutual orientations of the atomic fluctuations and, to a lesser extent, attributed to nonlinear correlations. Our generalized correlation measure overcomes these problems and, moreover, allows for an improved understanding of the conformational dynamics by separating linear and nonlinear contributions of the correlation.     

Side-chain dynamics and protein folding

The processes by which protein side chains reach equilibrium during a folding reaction are investigated using both lattice and all-atom simulations. We find that rates of side-chain relaxation exhibit a distribution over the protein structure, with the fastest relaxing side chains located in positions kinetically important for folding. Traversal of the major folding transition state corresponds to the freezing of a small number of side chains, belonging to the folding nucleus, whereas the rest of the protein proceeds toward equilibrium via backbone fluctuations around the native fold. The postnucleation processes by which side chains relax are characterized by very slow dynamics and many barrier crossings, and thus resemble the behavior of a glass.     

From helix-coil transitions to protein folding

An evolution of procedures to simulate protein structure and folding pathways is described. From an initial focus on the helix-coil transition and on hydrogen-bonding and hydrophobic interactions, our original attempts to determine protein structure and folding pathways were based on an experimental approach. Experiments on the oxidative folding of reduced bovine pancreatic ribonuclease A (RNase A) led to a mechanism by which the molecule folded to the native structure by a minimum of four different pathways. The experiments with RNase A were followed by development of a molecular mechanics approach, first, making use of global optimization procedures and then with molecular dynamics (MD), evolving from an all-atom to a united-residue model. This hierarchical MD approach facilitated probing of the folding trajectory to longer time scales than with all-atom MD, and hence led to the determination of complete folding trajectories, thus far for a protein containing as many as 75 amino acid residues. With increasing refinement of the computational procedures, the computed results are coming closer to experimental observations, providing an understanding as to how physics directs the folding process.     

Functional domain motions in proteins on the ~1-100 ns timescale: comparison of neutron spin-echo spectroscopy of phosphoglycerate kinase with molecular-dynamics simulation

Protein function often requires large-scale domain motion. An exciting new development in the experimental characterization of domain motions in proteins is the application of neutron spin-echo spectroscopy (NSE). NSE directly probes coherent (i.e., pair correlated) scattering on the ~1-100 ns timescale. Here, we report on all-atom molecular-dynamics (MD) simulation of a protein, phosphoglycerate kinase, from which we calculate small-angle neutron scattering (SANS) and NSE scattering properties. The simulation-derived and experimental-solution SANS results are in excellent agreement. The contributions of translational and rotational whole-molecule diffusion to the simulation-derived NSE and potential problems in their estimation are examined. Principal component analysis identifies types of domain motion that dominate the internal motion's contribution to the NSE signal, with the largest being classic hinge bending. The associated free-energy profiles are quasiharmonic and the frictional properties correspond to highly overdamped motion. The amplitudes of the motions derived by MD are smaller than those derived from the experimental analysis, and possible reasons for this difference are discussed. The MD results confirm that a significant component of the NSE arises from internal dynamics. They also demonstrate that the combination of NSE with MD is potentially useful for determining the forms, potentials of mean force, and time dependence of functional domain motions in proteins.     

Insights into internal dynamics of 6-phosphogluconolactonase from Trypanosoma brucei studied by nuclear magnetic resonance and molecular dynamics

Nuclear magnetic resonance is used to investigate the backbone dynamics in 6-phosphogluconolactonase from Trypanosoma brucei (Tb6PGL) with (holo-) and without (apo-) 6-phosphogluconic acid as ligand. Relaxation data were analyzed using the model-free approach and reduced spectral density mapping. Comparison of predictions, based on 77 ns molecular dynamics simulations, with the observed relaxation rates gives insight into dynamical properties of the protein and their alteration on ligand binding. Data indicate dynamics changes in the vicinity of the binding site. More interesting is the presence of perturbations located in remote regions of this well-structured globular protein in which no large-amplitude motions are involved. This suggests that delocalized changes in dynamics that occur upon binding could be a general feature of protein-target interactions.     

Application of biasing‐potential replica‐exchange simulations for loop modeling and refinement of proteins in explicit solvent

Comparative or homology modeling of a target protein based on sequence similarity to a protein with known structure is widely used to provide structural models of proteins. Depending on the target‐template similarity these model structures may contain regions of limited structural accuracy. In principle, molecular dynamics (MD) simulations can be used to refine protein model structures and also to model loop regions that connect structurally conserved regions but it is limited by the currently accessible simulation time scales. A recently developed biasing potential replica exchange (BP‐REMD) method was used to refine loops and complete decoy protein structures at atomic resolution including explicit solvent. In standard REMD simulations several replicas of a system are run in parallel at different temperatures allowing exchanges at preset time intervals. In a BP‐REMD simulation replicas are controlled by various levels of a biasing potential to reduce the energy barriers associated with peptide backbone dihedral transitions. The method requires much fewer replicas for efficient sampling compared with T‐REMD. Application of the approach to several protein loops indicated improved conformational sampling of backbone dihedral angle of loop residues compared to conventional MD simulations. BP‐REMD refinement simulations on several test cases starting from decoy structures deviating significantly from the native structure resulted in final structures in much closer agreement with experiment compared to conventional MD simulations. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Improving the accuracy of protein stability predictions with multistate design using a variety of backbone ensembles

Multistate computational protein design (MSD) with backbone ensembles approximating conformational flexibility can predict higher quality sequences than single‐state design with a single fixed backbone. However, it is currently unclear what characteristics of backbone ensembles are required for the accurate prediction of protein sequence stability. In this study, we aimed to improve the accuracy of protein stability predictions made with MSD by using a variety of backbone ensembles to recapitulate the experimentally measured stability of 85 Streptococcal protein G domain β1 sequences. Ensembles tested here include an NMR ensemble as well as those generated by molecular dynamics (MD) simulations, by Backrub motions, and by PertMin, a new method that we developed involving the perturbation of atomic coordinates followed by energy minimization. MSD with the PertMin ensembles resulted in the most accurate predictions by providing the highest number of stable sequences in the top 25, and by correctly binning sequences as stable or unstable with the highest success rate (≈90%) and the lowest number of false positives. The performance of PertMin ensembles is due to the fact that their members closely resemble the input crystal structure and have low potential energy. Conversely, the NMR ensemble as well as those generated by MD simulations at 500 or 1000 K reduced prediction accuracy due to their low structural similarity to the crystal structure. The ensembles tested herein thus represent on‐ or off‐target models of the native protein fold and could be used in future studies to design for desired properties other than stability. Proteins 2014; 82:771–784. © 2013 Wiley Periodicals, Inc.     

Simulating nanoscale functional motions of biomolecules

We are describing efficient dynamics simulation methods for the characterization of functional motion of biomolecules on the nanometer scale. Multivariate statistical methods are widely used to extract and enhance functional collective motions from molecular dynamics (MD) simulations. A dimension reduction in MD is often realized through a principal component analysis (PCA) or a singular value decomposition (SVD) of the trajectory. Normal mode analysis (NMA) is a related collective coordinate space approach, which involves the decomposition of the motion into vibration modes based on an elastic model. Using the myosin motor protein as an example we describe a hybrid technique termed amplified collective motions (ACM) that enhances sampling of conformational space through a combination of normal modes with atomic level MD. Unfortunately, the forced orthogonalization of modes in collective coordinate space leads to complex dependencies that are not necessarily consistent with the symmetry of biological macromolecules and assemblies. In many biological molecules, such as HIV-1 protease, reflective or rotational symmetries are present that are broken using standard orthogonal basis functions. We present a method to compute the plane of reflective symmetry or the axis of rotational symmetry from the trajectory frames. Moreover, we develop an SVD that best approximates the given trajectory while respecting the symmetry. Finally, we describe a local feature analysis (LFA) to construct a topographic representation of functional dynamics in terms of local features. The LFA representations are low-dimensional, and provide a reduced basis set for collective motions, but unlike global collective modes they are sparsely distributed and spatially localized. This yields a more reliable assignment of essential dynamics modes across different MD time windows.     

Protein folding pathways and state transitions described by classical equations of motion of an elastic network model

Protein topology defined by the matrix of residue contacts has proved to be a fruitful basis for the study of protein dynamics. The widely implemented coarse-grained elastic network model of backbone fluctuations has been used to describe crystallographic temperature factors, allosteric couplings, and some aspects of the folding pathway. In the present study, we develop a model of protein dynamics based on the classical equations of motion of a damped network model (DNM) that describes the folding path from a completely unfolded state to the native conformation through a single-well potential derived purely from the native conformation. The kinetic energy gained through the collapse of the protein chain is dissipated through a friction term in the equations of motion that models the water bath. This approach is completely general and sufficiently fast that it can be applied to large proteins. Folding pathways for various proteins of different classes are described and shown to correlate with experimental observations and molecular dynamics and Monte Carlo simulations. Allosteric transitions between alternative protein structures are also modeled within the DNM through an asymmetric double-well potential.     

The flexibility and dynamics of protein disulfide isomerase

We have studied the mobility of the multidomain folding catalyst, protein disulfide isomerase (PDI), by a coarse‐graining approach based on flexibility. We analyze our simulations of yeast PDI (yPDI) using measures of backbone movement, relative positions and orientations of domains, and distances between functional sites. We find that there is interdomain flexibility at every interdomain junction but these show very different characteristics. The extent of interdomain flexibility is such that yPDI's two active sites can approach much more closely than is found in crystal structures—and indeed hinge motion to bring these sites into proximity is the lowest energy normal mode of motion of the protein. The flexibility predicted for yPDI (based on one structure) includes the other known conformation of yPDI and is consistent with (i) the mobility observed experimentally for mammalian PDI and (ii) molecular dynamics. We also observe intradomain flexibility and clear differences between the domains in their propensity for internal motion. Our results suggest that PDI flexibility enables it to interact with many different partner molecules of widely different sizes and shapes, and highlights considerable similarities of yPDI and mammalian PDI. Proteins 2016; 84:1776–1785. © 2016 Wiley Periodicals, Inc.     

The effects of cosolutes on protein dynamics: the reversal of denaturant-induced protein fluctuations by trimethylamine N-oxide

The protein stabilizing effects of the small molecule osmolyte, trimethylamine N-oxide, against chemical denaturant was investigated by NMR spin-relaxation measurements and model-free analysis. In the presence of 0.7 M guanidine hydrochloride increased picosecond-nanosecond dynamics are observed in the protein ribonuclease A. These increased fluctuations occur throughout the protein, but the most significant increases in flexibility occur at positions believed to be the first to unfold. Addition of 0.35 M trimethylamine N-oxide to this destabilized form of ribonuclease results in significant rigidification of the protein backbone as assessed by (1)H-(15)N order parameters. Statistically, these order parameters are the same as those measured in native ribonuclease indicating that TMAO reduces the amplitude of backbone fluctuations in a destabilized protein. These data suggest that TMAO restricts the bond vector motions on the protein energy landscape to resemble those motions that occur in the native protein and points to a relation between stability and dynamics in this enzyme.     

Mapping backbone dynamics in solution with site-directed spin labeling: GCN4-58 bZip free and bound to DNA

In site-directed spin labeling, a nitroxide-containing side chain is introduced at selected sites in a protein. The EPR spectrum of the labeled protein encodes information about the motion of the nitroxide on the nanosecond time scale, which has contributions from the rotary diffusion of the protein, from internal motions in the side chain, and from backbone fluctuations. In the simplest model for the motion of noninteracting (surface) side chains, the contribution from the internal motion is sequence independent, as is that from protein rotary diffusion. Hence, differences in backbone motions should be revealed by comparing the sequence-dependent motions of nitroxides at structurally homologous sites. To examine this model, nitroxide side chains were introduced, one at a time, along the GCN4-58 bZip sequence, for which NMR (15)N relaxation experiments have identified a striking gradient of backbone mobility along the DNA-binding region [Bracken et al. (1999) J. Mol. Biol. 285, 2133]. Spectral simulation techniques and a simple line width measure were used to extract dynamical parameters from the EPR spectra, and the results reveal a mobility gradient similar to that observed in NMR relaxation, indicating that side chain motions mirror backbone motions. In addition, the sequence-dependent side chain dynamics were analyzed in the DNA/protein complex, which has not been previously investigated by NMR relaxation methods. As anticipated, the backbone motions are damped in the DNA-bound state, although a gradient of motion persists with residues at the DNA-binding site being the most highly ordered, similar to those of helices on globular proteins.     

Protein structural dynamics revealed by site-directed spin labeling and multifrequency EPR

Multifrequency electron paramagnetic resonance (EPR), combined with site-directed spin labeling, is a powerful spectroscopic tool to characterize protein dynamics. The lineshape of an EPR spectrum reflects combined rotational dynamics of the spin probe's local motion within a protein, reorientations of protein domains, and overall protein tumbling. All these motions can be restricted and anisotropic, and separation of these motions is important for thorough characterization of protein dynamics. Multifrequency EPR distinguishes between different motions of a spin-labeled protein, due to the frequency dependence of EPR resolution to fast and slow motion of a spin probe. This gives multifrequency EPR its unique capability to characterize protein dynamics in great detail. In this review, we analyze what makes multifrequency EPR sensitive to different rates of spin probe motion and discuss several examples of its usage to separate spin probe dynamics and overall protein dynamics, to characterize protein backbone dynamics, and to resolve protein conformational states.     

Conformational and dynamics changes induced by bile acids binding to chicken liver bile acid binding protein

The correlation between protein motions and function is a central problem in protein science. Several studies have demonstrated that ligand binding and protein dynamics are strongly correlated in intracellular lipid binding proteins (iLBPs), in which the high degree of flexibility, principally occurring at the level of helix-II, CD, and EF loops (the so-called portal area), is significantly reduced upon ligand binding. We have recently investigated by NMR the dynamic properties of a member of the iLBP family, chicken liver bile acid binding protein (cL-BABP), in its apo and holo form, as a complex with two bile salts molecules. Binding was found to be regulated by a dynamic process and a conformational rearrangement was associated with this event. We report here the results of molecular dynamics (MD) simulations performed on apo and holo cL-BABP with the aim of further characterizing the protein regions involved in motion propagation and of evaluating the main molecular interactions stabilizing bound ligands. Upon binding, the root mean square fluctuation values substantially decrease for CD and EF loops while increase for the helix-loop-helix region, thus indicating that the portal area is the region mostly affected by complex formation. These results nicely correlate with backbone dynamics data derived from NMR experiments. Essential dynamics analysis of the MD trajectories indicates that the major concerted motions involve the three contiguous structural elements of the portal area, which however are dynamically coupled in different ways whether in the presence or in the absence of the ligands. Motions of the EF loop and of the helical region are part of the essential space of both apo and holo-BABP and sample a much wider conformational space in the apo form. Together with NMR results, these data support the view that, in the apo protein, the flexible EF loop visits many conformational states including those typical of the holo state and that the ligand acts stabilizing one of these pre-existing conformations. The present results, in agreement with data reported for other iLBPs, sharpen our knowledge on the binding mechanism for this protein family.