Investigating Phosphorylation-Induced Conformational Changes in WNK1 Kinase by Molecular Dynamics Simulations.

The With-No-Lysine (WNK) kinase is considered to be a master regulator for various cation-chloride cotransporters involved in maintaining cell-volume and ion homeostasis. Here, we have investigated the phosphorylation-induced structural dynamics of the WNK1 kinase bound to an inhibitor via atomistic molecular dynamics simulations. Results from our simulations show that the phosphorylation at Ser382 could stabilize the otherwise flexible activation loop (A-loop). The intrahelix salt-bridge formed between Arg264 and Glu268 in the unphosphorylated system is disengaged after the phosphorylation, and Glu268 reorients itself and forms a stable salt-bridge with Arg348. The dynamic cross-correlation analysis shows that phosphorylation diminishes anticorrelated motions and increases correlated motions between different domains. Structural network analysis reveals that the phosphorylation causes structural rearrangements and shortens the communication path between the αC-helix and catalytic loop, making the binding pocket more suitable for accommodating the ligand. Overall, we have characterized the structural changes in the WNK kinase because of phosphorylation in the A-loop, which might help in designing rational drugs.

Introduction

Hypertension or high blood pressure is a common disorder affecting billion people worldwide.[1] It leads to various other disease conditions, such as myocardial infarction, renal failure, and congestive heart failure. Although available antihypertensive drugs help to control blood pressure, the emergence of the drug-resistant variant has created a global challenge in treating hypertension. Further, the molecular pathogenesis mechanism of hypertension is poorly understood. Gene mutations involved in the pathway of renal salt reabsorption could affect the blood pressure. The first link between the With-No-Lysine (WNK) signaling pathway and hypertension came into the picture with the discovery of an inherited form of hypertension, known as Pseudohypoaldosteronism type II (PHA II), or Gordon’s syndrome, caused by the WNK mutation.[2] Intron deletion in the WNK1 gene as well as the missense mutation in the WNK4 gene cause PHAII.[3] The WNK signaling pathway that is involved in regulating blood pressure participates in the activation of its downstream kinases, such as oxidative stress-responsive kinase 1 (OSR1) and STE20/SPS1-related proline/alanine-rich kinase (SPAK) being phosphorylated by WNK1 and WNK4, which in turn phosphorylate and activate the Na+/Cl– cotransporter (NCC).[4,5] The inhibiting function of WNK is considered as an effective approach for controlling blood pressure.

One of the largest groups of a gene family in eukaryotes involving cell growth and differentiation is the protein kinases (PKs), which carried out phosphorylation reaction by transferring γ-phosphate of adenosine triphosphate (ATP) to a peptide substrate.[6] The catalytic domain of these PKs is highly conserved[7,8] and forms a gene superfamily except few including WNK, sterile serine/threonine kinases (STE), and tyrosine kinase-like (TKL) family.[9] The name of WNK comes from the characteristic feature of the unique placement of catalytic lysine present in the β2 strand, instead of β3 strand as found in other kinases. The WNK kinase domain (KD) comprises an N-terminus (residue 218–483) and two coiled-coil domains in the C-terminus (residue 563–597 and 1814–1841).[9] Along with the KD, it has an autoinhibitory domain between residues 490 and 550 that suppresses kinase activity as well as other different protein-motif interactions. WNK undergoes autophosphorylation at Ser382 and becomes active. It is known as a major autophosphorylation site as the mutant (S382A) yields ∼1–3% of wild-type activity. Further, the residue Ser378 located in the activation-loop is also found to be autophosphorylated to enhance the kinase activity.[10] It is a minor phosphorylation site. The mutant S378A showed ∼50% of the activity of wild-type WNK1. Overall, the major phosphorylation site is Ser382, which is essential for the activity of WNK1.[11]

The overall architecture of WNK resembles other kinases consisting of dual domains: the small N-lobe comprises β strands and an α-helix (αC-helix), whereas the C-lobe comprises an α-helix and a hinge region that acts as a bridge between N- and C-lobes by accommodating ATP and ATP-competitive inhibitors[12] (see Figure ). There is a unique N-terminal domain consisting of six β-sheets forming a full β-barrel. On the contrary, five β-strands are seen in other kinases. This additional β0 forms a hydrogen bond with β1. The N-lobe comprises a glycine-rich loop (GXGXXKXV) from the first two β strands that place the unique catalytic lysine of WNK (Lys233). Various conserved sites near phosphate transfer including a glycine-rich loop (G-loop) and an αC-helix in N-lobe, and the activation loop (A-loop) and catalytic loop (C-loop) in the C-lobe perform the crucial catalytic function. A critical catalytic triad Asp-Leu-Gly (DLG) (residues 368–370) motif is present in the activation segment (residues 360–390) in the C-lobe.[13] The usual salt-bridge between Lys72 and Glu91 (residue numbers corresponding to protein kinase A, PKA) has been observed in many kinases, and this Lys72 interacts with α and β phosphates of ATP and involves in protonating ADP or related acid–base reactions.[14−16] However, in WNK1, Cys250 occupies this position in subdomain II (β3 strand) although its side chain is positioned away from the active site. Instead, Lys233 in β2 strand of subdomain I reaches the active site and plays the role of catalytic lysine (see Figure ). Previously, mutagenesis study has demonstrated that the mutation at Lys233 abolishes the catalytic activity of WNK1.[9] For most of the kinases, primary regions involved in controlling the activation process are the Asp-Phe-Gly (DFG) motif, A-loop, and the αC-helix.[17,18] For many kinases, such as IGF-1RK, it has been observed that the autophosphorylation of A-loop leads to an extended conformation for accommodating ATP and substrate and thus forms the catalytically active kinase.[19,20] It implies that the kinases interconvert between an active and an inactive conformation because of phosphorylation at the A-loop site.[21−23]

