An analysis of the influence of protein intrinsic dynamical properties on its thermal unfolding behavior

The influence of the protein topology-encoded dynamical properties on its thermal unfolding motions was studied in the present work. The intrinsic dynamics of protein topology was obtained by the anisotropic network model (ANM). The ANM has been largely used to investigate protein collective functional motions, but it is not well elucidated if this model can also reveal the preferred large-scale motions during protein unfolding. A small protein barnase is used as a typical case study to explore the relationship between protein topology-encoded dynamics and its unfolding motions. Three thermal unfolding simulations at 500 K were performed for barnase and the entire unfolding trajectories were sampled and partitioned into several windows. For each window, the preferred unfolding motions were investigated by essential dynamics analysis, and then associated with the intrinsic dynamical properties of the starting conformation in this window, which is detected by ANM. The results show that only a few slow normal modes imposed by protein structure are sufficient to give a significant overlap with the preferred unfolding motions. Especially, the large amplitude unfolding movements, which imply that the protein jumps out of a local energy basin, can be well described by a single or several ANM slow modes. Besides the global motions, it is also found that the local residual fluctuations encoded in protein structure are highly correlated with those in the protein unfolding process. Furthermore, we also investigated the relationship between protein intrinsic flexibility and its unfolding events. The results show that the intrinsic flexible regions tend to unfold early. Several early unfolding events can be predicted by analysis of protein structural flexibility. These results imply that protein structure-encoded dynamical properties have significant influences on protein unfolding motions.     

Increasing temperature accelerates protein unfolding without changing the pathway of unfolding

We have traditionally relied on extremely elevated temperatures (498K, 225 degrees C) to investigate the unfolding process of proteins within the timescale available to molecular dynamics simulations with explicit solvent. However, recent advances in computer hardware have allowed us to extend our thermal denaturation studies to much lower temperatures. Here we describe the results of simulations of chymotrypsin inhibitor 2 at seven temperatures, ranging from 298K to 498K. The simulation lengths vary from 94ns to 20ns, for a total simulation time of 344ns, or 0.34 micros. At 298K, the protein is very stable over the full 50ns simulation. At 348K, corresponding to the experimentally observed melting temperature of CI2, the protein unfolds over the first 25ns, explores partially unfolded conformations for 20ns, and then refolds over the last 35ns. Above its melting temperature, complete thermal denaturation occurs in an activated process. Early unfolding is characterized by sliding or breathing motions in the protein core, leading to an unfolding transition state with a weakened core and some loss of secondary structure. After the unfolding transition, the core contacts are rapidly lost as the protein passes on to the fully denatured ensemble. While the overall character and order of events in the unfolding process are well conserved across temperatures, there are substantial differences in the timescales over which these events take place. We conclude that 498K simulations are suitable for elucidating the details of protein unfolding at a minimum of computational expense.     

Molecular dynamics simulations of peptides and proteins with amplified collective motions

We present a novel method that uses the collective modes obtained with a coarse-grained model/anisotropic network model to guide the atomic-level simulations. Based on this model, local collective modes can be calculated according to a single configuration in the conformational space of the protein. In the molecular dynamics simulations, the motions along the slowest few modes are coupled to a higher temperature by the weak coupling method to amplify the collective motions. This amplified-collective-motion (ACM) method is applied to two test systems. One is an S-peptide analog. We realized the refolding of the denatured peptide in eight simulations out of 10 using the method. The other system is bacteriophage T4 lysozyme. Much more extensive domain motions between the N-terminal and C-terminal domain of T4 lysozyme are observed in the ACM simulation compared to a conventional simulation. The ACM method allows for extensive sampling in conformational space while still restricting the sampled configurations within low energy areas. The method can be applied in both explicit and implicit solvent simulations, and may be further applied to important biological problems, such as long timescale functional motions, protein folding/unfolding, and structure prediction.     

Dimerization affects collective dynamics of triosephosphate isomerase

