Role of Terahertz (THz) Fluctuations in the Allosteric Properties of the PDZ Domains.

With the aim of investigating the relationship between the fast fluctuations of proteins and their allosteric behavior, we perform molecular dynamics simulations of two model PDZ domains with differential allosteric responses. We focus on protein dynamics in the THz regime (0.1-3 THz) as opposed to lower frequencies. By characterizing the dynamic modulation of the protein backbone induced by ligand binding in terms of single residue and pairwise distance fluctuations, we identify a response nucleus modulated by the ligand that is visible only at THz frequencies. The residues of this nucleus undergo a significant stiffening and an increase in mutual coordination upon binding. Additionally, we find that the dynamic modulation is significantly more intense for the side chains, where it is also redistributed to distal regions not immediately in contact with the ligand allowing us to better define the response nucleus at THz frequencies. The overlap between the known allosterically responding residues of the investigated PDZ domains and the modulated region highlighted here suggests that fast THz dynamics could play a role in allosteric mechanisms.

Introduction

Allostery is defined as a process in which an event at one site of a protein impacts the structure and/or the dynamics of another site, regulating protein functions primarily via the modulation of affinity or binding modes for ligands. Such a process involves a long-range communication between the two sites that can be mediated by changes both in structure and in dynamics. While one can measure the allosteric coupling between two binding sites in terms of affinity change, i.e., global free energy shift induced for the binding event in the presence of the allosteric effector,[1] the events underlying the coupling effect typically involve a group of amino acids that define the responsive subset and can be investigated through a number of approaches. In the presence of structural changes, high resolution crystal structures of the apo and holo complexes can give a mechanistic insight into allostery. However, besides structural changes, the rearrangements that underlie allosteric functional regulation can also include dynamic modulation, which can be indirectly probed by mutagenesis,[2] sequence-based bioinformatics methods such as statistical coevolution analysis (SCA),[3,4] and computational approaches based on Gaussian network models[5] or molecular dynamics.[6] In particular, the role of protein dynamics in allostery has progressively gained attention since the seminal paper of Cooper and Dryden:[7,8] Over recent years, it has been widely investigated, with studies focusing on motions on different time scales from picosecond-to-nanosecond (ps–ns)[9,10] to microsecond-to-millisecond (μs–ms),[9,11−15] where the former are mainly dictated by backbone and side-chain fluctuations and the latter involve large-scale conformational changes. In this context, the potential relation of intramolecular vibrations in the terahertz (THz) frequency range to allosteric effects has just started to be investigated by the new emerging THz spectroscopy techniques.[16,17] These vibrations may play a key role in allosteric control by distorting the energy landscape upon ligand binding to enhance access to functional configurations.[16] For rhodopsin-like G protein-coupled receptors, a family of model allosteric membrane proteins, it has been proposed that a protein–ligand vibrational resonance in the THz frequency range creates an allosteric association of coupled fluctuations forming a coherent signal pathway from the orthosteric ligand-binding site to the activation region.[17] Recently, the importance of these fast fluctuations has become a topic of intense debate for several processes occurring in biomolecules, e.g., in directing biochemical reactions, in mediating efficient protein–ligand binding, in initiating and modulating slower dynamical processes, and in assisting transmembrane transport of small molecules.[18−20] In the framework of molecular recognition, shedding light on the role of the THz dynamics of proteins could provide significant insights into the fine determinants of dynamic allostery. Finally, allosteric mechanisms could be linked to a hierarchy of time scales in protein dynamics: Indeed, elucidating the relationships between slow and fast dynamics and the role they play in allostery still represents a challenging task.[21]

Here, by means of molecular dynamics (MD) simulation analyses, we investigate the issue of how protein dynamics in the THz frequency regime may contribute to allosteric mechanisms.

We consider two well-known model proteins, for which long-range modulation and allostery have been addressed and extensive dynamic, sequence-based, and mutagenesis data are available: the third PDZ domain of PSD95, PDZ3, and the second PDZ domain from tyrosine phosphatase, PDZ2.[4,22−25]

