Screening and Analysis of Possible Drugs Binding to PDGFRα: A Molecular Modeling Study

The platelet-derived growth factor receptor (PDGFR) is a membrane tyrosine kinase receptor involved in several metabolic pathways, not only physiological but also pathological, as in tumor progression, immune-mediated diseases, and viral diseases. Considering this macromolecule as a druggable target for modulation/inhibition of these conditions, the aim of this work was to find new ligands or new information to design novel effective drugs. We performed an initial interaction screening with the human intracellular PDGFRα of about 7200 drugs and natural compounds contained in 5 independent databases/libraries implemented in the MTiOpenScreen web server. After the selection of 27 compounds, a structural analysis of the obtained complexes was performed. Three-dimensional quantitative structure–activity relationship (3D-QSAR) and absorption, distribution, metabolism, excretion, and toxicity (ADMET) analyses were also performed to understand the physicochemical properties of identified compounds to increase affinity and selectivity for PDGFRα. Among these 27 compounds, the drugs Bafetinib, Radotinib, Flumatinib, and Imatinib showed higher affinity for this tyrosine kinase receptor, lying in the nanomolar order, while the natural products included in this group, such as curcumin, luteolin, and epigallocatechin gallate (EGCG), showed sub-micromolar affinities. Although experimental studies are mandatory to fully understand the mechanisms behind PDGFRα inhibitors, the structural information obtained through this study could provide useful insight into the future development of more effective and targeted treatments for PDGFRα-related diseases, such as cancer and fibrosis.


Introduction
Platelet-derived growth factor receptors alpha and beta (PDGFRα and PDGFRβ) are membrane receptors that play a key role in a variety of diseases, such as cancer [1][2][3][4][5], immune-mediated pathologies, such as systemic sclerosis (SSc) [6][7][8][9][10], and viral infection [11], by exerting function on the transduction of extracellular signals into the cell. PDGFRs are transmembrane glycoprotein dimer molecules consisting of an extracellular ligand-binding region divided into five Ig-like domains bound to an intracellular receptor tyrosine kinase (RTK) domain through a single transmembrane alpha helix. The tyrosine kinase domain is a typical type III RTK formed by a bilobal kinase, a juxtamembrane (JM) domain, and an activation loop (A-loop) [11,12]. The JM domain is inserted into the active site, while the A-loop (that starts with a conserved amino acid sequence DFG) controls access to it and partially blocks the ATP and substrate site in the inactive state of the enzyme. In the active state of the RTK, the DFG motif is located in a hydrophobic pocket close to the ATP binding site. When an inhibitor binds this pocket, the ATP binding site undergoes a conformational change that reversibly or irreversibly prevents its normal function [13]. Considering that the human PDGFRα has become one of the most important therapeutic targets for the above-mentioned diseases [4][5][6][7][8][14][15][16][17], its modulation/inhibition is a critical aspect that determines the investigation of new ligands or new information to design new effective drugs. As previously reported [18][19][20][21][22][23], inhibiting PDGFRα (overexpressed and/or abnormally activated) can reduce many types of cancer, including prostate, ovarian, breast, pancreatic, and liver cancers, by inactivating downstream signaling pathways that regulate cell proliferation, migration, and angiogenesis [24][25][26]. On the other hand, these signal pathways can induce proliferation and chemotaxis of myofibroblasts, collagen production, and adhesion in endothelial cells [27]. In pathological states, which are also correlated to viral infections, they are upregulated in patients with a progressive fibrotic disease such as SSc [28][29][30].
For these reasons, the structure-based drug design (SBDD), which promotes the development of novel drugs with a potential affinity for therapeutic targets, has turned out to be an essential tool for rapid and cost-efficient lead discovery. SBDD is a growing, iterative, and powerful approach that involves the structural evaluation of targets in the drug discovery process. It has the ability to reduce the time and cost of developing new drug lead molecules with potential therapeutic effects [31]. A virtual ligands database screening, followed by molecular docking, molecular dynamics (MD), and three-dimensional quantitative structure-activity relationship (3D-QSAR) analysis, as previously reported [2,3,5,32], provides an excellent workflow ( Figure 1) to assess an efficient SBDD.
close to the ATP binding site. When an inhibitor binds this pocket, the ATP binding site undergoes a conformational change that reversibly or irreversibly prevents its normal function [13].
Considering that the human PDGFRα has become one of the most important therapeutic targets for the above-mentioned diseases [4][5][6][7][8][14][15][16][17], its modulation/inhibition is a critical aspect that determines the investigation of new ligands or new information to design new effective drugs. As previously reported [18][19][20][21][22][23], inhibiting PDGFRα (overexpressed and/or abnormally activated) can reduce many types of cancer, including prostate, ovarian, breast, pancreatic, and liver cancers, by inactivating downstream signaling pathways that regulate cell proliferation, migration, and angiogenesis [24][25][26]. On the other hand, these signal pathways can induce proliferation and chemotaxis of myofibroblasts, collagen production, and adhesion in endothelial cells [27]. In pathological states, which are also correlated to viral infections, they are upregulated in patients with a progressive fibrotic disease such as SSc [28][29][30].
For these reasons, the structure-based drug design (SBDD), which promotes the development of novel drugs with a potential affinity for therapeutic targets, has turned out to be an essential tool for rapid and cost-efficient lead discovery. SBDD is a growing, iterative, and powerful approach that involves the structural evaluation of targets in the drug discovery process. It has the ability to reduce the time and cost of developing new drug lead molecules with potential therapeutic effects [31]. A virtual ligands database screening, followed by molecular docking, molecular dynamics (MD), and three-dimensional quantitative structure-activity relationship (3D-QSAR) analysis, as previously reported [2,3,5,32], provides an excellent workflow ( Figure 1) to assess an efficient SBDD. Here, we analyzed 27 molecules selected from the results of structure-based virtual screening following this workflow and characterizing the most important features of an ideal ligand/inhibitor for the intracellular RTK domain of human PDGFRα.