Molecular dynamics simulations (30-60 ns runs) are performed on free/apo triosephosphate isomerase (TIM) to determine any correlation between collective motions and loop 6 dynamics. Native TIM is reported to be active only as a homodimer even though cooperativity has not been observed between the two identical subunits. Both dimeric and monomeric (isolated from dimer) forms of TIM are simulated in explicit water at 300 K and 1 bar to inspect any differences between the structures in terms of fluctuation dynamics and functionally important loop 6 dynamics/closure. Significant cross-correlations between residue fluctuations are observed in the dimer, which result from the global counter-rotations of the two identical subunits in the essential modes of the dimer. Specifically, the first essential mode contributing to 34% of overall motion of the dimer is strongly coupled to the loop 6's closure over the active site. In contrast, such significant correlations cannot be observed in the monomeric structure, which maintains relatively localized motions of the loops in the essential modes. Thus, the onset of collective motions at ns time scale due to dimerization has functional implications as to the coordination of loop 6 closure.     

The essential dynamics of Cu, Zn superoxide dismutase: suggestion of intersubunit communication.

A 300-ps molecular dynamics simulation of the whole Cu, Zn superoxide dismutase dimer has been carried out in water, and the trajectory has been analyzed by the essential dynamics method. The results indicate that the motion is defined by few preferred directions identified by the first four to six eigenvectors and that the motion of the two monomers at each instant is not symmetrical. The vectors symmetrical to the eigenvectors are significantly sampled, suggesting that, on average, the motions of the two subunits will exchange. Large intra- and intersubunit motions involving different subdomains of the protein are observed. A mechanical coupling between the two subunits is also suggested, because displacements of the loops surrounding the active site in one monomer are correlated with the motion of parts of the second toward the intersubunit interface.     

The consistency of large concerted motions in proteins in molecular dynamics simulations.

A detailed investigation is presented into the effect of limited sampling time and small changes in the force field on molecular dynamics simulations of a protein. Thirteen independent simulations of the B1 IgG-binding domain of streptococcal protein G were performed, with small changes in the simulation parameters in each simulation. Parameters studied included temperature, bond constraints, cut-off radius for electrostatic interactions, and initial placement of hydrogen atoms. The essential dynamics technique was used to reveal dynamic differences between the simulations. Similar essential dynamics properties were found for all simulations, indicating that the large concerted motions found in the simulations are not particularly sensitive to small changes in the force field. A thorough investigation into the stability of the essential dynamics properties as derived from a molecular dynamics simulation of a few hundred picoseconds is provided. Although the definition of the essential modes of motion has not fully converged in these short simulations, the subspace in which these modes are confined is found to be reproducible.     

Molecular dynamics study of a hyperthermophilic and a mesophilic rubredoxin

In recent years, increased interest in the origin of protein thermal stability has gained attention both for its possible role in understanding the forces governing the folding of a protein and for the design of new highly stable engineered biocatalysts. To study the origin of thermostability, we have performed molecular dynamics simulations of two rubredoxins, from the mesophile Clostridium pasteurianum and from the hyperthermophile Pyrococcus furiosus. The simulations were carried out at two temperatures, 300 and 373 K, for each molecule. The length of the simulations was within the range of 6-7.2 ns. The rubredoxin from the hyperthermophilic organism was more flexible than its mesophilic counterpart at both temperatures; however, the overall flexibility of both molecules at their optimal growth temperature was the same, despite 59% sequence homology. The conformational space sampled by both molecules was larger at 300 K than at 373 K. The essential dynamics analysis showed that the principal overall motions of the two molecules are significantly different. On the contrary, each molecule showed similar directions of motion at both temperatures.     

Protein Collective Motions Coupled to Ligand Migration in Myoglobin

Ligand migration processes inside myoglobin and protein dynamics coupled to the migration were theoretically investigated with molecular dynamics simulations. Based on a linear response theory, we identified protein motions coupled to the transient migration of ligand, carbon monoxide (CO), through channels. The result indicates that the coupled protein motions involve collective motions extended over the entire protein correlated with local gating motions at the channels. Protein motions, coupled to opening of a channel from the distal pocket to a neighboring xenon site, were found to share the collective motion with experimentally observed protein motions coupled to a doming motion of the heme Fe atom upon photodissociation of the ligand. Analysis based on generalized Langevin dynamics elucidated slow and diffusive features of the protein response motions. Remarkably small transmission coefficients for rates of the CO migrations through myoglobin were found, suggesting that the CO migration dynamics are characterized as motions governed by the protein dynamics involving the collective motions, rather than as thermally activated transitions across energy barriers of well-structured channels.     

