2.1.1 Pharmacophore model generation

The main goal of the study was to generate a receptor structure-based pharmacophore model depicting the characteristic steric and electronic features of ligands, facilitating their binding to multiple binding sites on PLpro. Hence, we chose the recently reported 3D X-ray crystalline structures of PLpro receptor bound to potent inhibitors available in PDB ( https://www.rcsb.org). The PDB IDs of 3D structures of PLpro are: 7LBS (native ligand: XR8-24), 7LOS (native ligand: XR8-65), 7LLZ (native ligand: XR8-69), and 7LLF (native ligand: XR8-83) (Shen et al., 2022). These PDB structures were directly fetched into LigandScout 4.4.8. software using their codes (Wolber & Langer, 2005). LigandScout can be accessed through a valid license on Inte: Ligand [http:// www.inteligand.com/ligandscout/download.shtml]. All four receptor-ligand structures were visualized in structure-based modeling perspective mode, checked for missing atoms in the ligand structures, and four pharmacophore models were generated individually. To select a suitable pharmacophore model for virtual screening, all 4 models were evaluated against a set of known actives derived from literature (Shen et al., 2022). Known active compounds were carefully chosen based on their experimentally determined PLpro inhibition constant (IC₅₀) ranging from 0.1 to 5.7 μM provided possess structural characteristics necessary for exhibiting multiple binding site modalities. Pharmacophore models were evaluated against the actives using the screen pharmacophore mode under the structure-based perspective. During the process, every model was optimized manually by varying the degree of tolerance of pharmacophoric features under consideration. We varied the features by increasing or decreasing the tolerance, and weightage, adding or excluding exclusion volume spheres, and omitting a minimum number of features to increase the sensitivity of the models (Opo et al., 2021). The optimum pharmacophore model is the one that identified maximum hits among the actives, which was considered for further study (Luo et al., 2021).

2.1.2 Pharmacophore model validation

The pharmacophore model obtained by optimization in Section 2.1.1 has to be validated by evaluation of its ability to distinguish active compounds from decoys (Temml et al., 2017). The structures of 23 known actives identified from the literature were drawn in ChemDraw Ultra 10.0 and the smiles were copied to OpenBabelGUI ( https://openbabel.org/wiki/OpenBabelGUI) to convert into a single sdf file of actives (O’Boyle et al., 2013). Decoys for PLpro are available in the DEKOIS.2.0 database (Bauer et al., 2013). These property-matched decoys were recently generated and validated for PLpro binding (Ibrahim et al., 2020). The decoy set consisted of 720 compounds in a single sdf file prepared for screening. Actives and decoys sets were processed to ldb (library in local database) files in LigandScout using create screening database module, which is the only accepted format by the software. Actives and decoys sets were loaded in the ldb formats. The 7LBS model was optimized for sensitivity by increasing the tolerance of positive ionizable area, two of the hydrophobic interactions of alkyl groups to 0.3 Å, and decreasing the tolerance of hydrogen bond donor and acceptor features to 0.3 Å. Exclusion volume spheres were also added. The optimized pharmacophore model was copied to the screening perspective. The pharmacophore-fit score was set as the scoring function, match all query features was the screening mode, get the best matching conformation was the retrieval mode, and maximum number of omitted features was 2. Boolean expression was set to choose the optimized pharmacophore model depicted as {1} and the screening was initiated. After the run, results were analyzed in LigandScout by retrieving the receiver operating characteristic curve (ROC) constructed based on the ratio of actives and decoys identified by the pharmacophore model. The area under the curve (ROC-AUC) is used to estimate the detecting power of the model A pharmacophore model is said to be valid if it exhibits an AUC > 0.5 and will have excellent detective capacity if it is near 1. Early enrichment factors (EF) at 1, 5, 10, and 100 % also reveal the early detection of actives and support model validation (Kirchmair et al., 2008).

2.1.3 Pharmacophore model aided virtual screening

The CMNPD facilitates the download of chemical structures in sdf format ( https://www.cmnpd.org/). The entire database is available for download and all compounds were downloaded. The sdf file contained 47,451 marine natural compounds which were then converted to a ldb file in LigandScout during which 1,320 duplicates were removed. LigandScout generates 25 conformers for each compound. The iconFast mode for fast, high throughput generation of high-quality conformers was chosen to create the database for screening. The optimized pharmacophore model was fed to the screening perspective. CMNPD and decoy libraries were loaded to the screening mode and screening was conducted adopting the scoring, query, matching, and omitted features as mentioned in Section 2.1.2. The results were visualized in the library view, ROC was constructed, and the validity was assessed based on ROC-AUC. The data of resultant initial hits were saved as a sdf file and further refinement in OpenBabelGUI applying a molecular weight filter of ≤ 500 was carried out. The compounds with molecular weight ≤ 500, known from OpenBabelGUI results were chosen manually to align to the pharmacophore model in LigandScout. This refined hit library was stored in sdf and xls formats for the next stage of molecular docking studies. The sdf files can be saved to retain 2D or 3D coordinates of the selected conformation of the ligands. The minimum energy conformer was chosen for all ligands and the sdf file was written with 2D and 3D coordinates separately.

2.2.1 Processing of PLpro and ligand structures

The PLpro structure used for docking was 7LBS which contains a native ligand. MGL tools was used to prepare the target structure. The native ligand, water molecules, and other non-standard residues were removed. Gasteiger charges and polar hydrogens were added. The energy-minimized target structure was saved in pdbqt format. For the hits and decoys, the sdf file containing the 3D coordinates of the minimum energy structure obtained in LigandScout was used. These structures in sdf were then converted to individual pdbqt files applying OpenBabelGUI. The native ligand XR8-24 structure saved in PDB was converted to pdbqt in OpenBabelGUI and was included in the docking study for comparison with ligands.

2.2.2 Molecular docking

The binding site of PLpro was identified according to the literature by highlighting the residues Asp164, Arg166, Glu167, Pro248, Tyr264, Gly266, Tyr268, Tyr273, and Pro299 (Shen et al., 2022). This created a grid box with dimensions of 9 Å X 8 Å X 35 Å, and size of 30 Å in X, Y, and Z coordinates. This docking grid covered only the binding sites and not the whole surface of the target. Same docking coordinates were applied for docking in both docking engines. All other docking parameters were kept at default for both docking engines. Chiefly, in AutoDock the genetic algorithm runs were 10, and in AutoDock Vina the number of conformers generated was 10 at exhaustiveness of search 8 with an energy difference of 3 kcal/mol. The binding affinities in terms of binding energy (ΔG in kcal/mol) for the hits and decoys were predicted by both docking engines. The conformer with minimum binding energy, the one with the highest negative value was chosen for binding interaction analysis, individual, and consensus scoring. PLpro bound to the predicted best inhibitor was visualized for the binding site engagement in UCSF Chimera 1.13.1. ( https://www.cgl.ucsf.edu/chimera/download.html) (Pettersen et al., 2004). Binding interaction analysis was performed in Discovery Studio 2016 ( https://discover.3ds.com/discovery-studio-visualizer-download).

2.2.3 Benchmarking molecular docking

Benchmarking procedures to assess the predictive potential of docking engines based on individual scoring and consensus scoring was applied.

3.1.1 Pharmacophore model generation

We used LigandScout 4.4.8. software for structure-based pharmacophore model generation. LigandScout is advanced molecular design software that functions on a graphical user interface (GUI) and is robust in creating 3D and 2D structure-based and ligand-based pharmacophore models. LigandScout automatically generates 3D pharmacophore models applying unique algorithms that deduce and interpret ligands and their target structure environment from PDB files. It also enables visualization, alignment, and manual refining of pharmacophore models to understand and optimize models for improved prediction accuracy (Seidel et al., 2017).

Recent PDB depositions of PLpro 3D structures with IDs 7LBS, 7LOS, 7LLZ, and 7LLF were chosen for pharmacophore model generation. We chose these PLpro structures with ligands because the chemistry of ligands is proposed to encourage cooperative binding that strengthens the binding to multiple binding sites resulting in a potent inhibitory activity (Shen et al., 2022). Each structure has a ligand bound to it that guided the pharmacophore generation. Four pharmacophore models containing the essence of the corresponding ligand characteristics were generated by applying default settings. We hand-picked a set of 23 active compounds having the chemistry supporting cooperative binding to PLpro and experimental inhibition constant, IC₅₀ in the range of 0.1 to 5.7 μM from the literature to evaluate the selectivity index, i.e., the number of actives retrieved to the total number of actives (Luo et al., 2021; Shen et al., 2022). During this process of evaluation, the pharmacophoric features of all 4 models were manually optimized to improve the number of actives retrieved. After several attempts, the pharmacophore model derived from 7LBS was ranked 1 because it exhibited a better selectivity profile than other models, as shown in Table 1. The final optimized 7LBS model was obtained by choosing the screening mode to match all pharmacophoric features, decreasing the tolerance of HBD, and HBA by 0.3 Å, increasing the tolerance of two hydrophobic features of the alkyl groups by 0.3 Å, increasing the tolerance of PI group by 0.3 Å, and permission to omit a maximum of 2 features. We also added exclusion volume coats to the model. Exclusion volume spheres increase the sensitivity of pharmacophore models because these areas are forbidden for the ligand to interact and they pose a stereochemical constraint for ligand binding (Spitzer et al., 2010). Hence, the exclusion volume spheres in our refined pharmacophore model present a barrier to binding to non-specific ligands (Seidel et al., 2010). Fig. 2 represents the structure-based pharmacophore model derived from the structure of PLpro-7LBS and its complementary ligand XR8-24.

Table 1 Pharmacophore models and their rankings: blue (positive ionization centre); green (hydrogen bond donor); red (hydrogen bond acceptor); yellow (hydrophobic interactions).

Close

Fig. 2

A) 7LBS-XR8-24-based pharmacophore model at the PLpro active site B) Optimized pharmacophore model in 3D C) Optimized view of the model with exclusion volume coats D) 2D depiction of the model showing the interacting residues of PLpro.