Figure 1

(A) Structure of WNK1 KD, N-lobe with β-strands and αC-helix, and C-lobe exclusively with α-helices. Different structural parts are shown including glycine-rich loop (light blue), hinge region (orange), αC helix (blue), activation segment (red), catalytic loop (cyan), and DLG motif (olive green). Ligand (WNK463) in a ball and stick model (purple color). (B) C-spine and R-spine residues are shown with green and red, respectively.

Recently, Yamada and co-workers have resolved the crystal structure of unphosphorylated WNK1 complexed with an inhibitor, WNK463 (PDB code: 5DRB).[24] The 3-D structure of the inhibitor is shown in Figure S1 (Supporting Information). They have found that the A-loop and αC-helix lie away from the hinge region compared to apo-WNK1 to accommodate the inhibitor and adopts a DLG-in and an αC-helix-out (i.e., inactive) conformation. Previously, many experimental and theoretical studies have been conducted elucidating the effects of ligand binding and mutation on kinase activities, dynamics, and conformational changes.[25−28] However, the molecular mechanism and impact of autophosphorylation, that is, A-loop phosphorylation on its structure and dynamics, is not known from the crystal structure alone. Therefore, we conducted microsecond long unbiased MD simulations to elucidate the effect of phosphorylation on structural dynamics of WNK kinase and investigated the biophysical basis of binding of WNK463 to the unphosphorylated and phosphorylated WNK.

Results and Discussion

Structural Stability and Flexibility Analysis

First, root-mean-square deviations (rmsds) of Cα atoms from the initial structure were computed to find the overall stability and convergence of simulations. From Figure A, it is observed that rmsds of uWNK and pWNK fluctuate initially, and then become stable from 400 ns onward, indicating the trajectory convergence. It is evident from the figure that in the case of uWNK, the system drifts from its initial structure till 100 ns and again there is a sudden jump in rmsd at 360 ns, and then slowly it reaches a plateau after 480 ns. On the contrary, the pWNK system reaches equilibrium soon after 200 ns and remains stable after that. The average rmsd values of uWNK and pWNK are found to be 2.5 and 1.7 Å, respectively. An intermediate rmsd of 2.1 Å is obtained for 2pWNK (see Table ). However, a similar trend in Cα rmsd is observed in cases of u-Apo and p-Apo systems (see Figure S2).

Figure 2

(A) Time evolution of the rmsds of Cα atoms with respect to the initial structure for uWNK (red), pWNK (blue), and 2pWNK (gold). (B) Superimposition of uWNK and pWNK final structures from the simulations. (C) RMSFs of the backbone atoms of uWNK (red) and pWNK (blue) complexes during simulations. Also, the RMSF of the crystal structure (black) is included.

Table 1

Average C rmsds of the Protein, R-Spine, C-Spine, Radius of Gyration (Rg), and Average Distances for Ser378, Ser382, and Thr386 from the Simulations in Åa

systems	⟨protein⟩	⟨R-spine⟩	⟨C-spine⟩	⟨Rg⟩	D1b	D2c
uWNK	2.5 (0.1)	0.9 (0.2)	0.6 (0.1)	18.9 (0.1)	9.6 (1.2)	8.9 (1.4)
pWNK	1.7 (0.1)	0.6 (0.3)	0.5 (0.1)	18.8 (0.1)	12.5 (0.5)	11.7 (1.1)
2pWNK	2.1 (0.2)	0.8 (0.3)	0.5 (0.1)	18.7 (0.1)	8.9 (0.7)	11.2 (1.1)
u-Apo	1.8 (0.1)	0.6 (0.2)	0.8 (0.3)	18.7 (0.1)	6.2 (0.4)	7.3 (0.9)
p-Apo	1.8 (0.2)	0.8 (0.2)	0.9 (0.3)	18.8 (0.1)	10.4 (1.1)	10.7 (1.4)

Standard deviations are given in the parenthesis.

D1 = Avg. Cα distance of Ser378-Ser382.

D2 = Avg. Cα distance of Ser382-Thr386.

The structural compactness can be analyzed by calculating the radius of gyration (Rg) of the WNK complexes. The average Rg calculated for pWNK and 2pWNK were found to be the same (∼18.7 Å), whereas in the case of uWNK, a slightly higher value of 18.9 Å was obtained (see Table ). The frequency distribution of Rg for the complexes are shown in Figure S3. Further, we superimposed the final structural coordinates from the simulations of uWNK and pWNK using Superpose server[29] to observe the structural changes as shown in Figure b and found that major changes have occurred in the A-loop only. The atomic positional fluctuations or root-mean-square-fluctuations (RMSFs) describe the flexibility of an individual residue. RMSFs of protein Cα atoms of uWNK and pWNK during the simulations are shown in Figure c. It is evident from Figure c that the activation segment in uWNK is more flexible than in pWNK. Furthermore, loop regions between β4−β5 and β7−β8 show higher flexibility in uWNK compared to pWNK. A higher degree fluctuation regions in the complexes indicate that it may be involved in the inactivation mechanism and substrate binding.[30]

Structural Changes in the Activation Segment

