This review article provides an overview of structurally oriented experimental datasets that can be used to benchmark protein force fields, focusing on data generated by nuclear magnetic resonance (NMR) spectroscopy and room temperature (RT) protein crystallography. We discuss what the observables are, what they tell us about structure and dynamics, what makes them useful for assessing force field accuracy, and how they can be connected to molecular dynamics simulations carried out using the force field one wishes to benchmark. We also touch on statistical issues that arise when comparing simulations with experiment. We hope this article will be particularly useful to computational researchers and trainees who develop, benchmark, or use protein force fields for molecular simulations.

Stapled peptides are a promising class of molecules with potential as highly specific probes of protein-protein interactions and as therapeutics. Hydrocarbon stapling affects the peptide properties through the interplay of two factors: enhancing the overall hydrophobicity and constraining the conformational flexibility. By constructing a series of virtual peptides, we study the role of each factor in modulating the structural properties of a hydrocarbon-stapled peptide PM2, which has been shown to enter cells, engage its target Mouse Double Minute 2 (MDM2), and activate p53. Hamiltonian replica exchange molecular dynamics (HREMD) simulations suggest that hydrocarbon stapling favors helical populations of PM2 through a combination of the geometric constraints and the enhanced hydrophobicity of the peptide. To further understand the conformational landscape of the stapled peptides along the binding pathway, we performed HREMD simulations by restraining the peptide at different distances from MDM2. When the peptide approaches MDM2, the binding pocket undergoes dehydration which appears to be greater in the presence of the stapled peptide compared with the linear peptide. In the binding pocket, the helicity of the stapled peptide is increased due to the favorable interactions between the peptide residues as well as the staple and the microenvironment of the binding pocket, contributing to enhanced affinity. The dissection of the multifaceted mechanism of hydrocarbon stapling into individual factors not only deepens fundamental understanding of peptide stapling, but also provides guidelines for the design of new stapled peptides.

Computational modeling (CM) is a versatile scientific methodology used to examine the properties and behavior of complex systems, such as polymeric materials for biomedical bioengineering. CM has emerged as a primary tool for predicting, setting up, and interpreting experimental results. Integrating in silico and in vitro experiments accelerates scientific advancements, yielding quicker results at a reduced cost. While CM is a mature discipline, its use in biomedical engineering for biopolymer materials has only recently gained prominence. In biopolymer biomedical engineering, CM focuses on three key research areas: (A) Computer-aided design (CAD/CAM) utilizes specialized software to design and model biopolymers for various biomedical applications. This technology allows researchers to create precise three-dimensional models of biopolymers, taking into account their chemical, structural, and functional properties. These models can be used to enhance the structure of biopolymers and improve their effectiveness in specific medical applications. (B) Finite element analysis, a computational technique used to analyze and solve problems in engineering and physics. This approach divides the physical domain into small finite elements with simple geometric shapes. This computational technique enables the study and understanding of the mechanical and structural behavior of biopolymers in biomedical environments. (C) Molecular dynamics (MD) simulations involve using advanced computational techniques to study the behavior of biopolymers at the molecular and atomic levels. These simulations are fundamental for better understanding biological processes at the molecular level. Studying the wide-ranging uses of MD simulations in biopolymers involves examining the structural, functional, and evolutionary aspects of biomolecular systems over time. MD simulations solve Newton's equations of motion for all-atom systems, producing spatial trajectories for each atom. This provides valuable insights into properties such as water absorption on biopolymer surfaces and interactions with solid surfaces, which are crucial for assessing biomaterials. This review provides a comprehensive overview of the various applications of MD simulations in biopolymers. Additionally, it highlights the flexibility, robustness, and synergistic relationship between in silico and experimental techniques.

Proteins are molecular machines and to understand how they work, we need to understand how they move. New pump-probe time-resolved X-ray diffraction methods open up ways to initiate and observe protein motions with atomistic detail in crystals on biologically relevant timescales. However, practical limitations of these experiments demands parallel development of effective molecular dynamics approaches to accelerate progress and extract meaning. Here, we establish robust and accurate methods for simulating dynamics in protein crystals, a nontrivial process requiring careful attention to equilibration, environmental composition, and choice of force fields. With more than seven milliseconds of sampling of a single chain, we identify critical factors controlling agreement between simulation and experiments and show that simulated motions recapitulate ligand-induced conformational changes. This work enables a virtuous cycle between simulation and experiments for visualizing and understanding the basic functional motions of proteins.