PDZ domains function as modular protein–protein binding elements, for which evolutionary coupling between binding site and distal regions has been established.[4] It has been proven that under physiological conditions ligand binding is allosterically modulated for the Par6 member of the PDZ family[26,27] and more generally long-range modulation of ligand affinity has been observed, both through mutations and through post-translational modifications.[28−30] As such, PDZ domains constitute a popular reference system to test computational approaches for detecting allostery in single domain proteins.[31−33]

In particular, evidence indicates that in the PDZ3 and PDZ2 protein domains distal regions, which are evolutionarily coupled to the binding site, dynamically respond to the presence of the peptide ligand[4,23−25] and the modulated region extends to distal regions such as helix H1 and the opposite side of the binding site (see Figures and 2 for reference to the nomenclature of secondary structure elements). Overall, these systems appear suited for investigating the role of the THz dynamics in proteins with modulated allosteric responses.

Figure 1

PDZ3. Slow dynamics (t > 10 ps; E < 0.1 THz): (A) RMSF per residue, unbound (bound) state in black (orange), (B) protein structure with highlighted residues showing statistical significant increased (red) and decreased (blue) mobility upon binding. The binding site is represented by means of a wireframe gray surface. The secondary structure elements discussed in the text, loop L12, loop L23, and helices H1 and H2, are highlighted with labels using the nomenclature of ref (31) and correspond to the following residues: L12 res 318–322; L23 res 330–335; H1 res 345–351, H2 res 371–383. Fast THz dynamics (0.3–10 ps; 0.1–3 THz): (C) same as A and (D) same as B.

Figure 2

PDZ2. Details are the same as in Figure . The secondary structure elements discussed in the text: loop L12, loop L23, and helices H1 and H2 are highlighted with labels using the nomenclature of ref (31) and correspond to the following residues: L12 res 13–19; L23 res 24–34; H1 res 45–50; H2 res 70–79.

In this study, by comparing the variations upon ligand binding of THz dynamics (0.3–10 ps; 0.1–3 THz; ∼3–100 cm–1) with those at lower frequencies (t > 10 ps; E < 0.1 THz), we aim to clarify if and how patterns of coordination among residues that may be linked to allosteric couplings can be highlighted by filtering the frequency domain.

We first set out to evaluate whether a dynamic modulation upon ligand binding is observed at THz frequencies and if it differs from that at lower frequencies. We next address the question regarding whether the observed dynamic modulation can be related to the protein allosteric behavior. We identify a subset of modulated residues, named the response nucleus for binding (RNB), which we characterize dynamically in terms of single residue fluctuations, correlations, and pair rigidity, focusing on the THz frequencies as opposed to the lower frequency regime, and analyze it in comparison to the known allostery of PDZ domains[22−24] to highlight specific dynamic signatures of allosteric modulation. In particular, for a better characterization of the RNB we investigate the dynamics both of the backbone and of the side chains: due to their heterogeneity and high sensitivity to perturbation, the latter could reveal fast functional dynamics along the allosteric pathway that are invisible to the former.[34]

Methods

Constructs and Molecular Dynamics Trajectories

MD simulations of the PDZ3 domain from the synaptic protein PSD-95 were carried out in its free state (PDB structure: 1BFE) and in complex with the pentapeptide CRIPT (PDB sequence, 1BE9; peptide sequence, KQTSV). The simulations were performed with the GROMACS suite,[35] using the GROMOS G43A1 force field both for the protein and for the peptide ligand[36] and the SPC water model[37] for the solvent. Each system was solvated with ∼6000 explicit water molecules filling an octahedral box. For each of the two investigated states, equilibration was followed by 400 ns simulations at constant pressure and temperature with a time step of 2 fs, as described in detail in ref (31). The atomic positions and velocities, saved every 10 ps, were collected for the analyses of protein dynamics slower than 10 ps (i.e., frequencies <0.1 THz). In the following we will refer to these dynamics as slow dynamics or low frequencies. Additionally, 9 snapshots were selected from the long trajectory of bound and unbound states at given time interval (every 10 ns between 90 and 180 ns) in order to generate a set of distinct equilibrated starting structures in the vicinity of the native conformation: these snapshots were used as starting points for independent 100 ps microcanonical runs with a time step of 2 fs. The atomic positions and velocities, saved every 10 fs, were collected.