The time evolution of the Cα rmsds of A-loop from the initial structure is shown in Figure A. From Figure A, uWNK shows drifting in the initial period of simulation with a sudden jump at 380 ns, which then converges to an equilibrium state and reaches a plateau after 400 ns. In contrast, for pWNK, the rmsd reaches an equilibrium immediately after 150 ns and then stays stable thereafter. The average Cα rmsd of the A-loop for uWNK and pWNK are determined as 4.7 and 2.3 Å, respectively, that is, the average rmsd of uWNK is twice that of pWNK. rmsd values show a massive configurational change in the A-loop after phosphorylation. Similarly, in the 2pWNK system, the average Cα rmsd of A-loop is found to be 3.2 Å, which falls between that of uWNK and pWNK systems.

Figure 3

(A) Time evolution of Cα rmsd of the activation segment with respect to the initial structure for uWNK (red), pWNK (blue), and 2pWNK (gold) systems. Dominant conformations of the activation segment obtained from the cluster analysis for uWNK (B) and pWNK (C) systems. (D) Position of conserved Thr386 in uWNK (green cartoon) and pWNK (pink cartoon) and the average distance between Thr386-Ser382 in both complexes.

Furthermore, to extract the predicted structures from the simulated trajectories and to see structural variations in the A-loop, we have performed clustering analysis using k-means clustering algorithm.[31,32] Clustering was done based on the Cα rmsd distance metric. A representative frame in a cluster is defined as the smallest average distance from other frames that also belongs to that cluster. As the number of clusters in k-mean clustering is defined arbitrarily, and Figure A suggests three possible distinct states in uWNK, we have chosen three clusters for the k-means clustering analysis for both uWNK and pWNK. The size of the cluster with respect to the total trajectory represents the percentage of the dominant conformations of the A-loop out of three clusters for both the systems and is shown in Figure B,C for uWNK and pWNK, respectively. Clusters 1, 2, and 3 in uWNK correspond to 58.6, 27.8, and 13.6%, respectively. Similarly, for pWNK, clusters 1, 2, and 3 correspond to 82.4, 10, and 7.6%, respectively. It suggests that the unphosphorylated A-loop undergoes dynamic fluctuations, resulting in more distinct conformations; on the other hand, phosphorylation has stabilized the A-loop into a single distinct conformation. The dominant conformation of uWNK is outwardly displaced compared to the activation segment of pWNK. Overall, the phosphorylation at Ser382 stabilizes the A-loop significantly compared to uWNK and 2pWNK. However, a similar trend in Cα rmsd of the A-loop for u-Apo and p-Apo was observed as seen in Figure S4. A similar observation was made in other kinases, such as IRAK[30] and IGF-1RK,[33] suggesting that upon phosphorylation in A-loop, it may allow substrate access and activate the kinase.[34] Furthermore, the conformational changes in the A-loop are described by the free energy landscape (FEL) of rmsd versus Rg of the A-loop region and is shown in Figure S5. The FEL described the conformational space covered by the A-loop. It shows that in uWNK there are two major conformational states characterized by rmsd ∼4.5 Å/Rg ∼9.5 Å and rmsd ∼2.1 Å/Rg ∼10 Å, whereas for pWNK there is a single broad conformational space explored with rmsd ≈ 2.2 Å/Rg ≈ 9.7 Å, indicating that more deviation is observed in the A-loop of uWNK compared to pWNK. Also, another FEL as a function of rmsd of A-loop versus rmsd of αC-helix is drawn and shown in Figure S6. This again supports the above observation.

To further characterize the induced structural variation because of phosphorylation, the Cα–Cα distance between Ser378 and Ser382 is calculated for all systems and found to vary between 6.2 and 12.5 Å (see Table ). For u-Apo, an average distance of 6.2 Å is obtained, whereas 10.4 Å is obtained for the p-Apo system. Similarly, a relatively higher distance is obtained for pWNK (12.5 Å) compared to uWNK (9.6 Å). A similar trend is observed for the Cα–Cα distance between Ser382 and Thr386 (see Table ). Overall, this suggests that the A-loop adopts an extended conformation upon phosphorylation at Ser382. Furthermore, we have estimated the solvent-accessible surface area (SASA) of Ser378, Ser382, and Thr386 for all systems, and the corresponding average value is reported in Table S1 (Supporting Information). It is evident from Table S1 that the average SASA of Ser382 and Thr386 increases significantly because of phosphorylation at Ser382. This may explain the extended conformation of the A-loop because of phosphorylation.

Dynamics of Regulatory and Catalytic Spines

The regulatory spine or R-spine, also known as the hydrophobic spine, consists of hydrophobic residues between N- and C-lobes and forms a conserved motif including Leu299, Leu272, His347, and Leu369 from the β5 strand, αC-Helix of N-lobe, HRD, and DLG motifs from the C-lobe. This spine is least exposed to the solvent and arranged contiguously in the case of an active kinase, whereas in an inactive kinase, the structure disintegration occurs in part formed by the αC-Helix and A-loop.[35] In Figure S7, the probability distribution of Cα rmsd from the initial structure is shown. The average Cα rmsd for uWNK and pWNK over the simulations were estimated as ∼1.0 and ∼0.6 Å, respectively. An intermediate Cα rmsd of 0.8 Å was obtained for the 2pWNK system (see Table ). Again, a slightly lower rmsd was obtained for pWNK compared to uWNK and 2pWNK, suggesting that the unphosphorylation causes slight destabilization to the R-spine.