Export to PPT

3.1.2 Pharmacophore features analysis

Fig. 2A provides the complete view of native ligand XR8-24 bound to the active site of PLpro and the native pharmacophore model aligned to the XR8-24 structure. Fig. 2B is the optimized pharmacophore model in 3D and Fig. 2C shows the optimized model with the exclusive volume coats. Fig. 2D is the two-dimensional depiction of the optimized pharmacophore model with reference points, i.e., the interacting amino acid residues of PLpro. The model derived from the key structural elements of XR8-24 in a 3D chemical space consisted of 9 features. According to the model in Fig. 2B, a positive ionization centre, 2 hydrogen bond donors, 2 hydrogen bond acceptors, and 4 hydrophobic features in a particular geometry are essential for the ligand to recognize multiple binding sites on PLpro ensuring binding cooperativity and potent inhibition. Exclusion volume spheres in Fig. 2C contemplates the steric nature of the target’s binding site and confer additional stereochemical constraint for binding of ligands enhancing ligand selectivity by the target macromolecule. It is impressive to view the model in 2D because it is a simplified version of 3D. In the 2D depiction of the model, Fig. 2D shows the native ligand structure besides the deduced pharmacophoric characteristics. Typically, pharmacophore models consider the binding site characteristics and its geometry while mapping the pharmacophore features of a ligand and the 2D model proposes the interacting binding site residues. In the PLpro-XR8-24-based model, the interacting residues are Leu162, Asp164, Glu167, Tyr264, Tyr268, Gln269, Tyr273, and Thr301. A viewer can understand from the 2D diagram that the positive ionization centre is an ionized or protonated amine proposed to interact with Asp164 through a hydrogen bond and ionic bond. Any ligand to bind to PLpro may require keto oxygen providing the scope of hydrogen bond acceptor atom for two hydrogen bonds, as specified in the model. There is a need for a hydrogen bond donor group like –NH for any ligand to leverage binding to multiple sites on PLpro. Four hydrophobic areas, in particular, phenyl rings or alkyl groups in the ligand shall provide excellent binding and PLpro inhibition according to the model. Therefore, the generated model supports the study with knowledge about ligand’s structural elements crucial for the predicted interaction pattern with PLpro.

