Methods for molecular dynamics simulations of protein folding/unfolding in solution

All atom molecular dynamics simulations have become a standard method for mapping equilibrium protein dynamics and non-equilibrium events like folding and unfolding. Here, we present detailed methods for performing such simulations. Generic protocols for minimization, solvation, simulation, and analysis derived from previous studies are also presented. As a measure of validation, our water model is compared with experiment. An example of current applications of these methods, simulations of the ultrafast folding protein Engrailed Homeodomain are presented including the experimental evidence used to verify their results. Ultrafast folders are an invaluable tool for studying protein behavior as folding and unfolding events measured by experiment occur on timescales accessible with the high-resolution molecular dynamics methods we describe. Finally, to demonstrate the prospect of these methods for folding proteins, a temperature quench simulation of a thermal unfolding intermediate of the Engrailed Homeodomain is described.     

Elucidation of stable intermediates in urea-induced unfolding pathway of human carbonic anhydrase IX

Human carbonic anhydrase IX (CAIX) has evolved as a promising biomarker for cancer prognosis, due to its overexpression in various cancers and restricted expression in normal tissue. However, limited information is available on its biophysical behavior. The unfolding of CAIX in aqueous urea solution was studied using all-atom molecular dynamics simulation approach. The results of this study revealed a stable intermediate state along the unfolding pathway of CAIX. At intermediate concentrations of urea (2.0–4.0 M), the protein displays a native-like structure with a large population of its secondary structure and hydrophobic contacts remaining intact in addition to small confined overall motions. Beyond 4.0 M urea, the unfolding is more gradual and at 8.0 M urea the structure is largely collapsed due to the solvent effect. The hydrophobic contact analysis suggests that the contact in terminal α-helices is separated initially which propagates in the loss of contacts from centrally located β-sheets. The reduction of 60–65% tertiary contacts in 7.0–8.0 M urea suggested the presence of residual structure in unfolded state and is confirmed with structural snap shot. Free energy landscape analysis suggested that unfolding of CAIX exists through the different intermediate states.     

Effect of external stresses on protein conformation: a computer modelling study

The increasing use of digital technologies such as mobile phones has led to major health concerns about the effects of non-ionizing pulsed radiation exposure. We believe that the health implications of exposure to radiation cannot be fully understood without establishing the molecular mechanisms of biological effects of pulsed microwaves. We aim to establish methods for studying the molecular mechanisms of protein structural and energetic changes occurring due to external stresses related to non-ionizing radiation by using a combination of experimental and theoretical approaches. In this paper, we present the results from our fully atomistic simulation study of chemical and thermal stress response of a prototype protein, insulin. We performed a series of molecular dynamics simulations of insulin in solution under equilibrium conditions, under chemical stress (imitated by reducing the disulfide bonds in the protein molecule), and under short-lived thermal stress (imitated by increasing simulation temperature for up to 2 ns). The resultant protein conformational behaviour was analysed for various properties with the aim of establishing analysis routines for classification of protein unfolding pathways and associated molecular mechanisms.Submitted as a record of the 2002 Australian Biophysical Society meeting     

Direct observation of markovian behavior of the mechanical unfolding of individual proteins

Single-molecule force-clamp spectroscopy is a valuable tool to analyze unfolding kinetics of proteins. Previous force-clamp spectroscopy experiments have demonstrated that the mechanical unfolding of ubiquitin deviates from the generally assumed Markovian behavior and involves the features of glassy dynamics. Here we use single molecule force-clamp spectroscopy to study the unfolding kinetics of a computationally designed fast-folding mutant of the small protein GB1, which shares a similar β-grasp fold as ubiquitin. By treating the mechanical unfolding of polyproteins as the superposition of multiple identical Poisson processes, we developed a simple stochastic analysis approach to analyze the dwell time distribution of individual unfolding events in polyprotein unfolding trajectories. Our results unambiguously demonstrate that the mechanical unfolding of NuG2 fulfills all criteria of a memoryless Markovian process. This result, in contrast with the complex mechanical unfolding behaviors observed for ubiquitin, serves as a direct experimental demonstration of the Markovian behavior for the mechanical unfolding of a protein and reveals the complexity of the unfolding dynamics among structurally similar proteins. Furthermore, we extended our method into a robust and efficient pseudo-dwell-time analysis method, which allows one to make full use of all the unfolding events obtained in force-clamp experiments without categorizing the unfolding events. This method enabled us to measure the key parameters characterizing the mechanical unfolding energy landscape of NuG2 with improved precision. We anticipate that the methods demonstrated here will find broad applications in single-molecule force-clamp spectroscopy studies for a wide range of proteins.     