In order to specifically investigate protein dynamics in the 0.1–3 THz range, we have applied to these trajectories a Fourier filtering method[18,38] that enables us to analyze motions in this selected frequency window. These filtered trajectories were used for the analyses of the fast fluctuations: in the following we will refer to these dynamics as fast dynamics or THz frequencies.

The Fourier filtering method comprises three main steps:[18] (1) A Fourier transform operation is applied to the trajectory of the system, yielding the amplitude spectrum in the frequency domain. (2) A specific frequency window is selected and the amplitude spectrum is set to zero outside this window; this operation yields a reduced spectrum. (3) The reduced spectrum is transformed back to the time domain, so as to obtain the filtered trajectory. Here, we have implemented Fourier filtering thanks to an efficient program, developed by Turton and coauthors.[18]

In summary, for each of the two investigated states, i.e., PDZ3 in the unbound and the bound state, the main results of our study have been obtained from the following sets of trajectories: (a) a 400 ns trajectory (slow dynamics, low frequencies), (b) 9 unfiltered trajectories of 100 ps length, and (c) 9 filtered trajectories of 100 ps length (fast dynamics, THz frequencies). The quantities calculated on the groups b and c are averages over the full set of trajectories.

The same protocol was adopted for the PDZ2 domain from the tyrosine phosphatase PTP-1E, investigated in its free state (PDB structure: 3LNX) and in complex with the RAGEF2 C-terminal peptide, (sequence: ENEQVSAV) (PDB structure: 3LNY).

Communication Propensity

The analysis of the communication propensity CP between any two residues, introduced in ref (6), is based on evaluating the distance fluctuations between selected atom pairs. The distance fluctuations among residue pairs report on the relative rigidity of their position, hence the residue–residue coordination, and therefore give insight into the internal dynamics sampled along the trajectory.

The calculated quantity is CP, defined for any two residues as CP = ⟨(d – dij,ave)2⟩, where d = |r – r| is the distance between the Cα atoms of residue i and residue j.

Residues whose Cα–Cα distance fluctuates with a relatively small intensity during the trajectory, and hence with low CP values, are considered more efficiently coordinated than residues whose distance fluctuations are large.[6,39] The change of distance fluctuations for distal residue pairs, induced by ligand binding, can report on dynamic allosteric effects caused by the ligand.

Results

RMSF Fluctuations

In order to highlight changes in protein dynamics upon ligand binding in the two investigated frequency windows (set of trajectories a and c), we have first computed the root-mean-square fluctuations (RMSF) of PDZ3 and PDZ2 for the bound and the unbound state. The profiles, reported in Figures and 2, allow us to highlight regions undergoing a statistically significant increase or decrease of flexibility upon binding. The former (latter) are identified by the following criteria: RMSFbound–σbound > RMSFunbound+σunbound (RMSFbound+σbound < RMSFunbound–σunbound) and are represented by means of a red (blue) surface in the protein structure (panels B and C). Here, σ is the standard deviation that has been calculated over the set of independent trajectories c for the fast dynamics. For the long trajectory (set a), the standard deviation is estimated by means of a block average: The trajectory is divided into blocks of 10 ns, and for each block the RMSF per residue is calculated. The average RMSF is represented in the profile together with the standard deviation. σunbound (σbound) refers to the simulations of the free (complex) state.