Molecular dynamics (MD) simulations are routinely performed of biomolecules in solution, because this is their native environment. However, the structures used in such simulations are often obtained with X-ray crystallography, which provides the atomic coordinates of the biomolecule in a crystal environment. With the advent of free electron lasers and time-resolved techniques, X-ray crystallography can now also access metastable states that are intermediates in a biochemical process. Such experiments provide additional data, which can be used, for example, to optimize MD force fields. Doing so requires that the simulation of the biomolecule is also performed in the crystal environment. However, in contrast to simulations of biomolecules in solution, setting up a crystal is challenging. In particular, because not all solvent molecules are resolved in X-ray crystallography, adding a suitable number of solvent molecules, such that the properties of the crystallographic unit cell are preserved in the simulation, can be difficult and typically is a trial-and-error based procedure requiring manual interventions. Such interventions preclude high throughput applications. To overcome this bottleneck, we introduce gmXtal, a tool for setting up crystal simulations for MD simulations with GROMACS. With the information from the protein data bank (rcsb.org) gmXtal automatically (i) builds the crystallographic unit cell; (ii) sets the protonation of titratable residues; (iii) builds missing residues that were not resolved experimentally; and (iv) adds an appropriate number of solvent molecules to the system. gmXtal is available as a standalone tool https://gitlab.com/pbuslaev/gmxtal .

Proteins are molecular machines and to understand how they work, we need to understand how they move. New pump-probe time-resolved X-ray diffraction methods open up ways to initiate and observe protein motions with atomistic detail in crystals on biologically relevant timescales. However, practical limitations of these experiments demands parallel development of effective molecular dynamics approaches to accelerate progress and extract meaning. Here, we establish robust and accurate methods for simulating dynamics in protein crystals, a nontrivial process requiring careful attention to equilibration, environmental composition, and choice of force fields. With more than seven milliseconds of sampling of a single chain, we identify critical factors controlling agreement between simulation and experiments and show that simulated motions recapitulate ligand-induced conformational changes. This work enables a virtuous cycle between simulation and experiments for visualizing and understanding the basic functional motions of proteins.

DNA-protein interactions are pervasive in a number of biophysical processes ranging from transcription and gene expression to chromosome folding. To describe the structural and dynamic properties underlying these processes accurately, it is important to create transferable computational models. Toward this end, we introduce Coarse-grained Force Field for Energy Estimation, COFFEE, a robust framework for simulating DNA-protein complexes. To brew COFFEE, we integrated the energy function in the self-organized polymer model with side-chains for proteins and the three interaction site model for DNA in a modular fashion, without recalibrating any of the parameters in the original force-fields. A unique feature of COFFEE is that it describes sequence-specific DNA-protein interactions using a statistical potential (SP) derived from a data set of high-resolution crystal structures. The only parameter in COFFEE is the strength (λDNAPRO) of the DNA-protein contact potential. For an optimal choice of λDNAPRO, the crystallographic B-factors for DNA-protein complexes with varying sizes and topologies are quantitatively reproduced. Without any further readjustments to the force-field parameters, COFFEE predicts scattering profiles that are in quantitative agreement with small-angle X-ray scattering experiments, as well as chemical shifts that are consistent with NMR. We also show that COFFEE accurately describes the salt-induced unraveling of nucleosomes. Strikingly, our nucleosome simulations explain the destabilization effect of ARG to LYS mutations, which do not alter the balance of electrostatic interactions but affect chemical interactions in subtle ways. The range of applications attests to the transferability of COFFEE, and we anticipate that it would be a promising framework for simulating DNA-protein complexes at the molecular length-scale.

Molecular-dynamics (MD) simulations of protein crystals enable the prediction of structural and dynamical features of both the protein and the solvent components of macromolecular crystals, which can be validated against diffraction data from X-ray crystallographic experiments. The simulations have been useful for studying and predicting both Bragg and diffuse scattering in protein crystallography; however, the preparation is not yet automated and includes choices and tradeoffs that can impact the results. Here we examine some of the intricacies and consequences of the choices involved in setting up MD simulations of protein crystals for the study of diffraction data, and provide a recipe for preparing the simulations, packaged in an accompanying Jupyter notebook. This article and the accompanying notebook are intended to serve as practical resources for researchers wishing to put these models to work.

Some of our most detailed information about structure and dynamics of macromolecules comes from X-ray-diffraction studies in crystalline environments. More than 170,000 atomic models have been deposited in the Protein Data Bank, and the number of observations (typically of intensities of Bragg diffraction peaks) is generally quite large, when compared to other experimental methods. Nevertheless, the general agreement between calculated and observed intensities is far outside the experimental precision, and the majority of scattered photons fall between the sharp Bragg peaks, and are rarely taken into account. This chapter considers how molecular dynamics simulations can be used to explore the connections between microscopic behavior in a crystalline lattice and observed scattering intensities, and point the way to new atomic models that could more faithfully recapitulate Bragg intensities and extract useful information from the diffuse scattering that lies between those peaks.