Virtual Screening of the Small Compounds Database
The structure-based virtual screening of the databases (as descripted in the Methods section) against the human intracellular PDGFRα (hiPDGFRα) provided an early version of the best compounds list (Tables S1-S5). This was evaluated by taking the best ligands of each of the 5 drug-like chemical libraries, resulting in the 27 ligands reported in Table  1. Among these libraries, the Drugs-lib had the most involvement, showing ligands with higher affinities, although the natural products selected showed Kd values in the sub-micromolar range, indicating medium-high affinities. Here, we analyzed 27 molecules selected from the results of structure-based virtual screening following this workflow and characterizing the most important features of an ideal ligand/inhibitor for the intracellular RTK domain of human PDGFRα.

Virtual Screening of the Small Compounds Database
The structure-based virtual screening of the databases (as descripted in the Methods section) against the human intracellular PDGFRα (hiPDGFRα) provided an early version of the best compounds list (Tables S1-S5). This was evaluated by taking the best ligands of each of the 5 drug-like chemical libraries, resulting in the 27 ligands reported in Table 1. Among these libraries, the Drugs-lib had the most involvement, showing ligands with higher affinities, although the natural products selected showed K d values in the submicromolar range, indicating medium-high affinities.

Molecular Docking
The molecular docking analysis, which was used to check the virtual screening results, provided affinity values reported in Table 1 and binding geometries of the ligand/hiPDGFRα complexes. These values are in the nanomolar/sub-micromolar range and confirm the ranking obtained by the structure-based virtual screening ( Figure S1).
To validate the docking process, we compared the Imatinib/hiPDGFRα predicted complex with the available crystal structure, obtaining a ligand-ligand superimpose root mean square deviation (RMSD) value of 0.198 Å. This value is significantly lower than 1 Å, indicating a high reliability of the method.
All complexes involved the hiPDGFRα active site, although the docking box was set for the entire protein, indicating no other binding site. The four best ligands interacted with hiPDGFRα through the formation of five hydrogen bonds with Glu644, Thr674, Cys677, and Asp836, showing very similar geometries. As reported in Figure 2A, Bafetinib also interacted with the protein, forming six hydrophobic interactions with Val607, Glu644, Ile647, Thr674, Leu825, and Phe837, and one π-cation interaction with Lys627. This π-cation interaction was also formed in the Flumatinib/hiPDGFRα complex, which additionally showed nine hydrophobic interactions (with Leu599, Val607, Lys627, Glu644, Val658, Thr674, Leu825, Asp836, and Phe837) and a halogen bond between Glu644 and the trifluoromethyl group of the ligand ( Figure 2B). This π-cation interaction, the halogen bond, and the distances and number of hydrophobic interactions between the ligand and the hiPDGFRα seem to be crucial for binding affinity. Although the Imatinib/hiPDGFRα crystal structure shows seven hydrogen bonds (with Glu644, Met648, Thr674, Cys677, His815, Val816), a salt bridge (Asp836), two π-cation interactions (with His816 and Tyr676), and seven hydrophobic interactions (with Leu599, Lys627, Glu644, Val658, Ile672, Thr674, and Asp836), the predicted Imatinib/hiPDGFRα complex, after the molecular dynamics analysis, shows only five hydrogen bonds previously described and nine hydrophobic interactions (with Leu599, Val607, Lys627, Ile672, Thr674, Tyr676, Leu825, Asp836, and Phe837). This could be explained by the difference in temperature between the two types of tests (293 K for the crystallographic analysis and 300 K for the molecular dynamics).
Regarding the binding between natural ligands such as curcumin and quercetin and hiPDGFRα, this involves only a minor number of hydrogen bonds and hydrophobic interactions, even if their strength and number establish a medium-high affinity (Kd,pred = 0.2 µM and Kd,pred = 0.7 µM, respectively). This still determines an excellent competition against ATP that has an affinity for the wildtype enzyme in the order of sub-millimolar (Km = 179.6 µM) [12], underscoring the potency of their enzymatic inhibition capabilities. Although the Imatinib/hiPDGFRα crystal structure shows seven hydrogen bonds (with Glu644, Met648, Thr674, Cys677, His815, Val816), a salt bridge (Asp836), two π-cation interactions (with His816 and Tyr676), and seven hydrophobic interactions (with Leu599, Lys627, Glu644, Val658, Ile672, Thr674, and Asp836), the predicted Imatinib/hiPDGFRα complex, after the molecular dynamics analysis, shows only five hydrogen bonds previously described and nine hydrophobic interactions (with Leu599, Val607, Lys627, Ile672, Thr674, Tyr676, Leu825, Asp836, and Phe837). This could be explained by the difference in temperature between the two types of tests (293 K for the crystallographic analysis and 300 K for the molecular dynamics).
Regarding the binding between natural ligands such as curcumin and quercetin and hiPDGFRα, this involves only a minor number of hydrogen bonds and hydrophobic interactions, even if their strength and number establish a medium-high affinity (K d,pred = 0.2 µM and K d,pred = 0.7 µM, respectively). This still determines an excellent competition against ATP that has an affinity for the wildtype enzyme in the order of sub-millimolar (K m = 179.6 µM) [12], underscoring the potency of their enzymatic inhibition capabilities.

