Pharmacoinformatics and molecular dynamic simulation studies to identify potential small-molecule inhibitors of WNK-SPAK/OSR1 signaling that mimic the RFQV motifs of WNK kinases

2.1 Generation and validation of pharmacophore model

The crystal structure of C-terminal domain (CTD) of OSR1 in complex with RFQV (Arg-Phe-Asn-Val) peptide derived from WNK4 (PDB: 2V3S) was imported into LigandScout software from RCSB protein data bank (Villa et al., 2007). The structure-based pharmacophore was generated using automatic pharmacophore generating tool in LigandScout program (Wolber and Langer, 2005). The resulted pharmacophore model consists of the whole features involved in the binding of RFQV peptide residues to the primary pocket of OSR1-CTD. STOCK1S-50699, a known WNK and SPAK binding inhibitor, was docked (using Autodock vina) and mapped on the generated pharmacophore model using LigandScout program to obtain the final pharmacophore model (Mori et al., 2013). The final pharmacophore model was then validated using the receiver operating characteristic (ROC) curve, with LigandScout software, by screening the pharmacophore model against a set of active and inactive compounds to determine the ability of this pharmacophore to distinguish between these compounds. STOCK1S-50699 and STOCK2S-26016, known WNK-SPAK binding inhibitors, were used as active compounds. The two compounds were also used to generate the decoy set of 100 inactive compounds (50 compounds per each) using DUD-E webserver (http://dude.docking.org/generate) (Mysinger et al., 2012).

2.2 Pharmacophore-based virtual screening

In silico pharmacophore-based virtual screening was performed with “Asinex focused protein–protein interaction (PPI)” library having 11,870 small-molecules against the generated pharmacophore model using LignadScout software. The library contains non-macrocyclic compounds with a diversity of more than 500 scaffolds. The library was obtained from (https://www.asinex.com/ppi/) in sdf format and was converted into Idb using LigandScout library generation tool. The compounds that meet all pharmacophore features were considered as hit compounds and ranked based on their pharmacophore-fit scores which reflect to which degree the molecules fit the pharmacophore features.

2.3 Docking-based virtual screening

The retrieved hit compounds from pervious screening was subjected to docking-based virtual screening against the 3D structure of OSR1-CTD (PDB: 2V3S) using Autodcok Vina in PyRx 0.8 program (Dallakyan and Olson, 2015). Before docking, hit compounds were energy minimized and converted from sdf files into pdbqt files using Open Babel tool in PyRx 0.8 program (O'Boyle et al., 2011). The grid box was cantered to cover the amino acid residues involved in the topology of the primary pocket of OSR1-CTD. Prior screening, STOCK1S-50699 was added to the database as a control. Compounds that bind to OSR1-CTD with high binding affinities in comparison to STOCK1S-50699 were considered for further analysis.

2.4 In silico ADME-T analysis

pkCSM server was used to evaluate the absorption, distribution, metabolism and excretion- toxicity (ADME-T) parameters for the identified hit compounds (Pires et al., 2015). For the compound to be selected as a hit, it must be non-hepatotoxic and non-carcinogenic. SwissADME was used to assess other physiochemical properties of these hit compounds (Daina et al., 2017).

2.5 Calculation of ligands efficiency metrics and inhibition constants

The inhibition constants (K_i) of hit compounds were predicted from Autodock vina binding energy scores using Eq. (1) (Edwards and Price, 2010; Hopkins et al., 2004; Hopkins et al., 2014; Murray et al., 2014; Reynolds et al., 2007).

The ligand efficiency parameters were estimated using the following equations: (Edwards and Price, 2010; Hopkins et al., 2004, 2014; Murray et al., 2014; Reynolds et al., 2007).

In which LE, LE_scale, LLE, FQ and LELP are stand for Ligand Efficiency, Ligand Lipophilic Efficiency, Ligand Efficiency Scaled, Fit Quality and Ligand Efficiency Lipophilic Price, respectively.

2.6 Pharmacophore mapping of hit compounds

Hit compounds were mapped into the generated pharmacophore model using the alignment tap in LigandScout software.

2.7 Molecular docking study

Hit compounds fulfilling the pervious filters were docked against the 3D structure of OSR1-CTD using Autodock vina program (Trott and Olson, 2010). The protein structure (PDB: 2V3S) was obtained from RCSB protein data bank (Villa et al., 2007). The structure was solved at an X-ray resolution of 1.7 Å and it was composed of dimers of OSR1-CTD in complex with RFQV peptide. Discovery studio 4.5 (Accelrys, San Diego, CA, USA) was used to remove the unwanted water molecules and ligands as well as to generate the pdb files for the protein in monomer form. Autodock tools program was used to generate the pdbqt files and to prepare the gridbox for the docking configuration files (Sanner, 1999). The gridbox was centered to cover the primary pocket with the following parameters; the box size of x = 14 y = 14 z = 22 and the box center: x = 1.605 y = 11.139 z = 23.381. Discovery Studio 4.5 and PyMOL Molecular Graphics System 1.3 were used to visualized and analyzed the docking results.

2.8 Molecular dynamic (MD) simulation

The dynamic behavior of docked inhibitor-OSR1-CTD complexes was evaluated via all-atom MD simulation for 20 ns using GROMACS 2018.1 package (Hess et al., 2008). The topology files of all docked-inhibitors were obtained using SwissParem tool (Zoete et al., 2011). The OPLS-AA/L force field was applied to the system to carry out the MD simulation. A triclinic water box of TIP3P water model molecules (Jorgensen et al., 1983) with 1.0 nm distance from the edge of the box to protein was surrounded to each complex of protein-ligand system. A suitable numbers of counter ions were added to neutralize the system. The system was then equilibrated and energy minimized using steepest decent algorithm with tolerance value of 1000 kJ mol⁻¹ nm⁻¹ followed by equilibration using NVT and NPT ensemble for 100 ps. Bond lengths and electrostatic calculations were constrained using LINear Constraint Solver (LINCS) algorithm and particle mesh Ewald method, respectively (Essmann et al., 1995; Hess et al., 1997). MD simulation was carried out for 20 ns MD production with a time step of 2 fs (femto-second) at the constant pressure of 1 atm and constant temperature of 300 K and snapshots saves every 2 pico-second (ps). Several parameters included root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF) and radius of gyration (Rg) were analyzed using GROMACS to determine the conformational and performance stability of each complex system in the dynamic environment.

3.1 Pharmacophore model generation and validation

A structure-based pharmacophore model was generated based on the crystal structure of the RFQV (Arg-Phe-Asn-Val) peptide derived from WNK4 binding to the C-terminal domain (CTD) of OSR1 (Fig. 2) (Villa et al., 2007).

Close

Fig. 2 Crystal structure of OSR1-CTD (PDB: 2v3s): (A) Ribbon representation of OSR1-CTD in complex with RFQV peptide-derived from WNK4 (yellow). α- helices and β- sheets were shown in cyan and purple colors, respectively. (B) Molecular surface representation of OSR1-CTD (white) shown the binding mode of RFQV peptide-derived from WNK4 (yellow) to the primary pocket.

Export to PPT

Based on the interaction of RFQV motif with the primary pocket of OSR1-CTD, the extracted pharmacophore features consist of eleven hydrogen bond acceptors (HBA), ten hydrogen bond donors (HBD), two hydrophobic (HYD) and one positive ionizable features beside twenty-eight exclusion-volume spheres which are the essential regions that determine the overall shape of the binding pocket (Figs. 3A and 1B).

Close

Fig. 3 Generation of pharmacophore model. Pharmacophore model was generated based on the interaction of RFQV peptide derived from WNK4 kinase with the CTD of OSR1 kinase (PDB:2V3S). (A) Interaction of RFQV peptide with OSR1-CTD. (B) Structure-based Pharmacophore model that was generated using LigandScout program. (C) Interaction of STOCK1S-50699, WNK-SPAK binding inhibitor, with OSR1 C-terminal domain. (D) Mapping the RFQV (cyan) and STOCK1S-50699 (pink) onto the generated structure-based pharmacophore model. The HBA, HBD, H, and PI represent hydrogen bond acceptor, hydrogen bond donor, hydrophobic and pi-pi interaction, respectively. The pharmacophore features were defined in LigandScout by colour codes; red, green, yellow, blue and grey spheres which represent hydrogen bond acceptor, hydrogen bond donor, hydrophobic, positive ionizable group and exclusion volume, respectively.

Export to PPT

Since the using of the peptide-based pharmacophore model for screening of small-molecules library may result in identification of no hit compounds as the pharmacophore-based screening becomes inefficient with more than eight features, a rational method to reduce the number of pharmacophore features is needed (Jung et al., 2018). Therefore, a known WNK and SPAK binding inhibitor with a K_d value of 37 µM, namely STOCK1S-50699, was docked and mapped onto the generated pharmacophore model to identify the most important interaction features for optimal binding to OSR1-CTD (Figs. 3C and 1D) (Mori et al., 2013). Five pharmacophore features were identified; three hydrogen bond acceptors and two hydrophobic regions. This pharmacophore model beside all of the exclusion-volumes were considered as the final structure-based pharmacophore model (Fig. 4A). Notably, the distances between the pharmacophore features are large which reflect the size of typical protein-protein interaction binding pockets (Voet et al., 2013). To estimate the performance of pharmacophore model, the pharmacophore model was screened against a set of active and decoy compounds to determine its ability to correctly recognized a list of compounds as actives or inactive (decoys). The ROC analysis which is indicated by the area under the curve (AUC) as well as enrichment factor (EF) values showed that the pharmacophore yielded a ROC score of 0.86 which means that a randomly-selected-active compound has a higher score than a randomly-selected-decoy 8.6 times out of 10 (Fig. 4B). Since the value of AUC was beyond 0.5, the pharmacophore performed good in distinguishing between active and inactive samples within the screened dataset. The EF values of the set screened by the pharmacophore at 1%, 5%, 10%, and 100%, were 0.0, 10.2, 5.1, and 2.4, respectively. The pharmacophore retrieved the two active (100%) compounds from screen dataset, corroborating the ROC statistics.

Close

Fig. 4 The final pharmacophore model and its performance. (A) The final pharmacophore model. The pharmacophore features were shown in LigandScout by colour codes; red and yellow spheres which represent hydrogen bond acceptor and hydrophobic, respectively. (B) The performance of the pharmacophore model by ROC curves using the Directory of Useful Decoys (DUD) dataset. The ROC plot was generated using LigandScout.

Export to PPT

3.2 Combined virtual screening of protein–protein interaction (PPI) library

The validated structure-based pharmacophore model was used as a 3D query for the pharmacophore-based virtual screening of “Asinex focused protein-protein interaction library” having 11,870 small-molecules. Of these compounds, 150 hit compounds were found to meet the pharmacophore query features with a hit rate of 1.3%. The pharmacophore-fit scores of these compounds was in between 54.85 and 58.31. In the next step, these hit compounds were subjected to docking-based virtual screening against the 3D structure of OSR1-CTD to filter them further based on the free energy binding score. From this exercise, 31 hit compounds were identified to bind to the OSR1-CTD with high docking scores in comparison to the binding score of STOCK1S-50699, a known SPAK/OSR1 inhibitor, which was used as a control. The docking scores for these compounds were between −1.8 and −9 Kcal/mol.

3.3 ADME-T properties

The determination of the absorption, distribution, metabolism and excretion-toxicity) (ADME-T) parameters is a significant step in the early phases of drug discovery. Computational tools have provided useful and efficient measurements of these ADMET parameters in time- and cost effective manners (Hughes et al., 2011). The key properties were measured for the 31 compounds. Among the 31 compounds, only non-hepatotoxic and non-carcinogenic compounds were selected and three hit compounds were obtained (Fig. 5). The summary of their ADME-T properties is depicted in Table 1. Importantly, hit 2 has the probability of crossing the blood brain barrier (BBB+) while the other two hits (hit 1 and 3) were not. It was also observed that hit 2 has the highest probability of being absorbed by human intestine (94.219%), then hit 2 (89.239%) and then hit 3 (85.059%). The aqueous solubility is critical property for drug oral activity as well as for pharmaceutical preparation. The values of aqueous solubility for the identified hit compounds were within the standard range which should be between 1 and 5 (Tsaioun and Kates, 2011). The bioavailability scores for all hit compounds were 0.55 meaning that the compounds may have >10% bioavailability in rat (Martin, 2005). Therefore, these values indicated that the hit compounds may have good absorption and distribution properties. The cytochrome P450 (CYP2D6) is a key enzyme responsible for metabolism of >25% of currently available drugs (Pirmohamed and Park, 2003). Interestingly, the hit compounds were predicted to have no inhibitory effect on this enzyme. The predicted LD50 values for the hit compounds were 2.434, 2.666 and 3.595 mol/kg for hit1, 2 and 3, respectively. The LD50 values were expressed in mol/kg according to the standard practice of QSAR in which each 1/mol/kg is corresponding to ∼500 mg/kg (Raevsky et al., 2018). Accordingly, all hit compounds fall into class II labelled as ‘moderately toxic’ as they have LD50 values of ∼1200, 1300 and 1800 mg/kg for hit1, 2 and 3, respectively.

Close

Fig. 5 Chemical structures of candidate compounds. A, B and C are hit 1, 2 and 3 respectively.

Export to PPT

3.4 Physicochemical properties measurement and bioactivity prediction of the hit compounds

The physicochemical properties provide a deeper insight into the drug-likeness of a drug molecule. Lipinski’s rule of five is a famous method to evaluate the drugability of compounds which states: For a given molecule, the number of hydrogen-bond donors (HBD) and hydrogen-bond acceptors (HBA) must not be greater than 5 and 10, respectively, the molecular mass and logP should not exceed 500 g/mol and 5, respectively (Lipinski, 2016). The summary of the physicochemical properties of the identified hit compounds is illustrated in Table 2. The results indicated that all hit compounds obey the Lipinski’s rule of five. The number of rotatable bonds is known to modulate the bioavailability of compounds and it should not be more than 7 (Veber et al., 2002). All three hit compounds met this standard except hit 3 that has a value of 9 for number of rotatable bonds. The range for refractivity is another parameter of drug likeness which recommended to be between 40 and 130 (Ghose et al., 1999). None of the compounds violate this limit except hit 2 that has slightly high value of 132.72. Moreover, the three hit compounds have ideal polar surface area (PSA) values which is recommended to be less than 140 for optimal drug absorption and distribution (Cerqueira et al., 2015).

The binding energy scores of the three hit compounds were used to calculate the inhibitory constant (K_i) in micromolar concentration range (µM) Table 3. The K_i value determines the activity of compounds, typically it should be in a micromolar concentration range for a lead compound as well as in a low nanomolar concentration range for a drug (Hughes et al., 2011; Stevens, 2014). The calculated K_i values for the three compounds are 24.43-, 24.43- and 56.75 µM, respectively. Therefore, these compounds can be defined as lead compounds for discovery of WNK-SPAK/OSR1 signaling inhibitors.

In term of ligand efficiency metrics, a qualified hit should possess a threshold value of 0.3, 3, and 0.8 for the ligand Efficiency (LE), ligand lipophilic efficiency (LLE), and fit quality (FQ), respectively (Hopkins et al., 2004; Hopkins et al., 2014; Murray et al., 2014). Due to the nature of the primary pocket which has been showing to be a surface exposure binding site, the in silico binding energy of ligands are expected to be low (Villa et al., 2007). Consequently, the value of LE for the identified three hit compounds were lower than the recommended limit Table 3. This phenomenon was also observed with STOCK1S-50699 which has LE = 0.19. The FQ value is directly affected by LE value which results in lower value of FQ for the three hit compounds (Hopkins et al., 2004). In another hand, the values of ligand lipophilic efficiency (LLE) for the hit compounds were within the standard range (Hopkins et al., 2004). Furthermore, the Ligand Efficiency Lipophilic Price (LELP) measures the ligand efficiency in term of lipophilicity of compounds and it should be between −10 and 10 for a given lead compound (Hopkins et al., 2004). Obviously, the three hit compounds showed higher values of LELP due to the lower values of LE. The LELP value was also high for STOCK1S-50699 (LELP = ∼42).

3.5 Pharmacophore mapping

The hit compounds remarkably mapped well onto all the features of the pharmacophore model (Fig. 6). Table 3 showed the pharmacophore fit scores which represents how well a compound maps to the pharmacophore; the higher fit score indicates a better fit to the pharmacophore model and the molecules with high values should be active as WNK and SPAK/OSR1 binding inhibitors.

Close

Fig. 6 Mapping the three hit compounds onto the structure-based pharmacophore model. A, B and C are for hit 1, 2 and 3, respectively. Pharmacophore features were represented by colour codes; red and yellow spheres which indicated hydrogen bond acceptor and hydrophobic, respectively.

Export to PPT

3.6 Molecular interactions and binding modes

To analyze the binding modes as well as the type of interactions of the hit compounds with the CTD of OSR1, a molecular docking was performed using Autodock Vina (Trott and Olson, 2010). Hit 1 exhibits conventional and carbon hydrogen bonds with Arg451 and Ile450, respectively, Pi-alky bonds with Leu473, Ala 471 and Val 464 and Pi-anion interaction with Glu467 (Fig. 7A). The molecular interactions of hit 2 include conventional and carbon hydrogen bonds with Glu453 and Phe542, respectively and Pi-alky interactions with Leu473 and Ile450 (Fig. 7B). Hit 3 involves in conventional hydrogen bond with Arg451 and Pi-alky bonds with Ala471, Leu468 and Ile450 (Fig. 7C). Interestingly, all hit compounds adapt similar binding mode in the primary pocket of OSR1-CTD (Fig. 7D–F). Mutation and NMR binding studies, indicated that most of the residues involve in the interactions with the identified hits such as Leu473, Arg451, Ile450, Ala471 have been shown to be essential for the binding of RFQV peptide to the CTD of OSR1 (AlAmri et al., 2019, 2017; Villa et al., 2007; Vitari et al., 2006).

Close

Fig. 7 The molecular interactions and binding modes of hit compounds with the CTD of OSR1. Ribbon representation of the interaction of (A) hit1 (cyan) (B) hit 2 (pink), (C) hit 3 (blue) with OSR1-CTD. The type of interactions was illustrated in green, light green, pink and yellow which represent conventional hydrogen bond, carbon hydrogen bond, hydrophobic (Pi-Alkyl) and Pi-Anion type of interactions. Molecular surface representation of the binding of (D) hit 1 (cyan) (E) hit 2 (pink), (F) hit 3 (blue) to the primary pocket of OSR1-CTD. The co-crystal RFQV peptide-derived from WNK4 is shown in green.

Export to PPT

3.7 MD simulation

The dynamic stability and behavior of each docked-inhibitor-OSR1-CTD complex was explored through a 20 ns molecular dynamic simulation study. The root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF) and radius of gyration (Rg) were calculated for the protein backbone. The results of calculated RMSD of docked complexes showed that after sharp raised in the RMSD values at the beginning of the MD simulation all system reach equilibrium after ∼1.5 ns for hit 1 and 2 and after ∼2.5 ns for hit 3 with average RMSD values of 0.13 ± 0.02, 0.14 ± 0.02 and 0.14 ± 0.02 Å for hit 1, 2 and 3, respectively. (Fig. 8).