Local feature analysis: a statistical theory for reproducible essential dynamics of large macromolecules

Multivariate statistical methods are widely used to extract functional collective motions from macromolecular molecular dynamics (MD) simulations. In principal component analysis (PCA), a covariance matrix of positional fluctuations is diagonalized to obtain orthogonal eigenvectors and corresponding eigenvalues. The first few eigenvectors usually correspond to collective modes that approximate the functional motions in the protein. However, PCA representations are globally coherent by definition and, for a large biomolecular system, do not converge on the time scales accessible to MD. Also, the forced orthogonalization of modes leads to complex dependencies that are not necessarily consistent with the symmetry of biological macromolecules and assemblies. Here, we describe for the first time the application of local feature analysis (LFA) to construct a topographic representation of functional dynamics in terms of local features. The LFA representations are low dimensional, and like PCA provide a reduced basis set for collective motions, but they are sparsely distributed and spatially localized. This yields a more reliable assignment of essential dynamics modes across different MD time windows. Also, the intrinsic dynamics of local domains is more extensively sampled than that of globally coherent PCA modes.     

Structural and dynamic effects of alpha-helix deletion in Sso7d: implications for protein thermal stability

Sso7d is a 62-residue protein from the hyperthemophilic archaeon Sulfolobus solfataricus with a denaturation temperature close to 100 degrees C around neutral pH. An engineered form of Sso7d truncated at leucine 54 (L54Delta) is significantly less stable, with a denaturation temperature of 53 degrees C. Molecular dynamics (MD) studies of Sso7d and its truncated form at two different temperatures have been performed. The results of the MD simulations at 300 K indicate that: (1) the flexibility of Sso7d chain at 300 K agrees with that detected from X-ray and NMR structural studies; (2) L54Delta remains stable in the native folded conformation and possesses an overall dynamic behavior similar to that of the parent protein. MD simulations performed at 500 K, 10 ns long, indicate that, while Sso7d is in-silico resistant to high temperature, the truncated variant partially unfolds, revealing the early phases of the thermal unfolding pathway of the protein. Analysis of the trajectories of L54Delta suggests that the unzipping of the N-terminal and C-terminal beta-strands should be the first event of the unfolding pathway, and points out the regions more resistant to thermal unfolding. These findings allow one to understand the role played by specific interactions connecting the two ends of the chain for the high thermal stability of Sso7d, and support recent hypotheses on its folding mechanism emerged from site-directed mutagenesis studies.     

Thermal fluctuation spectroscopy of DNA thermal denaturation

We have developed the technique of thermal fluctuation spectroscopy to measure the thermal fluctuations in a system. This technique is particularly useful to study the denaturation dynamics of biomolecules like DNA. Here we present a study of the thermal fluctuations during the thermal denaturation (or melting) of double-stranded DNA. We find that the thermal denaturation of heteropolymeric DNA is accompanied by large, non-Gaussian thermal fluctuations. The thermal fluctuations show a two-peak structure as a function of temperature. Calculations of enthalpy exchanged show that the first peak comes from the denaturation of AT rich regions and the second peak from denaturation of GC rich regions. The large fluctuations are almost absent in homopolymeric DNA. We suggest that bubble formation and cooperative opening and closing dynamics of basepairs causes the additional fluctuation at the first peak and a large cooperative transition from a partially molten DNA to a completely denatured state causes the additional fluctuation at the second peak.     

Possibility for discriminating between two representative non two-state thermal unfolding models of proteins by DSC

Possible differences between two representative non two-state thermal unfolding mechanisms of protein are discussed concerning differential scanning calorimetry. Numerical simulations showed that, by DSC measurement, it is hard to discriminate between the independent model, which assumes independent unfolding domains in a protein, and the sequential model, which assumes intermediate(s) between native and denatured states, especially when values of molecular weight, denaturation enthalpy, and difference in denaturation temperature of each denaturation process are large. DSC curve analysis of Aspergillus niger glucoamylase based on these two models gave essentially the same thermodynamic parameters.     