DNA-protein interactions are pervasive in a number of biophysical processes ranging from transcription, gene expression, to chromosome folding. To describe the structural and dynamic properties underlying these processes accurately, it is important to create transferable computational models. Toward this end, we introduce Coarse grained force field for energy estimation, COFFEE, a robust framework for simulating DNA-protein complexes. To brew COFFEE, we integrated the energy function in the Self-Organized Polymer model with Side Chains for proteins and the Three Interaction Site model for DNA in a modular fashion, without re-calibrating any of the parameters in the original force-fields. A unique feature of COFFEE is that it describes sequence-specific DNA-protein interactions using a statistical potential (SP) derived from a dataset of high-resolution crystal structures. The only parameter in COFFEE is the strength (λ D N A P R O ) of the DNA-protein contact potential. For an optimal choice of λ D N A P R O , the crystallographic B-factors for DNA-protein complexes, with varying sizes and topologies, are quantitatively reproduced. Without any further readjustments to the force-field parameters, COFFEE predicts the scattering profiles that are in quantitative agreement with SAXS experiments as well as chemical shifts that are consistent with NMR. We also show that COFFEE accurately describes the salt-induced unraveling of nucleosomes. Strikingly, our nucleosome simulations explain the destabilization effect of ARG to LYS mutations, which does not alter the balance of electrostatic interactions, but affects chemical interactions in subtle ways. The range of applications attests to the transferability of COFFEE, and we anticipate that it would be a promising framework for simulating DNA-protein complexes at the molecular length-scale.

Accurate modeling of protein-water interactions in molecular dynamics (MD) simulations is important for understanding the molecular basis of protein function. Data from x-ray crystallography can be useful in assessing the accuracy of MD simulations, in particular, the locations of crystallographic water sites (CWS) coordinated by the protein. Such a comparison requires special methodological considerations that take into account the dynamic nature of proteins. However, existing methods for analyzing CWS in MD simulations rely on global alignment of the protein onto the crystal structure, which introduces substantial errors in the case of significant structural deviations. Here, we propose a method called local alignment for water sites (LAWS), which is based on multilateration-an algorithm widely used in GPS tracking. LAWS considers the contacts formed by CWS and protein atoms in the crystal structure and uses these interaction distances to track CWS in a simulation. We apply our method to simulations of a protein crystal and to simulations of the same protein in solution. Compared with existing methods, LAWS defines CWS characterized by more prominent water density peaks and a less-perturbed protein environment. In the crystal, we find that all high-confidence crystallographic waters are preserved. Using LAWS, we demonstrate the importance of crystal packing for the stability of CWS in the unit cell. Simulations of the protein in solution and in the crystal share a common set of preserved CWS that are located in pockets and coordinated by residues of the same domain, which suggests that the LAWS algorithm will also be useful in studying ordered waters and water networks in general.

Enzymes catalyze biochemical reactions through precise positioning of substrates, cofactors, and amino acids to modulate the transition-state free energy. However, the role of conformational dynamics remains poorly understood due to lack of experimental access. This shortcoming is evident with E. coli dihydrofolate reductase (DHFR), a model system for the role of protein dynamics in catalysis, for which it is unknown how the enzyme regulates the different active site environments required to facilitate proton and hydride transfer. Here, we present ligand-, temperature-, and electric-field-based perturbations during X-ray diffraction experiments that enable identification of coupled conformational changes in DHFR. We identify a global hinge motion and local networks of structural rearrangements that are engaged by substrate protonation to regulate solvent access and promote efficient catalysis. The resulting mechanism shows that DHFR's two-step catalytic mechanism is guided by a dynamic free energy landscape responsive to the state of the substrate.

There still is little treatment available for amyloid diseases, despite their significant impact on individuals and the social and economic implications for society. One reason for this is that the physical nature of amyloid formation is not understood sufficiently well. Therefore, fundamental research at the molecular level remains necessary to support the development of therapeutics. A few structures of short peptides from amyloid-forming proteins have been determined. These can in principle be used as scaffolds for designing aggregation inhibitors. Attempts to this end have often used the tools of computational chemistry, in particular molecular simulation. However, few simulation studies of these peptides in the crystal state have been presented so far. Hence, to validate the capability of common force fields (AMBER19SB, CHARMM36m, and OPLS-AA/M) to yield insight into the dynamics and structural stability of amyloid peptide aggregates, we have performed molecular dynamics simulations of twelve different peptide crystals at two different temperatures. From the simulations, we evaluate the hydrogen bonding patterns, the isotropic B-factors, the change in energy, the Ramachandran plots, and the unit cell parameters and compare the results with the crystal structures. Most crystals are stable in the simulations but for all force fields there is at least one that deviates from the experimental crystal, suggesting more work is needed on these models.

A molecular dynamics (MD)-based pipeline has been designed and implemented to emulate the entire process of collecting diffraction photographs and calculating crystallographic structures of proteins from them. Using a structure of lysozyme solved in-house, a supercell comprising 125 (5 × 5 × 5) crystal unit cells containing a total of 1000 protein molecules and explicit interstitial solvent was constructed. For this system, two 300 ns MD trajectories at 298 and 250 K were recorded. A series of snapshots from these trajectories were then used to simulate a fully realistic set of diffraction photographs, which were further fed into the standard pipeline for structure determination. The resulting structures show very good agreement with the underlying MD model not only in terms of coordinates but also in terms of B factors; they are also consistent with the original experimental structure. The developed methodology should find a range of applications, such as optimizing refinement protocols to solve crystal structures and extracting dynamics information from diffraction data or diffuse scattering.