Molecular dynamics simulation reveals a surface salt bridge forming a kinetic trap in unfolding of truncated Staphylococcal nuclease

Surface salt bridges are ubiquitous in globular proteins. Their contribution to protein stability has been extensively debated in the past decade. Here, molecular dynamics simulations are performed starting from a non-equilibrium state of Staphylococcal nuclease (SNase) with C-terminal truncation (SNaseDelta). The results indicate a key role in the unfolding of the surface salt bridge between arginine 105 and glutamate 135. Experimentally, SNaseDelta is known to be partially unfolded. However, in simulations over 1 ns at 300 K and over 500 ps at 400 K, SNaseDelta remains stable in the native-like folded conformation, the salt bridge hindering unfolding. When the potential function is altered so as to selectively weaken the salt bridge, which then breaks rapidly at 430 K, the protein starts to unfold. The results suggest that breaking of this salt bridge presents a significant barrier to the unfolding transition of SNaseDelta from a native-like state to the unfolded state. Potential of mean force calculations indicate that the barrier height for this transition is approximately 7 kcal/mol.     

Molecular dynamics simulations of lysozyme–lipid systems: probing the early steps of protein aggregation

Using the molecular dynamics simulation, the role of lipids in the lysozyme transition into the aggregation-competent conformation has been clarified. Analysis of the changes of lysozyme secondary structure upon its interactions with the model bilayer membranes composed of phosphatidylcholine and its mixtures with phosphatidylglycerol (10, 40, and 80 mol%) within the time interval of 100 ns showed that lipid-bound protein is characterized by the increased content of β-structures. Along with this, the formation of protein–lipid complexes was accompanied by the increase in the gyration radius and the decrease in RMSD of polypeptide chain. The results obtained were interpreted in terms of the partial unfolding of lysozyme molecule on the lipid matrix, with the magnitude of this effect being increased with increasing the fraction of anionic lipids. Based on the results of molecular dynamics simulation, a hypothetical model of the nucleation of lysozyme amyloid fibrils in a membrane environment was suggested.     

Non-equilibrium molecular dynamics study of nanoscale thermal contact resistance

Interfaces play an important role in microscale and nanoscale heat transfer processes with molecular dynamics (MD) simulations often used to study these interfacial phenomena. In this study, two models were used to simulate thermal conduction across micro contact points and the thermal contact resistance using non-equilibrium molecular dynamics simulations with consideration of the near field radiation. When the ratio of the length of the micro contact to the length of the conduction region is less than 0.125, the influence of the near field radiation should be considered; but when the ratio is larger than 0.2, it can be neglected. When the computational domain sizes are 8.50 × 10.62 × 8.50 nm and 10.62 × 10.62 × 10.62 nm, the MD results show that the thermal contact resistance exponentially increases with decreasing area of the micro contact point and increases with increasing micro contact layer thickness. The MD thermal contact resistances in nanoscale are much larger than that of the classical thermal analysis since the material thermal conductivity reduction is ignored in the classical model. The results also show that material defects increase the thermal resistance.     

Thermal unfolding simulations of apo-calmodulin using leap-dynamics

The simulation method leap-dynamics (LD) has been applied to protein thermal unfolding simulations to investigate domain-specific unfolding behavior. Thermal unfolding simulations of the 148-residue protein apo-calmodulin with implicit solvent were performed at temperatures 290 K, 325 K, and 360 K and compared with the corresponding molecular dynamics trajectories in terms of a number of calculated conformational parameters. The main experimental results of unfolding are reproduced in showing the lower stability of the C-domain: at 290 K, both the N- and C-domains are essentially stable; at 325 K, the C-domain unfolds, whereas the N-domain remains folded; and at 360 K, both domains unfold extensively. This behavior could not be reproduced by molecular dynamics simulations alone under the same conditions. These results show an encouraging degree of convergence between experiment and LD simulation. The simulations are able to describe the overall plasticity of the apo-calmodulin structure and to reveal details such as reversible folding/unfolding events within single helices. The results show that by using the combined application of a fast and efficient sampling routine with a detailed molecular dynamics force field, unfolding simulations of proteins at atomic resolution are within the scope of current computational power.     