There is another stretch of hydrophobic residues between N and C-lobes (such as Val235, Ala248, Leu313, Ile355, Ile357, Leu413, and Ala416), known as catalytic spine or C-spine, which forms an interaction with ATP. The distribution of Cα rmsd of the C-spine (see Figure S8) is quite stable in both the uWNK and pWNK complexes with an average rmsd of ∼0.5 Å. A similar average rmsd of 0.5 Å was obtained for 2pWNK (see Table ), which suggests that in contrast to the R-spine, the C-spine remains unaffected because of phosphorylation.

Salt-Bridge between the αC-Helix and Catalytic Loop

The WNK kinase does not have a usual Lys–Glu ion pair between the αC-helix and β3 strand, which is evident from Figure S9 (Supporting Information). The average distance between NZ atom of Lys233 and CD atom of Glu268 was found to be 17.2 and 15.7 Å for uWNK and pWNK, respectively. It is in agreement with a previous observation[24] where this distance was found to be ∼13 Å for uWNK. However, a different ion pair between Glu268 of the αC-helix and Arg348 of the catalytic loop was observed, and it was formed transiently in uWNK as shown in Figure a, which was also seen in other kinases, such as IGF-1RK33. On the contrary, this salt-bridge was quite stable throughout the simulations in case of pWNK. Average distances of 7.8 and 3.5 Å were found for uWNK and pWNK, respectively (see Figure a). Figure b shows the probability of occurrence of the salt-bridge throughout the simulations, and it was observed that pWNK has a higher probability of maintaining the salt-bridge compared to uWNK and 2pWNK. Hence, our simulation results suggest that phosphorylation at Ser382 only forms the stable salt-bridge and helps to activate the kinase. Also, the interactions between Glu268 and Arg348 for various atoms are shown in Figure S10, suggesting a close hydrogen bond interaction in case of the pWNK compared to the uWNK system.

Figure 4

(A) Salt-bridge distance calculated between CD of Glu-268 and CZ of Arg-348 from the simulations of the complexes uWNK (red), pWNK (blue), and 2pWNK (gold). (B) Probability distribution of the salt-bridge distances of the complexes.

Further, we have calculated the native contacts between the αC-helix and catalytic loop. Residues 265–271 from the αC-helix and residues 345–352 from the catalytic loop are considered for the calculation of native contacts and the results are shown in Figure S11. It is evident from Figure S11 that the number of native contacts obtained in the case of pWNK is relatively higher than in uWNK. Further, it is evident from Figure S11 that the fraction of native contacts between Glu268-Arg348 is also higher in pWNK (∼71%) than in uWNK (∼21%).

To complement the above results, we studied the effect of phosphorylation on the size and shape of the ATP binding pocket using CASTp,[36] and Chimera tool,[37] and the last configurations were used for this purpose. The calculated volumes for uWNK and pWNK were found to be 640.9 and 586.9 Å3, respectively. Overall, it signifies that the A-loop phosphorylation alters the ATP binding pocket (see Figure S12) by contracting the pocket volume, thus signifying that salt-bridge formation decreases the pocket size and alters the structure and function of the kinase.

Interactions Affecting the Dynamics of the A-Loop and αC-Helix

We have identified several interactions (see Table ) affecting the dynamic behavior of the A-loop and αC-helix. Critical electrostatic interactions between the phosphoserine residue of the A-loop with the adjacent charged residues affect the folding of A-loop and help in positioning substrate binding. Besides the Arg348-Glu268 salt-bridge interaction discussed above, other interactions that were affected after phosphorylation at Ser382 in the A-loop were studied. From the interaction index graph (see Figure S13), it can be suggested that in the case of pWNK, the interaction between pSer382 and charged residues of the HRD motif (Arg348 and Asp349) decreases significantly after phosphorylation. However, the interaction of pSer382 with the Arg255 and Lys256 increases after phosphorylation near the αC-helix (see Figure a,b). Also, different atoms involved in the interaction between pSer382 and Arg255 as well as other interactions from A-loop regions were calculated and are shown in Figure S14, correlating well with the observation that phosphorylation imparts electrostatic effects on nearby residues, resulting in a conformational change. Overall, it suggests that in the WNK kinase, the interaction between the A-loop and αC-helix contributes more toward the stability of its active conformation and influences the orientation and dynamics of the αC-helix. In addition to interactions with pSer382, other interactions such as the interaction between Glu268 and Arg264 also influence the orientation of the αC-helix. They form an intra-helix salt-bridge in uWNK (see Figure c). After phosphorylation, Arg264 disengaged from the Glu268, and Glu268 then reoriented toward the active site (see Figure d), and formed a salt-bridge with Arg348, leading to the proper location of catalytic residues and substrate for the phosphoryl transfer reaction. Further, it explained the association of the αC-helix with the A-loop in pWNK. Similar conserved interactions were observed in other kinases,[33,38,39] supporting the role of interactions of the αC-helix with the A-loop for the catalytic activity of the kinase.

Table 2

Average Interaction Distances between Selected Residue Pairs for Unphosphorylated and Phosphorylated Kinases in Åa