Close

Fig. 8 RMSD (Å) vs time (ns) of OSR1-CTD backbone obtained from complexes of OSR1-CTD-screened inhibitors.

Export to PPT

The RMSD results were all lower than 0.2 Å indicating a stable dynamic behavior for the last 18.5 ns for hit 1 and 2 and 17.5 ns for hit 3. Due to the critical role of individual amino acid in the stability of ligand inside the binding pocket, the RMSF value of was calculated to explore the flexibility of each residue (Fig. 9). The resultant RMSF values showed no significant fluctuations were observed at the ligand binding sites in the docked-inhibitor-OSR1-CTD complexes with fluctuations values ranged from 0.05 to 0.15 Å. The high RMSF peaks were observed in the loop regions formed by residue Glu507, Gly508, Ser509, Asp510 and Ile511. Notably, the latter effect was the least with compound hit 3 with RMSF value of 0.15 compared with RMSF values of 0.30 and 0.25 for hit 1 and 2, respectively (Fig. 9).

Close

Fig. 9 RMSF (Å) vs residue number of OSR1-CTD when bound to final screened inhibitors.

Export to PPT

To evaluate the structural stability and compactness of protein, the radius of gyration (Rg) was calculated (Fig. 10). The Rg of protein backbone for the docked-inhibitor-OSR1-CTD complexes showed stable behavior with values between 1.25 and 1.33 Å throughout the simulation at 20,000 ps. However, compound hit 2 showed sharp fluctuation in the Rg values after 10,000 ns and then remained stable through all the rest of simulation period. This could suggest that the hit 2 may adapt a new conformation within the binding pocket. The results of MD simulation indicate that the docked-inhibitor-OSR1-CTD complexes remained stable with favorable conformations throughout 20 ns suggesting that the identify inhibitors were stable at the active site of OSR1-CTD during the interactions.

Close

Fig. 10 Radius of gyration (Å) vs time (ns) obtained from complexes of OSR1-CTD-screened inhibitors.

Export to PPT