3.1.3 Validation of the pharmacophore model

Validation of the predicted pharmacophore model is critical to determine its effectiveness and reliability in distinguishing active compounds from decoys. Certain quality indicators like ROC, AUC-ROC, and EF are universally accepted for pharmacophore model validation (Kaserer et al., 2015; Kirchmair et al., 2008). LigandScout enables model validation in screening mode when a set of actives and decoys are loaded into the screening window. The 23 active compounds explained in Section 2.1.1 were utilized for model validation. DEKOIS 2.0 library provided 720 decoy structures for screening against the known actives. The actives and decoys libraries were screened using the pharmacophore template created in LigandScout. The predicted pharmacophore model was able to identify 20 out of 23 actives and 37 out of 720 decoys. All decoys had only one hydrophobic feature matched with the model and a pharmacophore-fit score ranging from 15.22 to 15.63. The native ligand XR8-24 was retrieved at rank 1 with all matched features and a pharmacophore fit score of 97.18 indicating the accuracy of the model. Table S1 has the details of the total actives and decoys, retrieved actives and decoys during model validation. ROC for model validation is a plot of specificity vs sensitivity. Specificity is plotted on the X-axis indicates the false positives or decoys retrieved, hence ROC close to zero on the X-axis indicates the high accuracy of the model. Sensitivity, on Y-axis, represents the true positives or actives retrieved, therefore, the higher the value on X-axis, the higher the accuracy. The dotted line in ROC represents random prediction. ROC closer to the axis-Y and the greater the AUC-ROC, the greater the prediction accuracy (Kirchmair et al., 2008; Luo et al., 2021). Fig. 3 is the ROC curve for pharmacophore model validation that reveals the model’s predictive ability is well above random prediction. AUC-ROC (AUC_1% = 1.00) and EF values (EF_1% = 24.3) on the curve indicate our pharmacophore model has high predictive accuracy, sensitivity, and specificity. Thus, it is a reliable model for the virtual screening of CMNPD.

