Elastic network models (ENMs) and constraint-based, topological rigidity analysis are two distinct, coarse-grained approaches to study conformational flexibility of macromolecules. In the two decades since their introduction, both have contributed significantly to insights into protein molecular mechanisms and function. However, despite a shared purpose of these approaches, the topological nature of rigidity analysis, and thereby the absence of motion modes, has impeded a direct comparison. Here, we present an alternative, kinematic approach to rigidity analysis, which circumvents these drawbacks. We introduce a novel protein hydrogen bond network spectral decomposition, which provides an orthonormal basis for collective motions modulated by noncovalent interactions, analogous to the eigenspectrum of normal modes. The zero modes decompose proteins into rigid clusters identical to those from topological rigidity, while nonzero modes rank protein motions by their hydrogen bond collective energy penalty. Our kinematic flexibility analysis bridges topological rigidity theory and ENM, enabling a detailed analysis of motion modes obtained from both approaches. Analysis of a large, structurally diverse data set revealed that collectivity of protein motions, reported by the Shannon entropy, is significantly reduced for rigidity theory compared to normal mode approaches. Strikingly, kinematic flexibility analysis suggests that the hydrogen bonding network encodes a protein-fold specific, spatial hierarchy of motions, which goes nearly undetected in ENM. This hierarchy reveals distinct motion regimes that rationalize experimental and simulated protein stiffness variations. Kinematic motion modes highly correlate with reported crystallographic B factors and molecular dynamics simulations of adenylate kinase. A formal expression for changes in free energy derived from the spectral decomposition indicates that motions across nearly 40% of modes obey enthalpy-entropy compensation. Taken together, our results suggest that hydrogen bond networks have evolved to modulate protein structure and dynamics, which can be efficiently probed by kinematic flexibility analysis.

Protein crystallography data collection at synchrotrons is routinely carried out at cryogenic temperatures to mitigate radiation damage. Although damage still takes place at 100 K and below, the immobilization of free radicals increases the lifetime of the crystals by approximately 100-fold. Recent studies have shown that flash-cooling decreases the heterogeneity of the conformational ensemble and can hide important functional mechanisms from observation. These discoveries have motivated increasing numbers of experiments to be carried out at room temperature. However, the trade-offs between increased risk of radiation damage and increased observation of alternative conformations at room temperature relative to cryogenic temperature have not been examined. A considerable amount of effort has previously been spent studying radiation damage at cryo-temperatures, but the relevance of these studies to room temperature diffraction is not well understood. Here, the effects of radiation damage on the conformational landscapes of three different proteins (T. danielli thaumatin, hen egg-white lysozyme and human cyclophilin A) at room (278 K) and cryogenic (100 K) temperatures are investigated. Increasingly damaged datasets were collected at each temperature, up to a maximum dose of the order of 10⁷ Gy at 100 K and 10⁵ Gy at 278 K. Although it was not possible to discern a clear trend between damage and multiple conformations at either temperature, it was observed that disorder, monitored by B-factor-dependent crystallographic order parameters, increased with higher absorbed dose for the three proteins at 100 K. At 278 K, however, the total increase in this disorder was only statistically significant for thaumatin. A correlation between specific radiation damage affecting side chains and the amount of disorder was not observed. This analysis suggests that elevated conformational heterogeneity in crystal structures at room temperature is observed despite radiation damage, and not as a result thereof.

