The 3D structure of FGF2 was selected in this study as a drug target and retrieved from the Protein Data Bank (PDB) by considering a specific four-digit code of 6L4O. In the given crystal structure, FGF2 is complexed with API5, both of which are highly expressed in different cancers including glioblastoma (Bong et al., 2020). The structure was determined by X-ray diffraction with a resolution value of 2.60 Å. Immediately, the structure was processed in UCSF Chimera v1.15 to remove the AP15 co-crystalized ligand (Kaliappan and Bombay, 2018). The apo-FGF2 was energy minimized via two algorithms: steepest descent and conjugate gradient algorithms. Each algorithm was run for a maximum of 1,500 steps, while the step size was set to 0.02 Å. During the process, the missing hydrogen atoms were added while appropriate charges were assigned. After minimization, the structure was saved into a .pdb format to make it ready to be used in virtual screening processes.

In this study, different libraries of Asinex databases were used against the 3D structure of FGF2. These libraries are the BioDesign library, elite library, synergy library, and gold and platinum library [https://www.asinex.com/screening-libraries-(all-libraries)]. The compounds in the mentioned libraries can be easily utilized in a rapid hit to lead optimization and have a high degree of drug-likeness. Collectively, these libraries contain ∼ 575302 compounds that can be easily accessed for experimental testing. The drug library was imported to PyRx 0.8 (Dallakyan and Olson, 2015), energy minimized via the MM2 force field, and subsequently converted into .pdbqt format (Halgren, 1996).

The virtual screening process is a highly efficient technique to filter hits that bind best to a receptor biomolecule. The structure-based virtual screening process docked the library compounds one by one to the receptor molecule with specified coordinates and ranked the best binders on the lowest binding energy (Lionta et al., 2014; Shaker et al., 2021). The lowest binding energy compounds achieved stable conformation and are, thus, used in further investigations. Herein, the coordinates used were of the following residues: Asn169, Arg223, Arg262, Thr263, Lys267, Lys271, and Lys277 (Bong et al., 2020). The number of iterations generated for each of the library compound was set to 100, and only the best eight stable iterations were shortlisted. The grid box set around the aforementioned active site residues were in 12.49 Å on the X-axis, −6.804 Å on the Y-axis, and −47.32 Å on the Z-axis. The grid box sizes along XYZ planes were 26.50 Å, 25.0 Å, and 28.75 Å, respectively. For validation of the docking protocol, the co-crystalized ligand AP15 was redocked with FGF2 at the same site, as reported in the crystal structure. The exact binding conformation achieved as produced in the crystal structure denotes accuracy of the virtual screening process. The top three docked complexes from this process were used further for additional investigations of binding conformational analysis and interaction analysis by UCSF Chimera 1.15 (Kaliappan and Bombay, 2018), Discovery Studio Visualizer v2021 (Biovia, 2017), and PDBsum generate tool (Laskowski, 2001).

The selected three complexes of compounds that showed a best binding affinity for FGF2 were opted further for dynamics studies. This was accomplished using the AMBER20 simulation package (Case et al., 2020). Initial processing of the complexes was performed using the Antechamber program (Wang et al., 2001). Both the FGF2 receptor and compounds were subjected to parameter generation using the AMBER tleap interface. The force field applied for the FGF2 receptor was FF14Sb, while for the compounds, the force field used was GAFF (Maier et al., 2015; Sprenger et al., 2015). Each complex was solvated into a TIP3P water model, followed by neutralization and parameter correction. Steric clashes of complexes were discarded by subjecting the complexes to the energy minimization process for a total of 3,000 steps (the first 1,500 steps were performed via the steepest descent algorithm, and the last 1,500 steps were accomplished through the conjugate gradient algorithm). Heating was then performed to heat the systems up to 310 K. Systems equilibration was accomplished by running it for 100 ps in the NVT ensemble. The simulation time length was set to 200 ns, where the SHAKE algorithm was used to restrain hydrogen bonds (Kräutler et al., 2001). Langevin dynamics were employed to obtain a stable temperature during a simulation production run (Izaguirre et al., 2001). The simulation trajectories were investigated for several structural parameters based on alpha-carbon atoms of the systems and performed using the CPPTRAJ module of AMBER (Roe and Cheatham, 2013). Visualization of different simulation frames was performed via visual molecular dynamic (VMD) v 1.93 software (Humphrey et al., 1996; Falsafi-Zadeh et al., 2012). Plots to see the structural stability of complexes were produced by XMGRACE v 5.1 (Tuccinardi, 2021).

Next, the system intermolecular binding free energies were calculated by the MMPBSA.py module (Miller et al., 2012). Two methods were used for this, MM-PBSA and MM-GBSA, to cross-validate the docking affinity of compounds for FGF2 (Tuccinardi, 2021; Wang et al., 2019). The MM-PB/GBSA estimates the energy difference between the docked receptor–ligand complex and the apo protein and ligand alone (Hou et al., 2011). For the calculation, 1,000 frames were selected for simulation trajectories at a regular time interval of 0.2 ns. The MM-PB/GBSA free energy estimation was performed via equation I. Furthermore, the residue-based binding free energy was estimated to highlight favorable contributing residues in anchoring and stabilizing docked compounds. Δ G b i n d i n g = E b i n d i n g + E e l + E v d w + G p o l + G n p − T S .

The radial distribution function (RDF) analysis was conducted for strong intermolecular interactions between FGF2 and compounds (Donohue, 1954). This analysis was vital to shed light on how key chemical interactions are vital for keeping the compound binding conformation stable with the receptor. The RDF demonstrates the density of interatomic radii interactions versus time. A high and stable interaction density plays a crucial role in long-term inhibition of the receptor and its binding with interacting partners. The RDF analysis was accomplished through the simulation length and performed via an AMBER CPPTRAJ module (Roe and Cheatham, 2013). Plotting of the interaction data was achieved through XMGRACE v 5.1 (Turner, 2005).

The entropy energy contribution to an overall net binding energy of the systems was estimated using AMBER normal mode calculations (Genheden et al., 2012). In this calculation, different entropy energies were generated such as translational, vibrational, rotational, and net entropy contribution. For each parameter, average, standard deviation, and errors were also calculated. As the normal mode entropy energy calculation is very expensive, only a limited number of frames from the simulation trajectories were processed.

The WaterSwap method is considered more refined than the conventional MM-PB/GBSA method (Woods et al., 2011; Woods et al., 2014). This method considered the contribution of water molecules in the interaction between the ligand and FGF2 active site residues that are commonly skipped in the MM-PB/GBSA method (Bergström and Larsson, 2018). During the procedure, the solvation box buffer was set at 10.0 Å. The iteration number used was 1,000, and the maximum threads were treated as Auto. The water monitor distance used was 7.00 Å. Three types of energies were calculated, namely, thermodynamic integration (TI), free energy perturbation (FEP), and Bennetts. The lowest energy implies a strong intermolecular docked stability of the complexes. The energy value difference among the algorithm must be less than 1 kcal/mol to get a well-converged system (Raza et al., 2019).

The drug-like properties and pharmacokinetics of the filtered predictions were investigated through SwissADME (Daina et al., 2017) and pkCSM (Pires et al., 2015).

AP15 and FGF2 are highly expressed in several cancer types and associated with poor prognosis. The binding conformation of these proteins with each other is given in Figure 2, while the chemical interactions they produced are provided in Figure 3. The intermolecular interaction produces seven salt bridges, five hydrogen bonds, and 48 non-bounded contacts. AP15 interacts with FGF2 through 10 interface residues, while nine residues from FGF2 contribute their role in binding with AP15.

The virtual screening process identified 81 compounds that docked well with the FGF2 receptor. The top stable binding conformation of these compounds with FGF2 can be categorized according to their binding energy value, as given in Figure 4. The top three best binding compounds were chosen: Asinex 51412372 ((E)-2-(3′-benzyl-2′, 5′-dioxo-2.2′, 3.3′, 4′, 5.5′, 6, 10a′, 11′, 14′, 14a′-dodecahydro-1′H-spiro[pyran-4,6′-pyrido[4,3-e][1,4]diazacyclododecin]-12′ (7′H,10′H, 13′H)-yl)-N-(4-phenoxyphenyl)acetamide), Asinex 51217461 (2-(2-(benzylsulfonyl)-3,3-dimethyl-2,8-diazaspiro [4.5]decan-8-yl)-N-(4-(piperidin-1-ylsulfonyl)phenyl)acetamide), and Asinex 51216586 (2-(3-benzyl-4-tosyl-1-oxa-4,8-diazaspiro [5.5]undecan-8-yl)-N-(5.6, 7, 8-tetrahydronaphthalen-2-yl)acetamide) with binding energy scores of –8.3 kcal/mol, −8.2 kcal/mol, and −7.8 kcal/mol, respectively, for further analysis. These compounds were prioritized as the best docked inhibitor upon a repeated docking run and compared to control. The control AP15 was found to have a binding energy of −12.31 kcal/mol.

Details about compound binding and interactions with the FGF2 are given in Figure 5. The compounds interact with the same interface reported in the crystal structure. The Asinex 51412372 N-(4-phenoxyphenyl)acetamide chemical structure formed a strong hydrogen bond with Asp183, while the rest of the compound structure ((E)-3′-benzyl-2, 3,3′, 4′,5.6, 10a′, 11′, 12′, 13′, 14′, 14a′-dodecahydro-1′H-spiro[pyran-4, 6′-pyrido[4, 3-e][1, 4]diazacyclododecine]-2′, 5′ (7′H, 10′H)-dione) produced several van der Waals bonds. Asinex 51217461 terminal structures [(hydrosulfonylmethyl) benzene and 1-hydrosulfonylpiperidine] generated hydrogen bonds with Arg223 and Arg175, while the central 1-hydrosulfonylpiperidine made several van der Waals interactions. Asinex 51216586 formed a hydrogen bond with Gln265, and its oxygen atom of 5.6, 7, 8-tetrahydronaphthalen-2-amine, while 1-hydrosulfonyl-4-methylbenzene is involved in the hydrogen bond with Arg175. Previously, the pan-fibroblast growth factor receptor (FGFR) inhibitor termed as LY2874455 was identified and is currently under phase I of a clinical trial investigation. This inhibitor is known to inhibit the same active pocket and interact with the same set of active residues, as reported here by the filtered molecules (Michael et al., 2017). In another study, Mahfuz et al. (2022) reported a novel compound (PubChem 137300327) to show stable interactions with FGF2. PubChem 137300327 was docked with several PDBS of FGF2 and reported a score of –10.3 kcal/mol with 561M, −11.3 kcal/mol with 564F, −11.2 kcal/mol with 564I, −11.0 kcal/mol with 564 M and −10.5 kcal/mol with 555 L, and −9.6 kcal/mol with 555 M.

The understanding of the molecular structure and dynamics is critical toward underpinning molecular function and biology. Molecular dynamic simulation is considered an integral part of modern computational drug discovery as it aids in the study of atomic-level dynamics and conformational flexibility of protein–drug complexes (Karplus, 2002). This analysis is important to decipher biological events of receptor–drug binding and interactions in theoretical simulations. It is now in regular use in CADD and provides detail and accurate estimation of receptor–drug thermodynamics and kinetics (Moradi et al., 2021). The simulation indicated all the docked complexes including the control relatively stable from the perspective of the structure. The FGF2_51412372 and FGF2_51217461 systems were reported to be more dynamically stable than FGF2_51216586 and control. The first 50 ns RMSD of all systems depicted constant dynamics, followed by divergence. The AP15 control and FGF2_51216586 afterward to around 150 ns demonstrated continuous small structural variations. These variations when dissected revealed regular conformational jumps by the ligand at the FGF2 docked site to attain more stable conformation, which can be witnessed in the last phase of simulation. In the latter part of simulation, conformational equilibrium can be seen for all systems. The maximum RMSD touched by the control and FGF2_51216586 is around 3 Å. The major RMSD jump for FGF2_51412372 noticed is 2 Å, while that for FGF2_51217461 is 1.5 Å. The structural deviations in the systems are due to loop dynamically flexible regions, which upon ligand pressure, behave more unstable. The RMSD plot of each system can be seen in Figure 6A. The study by Mahfuz et al. also revealed very stable dynamics of PubChem 137300327 with GFR2 in 100 ns of simulation time. The RMSD of the complex touches the maximum of 2 Å despite some random conformation variations between 20 and 60 ns (Mahfuz et al., 2022). Further observation on the validation of the system stability was gained through RMSD analysis to get an understanding whether the FGF2 residues are stable or unstable in the presence of compounds/control. A quite similar trend to that of RMSD was found in this study (Figure 6B). The maximum RMSF deviation for the systems was seen around residues 30–40 and 120–130. The maximum RMSF value was found for control and FGF2_51217461 that is 5 Å. Generally, all the systems were discovered to have low a RMSF value for most of the receptor residues. The hydrogen bonds involving residues of the FGF2 RMSF value are as follows: Arg15 (average 1.44 Å), Asp23 (average 1.68 Å), Arg63 (average 1.02 Å), and Gln105 (average 1.82 Å). Next, RoG analysis was performed to understand the compactness and the relax nature of the FGF2 receptor during simulation time (Figure 6C). A higher RoG reported a highly relaxed nature, and thus, the detachment of ligand is easy. The control system is the most unstable complex compared to other studied systems with an average RoG value of 69.52 Å. The RoG of PubChem 137300327 with GFR2 showed the same compact structure, as reported in this study.

Furthermore, RDF analysis of the key interactions between the compounds and FGF2 residues was performed to determine how the key intermolecular interaction density play a role in compound tight stability along the simulation length. It can be seen that Arg15 and Asp23 provided considerable contributions in holding the compounds at the FGF2 docked site. The maximum RDF value noticed for the Asinex 51216586–Arg15 interaction is at 1.8 Å with an RDF value of 0.14. For Asinex 51217461–Arg15, the maximum RDF value is noticed at 0.10 at a distance of 1.96 Å. Last, the maximum RDF observed for Asinex 51412372–Asp23 is 0.3 at a distance of 2 Å. The RDF plots are given in Figure 7.

The docking study results are often false positive, and error chances are high. It is usually preferred that the docking calculations must be revalidated by more sophisticated methods such as MM-GBSA and MM-PBSA (Genheden and Ryde, 2015). These methods are considered to be more accurate than scoring functions and are modest in the use of computational resources. They are a more preferred choice to be carried out than the docking algorithms and, thus, are applied in the current study to verify the binding affinity of the shortlisted compounds for FGF2. The main difference between MM-GBSA and MM-PBSA is the way ΔG _PB/GB term is estimated. In the former, this calculation is faster and provides more approximation of the GB model, while in the latter, this term calculation is time-consuming. By using these methods, it was concluded that systems are energy-wise very much stable. The control in each case was reported to produce robust binding free energies. The net energy of control was −86.72 kcal/mol in MM-GBSA and −85.06 kcal/mol in MM-PBSA. For the control system, both van der Waals and electrostatic energy terms dominated the overall system energy and provided a favorable contribution in compounds binding to FGF2. In MM-GBSA, the ranking of compound systems in terms of energy stability was in the following order: FGF2_51217461 (−66.42 kcal/mol) > FGF2_51412372 (−59.24 kcal/mol) > FGF2_51216586 (−42.74 kcal/mol). In MM-PBSA, FGF2_51217461 achieved the most optimal energy state with a net energy of −64.17 kcal/mol, followed by FGF2_51412372 (−56.79 kcal/mol) and FGF2_51216586 (−42.28 kcal/mol). Similarly, like control, for compounds, the van der Waals energy seems the most favorable in terms of intermolecular binding, followed by electrostatic energy. The polar energy in all systems provided unfavorable contributions. Each energy term calculated for the systems can be found in Table 1.

The entropy contribution to the net binding free energy of the systems is given in Table 2. The values demonstrated the presence of entropy energy in the systems due to the ligand pressure, as noticed in the simulation plots (Genheden et al., 2012). However, the entropy contribution is much less than the net binding energy of the systems and, thus, points to a stable docked nature of the complexes, hence the long-term inhibition of FGF2. This analysis again demonstrated the stable nature of complexes in the presence of compounds. All the simulation-based analysis proved that docked compounds are energy-wise very much stable with the receptor and enjoy the intermolecular interactions.

Although the MM-GBSA and MM-PBSA binding free energy calculations are reliable in predicting the binding affinity of compounds to FGF2, still they suffer from several limitations. In particular, the role of water molecules connects the ligand with FGF2 active site residues (Woods et al., 2014). The WaterSwap is a more reliable method and actually swaps the ligand and the surrounding water molecules. Three algorithms were employed to revalidate the net binding free energy of the systems. Herein, the control system was unveiled again to the most stable system with net energies of −30.25 kcal/mol (Bennetts), −31.43 kcal/mol (FEP), and −30.25 kcal/mol (TI). Among compound systems, FGF2_51216586 was found to be the most favorable system with a net energy of −27.85 kcal/mol, −27.66 kcal/mol, and −27.85 kcal/mol in Bennetts, FEP, and TI, respectively. All the systems were seen to be well converged as the energy difference among the algorithms is less than 1 kcal/mol (Ahmad et al., 2019). The WaterSwap energies of the control and compound systems are given in Figure 8.

Prediction of drug-like properties of the compounds and pharmacokinetics holds significant importance in drug development as they contribute majorly to reduce drug failure during trial testing (Daina et al., 2017). The oral bioavailability radar of the studied compounds is given in Figure 9. As the compounds fulfill most of the radar parameters, they are likely to be good drug molecules for additional structure optimization. As can be seen in the figure, most of the parameters of the radar are within the range of the compounds and are, thus, favorable candidates from the drug-likeness point of view. This was also demonstrated by Egan and Veber-like rules that indicated the compounds fulfilled all the parameters to be considered drug-like compounds and, hence, have higher chances to be successful in the experimental studies.

The compounds have a lower topological polar surface area (TPSA). For example, FGF2_51412372, FGF2_51217461, and FGF2_51216586 TPSA values are 109.00 Ų, 123.86 Ų, and 87.33 Ų, respectively. These values demonstrated the good ability of the compounds to cross the cell membrane and reached the target site, thus increasing the compound bioavailability for rapid therapeutic action (Veber et al., 2002). Similarly, the compounds are moderately water soluble, making them good candidates for oral formulation (Alex, 2003). The compounds have high gastrointestinal absorption and do not cross the blood–brain barrier (BBB). From a drug-likeness point of view, the compounds are considered drug-like by Veber and Egan’s drug rules (Veber et al., 2002; Ahmad et al., 2018). The compounds have a good synthetic accessibility score, which means that they can be easily synthesized in a laboratory to be used in in vitro and in vivo experiments. Additionally, the compounds are predicted to show no toxicity and can be cleared easily from the body. The different drug-likeness and pharmacokinetic properties of the compounds are tabulated in Table 3.