Molecular Dynamics
The stability of the best 4 predicted complexes, obtained through molecular docking, was analyzed using molecular dynamics simulation for 50 ns. The first four nanoseconds of these simulations are represented in Figure 3 since no further RMSD variations were established. The calculated RMSD of each ligand and complex trajectory with respect to the backbone indicates that the complexes achieve sufficient stability (RMSD < 0.17 nm, 1.7 Å) within less than 4 ns, except for Flumatinib ( Figure 3G), which shows RMSD values of about 0.22 nm.

Molecular Dynamics
The stability of the best 4 predicted complexes, obtained through molecular docking, was analyzed using molecular dynamics simulation for 50 ns. The first four nanoseconds of these simulations are represented in Figure 3 since no further RMSD variations were established. The calculated RMSD of each ligand and complex trajectory with respect to the backbone indicates that the complexes achieve sufficient stability (RMSD < 0.17 nm, 1.7 Å) within less than 4 ns, except for Flumatinib ( Figure 3G), which shows RMSD values of about 0.22 nm. The RMSD trend of Bafetinib ( Figure 3A) and Flumatinib ( Figure 3G) with respect to the other ligands seems to require a few more picoseconds and nanoseconds, respectively, to reach final stability. This behavior could be explained by the establishment of the πcation interaction.
The average short-range Coulombic and Lennard-Jones interaction energies calculated during the molecular dynamics analysis are reported in Table 2. As noted in the molecular docking section, the electrostatic contributions between Bafetinib and Flumatinib are essentially equivalent, while there is a significant difference in terms of intermolecular pair Lennard-Jones potential, indicating a more efficient binding of Bafetinib to the active site.  The RMSD trend of Bafetinib ( Figure 3A) and Flumatinib ( Figure 3G) with respect to the other ligands seems to require a few more picoseconds and nanoseconds, respectively, to reach final stability. This behavior could be explained by the establishment of the π-cation interaction.
The average short-range Coulombic and Lennard-Jones interaction energies calculated during the molecular dynamics analysis are reported in Table 2. As noted in the molecular docking section, the electrostatic contributions between Bafetinib and Flumatinib are essentially equivalent, while there is a significant difference in terms of intermolecular pair Lennard-Jones potential, indicating a more efficient binding of Bafetinib to the active site.