NMR relaxation dispersion measurements report on conformational changes occurring on the μs-ms timescale. Chemical shift information derived from relaxation dispersion can be used to generate structural models of weakly populated alternative conformational states. Current methods to obtain such models rely on determining the signs of chemical shift changes between the conformational states, which are difficult to obtain in many situations. Here, we use a "sample and select" method to generate relevant structural models of alternative conformations of the C-terminal-associated region of Escherichia coli dihydrofolate reductase (DHFR), using only unsigned chemical shift changes for backbone amides and carbonyls (¹H, ¹⁵N, and ¹³C'). We find that CS-Rosetta sampling with unsigned chemical shift changes generates a diversity of structures that are sufficient to characterize a minor conformational state of the C-terminal region of DHFR. The excited state differs from the ground state by a change in secondary structure, consistent with previous predictions from chemical shift hypersurfaces and validated by the x-ray structure of a partially humanized mutant of E. coli DHFR (N23PP/G51PEKN). The results demonstrate that the combination of fragment modeling with sparse chemical shift data can determine the structure of an alternative conformation of DHFR sampled on the μs-ms timescale. Such methods will be useful for characterizing alternative states, which can potentially be used for in silico drug screening, as well as contributing to understanding the role of minor states in biology and molecular evolution.

Detailed descriptions of atomic coordinates and motions are required for an understanding of protein dynamics and their relation to molecular recognition, catalytic function, and allostery. Historically, NMR relaxation measurements have played a dominant role in the determination of the amplitudes and timescales (picosecond-nanosecond) of bond vector fluctuations, whereas high-resolution X-ray diffraction experiments can reveal the presence of and provide atomic coordinates for multiple, weakly populated substates in the protein conformational ensemble. Here we report a hybrid NMR and X-ray crystallography analysis that provides a more complete dynamic picture and a more quantitative description of the timescale and amplitude of fluctuations in atomic coordinates than is obtainable from the individual methods alone. Order parameters (S(2)) were calculated from single-conformer and multiconformer models fitted to room temperature and cryogenic X-ray diffraction data for dihydrofolate reductase. Backbone and side-chain order parameters derived from NMR relaxation experiments are in excellent agreement with those calculated from the room-temperature single-conformer and multiconformer models, showing that the picosecond timescale motions observed in solution occur also in the crystalline state. These motions are quenched in the crystal at cryogenic temperatures. The combination of NMR and X-ray crystallography in iterative refinement promises to provide an atomic resolution description of the alternate conformational substates that are sampled through picosecond to nanosecond timescale fluctuations of the protein structure. The method also provides insights into the structural heterogeneity of nonmethyl side chains, aromatic residues, and ligands, which are less commonly analyzed by NMR relaxation measurements.

While protein conformational heterogeneity plays an important role in many aspects of biological function, including ligand binding, its impact has been difficult to quantify. Macromolecular X-ray diffraction is commonly interpreted with a static structure, but it can provide information on both the anharmonic and harmonic contributions to conformational heterogeneity. Here, through multiconformer modeling of time- and space-averaged electron density, we measure conformational heterogeneity of 743 stringently matched pairs of crystallographic datasets that reflect unbound/apo and ligand-bound/holo states. When comparing the conformational heterogeneity of side chains, we observe that when binding site residues become more rigid upon ligand binding, distant residues tend to become more flexible, especially in non-solvent-exposed regions. Among ligand properties, we observe increased protein flexibility as the number of hydrogen bonds decreases and relative hydrophobicity increases. Across a series of 13 inhibitor-bound structures of CDK2, we find that conformational heterogeneity is correlated with inhibitor features and identify how conformational changes propagate differences in conformational heterogeneity away from the binding site. Collectively, our findings agree with models emerging from nuclear magnetic resonance studies suggesting that residual side-chain entropy can modulate affinity and point to the need to integrate both static conformational changes and conformational heterogeneity in models of ligand binding.

Virtual screening is a widely used tool for drug discovery, but its predictive power can vary dramatically depending on how much structural data is available. In the best case, crystal structures of a ligand-bound protein can help find more potent ligands. However, virtual screens tend to be less predictive when only ligand-free crystal structures are available, and even less predictive if a homology model or other predicted structure must be used. Here, we explore the possibility that this situation can be improved by better accounting for protein dynamics, as simulations started from a single structure have a reasonable chance of sampling nearby structures that are more compatible with ligand binding. As a specific example, we consider the cancer drug target PPM1D/Wip1 phosphatase, a protein that lacks crystal structures. High-throughput screens have led to the discovery of several allosteric inhibitors of PPM1D, but their binding mode remains unknown. To enable further drug discovery efforts, we assessed the predictive power of an AlphaFold-predicted structure of PPM1D and a Markov state model (MSM) built from molecular dynamics simulations initiated from that structure. Our simulations reveal a cryptic pocket at the interface between two important structural elements, the flap and hinge regions. Using deep learning to predict the pose quality of each docked compound for the active site and cryptic pocket suggests that the inhibitors strongly prefer binding to the cryptic pocket, consistent with their allosteric effect. The predicted affinities for the dynamically uncovered cryptic pocket also recapitulate the relative potencies of the compounds (τ_b = 0.70) better than the predicted affinities for the static AlphaFold-predicted structure (τ_b = 0.42). Taken together, these results suggest that targeting the cryptic pocket is a good strategy for drugging PPM1D and, more generally, that conformations selected from simulation can improve virtual screening when limited structural data is available.

Structure-based, virtual High-Throughput Screening (vHTS) methods for predicting ligand activity in drug discovery are important when there are no or relatively few known compounds that interact with a therapeutic target of interest. State-of-the-art computational vHTS necessarily relies on effective methods for pose sampling and docking and generating an accurate affinity score from the docked poses. However, proteins are dynamic; in vivo ligands bind to a conformational ensemble. In silico docking to the single conformation represented by a crystal structure can adversely affect the pose quality. Here, we introduce AtomNet PoseRanker (ANPR), a graph convolutional network trained to identify and rerank crystal-like ligand poses from a sampled ensemble of protein conformations and ligand poses. In contrast to conventional vHTS methods that incorporate receptor flexibility, a deep learning approach can internalize valid cognate and noncognate binding modes corresponding to distinct receptor conformations, thereby learning to infer and account for receptor flexibility even on single conformations. ANPR significantly enriched pose quality in docking to cognate and noncognate receptors of the PDBbind v2019 data set. Improved pose rankings that better represent experimentally observed ligand binding modes improve hit rates in vHTS campaigns and thereby advance computational drug discovery, especially for novel therapeutic targets or novel binding sites.

Most macromolecular X-ray structures are determined from cryocooled crystals, but it is unclear whether cryocooling distorts functionally relevant flexibility. Here we compare independently acquired pairs of high-resolution data sets of a model Michaelis complex of dihydrofolate reductase (DHFR), collected by separate groups at both room and cryogenic temperatures. These data sets allow us to isolate the differences between experimental procedures and between temperatures. Our analyses of multiconformer models and time-averaged ensembles suggest that cryocooling suppresses and otherwise modifies side-chain and main-chain conformational heterogeneity, quenching dynamic contact networks. Despite some idiosyncratic differences, most changes from room temperature to cryogenic temperature are conserved and likely reflect temperature-dependent solvent remodeling. Both cryogenic data sets point to additional conformations not evident in the corresponding room temperature data sets, suggesting that cryocooling does not merely trap preexisting conformational heterogeneity. Our results demonstrate that crystal cryocooling consistently distorts the energy landscape of DHFR, a paragon for understanding functional protein dynamics.

Proteins and ligands sample a conformational ensemble that governs molecular recognition, activity, and dissociation. In structure-based drug design, access to this conformational ensemble is critical to understand the balance between entropy and enthalpy in lead optimization. However, ligand conformational heterogeneity is currently severely underreported in crystal structures in the Protein Data Bank, owing in part to a lack of automated and unbiased procedures to model an ensemble of protein-ligand states into X-ray data. Here, we designed a computational method, qFit-ligand , to automatically resolve conformationally averaged ligand heterogeneity in crystal structures, and applied it to a large set of protein receptor-ligand complexes. We found that up to 29 % of a dataset of protein crystal structures bound with drug-like molecules present evidence of unmodeled, averaged, relatively isoenergetic conformations in ligand-receptor interactions. In many retrospective cases, these alternate conformations were adventitiously exploited to guide compound design, resulting in improved potency or selectivity. Combining qFit-ligand with high-throughput screening or multi-temperature crystallography could therefore augment the structure-based drug design toolbox.