2.1 Virtual screening

2.2 Molecular dynamics stability study

The movements of the protein docked to the predicted top inhibitory compound were analyzed and compared with the apo form by MD simulation in the Schrodinger Desmond module (Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, 2021). The simple point charge model was used, and the force field OPLS2005 was applied. Counter ions were added to neutralize charges (42Na⁺ and 42Cl^- ions), and a 0.15 M NaCl level was maintained to simulate the physiological condition of the human system. At the NPT ensemble temperature of 300 k and 1,01325 bar of pressure strain, the MD simulation was run for 50 ns. The results of MD included RMSD, RMSF, and protein–ligand interactions (Ouassaf et al., 2021). The SigmaPlot (SigmaPlot Version 14.0, Systat Software, San Jose, CA. 2018) program was used for creating comparative graphs.

3.1 Virtual screening by benchmarking molecular docking

OS from the OliveNet^TM library underwent virtual screening for their potential to inhibit the SARS-CoV-2 PLpro (PDB 6WX4) along with benchmarking known active molecules and decoys from DEKOIS 2.0 database. The ligands were flexible throughout the docking run and moved freely inside the binding sites, and the protein was rigid. Three secoiridoids, Oleuropein pentamer, Oleuropein tetramer, and Oleuropein trimer (MW = 2684.67, 2150.11, and 1613.58 g/mol, respectively), were excluded from the study because of their unfavorable intramolecular interactions during energy minimization. ADV and AD were used for benchmarking the accuracy of prediction, wherein both programs apply empirical free energy scoring. They differ in the search algorithms, where AD uses the Lamarckian genetic algorithm while ADV uses an optimized gradient global and local conformational search (Chang et al., 2010). The binding energies predicted (kcal/mol) by the docking programs were used for active ligand rankings, as shown in Fig. 1.

Close

Fig. 1

Results of benchmarked docking of olive secoiridoids against SARS-CoV-2 PLpro. Predicted binding energy scores (kcal/mol) plotted against the ranks of the active molecules. (Set 1, 2 = AutoDock Vina, AutoDock). The rug plot above the graph shows the positions of active molecules.

Export to PPT

3.2 Performance assessment metrics

The rankings obtained from the docking results were used to construct ROC, EF, and PC curves from which the assessment metrics were calculated. Fig. 2A is the ROC curve at a maximum threshold with ADV ROC-AUC: 0.878; AD ROC-AUC: 0.847 indicating the successful application of scoring functions in differentiating the active molecules from decoys, and both programs have shown relatively equal capacity in rightly ranking the ligands (Empereur-Mot et al., 2015). Fair values of ADV BEDROC: 0.608; AV BEDROC: 0.510 indicate both programs are good in early recognition of active molecules, with ADV showing a marginal improvement in the results (Empereur-Mot et al., 2016). RIE is not reported because BEDROC is the normalization of RIE (Empereur-Mot et al., 2016). TG scores ADV TG: 0.511; AD TG: 0.495, that are above 0.2 suggest a significant relationship between the score variations and active molecules detection. EF_max calculated from ROC of Fig. 2A is 10.35 indicating satisfactory early recognition by both programs. Fig. 2B is the enrichment curve that enables the quantification of early recognition ability by visualization of a true positive fraction over a specific fraction (partial) of the entire ligand set. We analyzed the partial metrics at 1 %: pAUC ADV: 0.588; pAUC AD: 0.505, pTG ADV: 0.684; pTG AV: 0.527, EF_1% ADV: 6.47; EF_1% AV: 3.88 and partial metrics at 5 % pAUC ADV: 0.558; pAUC AD: 0.499, pTG ADV: 0.545; pTG AV: 0.418, EF_5% ADV: 4.66; EF_5% AV: 3.62 demonstrating that AutoDock Vina has performed better than AutoDock in predicting the actives. Fig. 2C, PC visually displays the better than random retrieval of actives by both programs, and ADV has a clear edge over AD, with a probability of activity score p(act) at 1 % threshold for ADV: 0.691; p(act) AD: 0545. Fig. 2D displays the ability of the programs to retrieve the known actives: ROC-AUC ADV: 0.905; ROC-AUC AD: 0.895, TG ADV: 0.676; TG AD: 0.655, BEDROC ADV: 0.627; BEDROC AD: 0.610 indicating both programs were equally excellent (Empereur-Mot et al., 2016). The ligand rankings by ADV and AD showed a moderate correlation with Kendall's correlation (τ) = 0.581 (Chang et al., 2010). The ligands placed in the top 1 % rank of the screened molecules were retrieved from Fig. 2B by choosing the set fraction at 0.01 on X-axis. Eight active compounds were retrieved from the 2 docking programs as follows: ADV: OS-3, known actives-2, decoys-3; AD: OS-3, known actives-1, decoys-4. Table 1 shows the enriched top 1 % OS, with two same molecules occupying the list, and Demethyloleuropein aglycone was identified by both software as the best compound with the lowest binding energy (ΔG: ADV = -10.5, AD = -9.9 kcal/mol), hence was chosen for MD analysis. Fig. 3 represents the active conformations of the top 1 % ranked OS bound to PLpro retrieved from AutoDock Vina docking.