Three-Dimensional Quantitative Structure-Activity Relationship (3D-QSAR) and ADMET Analyses
The three-dimensional QSAR analysis produced comparative molecular fields analysis (CoMFA) models considering the CoMFA potential or field, the best of which are summarized in Table 3 based on optimal principal components (OPC). Among these, MW, nHA, nHD, LogP, TPSA, rotatable bonds (RBs), molar refractivity (MR), length, and max length were taken into account and modulated to achieve the best model. Some of them are reported in Tables 4 and 5. Abbreviations: molecular weight (MW), logarithm of the n-octanol/water distribution coefficient (LogP), logarithm of the n-octanol/water distribution coefficients at pH = 7.4 (LogD), number of hydrogen bond acceptors (nHA), number of hydrogen bond donors (nHD), and topological polar surface area (TPSA). Figure 4 shows the relationships between experimental and calculated activities by CoMFA models. The fitting statistical results, also reported in Table 3, indicate that the steric model (STE) is the more consistent and robust one, even if the electrostatic one (ELE) shows a relatively good squared correlation coefficient r 2 and internal predictive coefficient q 2 . This would confirm the result of interaction energies obtained by Gromacs, which describes the best complex in terms of affinity as the one with the higher Lennard-Jones energy contribution, the forces involved in the steric effect. ADMET predictions, in particular, some properties and relative descriptors of the 27 selected molecules, are reported in Table 4. A major part of these compounds satisfies the conditions regarding the acceptance of the first two indexes (Lipinski and Pfizer) related to absorption, permeability, and toxicity. Only Ditercalinium does not satisfy the conditions of these indexes since it presents the heaviest and longest molecule and the highest MR and LogP, which are out of range. Several compounds could cause skin sensitization and respiratory toxicity-although no compound would cause acute toxicity during oral administration.  ADMET predictions, in particular, some properties and relative descriptors of the 27 selected molecules, are reported in Tables 4 and 5. A major part of these compounds satisfies the conditions regarding the acceptance of the first two indexes (Lipinski and Pfizer) related to absorption, permeability, and toxicity. Only Ditercalinium does not satisfy the conditions of these indexes since it presents the heaviest and longest molecule and the highest MR and LogP, which are out of range. Several compounds could cause skin sensitization and respiratory toxicity-although no compound would cause acute toxicity during oral administration.  Abbreviations: drug-likeness based on the concept of desirability (QED), Natural Product-likeness score (NPlikeness), human colon adenocarcinoma cell lines permeability (Caco-2), Madin−Darby Canine Kidney cells permeability (MDCK), skin sensitization (Skin Sen), eye irritation/corrosion (EI/EC) potential, respiratory toxicity (Respiratory), and acute toxicity during oral administration (LD50 oral).

Discussion
The PDGFRα is a critical factor that plays an essential role in regulating cell proliferation, survival, and chemotaxis [24][25][26][27]. Several studies have identified PDGFR as a potential target to treat pathologies, such as cancer and progressive fibrotic diseases [4][5][6][7][8][14][15][16][17], by searching for synthetic or natural compounds acting as ligands/inhibitors. Currently, however, no one has considered involving such large compound libraries for this purpose, such as those implemented in the MTiOpenScreen web server.
The details of binding interactions between human intracellular PDGFRα and the compounds selected after the structure-based virtual screening were revealed by molecular docking and the subsequent MD simulations. Molecular docking analysis, in addition to complexes' affinities, was able to provide important critical issues regarding the geometries. The predicted affinities for the 27 selected complexes vary from the nanomolar order (1.01 nM for the Bafetinib/hiPDGFRα complex) to the micromolar one, a typical range of values for a good drug inhibition and reliable QSAR analysis. Some of the most stable complexes show a π-cation possible interaction with Lys627, halogen bonds with Glu644, and numerous hydrophobic interactions.
MD analysis reveals that the complexes achieve their stability in a very short time and remain stable for all the MD simulations, demonstrating no substantial unfolding of the structures. The type of interaction energies calculated during this analysis shows that the steric contribution (the Lennard-Jones energy) is predominant in the most stable complexes. This behavior is confirmed by CoMFA models obtained through the three-dimensional QSAR. The steric model field (or at least a combination of the two models) produces a very consistent activity predictor model-suggesting the descriptors included in this field and the structural information obtained by molecular docking analysis-as a basis for designing new effective drugs. These predicted models and structural information confirm the molecule ranking reported in Table 1 and identify Bafetinib as the best candidate for a hiPDGFR ligand/inhibitor.
Although this analysis was carried out with bioinformatics tools and needs experimental evidence and confirmation, we can assert that it has the potential to create a platform for developing new drugs or repositioning already-known drugs. Though the predicted model is an excellent starting point for developing an ideal hiPDGFRα ligand/inhibitor, all 27 selected compounds have the potential to effectively bind and inhibit hiPDGFRα activity based on predicted geometries and affinities. This is partially supported by the fact that compounds such as Imatinib, Nilotinib, Crenolanib, and Nintedanib, which are included in this set, have already been experimentally shown to have activity against hiPDGFR [1,4,8,9,11,15]. Moreover, Bafetinib is already known to inhibit the Bcr/Abl fusion protein tyrosine kinase and the Src-family member Lyn tyrosine kinase for treating Bcr-Abl+ leukemias, including chronic myelogenous leukemia (CML) and Philadelphia+ acute lymphoblastic leukemia [33,34]. Further studies are necessary to determine whether these compounds should proceed to the pre-clinical and clinical phases of evaluation.
As reported in a previous study on Nintedanib [35], it will be crucial to carefully adjust the dosage of these new drugs in a real-life setting to ensure their optimal efficacy and safety for patients. This requires close monitoring of potential side effects and close collaboration between healthcare providers and patients to make necessary adjustments based on individual patient factors, such as age, weight, and comorbidities.