Comparative molecular dynamics simulation studies for determining factors contributing to the thermostability of chemotaxis protein “CheY”

Comparative molecular dynamics simulations of chemotaxis protein “CheY” from thermophilic origin Thermotoga maritima and its mesophilic counterpart Salmonella enterica have been performed for 10?ns each at 300 and 350?K, and 20?ns each at 400 and 450?K. The trajectories were analyzed in terms of different factors like root-mean-square deviation, root-mean-square fluctuation, radius of gyration, solvent accessible surface area, H-bonds, salt bridge content, and protein–solvent interactions which indicate distinct differences between the two of them. The two proteins also follow dissimilar unfolding pathways. The overall flexibility calculated by the trace of the diagonalized covariance matrix displays similar flexibility of both the proteins near their optimum growth temperatures. However, at higher temperatures mesophilic protein shows increased overall flexibility than its thermophilic counterpart. Principal component analysis also indicates that the essential subspaces explored by the simulations of two proteins at different temperatures are nonoverlapping and they show significantly different directions of motion. However, there are significant overlaps within the trajectories and similar direction of motions are observed for both proteins at 300?K. Overall, the mesophilic protein leads to increased conformational sampling of the phase space than its thermophilic counterpart. This is the first ever study of thermostability of CheY protein homologs by using protein dynamism as a main impact. Our study might be used as a model for studying the molecular basis of thermostability of two homologous proteins from two organisms living at different temperatures with less visible differences.     

Photochemical Reaction Dynamics of the Primary Event of Vision Studied by Means of a Hybrid Molecular Simulation

The photoisomerization reaction dynamics of a retinal chromophore in the visual receptor rhodopsin was investigated by means of hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations. The photoisomerization reaction of retinal constitutes the primary step of vision and is known as one of the fastest reactions in nature. To elucidate the molecular mechanism of the high efficiency of the reaction, we carried out hybrid ab initio QM/MM MD simulations of the complete reaction process from the vertically excited state to the photoproduct via electronic transition in the entire chromophore-protein complex. An ensemble of reaction trajectories reveal that the excited-state dynamics is dynamically homogeneous and synchronous even in the presence of thermal fluctuation of the protein, giving rise to the very fast formation of the photoproduct. The synchronous nature of the reaction dynamics in rhodopsin is found to originate from weak perturbation of the protein surroundings and from dynamic regulation of volume-conserving motions of the chromophore. The simulations also provide a detailed view of time-dependent modulations of hydrogen-out-of-plane vibrations during the reaction process, and identify molecular motions underlying the experimentally observed dynamic spectral modulations.     

Molecular dynamics simulations as a tool for improving protein stability

Haloalkane dehalogenase (DhlA) was used as a model protein to explore the possibility to use molecular dynamics (MD) simulations as a tool to identify flexible regions in proteins that can serve as a target for stability enhancement by introduction of a disulfide bond. DhlA consists of two domains: an alpha/beta-hydrolase fold main domain and a cap domain composed of five alpha-helices. MD simulations of DhlA showed high mobility in a helix-loop-helix region in the cap domain, involving residues 184-211. A disulfide cross-link was engineered between residue 201 of this flexible region and residue 16 of the main domain. The mutant enzyme showed substantial changes in both thermal and urea denaturation. The oxidized form of the mutant enzyme showed an increase of the apparent transition temperature from 47.5 to 52.5 degrees C, whereas the T(m,app) of the reduced mutant decreased by more than 8 degrees C compared to the wild-type enzyme. Urea denaturation results showed a similar trend. Measurement of the kinetic stability showed that the introduction of the disulfide bond caused a decrease in activation free energy of unfolding of 0.43 kcal mol(-1) compared to the wild-type enzyme and also indicated that the helix-loop-helix region was involved early in the unfolding process. The results show that MD simulations are capable of identifying mobile protein domains that can successfully be used as a target for stability enhancement by the introduction of a disulfide cross-link.     

Interdomain dynamics and ligand binding: molecular dynamics simulations of glutamine binding protein

