Protein Vibrations and Elastic Network Models

Molecular dynamics (MD) simulations were employed to investigate the structure, dynamics, and local base-pair step deformability of the free 16S ribosomal helix 44 from Thermus thermophilus and of a canonical A-RNA double helix. While helix 44 is bent in the crystal structure of the small ribosomal subunit, the simulated helix 44 is intrinsically straight. It shows, however, substantial instantaneous bends that are isotropic. The spontaneous motions seen in simulations achieve large degrees of bending seen in the X-ray structure and would be entirely sufficient to allow the dynamics of the upper part of helix 44 evidenced by cryo-electron microscopic studies. Analysis of local base-pair step deformability reveals a patch of flexible steps in the upper part of helix 44 and in the area proximal to the bulge bases, suggesting that the upper part of helix 44 has enhanced flexibility. The simulations identify two conformational substates of the second bulge area (bottom part of the helix) with distinct base pairing. In agreement with nuclear magnetic resonance (NMR) and X-ray studies, a flipped out conformational substate of conserved 1492A is seen in the first bulge area. Molecular dynamics (MD) simulations reveal a number of reversible alpha-gamma backbone flips that correspond to transitions between two known A-RNA backbone families. The flipped substates do not cumulate along the trajectory and lead to a modest transient reduction of helical twist with no significant influence on the overall geometry of the duplexes. Despite their considerable flexibility, the simulated structures are very stable with no indication of substantial force field inaccuracies.

The term molecular mechanics (MM) refers to the use of simple potential-energy functions (e.g., harmonic oscillator or Coulombic potentials) to model molecular systems. Molecular mechanics approaches are widely applied in molecular structure refinement, molecular dynamics (MD) simulations, Monte Carlo (MC) simulations, and ligand-docking simulations. Typically, molecular mechanics models consist of spherical atoms connected by springs which represent bonds. Internal forces experienced in the model structure are described using simple mathematical functions. For example, Hooke’s law is commonly used to describe bonded interactions, and the nonbonded atoms might be treated as inelastic hard spheres or may interact according to a Lennard-Jones potential. Using these simple models, a molecular dynamics simulation numerically solves Newton’s equations of motion, thus allowing structural fluctuations to be observed with respect to time. Dynamic simulation methods are widely used to obtain information on the time evolution of conformations of proteins and other biological macromolecules1–4 and also kinetic and thermodynamic information. Simulations can provide fine detail concerning the motions of individual particles as a function of time. They can be utilized to quantify the properties of a system at a precision and on a time scale that is otherwise inaccessible, and simulation is, therefore, a valuable tool in extending our understanding of model systems. Theoretical consideration of a system additionally allows one to investigate the specific contributions to a property through “computational alchemy”,5 that is, modifying the simulation in a way that is nonphysical but nonetheless allows a model’s characteristics to be probed. One particular example is the artificial conversion of the energy function from that representing one system to that of another during a simulation. This is an important technique in free-energy calculations.6 Thus, molecular dynamics simulations, along with a range of complementary computational approaches, have become valuable tools for investigating the basis of protein structure and function. This review offers an outline of the origin of molecular dynamics simulation for protein systems and how it has developed into a robust and trusted tool. This review then covers more recent advances in theory and an illustrative selection of practical studies in which it played a central role. The range of studies in which MD has played a considerable or pivotal role is immense, and this review cannot do justice to them; MD simulations of biomedical importance were recently reviewed.4 Particular emphasis will be placed on the study of dynamic aspects of protein recognition, an area where molecular dynamics has scope to provide broad and far-ranging insights. This review concludes with a brief discussion of the future potential offered to advancement of the biological and biochemical sciences and the remaining issues that must be overcome to allow the full extent of this potential to be realized. 1.1. Historical Background MD methods were originally conceived within the theoretical physics community during the 1950s. In 1957, Alder and Wainwright7 performed the earliest MD simulation using the so-called hard-sphere model, in which the atoms interacted only through perfect collisions. Rahman8 subsequently applied a smooth, continuous potential to mimic real atomic interactions. During the 1970s, as computers became more widespread, MD simulations were developed for more complex systems, culminating in 1976 with the first simulation of a protein9,10 using an empirical energy function constructed using physics-based first-principles assumptions. MD simulations are now widely and routinely applied and especially popular in the fields of materials science11,12 and biophysics. As will be discussed later in this review, a variety of experimental conditions may be simulated with modern theories and algorithms. The initial simulations only considered single molecules in vacuo. Over time, more realistic or at least biologically relevant simulations could be performed. This trend is continuing today. The initial protein MD simulation, of the small bovine pancreatic trypsin inhibitor (BPTI), covered only 9.2 ps of simulation time. Modern simulations routinely have so-called equilibration periods much longer than that, and production simulations of tens of nanoseconds are routine, with the first microsecond MD simulation being reported in 1998.13 In addition, the original BPTI simulation included only about 500 atoms rather than the 104-106 atoms that are common today. While much of this advancement results from an immense increase in availability of computing power, major theoretical and methodological developments also contribute significantly. The number of publications regarding MD theory and application of MD to biological systems is growing at an extraordinary pace. A single review cannot do justice to the recent applications of MD. Using data from ISI Web of Science, the authors estimate that during 2005 at least 800 articles will be published that discuss molecular dynamics and proteins. The historical counts are shown in Figure 1. Open in a separate window Figure 1 Articles matching ISI Web of Science query “TS=(protein) AND TS=(molecular dynamics)”.