Probing the dynamics of aromatic side chains provides important insights into the behavior of a protein because flips of aromatic rings in a protein's hydrophobic core report on breathing motion involving a large part of the protein. Inherently invisible to crystallography, aromatic motions have been primarily studied by solution NMR. The question how packing of proteins in crystals affects ring flips has, thus, remained largely unexplored. Here we apply magic-angle spinning NMR, advanced phenylalanine 1H-13C/2H isotope labeling and MD simulation to a protein in three different crystal packing environments to shed light onto possible impact of packing on ring flips. The flips of the two Phe residues in ubiquitin, both surface exposed, appear remarkably conserved in the different crystal forms, even though the intermolecular packing is quite different: Phe4 flips on a ca. 10-20 ns time scale, and Phe45 are broadened in all crystals, presumably due to µs motion. Our findings suggest that intramolecular influences are more important for ring flips than intermolecular (packing) effects.

Correlated motions in proteins arising from the collective movements of residues have long been proposed to be fundamentally important to key properties of proteins, from allostery and catalysis to evolvability. Recent breakthroughs in structural biology have made it possible to capture proteins undergoing complex conformational changes, yet intrinsic correlated motions within a conformation remain one of the least understood facets of protein structure. For many decades, the analysis of total X-ray scattering held the promise of animating crystal structures with correlated motions. With recent advances in both X-ray detectors and data interpretation methods, this long-held promise can now be met. In this Perspective, we will introduce how correlated motions are captured in total scattering and provide guidelines for the collection, interpretation, and validation of data. As structural biology continues to push the boundaries, we see an opportunity to gain atomistic insight into correlated motions using total scattering as a bridge between theory and experiment.

The determination of G protein-coupled receptor (GPCR) structures at atomic resolution has improved understanding of cellular signaling and will accelerate the development of new drug candidates. However, experimental structures still remain unavailable for a majority of the GPCR family. GPCR structures and their interactions with ligands can also be modelled computationally, but such predictions have limited accuracy. In this work, we explored if molecular dynamics (MD) simulations could be used to refine the accuracy of in silico models of receptor-ligand complexes that were submitted to a community-wide assessment of GPCR structure prediction (GPCR Dock). Two simulation protocols were used to refine 30 models of the D3 dopamine receptor (D3R) in complex with an antagonist. Close to 60 μs of simulation time was generated and the resulting MD refined models were compared to a D3R crystal structure. In the MD simulations, the receptor models generally drifted further away from the crystal structure conformation. However, MD refinement was able to improve the accuracy of the ligand binding mode. The best refinement protocol improved agreement with the experimentally observed ligand binding mode for a majority of the models. Receptor structures with improved virtual screening performance, which was assessed by molecular docking of ligands and decoys, could also be identified among the MD refined models. Application of weak restraints to the transmembrane helixes in the MD simulations further improved predictions of the ligand binding mode and second extracellular loop. These results provide guidelines for application of MD refinement in prediction of GPCR-ligand complexes and directions for further method development.

Conventional protein:ligand crystallographic refinement uses stereochemistry restraints coupled with a rudimentary energy functional to ensure the correct geometry of the model of the macromolecule-along with any bound ligand(s)-within the context of the experimental, X-ray density. These methods generally lack explicit terms for electrostatics, polarization, dispersion, hydrogen bonds, and other key interactions, and instead they use pre-determined parameters (e.g. bond lengths, angles, and torsions) to drive structural refinement. In order to address this deficiency and obtain a more complete and ultimately more accurate structure, we have developed an automated approach for macromolecular refinement based on a two layer, QM/MM (ONIOM) scheme as implemented within our DivCon Discovery Suite and "plugged in" to two mainstream crystallographic packages: PHENIX and BUSTER. This implementation is able to use one or more region layer(s), which is(are) characterized using linear-scaling, semi-empirical quantum mechanics, followed by a system layer which includes the balance of the model and which is described using a molecular mechanics functional. In this work, we applied our Phenix/DivCon refinement method-coupled with our XModeScore method for experimental tautomer/protomer state determination-to the characterization of structure sets relevant to structure-based drug design (SBDD). We then use these newly refined structures to show the impact of QM/MM X-ray refined structure on our understanding of function by exploring the influence of these improved structures on protein:ligand binding affinity prediction (and we likewise show how we use post-refinement scoring outliers to inform subsequent X-ray crystallographic efforts). Through this endeavor, we demonstrate a computational chemistry ↔ structural biology (X-ray crystallography) "feedback loop" which has utility in industrial and academic pharmaceutical research as well as other allied fields.

Recent technical advances have dramatically increased the power and scope of structural biology. New developments in high-resolution cryo-electron microscopy, serial X-ray crystallography, and electron diffraction have been especially transformative. Here we highlight some of the latest advances and current challenges at the frontiers of atomic resolution methods for elucidating the structures and dynamical properties of macromolecules and their complexes.

