Assessment of the binding of Cpd22 to ILK by surface plasmon resonance
Initially, to lay the groundwork for our subsequent simulations, we decided to test the hypothesis of direct ligand binding to ILK. In this regard, to gain both qualitative and quantitative insights, we carried out a surface plasmon resonance (SPR) assay to confirm Cpd22 ligation toward the protein (Fig. 2). Both Cpd22 and recombinant ILK were obtained from commercial sources. Using steady-state analysis of SPR measurements for the interaction between Cpd22 and ILK at room temperature, the dissociation constant (KD) was adjusted by kinetic fitting to 5.4 ± 1.65 µM. This evidence confirmed the ability of this protein to bind the molecule with moderate affinity in vitro. Interestingly, Cpd22 displayed a very similar KD to ATP (3.64 ± 0.49 µM), according to data reported in the literature [30]. This is not surprising as the ILK catalytic core has evolved, losing its catalytic activity, and those changes in the ILK sequence and its associated three-dimensional architecture can result in a reduction in nucleotide binding affinity, as reported for other ATP-binding pseudokianses such as STRADα or TYK2 JH2 [31].
Protein Energy Landscape Exploration Calculations Point To The Nucleotide Cleft
Once direct binding to ILK had been confirmed, we decided to study the molecular recognition phenomenon. As mentioned previously, the lack of information concerning the molecular interaction between Cpd22 and ILK prompted us to study this question using biomolecular modeling techniques. Given the ability to bind ATP, the nucleotide-binding site (herein pseudo-active site) is by far the most promising druggable pocket in ILK, similarly to kinases. As such, and given the reported driven hypothesis of Cpd22 development, our computational studies were conducted using a refined model of the kinase domain of ILK (henceforth ILK) generated from the PDB entry 3KMW [32]. This crystal contains the kinase domain of ILK bound to the N-terminal domain of α-parvin, with the pseudo-active site filled by a non-hydrolyzed ATP molecule and a Mg2+ cation.
Based on the aforementioned observations, the main hypothesis driving the computational work was that Cpd22 could be recognized by the pseudo-active site of ILK. To lay the groundwork for our subsequent calculations, we started searching for potential binding sites, considering the whole protein surface. To this end, we started by setting up several non-biased MD simulations to allow the ligand to diffuse and associate with ILK in the solvent. Thus, Cpd22 was randomly placed in the solvent around ILK and, after three independent replicas of 2 µs, the ligand was found not to reach any stable binding site on the protein (data not shown). Consequently, we decided to employ PELE (Protein Energy Landscape Surface), a relatively novel Monte Carlo algorithm, to investigate Cpd22 binding in a computationally affordable manner [33]. The PELE algorithm has been used successfully to determine the binding mode and location of several enzyme substrates, drugs and protein-protein inhibitors [34, 35]. For PELE calculations, seven independent simulations were carried out starting from four different starting ligand positions in order to explore the ILK surface (Fig.s S1 and S2).
This protocol allowed us to check and study the tendency of Cpd22 to associate with the ATP cleft or other regions on the protein. During these simulations, the ligand was protonated at the distal nitrogen atom of the piperazine ring given the estimated pKa for a realistic physiological environment (pH 7.4). A refined ILK model was employed by extracting the ATP molecule but maintaining the metallic cation and by placing Cpd22 in up to four different locations around ILK (Fig. S1).
Global exploration of the ILK surface revealed various Cpd22 minimal energy locations with strong interaction energies closer to the Mg2+ atom (Fig. 3). Interestingly, these calculations showed that at least one of the seven trajectories generated for each exploration simulation over the ILK surface reached the ATP pocket, and these tended to have the best energy values (Fig. 3 and S2). Furthermore, in the case of starting positions 1 and 3, two independent trajectories found the pseudo-active site of the protein to be most favorable region for Cpd22 binding (see Fig. S2). The remaining simulations did not converge on a consensus or well-defined site over ILK. Moreover, the resulting PELE interaction energies were worse than for the trajectories that reached the ATP cleft between the N- and C-terminal lobes. These findings suggest that Cpd22 tends to bind to a shallow cavity rather than a solvent exposed surface, as expected for a small molecule modulator that targets an enzymatic active site. Visual examination of the final protein-ligand complexes inside the ATP pocket showed that the resulting Cpd22 poses were not similar, thus resulting in different orientations of the molecule inside this cavity (Fig. S2). This is not surprising as the PELE protocol chosen was designed for a global exploration of the protein rather than ligand pose refinement.
As the initial approach was intended to roughly identify a plausible location for Cpd22 binding, more sophisticated methods were needed to propose a plausible and accurate binding mode. In this regard, we first carried out docking calculations considering a rigid receptor and flexible ligand during the conformational sampling step. Glide, which is a broadly validated docking method, in XP mode was selected for docking calculations [36, 37]. Its empirical scoring function was calibrated using kinase inhibitors for the training dataset and, as a result, Glide seems to be a very appropriate docking tool for the purpose of this study with a kinase-like protein such as ILK.
As mentioned above, PELE simulations converged on the ATP binding pocket and adjacent areas, therefore docking calculations were performed in this region of the protein. Initially, the docking protocol was validated by redocking the ATP molecule as positive control. This docking procedure was able to reproduce the binding mode (Fig. S3), and Cpd22 was docked.
Glide produced three different poses, which differed by less than 0.2 RMSD Å, thus resulting in a single well-defined binding mode (Fig. 4). Interestingly, the predicted binding mode was not one of those generated by PELE simulations. However, this was partly expected as the PELE global exploration was not intended to find an accurate binding mode. The Glide XP docking energies ranged from, -6.142 to -5.090 kcal/mol, thus supporting the plausibility of molecular recognition by the ILK pseudo-active site.
The molecule is accommodated inside the ATP pocket, close to the Mg2+ cation. Branch A of Cpd22 is oriented towards the outer region of the cleft, with the piperazine ring pointing towards Glu283. This facilitates an ionic interaction reinforced by a doubled hydrogen bond between the positively charged nitrogen atom of Cpd22 and the carboxylate group of Glu283 (Fig. 4). Given its strength, this is a key anchoring interaction and its presence also explains the experimental loss of activity when the piperazine ring is replaced by β-alanine or other heterocycles such as morpholine [28].
The second branch of the pyrazole ring, which corresponds to the N-methyl propanamide chain (branch B), is located in the region neighboring Mg2+. The oxygen atom from the carbonyl group is hydrogen bonded to Lys202 and might chelate the metal thanks to the lone electron pairs. A long-distance hydrogen bond (~ 4 Å) is observed between the Ala338 carbonyl and Cpd22 -NH amide groups. With regard to the diaryl moiety positioned at C5 of the pyrazole ring (branch C), docking simulations show that it is oriented towards a hydrophobic tunnel delimited by the side chains of amino acids Leu199, Leu207, Trp271 and Pro273. A detailed visualization of the solvent-accessible surfaces revealed that this is a narrow space, thus hindering the accommodation of bulky groups such as the phenanthrene ring, as reported by Liu and collaborators [28].
Molecular dynamics simulation studies
Based upon the results discussed above, our findings point to the pseudo-active site of ILK as being the target cavity for Cpd22 molecular recognition. Molecular simulations, especially all-atom molecular dynamics (MD) simulations, are an essential tool for the study of biomolecules in a time-dependent manner. In this regard, classical MD can provide a realistic picture, at an atomic level, of how molecular recognition between a protein and a ligand occurs. This permits a better understanding of the interaction between Cpd22 and ILK and its implication for the stability and conformational space explored by the protein-ligand complex.
Initially, we aimed to validate the binding mode generated by way of docking calculations. To this end, three replicas of a 1 µs classical MD with explicit water molecules were performed (ILK-Cpd22). To monitor ligand binding at the proposed site, we measured the root mean square deviation (RMSD) corresponding to the heavy atoms of Cpd22 as the observable variable, combined with a visual inspection of the protein-ligand trajectories. Furthermore, to compare and investigate the solution structure and flexibility of ILK in the absence of any ligand, an additional three MD simulations of 1 µs were carried out in explicit water with the protein alone (Apo-ILK).
The simulated systems were considered equilibrated beyond 100 ns, as shown by the RMSD convergence (Fig. S4). Furthermore, the protein remained stable during simulations with mean backbone RMSD values of approximately 2 Å. It is noteworthy that ILK-Cpd22 replicas presented lower RMSD values, thus suggesting stable and less flexible systems compared to the ILK kinase domain alone in solution. This observation seems to be related to the effect that Cpd22 can exert on the protein upon ligation.
With regard to Cpd22, during our MD simulations, only one metastable binding mode was observed across the three replicas. This is reflected in terms of low fluctuating RMSD values during the simulation, as shown in Fig.s 5A and S4. Furthermore, ligand conformational changes are also associated with a slight decrease in the solvent accessible surface area (SASA) values during the majority of the simulation, which reflects the tendency of the ligand to bury in the ATP pocket. Cpd22 only started to abandon the active site during the third replica, after 800 ns, probably initiating an unbinding event that was not completed in the remaining simulation time.
Moreover, a detailed visual examination of the MD simulation allowed us to observe how the initial binding mode (from docking calculations) suffered subtle but significant shifts that altered some of the initial interaction patterns. Thus, a few nanoseconds after the production phase, Cpd22 is accommodated inside the ATP pocket. As result, branch B orients its carbonyl group towards the Mg2+, thus chelating the cation during the rest of the simulation by displacing a solvent water molecule. In the original crystal structure (PDB: 3KMW) [32], Mg2+ is coordinated by the carboxylate functionality of Asp339 and the coordination sphere is completed by three oxygens atoms from the α-, β- and γ-phosphate groups of ATP, plus two water molecules.
Interestingly, in our simulations for ILK bound to Cpd22, the carbonyl group from ligand branch B penetrates into the Mg2+ coordination sphere by way of the oxygen atom (substituting the role of ATP γ-phosphate). Therefore, in the ILK-ATP and apo-ILK trajectories, it can be seen that the Asp339 carboxylate and Arg323 carbonyl participates in the coordination with the cation. Finally, four water molecules for apo-ILK and ILK-Cpd22 and only one molecule in the case of ILK-ATP are in charge of maintaining a quasi-octahedral coordination geometry across the simulations (Fig. 5B). As such, our results also support the role of Cpd22 mimicry of ATP to chelate Mg2+, thereby preventing its loss from the pseudo-active site [38].
This ligand rearrangement is also accompanied by a partial reorientation of the biphenyl moiety in the hydrophobic pocket to an outer region. This sub-pocket, which was previously occupied by branch C, corresponds to the location where the adenine ring of the ATP molecule is lodged (Fig. S5) and is a critical optimization region for kinase inhibitors. In light of these findings, it appears that the incorporation of a scaffold on Cpd22 that can mimic hydrogen bonding features will be beneficial for ligand enthalpic optimization [39, 40], a strategy that was not fully covered by the inhibitors used in the original study [28]. Moreover, the pyrazole ring is partially shifted, thus meaning that it can occasionally establish a new hydrogen bond between the pyrazole N2 and Arg323 sidechain (mean 10.58% lifetime occupancy across the three replicas).
Overall, these observations seem to indicate that the Cpd22 binding mode is mainly governed by the strong ionic interaction between the piperazine ring and Glu283 (which does not suffer significant changes). Additionally, this binding mode is also stabilized by the stable coordination bond with the Mg2+ cation, which also remained quite stable during MD simulations, with a mean distance between the oxygen and magnesium atom pairs of around 2 Å across replicas (Fig. S6).
To further understand and support this hypothetical binding mode, a simulated annealing protocol was carried out to characterize this protein-ligand conformation as a local or even global minimum in the conformational landscape. Simulated annealing techniques have previously been employed successfully to optimize structures from experimental assays [41, 42], characterize binding modes [43], and predict protein conformations [44, 45]. Using this protocol, it was discovered that the Arg323 sidechain changes its orientation towards the phenyl ring present in Cpd22 (branch A) for all ILK-Cpd22 complexes (Fig. 5C). As such, the cationic amino acid receives π density from the phenyl ring, thus contributing to Cpd22 anchoring. In addition, the N-methyl substituent (branch B) is able to establish a CH···π interaction with the Tyr351 aromatic ring. This finding may explain the experimental evidence that replacement of a methyl by an ethyl group preserves the desired activity [28].
Cpd22 binding induces changes in ILK structural plasticity
In light of the findings discussed above, we then studied whether Cpd22 binding can induce conformational changes to the ILK kinase domain. To this end, an extra complex corresponding to ILK bound to its cognate ligand ATP (ILK-ATP) was simulated in triplicate. This simulation also displayed a stable profile in terms of RMSD, especially for ATP (Fig. S4), the hydrogen bonding network for which anchors the molecule strongly during simulations. When comparing stability between ATP and Cpd22 molecules, we found that the former binds more tightly to ILK in terms of mean RMSD (1.64 ± 0.09 Å) than the synthetic ligand (1.95 ± 0.28 Å) (Fig. 6A). This is not surprising as the pseudo-active site conserves most of the molecular determinants to bind nucleotides, with ATP being the natural ligand of this pseudokinase [30].
We then focused on the influence that ligand binding could have on the molecular rigidity of ILK. An exhaustive RMSD analysis led us to conclude that the overall flexibility of the protein is significantly higher in the absence of any ligand (apo-ILK) compared with both holo proteins (Fig. 6B). In terms of RMSD, apo-ILK exhibited a mean value of 2.08 ± 0.25 Å (SD) while ILK-ATP showed a slightly lower value (1.96 ± 0.1 Å). This is in good agreement with recent studies relating ATP binding to a conformational stabilization of the ILK kinase domain, thus also validating our approach [46]. Moreover, our observation is in full agreement with experimental studies, which observed stabilization induced by ATP binding rather than gross conformational changes [7, 30, 46].
A similar picture was observed when comparing ILK-Cpd22, which resulted in an even more pronounced reduction in the mean RMSD value (1.86 ± 0.21 Å) when compared to ILK-apo and ILK-ATP. Indeed, this was also corroborated by the change in protein radius of gyration, the mean values of which decreased across the holo systems (see supporting info Table S1), thus indicating that ILK stabilization by pseudo-active-site ligands is accompanied by a slight increment in protein compactness. This effect may be of interest from the point of view of ILK protein-protein interactions with PINCH and parvin partner adaptor proteins, as the response elicited by Cpd22 on ILK flexibility as a whole could affect PPI assembly and/or stability given the lower conformational motions. This hypothesis is supported by the previously observed effect of ATP binding, which stabilizes the pseudokinase domain, thus facilitating a more effective interaction between the focal adhesion proteins PINCH and parvin [38, 46].
To study the flexibility in terms of single residues, RMSF values per Cα were calculated (supporting info Fig. S7) and inspected visually in the context of protein dynamics. Fluctuations in the nucleotide-binding cleft residues were observed in the case of ILK-Cpd22 simulations, as would be expected due to ligand accommodation inside the pocket, along with branch C shifting alongside the hinge region. Remarkably, we observed differences in the pseudo-active site loop, spanning residues 339 to 359. Thus, while for apo-ILK this remains with low mobility (according to RMSF values, Fig. S8), ATP induced a significant increase (up to approx. 6 Å). This is in contrast with the findings of Fukuda et al. for this loop, who suggested its rigidity based on the relatively low b-factors for the crystal structure (PDB ID: 3KMW) compared to the rest of the protein [32]. Interestingly, for ILK-Cpd22 simulations, higher RMSF values were also observed in the segment spanning residues His309 to Asn321, which is located in the C-lobe subdomain and flanking the pseudo-active binding site. Another less intense peak, corresponding to the glycine-rich P loop (residues 202–205), fluctuates to accommodate ATP and Cpd22 inside the cleft with values close to 2 Å (Fig. S6).
Analysis Of End-point Binding Free Energy
To gain further support for our binding mode hypothesis and further insights into binding thermodynamics, the approximate free energy of binding was calculated using the MM-ISMSA method [47]. In addition, this method was also employed for the ATP molecule in order to have a positive control as reference (Supporting Table S2). The ligand binding energy calculations obtained using this method are easily reproducible and computationally inexpensive compared to other more sophisticated methods, and provide good accuracy [43, 47–49].
The free energies averaged across replicas revealed how both ligands bind with affinity for the ATP pocket in ILK with mean values of -50.49 and − 48.93 kcal/mol for ATP and Cpd22 respectively. Thus, the predicted interaction affinity for ATP towards ILK is slightly more favorable than for Cpd22 (Table S2), which is expected as the former is the cognate ligand of this pseudo-kinase. However, the small difference between both affinities, in terms of binding free energy, is consistent with the experimentally measured values of KD, which also exhibited quite similar values for both molecules, as mentioned previously (see “Assessment of binding free energy of Cpd22 by surface plasmon resonance” section).
Indeed, whereas the main favorable contribution for ATP binding energy is provided by electrostatic interactions (coulombic term), in the case of Cpd22 it is given by van der Waals forces (12–6 Lennard–Jones potential) (Fig. 7A). This can be explained by the larger number of polar substituents present in the ATP molecule, especially phosphate groups. These moieties participate actively in the interaction with several polar and charged residues of the binding pocket, including Ser204, Lys220, Asn279, Arg323, Asp339, and Lys341, thereby exerting a strong influence on the binding pose and its stability during MD simulations (Fig. S5B). A detailed analysis of the per-residue contribution to the binding free energy for CPD22 may reveal the importance of some residues and their role in anchoring the molecule into the pseudo-active site cleft (Fig. 7B and S7). As observed previously in docking and MD simulations, the largest contribution comes from Glu283 (approx. -10 kcal/mol averaged over three independent replicas) and its ionic interactions with the
piperazine moiety (Fig. 7B). Interestingly, Arg323 is also considered to be an important residue for the binding free energy according to these calculations, although it depends on the side-chain orientation. During MD simulations, Arg323 can flip and accommodate its side-chain to establish a cation···π interaction between the guanidinium group and the ligand phenyl ring from branch A. Furthermore, Arg323 also interacts with this branch of Cpd22 via van der Waals interactions with the CH2 groups in its side chain, thus making this residue an important anchoring partner and stabilizing branch A.
Further contributions are provided by aliphatic and aromatic amino acids such as Leu199, Leu207, Tyr278, Trp271 and Tyr351, mainly as a result of non-polar contacts with Cpd22. It is worth noting that the stable coordination observed during MD simulations influences the binding free energy via the Mg2+ cation contribution (approx. 2 kcal/mol across replicas) to Cpd22 binding.