Close

Fig. 3

Pharmacophore model validation. ROC metrics indicate the accuracy and reliability of the model.

Export to PPT

3.1.4 Pharmacophore model-aided virtual screening

The prepared CMNPD for virtual screening was made up of 46,131 marine natural compounds. The decoy set containing 720 decoys was used against the CMNPD in virtual screening to evaluate the reliability of the results. So, the final dataset for virtual screening had a total of 46,851 compounds. LigandScout has an impressive protocol for screening large databases by applying the validated pharmacophore template. The screening protocol attempts to search ligands by matching all pharmacophoric features preferably, when it fails, it omits a maximum of two features and refines the search on each ligand until it gets a single pharmacophore feature match. The pharmacophore model-aided virtual screening of CMNPD in LigandScout yielded 66 hit compounds. Table S2 has information about 66 hit compounds including their pharmacophore-fit scores and matching pharmacophoric features. In a nutshell, virtual screening retrieved a total of 103 hits, which included 66 hits and 37 decoys. The AUC-ROC at 1, 5, 10, and 100 % are 1.00, 0.77, 0.70, and 0.57, respectively. The EF at 1, 5, 10, and 100 % are 24.3, 12.6, 6.2, and 3.3, respectively. AUC and EF values are good and acceptable. Therefore, virtual screening results are valid and the protocol was successful in retrieving hit compounds in alignment with the model.