Crystallographic data comes from a space-time average over all the unit cells within the crystal, so dynamic phenomena do not contribute significantly to the diffraction data. Many efforts have been made to reconstitute the movement of the macromolecules and explore the microstates that the confined proteins can adopt in the crystalline network. We explored different strategies to simulate a heart fatty acid binding protein (H-FABP) crystal by means of Molecular Dynamics (MD) simulations. We evaluate the effect of introducing restraints according to experimental isotropic B-factors and we analyzed the H-FABP motions in the crystal using Principal Component Analysis (PCA), isotropic and anisotropic B-factors. We compared the behavior of the protein simulated in the crystal confinement versus in solution, and we observed the effect of that confinement in the mobility of the protein residues. Restraining one-third of Cα atoms based on experimental B-factors produce lower B-factors than simulations without restraints, showing that the position restraint of the atoms with the lowest experimental B-factor is a good strategy to maintain the geometry of the crystal with an obvious decrease in the degrees of motion of the protein. PCA shows that, as position restraint reduces the conformational space explored by the system, the motion of the crystal is better recovered, for an essential subspace of the same size, in the simulations without restraints. Restraining only one Cα seems to be a good balance between giving flexibility to the system and preserving its structure. Communicated by Ramaswamy H. Sarma.

Crystal simulations provide useful tools, along with solution simulations, to test nucleic acid force fields, but should be interpreted with care owing to the difficulty of establishing the environmental conditions needed to reproduce experimental crystal packing. These challenges underscore the need to construct proper protocols for carrying out crystal simulations and analyzing results to identify the origin of deviations from crystallographic data. Toward this end, we introduce a novel framework for B-factor decomposition into additive intramolecular, rotational, and translational atomic fluctuation components and partitioning of each of these components into individual asymmetric unit and lattice contributions. We apply the framework to a benchmark set of A-DNA, Z-DNA, and B-DNA double helix systems of various chain lengths. Overall, the intramolecular deviations from the crystal were quite small (≤1.0 Å), suggesting high accuracy of the force field, whereas crystal packing was not well reproduced. The present work establishes a framework to conduct and analyze crystal simulations that ultimately take on issues of crystal packing and can provide insight into nucleic acid force fields.

The structures of biological macromolecules would not be known to their present extent without X-ray crystallography. Most simulations of globular proteins in solution begin by surrounding the crystal structure of the monomer in a bath of water molecules, but the standard simulation employing periodic boundary conditions is already close to a crystal lattice environment. With simple protocols, the same software and molecular models can perform simulations of the crystal lattice, including all asymmetric units and solvent to fill the box. Throughout the history of molecular dynamics, studies of crystal lattices have served to investigate the quality of the underlying force fields, correlate the simulated ensembles to experimental structure factors, and extrapolate the behavior in lattices to behavior in solution. Powerful new computers are enabling molecular simulations with greater realism and statistical convergence. Meanwhile, the advent of exciting new methods in crystallography, including femtosecond free-electron lasers and image reconstruction for time-resolved crystallography on slurries of small crystals, is expanding the range of structures accessible to X-ray diffraction. We review past fusions of simulations and crystallography, then look ahead to the ways that simulations of crystal structures will enhance structural biology in the future.

What factors favor protein folding? This is a textbook question. Parsing the experimental free energies of folding/unfolding into diverse enthalpic and entropic components of solute and solvent favoring or disfavoring folding is not an easy task. In this study, we present a computational protocol for estimating the free energy contributors to protein folding semi-quantitatively using ensembles of unfolded and native states generated via molecular dynamics simulations. We tested the methodology on 35 proteins with diverse structural motifs and sizes and found that the calculated free energies correlate well with experiment (correlation coefficient ∼ 0.85), enabling us to develop a consensus view of the energetics of folding. As a more sensitive test of the methodology, we also investigated the free energies of folding of an additional 33 single point mutants and obtained a correlation coefficient of 0.8. A notable observation is that the folding free energy components appear to carry signatures of the fold (SCOP classification) of the protein. © 2017 Wiley Periodicals, Inc.

An accurate potential energy model is crucial for biomolecular simulations. Despite many recent improvements of classical protein force fields, there are remaining key issues: much weaker temperature dependence of folding/unfolding equilibrium and overly collapsed unfolded or disordered states. For the latter problem, a new water model (TIP4P-D) has been proposed to correct the significantly underestimated water dispersion interactions. Here, using TIP4P-D, we reveal problems in current force fields through failures in folding model systems (a polyalanine peptide, Trp-cage, and the GB1 hairpin). By using residue-specific parameters to achieve better match between amino acid sequences and native structures and adding a small H-bond correction to partially compensate the missing many-body effects in α-helix formation, the new RSFF2+ force field with the TIP4P-D water model can excellently reproduce experimental melting curves of both α-helical and β-hairpin systems. The RSFF2+/TIP4P-D method also gives less collapsed unfolded structures and describes well folded proteins simultaneously.