interaction index	residue 1	residue 2	distance (uWNK)	distance (pWNK)
1	Ser382(OH)	Lys256(NZ)	19.8 (2.8)	9.3 (3.2)
2	Ser382(OH)	Arg255(CZ)	13.0 (2.6)	5.1 (2.0)
3	Asp349(CG)	Lys375(NZ)	3.6 (0.6)	4.1 (1.0)
4	Arg348(CZ)	Glu268(CD)	7.4 (1.1)	4.2 (0.6)
5	Glu268(CD)	Arg264(CZ)	4.8 (1.2)	8.0 (1.7)
6	Ser382(OH)	Arg348(CZ)	12.7 (1.9)	20.1 (1.5)
7	Ser382(OH)	Arg376(CZ)	6.2 (2.4)	10.4 (2.0)
8	Ser382(OH)	Asp349(CG)	11.4 (1.6)	15.8 (1.5)

Standard deviations are given in the parenthesis.

Figure 5

(A,B) Interaction distances calculated between S382(OH)-R255(CZ) (blue), S382(OH)-K256(NZ) (orange), R348(CZ)-E268(CD) (red), and E268(CD)-R264(CZ) (green) pairs for uWNK and pWNK, respectively, from the simulations. (C,D) 3-D view of the switch of the salt-bridge interactions between the Glu268(CD)-Arg264(CZ) and Glu268(CD)-Arg348(CZ) pairs for uWNK (grey) and pWNK (sea green), respectively.

Structural Network Analysis

Protein stability or protein function relies on interactions between residues. In protein contact network (PCN), amino acids are considered as nodes and edges between these nodes can be drawn in terms of distance or interaction energy that depicts the interactions between the residues. Previously, various structural network approaches have been reported that provide the structure–function relationship of proteins.[40−42] These topological studies revealed the important residues involved in protein stability,[43] protein dynamics,[44] allosteric regulations,[45] and signal transduction.[45] Different parameters such as shortest path distance in the network analysis of different proteins could depict the influence of binding of the inhibitor in the binding pocket of different systems and deduce the structural significance of the induced-fit mechanism in protein-inhibitors.[46] Similarly, in the WNK kinase, to study how the αC-helix and catalytic loop communicate in uWNK and pWNK and how the communication affects the different interactions were studied via network analysis of protein structures (NAPS).[47] The shortest distance path on the protein network structures was analyzed by considering the salt-bridge interaction between Glu268 (αC-helix) and Arg348 (at catalytic loop) and used as a terminal residue for calculation of the shortest path in NAPS and shown in Figure . The shortest path distance between Glu268 and Arg348 in uWNK spanned three residues (Pro343, Ile345, Ile346) (Figure a); however, in pWNK it spanned two residues (Gly370 and Ala372) (Figure b). This result suggests that the phosphorylation modifies the shortest communication path between the αC-helix and catalytic loop, making a region suitable for accommodating the ligand. Furthermore, the betweenness parameter was used to study the correlation between residues in the protein network, and the betweenness values of uWNK and pWNK are shown in Figure S15A. Residues showing higher betweenness centrality (BC) values in both systems are reported in Table . It is evident from Figure S15A that the betweenness of residues in pWNK show higher values compared to uWNK, especially in the N-terminal, αC-helix, and catalytic and A-loop, indicating that upon phosphorylation the binding pocket regions rearranged and made more related residues more connected. A similar observation was found in u-Apo and p-Apo (see Figure S15B) with higher residue–residue relevance in the case of p-Apo compared to u-Apo. Wilson et al.[48] suggested that drug binding leads to conformational changes via an induced-fit process. On the basis of the betweenness value, it reflects that the ligand interacts with the residues of the pre-organized pocket, and phosphorylation causes the stabilization of the binding pocket through the induced-fit mechanism.

Figure 6

Shortest communication path between the αC-helix and catalytic loop investigated via NAPS in (A) uWNK and (B) pWNK systems.

Table 3

Residues Showing High BC in Unphosphorylated, Phosphorylated, and Double Phosphorylated Systems Calculated from NAPS constructed from MD Simulations

system	residues
uWNK	V281, T301, C411, V247, N354
pWNK	A372, K351, E302, I280, D349, K375, V247, N279, M390, F232
2pWNK	A372, L374, D349, V247, E302, V281, L252

Cooperative and Anti-Cooperative Motions

To investigate how phosphorylation affects the domain motions, we generated cross-correlation maps using Cpptraj module(49) of AMBER tools from the conformations obtained by the molecular dynamics simulations. In the apo systems, Figure A,B), more correlated motions were observed in u-Apo compared to p-Apo and interestingly when the inhibitor binds the anti-correlated motions increase, whereas upon phosphorylation again it diminishes the anti-correlated motion and increases the correlated motions. In case of uWNK and pWNK (Figure C,D), the regions (i) β1−β3, (ii) hinge, (iii) αEF, (iv) αG-helix, and (v) C-terminal tail (Ala448-Gln480) are highly correlated. In uWNK, the interdomain correlated motions were observed in different regions, such as the R1 region in Figure C showed correlated residue interaction between the catalytic loop with the A-loop and near catalytic residue (HRD motif, residues 348–350) with a hinge region (residues 303–307). The anti-cooperative interactions were observed with the catalytic loop region (residues 348–350) with the A-loop (residues 376–384 in R2 region) as well as with the αC-helix region (R3 region) as supported above with the higher fluctuation and more considerable distance among these residues. When WNK KD was phosphorylated at Ser382 (Figure D), the observed highly anti-correlated motions seen in R2 and R3 regions of uWNK was diminished and transformed toward correlated motion. Overall, correlated as well as anti-correlated domain motions were modified after the phosphorylation at Ser382 of A-loop and affect the functional dynamics of WNK. In particular, all the interactions explained in the previous section were justified by the cross-correlation map such as the salt-bridge between the catalytic loop (Arg348) and αC-helix (Glu268), as well as the distance between different residues that control the stability within the αC-helix (Glu268-Arg264), and Ser382 with Arg255 and Lys256.