Close

Fig. 2

Quantification of docking performance for benchmarking virtual screening of olive secoiridoids as potential SARS-CoV-2 PLpro inhibitors. A) ROC curves B) Enrichment curves for all the screened molecules C) Predictiveness curves D) Enrichment curves for early identification of known active molecules. Red: AutoDock Vina; Green: AutoDock; Dashed lines: Random retrieval.

Export to PPT

Table 1 Structures of the top 1% ranked olive secoiridoids, predicted as SARS-CoV-2 PLpro inhibitors by benchmarked docking.

Olive secoiridoids predicted by AutoDock Vina in the top 1 % rank	Olive secoiridoids predicted by AutoDock in the top 1 % rank
Rank 1 Demethyloleuropein aglycone	Rank 4 Demethyloleuropein aglycone
Rank 3 Oleuroside	Rank 5 Oleuropein
Rank 6 Oleuroside-10-carboxylic acid	Rank 6 Oleuroside

Close

Fig. 3

AutoDock Vina retrieved docked poses of the top 1% ranked olive secoiridoids (shown as sticks, red: Demethyloleuropein aglycone, yellow: Oleuroside, magenta: Oleuropein, green: Oleuroside-10-carboxylic acid) inside the binding site of SARS-CoV-2 PLpro (Connolly surface).

Export to PPT

3.3 Molecular dynamics stability analysis

Results of MD analysis of Demethyloleuropein aglycone, the top-ranked ligand in complex with PLpro, are shown in Fig. 4, indicating high stability achieved after ligand binding. The RMSD graph (Fig. 4A) exhibits the ligand-bound protein, apo-protein, and ligand conformational changes. The ligand-bound protein reached a stable conformation with a global minimum equal to the apo-protein and maintained stability from 30 ns to 50 ns with RMSD < 2 Å and RMSF < 2.5 Å (Fig. 4B). Low RMSF indicates a less significant fluctuations of the PLpro residues. Fig. 5 is a detailed illustration of the Demethyloleuropein aglycone interactions with PLpro during MD. Arg166 is a positively charged residue located in binding site 2 of PLpro critical for inhibition. Demethyloleuropein aglycone was predicted to form non-covalent electrostatic interactions (67 % of the 50 ns) through the oxygen of phenolic hydroxy moiety, charged interactions through keto oxygen of ester via water bridge (62 %), and Pi-cation interactions (81 %) through the aromatic ring with Arg166, suggesting a potent non-covalent inhibition (Osipiuk et al., 2021). MD results also demonstrated that the ligand occupied site 3 comprising the BL2 loop with Tyr273 and Asp302. Asp302, in a negatively charged state, displayed two strong hydrogen bonds with the H atom of the two phenolic –OH groups for 98 % of the simulation time. Tyr273 exhibited Pi-Pi stacked hydrophobic interactions for 56 % of MD run time. The two phenolic hydroxy groups of Demethyloleuropein aglycone were engaged in hydrogen bonds with Asp302 and Arg166 for a long duration, suggesting that these groups have a significant influence on the nature of binding and stability of interactions with PLpro. Arg166 contributes to PLpro peptidase activity by cleaving the amide linkage between ubiquitin and ISG15 interfering with the human viral defense, thus its interactions with ligand imply its inhibition, delineating the preventative action of OS (Shen et al., 2022). Interestingly the results also demonstrate that the OS occupied the recently identified Pro248-BL2 groove (Shen et al., 2022). Engagement of multiple binding sites leads to cooperative binding resulting in better binding affinity than binding to a single site, justifying potent inhibition (Shen et al., 2022).

Close

Fig. 4

Molecular dynamics stability analysis of Demethyloleuropein aglycone bound SARS-CoV-2 PLpro. A) Comparative RMSD plot B) RMSF of bound PLpro (blue) and apo PLpro (black).

Export to PPT

Close

Fig. 5

Schematic illustration of Demethyloleuropein aglycone atom interactions with SARS-CoV-2 PLpro residues during molecular dynamics.

Export to PPT

3.4 Binding free energy calculation

The binding free energy (MM-GBSA, ΔG_binding, kcal/mol) of Demethyloleuropein aglycone bound PLpro complex was −94.54 ± 6.05 contributed by stabilizing energies: lipophilic energy (ΔG_lipo = -16.78 ± 1.23), Van der Waals energy (ΔG_vdW = -34.68 ± 2.17), coulomb energy (ΔG_coulomb = –22.56 ± 5.54), and destabilizing energies: generalized Born electrostatic solvation energy (ΔG_{Solv GB} = 22.45 ± 2.78), covalent binding energy (ΔG_covalent = 2.13 ± 1.90). Stabilizing energies contributed more than destabilizing energies suggesting a favorable bound conformation of the protein and ligand (Wang et al., 2019).