Periplasmic binding proteins from Gram-negative bacteria possess a common architecture, comprised of two domains linked by a hinge region, a fold which they share with the neurotransmitter-binding domains of ionotropic glutamate receptors (GluRs). Glutamine-binding protein (GlnBP) is one such protein, whose crystal structure has been solved in both open and closed forms. Multi-nanosecond molecular dynamics simulations have been used to explore motions about the hinge region and how they are altered by ligand binding. Glutamine binding is seen to significantly reduce inter-domain motions about the hinge region. Essential dynamics analysis of inter-domain motion revealed the presence of both hinge-bending and twisting motions, as has been reported for a related sugar-binding protein. Significantly, the influence of the ligand on GlnBP dynamics is similar to that previously observed in simulations of rat glutamate receptor (GluR2) ligand-binding domain. The essential dynamics analysis of GlnBP also revealed a third class of motion which suggests a mechanism for signal transmission in GluRs.     

Diffusing and colliding: the atomic level folding/unfolding pathway of a small helical protein

Proteins with ultra-fast folding/unfolding kinetics are excellent candidates for study by molecular dynamics. Here, we describe such simulations of a three helix bundle protein, the engrailed homeodomain (En-HD), which folds via the diffusion-collision model. The unfolding pathway of En-HD was characterized by seven simulations of the protein and 12 simulations of its helical fragments yielding over 1.1 micros of simulation time in water. Various conformational states along the unfolding pathway were identified. There is the compact native-like transition state, a U-shaped helical intermediate and an unfolded state with dynamic helical segments. Each of these states is in good agreement with experimental data. Examining these states as well as the transitions between them, we find the role of long-range tertiary contacts, specifically salt-bridges, important in the folding/unfolding pathway. In the folding direction, charged residues form long-range tertiary contacts before the hydrophobic core is formed. The formation of HII is assisted by a specific salt-bridge and by non-specific (fluctuating) tertiary contacts, which we call contact-assisted helix formation. Salt-bridges persist as the protein approaches the transition state, stabilizing HII until the hydrophobic core is formed. To complement this information, simulations of fragments of En-HD illustrate the helical propensities of the individual segments. By thermal denaturation, HII proved to be the least stable helix, unfolding in less than 450 ps at high temperature. We observed the low helical propensity of C-terminal residues from HIII in fragment simulations which, when compared to En-HD unfolding simulations, link the unraveling of HIII to the initial event that drives the unfolding of En-HD.     

Functionally important correlated motions in the single-stranded DNA-binding protein encoded by filamentous phage Pf3

To elucidate the interplay between different parts of dimeric single-stranded DNA-binding proteins we have studied the correlated motions in the protein encoded by filamentous phage Pf3 via the combined use of 15N-NMR relaxation experiments, molecular dynamics simulations and essential dynamics calculations. These studies provide insight into the mechanism underlying the protein-DNA binding reaction. The most important motions can be described by a few essential modes. Most outstanding is the correlated symmetric motion of the DNA-binding wings, which are far apart in the structure. This motion determines the access of DNA to the DNA-binding domain. A correlation between the motion of the DNA-binding wing and the complex loop is indicated to play a role in the cooperative binding of the protein to DNA. These motions are in the nanosecond regime in correspondence with the 15N-NMR relaxation experiments.     

Picosecond Internal Dynamics of Lysozyme as Affected by Thermal Unfolding in Nonaqueous Environment

A neutron-scattering investigation of the internal picosecond dynamics of lysozyme solvated in glycerol as a function of temperature in the range 200–410 K has been undertaken. The inelastic contribution to the measured intensity is characterized by the presence of a bump generally known as “boson peak”, clearly distinguishable at low temperature. When the temperature is increased the quasielastic component of the spectrum becomes more and more intrusive and progressively overwhelms the vibrational bump. This happens especially for T > 345 K when the protein goes through an unfolding process, which leads to the complete denaturation. The quasielastic term is the superposition of two components whose intensities and linewidths have been studied as a function of temperature. The slower component describes motions with characteristic times of ~4 ps corresponding to reorientations of polypeptide side chains. Both the intensity and linewidth of this kind of relaxations show two distinct regimes with a crossover in the temperature range where the melting process occurs, thus suggesting the presence of a dynamical transition correlated to the protein unfolding. Conversely the faster component might be ascribed to the local dynamics of hydrogen atoms caged by the nearest neighbors with characteristic time of ~0.3 ps.