All the computational studies were carried out on a PC running Windows 7 Ultimate with an Intel Core i3 microprocessor, 4.00 GB of RAM and a 64-bit operating system. AutoDock 4.2.1 (ADock) (https://autodock.scripps.edu/download-autodock4/) and AutoDock Vina 1.1.2 (Avina) (https://vina.scripps.edu/downloads/) were used to benchmark molecular docking. Discovery Studio Visualizer 3.5 (https://discover.3ds.com/discovery-studio-visualizer-download) and PyMol 2.5 (https://pymol.org/2/) were employed to analyse the binding conformations and inter-molecular interactions. OpenBabel 3.1.1 (https://sourceforge.net/projects/openbabel/) was used to create PDBQT files for all test compounds and decoys. To validate and compare the docking results, an online tool called Screening Explorer (http://stats.drugdesign.fr/) was used.

The 3D crystal structure of AR (PDB accession number ID: 2INE with the best resolution 1.90 Å) was downloaded from the RCSB PDB (http://www.rcsb.org) in PDB format. The native ligand, cofactors, metal ions and water were removed from the protein structure. The residues that constitute the binding pocket of AR were analysed using the Computed Atlas for Surface Topography of Proteins (CASTp), a web-based tool for binding site prediction (http://cast.engr.uic.edu). The 3D structures of 58 TFG compounds and Epalrestat were obtained from the PubChem database (http://pubchem.ncbi.nlm.nih.gov/pccompound). The SDF file of 3D structures of decoys for the AR enzyme was downloaded from the DEKOIS 2.0 database (http://www.dekois.com).

The target protein in PDB format obtained from the above step was further processed by adding polar hydrogens and Gasteiger charges in Avina. This energy-minimised AR enzyme structure was saved as a PDBQT file and was used for docking the ligands in ADock and Avina. All test compounds, including TFG compounds, Epalrestat and decoys in SDF format, were loaded into Open Babel for conversion to PDBQT, during which the molecules are energy minimised by adding polar hydrogen and Gasteiger charges. The binding site was defined by positioning a grid box with a centre of X, Y and Z dimensions 18.61 X -11.39 X 17.71 Å and a size of 60 X 60X 65 Å, respectively, with a spacing value of 0.375 Å. Both ADock and Avina docking were performed using the exact dimensions mentioned above. The genetic algorithm parameter was set at 10 in ADock, while all other parameters were fixed at default values. In Avina docking, ten conformers were generated at a 3 kcal/mol energy gradient, while the rest of the parameters chosen were at their default numbers. Both software allowed the ligand flexibility; their rotatable bonds could move freely during docking. The energy calculations of protein-ligand complexes were done through the Lamarckian genetic algorithm in both software. The protein-ligand complexes with the least negative binding energy (∆G) indicate robust binding and favourable conformations were selected for analysing the interaction between the protein-ligand complexes (Alshahrani et al. 2021). The 2D and 3D interactions of docked complexes were studied in Discovery Studio.

Individual and comparative scoring methods were used to compare the results of ADock and Avina docking software (Thangavel and Albratty 2022). For this purpose, Screening Explorer was employed. The binding energy results obtained from ADock and Avina, along with the ligand ID numbers for TFG compounds, standard drug and decoys, were loaded into Screening Explorer. The description of active compounds was 1, while decoys were 0, which was necessary to compare docking results. These were fed into Screening Explorer in CSV format. Screening Explorer then used the binding energy scores to rank the ligands. The ligand rankings were used to validate the scoring methods used by the docking software. The ligand’s binding energy and corresponding rankings were used to construct predictiveness, ROC and enrichment graphs in Screening Explorer. To assess the performance of ADock and Avina, widely-accepted metrics, such as auROC (Area Under Receiver Operating Characteristic Curve), BEDROC (Boltzmann-enhanced discrimination of ROC), EF (Enrichment factor), TG (Total gain) and RIE (Robust initial enhancement), were used. The results of benchmarking docking from individual, comparative scoring and binding interaction analysis were considered to choose the best TFG compound. This compound had the maximum potential to inhibit the AR enzyme and molecular dynamics further assessed the protein-ligand complex to ensure its stability.

DFT calculations were conducted to investigate the electronic properties of the best-forming TFG compounds. The calculation used Maestro’s Jaguar module with the B3LYP functional method and 6-31G* basis set. The electronic properties, such as the Highest Occupied Molecular Orbital (HOMO) energy, Lowest Unoccupied Molecular Orbital (LUMO) energy and its corresponding Band gap energy, were computed.

The GROMACS 4.3.1 package (https://www.gromacs.org/) with GROMOS43al force field was used to perform MD to evaluate the structural stability of protein-ligand complexes. The PRODRG online server (http://prodrg1.dyndns.org/) was utilised to create the ligand topology. The systems were equilibrated using a cubic box with dimensions of 10 × 10 × 10 nm³ and solvated with the Single Point Charge (SPC) water model. Appropriate counter ions (Na⁺ Cl-) were added to neutralise the systems. The steepest descent algorithm was employed for energy minimisation for 50,000 to eliminate weak van der Waals contacts. The Parrinello-Rahman barostat and modified Berendsen thermostat were applied to maintain constant pressure and temperature at 1 bar and 300 K, respectively, using the NPT and NVT ensemble. The MD simulations were conducted for 50 ns. The g_mmpbsa tool determined the binding free energy of the AR protein-Tigogenin and AR protein-Gitogenin complexes for the last 5 ns of the MD trajectories (Omolo et al. 2018). The expression used to calculate the MM/PBSA energy is:

∆G_bind = G_docked – (G_enzyme + G_compound),

wherein, G_docked is the binding energy of the docked complexes of TFG compounds with AR, G_enzyme and G_compound are the energies of AR and TFG compounds in a solvent-filled system.

The binding pocket of AR protein was analysed using the CASTp web server and the results are depicted in Fig. 1. A single binding site was identified in the AR enzyme. CASTp results indicated that Gly18, Thr19, Trp20, Lys21, Asp43, Val47, Tyr48, Lys77, Trp79, His110, Trp111, Phe122, Ser159, Asn160, Gln183, Tyr209, Ser210, Pro211, Leu212, Ser214, Asp216, Trp219, Leu228, Ala245, Ile260, Pro261, Lys262, Ser263, Val264, Thr265, Arg268, Thr268, Glu271, Asn272, Cys298, Ala299 and Leu300 are the amino acids surrounding the binding site. The predicted binding site residues correlated with the previous reports (Rondeau et al. 1992; El-Kabbani et al. 2003). The predicted binding site residues were used to draw the grid box in Auto Grid for performing docking analysis in AutoDock (ADock) and AutoDock Vina (Avina).

The benchmarking molecular docking protocol for assessing the AR inhibitory potential of TFG compounds utilising ADock and Avina involved 1200 decoys from the DEKOIS database (see Suppl. material 1). These decoys are designed explicitly for AR. These decoys matched the test ligands in less than 500 Da of molecular weight. The test set consisted of 58 TFG compounds and one standard drug, Epalrestat. As a result, 1259 molecules were in the ligand collection overall for benchmarking docking. The distribution of binding energy scores for TFG compounds, standard drug and decoys predicted by ADock and Avina are shown in Fig. 2. The rankings from Fig. 2 were used to compare the software’s predictive power by computing metrics through ROC, predictiveness and enrichment curves. The global and partial metrics computed from these graphs are presented in Table 1. Fig. 3 shows the ROC, predictiveness and enrichment curves derived from the compound rankings.

The performance of docking software in distinguishing active and inactive compounds can be evaluated using the auROC/pauROC. A higher pauROC score at the top 5% retrieval indicates better discrimination by Avina at an early stage. The TG metric reflects the software’s ability to prioritise active compounds effectively across various ranked compounds in which Avina performed better. RIE evaluates the efficiency of docking software in terms of information gain per compound. A higher RIE value for prediction by Avina suggests a more robust ranking of compounds at the initial stage. BEDROC emphasises docking software’s early recognition ability by quantifying the enrichment of top-ranked compounds. A higher BEDROC score indicates early enrichment and Avina performed better. The TPF and P(Act) scores indicate better predictive power of Avina than ADock for activity prediction.

Table 2 presents information regarding the structures and binding energy scores of 5 TFG compounds that Avina successfully retrieved within the top 5% level. Upon further analysis of the recovered compounds using two different types of software, Tigogenin and Gitogenin were identified as the top two enriched compounds in the top three ranks. Therefore, these two compounds were considered for further molecular interaction, density functional theory (DFT) calculations and molecular dynamic studies.

Fig. 4 shows the intermolecular interactions of Tigogenin, Gitogenin and the standard drug with AR’s binding site residues. Tigogenin forms hydrogen bonds with Trp111 and Cys298, a π-sigma interaction with Tyr209, alkyl interactions with Val47, Cys298, Leu300 and Leu301 and π-alkyl interactions with Trp20, Tyr48, Trp78, His110, Trp111, Phe122 and Trp219. Gitogenin forms hydrogen bonds with Trp111 and Cys298, a π-sigma interaction with Tyr209, alkyl interactions with Val47, Cys298, Leu300 and Leu301 and π-alkyl interactions with Trp20, Tyr48, Trp78, His110, Trp111, Phe122 and Trp219. The TGF compounds also interact through the van der Walls interactions with residues such as Trp20, Ala47, Tyr48, Lys77, Trp111, Gln183, Trp219, Ala299 and Leu301. These amino acid interactions align with earlier literature reports that interactions with Trp20, Tyr48, His110 and Trp111 located in the anion binding site of the AR are critical to its inhibition (Tanawattanasuntorn et al. 2021). The molecular docking study predicted that Tigogenin and Gitogenin can potentially inhibit AR.

We investigated the electronic and energetic states of the top two TFG compounds, Tigogenin and Gitogenin. Fig. 5 is the HOMO-LUMO distribution plots of Tigogenin and Gitogenin to analyse the atomic contribution of these orbitals. For Tigogenin and Gitogenin, the HOMO energy is -0.254 eV and -0.253 eV, respectively, while the LUMO energy values are -0.077 eV and -0.075 eV. The HOMO-LUMO energy gap (HLG) is -0.197 eV and 0.195 eV for Tigogenin and Gitogenin, respectively. The compounds have smaller HLG, which indicates high chemical reactivity. A higher HLG specifies higher kinetic stability and low chemical reactivity. Since the selected ligands reveal a lower HLG, they are more chemically reactive. The energy value of HOMO is higher than the LUMO energy for the ligands, demonstrating their ability to donate electrons to the binding site of the target.

Thus, the result suggests that orbitals of such ligands are favourable for forming hydrogen bond interactions with binding site amino acid residues of AR.

The Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF) and hydrogen bond interactions of the docked complexes were analysed and presented (Fig. 6). Fig. 6a shows the RMSD plots of apoprotein AR, AR with Tigogenin and AR with Gitogenin complexes. The RMSD of the apoprotein achieved stability after 15 ns of equilibration and maintained a stable RMSD value of 0.56 nm. The RMSD plot of AR with Tigogenin and Gitogenin complexes demonstrated stability, with no significant deviations observed during the simulation. The AR with Tigogenin complex maintained an acceptable RMSD value of 0.53 nm throughout the simulation. The AR with Gitogenin equilibrated in 5 ns and had a constant RMSD value of 0.58 nm throughout the simulation.

Fig. 6b shows the RMSF plot of AR with Tigogenin and AR with Gitogenin complexes during a 50 ns simulation period. The average RMSF value is 0.26 nm and most higher fluctuations occur in the loop rather than the structural region. Specifically, the loop region between the residues 110–135 and 217–223 experiences the highest level of fluctuation in each simulation. However, these fluctuations do not affect the structural conformation of the docked complexes (Antony and Vijayan 2015).

It is necessary to analyse the hydrogen bond interaction profile to determine the strength of the hydrogen bond between a ligand and a protein’s binding site. As illustrated in Fig. 6c, each of the complex’s hydrogen bond interaction profiles was examined during the simulation period. The simulation revealed that the AR with Tigogenin and AR with Gitogenin complexes had a maximum of five and a minimum of four H-bonds. The H-Bond interaction profiles showed that Tigogenin and Gitogenin have firmly occupied the catalytic pockets of the AR protein.

The binding energy of the docked complexes was calculated using the last 10 ns of the MD trajectories with the MM/PBSA method implemented using the g_mmpbsa tool of GROMACS (Kumari et al. 2014). The resulting Table 3 provides the binding energy (∆G_bind), van der Waals energy (E_vdw), electrostatic energy (E_elec), polar solvation energy (G_polar) and non-polar solvation energy (G_nonpolar) for each complex.

It is important to note that the binding energy is negatively impacted (stabilised) by E_vdw, E_elec and G_nonpolar while positively affected (destabilised) by G_polar. With both Tigogenin and Gitogenin, the binding energies with AR are below zero, indicating a strong affinity and interaction. In both complexes studied, van der Waals interactions contribute most to the binding energy.