X-ray scattering is uniquely suited to the study of disordered systems and thus has the potential to provide insight into dynamic processes where diffraction methods fail. In particular, while X-ray crystallography has been a staple of structural biology for more than half a century and will continue to remain so, a major limitation of this technique has been the lack of dynamic information. Solution X-ray scattering has become an invaluable tool in structural and mechanistic studies of biological macromolecules where large conformational changes are involved. Such systems include allosteric enzymes that play key roles in directing metabolic fluxes of biochemical pathways, as well as large, assembly-line type enzymes that synthesize secondary metabolites with pharmaceutical applications. Furthermore, crystallography has the potential to provide information on protein dynamics via the diffuse scattering patterns that are overlaid with Bragg diffraction. Historically, these patterns have been very difficult to interpret, but recent advances in X-ray detection have led to a renewed interest in diffuse scattering analysis as a way to probe correlated motions. Here, we will review X-ray scattering theory and highlight recent advances in scattering-based investigations of protein solutions and crystals, with a particular focus on complex enzymes.

Protein structures are essential in modern biology yet experimental methods are far from being able to catch up with the rapid increase in available genomic data. Computational protein structure prediction methods aim to fill the gap while the role of protein structure refinement is to take approximate initial template-based models and bring them closer to the true native structure. Current methods for computational structure refinement rely on molecular dynamics simulations, related sampling methods, or iterative structure optimization protocols. The best methods are able to achieve moderate degrees of refinement but consistent refinement that can reach near-experimental accuracy remains elusive. Key issues revolve around the accuracy of the energy function, the inability to reliably rank multiple models, and the use of restraints that keep sampling close to the native state but also limit the degree of possible refinement. A different aspect is the question of what exactly the target of high-resolution refinement should be as experimental structures are affected by experimental conditions and different biological questions require varying levels of accuracy. While improvement of the global protein structure is a difficult problem, high-resolution refinement methods that improves local structural quality such as favorable stereochemistry and the avoidance of atomic clashes are much more successful.

Biomolecular X-ray structures typically provide a static, time- and ensemble-averaged view of molecular ensembles in crystals. In the absence of rigid-body motions and lattice defects, B-factors are thought to accurately reflect the structural heterogeneity of such ensembles. In order to study the effects of averaging on B-factors, we employ molecular dynamics simulations to controllably manipulate microscopic heterogeneity of a crystal containing 216 copies of villin headpiece. Using average structure factors derived from simulation, we analyse how well this heterogeneity is captured by high-resolution molecular-replacement-based model refinement. We find that both isotropic and anisotropic refined B-factors often significantly deviate from their actual values known from simulation: even at high 1.0 Å resolution and Rfree of 5.9%, B-factors of some well-resolved atoms underestimate their actual values even sixfold. Our results suggest that conformational averaging and inadequate treatment of correlated motion considerably influence estimation of microscopic heterogeneity via B-factors, and invite caution in their interpretation.

During protein folding and as part of some conformational changes that regulate protein function, the polypeptide chain must traverse high-energy barriers that separate the commonly adopted low-energy conformations. How distortions in peptide geometry allow these barrier-crossing transitions is a fundamental open question. One such important transition involves the movement of a non-glycine residue between the left side of the Ramachandran plot (that is, ϕ < 0°) and the right side (that is, ϕ > 0°). We report that high-energy conformations with ϕ ~ 0°, normally expected to occur only as fleeting transition states, are stably trapped in certain highly resolved native protein structures and that an analysis of these residues provides a detailed, experimentally derived map of the bond angle distortions taking place along the transition path. This unanticipated information lays to rest any uncertainty about whether such transitions are possible and how they occur, and in doing so lays a firm foundation for theoretical studies to better understand the transitions between basins that have been little studied but are integrally involved in protein folding and function. Also, the context of one such residue shows that even a designed highly stable protein can harbor substantial unfavorable interactions.

Molecular mechanics is powerful for its speed in atomistic simulations, but an accurate force field is required. The Amber ff99SB force field improved protein secondary structure balance and dynamics from earlier force fields like ff99, but weaknesses in side chain rotamer and backbone secondary structure preferences have been identified. Here, we performed a complete refit of all amino acid side chain dihedral parameters, which had been carried over from ff94. The training set of conformations included multidimensional dihedral scans designed to improve transferability of the parameters. Improvement in all amino acids was obtained as compared to ff99SB. Parameters were also generated for alternate protonation states of ionizable side chains. Average errors in relative energies of pairs of conformations were under 1.0 kcal/mol as compared to QM, reduced 35% from ff99SB. We also took the opportunity to make empirical adjustments to the protein backbone dihedral parameters as compared to ff99SB. Multiple small adjustments of φ and ψ parameters were tested against NMR scalar coupling data and secondary structure content for short peptides. The best results were obtained from a physically motivated adjustment to the φ rotational profile that compensates for lack of ff99SB QM training data in the β-ppII transition region. Together, these backbone and side chain modifications (hereafter called ff14SB) not only better reproduced their benchmarks, but also improved secondary structure content in small peptides and reproduction of NMR χ1 scalar coupling measurements for proteins in solution. We also discuss the Amber ff12SB parameter set, a preliminary version of ff14SB that includes most of its improvements.

