PGs are biological active molecules with wide range of activities and they exist in many plant genus such as Veronica, Magnolia and Forsythia [15]. Out of 2370 natural products, we identify 3 PGs (Fig. 1) with high docking score which interact with amino acids at the binding site, including the catalytic residue His-41. Verbascoside had the highest docking score (-9.1317, RMSD = 1.51) followed by forsythiaside (-8.0889, RMSD = 2.10) and isoacteoside (-7.7302, RMSD = 1.68). The sugar part of forsythiaside A was able to bind to the catalytic residue His-41 and Asp-187 while one aromatic ring formed a hydrogen bond with Thr-190 and arene-H bond with Pro-168. Isoacteoside’s sugars were able to form many hydrogen bonds at the binding site and one aromatic ring formed a hydrogen bond with Glu-166 and the second aromatic ring formed arene-H bond with Gln-189. Verbascoside formed 3 hydrogen bonds with the Thr-24, Asn-142 and Phe-140 and one arene-H bond with Gln-189 (Fig. 2).
To analyze the selected 3 phenylethanoid glycosides properties, SwissADME was used (Table 1) [31–33]. All the three molecules are soluble in water and have low gastrointestinal absorption. We used Swiss Model (TargetPrediction) to predict the biological targets of the selected 3 PGs. Swiss model predicted the possible molecular targets along with the actual targets. Targets are ranked according to a score that combines both 2D and 3D similarity values with the most similar known active molecules to the query molecule [20, 21]. As shown in Fig. 3, PGs exhibit kinase and protease inhibitor activity according to the Swiss model. It has been reported that isoacteoside inhibits the expression of protease protein; matrix metalloproteinase 2, 12 and 13 with IC50 around 10 µM [34]. Additionally, forsythiaside and verbascoside inhibit protein kinase C alpha with IC50 1.9 µM and 68 nM respectively [35, 36].
Table 1
SwissADME results of top three selected compounds with ChemBL IDs.
Compound Names
(ChemBL ID)
|
Clinical Trail
|
Lipinsk’s Rules
|
Structural Alerts
|
Pharmacokinetics
|
Water Solubility
|
Forsythiaside A (CHEMBL504363)
|
-
|
3 violations MW > 500,
N/O > 10, NH/OH > 5
|
PAINS1: 1 alert (Catechol A), Brenk2: 2 alerts (catechol, Michael acceptor1)
|
Low GI absorption
|
Soluble; iLOGP3 2.15
|
Isoacteoside (CHEMBL504873)
|
-
|
3 violations MW > 500,
N/O > 10, NH/OH > 5
|
PAINS: 1 alert(Catechol A), Brenk: 2 alerts (catechol, Michael acceptor 1)
|
Low GI absorption
|
Soluble; iLOGP 2.15
|
Verbascoside (CHEMBL444478)
|
Phase II
|
3 violations MW > 500,
N/O > 10, NH/OH > 5
|
PAINS: 1 alert(Catechol A), Brenk: 2 alerts (catechol, Michael acceptor 1)
|
Low GI absorption
|
Soluble; iLOGP 2.15
|
1Pan-assay interference compounds, 2Brenk alerts to toxic groups, 3Octanol/water partition coefficient. |
A recent study reported that protein kinase C alpha and beta inhibitors have a pro-apoptotic effect in nucleated cells, and this creates a hostile environment for intracellular parasites including viruses. This can enhance the antiviral activity of antivirus drugs [37]. The dual activity (kinase/protease) of PGs could make them very promising as anit-COVID-19 agents.
From molecular dynamic simulations we assessed the stability of the protein-ligand bond. A label was associated to each different natural small molecule to improve the readability of the data analysis. Table 2 describes the PGs nomenclature adopted herein.
Table 2
Compounds associated labels.
Compounds
|
Associated Label
|
Forsythiaside A
|
L1
|
Isoacteoside
|
L2
|
Verbascoside
|
L3
|
As shown in Fig. 4, the root-mean-square deviation (RMSD) measurement was carried out on each of the three ligands. Interestingly each ligand reaches a steady configuration at different times within approximately 50 ns. However, L1 stands out for reaching a very steady configuration with respect to the other ligands. This implies a lower noise and mean value of the L1 RMSD. In the Supporting information we are showing that also the protein’s RMSD is associated to lower noise when interacting with L1 (Fig. S1).
To further characterize the system, we decided to compute the root mean square fluctuation (RMSF) of atomic positions in the trajectory. In Fig. 5 we consider a subset of residues in the active site corresponding to the “pocket” of interaction. Here, it is evident that different ligands interact in different characteristic ways. There are different moieties of the pocket, highlighted in Fig. 5, where the fluctuation of the atoms depends strongly on the binding ligand. However, the fluctuations in the backbone of the active sites are always limited to less than one Angstrom. Which means that none of the three scenarios implies a detrimental ligand-protease interaction.
The ligand themselves interact with the active site mainly by means of their aromatic moieties. In order to distinguish the two aromatic rings, characteristic of each ligand, we considered their bonding to the central THF heterocyclic ring of Fig. 1. Ring A is defined as the one connected to the central THF heterocyclic ring via a carboxylic group; whereas ring B is not connected via a carboxylic group. Recording the RMSF of the aromatic moieties, reported in Fig. 6, L1 turns out to be the most stable compound of the three. In Fact, L2 and L3 both show one of the aromatic rings to be fluctuating to twice the extent of the other which allows L2 and L3 to explore different configurations.
This speculation is backed up by the gyration radius data (Figure S2) of the ligands and the distance between the center of mass of ring A and B. Such a distance, shown in Fig. 7, is indicative of the L1 aromatic rings constraining. Nonetheless, Fig. 7 shows the relative angle between the aromatic rings. This clarifies that L1 does not explore as many configurations as the other ligands, owing to its enhanced bond stability. At this point, a mere analysis of the average number of contacts between the ligands and the active site. Figure S3, is not sufficient to determine which structure would be more reliable. Detailed contact map is presented in Figure S4. Binding free energy calculations allow the measurement of binding strength in protein–ligand or protein–protein complexes [38, 39].
Table 3
MMPBSA energy calculated from the last 400 ns of MD simulation.
|
Van der Waal energy (kJ/mol)
|
Electrostatic energy (kJ/mol)
|
Polar solvation energy (kJ/mol)
|
SASA energy (kJ/mol)
|
Binding energy (kJ/mol)
|
L1
|
-253 ± 24
|
-64 ± 25
|
221 ± 27
|
-25 ± 1
|
-121 ± 19
|
L2
|
-246 ± 24
|
-109 ± 22
|
279 ± 23
|
-25 ± 2
|
-101 ± 22
|
L3
|
-174 ± 25
|
-69 ± 35
|
211 ± 45
|
-21 ± 2
|
-53 ± 18
|
The molecular mechanics of the MM/PBSA method, even though is not as accurate as more computationally intensive methods such as thermodynamic integration [40], has been shown to be satisfactory for the work of both computational and experimental researchers. Last 100 ns of the MD trajectories was used for the binding free energy calculations for protein/ligand complexes and the result reported in Table 3. According to the result, L1 has a statistically significant lower binding energy in compare with other ligands.