Proteins are essential units of life that govern several functions. Understanding their behavior is closely related to their conformations, native folds, and change in conformations. Thus, the dynamic information of protein becomes essential to understand its properties at the molecular level. The molecular dynamics (MD) simulation approach provides atomistic-level dynamic information about proteins. However, more extended or complex MD simulations of protein are challenging to analyze and to gather meaningful confirmation from several snapshots of the dynamic system. To achieve it, i.e., analyzing MD simulation data, Markov State Model (MSM) is a powerful tool that has a statistical background. It represents the MD simulation system as a combination of finite memoryless states, i.e., states that are not dependent on prior states and transition probability among such states. MSM applications have grown from peptides to membrane protein simulations. The present book chapter sheds light on MD simulation’s role in protein dynamics and why MSM is required. The brief theoretical aspects of MSM techniques are demonstrated. Lastly, the chapter discusses the application of MSM in different protein folding and dynamics.

Selectivity and transport in aquaporins from molecular simulation studies.

Extensive molecular dynamics (MD) simulations have been used to investigate the intercalative binding of 9-hydroxyellipticine to the DNA oligonucleotide d(ATATATATATAT)2. Four independent simulations differing in the initial orientation of the drug at the intercalation site were carried out, and compared both with each other and a control simulation of the free DNA sequence. The structure of the latter was compared with structures obtained from x-ray crystallography and nmr spectroscopy, as well as the theoretically derived "alternating B-DNA" model [A. Klug et al. (1979), Journal of Molecular Biology, Vol. 131, p. 669]. The alternation of twist angles observed in experimental structures was reproduced in the simulation. All four independent simulations of the drug-DNA intercalation complex converged in placing the pyridine ring of the ellipticine chromophore in the major groove; in one case this involved a 180 degrees rotation of the drug at the intercalation site. At a more detailed level, the drug is seen to be capable of adopting several distinct orientations, each stable over a period of hundreds of pico-seconds. Despite the presence of several polar groups in the drug, however, no direct hydrogen bonding to the DNA occurs; instead, interactions between the methyl groups of the drug and the thymine bases at the intercalation site appear important in determining the orientational preferences of the drug. Comparison of the intercalation complexes with the free DNA sequence shows a degree of unwinding resulting from intercalation, in good agreement with experimental results, but spread over the three central base-pair steps, not confined to the intercalation site itself. Measurements of torsional rigidity indicate only a slight stiffening of the DNA restricted to the immediate site of intercalation. The structures obtained from the MD simulations were used to calculate theoretical CD spectra, with separate simulations giving very different results. This appears to indicate that given an accurate assignment of the main electronic transition dipole moment of the ellipticine chromophore, discrimination of the more realistic binding geometries may be possible. The relative merits of the various drug orientations observed in the simulations are discussed and a perpendicular orientation of the drug at the intercalation site is considered to be the most consistent with experimental data. While the simulations themselves represent a total of over 2 ns, however, the differences apparent between independent runs indicate that longer simulation times will be required before a complete, unequivocal view of DNA intercalation is obtained.