Molecular dynamics simulations of crystals can enlighten interpretation of experimental X-ray crystallography data and elucidate structural dynamics and heterogeneity in biomolecular crystals. Furthermore, because of the direct comparison against experimental data, they can inform assessment of molecular dynamics methods and force fields. We present microsecond scale results for triclinic hen egg-white lysozyme in a supercell consisting of 12 independent unit cells using four contemporary force fields (Amber ff99SB, ff14ipq, ff14SB, and CHARMM 36) in crystalline and solvated states (for ff14SB only). We find the crystal simulations consistent across multiple runs of the same force field and robust to various solvent equilibration schemes. However, convergence is slow compared with solvent simulations. All the tested force fields reproduce experimental structural and dynamic properties well, but Amber ff14SB maintains structure and reproduces fluctuations closest to the experimental model: its average backbone structure differs from the deposited structure by 0.37Å; by contrast, the average backbone structure in solution differs from the deposited by 0.65Å. All the simulations are affected by a small progressive deterioration of the crystal lattice, presumably due to imperfect modeling of hydrogen bonding and other crystal contact interactions; this artifact is smallest in ff14SB, with average lattice positions deviating by 0.20Å from ideal. Side-chain disorder is surprisingly low with fewer than 30% of the nonglycine or alanine residues exhibiting significantly populated alternate rotamers. Our results provide helpful insight into the methodology of biomolecular crystal simulations and indicate directions for future work to obtain more accurate energy models for molecular dynamics.

Currently, the best existing molecular dynamics (MD) force fields cannot accurately reproduce the global free-energy minimum which realizes the experimental protein structure. As a result, long MD trajectories tend to drift away from the starting coordinates (e.g., crystallographic structures). To address this problem, we have devised a new simulation strategy aimed at protein crystals. An MD simulation of protein crystal is essentially an ensemble simulation involving multiple protein molecules in a crystal unit cell (or a block of unit cells). To ensure that average protein coordinates remain correct during the simulation, we introduced crystallography-based restraints into the MD protocol. Because these restraints are aimed at the ensemble-average structure, they have only minimal impact on conformational dynamics of the individual protein molecules. So long as the average structure remains reasonable, the proteins move in a native-like fashion as dictated by the original force field. To validate this approach, we have used the data from solid-state NMR spectroscopy, which is the orthogonal experimental technique uniquely sensitive to protein local dynamics. The new method has been tested on the well-established model protein, ubiquitin. The ensemble-restrained MD simulations produced lower crystallographic R factors than conventional simulations; they also led to more accurate predictions for crystallographic temperature factors, solid-state chemical shifts, and backbone order parameters. The predictions for (15) N R1 relaxation rates are at least as accurate as those obtained from conventional simulations. Taken together, these results suggest that the presented trajectories may be among the most realistic protein MD simulations ever reported. In this context, the ensemble restraints based on high-resolution crystallographic data can be viewed as protein-specific empirical corrections to the standard force fields.

We report over 30 μs of unrestrained molecular dynamics simulations of six protein-RNA complexes in explicit solvent. We utilize the AMBER ff99bsc0χ(OL3) RNA force field combined with the ff99SB protein force field and its more recent ff12SB version with reparametrized side-chain dihedrals. The simulations show variable behavior, ranging from systems that are essentially stable to systems with progressive deviations from the experimental structure, which we could not stabilize anywhere close to the starting experimental structure. For some systems, microsecond-scale simulations are necessary to achieve stabilization after initial sizable structural perturbations. The results show that simulations of protein-RNA complexes are challenging and every system should be treated individually. The simulations are affected by numerous factors, including properties of the starting structures (the initially high force field potential energy, resolution limits, conformational averaging, crystal packing, etc.), force field imbalances, and real flexibility of the studied systems. These factors, and thus the simulation behavior, differ from system to system. The structural stability of simulated systems does not correlate with the size of buried interaction surface or experimentally determined binding affinities but reflects the type of protein-RNA recognition. Protein-RNA interfaces involving shape-specific recognition of RNA are more stable than those relying on sequence-specific RNA recognition. The differences between the protein force fields are considerably smaller than the uncertainties caused by sampling and starting structures. The ff12SB improves description of the tyrosine side-chain group, which eliminates some problems associated with tyrosine dynamics.

We introduce a labeling scheme for magic angle spinning (MAS) solid-state NMR that is based on deuteration in combination with dilution of the carbon spin system. The labeling strategy achieves spectral editing by simplification of the HαCα and aliphatic side chain spectral region. A reduction in both proton and carbon spin density in combination with fast spinning (≥50 kHz) is essential to retrieve artifact-free (13)C-R1 relaxation data for aliphatic carbons. We obtain good agreement between the NMR experimental data and order parameters extracted from a molecular dynamics (MD) trajectory, which indicates that carbon based relaxation parameters can yield complementary information on protein backbone as well as side chain dynamics.

