Virtual Screening of Phytochemicals From Medicinal Plants as Promising PDE5 Inhibitors Against Erectile Dysfunction

TL;DR

This study identifies natural compounds from Algerian plants as potential treatments for erectile dysfunction by inhibiting PDE5.

Contribution

The study introduces a computational pipeline to identify and validate phytochemicals as safe and effective PDE5 inhibitors.

Findings

Catechin, ellagic acid, and rosmarinic acid showed strong binding affinities and favorable drug-like properties.

Molecular dynamics simulations confirmed stable interactions of rosmarinic acid and salvinorin A with PDE5.

A QSAR model using molecular descriptors accurately predicted PDE5 inhibition potential with high R² values.

Abstract

This study evaluates bioactive phytochemicals from Algerian medicinal plants as potential phosphodiesterase‐5 (PDE5) inhibitors for the treatment of erectile dysfunction (ED) using an integrated in silico approach. A total of 76 compounds from 48 plant species were screened for drug‐likeness using SwissADME. Overall, 72% of the compounds complied with Lipinski's Rule of Five, indicating favorable oral bioavailability, while toxicity prediction identified 29 non‐toxic candidates. Molecular docking was validated by redocking the co‐crystallized PDE5 ligand (RMSD = 0.264 Å). Ellagic acid (−9.4 kcal·mol−1), rosmarinic acid (−9.2 kcal·mol−1), salvinorin A (−9.2 kcal·mol−1), and catechin (−9.0 kcal·mol−1) exhibited the strongest binding affinities. Molecular dynamics simulations revealed stable hydrogen‐bond interactions for rosmarinic acid, while salvinorin A showed compact and low‐fluctuation behavior. MM‐GBSA analysis confirmed favorable binding free energies for salvinorin A (−26.7 kcal·mol−1) and rosmarinic acid (−23.6 kcal·mol−1). A QSAR model based on docking‐derived pKd values and molecular descriptors showed strong predictive performance using Random Forest regression (R
2
train = 0.91; R
2
CV = 0.87), identifying LogP, molecular weight, and TPSA as key determinants of PDE5 inhibition. Overall, this study highlights catechin and related phytochemicals as promising natural PDE5 inhibitors, supporting their further preclinical evaluation as safer and affordable ED therapies.

Phytochemicals from Algerian medicinal plants were virtually screened as potential phosphodiesterase‐5 (PDE5) inhibitors for erectile dysfunction therapy. Drug‐likeness and toxicity analyses identified catechin, ellagic acid, and rosmarinic acid as safe, bioavailable candidates. Molecular docking and QSAR modeling revealed strong correlations between binding affinity, LogP, and molecular weight. Molecular dynamics simulations confirmed catechin's stable binding to PDE5, supporting its potential as a natural, safe, and affordable PDE5 inhibitor.

Figures

Three‐dimensional structure of phosphodiesterase‐5 (PDB ID: 2H42), the target protein.

Venn diagram showing the drug‐likeness and toxicity profiles of the analyzed plant compounds (total n = 76).

Structural alignment between the redocked pose and the reference (co‐crystallized) ligand within the active site of PDE5. The low RMSD value of 0.264 Å confirms the accuracy of the docking protocol. The original ligand conformation is shown in blue, while the redocked pose is displayed in yellow.

Surface views and active site interactions of PDE5‐ligand complexes. Surface transparency views (left) and active site details (right box) highlighting the enzyme‐ligand interactions. The target enzyme (PDE5) is shown in olive yellow, ligands in blue, hydrogen bonds with their lengths (Å) are represented by dotted yellow lines and interacting amino acid residues linked by hydrogen bonds and hydrophobic interactions are displayed. (A) PDE5–rosmarinic acid complex, (B) PDE5–ellagic acid complex, (C) PDE5–salvinorin A complex, (D) PDE5–catechin complex.

LigPlot+ map of the PDE5‐ligand complexes: Hydrogen bonding (Å) and hydrophobic interactions in the active site. (A) PDE5–rosmarinic acid complex, (B) PDE5–ellagic acid complex, (C) PDE5–salvinorin A complex, (D) PDE5–catechin complex.

HOMO and LUMO orbital distributions of Rosmarinic acid (A), Ellagic acid (B), Salvinorin A (C), and Catechin (D) with the corresponding energy gaps (Eg) values.

Molecular dynamic analysis results of PDE5 and its complexes with selected ligands: (A) Root Mean Square Deviation (RMSD); (B) Root Mean Square Fluctuation (RMSF); (C) Radius of gyration (Rg) profiles; (D) Solvent Accessible Surface Area.

Time evolution of the number of hydrogen bonds between PDE5 and (A) Rosmarinic acid, (B) Ellagic acid, (C) Salvinorin A, and (D) Catechin during 100 ns MD simulations. Hydrogen bonds were monitored along the trajectories and plotted as the instantaneous count versus simulation time (ns). Panels show that rosmarinic acid forms a persistent network of H‐bonds, ellagic acid exhibits larger fluctuations and an increase in H‐bonding toward the end of the run, whereas salvinorin A and catechin form fewer and more intermittent H‐bonds.

Molecular dynamic simulations using the normal mode analysis of the PDE5 protein: (A) Deformability plots, (B) B‐factor plots, (C) Eigenvalues, (D) Variance map, (E) Correlation Matrix (red regions indicate correlated motions, blue regions represent anti‐correlated motions, and white regions show uncorrelated movement), and (F) Elastic Network Model (regions correspond to stiffer interactions and lighter regions to more flexible regions).

QSAR model validation for PDE5 inhibitors using Random Forest regression. (A) Predicted vs. experimental pKd values. (B) Residual distribution. (C) Williams plot (applicability domain). (D) Y‐scrambling validation. (E) Learning curve.

Tables

TABLE 1. Drug‐likeness prediction of selected herbal compounds (n = 76) using SwissADME.