Figure 7

Cross-correlation matrices of the fluctuations of the coordinates for Cα atoms around their mean positions for Apo and complexes are shown in (A) u-Apo, (B) p-Apo, (C) uWNK, and (D) pWNK, respectively, from simulations. The extent of correlated and anticorrelated motions are color-coded (1 = highly correlated; 0 = no correlated; −1 = anticorrelated).

Principal Component Analysis

The principal component analysis (PCA) results in terms of eigenvalues (3N, 3 × 273 = 819) versus eigenvectors obtained after the diagonalization of the atomic fluctuations covariance matrix for uWNK and pWNK are shown in Figure S16. The corresponding eigenvalues of each eigenvector were plotted in decreasing order. The first few eigenvalues depict the collective motions of the localized fluctuations. Comparing the two systems, it indicates that the first few principal components that described the properties of motions were not the same, and in uWNK, the magnitudes of the eigenvalues were higher compared to pWNK. The first few 20 eigenvectors describe the collective motions approximately, and, the first 10 eigenvectors account for ∼85, and ∼79% overall motions in uWNK and pWNK, respectively. Similarly, the first two eigenvectors capture ∼47 and ∼32% of the total motions in uWNK and pWNK, respectively. Our result suggested the correlated motions were higher for uWNK in agreement with the cross-correlation results. FEL of the complexes for the first two principal components is represented in Figure . The pattern of a basin on the FEL is different in both the complexes, as is evident from Figure . In the case of uWNK, the conformational space showed a multiple dispersed basin, whereas in the case of pWNK, it showed a single broad and stable global minimum confined within a particular basin on the free energy surface. It suggests that for the structural disparity in both complexes, uWNK attains more conformations compared to pWNK.

Figure 8

FELs generated by projecting the principal components, PC1, and PC2 of uWNK (A) and pWNK (B) complexes in MD simulations at 300 K. The free energies are represented by −KBT ln PPC1,PC2 with PPC1,PC2 being the probability distribution.

Furthermore, the direction of motions and magnitude of selected eigenvectors can be visualized by using “porcupine plots”[50] to illustrate the effect of phosphorylation at the A-loop on the WNK kinase conformation. Porcupine plots, as shown in Figure , were generated using VMD,[51] corresponding to the first mode that shows the largest concerted motions of Cα atoms which had been mapped onto the average structure. For easy visualization, the cutoff was set to 4 Å. Figure shows that the overall direction of motions of uWNK changes after phosphorylation of the A-loop in pWNK, which agrees with the PCA plot (Figure ) and RMSF plot (Figure C). Overall, it is observed that the mobility is higher for uWNK compared to pWNK as can be ascertained from the length of the eigenvectors (see Figure ).

Figure 9

Porcupine plots showing prominent motions for (A) uWNK (blue cartoon) (B) pWNK (red cartoon). Yellow and green represent eigenvectors showing the direction of prominent motions for uWNK and pWNK, respectively. Length of the eigenvectors represents the magnitude of the motions.

Dihedral Principal Component Analysis

As compared to cartesian coordinates based PCA, which may not provide fully the correct internal and overall motions, hence to get more detail picture about the internal motions and dynamic conformations of the A-loop region after phosphorylation, we performed the dihedral principal component analysis (dPCA) (see Figure ). Figure depicts that in both complexes, dPCA exhibits a large number of minima showing the different specific conformations. In agreement with PCA results, a large number of local conformations separated by minima was observed for uWNK, whereas two major populated conformations were found in pWNK.

Figure 10

FEL of A-loop for (A) uWNK and (B) pWNK complexes generated from dPCA. Representative structures are also shown.

Energetics of Ligand Binding

To understand the biophysical basis of molecular recognition of WNK463 by the WNK kinase, all individual contributions to the total binding free energy were calculated for both uWNK and pWNK complexes using the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) scheme. A summary of the contributions of the binding interactions of the inhibitor (WNK463) with the kinase is shown graphically in Figure S17, and data are provided in Table . It is evident that the binding is favored by the van der Waals and intermolecular electrostatic interactions. This was already shown in our previous study.[52] Polar solvation energy and the configurational entropy mainly disfavor the complex formation. Results showed that the slightly higher binding affinity for pWNK compared to uWNK is due to lesser contributions from the disfavorable components of the binding free energy, such as polar solvation energy and entropy. Moreover, the overall polar contributions from the sum of electrostatic and polar solvation energy show positive results depicting that the van der Waals interaction is highly dominant for the complex formation. To ensure the above results, we have also calculated the H-bond occupancy of both the systems throughout the simulations, shown in Figure S18. It was in agreement that higher electrostatic interactions in uWNK compared with pWNK were supported by the higher percentage of H-bond occupancy between kinase inhibitor (KI) and particular residues of protein such as Met304@N-KI@N37, Lys233@N-KI@O42, and Asp368@N-KI@H312, respectively. Subsequently, we have also performed the decomposition of the free energy of individual residues for both the systems, shown in Figure S19. In uWNK, the remarkable residues contributing to the binding free energy comes from Lys233 (−3.0 kcal/mol), Phe283 (−2.4 kcal/mol), Met304 (−2.3 kcal/mol), Phe356 (−2.3 kcal/mol), Val235 (−2.2 kcal/mol), and Leu369 (−1.0 kcal/mol), and in pWNK major contributions come from the residues Phe283 (−2.3 kcal/mol), Met304 (−2.3 kcal/mol), Val235 (−2.2 kcal/mol), Phe356 (−2.1 kcal/mol), Lys233 (−1.6 kcal/mol), and Leu369 (−1.3 kcal/mol), respectively.

