PDZ (PSD-95/Dlg1/ZO-1) domains are highly abundant protein-protein interaction domains involved in regulating signaling pathways. [1][2][3][4][5][6] They play a critical role in many biological processes, such as managing cell polarity, regulating tissue growth and development, trafficking of membrane protein receptors and ion channels, and regulating cellular pathways. [7][8][9] So far, 268 PDZ domains have been identified in 151 unique human proteins. 10 Despite the broad function and relatively low sequence identity within PDZ domains, the secondary structure is highly conserved. The canonical PDZ domains contain six b-strands and two a-helices and have a single binding site in the hydrophobic groove between the aB helix and the bB strand, 11 as shown in Figure 1A. PDZ domains most commonly interact with the final three to five C-terminal residues of target proteins via the carboxylate binding loop that is defined by the conserved c-f-Gly-f motif, where c is any residue and f is any hydrophobic residue. 12 Various groups have revealed how these highly conserved protein-protein interactions propagate allosteric effects through the PDZ domain. 13,14,23- 27,15-22
The PDZ domain is considered to be a model system to study allostery within small modular domains. Allostery in the PDZ family was initially brought to the table when Lockless and Ranganathan 13 proposed a method to statistically predict allosteric residue networks using multiple sequence alignment. This method is based on networks of energetically coupled residues that are responsible for the propagation of allostery throughout the PDZ domain. This original work sparked a wide interest in studying allostery within the PDZ family. Many efforts have followed Lockless and Ranganathan's footsteps by applying various computational techniques, including direct coupling analysis, 28,29 deep coupling scan, 30 anisotropic thermal diffusion, 31,32 rigid-residue scan, 33 and interaction correlation via molecular dynamics simulations, 21,22,34 to reveal allosteric networks within the PDZ family. Furthermore, experimental groups have expanded our understanding of allostery in the PDZ family with applications of nuclear magnetic resonance (NMR) 16,35,36 and mutational analyses. 28,29,37 Despite the abundance of domains in the PDZ family, these efforts have primary focused on a few well-studied PDZ domains, including Par-6 PDZ, 38- 40 PSD-95 PDZ3, 17,20,23,36,37,[41][42][43][44][45] PTP-1E PDZ2, 16,20,21,35,41,46 PTP-BL PDZ. 17,47 To the best of our knowledge, little attention has yet been given to explore allostery of the PDZ domain in Protein Interacting with C Kinase-1 (PICK1).
PICK1 is a scaffolding protein involved in regulating the trafficking of various membrane proteins via endocytosis. [48][49][50] PICK1 is an especially unique PDZ protein as it is the only protein in the human proteome that is comprised of both a PDZ domain and a BAR (Bin/amphiphysin/Rvs) domain. [51][52][53] The PICK1 PDZ domain forms protein-protein interactions with a variety of integral membrane proteins, including the Dopamine Transporter (DAT) 54 and the GluR2 subunit of the AMPA receptor. 48 Widely accepted hypotheses suspect that such PDZ-protein interactions lead to a propagation of signals through PICK1 that alters its interdomain dynamics. 49,50 This global transduction of signal through PICK1 could be explained by allostery at the PICK1 PDZ domain. The presence of allostery at the PICK1 PDZ domain would have major implications in our understanding of the biological function of PICK1.
The purpose of this study is to use all-atom molecular dynamics (MD) simulations to reveal how the atomic-level interaction pattern affects the interaction mechanisms and dynamics between the PICK1 PDZ domain and two representative ligands. These ligands include the final five C-terminal residues of two natural ligands: DAT and AMPAR GluR2. The two systems of interest are shown in Figure 1B-C. Here, we see that both ligands induce dynamic allostery at the aA helix of the PICK1 PDZ domain. Furthermore, our results suggest that different ligands may trigger different dynamic changes to the PICK1 PDZ domain. Lastly, our work identifies that the hydrophobic core that is formed between the ligands and residue I35 may be key to inducing such dynamic allostery.
We studied two PICK1 PDZ systems: PICK1 PDZ-DAT complex and PICK1 PDZ-GluR2 complex. The DAT ligand refers to the final five C-terminal residues (HWLKV) of the Dopamine Transporter (DAT), and the GluR2 ligand refers to the final five C-terminal residues (ESVKI) of the carboxyl tail peptide of the AMPA receptor GluR2 subunit. Experimentally determined crystal structures of the complex systems were used to generate the starting structure for all all-atom molecular dynamics simulations. (PDB ID: 2LUI 55 and 2PKU, 56 respectively). The PDB file of the PICK1 PDZ-DAT complex (PDB ID: 2LUI) was manually edited by trimming terminal residues to ensure an identical sequence to the PICK1 PDZ-GluR2 system. Each starting structure is shown in Figure 1B-C. Each system was prepared using CHARMM-GUI. 57,58 The most recently developed CHARMM36m 59 force field with explicit solvent (TIP3P) was used in each simulation with the Groningen Machine for Chemical Simulations (GROMACS) package, [60][61][62] version 2020.4. Counter ions (Na + or Cl -) were added to neutralize the systems at 293 K. Steepest-descent minimization and 1-ns MD equilibrium simulations were carried out to generate equilibrated starting structures for the MD simulations. All bonds with hydrogen atoms were converted to constraints with the algorithm LINear Constraint Solver (LINCS) 63 . A Nose-Hoover temperature thermostat 64,65 was used in each simulation. The time step was set as 2 fs, and snapshots were taken every 100 ps. Each system was built in a 90 Å ´ 90 Å ´ 90 Å cubic water box. Each system (PICK1 PDZ-DAT and PICK1 PDZ-GluR2) had four replicates at 7 µs per trajectory, a total of 28 µs (4 ´ 7µs) per system.
The PICK1 PDZ-DAT and PICK1 PDZ-GluR2 complex systems had various dissociation events over the four trajectories (Figure S1). It is important to define a boundary that separates the bound states from the unbound states. Because the PICK1 PDZ-ligand complexes were very dynamic, we considered the distance distributions (Figure S2 andS3) of four key binding residue pairs that have been previously identified 55,56 between the PICK1 PDZ domain and the ligands. For the PICK1 PDZ-DAT and PICK1 PDZ-GluR2 complexes, residue pairs I37-L-2 and I37-V-2, respectively, display the clearest distinction on average between the bound state and unbound states. With these state-defining residue pairs, frames were classified bound or unbound. A bound state is defined as a distance less than 5.0 Å between any two atoms in I37 and L-2 for the PICK1 PDZ-DAT complex, and a distance less than 5.0 Å between any two atoms on I37-V-2 for the PICK1 PDZ-GluR2 complex. To test the accuracy of the defined cutoff, cluster analysis was performed over the bound state trajectories to reveal the most probable positions of DAT and GluR2 about the PIKC1 PDZ domain. In this way, we obtained the five most probable clusters of each ligand. Figure 2 shows the PICK1 PDZ domain in gray while the most probable positions of the DAT (A) and GluR2 (B) are shown by unique colors. Our results confirm that the ligands reside in the PICK1 PDZ binding pocket in the defined bound state trajectories.
The DFI metric estimates the resilience of residues within a given protein system. Being a residue specific metric, DFI calculates relative flexibility scores. 66 By incorporating Linear Response Theory (LRT) and Perturbation Response Scanning (PRS), 67 DFI calculates the response of a residue due to a perturbation on another residue normalized by the average response of all residues in the protein. 41 Position specific dynamics profiles are calculated by utilizing residue covariances.
The Hessian matrix, H, contains the second derivative of potentials. Residue covariances are calculated by taking the inverse of the Hessian matrix, H -1 . The Elastic Network Model (ENM) is commonly used to produce the Hessian matrix. However, to include explicit solvent and better estimate residue interactions, residue covariances can be gathered from an MD simulation production trajectory. In this study, we utilized the MD simulations to calculate residue covariances. ∆R is a response vector calculated by multiplying the covariance matrix with the force vector, F and contains the residue responses. The collection of DFI values calculated from this approach is further refined with a percentile ranking to normalize the scores. A residue with a DFI score less than 0.2 is considered a rigid location, while a position with a DFI score higher than 0.8 is considered a flexible residue. Rigid residues have been found to be important in protein stability and function. 68
Utilizing the same elemental principles as described above, the DCI metric captures the dynamic allosteric coupling of pair of residues in a protein. DCI calculates the response of a residue due to a Brownian force applied to another residue in the same system normalized by the average response of the same residue due to perturbations on the rest of the proteins. The magnitude of the response represents the strength of the dynamic allosteric coupling of a site to another residue being perturbed.
A DCI score applied on binding site residues can reveal other residues in the protein that are highly coupled, meaning a binding event or the dynamics of the residue upon binding will be highly affected. Notably, the DCI score is not an indicator of binding dynamics but rather how the binding dynamics are coupled to the rest of the protein. DCI metric can uncover long range allosteric communications related to the binding event. 66,69,70 Residues with a high DCI score indicate strong coupling with binding site and a position with a low DCI score is considered weakly coupled to the binding site.
Network analysis calculates the correlated movements between residues within a protein or protein complex by constructing residue-based and community-based weighted network graphs according to a trajectory. During the calculations, each residue is represented by a node in a network and the links between nodes are the cross-correlation values between these nodes. By using the algorithm developed by McCammon, A. J. and Harvey, S. C., 71 the displacement of the Ca atoms are used to assess the magnitude of all pairwise cross-correlation coefficients. If the correlation value is 1, the fluctuations of two Ca atoms are completely correlated. If the correlation value is -1, the fluctuations of two Ca atoms are completely anticorrelated (same period and opposite phase). Lastly, if the correlation value is 0, the fluctuations of two Ca atoms are not correlated. The analysis uses the calculated cross-correlation coefficients to return a community partition with the highest overall modularity value based on Girvan-Newman style clustering. 72 All the above analysis was carried out using the bio3d package [73][74][75] .
To quantify the degree of local frustration associated with the binding of different ligands to the PICK1 PDZ domain, the Frustratometer server (http://frustratometer.qb.fcen.uba.ar/) [76][77][78] was used to evaluate the two PDZ-ligand complexes investigated here. Default parameters were used when carrying out the assessments of local frustration, e.g., a 5Å radius cutoff value was applied. The PDB structures used in local frustration analysis contained only the PDZ domain, and the ligands have been removed.
Each trajectory experienced ligand dissociation events (Figure S1). These dissociation events present a unique opportunity to explore the switching of dynamic states at the aA helix in realtime. First, we reveal the unique and specific ligand-protein interactions related to the dissociation events by performing hydrogen bond analysis across the two complex systems. Hydrogen bond analysis reveals canonical Class II PDZ-ligand interactions with the carboxylate-binding loop in each system (Figure 3). These results are in good agreement with previous experimental work. 79 Additionally, we performed a statistical analysis to rank the probability of each hydrogen bond forming in the binding pocket (Figure S4-5). The PICK1 PDZ-DAT system has three hydrogen bonds that occur in at least 90% of the bound frames, including I37(N)-L-2(O), L-2(N)-I37(O) and V0(N)-I35(O) (Figure S4). The PICK1 PDZ-GluR2 system has three hydrogen bonds that occur in at least 70% of frames with the ligand bound, including I0(N)-I35(O), I37(N)-V-2(O) and V-2(N)-I37(O) (Figure S5). These three most probable pairs in each system are in agreement with each other. In both systems, the most probable hydrogen bonds occur between (1) I37 and the residue at position P-2 of the ligand and (2) I35 and the residue at position P0 of the ligand. These interactions are much more prevalent than interactions between G34-P0 and I33-P0. While the above analysis reveals the most probable hydrogen bonds within each complex, it is unclear if these interactions are simply essential to the stability of complex formation or, ultimately, if they effect the overall dynamics and subsequent dynamic allostery of the system. To connect the changes in protein-ligand hydrogen bonding interactions (particularly, as related to ligand dissociation) to protein dynamics, we explored the correlation between ligand dissociation and the dynamics of PICK1 PDZ domain by calculating the coupling of various residue-residue distance pairs over the first 3 ms of each trajectory. Five pairs were considered in the coupling calculation: I33-P0, G34-P-1, I35-P-2, S36-P-3 and I37-P-4. The five PICK1 PDZ residues were chosen because they comprise the bB strand which has been identified as a key player in ligand binding by previous work. 11 These pairs were selected to represent the overall interactions between PICK1 PDZ domain and ligand. Figure S6 lists the twenty residue-residue pairs for each system that are most strongly correlated with the distance changes between the five selected pairs. Having the relative highest rank in both systems, we consider the distance between I33 (bB strand) and A58 (aA helix) as directly dependent on the atomic-level interactions between the PICK1 PDZ domain and ligand. Interestingly, the I33-A58 distance can also be used to describe the overall distance between the bB strand and the aA helix. We explore the correlation between the PDZ-ligand interactions and the distance between the bB strand and the aA helix below. . First, we will consider PICK1 PDZ-DAT system, where the dissociation of the DAT is weakly correlated with the dynamics of the aA helix (Figure 4A-B). The distance between I37 of the PICK1 PDZ domain and L-2 of DAT was used to trace the dissociation as defined in the Methods section. At ~2.5 ms, the distance between I37 and L-2 spikes as the ligand dissociates from the binding pocket (Figure 4A, black). This dissociation is confirmed by hydrogen bond and surface area analysis. As DAT dissociates, the number of hydrogen bonds and the surface area between the PICK1 PDZ domain and DAT drops to zero (Figure 4A, blue and purple, respectively). The surface area between the PICK1 PDZ domain and DAT was calculated using solvent-accessible surface area. While the dissociation event does not clearly correlate with the RMSD of the aA helix (Figure 4A, red), it does result in a distinct increase in distance between aA helix and the bB strand (Figure 4A, green).
Next, we will consider the representative dissociation event for the PICK1 PDZ-GluR2 system (Figure 4C-D). As shown in Figure 4C, the dissociation of the GluR2 is directly correlated with the dynamics of the aA helix. The dissociation of GluR2 at ~2.0 µs is confirmed by a sharp distance increase between I37 of the PICK1 PDZ domain and V-2 of GluR2 (Figure 4C, black), a loss of hydrogen bonds between the PICK1 PDZ domain and GluR2 (Figure 4C, blue), and a loss of surface area contact between the PICK1 PDZ domain and GluR2 (Figure 4C, purple). Interestingly, the disruption of PICK1 PDZ-GluR2 interactions is correlated with dynamic changes at the aA helix. Figure 4C (red) shows that the RMSD of the aA helix increases with the dissociation of GluR2. Moreover, our analysis reveals a correlation between PICK1 PDZ-GluR2 interactions and the distance between the bB strand and the aA helix (Figure 4C, green). This distance separation may play a role in the destabilization of the aA helix.
Finally, we calculated the change in the dynamics flexibility index (DDFI) across the bound and unbound states of each system (Figure 4B and4D). DDFI reveals significant changes in dynamics of the PICK1 PDZ domain due to the dissociation of ligands. The important ligand binding regions, including the aB helix and bB strand, show enhanced flexibility upon ligand dissociation. When the interactions are disrupted, the key binding residues gain more conformational freedom, and the flexibility enhances. Thus, enhanced flexibility at the binding site is a direct indicator of a dissociation. More interestingly, DDFI also reveals unique changes to the aA helix upon dissociation of each unique ligand. As represented by the RMSD of the aA helix (Figure 4A, red), the dissociation of DAT does not enhance the flexibility of the aA helix (Figure 4B). Instead, the majority of the aA helix has little change in terms of flexibility while A59 shows enhanced rigidity (Figure 4B). Oppositely, there are significant changes in dynamics of the aA helix due to the dissociation of GluR2 (Figure 4D). Echoing the RMSD of the aA helix (Figure 4C, red) and the distance between I33 and A58 (Figure 4C, green), DFI analysis shows enhanced flexibility at the aA helix upon ligand dissociation (Figure 4D). As the I33-A58 distance increases, the interactions between the aA helix and the carboxylate-binding loop become weaker to allow more fluctuations. Advancing to a dynamically more flexible regime, the aA helix is observed be to allosterically being altered by the dissociation event.
To further explore the correlation between ligand binding and the dynamics at the aA helix, we performed protein network analysis. Protein network analysis can reveal the coupling of major movements by creating protein structure networks based on the primary motions of each residue. The analysis reveals the residues within the PICK1 PDZ domain that are most strongly coupled to the ligands' motion. The motions of DAT (Figure 5A) and GluR2 (Figure 5B) are both coupled to the motion of the distal aA helix and the bB-bC loop of the PICK1 PDZ domain. Interestingly, the motions of DAT are more strongly coupled to the bB and bC strands than are the motions of GluR2.
Dynamic Coupling Index (DCI) was applied to each system to explore the coupling of dynamics between binding site residues and the global protein. The DCI metric has previously been shown to capture allosteric coupling of distal site to critically important residues in a protein. Upon a binding event, the binding site residues experience exerted forces from the ligand so that the dynamics of the system may be affected. Notably, the force exerted by the ligand not only affects the dynamics of the binding site residues but may also affect the dynamics of the global protein due to allosteric communication. The DCI metric measures the coupling strength of a residue to a binding site. A highly coupled residue will experience the repercussions of binding more than weakly coupled residues. As shown in Figure 5C-D, DCI analysis on the PICK1 PDZ-DAT and PICK1 PDZ-GluR2 systems reveals a coupling trend that echoes results from network analysis at the aA helix. Both DAT and GluR2 binding residues observes strong coupling to the aA helix. Time-resolved force distribution analysis (TRFDA) 80 was performed to reveal the punctual stress on each PICK1 PDZ residue as a result of ligand binding. TRFDA was performed over each trajectory, and the per trajectory results were summed over each complex system. The summed results are shown in Figure S7. The ten PICK1 PDZ residues that experienced the greatest punctual stress for each system are listed in Figure S8. Both DAT and GluR2 induce the greatest punctual stress on the bB strand and aB helix, regions that directly interact with the ligands. In the PICK1 PDZ-DAT system, all six residues that experience the greatest punctual stress comprise the bB strand. Oppositely, GluR2 induces significant punctual stress on K83 of the aB helix. These results point to the different interaction patterns induced by different ligands binding.
Our analysis reveals that DAT and GluR2 can induce unique stresses on the PICK1 PDZ domain, but the specific residues and mechanisms through which dynamic allostery is propagated in the PICK1 PDZ domain remains in question. A recent review of allostery in the PDZ family 81 notes that A46 (aA helix) of PTP-BL PDZ2 and A347 (aA helix) of PSD-95 PDZ3 have been consistently identified as allosteric residues in a wide array of computational and experimental efforts. 15,18,85,31,34,35,41,46,[82][83][84] Furthermore, in a recent work exploring the interactions and dynamics between the PICK1 PDZ domain and the small molecule inhibitor BIO124, we propose that a structural alignment of PICK1 PDZ, PTP-BL PDZ2, and PSD-95 PDZ3 suggests that this allosteric alanine residue on the aA helix is evolutionarily conserved across all three PDZ domains. 86 This structural alignment also suggests that the interactions between BIO124 and I35 of the PICK1 PDZ domain may have a role in the propagation of signal to A58 of the aA helix. 86 Notably, A58 forms a van der Waals surface with I35, which is directly involved in ligand binding.
Here, our results support the importance of A58 as an allosteric residue in the PICK1 PDZ domain. Distance analysis reveals that I33-A58 distance is coupled with ligand binding, protein network analysis identifies A58 in the network of residues dynamically coupled to the ligand, and DCI analysis indicates A58 is strongly coupled to binding site residues. We suspect that interactions between natural ligands and I35 of the PICK1 PDZ domain may also have a role in the propagation of signal to the aA helix.
We explore the role of I35 in propagating allosteric signal to the aA helix of the PICK1 PDZ domain. Distance distribution and time-resolved force distribution analysis (TRFDA) are used to identify the degree of interactions between the ligands and I35. As shown in Figure 6A, distance distribution analysis was performed between the ligand and I35 for each system. Here, the distance is defined as the shortest distance between any two atoms in the ligand and I35. DAT (blue) and GluR2 (red) both form the close contact (~2 Å) with I35. In addition to exploring the distance distribution between ligands and I35, we also calculated the punctual stress on I35 induced by the ligand by using TRFDA. As shown in Figure S8, I35 is one of the top five residues that experiences the greatest punctual stress in each system. Figure 6B lists the punctual stress on I35 induced by DAT and GluR2. GluR2 induces a slightly greater punctual stress on I35 than DAT does. As demonstrated by Figure 4, GluR2 is more strongly coupled to the aA helix than DAT is. This stronger coupling between GluR2 and the aA helix may be a result of the strong punctual stress at I35. Together, distance distribution analysis and TRFDA point to the importance of interactions between the ligand and I35 in inducing dynamic allostery at the aA helix of the PICK1 PDZ domain. As discussed in previous work, the dynamic allostery can be closely related to the local conformational changes resulting from local frustrations. To explore the local frustration regions in PICK1 PDZ domains, the Frustratometer server was used. It can be seen from Fig. 7A that the aA helix is indeed a local high frustration region. Moreover, there are other local frustration regions, e.g., aB and bB-bC loop, which contain highly frustrated interactions. Interestingly, both of these two regions were identified in our network analysis (Fig. 5), showing their correlations with the ligands. The tight green lines at the center highlight that the major structural 'core' is conserved. The frustration projection on each residue is shown in Fig. 7B. The ligands are part of the core and, at the same time, trigger frustration on the protein surface.
The purpose of this work is to investigate the dynamic allostery in the PICK1 PDZ domain that can be induced by unique binding partners. We found that (1) the PICK1 PDZ domain exhibits dynamic allostery at the aA helix, (2) the unique interaction patterns between different binding partners and the PICK1 PDZ may induce unique dynamic changes to the PICK1 PDZ domain, and
(3) the hydrophobic core that is formed between the ligands and I35 may be key to inducing dynamic allostery at the aA helix.
Our results demonstrate that natural ligands DAT and GluR2 can induce dynamic allostery at the aA helix of the PICK1 PDZ domain. Protein structure network, DCI, TRFDA, and local frustration analysis show that both DAT and GluR2 are dynamically correlated with the aA helix. This dynamic correlation distant from the binding pocket points to the ability of DAT and GluR2 to induce dynamic allostery across the PICK1 PDZ domain. These results are in agreement with previous work which has identified the aA helix as an allosteric region within other PDZ domains, including Par-6 PDZ, PTP-1E PDZ2, PTP-BL PDZ1, and AF-6 PDZ. 16,[19][20][21]38,46,47 Furthermore, dissociation events captured during our simulations presented a unique opportunity to explore dynamic changes to the PICK1 PDZ domain in real time. GluR2 dissociation is directly coupled with increased fluctuations at the aA helix and increased distance between the aA helix and the bB strand. The distant shift of the aA helix and the bB strand agrees with secondary structure shifts seen in previously studied PDZ domains. 21,22 Notably, the dissociation of the PICK1 PDZ-DAT complex was not so clearly correlated to dynamic changes at the aA helix. These results suggest that different binding partners may induce different dynamic changes to the PICK1 PDZ domain.
Previous work on the PTP-BL PDZ2 domain 17,35 and the PSD-95 PDZ3 domain 13 has pointed to the importance of structural equivalents of I35 in propagating allosteric signal to the aA helix. Our work suggests that I35 may also be a key residue in propagating signals in the PICK1 PDZ domain.
Our results demonstrate that both DAT and GluR2 are dynamically coupled with the aA helix. Distance distribution analysis and TRFDA reveal that DAT and GluR2 form the close contact with and induce the strong punctual stress on I35. These results suggest that interactions between the ligand and I35 are key to inducing dynamic allostery at the aA helix in the PICK1 PDZ domain.
The release of the AlphaFold 2 provides a high-resolution solution 87,88 to compare PDZ domains across multiple species and different proteins, Our results identify dynamic allostery within the PICK1 PDZ domain. By comparing the responses of the PICK1 PDZ domain to the binding of different ligands, we see that the binding of different types of ligands may induce different dynamic changes to PICK1 PDZ domain. Our previous work on the PICK1 protein identified the aA helix of the PDZ domain as a key participant in interdomain PDZ-BAR and PDZ-linker interactions. 89 We suspect that the ligand-induced dynamic changes at the aA helix may affect interdomain interactions and ultimately explain the long hypothesized conformational change of PICK1 upon ligand binding. 49,50 An atomic-level resolution of the mechanism behind the PICK1 interdomain dynamics may greatly affect how we understand the PICK1 protein.