The mechanism of action of inosine-adenosine-guanosine nucleoside hydrolase (IAG-NH) has been investigated by long-term molecular dynamics (MD) simulation in TIP3P water using stochastic boundary conditions. Special attention has been given to the role of leaving group pocket residues and conformation of the bound substrate at the active site of IAG-NH. We also describe the positioning of the residues of an important flexible loop at the active site, which was previously unobservable by X-ray crystallography due to high B-factors. Five MD simulations have been performed with the Enzyme x Substrate complexes: Enzyme x anti-Adenosine with Asp40-COOH [E(40H) x Ade(a)], Enzyme x anti-Adenosine with Asp40-COO- [E(40) x Ade(a)], Enzyme x syn-Adenosine with Asp40-COOH [E(40H) x Ade(s)], Enzyme x syn-Adenosine with Asp40-COO- [E(40) x Ade(s)], and Enzyme x anti-Inosine with Asp40-COO- [E(40) x Ino(a)]. Overall, the structures generated from the MD simulation of E(40H) x Ade(s) preserve the catalytically important hydrogen bonds as well as electrostatic and hydrophobic interactions to provide a plausible catalytic structure. When deprotonated Asp40 (Asp4-COO-) is present, the active site is open to water solvent which interferes with the base stacking between Trp83 and nucleobase. A calculation using Poisson-Boltzmann equation module supports that Asp40 indeed has an elevated pK(app). Solvent accessible surface area (SASA) calculations on all the five MD structures shows that systems with protonated Asp40, namely, E(40H) x Ade(a) and E(40H) x Ade(s), have zero SASA. It is found that a water molecule is hydrogen-bonded to the N7 of the nucleobase and is probably the essential general acid to protonate the departing nucleobase anion. The N7-bonded water is in turn hydrogen-bonded to waters in a channel, held in place by the residues of the flexible loop, Tyr257, His247, and Cys245. Using normal-mode analysis with elastic network model, we find that the flexible loop explores a conformational space much larger than in the MD trajectory, leading to a "gating"-like motion with respect to the active site.

Biopolymers are versatile molecules that perform myriad cellular functions such as gene regulation and high-energy molecule synthesis. Alterations in biomolecular structure such as those arising from mutations can disrupt important cellular processes and lead to disease. Therefore, accurate biopolymer structure characterization plays an outsized role in improving human health. X-ray crystallography, cryogenic electron microscopy (Cryo-EM), and nuclear magnetic resonance (NMR) spectroscopy are the most powerful techniques for obtaining high-resolution, 3D structures.Valuable structural insight is also obtained using optical rotation, small angle X-ray scattering (SAXS), circular dichroism (CD) spectroscopy, and Forster resonance energy transfer (FRET) spectroscopy. Over the past few years, remarkable advances in computational power and algorithm development have revolutionized the field of biopolymer structure prediction. For example, leveraging machine learning (ML) and artificial intelligence (AI), AlphaFold and Open fold have demonstrated a remarkable ability to predict biopolymer structure from sequence. Although powerful, these techniques are limited in their ability to provide information about challenging systems where the structure is either unknown or highly flexible under physiological conditions. In this work, studies are presented in which structural information is gleaned from challenging biopolymer systems using novel computational and experimental methods. In Chapter 2, a new structure characterization methodology based on hydrogen-deuterium exchange-mass spectrometry (HDX-MS) in combination with molecular dynamics (MD) simulations is introduced. The approach is used to describe the structural motifs of the unstructured peptide, Nt17, a conserved sequence in Huntingtin protein and associated with Huntington’s Disease progression. The model is used to show that, although highly flexible, Nt17 is intrinsically primed for facile conversion to α-helical conformation upon binding with molecular partners. Having shown that HDX reactivity is predictive of the degree of structural flexibility and bias (propensity to form 2° structural elements such as α-helix, β-sheet, and β-turn) for intrinsically disordered regions (IDRs), the work in Chapter 3 explores an extension of structure characterization to protein-ligand systems. In this study human FK506 binding protein (FKBP12) structure is examined in the absence and presence of two, well-recognized ligand molecules (i.e., FK506 and rapamycin). The work provides a better understanding of binding/unbinding events, competitive binding, and folded/unfolded conformer heterogeneity. Native MS and HDX-MS were utilized in conjunction with MD simulations to understand this select intrinsically disordered region (IDR) and IDR-drug conformational landscape. Chapter 4 explores the functional property prediction of Proteorhodopsin, by integrating experimental studies with MD simulations validation. Here, a novel linear regression model accurately predicts the spectral