Virtual Screening of the Small Compounds Database
The structure-based virtual screening of compounds, such as drugs, food constituents, and natural products, against the human intracellular PDGFRα was performed using the MTiOpenScreen web server [36,37] that is based on Autodock Vina-based docking genetic algorithm and analyses 5 compounds libraries: iPPI-lib, diverse-lib, Drugs-lib, Food-lib, and NP-lib. These drug-like chemical libraries contain modulators of protein-protein interactions, diverse chemical compounds, purchasable approved drugs, food constituents, and natural products, respectively, currently for a total of about 7200 molecules. After the preparation of the receptor (pdbID: 6jol) [12] removing Imatinib and water molecules, the PDB file was uploaded on the server setting the grid center coordinates (64.868, 43.016, 0.324) and the size of the search space (30 × 28 × 36 Å), keeping all the other parameters as default.

Molecular Docking
Twenty-seven molecules were mainly selected on the basis of the best binding energies obtained from the results of the MTiOpenScreen web server. To check these results, we carried out a further docking analysis using the SwissDock web server, as previously reported [38], based on the EADock DSS algorithm [39,40], and setting the docking type as accurate without a definition of the region of interest to include the whole protein in the search space. The final geometry of the complexes obtained was rendered by PyMol software (The PyMOL Molecular Graphics System, Version 2.5.2, Schrödinger, LLC., Cambridge, MA, USA), and the 3D representation of the interaction was obtained by the Protein-Ligand Interaction Profiler (PLIP) web server [41].

Molecular Dynamics
The molecular dynamics simulation of the best 4 complexes obtained (including the Imatinib/receptor complex as reference) was performed using Gromacs software version 2022.4 [42], the CHARMM36 force field updated July 2022 [43], the SPC216 water model, and the CGenFF server for generating ligand topologies [43]. Each dynamics simulation started with solvation, neutralization adding three chlorine ions, and minimization of the system using the steepest descend algorithm for 50,000 steps. Next, the NVT and NPT equilibration stage was performed in a 100 ps run, stabilizing the temperature at 300 K (using V-rescale Berendsen thermostat) and the pressure at 1 bar (using C-rescale pressure coupling), respectively. Finally, MD simulation was run for 50 ns, and the rms tools of Gromacs were used to determine the root mean square deviation (RMSD) of each ligand trajectory in relation to the protein backbone and, also, to calculate the protein backbone RMSD relative to the energy minimized conformation. The Coulombic and Lennard-Jones interaction energies were calculated using the capability of Gromacs to decompose the short-range non-bonded energies.

Three-Dimensional Quantitative Structure-Activity Relationship (3D-QSAR) and ADMET Analyses
To simplify the design of selective human intracellular PDGFRα inhibitors, the threedimensional quantitative structure-activity relationship (3D-QSAR) analysis was carried out using the comparative molecular fields analysis (CoMFA) [44] available on the web portal www.3d-qsar.com [45]. To achieve this analysis, we first performed an alignment of molecules' conformations and a molecular interaction fields (MIFs) calculation of the 27 selected molecules. We set all parameters as default. The two CoMFA potentials obtained, called steric (STE) and electrostatic (ELE), were calculated by means of the Lennard-Jones and Coulomb law definitions, respectively.
The ADMET (absorption, distribution, metabolism, excretion, and toxicity) analysis was carried out using the screening section of ADMETlab 2.0 web server [46] and uploading the SMILES notation of the 27 selected molecules. In addition to the most common molecular descriptors, the results provide four acceptance indexes called Lipinski, Pfizer, GSK, and Golden Triangle. The Lipinski rule is related to absorption and permeability, while the Pfizer rule is related to toxicity. Compounds satisfying the GSK and Golden Triangle rule may have a more favorable ADMET profile.