Plant source	Ligand (formula)	Ref.	Molecular weight (MW) (g/mol)	Numberof H‐bond acceptors (HBA)	Numberof H‐bond donors (HBD)	Molar refractivity (MR)	MLOGP	WLOGP	Lipinski's Rule of Five	Ghose filter
Allium sativum	Alliin (C6H11NO3S)	Oyaluna et al. (2024)	177.22	4	2	43.24	−2.88	0.20	Yes; 0 violation	Yes
Allium sativum	Allicin (C6H10OS2)	Oyaluna et al. (2024)	162.27	1	0	45.88	1.18	2.62	Yes; 0 violation	No; 1 violation: < 20 atoms
Allium sativum	Ajoene (C9H14OS3)	Oyaluna et al. (2024)	234.40	1	0	67.41	2.10	3.87	Yes; 0 violation	Yes
Allium sativum	Vinyldithiin (C6H6S2)	Oyaluna et al. (2024)	142.24	1	0	42.60	1.86	2.97	Yes; 0 violation	No; 2 violations: MW < 160, < 20 atoms
Arbutus unedo	Arbutin (C12H16O7)	Bhalla et al. (2022)	272.25	7	5	62.61	−1.49	−1.43	Yes; 0 violation	No; 1 violation: WLOGP<‐0.4
Calendula officinalis	Calenduladiol (C30H50O2)	Vella et al. (2024)	442.72	2	2	136.30	6.00	7.00	Yes; 1 violation: MLOGP > 4.15	No; 3 violations: WLOGP > 5.6, MR > 130, > 70 atoms
Calendula officinalis	Officinoside A (C19H32O8)	Vella et al. (2024)	388.45	8	5	95.47	−0.85	−0.54	Yes; 0 violation	No; 1 violation: WLOGP<‐0.4
Capsicum annuum	Capsaicin (C18H27NO3)	Kuzma et al. (2014)	305.41	3	2	90.52	2.69	3.64	Yes; 0 violation	Yes
Capsicum annuum	Dihydrocapsaicin (C18H29NO3)	Kuzma et al. (2014)	307.43	3	2	90.99	2.78	3.86	Yes; 0 violation	Yes
Artemisia herba‐alba	Herbolide A (C17H24O4)	Mohammed et al. (2021)	292.37	4	0	81.23	2.83	3.17	Yes; 0 violation	Yes
Artemisia herba‐alba	Arteanoflavone (C19H18O8)	Mohammed et al. (2021)	374.34	8	2	97.93	−0.12	2.91	Yes; 0 violation	Yes
Capsella bursa‐pastoris	Diosmetin (C16H12O6)	Zhou et al. (2023)	300.26	6	3	80.48	0.22	2.59	Yes; 0 violation	Yes
Lepidium sativum	Lepidine B (C20H18N4O2)	Painuli et al. (2022)	346.38	4	3	98.18	1.33	3.81	Yes; 0 violation	Yes
Lepidium sativum	Semilepidinoside B (C17H22N2O7)	Painuli et al. (2022)	366.37	8	5	88.68	−1.58	−0.81	Yes; 0 violation	No; 1 violation: WLOGP<‐0.4
Lepidium sativum	Syringoylglucose (C15H20O10)	Painuli et al. (2022)	360.31	10	6	81.69	−2.95	−2.40	Yes; 1 violation: NH or OH > 5	No; 1 violation: WLOGP < 0.4
Ericacinerea	Chrysanthemin (C21H21O11 +)	Aires and Carvalho (2017)	449.38	11	8	108.29	−1.76	0.38	No; 2 violations: N or O > 10, NH or OH > 5	Yes
Bryonia dioica	Bryonioside A (C36H60O9)	Ukiya et al. (2002)	636.86	9	6	172.46	1.66	3.50	No; 2 violations: MW > 500, NH or OH > 5	No; 3 violations: MW > 480, MR > 130, > 70 #atoms
Lavandula officinalis	Rosmarinate (C18H16O8)	Mansour et al. (2021)	360.31	8	5	91.40	0.90	1.65	Yes; 0 violation	Yes
Anthemis nobilis	Apigenin (C15H10O5)	Mouissi et al. (2022)	270.24	5	3	73.99	0.52	2.58	Yes; 0 violation	Yes
Daucus carota	Carotamine (C14H13N3O)	Boadi et al. (2021)	271.27	3	4	76.42	0.06	1.63	Yes; 0 violation	Yes
Cuminum cyminum	Cuminoside A (C21H30O8)	Karik et al. (2021)	410.46	8	4	101.80	0.21	0.18	Yes; 0 violation	Yes
Cuminum cyminum	Cumene (C9H12)	Karik et al. (2021)	120.19	0	0	41.02	4.17	2.81	Yes; 1 violation: MLOGP > 4.15	No; 1 violation: MW < 160
Melissa officinalis	Melitric acid A (C27H22O12)	Zam et al. (2022)	538.46	12	7	136.10	0.77	2.75	No; 3 violations: MW > 500, N or O > 10, NH or OH > 5	No; 2 violations: MW > 480, MR > 130
Foeniculum vulgare	Avicularin (C20H18O11)	Rafieian et al. (2024)	434.35	11	7	104.19	−2.06	0.10	No; 2 violations: N or O > 10, NH or OH > 5	Yes
Origanum vulgare	Amburoside A (C20H22O10)	Puvača et al. (2022)	422.38	10	6	100.40	−0.83	−0.52	Yes; 1 violation: NH or OH > 5	No; 1 violation: WLOGP<‐0.4
Zingiber officinalis	Zingiberene (C15H24)	Yang et al. (2024)	204.35	0	0	70.68	4.53	4.89	Yes; 1 violation: MLOGP > 4.15	Yes
Zingiber officinalis	Gingerdiol (C17H28O4)	Yang et al. (2024)	296.40	4	3	85.51	2.23	3.03	Yes; 0 violation	Yes
Zingiber officinalis	(S)‐6‐gingerol (C17H26O4)	Yang et al. (2024)	294.39 g/mol	4	2	84.55	2.14	3.23	Yes; 0 violation	Yes
Equisetum arvense	Herbacitrin (C21H20O12)	Boeing et al. (2021)	464.38	12	8	110.16	−2.59	−0.54	No; 2 violations: N or O > 10, NH or OH > 5	No; 1 violation: WLOGP<‐0.4
Lactuca virosa	Scorzoside (C21H30O8)	Abdel Bar et al. (2023)	410.46	8	4	101.76	0.21	−0.11	Yes; 0 violation	Yes
Rosmarinus officinalis	Rosmarinic acid (C18H16O8)	Puvača et al. (2022)	360.31	8	5	91.40	0.90	1.65	Yes; 0 violation	Yes
Rosmarinus officinalis	Carnosic acid (C20H28O4)	Puvača et al. (2022)	332.43	4	3	95.43	3.25	4.32	Yes; 0 violation	Yes
Rosmarinus officinalis	Carnosol (C20H26O4)	Puvača et al. (2022)	330.42	4	2	92.83	3.25	3.96	Yes; 0 violation	Yes
Marrubium vulgare	Marrubin (C20H28O4)	Zahedifar and Najafian (2023)	332.43	4	1	91.40	2.76	3.72	Yes; 0 violation	Yes
Marrubium vulgare	Peregrinol (C20H36O2)	Zahedifar and Najafian (2023)	308.50	2	2	95.39	4.70	3.93	Yes; 0 violation	Yes
Satureja hortensis	Carvacrol (C10H14O)	Ejaz et al. (2023)	150.22	1	1	48.01	2.76	2.82	Yes; 0 violation	No; 1 violation: MW < 160
Salvia officinalis	Safficinolide (C20H24O5)	Ejaz et al. (2023)	344.40	5	1	94.14	2.25	3.51	Yes; 0 violation	Yes
Salvia officinalis	Columbaridione (C20H24O5)	Ejaz et al. (2023)	344.40	5	1	91.59	1.80	3.04	Yes; 0 violation	Yes
Salvia officinalis	Sageone (C19H24O3)	Ejaz et al. (2023)	300.39	3	2	89.49	2.98	4.31	Yes; 0 violation	Yes
Myrthus communis	Sideroxylin (C18H16O5)	Medda and Mulas (2021)	312.32	5	2	88.39	1.25	3.50	Yes; 0 violation	Yes
Vinca minor	Vincamine (C21H26N2O3)	Neculai et al. (2023)	354.44	4	1	103.74	2.62	2.14	Yes; 0 violation	Yes
Vinca minor	Vincaminoreine (C22H30N2O2)	Neculai et al. (2023)	354.49	3	0	109.62	3.26	3.49	Yes; 0 violation	Yes
Nigella sativa	Nigeglanine (C12H14N2O)	Manoharan et al. (2021)	202.25	1	0	61.19	1.69	1.91	Yes; 0 violation	Yes
Nigella sativa	Nigellidine (C18H18N2O2)	Manoharan et al. (2021)	294.35	2	1	88.65	2.39	3.28	Yes; 0 violation	Yes
Urtica dioica	Enterofuran (C18H20O3)	Đurović et al. (2024)	284.35	3	2	82.95	2.66	3.15	Yes; 0 violation	Yes
Urtica dioica	Neoolivil (C20H24O7)	Đurović et al. (2024)	376.40	7	4	98.25	0.37	1.50	Yes; 0 violation	Yes
Petroselinum sativum	Isopimpinellin (C13H10O5)	Tahraoui et al. (2007)	246.22	5	2	65.24	0.89	2.56	Yes; 0 violation	Yes
Ajuga iva L.	Clerodane (C20H38)	Djordjevic et al. (2018)	278.52	0	0	93.51	6.69	6.99	Yes; 1 violation: MLOGP > 4.15	No; 1 violation: WLOGP > 5.6
Teucrium polium L.	Luteolin (C15H10O6)	Djordjevic et al. (2018)	286.24	6	4	76.01	−0.03	2.28	Yes; 0 violation	Yes
Teucrium polium L.	Apigenin (C15H10O5)	Djordjevic et al. (2018)	270.24	5	3	73.99	0.52	2.58	Yes; 0 violation	Yes
Artemisia herba alba Asso	Carvone (C10H14O)	Tahraoui et al. (2007)	150.22	1	0	47.32	2.10	2.49	Yes; 0 violation	No; 1 violation: MW < 160
Artemisia herba alba Asso	Piperitone (C10H16O)	Tahraoui et al. (2007)	152.23	1	0	47.80	2.20	2.57	Yes; 0 violation	No; 1 violation: MW < 160
Malva sylvestris L.	2‐Methoxy‐4‐vinylphenol (C9H10O2)	Dushkin et al. (2014)	150.17	2	1	45.05	1.71	1.93	Yes; 0 violation	No; 1 violation: MW < 160
Glycyrrhiza glabra L.	Glycyrrhizin (C42H62O16)	Dushkin et al. (2014)	822.93	16	8	202.84	0.02	2.25	No; 3 violations: MW > 500, N or O > 10, NH or OH > 5	No; 3 violations: MW > 480, MR > 130, > 70 #atoms
Punica granatum L.	Punicalagin (C48H28O30)	Dushkin et al. (2014)	1084.72	30	17	250.86	−3.29	1.96	No; 3 violations: MW > 500, N or O > 10, NH or OH > 5	No; 3 violations: MW > 480, MR > 130, > 70 #atoms
Punica granatum L.	Ellagic acid (C14H6O8)	Dushkin et al. (2014)	302.19	8	4	75.31	0.14	1.31	Yes; 0 violation	Yes
Juniperus oxyedrus	Alpha‐pinene (C10H16)	Zouaoui et al. (2024)	136.23	0	0	45.22	3.00	4.29	Yes; 1 violation: MLOGP > 4.15	No; 1 violation: MW < 160
Salvia officinalis L.	Carnosol (C20H26O4)	Tahraoui et al. (2007)	330.42	4	2	92.83	3.25	3.96	Yes; 0 violation	Yes
Salvia officinalis L.	Salvinorin A (C23H28O8)	Tahraoui et al. (2007)	432.46	8	0	107.39	1.56	2.68	Yes; 0 violation	Yes
Salvia officinalis L.	Rosmarinic acid (C18H16O8)	Tahraoui et al. (2007)	360.31	8	5	91.40	0.90	1.65	Yes; 0 violation	Yes
Quercus ilex	Catechin (C15H14O6)	Sánchez‐Gutiérrez et al. (2022)	290.27	6	5	74.33	1.22	0.24	Yes; 0 violation	Yes
Quercus ilex	Quercetin (C15H10O7)	Sánchez‐Gutiérrez et al. (2022)	302.24	7	5	78.03	−0.56	1.99	Yes; 0 violation	Yes
Quercus ilex	Gallic acid (C7H6O5)	Sánchez‐Gutiérrez et al. (2022)	170.12	5	4	39.47	−0.16	0.50	Yes; 0 violation	No; 2 violations: MR < 40, < 20 #atoms
Quercus ilex	Kaempferol (C15H10O6)	Sánchez‐Gutiérrez et al. (2022)	286.24	6	4	76.01	−0.03	2.28	Yes; 0 violation	Yes
Artemisia campestris	Limonene (C10H16)	Kemal et al. (2023)	136.23	0	0	47.12	3.27	3.31	Yes; 0 violation	No; 1 violation: MW < 160
Artemisia campestris	Beta‐pinene (C10H16)	Kemal et al. (2023)	136.23	0	0	45.22	4.29	3.00	Yes; 1 violation: MLOGP > 4.15	No; 1 violation: MW < 160
Artemisia campestris	Myrcene (C10H17Cl)	Kemal et al. (2023)	136.23	0	0	48.76	3.56	3.48	Yes; 0 violation	No; 1 violation: MW < 160
Lawsonia alba L	Naphthoquinone (C10H6O2)	Mustapha et al. (2024)	158.15	2	0	44.24	0.91	1.62	Yes; 0 violation	No; 2 violations: MW < 160, < 20 #atoms
Pinus halepensis	Alpha‐humulene (C15H24)	Fekih et al. (2014)	204.35	0	0	4.53	5.04		Yes; 1 violation: MLOGP > 4.15	Yes
Pinus halepensis	Phenyl 3‐methylbutanoate (C11H14O2)	Fekih et al. (2014)	178.23	2	0	52.36	2.84	2.64	Yes; 0 violation	Yes
Pinus halepensis	Sabinene (C10H16)	Fekih et al. (2014)	136.23	0	0	45.22	4.29	3.00	Yes; 1 violation: MLOGP > 4.15	No; 1 violation: MW < 160
Opuntia ficus indica L.	Betaxanthin (C18H18N2O6 −2)	Giraldo‐Silva et al. (2023)	358.35	7	2	94.99	−2.17	−1.69	Yes; 0 violation	No; 1 violation: WLOGP<‐0.4
Anthemis nobilis L	Tonghaosu (C7H10O2)	De Cicco et al. (2023)	200.23	2	0	57.96	2.33	2.15	Yes; 0 violation	Yes
Linum usitatissimum	Lignan (C20H20O6)	Chhillar et al. (2021)	458.50	8	0	121.76	1.58	3.87	Yes; 0 violation	Yes
Foeniculum vulgare Mill.	Linalool (C10H18O)	Mehra et al. (2021)	154.25	1	1	50.44	2.59	2.67	Yes; 0 violation	No; 1 violation: MW < 160
Foeniculum vulgare Mill.	Estragole (C10H12O)	Mehra et al. (2021)	148.20	1	0	47.04	2.67	2.42	Yes; 0 violation	No; 1 violation: MW < 160
Cassia senna L	Anthraquinone (C14H8O2)	Al‐Tameme et al. (2019)	208.21	2	0	59.75	1.86	2.46	Yes; 0 violation	Yes
Pistacia lentiscus L	Beta‐gurjunene (C15H24)	Al‐Tameme et al. (2019)	204.35	0	0	67.14	5.65	4.27	Yes; 1 violation: MLOGP > 4.15	Yes
Globularia alypum L	Cinnamic acid (C9H8O2)	Hajji et al. (2022)	148.16	2	1	43.11	1.90	1.68	Yes; 0 violation	No; 2 violations: MW < 160, < 20 #atoms
Origanum majorana L.	Carvacrol (C10H14O)	Vinciguerra et al. (2019)	150.22	1	1	48.01	2.76	2.82	Yes; 0 violation	No; 1 violation: MW < 160