For both the PDZ3 and PDZ2 systems, our analysis reveals that the frequency range is critical to observe the specific effects of the ligand, e.g., in the binding site (here defined by the residues of the atoms that lie within 3 Å of the peptide and represented by means of a wireframe surface in the protein structure in Figures and 2; for PDZ3, residues 323–325, 327, 339, 372; for the PDZ2, residues 17–20, 22–25, 79).

For the PDZ3 system, several regions undergoing an increase of flexibility upon binding (residues around 320, 345, 365, 375–383, 412–415) show enhanced dynamics in both frequency ranges. In more detail, a statistically significant increase of flexibility upon binding is observed for residues 319, 320, 321, 344, 347, 348, 380 at lower frequencies (Figure A) and for residues 320, 344, 345, 363, 389, 395, 414, 415 at THz frequencies (Figure C): these residues are highlighted by means of red surfaces in the protein structure (Figure B,D). In contrast, for regions showing reduced mobility upon binding (blue surfaces in Figure B–D), there is qualitative agreement in the two frequency ranges only in the C terminal part around residue 405, corresponding to helix H3 and the terminal β sheet. In the THz regime, mobility decrease is observed for regions around residues 327 and 372 (more in detail: res 312, 324, 326, 327, 328, 338, 371, 372, 404, 405, 406). This set includes a variety of residue types, including some residues of the binding site (res Gly324, Ile327, His372) that are in contact with the peptide. Another region with reduced flexibility upon ligand binding, but not directly interacting with the ligand, is located around residue Arg312. In contrast, no significant impact of the ligand on the fluctuations of the binding site is observed at low frequencies, since on time scale longer than 10 ps, the binding site appears rather equally rigid both in the bound and unbound systems.[31]

In the PDZ2 system, analogously to the results found for PDZ3, we find that only in the THz frequency range a set of residues undergoes a statistically significant stiffening response upon binding (res 19, 20, 21, 22, 39, 40, 78). Once again, this set includes some residues of the binding site (res Gly19, Ile20, Val22). A consistent increase of flexibility upon binding is shown by residues located at the loop region of the protein structure (near helix H1 at residue Ala46 at THz frequencies and residues Ser17 and Ser65 at low frequencies).

Overall, we observe that the stiffening response of the binding site region is visible only at THz frequencies. This supports the importance of filtering over multiple frequency ranges, and in particular of carefully considering the THz regime, as a tool to highlight the modulations of protein dynamical responses to a ligand at the local level. On the other hand, for the PDZ3 system, our analysis suggests some overlap between high and low frequency ranges for the regions showing increased mobility. On this basis, we define the residues undergoing a consistent mobility decrease on the THz scale as the ligand responding subset, or response nucleus for binding (RNB). In parallel, the regions increasing their fluctuations, that are located at the periphery of the protein structure and mainly in protein loops, constitute a second responding nucleus.

We have further characterized the THz dynamics of the two responding nuclei defined by means of the RMSF analysis by calculating the vibrational density of states (VDOS) (on the set of trajectories b), as described in detail in the Supporting Information. In proteins the main feature of the VDOS is a broad peak centered at 2–2.5 THz, which has been attributed to bending and stretching motions of hydrogen bonds that connect groups of atoms in the side chains and the backbone, involving delocalized backbone torsional motions and librations of the side chains, and that partly reflects intramolecular collective short-wavelength motions.[19,40,41] Frequency shifts in the VDOS for a subset of protein residues due to an event, e.g., the binding of a ligand, give us precise and quantitative information regarding the change in the number of vibrational modes at a given frequency for the selected atoms. First, we have calculated the VDOS of the whole molecule (Figure S1), and we have found that there is no appreciable difference between the VDOS of the bound and the unbound system: this finding is in agreement with a recent study that suggested that the contribution of specific modes involved in substrate binding might be more important than the net difference between the bound and unbound state VDOS.[42] Hence, we decided to quantitatively estimate the local and delocalized density of vibrational modes involving the atoms belonging to the two response nuclei revealed by the RMSF analysis. Both for PDZ2 and PDZ3, we find that the ligand has an impact on the THz vibrational modes of the RNB, causing a blue-shift of its vibrational frequencies (cf. SI, Figures S2 and S3). This finding, together with the results of the RMSF analysis, points to a stiffening and a reduction in flexibility of the RNB at THz frequencies. Vice versa, the second responding set shows an opposite behavior, i.e., its vibrational frequencies are red-shifted in the presence of the ligand (cf. SI, Figures S2 and S3). Hence, the two responding sets highlighted here could be involved in specific modes contributing to substrate binding.