Physicochemical properties of chemical compounds influence their biological activity. Molecular weight, lipophilicity, rotational bonds, and number of hydrogen bond donors or acceptors are important molecular descriptors that significantly affect ionization, absorption, distribution, metabolism, excretion, toxicity, and binding affinity properties of drugs. Lipinski’s rule of 5 is a key to deciding the oral bioavailability of drugs and an MW ˃ 500 may result in poor oral absorption (Bunally et al., 2019). Accordingly, to downsize the number of hits, we filtered the initial 66-hit library through a molecular weight filter of ≤ 500, which yielded 50 hits. When we attempted to apply a rotational bond filter and ADME filter, the number of hits was reduced but the results obtained in sequential molecular docking were not fruitful. Hence, we proceeded with 50 hit compounds having a molecular weight of less than 500. Table S3 contains information about 50 hit compounds and their MW. Fig. 4 shows 50 hits of MW less than 500, aligned to the pharmacophore model.

Close

Fig. 4

Virtual screening results show 50 hit compounds aligned to the pharmacophore model.

Export to PPT

3.2.1 Benchmarking molecular docking

The hit library for molecular docking, generated by pharmacophore-based virtual screening contained marine natural compounds that are relatively small, as it is a library that has been filtered based on MW. The decoys from DEKOIS 2.0 were used for benchmarking. Therefore, 50 hits from virtual screening and 720 decoys were docked to 7LBS using both docking engines. The performance of AutoDock and AutoDock Vina were compared for their predictive ability and their accuracy to select active compounds at top ranks. Conventional evaluation metrices of individual and consensus scoring were applied for performance evaluation (Empereur-Mot et al., 2015; Empereur-Mot et al., 2016).

3.2.1.1

3.2.1.1 Individual scoring

The binding energy scores obtained from both docking engines were almost similar in the range of −9.9 to −6.4 kcal/mol for active compounds. Fig. S1 shows the distribution of scores for all actives and decoys and their corresponding ranks allotted by the two docking engines. Fig. 5 shows the predictiveness curve, ROC, and enrichment curves obtained from individual scoring. The predictiveness curve provides the probability of activity [P(act)] based on the scores and corresponding ranks of actives. The P(act) for AutoDock was 0.573 and AutoDock Vina was 0.715 at a set threshold of 0.01. Both docking engines identified the native ligand XR8-24 at the first position in the active compound list. None of the hit compounds had a score greater than the native ligand. Table 2, is the top 1 % rank list retrieved from two docking engines, has 6 actives each, out of which 3 compounds were enriched, but with different ranks. Table 3 shows the structures of 3 compounds that enrich the top 1 % rank list from AutoDock and AutoDock Vina. From Fig. 5 it is understood that both docking engines retrieved active compounds with early enrichment and possess good predictive power, AutoDock Vina scored a better position though. The metrices calculated from Fig. 5A, 5B, and 5C are reported in Table 4 under individual scoring. To reach an agreement between the results and to improve active compound selection, we attempted consensus scoring.

Close

Fig. 5

Comparative molecular docking results based on individual scoring. Red. AutoDock, green. AutoDock Vina. A) Predictiveness curves B) ROC curves C) Enrichment curves.

Export to PPT

Table 2 Top 1% compounds retrieved by comparative molecular docking.

AutoDock		AutoDock Vina
Compound retrieved	Predicted binding energy (ΔG, kcal/mol)	Compound retrieved	Predicted binding energy (ΔG, kcal/mol)
XR8-24	−9.8	XR8-24	−9.9
CMNPD23373	−9.8	CMNPD28766	−9.9
CMNPD30040	−9.8	CMNPD13919	−9.8
CMNPD28766	−9.7	CMNPD7007	−9.8
CMNPD13919	−9.6	CMNPD24799	−9.7
CMNPD1441	−9.6	CMNPD23373	−9.7

Table 3 Structures of compounds that enrich the top 1% rank list retrieved from AutoDock and AutoDock Vina.

Table 4 Individual and consensus scoring metrices of AutoDock and AutoDock Vina.