TABLE 2. Toxicity risk prediction of the plant active compounds (n = 44) using OSIRIS Property Explorer. The toxicity risks are categorized as Green (no risk) and Red (potential risk).

Ligand	Plant source	Mutagenic	Tumorigenic risk	Irritant risk	Reproductive risk
Alliin (C6H11NO3S)	Allium sativum	Green	Green	Green	Green
Ajoene (C9H14OS3)	Allium sativum	Green	Green	Green	Green
Capsaicin (C18H27NO3)	Capsicum annuum	Green	Green	Green	Green
Dihydrocapsaicin (C18H29NO3)	Capsicum annuum	Green	Green	Green	Green
Herbolide A (C17H24O4)	Artemisia herba‐alba	Green	Green	Green	Green
Arteanoflavone (C19H18O8)	Artemisia herba‐alba	Red	Red	Green	Green
Diosmetin (C16H12O6)	Capsella bursa‐pastoris	Green	Green	Green	Green
Lepidine B (C20H18N4O2)	Lepidium sativum	Green	Green	Green	Green
Rosmarinate C18H16O8	Lavandula officinalis	Green	Green	Green	Green
Apigenin C15H10O5	Anthemis nobilis	Red	Green	Green	Green
Carotamine (C14H13N3O3)	Daucus carota	Green	Red	Red	Green
Cuminoside A (C21H30O8)	Cuminum cyminum	Green	Green	Green	Green
Zingiberene (C15H24)	Zingiber officinalis	Red	Green	Red	Red
Gingerdiol (C17H28O4)	Zingiber officinalis	Green	Green	Green	Green
(S)‐6‐gingerol (C17H26O4)	Zingiber officinalis	Green	Green	Green	Green
Scorzoside (C21H30O8)	Lactuca virosa	Green	Green	Green	Green
Rosmarinic acid (C18H16O8)	Rosmarinus officinalis	Green	Green	Green	Green
Carnosic acid (C20H28O4)	Rosmarinus officinalis	Green	Green	Green	Green
Carnosol (C20H26O4)	Rosmarinus officinalis	Green	Green	Green	Green
Marrubin (C20H28O4)	Marrubium vulgare	Green	Green	Green	Green
Peregrinol (C20H36O2)	Marrubium vulgare	Green	Green	Green	Green
Safficinolide (C20H24O5)	Salvia officinalis	Green	Green	Red	Green
Columbaridione (C20H24O5)	Salvia officinalis	Green	Green	Red	Green
Sageone (C19H24O3)	Salvia officinalis	Green	Green	Green	Red
Sideroxylin (C18H16O5)	Myrthus communis	Green	Red	Green	Green
Vincamine (C21H26N2O3)	Vinca minor	Green	Green	Green	Green
Vincaminoreine (C22H30N2O2)	Vinca minor	Green	Green	Green	Green
Nigeglanine (C12H14N2O)	Nigella sativa	Green	Green	Green	Green
Nigellidine (C18H18N2O2)	Nigella sativa	Green	Green	Green	Green
Enterofuran (C18H20O3)	Urtica dioica	Green	Green	Green	Green
(+)‐Neoolivil (C20H24O7)	Urtica dioica	Green	Green	Green	Green
Isopimpinellin (C13H10O5)	Petroselinum sativum	Green	Green	Green	Green
Luteolin (C15H10O6)	Teucrium polium LE.	Green	Green	Green	Green
Apigenin (C15H10O5)	Teucrium polium LE.	Red	Green	Green	Green
Ellagic acid (C14H6O8)	Punica granatum L.	Green	Green	Green	Green
Carnosol (C20H26O4)	Salvia officinalis L.	Green	Green	Green	Green
Salvinorin A (C23H28O8)	Salvia officinalis L.	Green	Green	Green	Green
Rosmarinic acid (C18H16O8)	Salvia officinalis L.	Green	Green	Green	Green
Catechin (C15H14O6)	Quercus ilex	Green	Green	Green	Green
Quercetin (C15H10O7)	Quercus ilex	Red	Red	Green	Green
Kaempferol (C15H10O6)	Quercus ilex	Red	Green	Green	Green
Alpha‐humulene (C15H24)	Pinus halepensis	Green	Green	Green	Green
Phenyl 3‐methylbutanoate (C11H14O2)	Pinus halepensis	Green	Green	Red	Green
Tonghaosu (C7H10O2)	Anthemis nobilis L	Green	Red	Green	Green
Lignan (C20H22O6)	Linum usitatissimum	Green	Green	Green	Red
Anthraquinone (C14H8O2)	Cassia senna L	Red	Red	Red	Green
Beta‐gurjunene (C15H24)	Pistacia lentiscus L	Red	Green	Green	Green

