Ligand binding remodels protein side chain conformational heterogeneity

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 decrease 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 NMR 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.


INTRODUCTION
Ligand binding is essential for many protein functions, including enzyme catalysis, receptor activation, and drug response (Mobley and Dill, 2009). Ligand binding reshapes the protein conformational ensemble between the ligand-bound (holo) and unbound (apo) states, alternative conformations, or "harmonic disorder", which can be modeled by B-factors/atomic displacement parameters ( Figure 1A). Molecular dynamics experiments have demonstrated that if alternative conformations are not modeled correctly and consistently, then B-factors take on values that are not representative of the underlying conformational heterogeneity (Kuriyan et al., 1986;Kuzmanic et al., 2014). Moreover, B-factors incorporate many effects, including the biases and restraints of the refinement programs, modeling errors, crystal lattice defects, and occupancy changes of atoms. Therefore, consistently modeling X-ray structures as multiconformer models, with alternative side chain and backbone conformations, along with Bfactors, may better complement the view emerging from NMR and improve our understanding of the energetics of binding (van den . 70 Here, we examine how protein side chain conformational heterogeneity changes upon ligand binding by assembling a large, high-quality dataset of matched holo and apo X-ray crystallography structures. To integrate both harmonic and anharmonic disorder, we use a consistent multiconformer modeling procedure, qFit (Riley et al., 2021), and crystallographic order parameters (Fenwick et al., 2014). We test the hypothesis that ligand binding narrows the conformational ensemble, resulting in a decrease in heterogeneity of side chains in the holo structure compared with the apo structure. Our analysis reveals complex patterns of conformational heterogeneity that vary between and within proteins upon ligand binding.
Specifically, in proteins where binding site residues become more rigid upon ligand binding, 80 distant residues tend to become less rigid. This observation suggests that both natural and artificial ligands can modulate the natural composition of the protein conformational heterogeneity across the entire receptor to modulate the free energy of binding.

Figure 1. Representing structural data as multiconformer models. (A)
The grey outlines represent snapshots of the true underlying ensemble of the phenylalanine residue. The orange stick represents the residue modeled as a single conformer. The teal sticks represent the residue modeled as alternative conformers. The single conformer accounts for all heterogeneity in the B-factor, increasing the B-factor and reducing our ability to determine harmonic versus anharmonic motion. When a residue is modeled using alternative conformers, this heterogeneity is divided between harmonic heterogeneity, captured by the B-factors of each alternative 90 conformation and the anharmonic heterogeneity, captured by spread in coordinates between the alternative conformations. (B) To quantify the conformational heterogeneity of each residue, we used multi-conformer order parameters (Fenwick et al., 2014), which are the products of the ortho order parameter, which captures the harmonic or B-factor portion of each conformation and the angular order parameter, which captures the anharmonic portion or the displacement between alternative conformers. These are multiplied to produce the final order parameter (Methods). (C) The change in the number of alternative conformers (holo-apo) in binding site residues. In the re-refined dataset (orange), the majority structures have the same number of alternative conformers in the binding site, with the second most popular category gaining alternative conformers in the holo structure. In the qfit dataset (teal), the majority of structures lose an alternative conformer in the holo structure, with the second most common category being gaining an alternative conformer. (D) The differences in B-factors 100 (holo-apo) in the re-refined (orange) and qFit (teal) datasets. Overall, there was no significant difference in Bfactors between holo and apo structures in both the re-refined and qFit datasets.

RESULTS
Assembling the dataset To assess the differences in conformational heterogeneity upon ligand binding, we identified high quality, high resolution (2Å resolution or better) X-ray crystallography datasets from the PDB (Berman et al., 2000). We classified structures as holo if they had a ligand with 10 or more heavy atoms, excluding common crystallographic additives (Supplementary Table 1, Supplementary Figure 1). Structures without ligands, excluding common crystallographic additives, were classified as apo (Supplementary Figure 1A). We identified holo/apo matched 110 pairs by requiring the same sequence and near-isomorphous crystallographic parameters. Furthermore, we required the resolution difference between holo and apo pairs to be 0.1 Å or less, selecting representative apo structures to minimize the difference in resolution (Supplementary Figure 1B). This stringently matched ligand holo-apo dataset contained 1,205 pairs (Supplementary Table 2). We also used identical selection criteria to create a control dataset of 293 apo-apo pairs, taken from the set of holo/apo pairs (Supplementary Table 3).
Re-refining and qFit modeling of apo/holo pairs To minimize biases resulting from different model refinement protocols, we re-refined all structures using the deposited structure factors and phenix.refine (Liebschner et al., 2019). 120 The majority of structures in our re-refined dataset had less than 2% of residues modeled with alternative conformations, likely reflecting undermodeling of conformational heterogeneity represented in the PDB, based on prior literature (Lang et al., 2010). To more consistently assess conformational heterogeneity, we rebuilt all structures using qFit, an automated multiconformer modeling algorithm (Keedy et al., 2015;Riley et al., 2021) with subsequent refinement using phenix.refine (Liebschner et al., 2019). While qFit has biases, running all models through a consistent protocol will avoid manual biases that could creep into the holo or apo structures specifically. Additionally, by re-building each model as a multiconformer model, we were able to better distinguish the contributions of harmonic and anharmonic conformational heterogeneity across the structure (Figure 1A, B). All models went through 130 additional quality control, removing structures that resulted in large increases in R-free at each refinement step, high clashscores, or large RMSD between the pairs (Methods). This procedure resulted in 743 pairs (Supplementary Figure 2). Due to apo datasets serving as the reference state for multiple ligand bound structures, our dataset consists of 743 unique holo structures and 432 unique apo structures.
Properties of the apo/holo pairs The median resolution across our dataset was 1.6Å with a small trend towards improved (higher) resolution in the apo structure (0.01Å, median improvement (holo-apo); p=3.8x10 -20 , Wilcoxon signed rank test; Supplementary Figure 3). Across the dataset, 546 unique ligands were present in the structures, with 134 of these (e.g. NAG, AMP, etc) appearing in multiple 140 structures (Supplementary Figure 4A). The median number of ligand heavy atoms was 19, with only 10 very large ligands (>50 heavy atoms, e.g. Atazanavir; Supplementary Figure  4B). The proteins in the dataset represent 315 unique Uniprot IDs, with a bias towards enzymes that have been used for model systems for structural biology, including: Endothiopepsin (n=73 pairs), Lysozyme (n=62 pairs), Trypsin (n=48 pairs), and Carbonic Anhydrase 2 (n=46 pairs; Supplementary Figure 4C, D).

Conformational Heterogeneity across the Re-refined and qFit dataset
To determine the differences in conformational heterogeneity upon ligand binding in both the re-refined and qFit models, we assessed four commonly used metrics: the number of alternative conformers, B-factors (atomic displacement parameter), root-mean-square 150 fluctuations (RMSF), and rotamer changes.

Number of Alternative Conformations
Alternative conformations were modeled at low frequency in the re-refined dataset compared to the qFit modeled structures (1.7% vs. 47.8% of residues). In the re-refined dataset, there is a bias to increased modeling of alternative conformations in the holo dataset (50.5% gain vs. 29.8% loss), whereas more even representation was observed in the qFit dataset (44.3% gain vs. 54.8% loss; Supplementary Figure 5A). These results suggest that the trend of increased side chain conformational heterogeneity in PDB deposited structures may have its origin in human bias with more careful human attention to careful model building of binding site 160 residues in holo structures.
We next focused our analysis on binding site residues, defined as any residue with a heavy atom within 5Å of any ligand heavy atom. In the re-refined dataset, 23.9% of the matched pairs had a gain in alternative conformations in the holo model compared to only 19.3% losing an alternative conformer in the holo model, suggesting, counter-intuitively, that ligand binding increases local side chain mobility ( Figure 1C). However, in the qFit dataset, holo models tend to lose alternative conformations in the binding site residues (39.7% gain vs. 51.5% loss; Figure 1C).

B-factors 170
Next we explored the harmonic contribution to conformational heterogeneity as modeled by Bfactors on a pairwise, residue by residue basis. Across all residues in the re-refined dataset, Bfactors were slightly higher in holo models (0.31Å 2 , median difference (holo-apo); p=4.4x10 -208 , Wilcoxon signed rank test; Supplementary Figure 6A). In the qFit dataset, similar to the re-refined structures, holo residues had slightly higher B-factors (0.34Å 2 , median difference (holo-apo); p=5.6x10 -264 , Wilcoxon signed rank test; Supplementary Figure 5B, Supplementary Figure 6B). Of note, the B-factors in the qFit dataset are slightly smaller than the re-refined dataset (13.41Å 2 vs. 13.94Å 2 , average B-factors) reflecting the tendency for alternative conformation effects to be modeled as increased B-factors. When examining the binding site residues, there was no significant difference in B-factors between the holo and apo 180 models in both the re-refined (0.01Å 2 , median difference in B-factors; p=0.34, Wilcoxon signed rank test; Figure 1D) and qFit datasets (0.06Å 2 , median difference in B-factors; p=0.7, Wilcoxon signed rank test; Figure 1D, Supplementary Figure 6B). The lack of change in Bfactors close to ligands between the holo and apo models indicate that changes between the holo and apo B-factors are driven by signals distant from the binding site.

Conformational differences incorporating alternative conformations
Because of the low number of alternative conformers in the re-refined dataset, we only explored the anharmonic differences for side chains between the holo and apo models in the qFit dataset. First, to determine the extent of conformational change of alternative conformations, we compared the rotameric distribution of side chains. Side chain rotamer 190 changes between apo and holo structures have been reported to be very prevalent in single conformer models, with 90% of binding sites having at least one residue changing rotamers upon ligand binding (Gaudreault et al., 2012). To accommodate multiconformer models, we assigned all conformations to distinct rotamers using phenix.rotalyze. We classified each residue as having "no change" in rotamers if the set of rotamer assignments matched for the holo and apo residue. In binding sites, we also observed that "no change" was the most common outcome for residues (78.6%; Figure 2A). In the second largest category, "distinct", the holo and apo residue shared no rotamer assignments (15.5% of residues; Figure 2B).
A more complicated situation occurs when some, but not all, of the rotamer assignments are 200 shared across apo and holo residue. We classified 2.6% of residues as "remodeled -holo loss" ( Figure 2C) if distinct, additional rotameric conformations were populated in the apo residue only and 3.8% of residues as "remodeled -holo gain" (Figure 2D) if distinct, additional rotameric conformations were populated in the holo residue only. These results suggest a counterintuitive interpretation of binding site residues increasing their conformational heterogeneity upon ligand binding. However, a major potential confounder is that holo structures reflect an ensemble average of two compositional states (apo and holo) with alternative conformations representing the apo state at reduced occupancy, which we examined by subsetting the ligands based on relative B-factors (see below). A potential for a third category of remodeling, where both apo and holo residues share at least one 210 conformation and each have at least one additional conformation, did not occur in our dataset.
Across apo-holo pairs there was a large range of the percentages of binding site residues with the same rotamer classification in the pairs (23.2% to 100.0%), indicating that side chain remodeling can be quite variable ( Figure 2E). We found 11% of binding sites had all residues classified as "same" between pairs, consistent with a previous study that used single conformer models (Gaudreault et al., 2012). As an example of such a "pre-organized" binding site is Galectin-3 bound to thiodigalactoside (PDB: 5NFC, 4JC1; Figure 2F). In contrast, 67% of binding site residues have a rotamer status difference in transthyretin (PDB: 1CZR, 3CFN; Figure 2G), including a rotamer change in Leu101 to avoid a clash with the ligand. To compare the magnitude of fluctuations between alternative conformations, we calculated RMSF for all residues. This analysis suggested that, on average, apo residues have slightly greater conformational heterogeneity than holo residues (-0.006Å, mean difference of RMSF(holo-apo); p=3.7x10 -8 , Wilcoxon signed rank test; Supplementary Figure 6C). This trend was somewhat stronger in binding site residues (-0.02 Å, mean difference of RMSF(holoapo); p=4.5x10 -29 , Wilcoxon signed rank test; Supplementary Figure 6D). Our RMSF results 240 suggest that, on average, there is a slight decrease in heterogeneity upon ligand binding and that this reduction is most prevalent at residues distant from the binding site.
Collectively, these results do not conform to a simple model. There is a large amount of variability in the response across datasets and the median responses reveal only small biases. Nonetheless, considering those average responses, upon binding a ligand, the RMSF analysis suggests decreases in heterogeneity at the binding site, whereas the rotamer comparison has a slight bias to increased heterogeneity at the binding site, and B-factors only change at distant sites. One interpretation is that heterogeneity is reduced in binding site residues by a small number of anharmonic conformational changes, as observed by the RMSF reduction, paired 250 with an increase in harmonic fluctuations far away, as observed by an increase in the Bfactors. However, it is difficult to interpret these changes separately, as conformational heterogeneity is a combination of both harmonic and anharmonic motion and there is potential degeneracy in modeling alternative conformations, even with qFit (van den Bedem et al., 2009). Therefore, we moved to using an integrated measurement of order parameters that can account for these complications (Fenwick et al., 2014).
Order parameters integrate both harmonic and anharmonic conformational heterogeneity.
To integrate the anharmonic fluctuations between alternative conformers with the harmonic fluctuations modeled by B-factors (Kuzmanic et al., 2014), we used a crystallographic order 260 parameter ( Figure 1A) (Fenwick et al., 2014). While order parameters are traditionally used in NMR or molecular dynamic simulations, they can be calculated for multiconformer X-ray models and, in some cases, show reasonable agreement with solution measures (Fenwick et al., 2014). We focused on the order parameters of the first torsion angle ( 1) of every sidechain for all residues except for glycine and proline. Order parameters are measured on a scale of 0 to 1, with 1 representing a fully rigid residue and 0 representing a fully flexible residue. Below, we analyze the differences in normalized order parameters between paired residues (Methods).
As an additional control, we compared our apo/holo dataset to a dataset of apo/apo pairs. In 270 examining the differences in order parameters, both in the holo/apo pairs and the apo/apo pairs, there are no large differences in conformational heterogeneity, as indicated by a median difference in order parameters of approximately 0. However, in the holo/apo pairs there is a much wider range of order parameter differences, indicating that ligand binding impacts conformational heterogeneity beyond experimental variability (p=3.4x10 -17 , individual Mann-Whitney U test; Figure 3A).
Next, to examine whether different regions of the protein were driving this higher variability, we compared the differences in order parameters among binding site residues, within 5Å of any ligand heavy atom, compared to a control dataset which matched the number of, type and 280 solvent exposure within the protein for each binding site residue. In binding site residues, the holo structures had a slightly, but significantly, increased order parameters, suggesting reduced conformational heterogeneity compared to the control dataset (0.034 vs. 0, median difference (holo-apo) order parameter; p=3.4x10 -7 , individual Mann-Whitney U test; Figure 3B). While there is a larger range of responses, this indicated that, in general, binding site residues become more rigid upon ligand binding.

Spatial distribution of conformational heterogeneity changes
Based on the large range of order parameter differences we observed across the protein, along with the decrease in heterogeneity localized to binding site residues, we next explored the relationship between changes in heterogeneity in binding site residues and the rest of the protein. The difference in order parameters between the holo and apo models were correlated 310 in both the binding site and distant residues (Supplementary Figure 7A), indicating that ligand binding generally caused global changes to flexibility. Given the average rigidification of the binding site residues (Figure 3B), these results predict a general trend of decreased conformational heterogeneity in the ligand binding site would be associated with a relative increase in conformational heterogeneity at distant sites in the protein. This pattern suggests that the residual change in heterogeneity (the difference between the average order parameter of the distant residues and the average order parameter of the binding site residues) should be inversely related to the change in the binding site residues: more rigidified binding sites will have more flexible than expected distant sites, and vice versa. Indeed, on a protein-by-protein basis, the relationship between binding site residues and residual changes at distant sites 320 follows this trend (Supplementary Figure 7B). Consistent with studies suggesting significant residual conformational heterogeneity in folded buried residues (Wong and Daggett, 1998) and the potential for those buried residues to change heterogeneity upon ligand binding (Moorman et al., 2012), this trend is even stronger in residues that were more than 10Å away from any heavy atom in the ligand and less than 20% solvent exposed (slope=-0.41, r 2 =0.46; p=5.1x10 -50 , two-sided t-test; Figure 3C). This indicates that proteins that lose conformational heterogeneity in the binding site are associated with a relative increase in conformational heterogeneity in distant, non-solvent exposed residues.
There are three likely origins of this effect. First, this may reflect a feature of the distribution of 330 order parameters around the mean value within each protein. Second, this may reflect a topological feature of protein packing, whereby packing optimization of certain areas of a protein decreases the optimization of other parts of the protein (Bromberg and Dill, 1994). Third, this may reflect the stabilization of certain conformations in a ligand bound protein. As a control for these effects, we compared the residual order parameter differences between the buried, non-solvent exposed residues and the binding site residues in apo-apo pairs. Globally the trends were similar, but weaker in both correlation and magnitude (slope=-0.28, r 2 =0.20; Supplementary Figure 7C). Therefore, we interpret the trend we observe as mainly based on protein topology, specifically that proteins have areas where there are less efficiently packed alternative conformers, likely to enable entropic compensation across the protein during 340 various functions, including ligand binding. We interpret that the stronger signal we observed in the holo-apo dataset is due to the ligand perturbation, which is also reflected in the median rigidification of binding site residues ( Figure 3B). We hypothesize that we are observing this innate protein property being used, specifically optimizing the binding site residues to bind a ligand, while decreasing the optimization elsewhere in the protein.
As an example to visualize this trend, we mapped the change in order parameters onto the structure of the human ATAD2 bromodomain (PDB ID:5A5N). In ATAD2, the binding site residues rigidify upon ligand binding whereas the majority of distant residues are more heterogeneous compared to the binding site residues ( Figure 3D). Specifically, this difference 350 is greatest between binding residues and non-solvent exposed residues. However, as in the global analysis, this example demonstrates there is a large range of changes in binding site order parameters, consistent with NMR examples that show a heterogeneous response both close to and far from ligands (Caro et al., 2021;Moorman et al., 2012).

Ligand properties influence conformational heterogeneity
Next we investigated how ligand properties impact the conformational heterogeneity of binding site residues. For ligand properties dictated by the size of the ligand (number of rotatable bonds and number of hydrogen bonds) we normalized these metrics by the molecular weight of the ligand. For each property, we compared the highest and lowest quartiles by both the absolute order parameters of the holo structure and the order parameters differences between 360 holo and apo pairs. No significant associations existed when comparing the differences between holo and apo order parameters, but the characteristics of the holo binding site and the rotamer changes were correlated with ligand properties in several cases.
We hypothesized that ligand properties associated with increase ligand dynamics, including more rotatable bonds, higher lipophilicity (logP), fewer hydrogen bonds, and more heavy atoms would be associated with increased conformational heterogeneity (an increase in absolute order parameters or a smaller difference between the apo and holo order parameters; Wicker and Cooper, 2015). While molecules with fewer rotatable bonds (lower quartile: <2 (n=134) vs. upper quartile: >6 (n=134)) were indeed associated with more rigid binding sites 370 (lower quartile: 0.83 vs. upper quartile: 0.81, individual Mann Whitney U test), this was not significant. However, higher numbers of rotatable bonds were associated with a lower number of same rotamers between the apo and holo binding site residues (88% vs. 80%, percentage same rotamer; p=6.0x10 -6 , individual Mann Whitney U test; Supplementary Figure 8A). Increased lipophilicity (logP, upper quartile: <0.04 (n=134) vs. lower quartile: >2.69 (n=134)), was significantly associated with a more flexible binding site (0.84 vs. 0.79, median order parameters; p=7.5x10 -6 , individual Mann Whitney U test; Figure 4A). Previous studies have indicated that increased lipophilicity generates more nonspecific binding interactions (Olsson et al., 2008). Larger compounds (upper quartile: >26 heavy atoms (n=134) vs. lower quartile: <13 heavy atoms (n=134)) are also associated with more flexible binding sites (0.79 vs. 0.83, 380 median order parameter; p=0.0001, individual Mann Whitney U test; Figure 4B). Large compounds, thus a larger ligand surface area, are associated with more nonspecific binding interactions, which is compatible with increased protein conformational heterogeneity. Finally, more total hydrogen bonds per heavy atom(upper quartile: >0.47 (n=134) vs. lower quartile: <0.25 (n=134)) are associated with more rigid binding sites (0.84 vs. 0.79, median order parameter; p=5.9x10 -5 , individual Mann Whitney U test; Figure 4C). This trend holds even when examining hydrogen bond donors or acceptors separately.
From these results, an intuitive general picture emerges where more specific interactions, such as hydrogen bonds, are correlated with more rigid binding site residues, whereas the more 390 non-specific interactions are correlated with more flexible binding site residues. There is also a wide range of deviation from this general picture, likely reflecting that natural and artificial optimization of ligands is based on free energy, not any specific thermodynamic component or interaction type. These trends emphasize the need to monitor both the impacts of ligands on specific interactions with the protein along with conformational heterogeneity of the protein.
Additionally, these results suggest that specific interactions can be tuned to rigidify a binding site. Paired with our findings of the relationship between order parameters in binding site and distant residues, ligand impacts are likely propagated throughout the protein. Ligands with more specific interactions, thus a less flexible binding site, will likely have a corresponding increase in conformational heterogeneity far away from the binding site. 400 Reduced ligand occupancy and conformational heterogeneity One potential confounder for quantifying the change in conformational heterogeneity of binding site residues is that the ligands may not be fully occupied in the crystal. There were 193 structures with ligands with alternative conformations or partially occupied ligands in our datasets ( Figure 4D). Of these 193, 125 ligands had less than full occupancy, whereas 68 had alternative conformations that amounted to full occupancy. The vast majority of ligands (n=425) were modeled originally with full occupancy. Fully occupied ligands were associated with more rigid binding sites than partially occupied ligands or ligands with alternative conformers (0.84 vs. 0.79, mean order parameters of binding site residues; p=2.9x10 -7 , individual Mann Whitney U test; Figure 4D). There was no difference observed between the 410 partially occupied ligands and ligands with alternative conformers (p=0.15, individual Mann Whitney U test; Supplementary Figure 8B). We also explored if partially occupied ligands were associated with more rotamer changes between holo and apo pairs, but no significant difference existed (80% vs. 85%, median percentage of the same rotamer; Supplementary Figure 8B).
While the scattering contributions of B-factor and occupancy changes are subtle (but distinct), most models likely include true occupancy changes as elevated B-factors. We observed a wide range of average ligand B-factors and, as expected, a lack of correlation between the ligand B-factors and ligand occupancy (Bhat, 1989;Carugo, 1999;van Zundert et al., 2018). 420 As a proxy for likely partially occupied ligands, we normalized the ligand B-factor by the mean C-alpha B-factor to identify ligands with higher B-factors than expected (Supplementary Figure 8C). We examined the outer two quartiles of the normalized ligand B-Factors (>0.016 vs. <0.005, median normalized B-factor). In these "likely partially occupied" ligands, we observed greater conformational heterogeneity (0.86 vs 0.80, mean order parameter; p=1.6x10 -11 ; individual Mann Whitney U test, Figure 4E). In structures with modeled partially occupied ligands and likely partially occupied ligands, we learned that binding site residues tend to have more apparent conformational heterogeneity, likely due a combination of compositional and conformational heterogeneity.

Conformational heterogeneity for multiple ligands to CDK2
To better understand our findings in the context of multiple, diverse ligands binding to a single receptor, we examined Cyclin-dependent Kinase 2 (CDK2), a cyclin kinase family that regulates the G1 to S transition in the cell cycle. Our dataset contains 13 protein-inhibitor complexes, including both type I and type II inhibitors, all of which share the same apo model (PDB ID: 1PW2). We hierarchically clustered the residues and ligands by difference in order parameters between the holo and apo models, identifying three distinct clusters of residues. The first cluster (blue, Figure 5A,B), consisting of 13 residues, are rigidified upon ligand binding. This cluster included residues scattered throughout both the N-and C-lobes of CDK2 460 that rigidified upon ligand binding. This dispersed pattern is similar to the trend of rigidification that is observed by NMR in PKA upon substrate binding (Kim et al., 2017). Two notable residues in this cluster, Glu127 and Val18, contact the inhibitors. Upon ligand binding, Val18 transitions from multiple conformers to a single conformation. Glu127 has a similar conformation in the apo and type II structures of two distinct alternative side chain rotamers, whereas in the type I inhibitor structure, the alternative conformers cluster around the same rotamer (Supplementary Figure 9A,B).

A B C D P-Loop
The second cluster (salmon, Figure 5A,B), consists of 14 residues that increase flexibility upon ligand binding. The majority of these residues connect the p-loop and the activation loop 490 ( Figure 5B). The electron density is very weak for many of these residues in most of the holo structures, driving their modeling in multiple conformations and elevated B-factors (Supplementary Figure 9D). The third cluster (dark red, Figure 5A,B) is comprised of five residues that became more flexible in all, but two holo datasets, which are the only type II inhibitors in the dataset. These were all located on the activation loop of the kinase ( Figure  6D). As type II inhibitors, the two molecules [PDB: 1PXI (ligand: CK1) and PDB: 3QQL (ligand: X03)] bind the DFG out conformation present in the apo dataset (PDB: 1PW2) and do not have as drastic of a rigidifying effect as the type I inhibitors. Notably, these two inhibitors were also smaller than the type I inhibitors and the reduced contacts may also drive some of this effect. 500 The multiconformer models also provide a structural rationale for these changes. The differences in DFG conformation change the contacts with the P-loop, which allow for greater side chain flexibility in the "up" form compatible with type I inhibitors. The interface between the P-loop and the activation loop is weaker and residues such as Tyr155 adopt multiple conformations. At the base of the activation loop, Thr161, a critical phosphorylation site, changes conformation, with a rigidifying effect common to both type I and II inhibitors ( Figure 5C). The conformation of Thr161 found in the type II inhibitors overlap, with one of the conformations populated in the multiconformer apo model. In contrast, the type I inhibitors adopt a distinct conformation. This case study highlights how modeling information present in the density can reveal changes beyond those in single conformer structures.

DISCUSSION
By creating a large dataset of stringent matched pairs of apo and holo multiconformer models, we identified a pattern of conformational heterogeneity consistent with smaller scale studies of individual proteins (Caro et al., 2021). In general, we found that binding site residues tend to become more rigid upon ligand binding. However there was a large range of effects, including many sites becoming more flexible when bound to a ligand. The trends suggest that disorderorder transitions between binding site residues and distant residues are common and potentially a selected property of many proteins (Kim et al., 2017). Specifically, our data suggests that some of the entropy lost by the rigidification incurred by binding site residues 520 upon ligand binding can be compensated with an increase in disorder in distant residues. This finding generalizes the phenomenon has been observed in single protein analyses with X-ray, NMR, and MD simulation (Gohlke et al., 2004;Moorman et al., 2012;Wang et al., 2019). Both theoretical and experimental analyses suggest that the relationship between local packing optimization and small voids that permit alternative conformations will be key to predictably mapping this relationship (Bromberg and Dill, 1994;Caro et al., 2021). Using temperature or pressure as perturbations during X-ray data collection can help to further map the connection between packing "quality" and side chain conformational heterogeneity in greater detail (Caro et al., 2021). While NMR order parameter studies only take into account movement that is shorter than the tumbling time for the protein (Hoffmann et al., 2021;Gangé et al., 1998), our 530 results are insensitive to timescale. In addition, it is quite likely that our use of cyro-cooled structures causes an underestimate of the heterogeneity occurring in these datasets (Fraser et al., 2011).
We observed a complex interplay between conformational change and dynamics in our analysis of 13 inhibitor-bound datasets of the kinase CDK2, in the same crystal form and space group. The ability to explore one protein with multiple ligands highlights the utility of crystal systems amenable to isomorphous soaking or co-crystallization (Steuber et al., 2006). We identified differences in conformational heterogeneity between type I and type II inhibitors that can be classified along with well-known changes, such as differences in the DFG motif. 540 Tuning distal site dynamics may be a viable strategy for modulating the affinity of kinase inhibitors and affect the pattern of protein-protein interactions on distal surfaces, which is of critical importance in CDK inhibitor development (Jhaveri et al., 2021;Wood and Endicott, 2018).
We note that our work is not sensitive to many facets of the complex changes associated with ligand binding (Mobley and Dill, 2009). Our stringent resolution matching criterion may also render us blind to the most severe effects on conformational heterogeneity, whereby ligand binding causes a more widespread change leading to a loss or gain of diffraction power. In addition, water molecules play an important role in ligand binding, both in the release of 550 ordered water molecules contributing to binding through entropy and in forming specific interactions (Breiten et al., 2013;Verteramo et al., 2019). Additionally, ligand conformational heterogeneity has been highlighted by several recent studies (Caldararu et al., 2021;Jain et al., 2020;van Zundert et al., 2018). Another caveat in our analysis is the limitations of qFit modeling for modeling extensive backbone heterogeneity into weak electron density. Ensemble modeling methods, which leverage molecular dynamics for sampling and use a different model representation may be more appropriate for examining these systems (Burnley et al., 2012;Eshun-Wilson et al., 2019). Future work, integrating the conformational heterogeneity of the protein, ligand, and water molecules will create better predictions and explanations of the energetics of binding. In addition, this would allow us to interpret the impact 560 of specific interactions and alterations on both the entropy and enthalpy of all components of the system.
Our study, as well as previous NMR studies (Caro et al., 2017;Frederick et al., 2007), only leverage a limited set of side chain dihedral angles. However, comparisons with molecular dynamics simulations suggest that small sets of side chain dihedrals alone may be representative of the overall changes in dynamics of the system Sharp, 2018, Fleck et al., 2018). What is the thermodynamic impact of restricting side chain conformational heterogeneity? Protein folding studies and theory indicate that restricting the rotamer of even a single side chain can incur an entropic penalty of binding of ~0.5kcal/mol (Doig and Sternberg,570 1995). While we observe many such restrictions in binding sites due to specific interactions with ligands, our data point to corresponding changes away from the binding site that help balance this cost. Overall, the median increase in rigidity we observe in binding site residues (0.03 order parameter increase) would create an energetic penalty of approximately ~0.1-0.5kcal/mol (Wand and Sharp, 2018), with outliers having even larger thermodynamic consequences. Given the constraints of maintaining a folded conformational ensemble upon ligand binding, it is likely that ligand binding generally acts to restrict degrees of freedom locally and that protein topological constraints favor increased motion in distal regions (Bromberg and Dill, 1994). This overall effect likely manifests because optimizing affinity is desirable for medicinal chemistry and for the selective pressures experienced by many proteins. Such 580 optimization is insensitive as to whether the free energy is driven enthalpically or entropically. However, given the attention paid to designing and optimizing enthalpic interactions, there is likely unleveraged potential in optimizing the entropic component as well. As more sophisticated models of conformational heterogeneity are created and validated (Rosenbaum et al., 2021) the strategy of rationally tuning conformational heterogeneity to improve binding affinity may be an attainable design strategy.