Individual scoring - global metrices					Individual scoring- partial metrices
Docking engine	ROC-AUC	TG	RIE	BEDROC	EF1%	EF5%	EF10%
AD	0.957	0.680	6.99	0.630	11.34	7.56	7.36
AV	0.993	0.918	9.98	0.901	13.23	13.13	9.82
Consensus scoring - global metrices					Consensus scoring - partial metrices
Metrices	ROC-AUC		RIE	BEDROC	EF1%	EF5%	EF10%
Min of ranks	0.974		7.74	0.699	5.67	8.75	8.64
Average of ranks	0.991		9.93	0.896	11.34	11.34	9.82
Max of Z scores	0.995		10.32	0.931	13.53	13.13	9.82
Average of Z scores	0.981		9.02	0.814	11.34	11.14	8.83
Max of normalized scores	0.974		8.08	0.729	9.45	9.15	8.85
Average of normalized scores	0.992		9.91	0.894	11.34	12.73	9.82

3.2.1.2

3.2.1.2 Consensus scoring

Consensus scoring ensures fairness between scores of two docking engines and reduces the discrepancies due to outliers. Screening Explorer applies i) normalization of scores ii) standardization of scores, and iii) ranking-dependent methods for consensus scoring. For both docking engines, all 3 of the consensus scoring were carried out. A consensus selection of the best active compound arising from different scoring functions is possible based on the maximum of Z and normalized scores, and the minimum of ranks. The dominance or shortfall of scoring functions, if any, are balanced by implementing averages of ranks, Z, and score values. Fig. 6 shows the ROC and enrichment curves from consensus scoring.

Close

Fig. 6

Consensus scoring of AutoDock and AutoDock Vina results. A) ROC curves, B) Enrichment curves, blue (minimum of ranks), red (average of ranks), purple (maximum of z scores), green (average of Z scores), yellow (maximum of normalized scores), and cyan (average of normalized scores).

Export to PPT

3.2.1.4

3.2.1.4 Active compound properties and pharmacophore model alignment

CMNPD28766, is aspergillipeptide F. SwissADME predicted molecular properties of aspergillipeptide F: MW – 450.53 g/mol, HBA – 7, HBD – 6, rotatable bonds – 14, lipophilicity (Log P) – 2.41, pharmacokinetic properties: GI absorption low, BB permeation low, the substrate of p-glycoprotein, not an inhibitor of CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4, bioavailability score 0.55, and drug-likeness properties: Lipinski’s violation – 1 (HBA, HBD > 5), lead likeness violation – 2 (MW > 350, rotatable bonds > 7), synthetic accessibility – 3.98 (easy). No PAIN or structural alert from ADME. Collectively, predicted molecular and ADME properties indicate that aspergillipeptide F has good drug-likeness (up to 1 Lipinski’s violation is acceptable) and is less toxic (Daina et al., 2017). Chemically aspergillipeptide F is described in the literature as a linear tetrapeptide (Carroll et al., 2019; Liu, 2017). It was reported from the fungi Aspergillus sp. of the South China sea. No biological activities have been reported for aspergillipeptide F. CMNPD28766 pharmacophore-fit score is 75.916 in LigandScout pharmacophore model-aided virtual screening. Interestingly, this compound topped the pharmacophore-fit score list of compounds possessing MW less than 500 g/mol. CMNPD28766 had 7 out of 9 pharmacophoric features matched with the pharmacophore model. Fig. 7A and 7B show the pharmacophore features alignment of aspergillipeptide F (CMNPD28766) with the 3D and 2D pharmacophore models respectively. Aspergillipeptide F possesses a positive ionization centre, 2 HBA, 1 HBD, and 3 hydrophobic features aligned to the structure-based pharmacophore model. When compared to the developed pharmacophore model, 1 HBD and 1 hydrophobic feature are absent in the aspergillipeptide F derived model.

Close

Fig. 7

A) CMNPD28766, aspergillipeptide F aligned to 3D pharmacophore model B) 2D depiction of aspergillipeptide F pharmacophore model interaction features.

Export to PPT

3.2.1.5

3.2.1.5 Binding interaction analysis