Anisotropic network model (ANM) is used to analyze the collective motions of restriction enzyme EcoRI in free form and in complex with DNA. For comparison, three independent molecular dynamics (MD) simulations, each of 1.5 ns duration, are also performed for the EcoRI-DNA complex in explicit water. Although high mobility (equilibrium fluctuations) of inner and outer loops that surround the DNA is consistent in both methods and experiments, MD runs sample different conformational subspaces from which reliable collective dynamics cannot be extracted. However, ANM employed on different conformations from MD simulations indicates very similar collective motions. The stems of the inner loops are quite immobile even in the free enzyme and form a large, almost fixed, pocket for DNA binding. As a result, the residues that make specific and non-specific interactions with the DNA exhibit very low fluctuations in the free enzyme. The vibrational entropy difference between the EcoRI complex and free protein + unkinked DNA is positive (favorable), which may partially counteract the unfavorable enthalpy difference of DNA kink formation. Dynamic domains in EcoRI complex and cross-correlations between residue fluctuations indicate possible means of communication between the distal active sites.

Understanding the function of a protein requires not only knowledge of its tertiary structure but also an understanding of its conformational dynamics. Nuclear magnetic resonance (NMR) spectroscopy, polarization-resolved fluorescence spectroscopy and molecular dynamics (MD) simulations are powerful methods to provide detailed insight into protein dynamics on multiple time scales by monitoring global rotational diffusion and local flexibility (order parameters) that are sensitive to inter- and intramolecular interactions, respectively. We present an integrated approach where data from these techniques are analyzed and interpreted within a joint theoretical description of depolarization and diffusion, demonstrating their conceptual similarities. This integrated approach is then applied to the autophagy-related protein GABARAP in its cytosolic form, elucidating its dynamics on the pico- to nanosecond time scale and its rotational and translational diffusion for protein concentrations spanning 9 orders of magnitude. We compare the dynamics of GABARAP as monitored by 15N spin relaxation of the backbone amide groups, fluorescence anisotropy decays and fluorescence correlation spectroscopy of side chains labeled with BODIPY FL, and molecular movies of the protein from MD simulations. The recovered parameters agree very well between the distinct techniques if the different measurement conditions (probe localization, sample concentration) are taken into account. Moreover, we propose a method that compares the order parameters of the backbone and side chains to identify potential hinges for large-scale, functionally relevant intradomain motions, such as residues 27/28 at the interface between the two subdomains of GABARAP. In conclusion, the integrated concept of cross-fertilizing techniques presented here is fundamental to obtaining a comprehensive quantitative picture of multiscale protein dynamics and solvation. The possibility to employ these validated techniques under cellular conditions and combine them with fluorescence imaging opens up the perspective of studying the functional dynamics of GABARAP or other proteins in live cells.