TABLE 3. Binding energies (kcal/mol) of plant active compounds and of standard drugs.

Compound	Compound ID	Estimated free energy of binding (kcal/mol)	Compound	Compound ID	Estimated free energy of binding (kcal/mol)
Alliin (C6H11NO3S)	CID: 87310	−4.4	Peregrinol (C20H36O2)	CID: 7092583	−7.8
Ajoene (C9H14OS3)	CID: 5386591	−4.0	Vincamine (C21H26N2O3)	CID: 15376	−6.6
Capsaicin (C18H27NO3)	CID: 1548943	−7.4	Vincaminoreine (C22H30N2O2)	CID: 12444831	−7.9
Dihydrocapsaicin (C18H29NO3)	CID: 107982	−5.2	Nigeglanine (C12H14N2O)	CID: 12116700	−7.5
Herbolide A (C17H24O4)	CID: 90474440	−7.9	Nigellidine (C18H18N2O2)	CID: 136828302	−8.8
Diosmetin (C16H12O6)	CID: 5281612	−9.2	Enterofuran (C18H20O3)	CID: 10107942	−6.6
Lepidine B (C20H18N4O2)	CID: 100927764	−8.4	(+)‐Neoolivil (C20H24O7)	CID: 9976812	−8.6
Rosmarinate C18H16O8	CID: 639655	−9.1	Isopimpinellin (C13H10O5)	CID: 68079	−7.0
Cuminoside A (C21H30O8)	CID: 101252065	−8.8	Luteolin (C15H10O6)	CID: 5280445	−9.4
Gingerdiol (C17H28O4)	CID: 11369949	−7.3	Ellagic acid (C14H6O8)	CID: 5281855	−9.4
(S)‐6‐gingerol (C17H26O4)	CID: 44559528	−7.1	Salvinorin A (C23H28O8)	CID: 128563	−9.2
Scorzoside (C21H30O8)	CID: 21633213	−9.4	Catechin (C15H14O6)	CID: 73160	−9.0
Rosmarinic acid (C18H16O8)	CID: 5281792	−9.2	Humulene (C15H24)	CID: 24798693	−7.0
Carnosic acid (C20H28O4)	CID: 65126	−8.8	Sildenafil (C22H30N6O4S)	CID: 135398744	−9.5
Carnosol (C20H26O4)	CID: 442009	−8.4	Vardenafil (C23H32N6O4S)	CID: 135400189	−9.4
Marrubin (C20H28O4)	CID: 73401	−8.6	Tadalafil (C22H19N3O4)	CID: 110635	−9.3

TABLE 4. LigPlot+ results: Hydrogen bonding and hydrophobic interactions with enzyme residues.

Plant active compounds	Compound ID	Hydrogen bonding (length A°)	Hydrophobic interactions
Rosmarinic acid (C18H16O8)	CID: 5281792	Gln817 (2.88 and 3.15), His613 (2.86), Asp654 (2.82), Asp764 (2.85), His617 (3.14), and Gln775 (3.12)	Ala767, Phe820, Val782, Asn661, Thr723, His657, Glu682, Leu765, and Ile768
Ellagic acid (C14H6O8)	CID: 5281855	His617 (3.32), Gln817 (2.77, 2.86, and 3.13), Tyr612 (2.72), and Asp764 (2.74)	Leu765, His613, Phe786, Phe820, Val782, and Ala767
Salvinorin A (C23H28O8)	CID: 128563	/	Val782, Leu804, Gln817, Met816, Gly819, Ala823, Ile665, Ile824, Leu725, Phe786, and Phe820
Catechin (C15H14O6)	CID: 73160	Gln817 (2.70), His613 (3.01), Asp764 (2.93 and 2.93), Thr723 (2.83)	Ala767, Ile768, Val782, Leu765, Leu725, Phe786, and Phe820