Molecular dynamics simulation integrating the inhibition kinetics of hydroxysafflor yellow A on α-glucosidase

Inhibition of α-glucosidase has attracted the attention of researchers due to its connection to type-2 diabetes. Hydroxysafflor yellow A (HSYA) extracted from Carthamus tinctorius L. is a natural antioxidant used in traditional Chinese medicine. In this study, the effect of HSYA on α-glucosidase was evaluated using inhibitory kinetics based on the antioxidant properties of HSYA and by performing computational simulation integration methods. HSYA reversibly inhibited α-glucosidase in a competitive inhibition manner and the evaluated kinetic parameters were IC₅₀ = 1.1 ± 0.22 mM and K_i = 1.04 ± 0.23 mM, respectively. The results of spectrofluorimetry showed that the inner hydrophobic regions of α-glucosidase, which are mostly in the active site, were exposed to the surface with increasing HSYA concentrations, indicating that the inactivation of α-glucosidase by HSYA was accompanied by regional unfolding. The molecular dynamics simulations indicated that the four rings of HSYA interact with four residues such as G217, A278, H279, and G280 at the entrance of the active site. Our study provides insight into the inhibition of α-glucosidase and the accompanying structural changes by HSYA. Based on its α-glucosidase-inhibiting effect and its potential as a natural antioxidant, HSYA is a potential agent for treating α-glucosidase-associated type-2 diabetes.     

The kinetic pathway of folding of barnase

To search for folding intermediates, we have examined the folding and unfolding kinetics of wild-type barnase and four representative mutants under a wide range of conditions that span two-state and multi-state kinetics. The choice of mutants and conditions provided in-built controls for artifacts that might distort the interpretation of kinetics, such as the non-linearity of kinetic and equilibrium data with concentration of denaturant. We measured unfolding rate constants over a complete range of denaturant concentration by using by 1H/2H-exchange kinetics under conditions that favour folding, conventional stopped-flow methods at higher denaturant concentrations and continuous flow. Under conditions that favour multi-state kinetics, plots of the rate constants for unfolding against denaturant concentration fitted quantitatively to the equation for three-state kinetics, with a sigmoid component for a change of rate determining step, as did the refolding kinetics. The position of the transition state on the reaction pathway, as measured by solvent exposure (the Tanford beta value) also moved with denaturant concentration, fitting quantitatively to the same equations with a change of rate determining step. The sigmoid behaviour disappeared under conditions that favoured two-state kinetics. Those data combined with direct structural observations and simulation support a minimal reaction pathway for the folding of barnase that involves two detectable folding intermediates. The first intermediate, I(1), is the denatured state under physiological conditions, D(Phys), which has native-like topology, is lower in energy than the random-flight denatured state U and is suggested by molecular dynamics simulation of unfolding to be on-pathway. The second intermediate, I(2), is high energy, and is proven by the change in rate determining step in the unfolding kinetics to be on-pathway. The change in rate determining step in unfolding with structure or environment reflects the change in partitioning of this intermediate to products or starting materials.     

Theoretical characterization of excitation energy transfer in chlorosome light-harvesting antennae from green sulfur bacteria

We present a theoretical study of excitation dynamics in the chlorosome antenna complex of green photosynthetic bacteria based on a recently proposed model for the molecular assembly. Our model for the excitation energy transfer (EET) throughout the antenna combines a stochastic time propagation of the excitonic wave function with molecular dynamics simulations of the supramolecular structure and electronic structure calculations of the excited states. We characterized the optical properties of the chlorosome with absorption, circular dichroism and fluorescence polarization anisotropy decay spectra. The simulation results for the excitation dynamics reveal a detailed picture of the EET in the chlorosome. Coherent energy transfer is significant only for the first 50 fs after the initial excitation, and the wavelike motion of the exciton is completely damped at 100 fs. Characteristic time constants of incoherent energy transfer, subsequently, vary from 1 ps to several tens of ps. We assign the time scales of the EET to specific physical processes by comparing our results with the data obtained from time-resolved spectroscopy experiments.     