Table 4

Energetic Components of the Binding Free Energy in Unphosphorylated and Phosphorylated Systems with WNK463 (Ligand) in kcal/mola

component	uWNK	pWNK
ΔEvdW	–66.9 (0.1)	–65.5 (0.1)
ΔEele	–29.6 (0.2)	–26.1 (0.2)
ΔGpol	59.3 (0.1)	55.0 (0.2)
ΔGnp	–5.2 (0.0)	–5.1 (0.0)
ΔGb	29.7 (0.1)	28.9 (0.2)
–TΔS	26.0 (1.2)	24.9 (1.5)
ΔGbind	–16.3 (0.1)	–16.7 (0.1)

Standard errors of the mean are provided in parentheses.

ΔGpol,ele = ΔGpol + ΔEele.

WNK–Ligand Interactions

One of the important interactions between a ligand and protein is the hydrogen bond that gives specificity for ligand recognition by proteins. For hydrogen-bond analysis, the acceptor to donor distance cutoff was taken as 3 Å, and angle cutoff was set to 135°. To understand these interactions between kinase and the ligand (WNK463), the lowest energy structures obtained from the MD simulations were considered. Kuenemann and Fourches[53] had observed various hydrogen bond and hydrophobic interactions such as between imidazole nitrogen of ligand and backbone Met304 and with Asp368 and Leu369, as well as multiple hydrophobic interactions such as Val235, Ala248, Phe256, Leu272, Leu299, Leu371, and Phe356. We also have observed the similar interactions between imidazole nitrogen and the backbone Met304. Interestingly, we also found the hydrogen bond between amide nitrogen of Lys233 with the oxygen of WNK463 from our MD simulations, suggesting that conservation of these hydrogen bonds is crucial for the binding of ligand and designing of WNK463 analogs against WNK could be reinforced by increasing the conservation of the interactions with the catalytic residue Lys233 and gatekeeper residue Met304. Common hydrophobic interactions observed in complexes as shown in Figure S20 were calculated from the LigPlot,[54] which include Phe283, Leu369, Gly367, Val281, Asp368, Glu302, Val235, Ala248, Gly228, Thr301, and Phe356. Also, we have found the conserved π–π stacking interaction (see Figure S21) between Phe283 and oxadiazole group of WNK463 throughout our simulations as shown in Figure S22 and were seen as constant in both systems with an average distance of ∼4.3 Å in accordance with the experimental structure.[24]

Conclusions

In this work, 1 μs long MD simulations were performed to investigate the effect of phosphorylation on the structural dynamics of the KD by phosphorylating Ser382 of the A-loop. Simulation results suggest that the phosphorylation leads to the stabilization of the flexible A-loop by the formation of multiple side-chain interactions and affects the protein motions in different segments of the KD, such as the αC-helix, catalytic loop, and A-loop. Stable salt-bridge interaction between Glu268 and Arg348 was observed in pWNK. Other interactions such as between Glu268 and Arg264 was formed only in the case of uWNK, which after phosphorylation, Arg264 disengaged from Glu268, and the side chain of Glu268 reoriented toward the active site, and formed a salt-bridge with Arg348, for the appropriate phosphoryl transfer reaction. Structural network analysis reveals that the phosphorylation causes structural rearrangements and shortens the communication path between the αC-helix and catalytic loop, making binding pocket more suitable to accommodate ligand. Free energy basin clearly showed there was no conformation selection for uWNK compared to pWNK. Our simulations’ study also predicts that the monophosphorylation at Ser382 significantly stabilizes the structure compare with the other systems and makes the WNK kinase active in agreement with the experimental findings.[10,55]

We have also computed and compared the binding affinity for the pWNK and uWNK systems. Electrostatic energy and van der Waals interactions, as well as nonpolar components, always favor the binding, on the other hand, polar solvation energy and entropy component always disfavor the binding. Here, the higher binding affinity of pWNK than uWNK is due to the fewer contributions from the disfavor components of binding free energy such as polar solvation energy and entropy energy. Moreover, the overall polar contributions from the sum of electrostatic and polar solvation energy show positive results, depicting that the main driving force for the complex formation is provided by the van der Waals interaction. Overall, this kind of computational study helps us to understand the structure to function relationship of proteins to know more about the mechanisms that aid in drug development.

Materials and Methods

Summary of the WNK Systems

In the present study, WNK1 KD with WNK463 was taken from the crystallographic structure (PDB ID: 5DRB),[24] while keeping all the crystal water molecules. The following systems were investigated: (i) unphosphorylated complex, designated as uWNK (ii) phosphorylated complex, where Ser382 was modified to phosphorylated form, designated as pWNK, and (iii) double-phosphorylated system, where Ser382 and Ser378 were phosphorylated and designated as 2pWNK. Also, Apo-unphosphorylated (u-Apo) and Apo-phosphorylated (p-Apo) systems were studied. To understand the phosphorylation-induced conformational changes, we have simulated the above systems for 1 μs and also investigated the ligand-binding onto the structural conformation.

Molecular Dynamics Simulations

