CC-BY 4.0
1. Introduction
Pulmonary arterial hypertension (PAH) is a serious and progressive disorder marked by increased pressure in the pulmonary arteries, potentially resulting in right ventricular hypertrophy and ultimately heart failure [1, 2, 3]. The condition may be categorized as idiopathic, lacking an identifiable underlying cause, or associated with other diseases, including chronic obstructive pulmonary disease (COPD), congenital heart defects, and connective tissue disorders [4, 5]. The precise etiology of PAH is intricate, frequently encompassing a blend of genetic predispositions and environmental influences that lead to endothelial dysfunction, inflammation, and vascular remodeling [6, 7, 8].
The pathophysiological characteristics of PAH involve significant alterations in the pulmonary vasculature. Endothelial dysfunction is a critical characteristic that disrupts the equilibrium between vasodilatory and vasoconstrictive factors, resulting in the excessive proliferation of vascular smooth muscle cells [1, 2, 3]. Inflammation exacerbates this process, playing a central role in disease progression, alongside a reduction in programmed cell death (apoptosis) and the subsequent vascular remodeling that characterizes advanced disease stages. These processes result in the constriction of the pulmonary arteries, heightened pulmonary vascular resistance, and increased pressures in the pulmonary circulation. The right ventricle exerts progressively greater effort to circulate blood through the lungs, resulting in hypertrophy and potential heart failure if not addressed [9, 10].
Recent advances in the management of PAH have led to the development of various therapeutic strategies designed to alleviate symptoms, enhance functional capacity, and decelerate disease progression. Phosphodiesterase type 5 (PDE5) inhibitors represent a significant advancement in the therapy for PAH. Sildenafil (SIL), tadalafil (TAD) and avanafil (AVA) (Figure 1) are pharmacological agents acting on the nitric oxide (NO)–cyclic guanosine monophosphate (cGMP) pathway, a crucial vascular tone regulator [11, 12, 13, 14]. These drugs inhibit PDE5, an enzyme responsible for the degradation of cGMP, leading to elevated intracellular cGMP levels. This elevation results in the relaxation of pulmonary artery smooth muscle cells and a decrease in pulmonary vascular resistance [15]. First PDE5 inhibitor approved for PAH treatment, SIL has demonstrated significant clinical efficacy, though its relatively short half-life necessitates frequent dosing to maintain therapeutic effects [11]. Another PDE5 inhibitor, TAD, has a longer half-life, which allows less frequent dosing and makes it easier for patients to follow their treatment regimen [12]. More recently, AVA has been introduced as a highly selective PDE5 inhibitor, offering greater selectivity and a more favorable side-effect profile compared to earlier drugs, which reduces the risk of cardiovascular complications [13].
Chemical structures of the investigated compounds.
In addition to the synthetic PDE5 inhibitors, there has been growing interest in the potential of naturally occurring compounds as alternative or complementary therapies for PAH. Usnic acid (UA, Figure 1), a lichen-derived dibenzofuran, has emerged as a compound of interest due to its diverse biological activities [16, 17, 18], including antioxidant [19, 20], anti-inflammatory [21], and antimicrobial properties [22]. Although UA has not been directly investigated for its inhibitory activity on PDE5, its mechanisms of action in other biological contexts—such as reducing oxidative stress and inflammation—suggest that it warrants further exploration as a potential treatment option for PAH. Given that oxidative stress and inflammation are critical processes implicated in the pathogenesis of PAH, UA may still offer therapeutic potential through these pathways, possibly complementing existing treatments targeting the NO–cGMP axis [15]. This investigation into UA activity adds to the growing body of research that seeks to identify natural compounds with pharmacological effects that may rival or complement synthetic drugs. Its role as a potential PDE5 inhibitor, or as a compound influencing alternative pathways involved in pulmonary vasodilation, makes it a subject of considerable scientific interest.
This study aims to evaluate the inhibitory potential of usnic acid (UA) in comparison to clinically used PDE5 inhibitors AVA, SIL and TAD through a comprehensive in silico approach. ADMET profiling was performed to assess UA’s pharmacokinetic and toxicity-related properties relative to reference compounds. Molecular docking and molecular dynamic simulations explored the binding affinity, interaction patterns, and stability of enzyme–ligand complexes within the PDE5 active site. This integrative analysis seeks to elucidate the molecular basis of UA’s activity and its potential as a novel therapeutic candidate for the treatment of PAH.
2. Materials and methods
2.1. ADMET analysis
The ADMETlab 3.0 web server [23] facilitated a thorough assessment of the toxicity profiles and pharmacokinetic properties of compounds. A comprehensive ADMET analysis (Absorption, Distribution, Metabolism, Excretion, and Toxicity) is crucial for assessing the safety and biological efficacy of compounds. The analysis presented essential data regarding absorption, tissue distribution, metabolic stability, elimination pathways, and potential toxicities, facilitating predictions about behavior within a biological system.
2.2. DFT method
Computational analyses were performed using the Gaussian16 software suite [24], while molecular visualizations were executed with GaussView 6.0.16 [25]. Quantum chemical calculations were conducted to optimize molecular geometries using density functional theory (DFT) with the B3LYP-D3BJ functional and the 6-311++G(d,p) basis set [26, 27]. This combination utilizes polarization and diffuse functions to improve computational accuracy. The selected theoretical model was based on its proven capacity to closely match experimental geometrical parameters, as evidenced by multiple studies [28, 29, 30].
2.3. Molecular docking simulation
The computational analysis of PDE5 inhibition involved molecular docking simulations to evaluate the binding affinities and interactions of SIL, TAD, AVA, and UA. The docking studies utilized AutoDock 4.2 software [31], employing the Lamarckian genetic algorithm [32] to optimize binding poses and assess interaction energies. In this study, parameters were defined to ensure precise and dependable predictions: a maximum of 250 000 energy evaluations were established, and 27 000 generations were permitted for the evolutionary search. The mutation rate was established at 0.02, while the crossover rate was determined to be 0.8. These parameters maintain a balance between computational efficiency and comprehensive exploration of potential ligand conformations. Kollmann partial charges were assigned to the ligands, and polar hydrogens were added using AutoDockTools to enhance the accuracy of electrostatic interaction predictions in the docking simulations. In the docking procedure, ligands, whose structures were optimized according to the previously described method, were treated as flexible entities, enabling them to assume multiple conformations within the PDE5 active site, whereas the protein structure was maintained in a rigid state. This method facilitates the investigation of ligand flexibility, crucial for comprehending the dynamic characteristics of ligand–protein interactions. The protein structure of PDE5 was sourced from the Protein Data Bank (PDB code: 6L6E, accessed on 18.12.2023) [33]. Preprocessing was conducted using Discovery Studio 4.0, during which co-crystallized water molecules and ligands were eliminated to facilitate the preparation of the binding site. The docking grid box was centered at coordinates −2.479 × −26.083 × −27.632, featuring a grid size of 53 × 69 × 56 points and a grid spacing of 0.375 Å. The selected settings accurately define the active site of PDE5, ensuring comprehensive coverage of the region involved in ligand binding.
2.4. Molecular dynamic simulation
All molecular dynamics simulations were conducted utilizing the AMBER22 software suite [34], while AmberTools23 was used for system preparation [35]. Protein–ligand complexes were constructed utilizing the tleap module, employing the ff19SB force field for the protein and GAFF2 (General Amber Force Field 2) for the ligands [36]. The initial protein–ligand complexes were obtained from previous molecular docking simulations. Each complex was solubilized in a truncated octahedral box of TIP3P water molecules, maintaining a 12 Å buffer surrounding the solute. Counterions were added to neutralize the charge of the system. Energy minimization occurred in two phases: the initial phase involved the application of 500 kcal/mol/Å2 restraints to the solute, while the solvent and ions were minimized for 5000 steps to resolve potential steric clashes. During the second phase, all restraints were lifted, and the complete system underwent an additional 10 000 steps of full minimization to achieve structural relaxation. After minimization, the system was incrementally heated from 0 to 300 K in 100 ps, utilizing a Langevin thermostat with a collision frequency of 2 ps−1. In this phase, solute atoms were constrained using a force constant of 10 kcal/mol/Å2. A 1 ns equilibration run was conducted at a constant pressure of 1 atm and a temperature of 300 K following heating, during which the restraints were gradually diminished. The molecular dynamics simulations, lasting 100 ns, were performed utilizing the CUDA-enabled pmemd engine in AMBER22, facilitating efficient GPU acceleration [34, 35, 36]. These simulations employed the NPT ensemble, ensuring a constant number of particles, pressure, and temperature, while applying periodic boundary conditions. The SHAKE algorithm [37] constrained hydrogen-involved bonds, facilitating the application of a 2 fs time step. Long-range electrostatic interactions were handled using the particle mesh Ewald (PME) method [38], with a 10 Å cutoff applied for nonbonded interactions.
2.4.1. Molecular mechanics generalized Born surface area
The binding energies of the protein–ligand complexes were estimated using the molecular mechanics generalized Born surface area (MM-GBSA), as implemented in the MMPBSA.py module of AMBER [35, 39]. A total of 1000 snapshots, derived from the final 20 ns of the 100 ns molecular dynamics simulations, were utilized to confirm system equilibrium. The binding free energy (ΔGbind) is determined by the difference between the free energy of the complex and the sum of free energies of the receptor and ligand in their unbound states (Equation (1)):
| \begin {equation}\label {eq1} \Delta G_{\mathrm {bind}} = G_{\mathrm {complex}} - (G_{\mathrm {receptor}} + G_{\mathrm {ligand}}) \end {equation} | (1) |
The free energy of each system was partitioned into gas-phase and solvation components. The gas-phase energy (ΔGgas) was calculated by summing the van der Waals and electrostatic interactions between protein and ligand:
| \begin {equation}\label {eq2} \Delta G_{\mathrm {gas}} = E_{\mathrm {VDW}} + E_{\mathrm {ele}} \end {equation} | (2) |
The solvation free energy (Gsolv) was divided into polar and nonpolar components. In MM-GBSA, the polar solvation energy (Gpolar) is determined using the generalized Born model, which approximates the electrostatic component of solvation. The nonpolar solvation energy (Gnonpolar) is estimated based on the solvent-accessible surface area (SASA) utilizing the linear combination of pairwise overlap (LCPO) method. The nonpolar solvation energy was calculated as follows:
| \begin {equation}\label {eq3} G_{\mathrm {nonpolar}} = \gamma \times \mathrm {SASA} + b \end {equation} | (3) |
The total solvation free energy (ΔGsolv) was then determined as the sum of the polar and nonpolar components:
| \begin {equation}\label {eq4} \Delta G_{\mathrm {solv}} = G_{\mathrm {polar}} + G_{\mathrm {nonpolar}}. \end {equation} | (4) |
| \begin {equation}\label {eq5} \Delta G_{\mathrm {bind}} = \Delta G_{\mathrm {gas}} + \Delta G_{\mathrm {solv}}. \end {equation} | (5) |
The decomposing of binding energy into van der Waals, electrostatic, and solvation components facilitated a thorough examination of the factors influencing ligand binding to the PDE5 enzyme. The evaluation of these contributions facilitated a comprehensive understanding of the factors affecting the stability and affinity of the ligand–protein complex.
3. Results and discussion
3.1. ADMET analysis of investigated compounds
Standing for “Absorption, Distribution, Metabolism, Elimination, Toxicity”, ADMET analysis is essential in drug discovery, evaluating compounds for their pharmacokinetic and pharmacodynamic suitability. The bioavailability radar plots in Figure 2 depict 13 physicochemical parameters for SIL, TAD, AVA, and UA, with the yellow zone representing optimal bioavailability ranges. Comparatively, UA has a lower molecular weight than SIL, TAD, and AVA, which can enhance membrane permeability and systemic absorption. Its moderate lipophilicity (logP) balances solubility and distribution, reducing the risk of tissue accumulation. The higher topological polar surface area (TPSA) of UA suggests limited central nervous system (CNS) penetration, potentially lowering CNS side effects compared to PDE5 inhibitors.
Bioavailability radar chart of (a) SIL, (b) TAD, (c) AVA and (d) UA with the following parameters: molecular weight (MW), number of rigid bonds (nRig), formal charge (fChar), number of heteroatoms (nHet), number of atoms in the biggest ring (MaxRing), number of rings (nRing), number of rotable bonds (nRot), topological polar surface area (TPSA), number of hydrogen bond donors (nHD), number of hydrogen bond acceptors (nHA), log of the octanol/water partition coefficient (logP), logP at physiological pH (logD), and log of aqueous solubility (logS).
Regarding hydrogen bonding, UA shows increased H-bond donors and acceptors, which may enhance water solubility and gastrointestinal absorption. The rigidity of UA’s structure, with fewer rotatable bonds, may improve stability in binding but limits conformational flexibility. Finally, UA’s favorable water solubility (logS) supports effective oral absorption, while structural heteroatoms and a mild formal charge promote stable binding in diverse environments. Overall, UA’s unique pharmacokinetic profile—characterized by favorable solubility, balanced lipophilicity, and reduced CNS impact—presents a promising therapeutic profile compared to synthetic PDE5 inhibitors.
Absorption of the potential drug implies the transition into the organism’s systemic circulation. The evaluated parameters encompass Caco-2 permeability (a model for intestinal absorption using human colon carcinoma cells) [40, 41], MDCK permeability (Madin–Darby canine kidney cells used to assess renal absorption) [42], human intestinal absorption (HIA) [43], PAMPA permeability (parallel artificial membrane permeability assay, a model for passive diffusion across lipid membranes) [44], the biopharmaceutical absorption model (F20%, representing the fraction of drug absorbed) [45], plasma protein binding (PPB) [46], the fraction of unbound drugs (Fu%, indicating the portion of the drug that remains unbound in the plasma) [47], volume of distribution (VD, representing the extent of drug distribution into tissues) [48], and interactions with OATP1B1 and OATP1B3 transporters (organic anion transporting polypeptides, important for drug transport into liver cells) [49] (Table 1).
Pharmacokinetic properties of investigated compounds based on ADMET analysis: absorption and distribution parameters of investigated compounds
| Comp. | Absorption | Distribution | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Caco-2p (cm⋅s−1) | MDCK (cm⋅s−1) | HIA | PAMPA | F20% | PPB% | Fu% | VD (L⋅kg−1) | OATP1B1/OATP1B3 | |
| SIL | −4.812 | −4.771 | 0.003 | 0.171 | 0.049 | 96.845 | 2.745 | 0.153 | 0.017/0.005 |
| TAD | −5.138 | −4.401 | 0.000 | 0.007 | 1.000 | 96.596 | 2.968 | −0.084 | 0.092/0.000 |
| AVA | −5.453 | −5.054 | 0.001 | 0.787 | 0.013 | 96.077 | 3.136 | −0.130 | 0.998/0.998 |
| UA | −5.197 | −4.787 | 0.003 | 0.899 | 0.008 | 97.162 | 3.164 | −0.465 | 0.724/0.499 |
The absorption parameters indicate distinct profiles for the compounds under investigation. The Caco-2 and MDCK models indicate low permeability for UA (−5.197 and −4.787 cm/s), aligning with the values for TAD and AVA, whereas SIL demonstrates marginally improved permeability. The PAMPA model demonstrates that UA exhibits significantly higher permeability (0.899) than SIL, TAD, and AVA, suggesting enhanced passive diffusion across lipid membranes for UA. The HIA parameter for UA (0.003) aligns closely with SIL and AVA, indicating restricted absorption via the human intestinal system. The absorption fraction, denoted as F20%, shows nearly complete absorption for TAD (1.000), whereas UA exhibits a markedly lower value (0.008). The results indicate that, while UA demonstrates significant permeability in the PAMPA model, its overall absorption is restricted relative to other compounds, especially in comparison to TAD.
Plasma protein binding (PPB) indicates that UA exhibits a high binding rate of 97.162%, consistent with the values observed for other drugs. Although this may decrease the availability of the unbound drug in plasma, high binding may also extend the half-life in circulation. The unbound drug fraction (Fu%) for UA is 3.164%, comparable to the values for AVA. In contrast, SIL and TAD exhibit slightly lower values, suggesting that UA may persist in its unbound form in plasma for a longer duration than SIL. The volume of distribution (VD) for UA (3.164 L/kg) exceeds that of SIL (2.745 L/kg) while being comparable to AVA, suggesting considerable tissue distribution (Table 1). This parameter indicates that UA may be widely distributed in tissues, which could have implications for its therapeutic efficacy and potential toxicity.
The interaction between OATP1B1 and OATP1B3 transporters is essential for the hepatic distribution of drugs. The significant interaction of UA with OATP1B1/B3 (0.724/0.499) markedly exceeds the interactions observed with SIL and TAD. This may suggest enhanced uptake of UA by hepatocytes, potentially influencing its distribution and elimination. Minimal interaction occurs between SIL and OATP transporters, indicating a reduced level of hepatic uptake.
These results indicate that UA exhibits distribution characteristics akin to those of AVA, yet it possesses a unique absorption profile. Although UA exhibits significant permeability across lipid membranes (PAMPA), its overall absorption remains constrained, potentially impacting its bioavailability. The interaction with OATP1B1/B3 transporters indicates a potential for enhanced hepatic uptake, which could affect its metabolism and elimination.
Table 2 presents the estimated parameters of metabolism and elimination for the investigated compounds in the ADMET analysis, providing key insights. The metabolism of the compounds is assessed through interactions with cytochrome P450 enzymes (CYP1A2, CYP2D6, and CYP2C9), with values indicating their potential as inhibitors (I) or substrates (S) of these enzymes [50]. Elimination is evaluated by measuring the half-life (T1/2), which indicates how long the compound remains active in the body [51], and clearance (CL), representing the rate at which the compound is removed from the bloodstream [52].
Pharmacokinetic properties of investigated compounds based on ADMET analysis: metabolism and elimination parameters of investigated compounds
| Compounds | Metabolism; P450 | Elimination | |||||||
|---|---|---|---|---|---|---|---|---|---|
| 1A2 | 2D6 | 2C9 | HLM | T1/2 (h) | CL (mL⋅min−1⋅kg−1) | ||||
| I | S | I | S | I | S | ||||
| SIL | 0.000 | 0.167 | 0.000 | 0.000 | 0.003 | 0.999 | 0.999 | 0.400 | 8.757 |
| TAD | 0.000 | 1.000 | 0.000 | 0.000 | 0.948 | 0.000 | 0.008 | 1.023 | 2.930 |
| AVA | 0.966 | 0.374 | 0.045 | 0.200 | 0.791 | 0.984 | 0.641 | 0.515 | 7.078 |
| UA | 0.479 | 0.998 | 0.000 | 0.000 | 0.006 | 0.001 | 0.258 | 1.446 | 1.413 |
Regarding interaction with the 1A2 enzyme, AVA exhibits the highest level of inhibition (I = 0.966), whereas UA displays moderate inhibition (I = 0.479). SIL and TAD exhibit no inhibitory effects on this enzyme, whereas TAD serves as the primary substrate (S = 1.000) for 1A2. There is considerable potential for UA as a substrate for this enzyme (S = 0.998), highlighting the significant involvement of 1A2 in its metabolism. In the analysis of interactions with the 2D6 enzyme, SIL and UA exhibit low or negligible inhibition. They do not serve as substrates for this enzyme, whereas AVA exhibits low inhibition (I = 0.045) and moderate substrate activity (S = 0.200). The 2C9 enzyme should be significantly inhibited by AVA (I = 0.791), whereas UA shows minimal inhibition (I = 0.006). The values indicate distinct metabolic pathways for these compounds within the body. The elimination parameters indicate that UA possesses the longest half-life (T1/2 = 1.446 h), which is markedly longer than that of SIL (T1/2 = 0.400 h) and AVA (T1/2 = 0.515 h). This suggests that UA has an extended duration in the body, which may influence its therapeutic profile and cumulative effects. The clearance (CL) for UA is the lowest at 1.413 mL/min/kg, indicating a slower elimination rate compared to other compounds studied, particularly SIL at 8.757 mL/min/kg and AVA at 7.078 mL/min/kg. TAD exhibits a relatively low clearance rate of 2.930 mL/min/kg, although it remains faster than UA. The findings indicate that UA demonstrates distinct metabolic and elimination properties, exhibiting significant inhibition and substrate activities for the 1A2 enzyme, while displaying limited interaction with 2D6 and 2C9. The longer half-life and slower clearance of UA relative to other compounds suggests its potential for prolonged presence in the body, which may be significant for its pharmacokinetic and pharmacodynamic characterization.
Table 3 provides estimated toxicity parameters from the ADMET analysis of the investigated compounds, offering key insights into potential toxic effects. The analysis included parameters such as hERG channel inhibition (assessing the risk of cardiotoxicity by evaluating the potential to block the human ether-à-go-go-related gene potassium channel) [53, 54], drug-induced liver injury (DILI) [55], human hepatotoxicity (the potential to cause liver damage) [56], the AMES test (evaluating mutagenicity or the potential to cause genetic mutations) [57, 58], hematotoxicity (toxicity to blood cells), carcinogenicity (potential to cause cancer), nephrotoxicity (toxicity to the kidneys), and neurotoxicity (toxicity to the nervous system) for SIL, TAD, AVA, and UA.
Toxicological parameters of investigated compounds based on ADMET analysis
| Compounds | hERG | DILI | Human hepatotoxicity | AMES | Hematotoxicity | Carcinogenicity | Nephrotoxicity | Neurotoxicity |
|---|---|---|---|---|---|---|---|---|
| SIL | 0.570 | 1.000 | 0.966 | 0.507 | 0.720 | 0.894 | 0.968 | 0.935 |
| TAD | 0.372 | 0.995 | 0.984 | 0.125 | 0.776 | 0.231 | 0.988 | 0.979 |
| AVA | 0.821 | 0.985 | 0.979 | 0.562 | 0.183 | 0.691 | 0.946 | 0.970 |
| UA | 0.003 | 0.739 | 0.339 | 0.252 | 0.675 | 0.340 | 0.071 | 0.099 |
The lowest hERG channel inhibition value is observed for UA at 0.003, indicating a reduced risk of cardiotoxicity in comparison to SIL (0.570) and AVA (0.821), which exhibit a greater risk of channel inhibition and are potentially associated with cardiovascular complications. About the DILI parameter, UA exhibits a lower value of 0.739 in contrast to SIL at 1.000 and TAD at 0.995, suggesting a reduced risk of drug-induced liver injury. The human hepatotoxicity value for UA is 0.339, categorizing it as medium risk (yellow category). In contrast, SIL (0.966), TAD (0.984), and AVA (0.979) are classified as high risk for hepatotoxicity (red category).
The Ames test for mutagenicity indicates a low value for UA (0.252), comparable to TAD (0.125), while AVA and SIL exhibit higher risks than UA (0.562 and 0.507, respectively). The hematotoxicity of UA (0.675) is comparable to that of SIL, whereas AVA demonstrates markedly lower hematotoxicity (0.183). The carcinogenicity of UA is 0.340, which is lower than that of SIL at 0.894 and AVA at 0.691, but higher than TAD, which has a value of 0.231. Nephrotoxicity for UA (0.071) is significantly lower than that of the other compounds, indicating a diminished risk of kidney damage. The neurotoxicity of UA (0.099) is significantly lower than that of SIL (0.935), TAD (0.979), and AVA (0.970), suggesting a reduced risk of neurological complications. The data indicate that UA presents a low risk of cardiotoxicity and neurotoxicity, a medium risk of hepatotoxicity, and a reduced risk of nephrotoxicity, which may be pertinent to its long-term use.
3.2. Inhibitory activity toward PDE5 enzyme—molecular docking study
Before analyzing the docking results, the optimized geometries of the compounds—SIL, TAD, AVA, and UA—were determined using DFT methods (Figure 3). The rigid structure of UA is stabilized by intramolecular hydrogen bonds and van der Waals interactions, with key bond distances (1.605 and 1.709 Å) contributing to its stability. In contrast, PDE5 inhibitors such as SIL and TAD exhibit flexibility in their aromatic rings, enhancing their adaptability. AVA, with hydrogen bond lengths of 1.838 and 2.614 Å, is further stabilized through these key hydrogen bonding interactions, contributing to the overall stability of the molecule.
Optimized geometries of investigated compounds at the B3LYP-D3BJ/6-311++G(d,p) level of theory. Legend: grey—carbon atoms, red—oxygen atoms, blue—nitrogen atoms, white—hydrogen atoms, yellow—sulfur atom (in SIL), green—chlorine atom (in AVA). The dashed line represents a hydrogen bond and the number indicates the bond length (Å).
The molecular docking study detailed in Table 4 elucidates the thermodynamic parameters associated with the binding of investigated compounds to PDE5. While not a traditional PDE5 inhibitor, UA is highlighted as a natural compound that demonstrates considerable potential relative to established synthetic inhibitors. The binding free energy (ΔGbind) of UA at −10.12 kcal/mol is marginally higher than that of AVA (−11.39 kcal/mol) and SIL (−11.01 kcal/mol) yet still signifies a stable interaction with PDE5. The Ki value for UA is 3.8 × 10−2 μM, which exceeds that of AVA (4.5 × 10−3 μM) and SIL (8.5 × 10−3 μM), suggesting that a greater concentration of UA is necessary to attain half-maximal enzyme inhibition. The natural origin of UA enhances its importance, as natural compounds typically exhibit favorable safety profiles and possess additional biological effects, including antioxidant and anti-inflammatory properties. The energetic components of UA binding, including ΔGvdw+hbond+desolv (−9.69 kcal/mol) and ΔGelec (−0.10 kcal/mol), indicate stable intermolecular interactions. However, its ΔGtotal energy (−2.03 kcal/mol) is relatively lower than that of the other inhibitors. The findings indicate the potential of UA in the therapy of PAH, functioning as a PDE5 inhibitor and possessing natural properties that may aid in reducing oxidative stress and vascular remodeling, which are critical factors in the disease’s pathogenesis.
Important thermodynamic parameters (ΔGbind: binding free energy, Ki: inhibition constant, ΔGtotal: total internal energy, ΔGtor: torsional free energy, ΔGunb: energy of the unbound system, ΔGelec: electrostatic energy, and ΔGvdw+hbond+desolv: sum of van der Waals interactions, hydrogen bonds, and desolvation energies, in kcal/mol) for the most stable conformations of the investigated compounds within the PDE5 enzyme active site, obtained through molecular docking simulations
| Compound | ΔGbind | Ki (μM) | ΔGinter | ΔGvdw+hbond+desolv | ΔGelec | ΔGtotal | ΔGtor | ΔGunb |
|---|---|---|---|---|---|---|---|---|
| PDE5–SIL | −11.01 | 8.5 × 10−3 | −12.36 | −11.94 | −0.42 | −1.92 | 2.09 | −1.18 |
| PDE5–TAD | −10.36 | 2.5 × 10−2 | −10.83 | −10.77 | −0.06 | −0.61 | 0.30 | −0.78 |
| PDE5–AVA | −11.39 | 4.5 × 10−3 | −13.34 | −12.19 | −1.15 | −1.86 | 2.68 | −1.13 |
| PDE5–UA | −10.12 | 3.8 × 10−2 | −9.79 | −9.69 | −0.10 | −2.03 | 1.49 | −0.21 |
Figure 4 presents the crystal structure of the PDE5 enzyme in a complex with AVA (PDB: 6L6E), which features a well-defined active site [31]. This active site served as the basis for subsequent studies, where other investigated compounds, including SIL, TAD, and UA, were also examined. The structure of the active site allowed for a detailed analysis of the binding modes of various inhibitors, providing insights into the key factors influencing their selectivity and inhibitory potential.
Crystal structure of PDE5 enzyme in complex with AVA (PDB: 6L6E). The red α-helices represent the helical regions of the PDE5 enzyme, contributing to the overall stability and structural framework of the protein. The green loops and β-turns depict flexible regions of the protein, which allow conformational changes critical for ligand binding. The ligand AVA (avanafil) is located in the active site of PDE5 and is displayed near the center of the image. Carbon atoms are grey, nitrogen atoms blue, oxygen atoms red, hydrogen atoms white, and chlorine atoms in green. Water molecules (shown as small red dots) play a role in the stabilization of the ligand and protein structure.
Comparing the intermolecular interactions between inhibitors and the PDE5 enzyme is essential for elucidating their inhibitory potency and selectivity. The interactions determine the binding affinity of an inhibitor to the enzyme, thereby affecting its therapeutic efficacy. Figure 5 presents the most stable conformations of the investigated compounds within the active site of the protein, while Figure 6 highlights the key intermolecular interactions between the investigated compounds and amino acid residues in the active site.
The most stable conformations of the investigated compounds: SIL (a), TAD (b), AVA (c) and UA (d) in the active site of the PDE5 enzyme. The red α-helices represent the helical regions of the PDE5 enzyme, contributing to the overall stability and structural framework of the protein. The green loops and β-turns depict flexible regions of the protein, which allow conformational changes critical for ligand binding. Carbon atoms are grey, nitrogen atoms blue, oxygen atoms red, hydrogen atoms white, and chlorine atoms green.
Two-dimensional representation of significant intermolecular interactions between the PDE5 enzyme (PDB: 6L6E) and the inhibitors SIL, TAD, AVA, and UA, displaying interatomic distances (Å) obtained after molecular docking studies. Different colors represent various types of interactions, as detailed in the legend.
The analysis of intermolecular interactions between the inhibitors and the PDE5 enzyme indicates the presence of significant hydrogen bonds, π–π interactions, and alkyl and π–cation interactions. These interactions collectively enhance the stabilization of the complexes and suggest potential inhibition (Figure 6).
In protein–ligand complexes, conventional hydrogen bonds are crucial for the proper orientation of inhibitor molecules, facilitating optimal positioning within the enzyme’s binding pocket. A conventional hydrogen bond is formed between SIL and A:TYR 612 through its carbonyl oxygen, facilitating its anchoring in the active site. A hydrogen bond is formed between TAD and A:GLN 817 through the –NH– group of its indole ring, contributing to the molecule’s stabilization. Hydrogen bonds are formed between AVA and amino acids such as A:TYR 612, A:THR 723, and A:LEU 725, thereby ensuring stable binding within the enzyme’s active site and enhancing its affinity for PDE5. A conventional hydrogen bond is established between UA and A:LEU 725, a critical interaction that facilitates the correct positioning and stabilization of UA within the active site.
Besides hydrogen bonds, π–π interactions significantly contribute to the stabilization of the inhibitors. All studied compounds engage in π–π interactions with the aromatic amino acid residue A:PHE 820, which plays a significant role in stabilizing the ligand within the binding pocket of the PDE5 enzyme. The TAD and AVA molecules also form additional π–π interactions with A:PHE 786, which further anchors them in the hydrophobic region of the active site, thus increasing the inhibitors’ affinity for the enzyme.
Furthermore, the π–cation interactions observed in inhibitors such as AVA and SIL, particularly with A:HIS 613, highlight the critical role of electrostatic interactions in stabilizing the complexes. These interactions between positively charged residues and the ligand’s π-system significantly strengthen the inhibitors’ binding affinity to the PDE5 enzyme.
For the stabilization of complexes, π–Sigma contacts are important. All compounds established π–sigma interactions with the A:VAL 782 residue, thereby enhancing the stability of the ligand within the enzyme’s binding site. Furthermore, alkyl/π–alkyl interactions and carbon–hydrogen bonds/π-donor hydrogen bonds between the inhibitors and residues including A:LEU 765, A:ALA 767, A:PHE 786, A:LEU 804, A:ILE 813, A:MET 816, and A:GLN 817 play a crucial role in enhancing the stability of the complexes. The results are supported by thermodynamic parameters, with ΔGinter and ΔGvdw+hbond+desolv indicating the contributions of hydrophobic, van der Waals, and solvation interactions to the stability of the complexes (Table 4). The lowest values for these parameters are obtained for AVA, suggesting strong hydrophobic stabilization due to a high number of π–π, π–sigma, and alkyl/π–alkyl interactions. These intermolecular interactions collectively enhance the stability of the inhibitor complexes with the PDE5 enzyme, allowing each inhibitor to leverage a unique combination of interactions to achieve high inhibitory potency.
3.3. Inhibitory activity toward PDE5 enzyme—molecular dynamic study
This study utilized molecular dynamics (MD) simulations to examine the conformational changes and stability of a receptor in the presence of various ligands. Key parameters, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (Rg), were analyzed to evaluate the impact of ligands on the stability and flexibility of the protein structure. Special emphasis was placed on UA and its inhibitory potential compared to other tested ligands, including AVA, SIL, and TAD.
Analysis of RMSD indicated that all systems attained stable conformations throughout the simulations (Figure 7a). The receptor devoid of a ligand exhibited an average RMSD value of 5.631 nm, whereas the systems with ligands AVA, SIL, TAD, and UA fell within a narrow range of 5.625–5.711 nm. Despite the absence of significant deviations between the systems, UA demonstrated a stable RMSD value of 5.668 nm, comparable to TAD’s 5.625 nm, suggesting that UA does not cause notable conformational changes in the protein.
Comparative analysis of 100 ns molecular dynamics (MD) simulations showing (a) RMSD (root mean square deviation) of the C–Cα–N backbone, (b) Rg (radius of gyration), and (c) RMSF (root mean square fluctuation) for the PDE5 enzyme with and without various ligands (AVA, SIL, TAD, and UA).
The Rg values, indicating the compactness of the protein structure, exhibited minimal variations across the systems (Figure 7b). The receptor in the absence of a ligand exhibited an average Rg value of 4.003 nm, and the presence of ligands did not result in notable alterations in the protein’s overall compactness. The lowest Rg values recorded were for SIL at 3.989 nm and TAD at 3.991 nm, whereas UA exhibited an Rg value of 3.997 nm. Despite the similarity in values, UA exhibits a stable, compact structure that corresponds with its inhibitory potential. The compactness of the protein in the presence of UA may enhance structural stabilization and limit access to active sites, which is essential for inhibitor efficacy. The lack of notable variations in Rg values indicates that UA, similar to TAD and SIL, preserves a stable protein conformation with minimal changes in compactness, thereby reinforcing its inhibitory potential.
The values of RMSF revealed notable differences in the fluctuations of the proteins in the presence of the ligands (Figure 7c). The receptor without a ligand exhibited the highest average RMSF value (1.454 nm), indicating greater flexibility. The ligands SIL and UA significantly reduced protein fluctuations (0.973 nm and 0.829 nm, respectively), suggesting that these ligands stabilize key parts of the protein structure. The ligand that most effectively diminishes protein flexibility is UA, exhibiting the lowest RMSF value of 0.829 nm. This indicates that UA establishes more stable interactions with the protein than the other ligands, specifically AVA (1.403 nm) and SIL (0.973 nm), which is essential for inhibiting protein function. The diminished flexibility of protein chains in the presence of UA may hinder the dynamic alterations required for protein activity, significantly increasing its inhibitory potential.
In comparison to other ligands, UA demonstrates significant efficacy as a protein stabilizer, as indicated by its lowest RMSF values and stable RMSD. This combination suggests that UA reduces protein flexibility while preserving a stable conformation, characterized by minimal fluctuations and a lack of destabilization. Based on these characteristics, it can be concluded that UA exhibits significant potential as an inhibitor, especially when compared to AVA and SIL, which demonstrate greater fluctuations and relatively less stable RMSD values. Additionally, UA exhibits comparable Rg values and stability in conformational behavior when contrasted with TAD. The UA ligand demonstrates a marginally superior capacity to stabilize the protein, positioning it as a more viable candidate for subsequent research regarding inhibitory potential.
The energy analysis of the binding of SIL, TAD, AVA, and UA to the PDE5 enzyme reveals important trends in complex stability and the contributions of different energetic components. Van der Waals interactions play a pivotal role in stabilizing all complexes. The strongest van der Waals interactions (EVDW) are obtained for TAD and AVA, with values of −48.23 and −48.61 kcal/mol, respectively. In comparison, SIL and UA demonstrate favorable, but less pronounced, interactions of −38.96 and −22.18 kcal/mol, respectively (Table 5). These results suggest that, although somewhat less potent, UA still effectively interacts with the enzyme’s hydrophobic pocket.
Free energy decomposition analysis (MM-GBSA) for PDE5 enzyme complexed with different ligands (SIL, TAD, AVA, and UA): van der Waals (EVDW), electrostatic (Eele), polar and nonpolar solvation energies (Gpolar and Gnonpolar), gas-phase binding energy (ΔGgas), solvation free energy difference (ΔGsolv), and binding free energy (ΔGbind) in kcal/mol
| Compound | EVDW | Eele | Gpolar | Gnonpolar | ΔGgas | ΔGsolv | ΔGbind |
|---|---|---|---|---|---|---|---|
| PDE5-SIL | −38.96 | −18.89 | 51.78 | −5.49 | −57.85 | 46.29 | −11.56 |
| PDE5-TAD | −48.23 | −12.71 | 35.84 | −6.01 | −60.95 | 29.82 | −31.12 |
| PDE5-AVA | −48.61 | −9.76 | 50.05 | −6.79 | −58.36 | 43.25 | −15.11 |
| PDE5-UA | −22.18 | 3.69 | 10.75 | −3.27 | −18.49 | 7.48 | −11.01 |
Electrostatic interactions further contribute to the overall stability of the complexes. The most favorable negative energy electrostatic values (Eele) are obtained for SIL and TAD, indicating strong electrostatic interactions with amino acids in the active site. In contrast, UA shows a positive electrostatic energy value of +3.69 kcal/mol, suggesting an alternative interaction mechanism that relies less on electrostatic stabilization and more on other interaction types within the complex. All compounds display expected unfavorable polar solvation energies (Gpolar), attributed to the exposure of polar groups to the solvent during complex formation. Polar solvation energies are similar for AVA and SIL, around 50 kcal/mol, while TAD exhibits a more favorable value of 35.84 kcal/mol. Notably, UA shows a much lower polar solvation energy (10.75 kcal/mol), suggesting its binding is less impacted by solvent exposure compared to the other compounds (Table 5). This implies that, despite lacking strong electrostatic interactions, UA’s interaction with the solvent enhances its binding affinity.
The total binding energies (ΔGbind) provide a comprehensive evaluation of the stability of each complex. The lowest binding energy is obtained for TAD (−31.12 kcal/mol), indicating the strongest interaction with PDE5. In comparison, AVA and SIL exhibit slightly higher binding energies of −15.11 and −11.56 kcal/mol, respectively. The UA ligand shows a competitive binding energy of −11.01 kcal/mol, suggesting that, despite its unique interaction profile, it remains an effective PDE5 inhibitor relative to established pharmaceutical compounds. The ability of UA to form a stable complex with PDE5 highlights its potential for further exploration, particularly in terms of modifications that could enhance its inhibitory effectiveness (Table 5).
4. Conclusion
This study presents an analysis of the inhibitory potential of usnic acid (UA) in comparison to established PDE5 inhibitors, including sildenafil (SIL), tadalafil (TAD), and avanafil (AVA), for the treatment of pulmonary arterial hypertension (PAH). The ADMET analysis indicated that UA possesses a distinct pharmacokinetic profile characterized by high plasma protein binding and an extended half-life, potentially facilitating its prolonged retention in the body. Its absorption is limited relative to synthetic inhibitors, yet it exhibits considerable permeability across lipid membranes. Crucially for its potential therapeutic application, UA demonstrated low toxicity, exhibiting minimal risks of cardiotoxicity and neurotoxicity.
The combination of natural compounds such as UA with synthetic PDE5 inhibitors offers a promising approach for achieving synergistic effects in the treatment of PAH. This combination may improve the effectiveness of standard PDE5 inhibitors and potentially lower the necessary dosage of synthetic drugs, thereby reducing side effects. The antioxidant properties exhibited by UA may address oxidative stress and inflammation, both of which are crucial in the progression of PAH, thereby offering a multi-targeted therapeutic strategy. Exploring the synergistic potential of UA in conjunction with SIL or TAD may lead to the development of novel combination therapies that enhance both safety and efficacy.
Molecular docking results indicate that UA establishes stable interactions with the PDE5 enzyme, exhibiting a binding free energy of −10.12 kcal/mol, which is comparable to SIL at −11.01 kcal/mol and AVA at −11.39 kcal/mol. Significant intermolecular hydrogen bonds and π–π interactions are established by UA, which are essential for the stabilization of the complex. Molecular dynamics simulations indicated that UA decreases the flexibility of the PDE5 enzyme, as evidenced by low RMSF values, suggesting its potential role as an enzyme stabilizer.
The MM-GBSA analysis revealed further details regarding the binding energy, indicating that UA exhibited a total binding energy of −11.01 kcal/mol, which is comparable to that of SIL and AVA. Van der Waals interactions were crucial for the stability of the complex with UA, and favorable solvation interactions further enhanced this stability. While electrostatic interactions were less significant than those of synthetic inhibitors, UA demonstrated considerable potential for stable interaction with PDE5 through alternative mechanisms.
The findings indicate that UA possesses significant potential as a natural PDE5 inhibitor, along with antioxidant and anti-inflammatory properties that may enhance its therapeutic application in PAH. The results of this study indicate that UA may serve as a complementary or adjunctive therapy in the management of PAH. Future research must prioritize in vivo studies to evaluate the bioavailability, efficacy, and safety of UA in clinical contexts. Long-term studies are essential to assess the potential for chronic toxicity and interactions with current PAH medications. Furthermore, clinical trials that integrate UA with synthetic PDE5 inhibitors could investigate the potential of UA’s antioxidant and anti-inflammatory properties to enhance treatment outcomes and improve patients’ quality of life.
Acknowledgement
The authors wish to acknowledge the Ministry of Science, Technological Development and Innovation of the Republic of Serbia and the Ministry of Education of the Republic of Serbia (Agreements No. 451-03-66/2024-03/200378) for financial support.
Declaration of interests
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.
Supplementary materials
Supporting information for this article is available on the journal’s website under https://doi.org/10.5802/crchim.413 or from the author.