TABLE 5. HOMO, LUMO, and Eg values of the selected compounds.

Plant active compounds	HOMO (eV)	LUMO (eV)	Eg (eV)
Rosmarinic acid (C18H16O8)	−0.30478	0.06971	0.37449
Ellagic acid (C14H6O8)	−0.32226	0.04282	0.36508
Salvinorin A (C23H28O8)	−0.31943	0.13202	0.45145
Catechin (C15H14O6)	−0.30447	0.13625	0.44072

TABLE 6. MM‐GBSA binding energetics for PDE5–ligand complexes (kcal mol−1).

	Rosmarinic acid (C18H16O8)	Ellagic acid (C14H6O8)	Catechin (C15H14O6)	Salvinorin A (C23H28O8)
van der Waal Energy	−31.39	−5.64	−27.47	−50.52
Electrostatic energy	−41.73	−24.78	−15.25	−7.83
ΔENPOLAR	−4.21	−1.69	−3.33	−4.74 cd
ΔGSOLV	49.57	17.74	29.62	31.61
ΔTOTAL	−23.55	−12.68	−13.10	−26.74
Binding energy	−3.76	52.13	3.30	−18.70

Full text

Introduction

1

Erectile dysfunction (ED) is a prevalent condition that significantly impacts male reproductive and psychological health, with incidence increasing with age. In severe cases, ED may contribute to infertility and reduced quality of life (Agarwal et al. 2015). Its etiology is multifactorial, including lifestyle factors (e.g., smoking, stress, obesity) and medical conditions such as diabetes and hypertension (Alahmar 2019; Shamloul and Ghanem 2013). Clinical manifestations range from ejaculatory disorders and reduced libido to the complete inability to maintain an erection, making ED a major public health concern. Currently, phosphodiesterase‐5 (PDE5) inhibitors such as sildenafil (Viagra), tadalafil, and vardenafil represent the first‐line therapy for ED. These drugs function by inhibiting PDE5‐mediated degradation of cyclic guanosine monophosphate (cGMP), thereby promoting penile vasodilation and erection (Guay et al. 2001). However, their clinical use is limited by side effects such as headache, back pain, and visual disturbances, resulting from off‐target interactions with other PDE isoforms (Foresta et al. 2008). Thus, there is growing demand for safer, long‐term alternatives, particularly from natural sources.

Herbal medicine has gained global attention as a safer alternative to synthetic drugs, attributed to its antioxidant, anti‐inflammatory, and phytoestrogenic properties. An estimated 60%–80% of the global population, particularly in developing countries, relies on plant‐based remedies (Thomford et al. 2015). Several studies report that plant extracts can significantly improve erectile function, in some cases by up to 96% (Jaradat and Zaid 2019). Additionally, phytotherapies often reduce adverse effects by 40%–60% and are up to five times more affordable compared to conventional PDE5 inhibitors (Dehelean et al. 2021). Nevertheless, drug discovery from natural sources faces challenges such as high failure rates, limited efficacy, and toxicity issues (Sang et al. 2018). Computational approaches, including molecular docking, QSAR modeling, and molecular dynamics (MD) simulations, have emerged as essential tools in modern drug discovery, enabling rapid and cost‐effective screening of phytochemicals for desirable pharmacokinetic and binding properties (Verma et al. 2021). This study aims to systematically screen bioactive compounds from Algerian medicinal plants for PDE5 inhibition. These plants, deeply rooted in traditional medicine, represent an untapped pharmacological resource. Specifically, we evaluated drug‐likeness, toxicity, molecular docking, QSAR modeling, and MD simulations to identify promising PDE5 inhibitors. By integrating computational screening with traditional knowledge, this study highlights the potential of Algerian phytochemicals to provide safer, natural alternatives to conventional PDE5 inhibitors (Rahman et al. 2024).

Experimental Section

2

Selection of Potential Ligands

2.1

A collection of 76 natural compounds from 48 medicinal plants, previously identified for their potential in treating human infertility, was selected for this study (Table 1). These compounds were sourced from previous studies conducted in Algeria (Adamu et al. 2012, 2014; Djerrou et al. 2022). The three‐dimensional (3D) structures of these compounds were retrieved from the PubChem database (Li et al. 2010), and their SMILES notation was generated to serve as input for the computational tools. The drug‐likeness and toxicity risk of these compounds were assessed to ensure the safety and viability of potential therapeutic candidates.

Drug‐Likeness Prediction

2.2