A molecular dynamics simulation of SNase and its hydration shell at high temperature and high pressure

Temperature- and pressure-induced unfolding of staphylococcal nuclease (SNase) was studied by Royer, Winter et al. using a variety of experimental techniques (SAXS, FT-IR and fluorescence spectroscopy, DSC, PPC, densimetry). For a more detailed understanding of the underlying mechanistic processes of the different unfolding scenarios, we have carried out a series of molecular dynamics (MD) computer simulations on SNase. We investigated the initial changes of the structure of the protein upon application of pressure (up to 5 kbar) and discuss volumetric and structural differences between the native and pressure pre-denatured state. Additionally, we have obtained the compressibility of the protein and hydration water and compare these data with experimental results. As water plays a crucial role in determining the structure, dynamics and function of proteins, we undertook a detailed analysis of the structure of the interfacial water and the protein-solvent H-bond network as well. Moreover, we report here also MD results on the temperature-induced unfolding of SNase. The time evolution of the protein volume and solvent accessible surface area during thermal unfolding have been investigated, and we present a detailed discussion of the temperature-induced unfolding pathway of SNase in terms of secondary and tertiary structural changes.     

Comparative analysis of thermal unfolding simulations of RNA recognition motifs (RRMs) of TAR DNA-binding protein 43 (TDP-43)

TAR DNA-binding protein 43 (TDP-43) inclusions have been found in Amyotrophic lateral sclerosis (ALS) and several other neurodegenerative diseases. Many studies suggest the involvement of RNA recognition motifs (RRMs) in TDP-43 proteinopathy. To elucidate the structural stability and the unfolding dynamics of RRMs, we have carried out atomistic molecular dynamics simulations at two different temperatures (300 and 500 K). The simulations results indicate that there are distinct structural differences in the unfolding pathway between the two domains and RRM1 unfolds faster than RRM2 in accordance with the lower thermal stability found experimentally. The unfolding behaviors of secondary structures showed that the α-helix was more stable than β-sheet and structural rearrangements of β-sheets results in formation of additional α-helices. At higher temperature, RRM1 exhibit increased overall flexibility and unfolding than RRM2. The temperature-dependent free energy landscapes consist of multiple metastable states stabilized by non-native contacts and hydrogen bonds in RRM2, thus rendering the RRM2 more prone to misfolding. The structural rearrangements of RRM2 could lead to aberrant protein–protein interactions that may account for enhanced aggregation and toxicity of TDP-43. Our analysis, thus identify the structural and thermodynamic characteristics of the RRMs of TDP-43, which will serve to uncover molecular mechanisms and driving forces in TDP-43 misfolding and aggregation.     

Atomistic simulations of liquid–liquid coexistence in confinement: comparison of thermodynamics and kinetics with bulk

We review a few simulation methods and results related to the structure and non-equilibrium dynamics in the coexistence region of immiscible symmetric binary fluids, in bulk as well as under confinement, with special emphasis on the latter. Monte Carlo methods to estimate interfacial tensions for flat and curved interfaces have been discussed. The latter, combined with a thermodynamic integration technique, provides contact angles for coexisting fluids attached to the wall. For such three-phase coexistence, results for the line tension are also presented. For the kinetics of phase separation, various mechanisms and corresponding theoretical expectations have been discussed. A comparative picture between the domain growth in bulk and confinement (including thin-film and semi-infinite geometry) has been presented from molecular dynamics simulations. Applications of finite-size scaling technique have been discussed in both equilibrium and non-equilibrium contexts.     

Subtilase from Beauveria sp.: conformational and functional investigation of unusual stability