This study examines the thermodynamic and conformational dynamics of GTRNase and lysozyme proteins under varying temperature conditions using molecular dynamics (MD) simulations. The objective is to evaluate their structural stability, folding behavior, and thermodynamic properties to understand their responses to thermal fluctuations. Protein structures were retrieved from the Protein Data Bank (PDB), refined to remove extraneous molecules, and simulated using GROMACS 2022.2 under NVT and NPT ensembles, spanning temperatures from 270 K to 380 K. The results revealed distinct behaviors for the two proteins. GTRNase exhibited a slight escalation in radius of gyration (Rg) from 1.434 nm at 270 K to 1.445 nm at 380 K, suggesting a marginal conformational expansion. In contrast, lysozyme maintained a consistent Rg of 1.38 nm over the same temperature range, indicating structural compactness. Root mean square deviation (RMSD) data demonstrated increased flexibility in both proteins, with GTRNase escalating from 0.115 nm at 270 K to 0.179 nm at 380 K and lysozyme rising from 0.102 nm to 0.142 nm across the temperature range. Solvent-accessible surface area (SASA) for GTRNase fluctuated between 69 nm² and 73 nm², with the lowest value observed at 300 K and the highest at 370 K. These findings highlight that GTRNase is more susceptible to thermal perturbations than lysozyme, showing greater conformational flexibility and expansion. This research underscores the utility of MD simulations in exploring protein behavior and provides valuable insights for applications in protein engineering and drug design.

Protein dynamics is a manifestation of the complex trajectories of these biomolecules on a multidimensional rugged potential energy surface (PES) driven by thermal energy. At present, computational methods such as atomistic molecular dynamics (MD) simulations can describe thermal protein conformational changes in fully solvated environments over millisecond timescales. Despite these advances, a quantitative assessment of protein dynamics remains a complicated topic, intricately linked to issues such as sampling convergence and the identification of appropriate reaction coordinates/structural features to describe protein conformational states and motions. Here, we present the cumulative variance of atomic coordinate fluctuations (CVCF) along trajectories as an intuitive PES sensitive metric to assess both the extent of sampling and protein dynamics captured in MD simulations. We first examine the sampling problem in model one- (1D) and two-dimensional (2D) PES to demonstrate that the CVCF when traced as a function of the sampling variable (time in MD simulations) can identify local and global equilibria. Further, even far from global equilibrium, a situation representative of standard MD trajectories of proteins, the CVCF can distinguish different PES and therefore resolve the resultant protein dynamics. We demonstrate the utility of our CVCF analysis by applying it to distinguish the dynamics of structurally homologous proteins from the ubiquitin family (ubiquitin, SUMO1, SUMO2) and ubiquitin protein-protein interactions. Our CVCF analysis reveals that differential side-chain dynamics from the structured part of the protein (the conserved β-grasp fold) present distinct protein PES to distinguish ubiquitin from SUMO isoforms. Upon binding to two functionally distinct protein partners (UBCH5A and UEV), intrinsic ubiquitin dynamics changes to reflect the binding context even though the two proteins have similar binding modes, which lead to negligible (sub-angstrom scale) structural changes.

Calcium-dependent lectin I from Pseudomonas aeruginosa (PA-IL) binds specifically to oligosaccharides presenting an alpha-galactose residue at their nonreducing end, such as the disaccharides alphaGal1-2betaGalOMe, alphaGal1-3betaGalOMe, and alphaGal1-4betaGalOMe. This provides a unique model for studying the effect of the glycosidic linkage of the ligands on structure and thermodynamics of the complexes by means of experimental and theoretical tools. The structural features of PA-IL in complex with the three disaccharides were established by docking and molecular dynamics simulations and compared with those observed in available crystal structures, including PA-IL.alphaGal1-2betaGalOMe complex, which was solved at 2.4 A resolution and reported herein. The role of a structural bridge water molecule in the binding site of PA-IL was also elucidated through molecular dynamics simulations and free energy calculations. This water molecule establishes three very stable hydrogen bonds with O6 of nonreducing galactose, oxygen from Pro-51 main chain, and nitrogen from Gln-53 main chain of the lectin binding site. Binding free energies for PA-IL in complex with the three disaccharides were investigated, and the results were compared with the experimental data determined by titration microcalorimetry. When the bridge water molecule was included in the free energy calculations, the simulations predicted the correct binding affinity trends with the 1-2-linked disaccharide presenting three times stronger affinity ligand than the other two. These results highlight the role of the water molecule in the binding site of PA-IL and indicate that it should be taken into account when designing glycoderivatives active against P. aeruginosa adhesion.