Dataset:
Our dataset was compiled using a snapshot of the PDB (Berman, 2002) in September 2019, containing 156,187 structures. We then selected structures that had a resolution better or equal to 2Å (n= 64,557). 590 We also excluded any structure that contained nucleic acids (n=2,280) or covalently bound ligands (n=1,030). We identified holo structures(n=30,530), defined as those that contained at least one ligand, defined as any HETATM residue with 10 or more heavy atoms, excluding common crystallographic additives.
To create apo/holo pairs, we took each holo structure and compared it to each potential apo structure (n=30,717), defined as structures without a ligand bound. A pair was defined according to the following criteria: -same space group -exact sequence or exact sequence after removing the first or last five base pairs of either structure 600 -a resolution difference between the two structures less than 0.1Å -dimensions of unit cells do not differ by more than 1Å -angles of the unit cells do not differ by more than 1 degree This gave us 15,214 pairs. We then subsetted this list down to provide only one apo structure per holo structure, based on the smallest resolution difference. This produced a final pair set of 1,205 with 1,143 unique structures (Supplementary Table 2).
We also created a pairset with 458 unique apo/apo pairs using the same criteria as the apo/holo pairset (Supplementary Table 3). 610

Refinement:
We re-refined all structures using phenix.refine (https://www.phenix-online.org/documentation/ reference/refinement.html). This was done using phenix version 1.17.1-3660. We performed anisotropic refinement on all pairs where both PDBs had a resolution better than 1.5Å. All other refinement was run isotropically. Refinement used the following parameters: • Refine strategy: individual sites + individual adp + occupancies • Number of macro cycles: 8 • NQH flips: True • Optimize xyz weight: True • Optimize adp weight: True 620 • Hydrogen refine: Riding We removed 102 structures because of incompatibility with our re-refinement pipeline due to breaks in the protein chain or ligand incompatibility. We removed 88 structures where the R-free increased by >2.5% compared to the value reported in the PDB header (Supplementary Table 1; Supplementary  Figure 2A).

Running qFit:
qFit-3.0 (Riley et al., 2021) (version 3.2.0) was run using a composite omit map and the re-refined structure on the default parameters(https://github.com/ExcitedStates/qfit-3.0/). We ran qFit on Amazon Web Services (AWS). We used an auto scaling cluster of images controlled by the scheduler via 630 ParallelCluser. Please see the qFit github for a script that will install qFit on AWS's default OS image, using conda to install its dependencies.
After qFit, we re-ran refinement as suggested by qFit-3.0. Briefly, this involves three rounds of refinement. The first refines coordinates only, the second goes through a cyclical round of refinement until the majority of low occupancy conformers are removed, and the third refinement polishes the structure, including hydrogens. The script used for post qFit refinement can be found here: https://github.com/ExcitedStates/qfit-3.0/blob/master/scripts/post/qfit_final_refine_xray.sh. We removed 100 structures because of incompatibility with refinement after qFit rebuilding.
Quality Control: 640 From our original dataset (n=1,205 pairs), we removed 28 apo structures that had a crystallographic additive or amino acid in the binding site that partially overlaid with the holo structure. We set a minimum ligand occupancy threshold of 0.15, which did not remove any pairs from our dataset. Chains were renamed according to the difference in distance between the two chains. We also re-numbered each chain based on the apo structure. We then superimposed the two structures using the pymol align function. We measured the alpha carbon root mean squared difference (RMSD) between the two structures as well as the difference in just binding site residues. Structures were removed if the mean RMSD of the entire structure was greater than 1Å or if the mean RMSD in the binding site residues were greater than 0.5Å. We removed two pairs based on these RMSD criteria. 650 We also assessed the difference in R-free values for each refinement step (before/after qFit). If the post refinement R-free value was 2.5% larger than the pre refinement R-value, then the structure was removed (n=85, 77 structures removed; Supplementary Figure 2A,B) . Additionally, we compared the final R-free values between apo and holo pairs, removing pairs with R-free values with more than a 5% difference (n=16 pairs removed; Supplementary Figure 2C). We ran the clashscore function out of Molprobity (Williams et al., 2018) to identify severe clashes in our dataset. We removed any structures with a clashscore greater than 15, removing 52 structures. After filtering out pairs that failed our quality checks, our dataset contained 743 matched apo/holo pairs.

Alternative Conformations:
Side chains were considered alternative conformers if there was at least one atom that was modeled 660 with an alternative conformer. Our re-refinement procedure changes the occupancy, coordinates, and B-factors of these conformations, but it does not add or delete conformations.

B-factors:
B-factors were assessed on a residue basis by averaging the B-factors of all heavy atoms for each residue. For residues with multiple conformations, we took the mean B-factor for all heavy atoms in each side chain, weighted by occupancy. For structures modeled anisotropically, we used the isotropic equivalent B-factor from phenix.

Root Mean Squared Fluctuation (RMSF):
670 RMSF was chosen over root mean squared deviation as many alternative conformers were predicted to have the same occupancy, thus making it difficult to define which was the main conformer. RMSF was measured for each residue based on all side chain heavy atoms. RMSF finds the geometric center of each atom in all alternative conformers. It then takes the distance between the geometric mean of each conformer's side chain heavy atoms and the overall geometric center. It then takes the squared mean of all of those distances, weighted by occupancy.
Order Parameters: Order parameters were measured for each residue (except proline and glycine) by taking into account both the angle of alternative conformers (s2angle), by measuring the chi1 angle, and the B-factors of 680 alpha or beta carbons along with an attached hydrogen(s2ortho) (Fenwick et al., 2014). To account for differences in B-factors as resolution changes, we investigated the correlation between order parameters in 32 apo lysozyme structures ranging in resolution from 1.1 to 2Å. We optimized the s2ortho parameter by looking for the normalization that would provide a slope closest to one and have the smallest root mean squared error. We normalized the s2ortho portion using the following equation: 2 ℎ *+,-./0123 = 2 ℎ * ./:;. =.,>+* /10( ) The final order parameter reported in the paper is: 2 = 2 ℎ *+,-./0123 * 2 Rotamer Analysis: Rotamers were determined using phenix.rotalyze (Williams et al., 2018) with manually relaxing the 690 outlier criteria to 0.1%. Each alternative conformation has its own rotamer state. Rotamers were compared on a residue by residue basis between the holo and apo, taking into account each rotamer state for each alternative conformation. Residues were classified as "no change" if rotamers matched across the apo and holo residue, "distinct" if the matched residue shared no rotamer assignments. Residues were classified as "remodeled-holo loss" if distinct, additional rotameric conformations were populated in the apo residue only, and "remodeled -holo gain" if distinct, additional rotameric conformations were populated in the holo residue only.
an alpha of 0.05, as we were testing 10 hypotheses leaving us with a corrected significance value of 0.005.
Occupancy of the ligands were taken directly from the PDB file and correspond with the ligand occupancy from the deposited structure. Ligand B-factors were normalized using the mean alpha carbon B-factor of all residues in the structure. 710 If there were multiple ligands of interest in a structure, we looked at the properties of the ligand and surrounding protein residues in chain A or in the lowest alphabetical chain.

Protein Type Analysis:
Protein names and enzyme names were extracted from Uniprot(UniProt Consortium, 2015). Names and properties were connected using PDB IDs.

Statistics:
Paired Wilcoxen test was used for all matched properties (comparing holo v. apo matched residues or 720 structures). Individual Mann-Whitney U test was used for all non-match properties, including ligand properties. Two-sided t-test was used to compare the significance of the slopes. Code: Code can be found in the following repositories: -