Retention of total activity of the subtilisin-like serine protease from Beauveria sp. MTCC 5184 (Bprot) in the vicinity of (1) 3 M GdnHCl for 12 h, (2) 50 % methanol and dimethyl sulfoxide each for 24 h, and (3) proteolytic enzymes (trypsin, chymotrypsin, and proteinase K) for 48 h led to expect the enzyme to be a kinetically stable protein. Also, the structure of the protein was stable at pH 2.0. Biophysical characterization and conformational transitions were monitored using steady-state and time-resolved fluorescence, FTIR, and CD spectroscopy. Single tryptophan in the protein exists as two conformers, in hydrophobic and polar environment. The secondary structure of Bprot was stable in 3 M GdnHCl as seen in far-UV CD spectra. The active fraction of Bprot obtained from size-exclusion chromatography in the presence of GdnHCl (1.0–3.0 M) eluted at reduced retention time. The peak area of inactive or denatured protein with the same retention time as that of native protein increased with increasing concentration of denaturant (1.0–4.0 M GdnHCl). However, the kinetics of GdnHCl-induced unfolding as studied from intrinsic fluorescence revealed k _unf of native protein to be 5.407 × 10^?5 s^?1 and a half-life of 3.56 h. The enzyme is thermodynamically stable in spite of being resistant to the denaturant, which could be due to the effect of GdnHCl imparting rigidity to the active fraction and simultaneously unfolding the partially unfolded protein that exists in equilibrium with the folded active protein. Thermal and pH denaturation of Bprot exhibited interesting structural transitions.     

1H, 13C,and 15N backbone resonance assignments of the porcine pepsin and porcine pepsin complexed with pepstatin

Pepsin is formed as the zymogen, pepsinogen, which includes an additional 44 residue prosegment (PS) on the N-terminus. Upon acidification (pH <3) the PS is removed, yielding active pepsin. The PS is critical to such processes as the initiation of correct folding and protein stability. In the present study, the NMR assignments of the 34.6 kDa native porcine pepsin and porcine pepsin complexed with pepstatin are reported in order to obtain structural information regarding PS-catalyzed protein folding. Such information would contribute to a better understanding of the nature of folding/unfolding energy barrier of pepsin and other aspartic proteases.     

Farseer-NMR: automatic treatment,analysis and plotting of large,multi-variable NMR data

Molecular dynamics play a significant role in how molecules perform their function. A critical method that provides information on dynamics, at the atomic level, is NMR-based relaxation dispersion (RD) experiments. RD experiments have been utilized for understanding multiple biological processes occurring at micro-to-millisecond time, such as enzyme catalysis, molecular recognition, ligand binding and protein folding. Here, we applied the recently developed high-power RD concept to the Carr–Purcell–Meiboom–Gill sequence (extreme CPMG; E-CPMG) for the simultaneous detection of fast and slow dynamics. Using a fast folding protein, gpW, we have shown that previously inaccessible kinetics can be accessed with the improved precision and efficiency of the measurement by using this experiment.     

Computational analysis on conformational dynamics of bone morphogenetic protein-2 (BMP-2)

BMP-2 is widely used for bone regeneration because of its ability to induce osteoblast differentiation and proliferation. The pharmaceutical application of BMP-2 as bone implant makes the studies on stability and conformational dynamics very relevant as proteins are functional only in their native three-dimensional state. Knowing the factors affecting BMP-2 structure becomes essential for designing bone implants activated by BMP-2. In order to explore the influence of temperature and hydration on protein conformation, we have performed the molecular dynamics (MD) simulations at the time scale of 100 ns with two different force fields. We have examined the dynamic behaviour of BMP-2 monomer and dimer in aqueous medium as well as in vacuum at four different temperatures (300, 350, 400 and 450 K). MD simulation of BMP-2 monomer and dimer in water and vacuum environments shows the major contribution of water in structure stabilization. Temperature of the system affects the secondary structure differently in case of monomer and dimer simulation and the dynamics also depends on the environment viz. vacuum and aqueous. Vacuum simulations show very early loss of the major secondary structure content. On the other hand, BMP-2 monomer and dimer in aqueous environment show the unfolding of α-helix with increasing temperature. This unfolded α-helix is converted into β-sheet at 400 K in monomer of BMP-2. Contrary to this, we did not observe β-sheet formation in dimer BMP-2 even at 450 K indicating that monomers are more aggregation prone entity as compared to dimers of BMP-2.