Overall, the result shows that although the total VDOS or its side chains and backbone contributions[18] may show no appreciable difference between the bound and unbound state, some residues of the proteins, e.g., as in the case of the two responding sets highlighted here, do indeed show changes in some vibrational modes upon ligand binding.

Binding Induced Modulation of the Internal Dynamics

Prompted by the dynamic modulation of single residues evidenced in the RMSF profile, we set out to evaluate the dynamic coupling among residues by measuring their rigid coordination. The aim of this analysis is to detect the impact of ligand binding on the overall structural dynamics of the protein and particularly to elucidate the connection between the enhanced rigidity of the binding site observed at THz frequencies (Figures C and 2C) and the modulation of flexibility throughout the structure, also in regions distal from the binding site, focusing both on the backbone and on the side chains.

In order to elucidate the mechanisms of signal propagation from the binding site and its possible dependence on the frequency range, we applied the method of communication propensity (CP) analysis to our simulations.[6] This method is based on the hypothesis that signal transduction events in proteins can be described in terms of the modulation of rigidly coordinated subunits that affect the overall structure and dynamics (see Methods section for details).

We performed the CP analysis for the Cα atoms of PDZ3 and PDZ2 in the bound and unbound state in the two investigated frequency ranges (set of trajectories a and c). We calculated the difference between the CP of the two states (bound–unbound), hereafter referred to as deltaCP matrix: in this matrix blue (red) stripes describe residue pairs that are more (less) coordinated following the binding event.

The results of the analysis on the Cα atoms of PDZ3 and PDZ2 are shown in Figures A and 4A (frequencies < THz) and 3D and 4D (THz frequencies).

Figure 3

PDZ3. Results of the communication propensity (CP) analysis for the slow dynamics (upper panel) and fast THz dynamics (lower panel). The difference between the CP of the two states (bound–unbound) for the Cα (first column), Cβ (second column), and methyl C atoms (third column) is shown. Positive (negative) values indicate increased (decreased) distance fluctuations among residue pairs, corresponding to a decreased (increased) coordination upon binding. The axis labels for panels B, C, E, and F refer to the numbering given in the Supporting Information.

Figure 4

PDZ2. Caption details as in Figure .

From a general point of view, it is readily seen that, as in the case of RMSF (Figures and 2), the ligand effect on the pairwise distance fluctuations shows a dependence on the frequency range (cf. upper and lower panels in Figures and 4).

For the PDZ3 system, in the low frequency regime (Figure A) the regions showing increased mobility upon binding, highlighted in Figure B, determine a perturbation in the pairwise coordination that is visible as red stripes throughout the matrix (i.e., a decreased coordination). Loop L12 (res 320) and helix H1 (res 345) markedly move with respect to the protein bulk. Also, H2 (res 365) moves relative to the protein. On the other hand, stiffening of the chain at the C terminal end (res 405) is reflected in higher coordination relative to the rest of the molecule, and visible as set of parallel blue stripes. The rest of the binding site region however does not emerge in this representation.

