Virtual screening of small molecules against protein targets is an important aspect of drug development, it predict the binding orientation of molecules to the active site of the target. To ascertain the reliability and reproducibility of the docking procedure, the co-crystalized ligand was extracted from the binding site and re-docked back into the binding domain. Figure 3 revealed the superimposition of the re-docked ligand with slight deviation from the original geometry with root mean square deviation (RMSD) of 1.423 Å, which was observed to be less than acceptable 2.00 Å (Omoboyowa et al., 2022b).
Post-docking analysis and binding free energy estimation
Molecular docking is a computational approach that predicts the interaction between receptor-ligand complexes through binding orientation of the ligands in the active site of the protein and ranking score of the interactions (Balogun et al., 2021). To predict potent inhibitors of PDE-5 for the management of erectile dysfunction, small molecules from Aframomum melegueta were screened against the binding site of PDE-5 with the binding affinity estimated using three glide filtering procedures. The binding affinity obtained from the extra precision (XP) mode is presented in Fig. 4. The nine (9) top scored compounds showed scores between − 11.522 kcal/mol to -8.152 kcal/mol compared with the binding affinity of the co-crystalized ligand (sildenafil = -11.872 kcal/mol). Among the hit compounds, 1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate have the best binding affinity of -11.522 kcal/mol which is very close to the co-crystalized ligand. The binding affinity of sildenafil obtained in this study agrees with the value reported by Iwaloye et al. (Iwaloye et al., 2020); Khalid et al. (2022) also showed high binding affinity of sildenafil-PDE-5 complex (-10.01 kcal/mol) using autodock vina.
Table 1
Hydrogen bond interaction of Hits with Phosphodiesterase-5
Compounds
|
Compound Name
|
No of H-bond
|
Interacting residues
|
Co-ligand
|
Sildenafil
|
2
|
GLN 817 (1.82 Å); GLN 817 (4.98 Å)
|
CMPD 1
|
1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate
|
3
|
ASP 724 (2.38 Å); SER 663 (1.69 Å); SER 663 (2.02 Å)
|
CMPD 2
|
1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate
|
3
|
GLN 817 (5.24 Å); GLN 817 (2.49 Å); ASN 662 (2.00 Å)
|
CMPD 3
|
Apigenin
|
1
|
SER 663 (1.74 Å)
|
CMPD 4
|
quercetin-3,7,3,4_-tetramethylether
|
1
|
TYR 612 (2.29 Å)
|
CMPD 5
|
letestuianin A
|
2
|
GLN 775 (2.13 Å); ASP 724 (1.74 Å)
|
CMPD 6
|
kaempferol-3,7,4-trimethylether
|
Nil
|
Nil
|
CMPD 7
|
Buplerol
|
2
|
SER 663 (2.16 Å); GLN 817 (2.82 Å)
|
CMPD 8
|
letestuianin C
|
3
|
GLN 817 (2.00 Å); LEU 725 (2.65 Å); GLU 682 (2.27 Å)
|
CMPD 9
|
5-hydroxy-7-methoxyflavone
|
3
|
ASP 764 (3.39 Å); TYR 612 (2.04 Å); TYR 612 (2.74 Å)
|
E-Pharmacophore Hypothesis and fitness score of hit compounds
Generation of pharmacophore modeling based on small molecules formed an integral component of computational drug development. The hypothesis reports important features that are necessary for the inhibitory activity of active compounds (Miler and Roitberg, 2013). In this study, the ensemble of steric and electronic characters for phosphodiesterase-5 complex with sildenafil was generated using four partitioning features for the identification of the hypothesis model. The hypothesis of the target with sildenafil is shown in Fig. 6. One hydrogen bond acceptor, two hydrogen bond donors, hydrophobic interaction and one aromatic ring were the features selected for the prediction of model with best fitness score. The screening of the hit compounds form A. melegueta was performed according to these features. The fitness scores of the nine (8) hit compounds and sildenafil are presented in Table 2. The co-crystalized ligand was observed to have the highest fitness score (2.605) compared with the hit compounds. Among the hit compounds, 5-hydroxy-7-methoxyflavone possesses the best fitness score of 1.986 while other compounds showed fitness score of less than 1.000 except 1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate with fitness score of 0.754. Fitness score is the quantitative measure of the fitness of ligand into the target with the considering the ligand interaction energies (Natarajan et al., 2015). Therefore, the sildenafil, 5-hydroxy-7-methoxyflavone, Apigenin and quercetin-3,7,3,4-tetramethylether might be predicted to fit well into the binding site of PDE-5 than other bioactive compounds. All the bioactive molecules showed positive fitness score, these showed that they fit into the target.
Table 2
Hit compounds fitness score via pharmacophore model
Compound Label
|
Compound Name
|
Fitness Score
|
Co-ligand
|
Sildenafil
|
2.605
|
CMPD 1
|
1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate
|
1.187
|
CMPD 2
|
1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate
|
0.754
|
CMPD 3
|
Apigenin
|
1.908
|
CMPD 4
|
quercetin-3,7,3,4-tetramethylether
|
1.904
|
CMPD 5
|
letestuianin A
|
1.483
|
CMPD 6
|
kaempferol-3,7,4-trimethylether
|
1.355
|
CMPD 7
|
Buplerol
|
1.247
|
CMPD 8
|
letestuianin C
|
1.312
|
CMPD 9
|
5-hydroxy-7-methoxyflavone
|
1.986
|
AutoQSAR modeling and pIC50 Prediction of hit compounds
Quantitative structure-activity relationship (QSAR) is an automated machine learning tool for validating, building and deploying structural-activity relationship models (Dixion et al., 2016). AutoQSAR is a computational tool that reveals the relationship between the biological potency and structures of small molecules which is crucial in computational drug development (Kwon et al., 2019). In this study, kpls_desc_6 was selected as the best partial least square regression (kpls) analysis predicted from the QSAR model. The kpls was categorized into 25% test and 75% train set (table S1). The predicted activity against observed activity plot is revealed in Fig. 7, it was observed that more train set (blue colour) were present on the plot compared with test set (red colour) which correlate strongly with table S1. The analysis of the selected model parameters showed the standard deviation (SD) as 0.5804, R2 as 0.8833, RMSE as 0.7401and Q2 as 0.7737 as shown in Table 3. The generated autoQSAR model was use to predict the pIC50 of the hit molecules and presented in Table 4, all the hit ligands and sildenafil were observed to have high pIC50 more than 5.00 µM except 1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate (pIC50 = 3.835). The co-crystalized ligand (7.976 µM) showed the highest pIC50 than the hit compounds. The value of negative logarithm of IC50 is pIC50 and is necessary to compare the potency of small molecules at the same molar levels. High values of pIC50 indicate exponentially more potent inhibitors (Stewart and Watson, 1983). Hence, the hit compounds and sildenafil were suggested as potent inhibitor of phosphodiesterase-5.
Table 3
Parameters of the predicted model
Model code
|
S.D
|
R2
|
RMSE
|
Q2
|
kpls_desc_6
|
0.5804
|
0.8833
|
0.7401
|
0.7737
|
Table 4
Predicted pIC50 for the hit molecules via QSAR model
Compound Label
|
Compound Name
|
Predicted pIC50 (µM)
|
Co-ligand
|
Sildenafil
|
7.976
|
CMPD 1
|
1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate
|
5.592
|
CMPD 2
|
1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate
|
3.835
|
CMPD 3
|
Apigenin
|
4.765
|
CMPD 4
|
Quercetin-3,7,3,4_-tetramethylether
|
6.536
|
CMPD 5
|
Letestuianin A
|
6.167
|
CMPD 6
|
Kaempferol-3,7,4-trimethylether
|
5.742
|
CMPD 7
|
Buplerol
|
5.706
|
CMPD 8
|
Letestuianin C
|
5.727
|
CMPD 9
|
5-hydroxy-7-methoxyflavone
|
5.459
|
Pharmacokinetic and drug-likeness evaluation of small molecules as drug candidates is vital in their prediction as novel therapeutic agents. Hence, high percentage of potential drugs failed the stage of clinical trial as a result of poor pharmacokinetic profile with high toxicity rate (Omoboyowa et al., 2021). Due to the time consuming and cost of in vivo toxicity study, in silico model of pharmacokinetic screening remains the fastest and cost-effective approach. In this study, selected pharmacokinetic profiles such as IC50 value for the blockage of HERG K (QPLogHERG), apparent calcium carbonate (QPCaco), madin-Darby canine kidney (MDCK) cell, brain/blood partition coefficient (QPlogBB) and percentage oral absorption were predicted using the Qikprop tool of Schrodinger suite.
The physicochemical profile of the hit compounds presented in Table 5 revealed that, compound 2 [1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate] exhibit poor Caco cell permeability value below 25 nm/s, blood brain partition coefficient below − 3.0, MDCK value and percentage oral absorption. The Caco cell lines remain a model for human intestinal absorption prediction of drug candidates (Van-Breemen and Li, 2005). This blood-brain partition coefficient is vital in the prediction of drug molecules that can penetrate the central nervous system and MDCK cells are utilized in investigating the rate of drug efflux during drug design (Omoboyowa, 2022b). The QPLogHERG value for the co-crystalized ligand and letestuianin A (-6.33 and − 6.31 respectively) were observed to deviate from the reference range while all other molecules were predicted with pharmacokinetic indices within the reference range with high percentage oral absorption.
The result of the predicted drugability of the hit molecules was represented in Table 6, except 1-(3,4-dihydroxy-5-methoxyphenyl)-7-(3,4-dihdroxyphenyl)heptane-3,5-diyldiacetate which violated one Lipinski’s rule of five, The hit molecules obeys all the Lipinski’s rule of five with molecular weight less than 50Da, number of hydrogen bond donor less than 5 but compound 2 has 6, number of hydrogen bond acceptor less than 10 (Lipinski et al., 2001), but the co-crystalized ligand has 11.75. Compounds 1 and 2 were observed to exhibit a high topological surface area of 169.85 and 190.70 respectively.
Table 5
Pharmacokinetic prediction of Hit molecules
Entry Name
|
QPlogHERG
|
QPPCaco
|
QPlogBB
|
QPPMDCK
|
% ORAL ABS
|
Co-Ligand
|
-6.33
|
129.73
|
-0.96
|
60.19
|
77.47
|
Compound 1
|
-5.08
|
48.95
|
-2.90
|
18.97
|
74.61
|
Compound 2
|
-1.58
|
0.21
|
-4.13
|
0.08
|
17.04
|
Compound 3
|
-5.12
|
116.16
|
-1.44
|
48.28
|
73.34
|
Compound 4
|
-5.05
|
504.86
|
-0.74
|
236.33
|
87.86
|
Compound 5
|
-6.31
|
284.41
|
-1.82
|
127.09
|
93.91
|
Compound 6
|
-4.95
|
9906.04
|
-0.38
|
5899.29
|
100.00
|
Compound 7
|
-4.57
|
889.32
|
-0.91
|
435.80
|
100.00
|
Compound 8
|
-5.83
|
201.91
|
-1.81
|
87.76
|
83.64
|
Compound 9
|
-5.28
|
1235.16
|
-0.41
|
621.57
|
100.00
|
Reference values: QplogHERG = below − 5; QPlogBB = -3.0 to 1.2; QPPCaco = < 25 poor, > 500 great |
Table 6
Drug likeness prediction of hit molecules
Entry Name
|
MW
|
HBD
|
HBA
|
TPSA
|
RO5
|
Co-Ligand
|
474.58
|
1.00
|
11.75
|
119.09
|
0
|
Compound 1
|
492.52
|
4.00
|
8.50
|
169.85
|
0
|
Compound 2
|
462.50
|
6.00
|
7.75
|
190.70
|
1
|
Compound 3
|
270.24
|
2.00
|
3.75
|
100.03
|
0
|
Compound 4
|
296.28
|
0.00
|
4.75
|
89.04
|
0
|
Compound 5
|
324.38
|
2.00
|
3.25
|
87.36
|
0
|
Compound 6
|
324.38
|
0.00
|
2.75
|
28.35
|
0
|
Compound 7
|
372.42
|
1.00
|
6.00
|
85.35
|
0
|
Compound 8
|
312.37
|
2.00
|
5.50
|
95.15
|
0
|
Compound 9
|
268.27
|
0.00
|
3.00
|
63.25
|
0
|
Molecular weight (MW), Hydrogen bond donor (HBD), Hydrogen bond acceptor (HBA), topological surface area (TPSA), Rule of five (RO5)
Result of the MD simulation
The result of docking and post-docking study the bioactive compounds from Aframemum melegueta against phosphodiesterase-5 revealed nine (9) compounds as potential inhibitor of the target. From these compounds, 1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate possess the highest binding affinity compared with the standard ligand. Therefore, conformational stability of 1,7-bis(3,4-dihydroxy-5-methoxyphenyl)heptane-3,5-diyldiacetate and sildenafil complex with PDE-5 were validated using 100 ns MD simulation.
Root-mean-square-deviation Analysis
The analysis of the stability of docked complex is a vital index for predicting the mechanism of inhibition of ligand-protein complex during virtual screening. The root-mean-square-deviation (RMSD) gives detail concerning the structural backbone of the target protein (Omboyowa et al., 2023). From Fig. 8, the RMSD plot of the target (Cα) and ligands fit protein over 100 ns MD trajectory revealed initial deviation of less than 1.2 Å of the compound-protein complex over 20 ns between the protein and the ligand with stability observed from 20–100 ns of simulation time. The stability observed in the active complex simulation was similar to the observed stability in the co-ligand-protein complex, although there was no initial deviation in the co-ligand-protein complex but minimal fluctuation of less than 2.0 Å. The deviation that was observed in the active compound-protein complex was acceptable for small globular protein. The stability observed in the RMSD plot of the compound and co-ligand-protein complex is consistent with the close binding affinity observed in the molecular docking and post-docking analysis.
Result of Root-mean-square-fluctuation (RMSF)
The computation of the atomic fluctuation of the position of residues of protein backbone is use to evaluate the dynamic behavior protein–ligand complexes. The positional changes in the amino acid of phosphodiesterase-5 were analyzed using root-mean-square-fluctuation within 100 ns of simulation. Figure 9 shows that, the β-strands and α-helices of the protein crystal structure was observed to oscillate between 0.6–5.4 Å for the lead compound and 0.5–3.5 Å for the co-crystalized ligand (sildenafil). There was less fluctuation in the target’s secondary parameters (β-strands and α-helices) due to their stability compared with the unstructured part of the target. There was high level of fluctuation at the loop region of the protein up to 5.4 Å for the co-ligand but 3.5 Å for the lead compounds.
Protein-Ligand complex interaction mapping
The docking orientation of phosphodiesterase-5-ligand complex was validated for intermolecular contact using the complex contact mapping from individual frame of the trajectory for 100 ns simulation time as shown in Fig. 10. The co-crystalized ligand and compound 1 interacted with the amino acid residue at the binding site of phosphodiesterase-5 with hydrogen bond, hydrophobic bond and water bridges. Glutamine 817 (GLN 817) was observed in both complexes to possess higher interaction fraction of hydrogen bond. Compound-1-PDE-5 complex was predicted to possess more hydrogen bond which might contribute to the stability of the complex than the co-crystalized-PDE-5 complex.