Native ligand XR8-24 was included along with hit compounds for docking and its active conformation generated by the developed docking protocol was aligned to the native PDB conformer. Fig. 8A shows the superposed native and redocked conformers of XR8-24 that exhibited an RMSD of 0.764 Å. The active conformations of native ligand XR8-24 and aspergillipeptide F in the active site of PLpro are shown in 3D in Fig. 8B. The RMSD between these two conformers was 1.056 Å indicating an acceptable fit of aspergillipeptide F into the same binding site of XR8-24. Fig. 8C is the 3D illustration of the binding sites occupied by aspergillipeptide F. Fig. 9A and 9B represent the 2D interaction diagrams of XR8-24 and aspergillipeptide, respectively. The SARS-CoV-2 PLpro catalytic binding site is a triad of Cys111, His272, and Asp286 residues. GRL0617 is the prototype of SARS-CoV-2 PLpro non-covalent inhibitors and is a peptide (Osipiuk et al., 2021; Shen et al., 2022). Shen et al. (2022) aimed to develop potent PLpro inhibitors using GRL0617 as the lead and demonstrated that ligand XR8-24 undergoes non-covalent binding to multiple binding sites implementing binding cooperativity leading to enhanced in-vitro anti-viral potency than GRL0617. Shen et al. (2022) concluded that there are five binding sites far from the active catalytic triad residues of PLpro that needs to be occupied for potent inhibition. Site I, Glu167, and Asp164 are crucial for PLpro binding to its substrate ubiquitin/ISG15. Site II is surrounded by Arg166 which interacts with ubiquitin. The compounds designed by this research group failed to interact with site II, except Zn-3–56 which lacked the characteristic azetidine of other compounds and ZN-3–56 also exhibited improved in-vitro activity than GRL0617. Site III, the BL2 loop featuring Gln269 and Tyr268 also plays a significant role in substrate binding. Site IV is hydrophobic due to Leu162, Tyr264, and Tyr273. In addition, a new binding site V, named BL2 groove formed by the hydrophobic residues Pro248 and Pro299 was identified. Non-covalent inhibitors are capable of PLpro inhibition by engaging site V without accessing the active catalytic triad site. Engaging site V blocks the access of the catalytic site of PLpro to its substrate ubiquitin/ISG.

Close

Fig. 8

A) Superposed native and redocked conformations of XR8-24, red (native), green (redocked), B) Superposed active conformations of XR8-24 (stick) and aspergillipeptide F (ball and stick) inside the binding pocket of PLpro, C) Aspergillipeptide F conformation engaging all five binding sites of PLpro.

Export to PPT

Close

Fig. 9

2D interaction diagrams of A) XR8-24B) aspergillipeptide F with PLpro active site residues.

Export to PPT

CMNPD23373 and CMNPD13919 were unable to bind to all 5 binding sites. We compared the binding interactions of XR8-24 and aspergillipeptide F as shown in Fig. 9A and 9B. XR8-24 has 5 rings and interacts in a double ‘U’ conformation and aspergillipeptide F interacts in a single ‘U’ shaped conformation as it has only one ring. The most significant observation is that the LigandScout-predicted 7LBS-XR8-24-based model in Fig. 2, proposed electrostatic interactions through the positive ionization centre and HBD in the azetidine ring of XR8-24 with Glu167 of site I, but the active conformation of XR8-24 lacked electrostatic interactions and exhibited van der Waal’s interactions. Similarly, aspergillipeptide F also interacted with Glu167 through van der Waal’s interactions. The fascinating phenomenon was that aspergillipeptide F behaved according to the predicted pharmacophore model and the ionized primary amine in the side chain established an ionic interaction with Asp164 and a hydrogen bond with Ala246. Ionic hydrogen bonds play an important role in biomolecule recognition by receptors; hence it is suggested that the pharmacophore model was effective in predicting a suitable pharmacophore for molecular recognition by SARS-CoV-2 PLpro (Meot-Ner, 2005). Asp164 of site I binds through a Pi-anion bond to XR8-24.