Dynamical structures of intrinsically disordered proteins (IDPs) and multi-domain proteins that include large ID regions between the domains are unable to be determined by such conventional methods as X-ray crystallography and electron microscopy. Small-angle X-ray scattering (SAXS) is suitable to determine low-resolution structures of proteins and protein complexes in solution, but the structural data on protein dynamics are averaged over the structural ensemble in protein solution. To overcome this problem, we have developed a novel method, named MD-SAXS, of the combined use of SAXS and molecular dynamics (MD) simulation to analyze protein dynamics in solution of multi-subunit protein complexes and multi-domain proteins toward the structural characterization of IDPs. Here we show validity of the method through the structural characterization of restriction Endonuclease EcoO109I.

Human estrogen receptor alfa (ER alfa) is one of the most studied targets for in silico sc reening of bioactive compounds. The estrogenic activity of a vast number of chemicals has been studied for their potentially adverse effects on the hormone regulation of the endocrine system. The commonly accepted presentation of the ER alfa agonist pharmacophore includes terminal phenolic groups and a hydrophobic core with a rigid framework. In this study we report on molecular dynamic (MD) simulations of ER alfa to get a deeper structural insight into the agonist-receptor interactions and the pharmacophore pattern of compounds with agonistic activity. We rely on a crystallographic structure of a complex of ER alfa with an agonist of picomolar affinity . As the X-ray structure has mutation in the key structural element for ER alfa agonistic activity (helix 12, Y537S), a series of MD simulations have been performed on the wild type receptor to prove the stability of the agonist-receptor interactions. The results suggest that the pharmacophore of compounds with ER alfa a gonistic activity can be extended by a feature that occupies a free hydrophobic region of the binding pocket. The results imply also that MD simulations are a powerful in silico tool for both protein dynamics and structure investigation, especially when mutations are available that can potentially disturb the protein structure and functions.

Molecular dynamics (MD) simulations have been performed on adipocyte lipid-binding protein, using the apo and holo forms, bound with stearic and oleic acid. The contribution of electrostatics to protein dynamics and ligand stabilization was assayed by perturbing the electrostatic charge of Arg106 and Arg126 (positive-->neutral) and the fatty acid (132H) headgroup (negative-->neutral). MD simulations for charged holo forms demonstrated significantly greater electrostatic binding energy and a more stabilized hydrogen bond network than simulations performed using neutral forms. Electrostatics, however, appeared to have little effect on fatty acid behavior, e.g., fluctuation of the dihedral head group; number of dihedral transitions within the acyl chain; and change in the end-to-end distance for fatty acid. Instead, fatty acid behavior appeared to be dictated by the presence or absence of an unsaturated bond within the acyl chain. A significantly greater number of transitions were observed during MD simulations in oleic than stearic acid. In addition, significantly greater fluctuation was observed for oleic acid, within the C2 headgroup and C9 and C11 dihedrals (which lie adjacent to the olefin bond of oleic acid). The dynamic behavior of the acyl chain may thereby be more a property of van der Waals contact, and the degree of acyl chain unsaturation, than a function of electrostatics. In the absence of fatty acid, an increase in distance between guanidino carbon centered atoms of Arg126 and Arg106 was observed during MD simulations of the charged apo form. This effect not observed with the neutral apo form or in any of the holo complexes and, presumbably, was a result of repulsion between the negatively charged arginine sidechains. Conserved waters reflected substantially lower mean-square displacement (msd) in all simulations, except the neutral apo form. This suggests that the presence of either charged amino acids or lipid provides increased order for water within the binding pocket. These results provide a dynamic perspective of the interactive nature within the FABP binding pocket regulated in a complex manner by the electrostatics within the binding cavity, acyl chain structure and behavior, and water energetics.