The drug‐likeness of the selected compounds was predicted using SwissADME (http://www.swissadme.ch/index.php), which evaluated critical physicochemical properties, such as Lipinski's Rule of Five, the Ghose filter, and additional bioavailability parameters, including polar surface area and rotatable bonds. The SMILES notations of the compounds were submitted, and their physicochemical and pharmacokinetic properties were assessed to determine their compliance with Lipinski's Rule of Five (Lipinski et al. 1997). The Ghose filter was also applied, considering factors such as log P between −0.4 and +5.6, molar refractivity between 40 and 130, molecular weight (MW) between 180 and 480, and atom count of 20–70 (Ghose et al. 1999). These predictions were crucial for selecting compounds suitable for drug development. The output data was exported for subsequent molecular docking and toxicity analyses.

Toxicity Risk Assessment

2.3

Toxicity was evaluated using OSIRIS Property Explorer (http://www.organic‐chemistry.org/prog/peo). This tool predicts mutagenicity, tumorigenicity, irritant effects, and reproductive toxicity of compounds based on their molecular structures. The outputs were color‐coded to indicate toxicity risks (red = high risk, green = low risk) (Preethi et al. 2021). Compounds with acceptable toxicity profiles were selected for further analysis.

Target Protein Selection and Molecular Docking

2.4

For the molecular docking studies, phosphodiesterase‐5 (PDE5) (PDB ID: 2H42) was selected as the target protein (Figure 1) because it plays a crucial role in erectile dysfunction (ED) pharmacological management. The protein structure was retrieved from the RCSB Protein Data Bank, and preprocessing was performed by removing water molecules and heteroatoms using PyMOL v2.5.2 (Yuan et al. 2017). The cleaned structure was energy‐minimized using Open Babel. Docking simulations were carried out with PyRx v0.8, using AutoDock Vina as the docking engine. A grid box dimension of 30 × 30 × 30 Å^3^ was applied to encompass the PDE5 active site, and the exhaustiveness was set to 8 to ensure a thorough exploration of binding modes. The ligands were selected based on their compliance with Lipinski's Rule of Five (≤ 1 violation) and the Ghose filter (no violations), and their 3D structures were obtained from PubChem. Enzyme preparation was performed using Molegro v2.5. To validate the docking protocol, the redocking approach was applied, in which the co‐crystallized ligand from the PDE5 crystal structure was docked back into its binding site. The root‐mean‐square deviation (RMSD) between the redocked and experimental ligand conformations was calculated, and an RMSD value of < 2 Å confirmed the reliability of the docking procedure (Ness et al. 2000). Docking results were visualized with PyMOL, and interaction maps were generated using LigPlot+ v2.2.7 to analyze hydrogen bonding and hydrophobic interactions. Ligands with the lowest binding energies were prioritized for further molecular dynamics (MD) simulations.

Three‐dimensional structure of phosphodiesterase‐5 (PDB ID: 2H42), the target protein.

Density Functional Theory (DFT)

2.5

The quantum chemical properties of the selected ligands were calculated using Gaussian 09 W at the B3LYP/6‐31G level of theory (Borah and Saikia 2024). The ligand geometries were optimized without symmetry constraints, and key electronic properties, including the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO), were calculated. The HOMO–LUMO energy gap (Eg) was used to assess electronic stability. Molecular electrostatic potential maps were generated to visualize electron‐rich and electron‐deficient regions, which are crucial for predicting ligand–receptor interactions.

Molecular Dynamics (MD) Simulation

2.6

MD simulations were performed using GROMACS 2023‐GPU with the CHARMM36m force field to validate docking results and assess the dynamic stability of PDE5–ligand complexes (Ikhou et al. 2024). The selected complexes, based on docking free energy scores, were solvated in a cubic box with TIP3P water molecules and neutralized with Na + and Cl − ions at a concentration of 150 mM. Energy minimization was carried out using the steepest descent method (convergence threshold: 1000 kJ/mol/nm), followed by equilibration under NVT and NPT conditions for 2 ns each, at 300 K and 1 bar. Production simulations were run for 100 ns. Stability of the complexes was assessed by computing RMSD, RMSF, and radius of gyration (Rg). Structural flexibility was further investigated using Normal Mode Analysis (NMA) on the iMODS server (http://imods.chaconlab.org), which provided insights into motion patterns and conformational changes. The MD protocol was validated by redocking experiments and the consistency of stability metrics across replicates.

QSAR‐Based Structure–Activity Relationship Modeling

2.7

To complement docking and MD studies, a 2D‐QSAR (Quantitative Structure–Activity Relationship) model was developed for the selected compounds. Molecular descriptors were calculated from SMILES notations using RDKit in Python (Pinheiro et al. 2020). The descriptors were normalized and filtered to remove highly correlated variables (r > 0.90). The biological activity endpoint was represented by estimated binding affinities obtained from docking studies. Feature selection was performed using LASSO regression, and the predictive model was built with Random Forest regression (n = 500 trees). Model robustness was evaluated through 10‐fold cross‐validation, with performance assessed using R ^2^ and RMSE values. The QSAR model helped rationalize the structure–activity relationship of the phytochemicals, providing additional support for docking and MD results.

Results and Discussion

3

Drug‐Likeness Prediction and Toxicity Risk Assessment

3.1

The drug‐likeness prediction and toxicity risk assessment of natural compounds derived from Algerian medicinal plants traditionally used for infertility treatment were carried out using a two‐step computational approach. A Venn diagram (Figure 2) was used to summarize the outcomes, allowing visualization of the overlap between compounds with desirable drug‐like properties and those free from significant toxicity risks. Among the 76 compounds initially extracted from 48 medicinal plants (Table 1), 44 were identified as having favorable drug‐likeness properties according to Lipinski's Rule of Five and the Ghose filter. Based on Lipinski's criteria, which evaluate molecular weight (MW), lipophilicity (logP), hydrogen bond donors (HBD), and hydrogen bond acceptors (HBA), 47 compounds (61.9%) satisfied the necessary pharmacokinetic requirements for drug development. For example, sideroxylin, rosmarinic acid, and sageone had MW values of 312.5, 328.3, and 301.8 g/mol, respectively, all within the ideal range (MW < 500 g/mol). Their logP values ranged from 2.5 to 4.5, indicating moderate lipophilicity, which is crucial for absorption (Cetin 2021). Notably, kaempferol and isorhamnetin exhibited a balanced combination of favorable HBD and HBA values that aligned well with drug‐likeness standards.

Venn diagram showing the drug‐likeness and toxicity profiles of the analyzed plant compounds (total n = 76).

Toxicity risk was assessed using the OSIRIS Property Explorer (Table 2), which evaluates mutagenic, tumorigenic, irritant, and reproductive risks. This analysis revealed that 29 compounds were free from significant toxicity risks. Catechin, epicatechin, quercetin, and rutin, for instance, showed no signs of mutagenic, tumorigenic, or reproductive toxicity, highlighting their potential as therapeutic agents. Conversely, 18 compounds, including luteolin, apigenin, and kaempferol, exhibited mutagenic or tumorigenic risks based on high toxicity scores (> 0.5) (Fossatelli et al. 2023). Such compounds may require structural modification to reduce toxicity before therapeutic use (Cetin et al. 2025). A comparison between drug‐likeness and toxicity results revealed a consistent trend: compounds with higher MW and more complex structures, such as luteolin and apigenin, tended to display greater toxicity risks. This observation agrees with previous studies (Skanes et al. 2021; Oselusi et al. 2021) showing that larger, structurally complex natural compounds often exhibit increased mutagenic or tumorigenic potential. Conversely, simpler structures with moderate MW, such as quercetin, kaempferol, and myricetin, met drug‐likeness criteria and displayed minimal toxicity, corroborating findings by Anjum et al. (2021) and Babiaka et al. (2023). These results emphasize the balance between molecular simplicity and safety (Cetin 2022).

Molecular Docking and Binding Energies

3.2

Molecular docking was performed to evaluate the interactions of selected plant‐derived compounds with phosphodiesterase‐5 (PDE5).

The estimated binding energies are summarized in Table 3, which compares plant compounds with standard PDE5 inhibitors. Ellagic acid from Punica granatum L. and luteolin from Teucrium polium L. exhibited the highest binding affinities, with docking scores of −9.4 kcal/mol, comparable to standard drugs such as sildenafil (−9.5 kcal/mol), vardenafil (−9.4 kcal/mol), and tadalafil (−9.3 kcal/mol). Rosmarinic acid and salvinorin A also showed strong binding energies (−9.2 kcal/mol), while catechin followed closely with −9.0 kcal/mol. These results suggest that several plant compounds can bind PDE5 as effectively as, or even better than, established pharmaceutical inhibitors. The accuracy of the docking protocol was confirmed by redocking the co‐crystallized ligand in PDE5 (PDB ID: 2H42), yielding a root‐mean‐square deviation (RMSD) of 0.264 Å (Figure 3). Figure 4 illustrates the surface views and active site interactions of the top‐scoring ligands, while Figure 5 presents LigPlot+ interaction maps showing hydrogen bonding and hydrophobic contacts.

Structural alignment between the redocked pose and the reference (co‐crystallized) ligand within the active site of PDE5. The low RMSD value of 0.264 Å confirms the accuracy of the docking protocol. The original ligand conformation is shown in blue, while the redocked pose is displayed in yellow.

Surface views and active site interactions of PDE5‐ligand complexes. Surface transparency views (left) and active site details (right box) highlighting the enzyme‐ligand interactions. The target enzyme (PDE5) is shown in olive yellow, ligands in blue, hydrogen bonds with their lengths (Å) are represented by dotted yellow lines and interacting amino acid residues linked by hydrogen bonds and hydrophobic interactions are displayed. (A) PDE5–rosmarinic acid complex, (B) PDE5–ellagic acid complex, (C) PDE5–salvinorin A complex, (D) PDE5–catechin complex.

LigPlot+ map of the PDE5‐ligand complexes: Hydrogen bonding (Å) and hydrophobic interactions in the active site. (A) PDE5–rosmarinic acid complex, (B) PDE5–ellagic acid complex, (C) PDE5–salvinorin A complex, (D) PDE5–catechin complex.

LigPlot+ analysis (Table 4, Figure 5) shows the detailed interactions for the highest‐affinity compounds. Ellagic acid formed hydrogen bonds with His617, Gln817, Tyr612, and Asp764, and hydrophobic interactions with Leu765, His613, Phe786, Phe820, Val782, and Ala767. Rosmarinic acid interacted with Gln817, His613, Asp654, Asp764, His617, and Gln775 via hydrogen bonding, and with Ala767, Phe820, Val782, Asn661, Thr723, His657, Glu682, Leu765, and Ile768 via hydrophobic contacts. Catechin formed hydrogen bonds with Gln817, His613, Asp764, and Thr723, alongside hydrophobic interactions with Ala767, Ile768, Val782, Leu765, Leu725, Phe786, and Phe820. Salvinorin A primarily engaged hydrophobic interactions, including Val782, Leu804, Gln817, Met816, Gly819, Ala823, Ile665, Ile824, Leu725, Phe786, and Phe820. These interactions are consistent with prior studies highlighting the importance of residues Gln817 and His613 in PDE5 inhibition (Palanichamy et al. 2022; Li et al. 2010).

Molecular Orbital Analysis

3.3

Molecular orbital properties, including HOMO, LUMO, and energy gap (Eg), were evaluated to assess the reactivity and stability of selected compounds (Table 5). The energy gap reflects chemical reactivity, with smaller values indicating higher reactivity. Ellagic acid exhibited the lowest Eg (0.365 eV), suggesting high reactivity and strong potential for biological interactions. Salvinorin A had the highest Eg (0.451 eV), indicating greater stability but reduced reactivity. Rosmarinic acid and catechin displayed intermediate Eg values (0.374 and 0.441 eV, respectively). HOMO values reflect electron‐donating potential, with rosmarinic acid and catechin showing the highest values, suggesting strong antioxidant and enzyme‐inhibiting capabilities. Conversely, ellagic acid exhibited the lowest LUMO energy, supporting its role as an efficient electron acceptor. Figure 6 illustrates the orbital distributions. These results indicate that ellagic acid may act as a highly reactive inhibitor, while rosmarinic acid and catechin provide a balance between stability and reactivity. Such orbital characteristics reinforce the docking findings and highlight their potential pharmacological applications.

HOMO and LUMO orbital distributions of Rosmarinic acid (A), Ellagic acid (B), Salvinorin A (C), and Catechin (D) with the corresponding energy gaps (Eg) values.

Molecular Dynamics

3.4

MD simulations of PDE5 and its complexes with rosmarinic acid, ellagic acid, salvinorin A, and catechin provide crucial insights into their structural behavior, flexibility, and stability. The RMSD analysis (Figure 7A) revealed that all complexes achieved structural stability after an initial equilibration phase, RMSD values ranging from 0.2 to 0.4 nm. The PDE5‐ellagic acid complex showed slightly higher RMSD values, suggesting increased structural fluctuations compared with the other complexes. This observation indicates that ellagic acid may induce more dynamic conformational changes in the PDE5 structure.

Molecular dynamic analysis results of PDE5 and its complexes with selected ligands: (A) Root Mean Square Deviation (RMSD); (B) Root Mean Square Fluctuation (RMSF); (C) Radius of gyration (Rg) profiles; (D) Solvent Accessible Surface Area.

The RMSF analysis (Figure 7B) provides residue‐specific information on flexibility, with noticeable peaks indicating regions of higher flexibility. Ellagic acid induced greater fluctuations in some regions, indicating more significant movement in specific residues. ThePDE5‐catechin and PDE5‐salvinorin A complexes displayed similar profiles, with overall lower RMSF values, indicating better stability and less fluctuation. This suggests that catechin and salvinorin A promote greater structural rigidity, which could influence their binding stability to PDE5.

The Rg profiles (Figure 7C) further support these findings. The complexes of PDE5 with rosmarinic acid, salvinorin A, and catechin maintained a stable and compact structure (Rg values ranging from 2.0 to 2.2 nm). However, the PDE5‐ellagic acid complex displayed consistently lower Rg values (~1.6–1.7 nm), indicating significant structural compactness (Wu et al. 2022). These differences in Rg values highlight the unique structural behaviors of the ligands. Specifically, ellagic acid induces a more compact configuration and may impact its binding affinity to PDE5.

Concomitantly, the solvent‐accessible surface area (SASA) of the ellagic‐acid complex shrank by ~8% relative to the apo protein and remained the lowest throughout the trajectory (Figure 7D), corroborating the Rg data and suggesting burial of charged moieties that penalize solvation. MM‐GBSA post‐processing (Table 6) rationalizes these observations. Salvinorin A delivers the most favorable binding free energy (ΔGbind = −26.7 kcal mol^−1^) driven by strong van‐der‐Waals contacts (−50.5 kcal mol^−1^). Rosmarinic acid is electrostatically driven (ΔEele −41.7 kcal mol^−1^) and yields ΔGbind = −23.6 kcal mol^−1^. Catechin gives marginal affinity (ΔGbind = −13.1 kcal mol^−1^). Ellagic acid presents attractive vacuum terms (ΔEvdW + ΔEele = −30.4 kcal mol^−1^) but the largest desolvation penalty (+17.7 kcal mol^−1^), flipping the net balance to an unfavorable +52.1 kcal mol^−1^. Thus, is fully consistent with its elevated RMSD/RMSF and reduced SASA.

Hydrogen‐bonding analysis (Figure 8) further clarifies the interaction patterns that underlie the stability differences between complexes. Time traces of the number of hydrogen bonds formed between PDE5 and each ligand across the 100‐ns trajectories. Rosmarinic acid (Figure 8A) maintained a relatively high and persistent hydrogen‐bond network throughout the simulation, with an estimated average of ~3.2 ± 1.0 H‐bonds and repeated transient increases to 4–6 bonds. Ellagic acid (Figure 8B) exhibited larger temporal fluctuations and a rising trend in the final 20 ns (estimated average ~3.8 ± 1.4 H‐bonds), consistent with its higher RMSD and RMSF values and suggesting ellagic acid both forms multiple H‐bonds and induces dynamic conformational rearrangements of the binding site. In contrast, salvinorin A (Figure 8C) and catechin (Figure 8D) formed fewer hydrogen bonds (estimated averages ~0.9 ± 0.4 and ~1.0 ± 0.5 H‐bonds, respectively) that were mostly intermittent (predominantly 0–1 H‐bond with occasional 2‐bond spikes). The lower number of H‐bonds for salvinorin A and catechin indicates these complexes rely more on non‐polar contacts and steric complementarity to maintain stability, which matches their lower RMSF and compact Rg profiles. Overall, the hydrogen‐bonding patterns reinforce the previous observations: ellagic acid induces greater structural fluctuation (more and more variable H‐bonding), rosmarinic acid provides a consistent hydrogen‐bonding network associated with maintained stability, while salvinorin A and catechin form only transient H‐bonds but still stabilize PDE5 through other interactions.

Time evolution of the number of hydrogen bonds between PDE5 and (A) Rosmarinic acid, (B) Ellagic acid, (C) Salvinorin A, and (D) Catechin during 100 ns MD simulations. Hydrogen bonds were monitored along the trajectories and plotted as the instantaneous count versus simulation time (ns). Panels show that rosmarinic acid forms a persistent network of H‐bonds, ellagic acid exhibits larger fluctuations and an increase in H‐bonding toward the end of the run, whereas salvinorin A and catechin form fewer and more intermittent H‐bonds.

In addition, NMA offered valuable insights into PDE5 flexibility and functional dynamics. The deformability plot (Figure 9A) revealed regions of varying flexibility within PDE5, with deformability values ranging from 1.33E‐03 (low flexibility) to 8.20 × 10^−3^ (high flexibility). Notably, the indices 1 to 6 exhibited higher deformability, indicating flexible regions that are crucial for functional movements. Conversely, the indices 12–20 showed lower deformability, suggesting greater stability in regions vital for PDE5 structural integrity. Regions between indices 60 and 70, which demonstrate fluctuations, may play a significant role in PDE5‐ligand interactions.

Molecular dynamic simulations using the normal mode analysis of the PDE5 protein: (A) Deformability plots, (B) B‐factor plots, (C) Eigenvalues, (D) Variance map, (E) Correlation Matrix (red regions indicate correlated motions, blue regions represent anti‐correlated motions, and white regions show uncorrelated movement), and (F) Elastic Network Model (regions correspond to stiffer interactions and lighter regions to more flexible regions).

The B‐factor diagram (Figure 9B) provides additional insights into the thermal vibrations of the PDE5‐ligand complexes. Lower B‐factors correspond to more rigid regions, while higher B‐factors indicate flexibility. The strong correlation (Pearson coefficient = 0.85) between experimental and calculated B‐factors further confirmed the model reliability in capturing the dynamic behavior of the complexes. Residues 45–60, which exhibit high B‐factor values (> 40 Å^2^), correlated with PDE5 regions involved in conformational changes, highlighting their dynamic nature.

The Eigenvalue analysis (Figure 9C) revealed the vibration modes of the complexes with calculated eigenvalue of 3.597134 × 10^−5^ indicated significant flexibility (Sumera et al. 2022). Lower eigenvalues reflect easier deformation, which aligns with regions exhibiting mobility and suggests that minimal energy is required for conformational changes.

The variance map (Figure 9D) further emphasized regions of atomic fluctuation, with the cumulative variance reaching 100% by mode index 14. Low‐frequency modes, particularly modes 2–12, contribute significantly to atomic motion, reflecting greater structural flexibility in these regions.

The covariance matrix (Figure 9E) provided a detailed view of residue‐residue dynamic coupling. Red regions indicate correlated motions, blue regions represent anti‐correlated motions, and white regions show uncorrelated movements. This matrix revealed the intricate dynamics of the PDE5‐ligand complexes, emphasizing areas where residues move coordinated, potentially influencing ligand binding and protein function. The elastic network plot (Figure 10F) visualized the atom‐to‐atom connections within the complex, where darker regions correspond to stiffer interactions and lighter regions to more flexible regions. In summary, the MD simulations, supported by the NMA results, offer an in‐depth analysis of structural stability, flexibility, and interactions within the PDE5‐ligand complexes. The results emphasize the dynamic interplay of rigidity and flexibility, which is critical for understanding the ligand binding affinity and the functional adaptability of PDE5. The detailed MD findings highlight how each ligand affects PDE5 structural integrity and flexibility. Ellagic acid induced more fluctuations, while rosmarinic acid, salvinorin A, and catechin promoted greater stability and moderate flexibility, making them versatile candidates for therapeutic applications.

QSAR model validation for PDE5 inhibitors using Random Forest regression. (A) Predicted vs. experimental pKd values. (B) Residual distribution. (C) Williams plot (applicability domain). (D) Y‐scrambling validation. (E) Learning curve.

QSAR And Structure–Activity Relationship Analysis

3.5

To complement docking and MD simulations, a QSAR model was developed to quantitatively assess the structural determinants of PDE5 inhibition. Docking‐derived binding energies (ΔG) were first converted into predicted pKd values and then modeled against a set of molecular descriptors using a Random Forest (RF) regression approach. The RF model demonstrated strong predictive performance, with high correlation between experimental and predicted pKd values (R ^2^_train = 0.91, R ^2^_CV = 0.87) (Figure 10A). The residual distribution plot confirmed that prediction errors were symmetrically centered around zero (Figure 10B), indicating unbiased estimates. The Williams plot showed that the majority of compounds fell within the model's applicability domain (Figure 10C). Furthermore, Y‐scrambling validation produced markedly reduced R ^2^ and Q ^2^ values, ruling out chance correlations and confirming model robustness (Figure 10D). The learning curve revealed convergence between training and cross‐validation performance, supporting both dataset sufficiency and generalization capability (Figure 9E). QSAR analysis highlighted the key molecular features underlying PDE5 inhibition. Catechin exhibited strong activity due to its multiple hydrogen bond donors (HBD = 4), moderate polar surface area (TPSA = 110 Å^2^), and favorable lipophilicity (logP = 1.23). Ellagic acid benefited from its conjugated aromatic scaffold and high hydrogen donor count, while rosmarinic acid displayed an optimal balance of molecular weight (360 g/mol), TPSA (125 Å^2^), and multiple donors (HBD = 4). In contrast, Salvinorin A, despite lacking hydrogen bond donors, compensated with its hydrophobic aromatic skeleton (logP≈2.8), achieving notable affinity. Polyphenolic flavonoids such as apigenin, quercetin, kaempferol, and naringenin also displayed strong binding affinities (−8.5 to −8.9 kcal/mol), attributable to shared hydroxylation patterns and moderate TPSA values, consistent with previous QSAR studies on flavonoid bioactivity. Overall, these findings demonstrate that potent PDE5 inhibition arises from a fine balance between hydrogen‐bonding capacity, moderate polarity, and hydrophobic aromatic scaffolds. Flavonoids and polyphenols therefore represent promising structural frameworks for the design and development of natural PDE5 inhibitors.

Conclusion

4

This study provides a comprehensive in silico evaluation of 76 natural compounds derived from Algerian medicinal plants as potential PDE5 inhibitors for the treatment of erectile dysfunction (ED). Drug‐likeness and toxicity predictions revealed that the majority of compounds complied with pharmacokinetic rules and exhibited low toxicity, underscoring their suitability as safe therapeutic candidates. Molecular docking identified several compounds, including ellagic acid, rosmarinic acid, catechin, and salvinorin A, with binding affinities comparable to sildenafil, the reference PDE5 inhibitor. MD simulations further validated these interactions, demonstrating that catechin, rosmarinic acid, and salvinorin A formed stable complexes with PDE5, maintaining compact protein structures and stable binding throughout the 100 ns trajectories. In contrast, ellagic acid displayed higher fluctuations, suggesting reduced conformational stability. Frontier molecular orbital analysis confirmed the electronic reactivity of catechin and rosmarinic acid, supporting their strong interaction potential with PDE5. Importantly, the integrated QSAR analysis provided quantitative confirmation of these findings, revealing that effective PDE5 inhibition is driven by an optimal balance of hydrogen‐bond donor capacity, moderate polarity, and hydrophobic aromatic scaffolds. Flavonoids and polyphenols (i.e., catechin, quercetin, apigenin, and rosmarinic acid) emerged as the most promising scaffolds, with predictive modeling reinforcing their inhibitory potential and structural suitability. Taken together, these results highlight catechin, rosmarinic acid, and related flavonoids as promising natural PDE5 inhibitors with both stability and safety advantages over synthetic alternatives. This integrated workflow, combining drug‐likeness filtering, docking, MD simulations, and QSAR modeling, provides a powerful strategy for identifying natural leads in ED management and lays the foundation for future preclinical and clinical validation.

Author Contributions

Farouk Boudou designed the study, supervised the project, managed administration, acquired funding, and wrote and critically reviewed the manuscript. Alaeddine Berkane contributed equally by developing software, performing data curation, conducting molecular docking and virtual screening, and validating results. Amal Belakredar contributed equally to conceptualization, software development, validation, and visualization, and drafted sections of the manuscript. Ahcene Keziz contributed equally by developing methodology, performing computational modeling, and analyzing results. Huda Alsaeedi selected phytochemicals, managed databases, conducted literature review, and reviewed and edited the manuscript. Brian A. Murray evaluated the study critically, provided scientific guidance, and edited the manuscript. Mikhael Bechelany supported methodology, provided computational resources, and validated results. Ahmed Barhoum supervised the project, coordinated administration, acquired funding, and edited the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.