In contrast, THz frequencies (Figure D) are characterized by a pattern of locally increased coordination, in agreement with the RMSF analysis, surrounding the binding site (around res 327 and res 372) (Figure D) and generating a rigidly coordinated unit, as shown by the blue stripes in the delta CP matrix. A higher coordination is observed also here for the C terminal end. On the other hand and similarly to lower frequencies, increased distance fluctuations are observed both for loop L12 (res 320) and for helix H1 (res 345) relative to the rest of the protein (yellow-red stripes). A spot of local mobility increase is observed in the neighborhood of residue 365, just before helix H2, reminiscent of the motion of this region observed on longer time scales (red stripes in the interval 365–370 in Figure A). Overall, we can conclude that at THz frequencies the regions with increased (decreased) pairwise coordination are positioned as those with decreased (increased) mobility revealed by the single residue RMSF analysis. In agreement with the RMSF analysis, we observe that the increased mutual coordination of the binding site in the bound PDZ3 system with respect to the unbound system is revealed only at THz frequencies.

Notably, Fuentes and coauthors have shown that in PDZ2 the most significant dynamical changes on the ps–ns time scale upon binding are observed at the side-chain level, suggesting that allosteric behavior can be detected using side-chain methyl dynamics measurements on this time scale.[23] In order to highlight a possible differential modulation of main chain and side-chain response on the two time scales considered here, the CP analysis was then repeated on two selected sets of side-chain atoms of PDZ3, namely, on Cβ atoms and on the carbon atoms of the methyl groups. For residues containing multiple methyl groups, all the methyl groups have been taken into account in the analysis. Tables S1–S4 in the Supporting Information show the progressive numbering of the Cβ and methyl groups, associated with that of the corresponding amino acids.

The low frequency regime (Figure B) is characterized by a rather similar modulation pattern for Cα and Cβ atoms, with the exception of an increased rigidity observed for single residues (Cβ 20 and 37, corresponding to residues 328 and 352, cf. Table S1 in the SI). For residue 328, which is located in the vicinity of the binding site, the observed rigidity is due to side-chain interactions with the ligand. The same trend is observed for the methyl C atoms (Figure C). Therefore, on time scales slower than 10 ps, we conclude that the backbone and the side chains seem to have consistent relative motions.