Two 250 ns molecular simulations have been carried out to study the structure and dynamics of crystal toxin protein II from the scorpion Androctonus australis Hector employing the polarized protein-specific charge (PPC), as well as the standard AMBER99SB force field, to investigate the electrostatic polarization on the simulated crystal stability. Results show that under PPC, the monomers in unit cell as well as the lattice in supercell are more stable with smaller root-mean-square deviations and more accurate lattice atomic fluctuations compared with the crystallographic B-factors than under AMBER99SB force field. Most of the interactions at interfaces in the X-ray structure are quite well-preserved, underscoring the important effect of polarization on maintaining the crystal stability. However, the results also show that the hydrogen bond between Asp53 and Gln37 and the cation-π interaction between Arg56 and His64 are not stable, indicating that further optimization of force field, especially the van der Waals interaction parameters, is desired.

Intrinsically disordered proteins (IDPs) often rely on electrostatic interactions to bind their structured targets. To obtain insight into the mechanism of formation of the electrostatic encounter complex, we investigated the binding of the peptide Sos (PPPVPPRRRR), which serves as a minimal model for an IDP, to the c-Crk N-terminal SH3 domain. Initially, we measured ¹⁵N relaxation rates at two magnetic field strengths and determined the binding shifts for the complex of Sos with wild-type SH3. We have also recorded a 3 μs molecular dynamics (MD) trajectory of this complex using the Amber ff99SB*-ILDN force field. The comparison of the experimental and simulated data shows that MD simulation consistently overestimates the strength of salt bridge interactions at the binding interface. The series of simulations using other advanced force fields also failed to produce any satisfactory results. To address this issue, we have devised an empirical correction to the Amber ff99SB*-ILDN force field whereby the Lennard-Jones equilibrium distance for the nitrogen-oxygen pair across the Arg-to-Asp and Arg-to-Glu salt bridges has been increased by 3%. Implementing this correction resulted in a good agreement between the simulations and the experiment. Adjusting the strength of salt bridge interactions removed a certain amount of strain contained in the original MD model, thus improving the binding of the hydrophobic N-terminal portion of the peptide. The arginine-rich C-terminal portion of the peptide, freed from the effect of the overstabilized salt bridges, was found to interconvert more rapidly between its multiple conformational states. The modified MD protocol has also been successfully used to simulate the entire binding process. In doing so, the peptide was initially placed high above the protein surface. It then arrived at the correct bound pose within ∼2 Å of the crystallographic coordinates. This simulation allowed us to analyze the details of the dynamic binding intermediate, i.e., the electrostatic encounter complex. However, an experimental characterization of this transient, weakly populated state remains out of reach. To overcome this problem, we designed the double mutant of c-Crk N-SH3 in which mutations Y186L and W169F abrogate tight Sos binding and shift the equilibrium toward the intermediate state resembling the electrostatic encounter complex. The results of the combined NMR and MD study of this engineered system will be reported in the next part of this paper.

Molecular dynamics simulations can complement experimental measures of structure and dynamics of biomolecules. The quality of such simulations can be tested by comparisons to models refined against experimental crystallographic data.

We report simulations of DNA and RNA duplexes in their crystalline environment. The calculations mimic the conditions for PDB entries 1D23 [d(CGATCGATCG)2] and 1RNA [(UUAUAUAUAUAUAA)2], and contain 8 unit cells, each with 4 copies of the Watson-Crick duplex; this yields in aggregate 64μs of duplex sampling for DNA and 16μs for RNA.

The duplex structures conform much more closely to the average structure seen in the crystal than do structures extracted from a solution simulation with the same force field. Sequence-dependent variations in helical parameters, and in groove widths, are largely maintained in the crystal structure, but are smoothed out in solution. However, the integrity of the crystal lattice is slowly degraded in both simulations, with the result that the interfaces between chains become heterogeneous. This problem is more severe for the DNA crystal, which has fewer inter-chain hydrogen bond contacts than does the RNA crystal.

Crystal simulations using current force fields reproduce many features of observed crystal structures, but suffer from a gradual degradation of the integrity of the crystal lattice.

The results offer insights into force-field simulations that test their ability to preserve weak interactions between chains, which will be of importance also in non-crystalline applications that involve binding and recognition. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

An overview of the current capabilities and limitations of molecular simulation of biomolecular complexes in the context of computer-aided drug design is provided. Steady improvements in computer hardware coupled with more refined representations of energetics are leading to a new appreciation of the driving forces of molecular recognition. Molecular simulations are poised to more frequently guide the interpretation of biophysical measurements of biomolecular complexes. Ligand design strategies emerge from detailed analyses of computed structural ensembles. The feasibility of routine applications to ligand optimization problems hinges upon successful extensive large scale validation studies and the development of protocols to intelligently automate computations.