Molecular dynamics simulations of all the systems were carried out using Amber16 MD package.[56] All simulations were performed using GTX-1070 GPU. Amber ff14SB force field[57] was used to describe proteins. The phosphorylated serine was described using phosaa10 force field.[58] Generalized Amber force field[59] was used for the ligand, and AM1-BCC[60] atomic charges were calculated using the antechamber[61] module of Amber.[56] Each system was placed in an octahedron periodic box with a distance of 10 Å from the walls and solvated using a TIP3P[62] water model. Then, the systems were neutralized by adding an appropriate number of Cl–/Na+ ions. SHAKE[63] algorithm was used to constrain the hydrogen atoms. The system’s temperature was maintained at 300 K using a Langevin thermostat.[64] The long-range electrostatic interactions were calculated using the particle-mesh Ewald[65] scheme involving nonbonding interactions cutoff at 10 Å. Simulations of all the systems were performed as follows: first, minimization was run for 500 steps with the steepest descent and then 500 steps of a conjugate gradient algorithm restraining the solute. Then, the second minimization was run without any restraint for 100 steps using the steepest descent followed by another 900 steps using the conjugate gradient method. Further, equilibration was performed with positional restraint of 5 kcal mol–1 Å–2 under NVT ensemble at 300 K for 50 ps. Then, under NPT ensemble simulations were performed for 50 ps at 1 atm using Berendsen barostat.[66] Subsequently, without any restraint, the systems were equilibrated for 1 ns. Finally, the production run was performed for 1 μs under NPT at 300 K and 1 atm without any restraint. We have used a time-step of 2 fs. Coordinates were saved every 10 ps, leading to 100 000 configurations.

MM-PBSA Calculation

The binding affinity in terms of binding free energy of the complexes was evaluated in explicit water using the MM-PBSA method.[67−73] For the free energy calculation, 10 000 frames from the last 200 ns of simulations were taken for the calculation of binding free energy using the following equationwhere Gcomplex, Greceptor, and Gligand in eq represent the average Gibbs free energy for the complex, receptor, and ligand, respectively. Summation of the molecular mechanic’s free energy (ΔEMM), the solvation free energy (ΔGsol), and the entropy (TΔS) describe the Gibbs free energy as

Molecular mechanic’s energy (EMM) can further be decomposed intowhere Ecov, Eelec, and EvdW represent the covalent, electrostatic, and van der Waals interactions, respectively. Further, the covalent component includes the bond energy (Ebond), the angular vibrational energy (Eangle), and the dihedral angle energies (Edihedral). The solvation free energy (ΔGsol)[71,74] comprises polar as well the nonpolar free-energy component given by

The polar contributions (Gpol) are obtained from the Poisson–Boltzmann (PB) equation with dielectric constants of 1.0 and 80.0 for the solute and solvent, respectively. The nonpolar solvation free energy Gnp is estimated from the solvent-accessible surface area or SASA as given bywhere γ = 0.00542 kcal·mol–1 Å and b = 0. SASA was calculated with a probe radius of 1.4 Å by using a fast–linear combination of pairwise overlap algorithm.[75] The normal mode analysis[76] was used to calculate the entropy (SMM) part and for which 20 configurations were taken from the last 200 ns simulations.

Structural Network Analysis

The structural network was constructed via the NAPS portal.[47] The last trajectory of the simulations was used for the representative structures. NAPS provides an interactive visualization of the PCNs. The construction of networks can be done based on Cα, Cβ, atom pairs, centroid network, or interaction-strength network. Here, in the network, the nodes are represented by the Cα atom of the amino acid residues and the edge between Cα–Cα pair of residues is drawn with a threshold distance taken as Rc (∼7 Å). Floyd–Warshall algorithm[77] was used to calculate the shortest path between the residues in the network. The parameter of the shortest path is considered to describe the long-range interactions in proteins,[78,79] whereas the shortest path is defined by the minimum number of nodes required to reach from one node to other. The betweenness of a node represents the global/local centrality of the node represented by the number of shortest paths passing through the given node in a network. The residue connectivity for uWNK, pWNK, as well as for apo systems is shown by the betweenness parameter to show the structural rearrangements of the binding pocket induced upon phosphorylation as well as by the inhibitor binding.

PCA[80] describes the functional significance of correlated motions. The covariance matrix from the atomic fluctuations of protein is diagonalized and represented in terms of eigenvectors and eigenvalues that represent the sets of principal components which could be used to describe the motions.[81] The direction of motions is represented by the eigenvectors, and the amplitude of motions is given by the corresponding eigenvalues. Our analysis yielded 3 × N, 3 × 273 = 819 eigenvectors.

Dihedral Principal Component Analysis

It is a modified form of PCA which employs internal coordinates, such as dihedral angles (φ, ψ) instead of cartesian coordinates.[82] For the biomolecular study, considerations of dihedral angles are more appealing than bond angles and bond lengths as they do not undergo much change. Previously, Mu et al.[82] had shown that the Cartesian coordinate based PCA would not provide an appropriate FEL as it would not differentiate between internal and overall motions. On the other hand, dPCA considers the dihedral angles of proteins that help to distinguish the motions between the internal and overall dynamic motions and would construct the FEL involving large structural changes. The free energy surface is described with the first two eigenvectors (v1,v2) aswhere P(v1,v2) represents the probability distribution obtained from the MD simulations and Pmax denotes the maximum of distribution, kB is the Boltzmann constant, and T is the temperature (300 K).