XR8-24 established 3 hydrogen bonds whereas aspergillipeptide formed 5 hydrogen bonds (4 HBA and 1 HBD atoms). The developed pharmacophore model in Fig. 2D proposed four hydrogen bonds taking into account 2 HBA and 2 HBD groups of ligands. Interestingly, aspergillipeptide F interacted through its HBA ‘O’ atoms of the amide chain to form 2 hydrogen bonds with Arg166 of site II. Designing compounds that can block site I and site II simultaneously was reported to be difficult and harnessing natural compounds proved effective (Shen et al., 2022). XR8-24 lacks interaction with Arg166 which suggests the significance of the structural motif of aspergillipeptide F. XR8-24 and aspergillipeptide F (HBA ‘O’ atom of –COOH), both formed a hydrogen bond with Site III residue Gln269. Site IV Tyr264 interacted by forming a hydrogen bond with XR8-24 whereas Tyr273 of site IV interacted through a hydrogen bond with aspergillipeptide F (HBA ‘O’ atom of –COOH).

The predicted pharmacophore model in Fig. 2D proposed 4 hydrophobic interactions and in Fig. 7B, the model for aspergillipeptide F predicted 3 hydrophobic interactions. The chemical features of active poses of XR8-24 and aspergillipeptide F involved in binding are compared in Table 5. Molecular docking has predicted 2 hydrophobic interactions, besides van der Waal’s interactions. As per the developed model, methyl and 2 isopropyl groups of aspergillipeptide F were predicted to involve in hydrophobic interactions. Molecular docking also predicted the contributions of methyl and 1 isopropyl group of aspergillipeptide F towards hydrophobic interactions. However, the phenyl ring of aspergillipeptide F contributes significantly to hydrophobic interactions. Hydrophobic interactions were almost similar in XR8-24 and aspergillipeptide F, the most significant interaction being the Pi-Pi T-shaped interaction with Tyr268. Interactions of XR8-24 and aspergillipeptide F with Gln269 and Tyr268 lead to BL2 loop closure. Also, van der Waal’s interaction with Asn267 supports BL2 loop closure. The closure of the BL2 loop blocks the access of substrates to the active site, thus blocking the PLpro protease activity resulting in inhibition of virus replication. Pi-alkyl interaction of aspergillipeptide F with Pro248 of site V closes the BL2 groove, blocks the access of PLpro active site to ubiquitin/ISG15, and inhibits the deubiquitinating activity of PLpro, thus improving host anti-viral immunity.

Table 5 Comparison of chemical features of the active binding poses of XR8-24 and Aspergillipeptide F.

Nature of interactions	Chemical features
	XR8-24	Aspergillipeptide F
Hydrogen bond	Oxygen (acceptor) and nitrogen (donor) of the amide	Carboxylic acid oxygens (acceptor), keto oxygens (acceptor)of peptide, and nitrogen (donor) of amine
Electrostatic π-anion	Phenyl	Absent
π-sulfur	Sulfur of thiophene	Absent
Van der Waal’s	Azetidine, phenyl, pyrrole, thiophene, and aliphatic methyl	Phenyl, aliphatic isopropyl, and methyl
Hydrophobic	Phenyl, pyrrole, and methyl	Phenyl
Attractive charge (ionic)	Absent	Protonated amine

The interatomic distances for the charged interactions of aspergillipeptide F with PLpro range between 1.82 and 2.96 Å and hydrophobic interactions occurred at a distance of 2.44–4.46 Å suggesting stable interactions. Amide or peptide linkages have been associated with greater PLpro binding than other moieties (De Vita et al., 2020; Osipiuk et al., 2021; Pang et al., 2021; Shen et al., 2022). Molecular docking results are in favour of peptide linkages in aspergillipeptide F. Together, pharmacophore model-aided virtual screening along with consensus molecular docking was successful in retrieving a natural compound of marine origin with the potential of SARS-CoV-2 PLpro inhibition similar to synthetic potent non-covalent inhibitor XR8-24.

Thus, explored in-silico binding interactions delineate the mode of action of aspergillipeptide F as an anti-SARS-CoV-2 compound. Consequently, we proceeded to verify the stability of the aspergillipeptide F-PLpro complex and essential dynamics of PLpro during aspergillipeptide F binding by molecular dynamics simulation studies.