Many proteins experience widespread molecular motions that occur on the s-ms timescale that are important for their function. [1][2][3][4][5] These motions can include protein sub-domains, active site loops, and folding-unfolding transitions. [6] In addition, the phenomenon of allostery [7,8] that relays signals between two distinct protein sites (classically, ligand binding sites) often relies on networks of flexible amino acid residues to enable this communication. [9][10][11][12][13] Specific examples include disease-causing mutations that can alter the biological function of proteins and enzymes by disrupting the necessary global motions. [10,14] Despite the biological importance and sustained study a unified mechanistic description of how allosteric signals are relayed through the protein matrix remains in flux and is likely dissimilar for different proteins. [15] One experimental method that can characterize molecular motions and their role in conformational changes and allostery is solution NMR spectroscopy. Solution NMR spectroscopy can directly detect subtle, atomic-level structural and dynamical perturbations that mediate conformational changes including allosteric processes. [16][17][18][19] These structural fluctuations can be detected using Carr-Purcell-Meiboom-Gill (CPMG) relaxation dispersion experiments, [20,21] which provide residue-specific relaxation dispersion curves that are used to quantify kinetics (exchange constant, kex) and thermodynamic (populations, pA and pA') parameters for exchange processes occurring on the s-ms timescale. [3] The residue-specific nature of the information obtained from relaxation dispersion experiments also presents problems in interpretation. One would like to know whether multiple residues are moving together (in concert) or independently. Historically, this was assessed by comparing kex values and using some (often arbitrary) criteria, i.e., determining whether kex values are similar enough to warrant grouping into a single motional process. Furthermore, when many residues show promising relaxation dispersion curves, grouping them by similar kex values is a difficult task and often very tedious to perform manually. Typically, dynamic parameters are determined for individual sites with high precision and then compared to suggest potential correlations. However, such manual analysis is prone to bias of the individual researcher.
Here, we propose an alternative method to automatically cluster residues undergoing synchronous motions (motions with the same kex values and populations) in proteins. This method builds upon a previous approach we described to identify the collective dynamic response for the catalytic subunit of the cAMP-dependent protein kinase A (PKA-C) to ligand binding. [22] In this revised method, we fit the CPMG relaxation dispersion data (either backbone amides or sidechain methyl groups) for each residue pair in the protein using the Bloch-McConnell equations, [23] and utilize the goodness of the pairwise fitting as a metric to build two-dimensional dynamic correlation (SyncDyn) maps and cluster residues whose motions are fluctuating synchronously.
We applied our method to two proteins with well-characterized millisecond motions, PKA-C in different ligated forms and the T17A mutant of the ribonuclease A (RNase A). We show how ligand binding or single mutations modulate the organization of the internal dynamic network of these enzymes.
To fit the CPMG dispersion curves, we consider a single residue undergoing a two-site exchange between a ground (A) and higher energy (A') state:
where kf and kr are the forward and reverse rate constants of the exchange process, respectively.
The evolution of the magnetization for this system can be written as
where MA and MA' are the transverse magnetizations for A and A', respectively, 𝑅 2𝐴 0 and 𝑅 2𝐴′ 0 are the intrinsic transverse relaxation rates for A and A', and ∆𝜔 is the chemical shift difference between sites A and A' expressed in angular frequency units (s -1 ). If pA >> pA' and we if we also assume that 𝑅 2𝐴 0 = 𝑅 2𝐴′ 0 = 𝑅 2 0 , the solution of equation ( 1) is
For the constant relaxation time CPMG relaxation dispersion (RD) experiment, the transverse relaxation rates (R 2,Eff ) are determined from the observed resonance intensities using the following expression:
where TCPMG is the total CPMG relaxation delay, 𝐼 𝐴 is the observed peak intensity after the CPMG sequence (TCPMG ≠ 0) and 𝐼 𝐴 0 is the intensity of the reference spectrum (TCPMG = 0). The RD profile results from a series of 𝑅 2,𝐸𝑓𝑓 at different values as a function of the pulse frequency (CPMG)
where CPMG = 1/(2CPMG), in which CPMG is the time interval between 180º pulses in the CPMG element. To construct the SyncDyn map, we begin with equation (3) and use the following steps:
Step 1. 𝑅 2,𝐸𝑓𝑓 calculation for individual residues. After collecting the CPMG RD data at a single or multiple B0 fields, we calculate the 𝑅 2,𝐸𝑓𝑓 values using equation 3 for different CPMG values.
Step 2. Identification of all residue pairs. We identify all non-equivalent residues pairs, which are N × (N -1) / 2 where N is the number residues in the protein.
Step 3. Pairwise fitting of all RD profiles. For a pair of protein residues, 𝑖 and 𝑗, the two RD profiles at a single B0 field, 𝑅𝐷 𝑖 and 𝑅𝐷 𝑗 , depend on 𝑘 𝑓 𝑖 , 𝑘 𝑟 𝑖 , 𝑅 2 𝑂𝑖 , ∆𝜔 𝑖 and 𝑘 𝑓 𝑗 , 𝑘 𝑟 𝑗 , 𝑅 2 𝑂𝑗 , ∆𝜔 𝑗 , respectively. To fit the RD profiles of the residue pair simultaneously, we assume that the kinetic constants, kf and kr, are the same for both residues: 𝑘 𝑓 𝑖 = 𝑘 𝑓 𝑗 = 𝑘 𝑓 (4) 𝑘 𝑟 𝑖 = 𝑘 𝑟 𝑗 = 𝑘 𝑟 (5) Therefore, 𝑅𝐷 𝑖 and 𝑅𝐷 𝑗 can be fitted using the Bloch-McConnell equation, combining all the parameters , 𝑘 𝑓 , 𝑘 𝑟 , 𝑅 2 0 𝑖 , 𝑅 2 0 𝑗 , ∆𝜔 𝑖 , ∆𝜔 𝑗 , into a single set. The resulting fitted RD profiles, 𝑅𝐷 𝑓𝑖𝑡 𝑖 and 𝑅𝐷 𝑓𝑖𝑡 𝑗 , are then used to assess the goodness of the fit with the experimental RD data for residues i and j.
Step 4. Calculating the goodness of the pairwise fitting using R-squared. To assess the goodness of the pairwise fitting, we utilized the R-squared value between the experimental and fitted RD profiles using the following equation [24]:
Where
and
SSR is the regression sum-of-squares and SST total sum-of-squares, n is the number of CPMG points, 𝑅𝐷 and 𝑅𝐷 𝑓𝑖𝑡 are experimental and fitted dispersion profiles. We calculated the R-squared values for each residue across the various Bo fields and identified the minimum R-squared value among all the dispersion curves and all residue pairs. . The R-squared ranges from 0 to 1. A high R-squared value for a residue pair indicates that the initial assumption is highly probable, i.e., the residues are likely to exchange between two states with similar frequency. A low value of Rsquared indicates that the motions of the two residues occur on a different timescale (Fig. 1). The final SyncDyn map consists of a matrix of R-squared values, where each element (i, j) represents the R-squared value between the residue pairs.
To validate this approach, we considered two residues, X and Y, exchanging with their respective minor states X' and Y'. We simulated two hundred synthetic CPMG dispersion curves for X, assuming exchange rates ranging between 0 and 5000 s -1 . For each dispersion curve of X, we simulated 120 CPMG dispersion curves for residue Y. The exchange rate of Y was selected in such a way that the difference in kex between X and Y would vary uniformly from 0 to 3000 s -1 .
The population of X' and Y' were assumed to be 0.02 for all the simulations (Fig. 2).
We then calculated the pairwise R-squared values for the 200 x 120 pairs and plotted as a function of the exchange rate for X (𝑘 𝑒𝑥 𝑋 ) on the horizontal axis and the difference between the exchange rates of X and Y (∆𝑘 𝑒𝑥 = 𝑘 𝑒𝑥 𝑋 -𝑘 𝑒𝑥 𝑌 ) on the vertical axis (Fig. 2). We found that the Rsquared values decrease with increasing ∆𝑘 𝑒𝑥 , with the most reliable correlations occurring for exchange kinetics between 200 and 4000 s -1 , as expected for typical CPMG experiments. [3] On the other hand, faster exchange rates give rise to flat dispersion profiles, which may result in artifactual correlations in the SyncDyn maps, especially in the presence of experimental noise.
To test this, we compared two simulated sets of data using random noise of 2 and 5% (Fig. 2a and 2b). We found that as the noise amplitude increases, the R-squared values significantly decrease, indicating that for reliable correlations it is critical to have high-quality relaxation dispersion data with low noise levels, especially for exchange rates near the fast limit as suggested previously [25].
Validation of the protocol using a synthetic CPMG dataset. To test our clustering algorithm, we simulated a set of synthetic CPMG data for a hypothetical twenty amino acid peptide with two clusters of three residues, i.e., cluster 1 with kf and kr of 25 and 500 s -1 , and cluster 2 with kf and kr of 100 and 2000 s -1 (Fig. 3a). As shown in Fig. 3b, our algorithm identifies the correlations between the two clusters of residues that share similar exchange kinetics as indicated by the off-diagonal elements in the 2D SyncDyn map and illustrated in the dendrogram (Fig. 3c).
We applied our approach to analyze two dynamic enzymes (Fig. S1), PKA-C and RNAse A.
PKA-C is a ubiquitous kinase responsible for phosphorylating numerous cell substrates, [26,27] while RNAse A catalyzes the transphosphorylation and hydrolysis of single-stranded RNA. [28] Both enzymes exhibit critical conformational fluctuations in the millisecond time scale during their catalytic cycles. We reevaluated their NMR relaxation dispersion data to identify potential clusters of residues that undergo synchronous motions.
Clustering synchronous residues of PKA-C. We first examined the CPMG relaxation dispersion data for the methyl group of the aliphatic side chains of PKA-C. [22,[29][30][31] This bilobate enzyme undergoes open-to-close transitions through its catalytic cycle. [32] The apo enzyme binds the nucleotide and adopts a partially closed conformation with a higher affinity for the substrate (positive cooperativity). [33] After substrate binding, the chemical step (phosphoryl transfer) is relatively fast (500 s -1 ), and the product release constitutes the rate-determining step of the catalytic cycle. [34][35][36] We previously found that synchronous and asynchronous fluctuations are hallmarks for the formation of the dynamically committed state and for product release, respectively. [22,36] The motions of this enzyme are exquisitely tuned to kinase function, as allosteric mutations of PKA-C within the substrate recognition sequence affect the enzyme's opening and closing kinetics, influencing its catalytic efficiency [56][57][58], while the binding of inhibitors quenches the motions in the s-ms time scale. [33,35] The SyncDyn map of the apoenzyme shows that the s-ms motions are located into two separate clusters, spanning residues 5 -100 and 150 -200 (Fig. 4a and Fig. S2). Additionally, the SyncDyn map displays a few correlations among residues at the C-terminal tail, which wraps around the two lobes, and these two well-defined dynamic clusters.
A closer inspection of the overall and range-specific SyncDyn maps reveals that the apo enzyme displays residues spanning from the small to the large lobe, i.e., V57, L82, I94, L95, L103, V104, I150, and L172 (Figs. S2a-b). Among those residues, V57 is part of the kinase's C spine, a continuous array of hydrophobic residues necessary for its activation. [37,38] L103 and V104 bridge the C spine with the regulatory (R) spine, formed upon phosphorylation of the activation segment. I150 is another essential residue within the core of the enzyme, and when mutated into alanine, PKA-C does not undergo autophosphorylation at T197 and S338, rendering the enzyme inactive. [31]SyncDynSyncDyn
Upon nucleotide binding, the density of the R-squared correlations from residue 50 to 200 increases, which is consistent with the structural and dynamic role of the nucleotide acting as an allosteric effector. [37,38] Notably, the distribution of the exchange rates changes from an average of 800 to 2000 s --1 , suggesting that the ATPN binding shifts the frequency of the opening and closing of the enzyme's active site (Figs. 4a-b). Of note, this change in motion corresponds to an increase of binding affinity between the kinase substrates or pseudo-substrate inhibitor PKI. [33] The SyncDyn map, along with the residues mapped onto the X-ray structures (Figs. completely disrupts kinase's activity. [31] Motions in the same timescale is identified for residues that bridge the E-helix in the C-lobe and the C-helix in the N-lobe, including I94, L160, L162, and L157. The network of synchronous residues crosses the enzyme's structure horizontally, at the hydrophobic interface between the two lobes harboring the substrate binding pocket and radiates to the small lobe via the hydrophobic spines.(Fig. S5) Residues linking the B and C helices also fluctuates in a similar time scale, which agrees with the critical role of these segments in assembling the kinase's active conformation. [32] We then tested the two inhibited forms of the kinase. Both inhibitors quenched the internal dynamics in the s-ms time scale. [39] For both ligated forms of PKA-C, we observe a partial (balanol) and a complete (H89) disruption of the dynamic networks, showing that the change in frequency of the enzyme's fluctuations may affect its affinity for substrate binding (Fig. 4). The populations of the conformationally excited states as determined from the fitting of the CPMG relaxation dispersion curves are summarized for all the kinase forms in Fig. S4 and reflect the dynamic activation of the nucleotide and the suppression of the conformationally excited states by balanol, [39] a phenomenon previously observed in the catabolite activator protein. [40] Note that out of the 102 methyl probes assigned and analyzed with SyncDyn and using the GUARDD software, [41] eleven residues could not be fitted using a two-state model and are consequently omitted from the SyncDyn plots These residues are distributed throughout the enzyme and are illustrated in Fig. S6. Additional experiments at various field strengths will be necessary to accurately quantify the kinetics and populations of these residues, a task that lies beyond the scope of this study.
Overall, the SyncDyn maps highlighted the clusters of residues that fluctuate at the same frequency both in the free and ligated forms. This information is complementary to the CHESCA maps, [22] and provide a starting point to interrogate the allosteric phenomena occurring within the kinase upon ligand binding and mutations. [42][43][44] Analysis of the slow motions in the T17A mutant of RNase A. As a second example, we applied our analysis method to the 15
A. This well-studied enzyme catalyzes the cleavage of single-stranded RNA without metals or cofactors. The RNase A architecture consists of three -helices and seven -strands organized into two lobes with the active site positioned between the lobes (Fig. S1b). Motions play a significant role in the catalytic cycle of RNase A as the rate-determining step (product release) is coupled with a conformational change involving residues throughout the entire protein structure. [45][46][47] NMR relaxation dispersion data revealed that the residues of the RNase can be grouped into two dynamic clusters. [48][49][50] Specifically, loops 1 and 4 positioned at ~20 Å away from the active site impart substrate fidelity, recognizing purine residues at the 5' end of the phosphoester cleavage site (Fig. S1b). These loops undergo a conformational change of 2-3 Å upon ligand binding, defining the closed state of the enzyme. After the chemical step, RNaseA requires an additional conformational change on the millisecond timescale for product release. Mutagenesis and structural dynamic studies have shown that the H48 side chain forms critical H-bond interactions with the surrounding residues affecting the motions of remote structural elements, such as loops 1, helices 1 and 2, as well as 1, 2, and 4 strands (Fig. S1b). In particular, the H48A mutation disrupts these H-bond interactions dramatically affecting the enzyme's activity. In contrast, mutations at T17 (also involved in H-bond interactions) slightly accelerate the product off rate by increasing the motion of loop 1. [49] Because of the extensive data available for this mutant, we chose it for our SyncDyn analysis. The SyncDyn map generated with these data and the distribution of the kex values are reported in Figs. 5a-b and Fig. S7. The SyncDyn map shows that several residues intersperse throughout the protein experience a motion in the same timescale, i.e., residues located around loop 4 connecting strands 2 and 3 (A64, Q69, T70, and N71). These residues fluctuate in the same timescale as those belonging to the adjacent -strand and distal residues located in the -hairpin connecting strands 4 and 5 and helices 1 and 2. Finally, C95
fluctuates with a kex similar to residues located in the opposite dynamic cluster, i.e., I106, A109, E111, H119, and S123. Note that the comprehensive SyncDyn map contains the information for a broad range of kex, spanning from 0 to 5000 s -1 . However, the distribution of the kex values indicates that structural fluctuations with R-squared of 0.85 occur between 0 and 5000 s -1 (Fig. 5b), with most of the fluctuations occurring around 3000 s -1 , while only a few residues fluctuate at frequencies between 0 -1000 and 4000 -5000 s -1 range. Overall, the SyncDyn map for
RNAse A shows that many residues belonging to the two dynamic clusters move synchronously.
However, it should be noted that the motions between the two lobes of the enzyme are not coupled, although the fluctuation of the individual residues occurs in the same time scale [50].
Synchronous motions have been suggested to favor allosteric transmission from the ligandbinding site to remote sites of proteins. [51] Recent studies also infer that these structural fluctuations must be correlated and occur within a defined frequency range to modulate allosteric phenomena efficiently. [52,53] NMR studies of various enzymes have shown that V-type and K-type allosteric enzymes' regulation occurs in the s-ms time scale, [1,2,[54][55][56][57] a window of motions exquisitely captured by the CPMG dispersion experiments. [3,[58][59][60] The dispersion data obtained from CPMG experiments at multiple frequencies are used to extract kinetic and thermodynamic parameters on the conformational equilibrium of biomacromolecules. The most common analysis method involves a global fitting of the flexible residues using the Bloch-McConnell equations under the assumption of a two-site exchange equilibrium. [23] When a global kex is obtained, subsets of residues are selected to identify clusters of residues that fluctuate at a similar frequency and may be involved in allosteric networks. The latter comes at the expense of the accuracy of individual curve fitting. [61] In contrast, our method utilizes a geometrical approach to evaluate pairwise synchronicity by measuring the distance between residues in the kf/kr space. If residues share similar kinetic parameters, they cluster in the kf/kr space, displaying correlations in a 2D SyncDyn map. Although computationally demanding, this method does not rely on an arbitrary selection of residues as it systematically analyzes all residue pairs, and it objectively provides a more accurate fit of individual dispersion curves than the standard global fitting approach.
As expected, we found that pairwise dynamic correlations critically depend on the quality of the relaxation dispersion data. Specifically, spectral noise levels may cause false positives (or negatives) that affect the fidelity of the SyncDyn maps. To avoid this, the new algorithm utilizes pair-fitting procedures stricter than individual fitting, and the goodness of fit (R -squared) avoids the overestimation of pairwise correlations. Note that our method only clusters residues that fluctuate at similar frequencies (kf and kr) and does not identify whether these fluctuations are coupled or correlated (e.g., RNAse A). Indeed correlated motions have been proposed to be a critical property of allosteric proteins [62]; however, determining whether motions are correlated requires the support of computational approaches, [63] other spectroscopic methods [62], as well as experimental validations via mutagenesis. Finally, it should be noted that the dynamic networks of residues reported by the SyncDyn maps are underestimated as it depends on the number of NMR probes available and the window of motion analyzed, i.e., residues manifesting relaxation dispersion. A possible solution is to expand our analysis to other range of motion using T1 measurements [3] or adiabatic relaxation dispersion experiments. [64,65] Enzyme function involves both structural and dynamical changes in the biomacromolecular matrix. The analysis of ligand titrations on protein fingerprints has been instrumental in characterizing orthosteric and allosteric effects in proteins. [66] Also, chemical shifts of both backbone amide and methyl-group side chains can be used to describe conformational shifts upon ligand binding [67] and identify networks of distally connected residues. [68][69][70][71][72][73]. The SyncDyn maps complement these approaches to fully characterize the structural and dynamic response to ligand binding. The synchronicity of motions adds another layer of complexity to global and sub-global motions, though it is not universally applicable to all systems. However, as we gain access to dynamic data on larger biomacromolecules, it will be intriguing to determine whether synchronous motions play a role in biological activity and whether they are a common characteristic of cooperative systems. Additionally, this approach could be useful for exploring the organization of long-range dynamics in intrinsically disordered proteins [74] and polypeptides that exhibit multiple conformational states. and polypeptides with multiple conformational states. [75][76][77] In conclusion, we present a method to identify networks of residues that fluctuate with similar frequencies (synchronous) that may mediate allostery in certain proteins. Our method largely removes experimenter bias at the early stages of relaxation dispersion analysis and will allow a standardized protocol for assessment of important protein motions. This approach provides 2D correlation maps that, when integrated with chemical shift covariance (CHESCA) maps for fast exchanging systems, offer structural and dynamic canvas as a starting point to interrogate allosteric phenomena in proteins and to guide future enzyme engineering or drug discovery efforts.
Using Matlab®, we generated 15 N CPMG relaxation dispersion data for a synthetic twenty-aminoacid peptide, using the Bloch McConnell equations. We created two clusters each consisting of six residues that share the same kinetic parameters. Cluster 1 with kf = 25 s -1 , kr = 500 s -1 , Δ = 1 ppm, and R20 = 10 s -1 ; and Cluster 2 with kf = 100 s -1 , kr = 2000 s -1 , Δ = 1 ppm, and R20 = 10 s -1 . The data were generated for two magnetic fields with 1 H Larmor frequencies of 600 and 800
MHz. We set the TCPMG at 20 ms and generated synthetic data for CPMG = 50, 100, 150, 200, 250, 300, 350, 400, 500, 600, 700, 800, 900, and 1000 s -1 in which CPMG is the rate at which pulses are applied during the constant time relaxation, TCPMG. Finally, we included 5% random noise into each dispersion profile. For the remaining residues, we randomly assigned values for kf and kr within a range of 0 to 5000 s -1 .
The 13 C-ILV labeled catalytic subunit of protein kinase A (PKA-C) was expressed and purified as previously described. [22] Methyl-TROSY relaxation dispersion data, previously published, [22] were acquired on a 230 µM PKA-C sample at two magnetic fields (16.4 T and 19.97 T) using a constant relaxation time of 40 ms and inter-pulse delays frequencies CPMG = 50 (× 2), 100, 150, 200, 250, 300, 400, 500 (× 2), 600, 800,1000 (× 2) s -1 . All the experiments were carried out at 300 K. The spectra were acquired in interleaved mode as a pseudo-3D with two-dimensional planes of 80 x 1024 complex points. The peak intensities were extracted using the FuDA (https://www.ucl.ac.uk/hansen-lab/fuda/).
For the RNase A, experimental CPMG relaxation dispersion data, published previously, [50] were acquired at static magnetic fields of 11.7, 14.1, and 18.8 T on a 773 M 15 N-labeled T17A
RNAse A mutant and a 600 M 15 N-labeled wild-type (WT) RNAse A sample at pH = 6.4 and 298 K expressed and purified as previously described. [47,49] The NMR relaxation data were typically acquired with spectral widths in the t1 and t2 dimensions of 2700 and 10000 Hz with a digital resolution of 0.1 and 0.2 points/Hz respectively. The proton carrier frequency was set to the H2O resonance, and the 15 N frequency was placed in the middle of the amide region at 120 ppm. Spinrelaxation rate constants at each CPMG pulse repetition time (cp) were acquired using a constant time [78] total relaxation delay of 40 ms. Transverse relaxation rates were determined for CPMG frequencies CPMG = 50, 100, 150, 200 (× 2), 250, 300, 400, 500, 700 (× 2), 800 s -1 . The individual spectra were acquired as an interleaved three-dimensional experiment in which the two-dimensional planes were extracted, and peak intensities were determined from a 3 x 3 grid using inhouse written software. from 0 to 5000 s -1 and (∆𝑘 𝑒𝑥 = 𝑘 𝑒𝑥 𝑋 -𝑘 𝑒𝑥 𝑌 ranging from 0 to 3000 s -1 ). The response was simulated with noise amplitudes of 2% (a) and 5% (b) of the intrinsic relaxation rate (𝑅 2 0 = 10 s -1 ).