On the other hand we find that in the THz regime (Figure E for Cβ and Figure F for methyl carbon atoms) the intensity of the modulation, and specifically the stiffening effect determined by the ligand and propagated through the response nucleus, appears more intense in the side chains than in the backbone, and particularly in the methyl C atoms. Rigidified methyl carbon atoms, identified by selecting the most negative contributions to the CP difference matrix, are localized in the binding site (327) but also in core regions away from it (353, 388), before and after helix H1 (342, 343, 349), and in the terminal β strand (406, 412) (cf. Table S2 in SI → progressive # in Figure F are 12; 29, 54; 22, 24, 26, 27; 60, 62). The residues containing these rigidified methyl carbon atoms are mapped on the protein structure in Figure (left panel). In contrast, decreased coordination is observed locally in the N terminal end (residue 307), in the segment preceding helix H2 (360, 362, 365; progressive # 32–37) and at the C terminal part of helix H2 (residues 377 379; progressive # 44, 46, 47).

Figure 5

PDZ3 (left) and PDZ2 (right) protein structures with highlighted residues (cyan) containing methyl carbon atoms with increased coordination upon binding in the THz regime, as revealed by the CP method (see text for details). Bottom panel: amino acid sequence alignment of PDZ3 (PDB structure: 1BFE) and PDZ2 (PDB structure: 3LNX). The residues of the RNB are colored in blue (identified by means of the RMSF analysis) and cyan (identified by means of the CP analysis). The residue of the second responding nucleus (i.e., undergoing an increased mobility upon binding at THz frequencies as determined by the RMSF analysis) is colored in red.

Focusing on the CP analysis for the Cα atoms of the PDZ2 system, we notice, in agreement with the PDZ3 case, that the coordination pattern is modulated upon binding in a frequency-dependent manner.

A clear difference between PDZ2 and PDZ3 emerges at low frequencies (cf. Figures and 4, upper panels): the changes of the distance fluctuations due to the binding are mainly due to increased pairwise coordinations in the PDZ2 system (blue stripes) and decreased pairwise coordinations in the PDZ3 system (red stripes). The former result from structural changes or motions occurring on the nanosecond scale, such the previously observed local unfolding of H1 or motion of loop L23 that is damped by the ligand,[31] is not observed in PDZ3. In more detail, for PDZ2 in the low frequency regime (see Figure A), some residues belonging to the binding site region (around res 20), to L23 (28–31), to H1 (45–48), and to H2 (70–80) increase their mutual coordination (blue stripes along the matrix) upon binding, while nearby those around res 25 show in contrast a coordination decrease (red stripe). The pattern does not change when going from Cα to the side-chain atoms, which again suggests that this analysis mainly captures residue motions that are affected by the global conformational dynamics. In contrast, when focusing on the THz frequencies, the analysis reveals a modulation between the main chain and side chains. Similarly to the PDZ3 case, here a stronger coordination increase in the bound system with respect to the unbound one emerges for the methyl carbon atoms in the THz regime, involving an extended region: the binding site (specifically, residues 20–22–26) and the nearby loop L23 (28–30) and stability core (35) as well as distal residues near helix H1, 39–40–41, and 66, 87–89. These rigidified methyl carbon atoms are shown in the right panel of Figure (cf. Table S4 in SI: progressive # 5, 8, 9, 12, 15, 17, 18, 19, 24, 28, 42, 43, 50, 61, 62 of Figure F, corresponding to residues 11, 18, 20, 22, 26, 28, 30, 39, 40, 41, 66, 78, 89). Coordination decrease at the level of the side chains is only observed at isolated positions, such as the N terminal β strand containing residues 6, 11, 12, and around residue 61 before helix H2.

In summary, our analysis of the pair rigidity of residues based on the distance fluctuations between Cα, Cβ, and methyl C atoms reveals a consistent increase of the binding site coordination in the presence of the ligand at THz frequencies for PDZ2 and PDZ3, while at lower frequencies the coordination response differs for the two systems, as a consequence of motions occurring on longer time scales.

Moreover, the modulation induced by ligand binding at THz frequencies is significantly more intense for the side-chain atoms, where it is also redistributed to distal regions not immediately in contact with the ligands. This allows us to better characterize the RNB.

The increase of coordination involving an extended responding region, which includes not only the binding site residues but also core and distal residues, is observed in both proteins, and is specifically highlighted by focusing on the methyl carbon atoms, which provide a stronger modulation on the considered time scales.

Discussion

In this paper we have investigated whether a hierarchy of motions on different time scales can be related to the response to binding events and allosteric behavior of model systems. Our findings identify subsets of residues with similar dynamics on different regimes considered, which can be directly reconnected to a number of experimental results discussing allostery in PDZ domains.

According to the SCA approach,[22] while most residues in proteins evolve rather independently from each other, about 20% of amino acids form spatially contiguous networks of amino acids, named sectors, whose evolution is pairwise correlated. In the PDZ family of proteins the sector has been suggested to be the functional unit responsible for the protein allosteric behavior. Along the same line, from the dynamical point of view, NMR relaxation measurements have identified a subset of residues undergoing dynamic modulation upon binding.[23,24] Therefore, allostery in PDZ2 and PDZ3 domains has been mainly defined as dynamic modulation of a subgroup of coevolving residues identified by SCA.

In the PDZ3 system the sector of coevolving amino acids identified in the PDZ family by Lockless and Ranganathan corresponds to the following residues: 322, 325, 329, 330, 347, 353, 372, 376, 380.[4] A more recent study by the same group focusing on the impact of mutations in modulating the ligand-binding affinity[22] identified a region that comprises the sector and significantly overlaps with the dynamic response nucleus we highlight through our THz analysis, in particular with the subset of rigidified residues defined above. The RNB defined here is highlighted in the protein sequence (Figure ). It is worth noting that in PDZ3 the terminal helix α3 (residues 393–409), which we include in the RNB, has been found to be essential in modulating the affinity for the ligand properties.[25,43] On the other hand, the NMR relaxation pattern measured for methyl groups and shown by Petit et al.[25] highlights the increased mobility of helix H1 (residues 344–345), which we observe in the RMSF analysis.

Considering the PDZ2 system, the SCA analysis indicates residues 17, 20, 24, 25, 34, 46, 52, 61, 75, 79, 85 as a coevolving sector. With a different approach, Fuentes and coauthors[23,24] have probed the allosterically responsive region of PDZ2 using side-chain methyl dynamics NMR relaxation measurements, highlighting residues 18, 20, 22, 26, 30, 39, 40, 61, 64, 66, 69, 81, 85. Therefore, in particular, two distal regions were identified that show modulation of their ps–ns dynamics upon ligand binding, without being directly in contact with the peptide. As found for PDZ3, and also for PDZ2, this region overlaps well with the RNB defined here by our THz analysis, obtained by considering the RMSF analysis and the coordination analysis of CH3 (cf. Figure ). Moreover, other regions display increased flexibility on the THz range, confirming what observed in PZD3: these are the area around residue 61 in the loop preceding helix H2, and helix H1 (res 45, 46). The increased flexibility of helix H1 in PDZ2 and PDZ3 might be connected to the known functional properties of the PDZ region involving helix H1 and the nearby loop, as a site for protein–protein interactions that have been hypothesized to allosterically regulate the ligand affinity at the orthosteric binding site.[44,45] Overall, the analysis in the THz range proves effective and computationally efficient in highlighting the common features of the local dynamic modulation in two PDZ homologues, and connect it to the current knowledge about allostery in PDZ domains. Moreover, given the change in the vibrational modes of the response nucleus following the binding event (cf. Figures S2 and S3) and the good overlap between this nucleus and the allosteric region, it is tempting to speculate that these specific THz frequency vibrational modes involved in substrate binding might play a role in the allosteric mechanism.

Conclusions

The methodology adopted in this paper, i.e., the analysis of Fourier-filtered trajectories for the selected range of THz frequencies, has been found to be particularly advantageous to investigate fast protein motions. To the best of our knowledge, here it has been applied for the first time to the study of dynamic allostery. Computational studies of protein dynamics and dynamical allostery have often focused on slow global protein motions. Therefore, faster dynamics are only described by higher order components, mainly as constrained atomic motions. With a complementary approach, here we have focused on fast protein fluctuations, with the aim of reconnecting these motions to the response to a bound ligand. The good agreement between the dynamic modulation observed here in the THz regime and the known allostery in PDZ domains suggests that fast THz dynamics could play a role in protein allostery and long-range communication within a protein and should be taken into account for a proper description of these mechanisms. Our findings suggest that THz dynamics is part of the hierarchy of time scales in protein dynamics that are linked to dynamical allostery.[21,34] Importantly, we have shown that an increase of coordination in the presence of the ligand is more visible at THz frequencies than at lower frequencies and that in the THz regime this modulation is particularly intense for the side-chain carbon atoms of both the binding site residues and distal residues, thus allowing us to better define the response nucleus for binding and localize the allosteric region. This result supports the view that the characterization of side-chains dynamics can provide additional insight about the functional dynamics, which would not be inferable solely from the motions of the protein backbone.[34] In the case of PDZ2 and PDZ3, inspecting THz-scale motions allows us to highlight common dynamic responses, in agreement with the current knowledge of PDZ allostery, which however can coexist with variable dynamic properties on longer time scales, as well as different conformational changes.

We can generally conclude that the characterization of THz dynamics by means of the proposed approach might provide a robust basis for the interpretation of the experimental results obtained by means of the emerging THz spectroscopy techniques.