Molecular docking of potential SARS-CoV-2 papain-like protease inhibitors.

SARS-CoV-2 papain-like protease is considered as an important potential target for anti-SARS-CoV-2 drug discovery due to its crucial roles in viral spread and innate immunity. Here, we have utilized an in silico molecular docking approach to identify the possible inhibitors of the SARS-CoV-2 papain-like protease, by screening 21 antiviral, antifungal and anticancer compounds. Among them, Neobavaisoflavone has the highest binding energy for SARS-CoV-2 papain-like protease. These molecules could bind near the SARS-CoV-2 papain-like protease crucial catalytic triad, ubiquitination and ISGylation residues: Trp106, Asn109, Cys111, Met208, Lys232, Pro247, Tyr268, Gln269, His272, Asp286 and Thr301. Because blocking the papain-like protease is an important strategy in fighting against viruses, these compounds might be promising candidates for therapeutic intervention against COVID-19.

Introduction

The coronavirus disease 2019 (COVID-19) pandemic caused by the new severe acute respiratory syndrome coronavirus type 2 (SARS-CoV-2) poses an unprecedented challenge [1,2]. Globally, as of October 21, 2020, more than 40.6 million cases have been reported to the World Health Organization, and 1,121,843 deaths have occurred in more than 235 countries. The scale of the pandemic and the rapidly spread make COVID-19 one of the most serious diseases facing mankind so far.

At present, there is no vaccine or approved antiviral treatment for SARS-CoV-2. In the extensive research on the development of new antiviral drugs for COVID-19, main protease of SARS-CoV-2, which plays an important role in blocking viral infection, has been identified as a potential target [3]. Another attractive antiviral target is papain-like protease (PLpro) non-structural protein 3 (nsP3), an essential component of the replicase-transcriptase complex [4]. Its topological structure is divided into ubiquitin-like (UBL), thumb, palm, and fingers [5]. It process the viral polyprotein by recognizing the tetrapeptide Leu-X-Gly-Gly motif, and functions as de-ubiquitinating and de-ISGylating enzyme to suppresses inflammation and antiviral signaling [6,7]. Therefore, inhibition of PLpro activity can prevent virus replication and destroy its role in host immune response evasion, making it a promising target for antiviral drugs. Recently, an inhibitor for PLpro GRL-0617 has been proved to be an effective therapeutic agent for SARS-CoV PLpro [8,9].

In current study, the inhibitory effects towards PLpro for 21 compounds, including azvudine, aminoethyl, bavachinin, corylifol A, chromen, darunavir, disulfiram, ebselen, efavirenz, famciclovir, GRL-0617, isobavachalcone, mercaptopurine, neobavaisoflavone, oseltamivir, papyriflavonol A, psoralidin, ribavirin, sofosbuvir, tioguanine and 4′-O-methylbavachalcone were investigated by molecular docking. These compounds are antiviral, antifungal or anticancer drugs. They have been proved to be highly effective against HIV, HCV, influenza virus, coronavirus, antifungal and anticancer agents [[10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21]]. The focus of our study is on the repurposing and identification of compounds, aiming to accelerate the identification of potential drugs for COVID-19 treatment.

Materials and methods

Optimization of the 3D coordinates of the ligands

The 3D coordinates of the 21 drug ligands, azvudine (CSID: 24717759), aminoethyl (CSID: 29361441), bavachinin (CSID: 8512670), corylifol A (CSID: 22913562), chromen (CSID: 57578897), darunavir (CSID: 184733), disulfiram (CSID: 3005), ebselen (CSID: 3082), efavirenz (CSID: 57715), famciclovir (CSID: 3207), GRL-0617 (Chemical Book ID: CB74667969), isobavachalcone (CSID: 4444667), mercaptopurine (CSID: 580869), neobavaisoflavone (CSID: 4478223), oseltamivir (CSID: 58540), papyriflavonol A (CSID: 8518529), psoralidin (CSID: 4445118), ribavirin (CSID: 34439), sofosbuvir (CSID: 26286922), tioguanine (CSID: 2005804), and 4′-O-methylbavachalcone (CSID: 24845874) were obtained from ChemSpider (http://www.chemspider.com/) and Chemical Book (https://www.chemicalbook.com/) as a .mol file. The code within the parenthesis signifies the ChemSpider ID and ChemicalBook Number of the respective molecules. The secondary structure shown is the predicted by DSSP for SARS-CoV-2 (PDB: 7JIW), Rhinolophus affinis (GenBank: KF569996.1), SARS (PDB: 4OW0), MERS (PDB: 4RF1) and MHV (PDB: 5WFI) PLpro. Similarity and alignment calculations were performed using MEGA-X and ClustalW.

Drug-like properties of the ligands

The drug likeliness of a molecule is indicated by the Lipinski’s rule of five parameters (molecular weight ranges from 152.18 Da to 547.66 Da, no more than 5 hydrogen bond donors, no. of hydrogen bond acceptors should be less than 10 and logP should not be greater than 5). The Lipinski’s rule of five parameters was obtained from the SWISSADME server (www.swissadme.ch/index.php). The chemical structures, chemical formula and the Lipinski’s rule parameters of the ligands are listed in Table 1 .

Table 1

Candidate drugs properties screened and dockscore through LibDock in the study.

No.	Compounds	Structure	Treatment	Mass	Dockscore
1	Azvudine	Image 1	HIV	286.22	83.52
2	Aminoethyl	Image 2	HIV	312.41	65.25
3	Bavachinin	Image 3	SARS	338.40	82.65
4	Corylifol A	Image 4	SARS	390.47	89.15
5	Chromen	Image 5	H1N1, H3N2	390.47	88.00
6	Darunavir	Image 6	HIV	547.66	99.38
7	Disulfiram	Image 7	HIV, MERS, SARS	296.54	36.89
8	Ebselen	Image 8	HIV	274.18	69.10
9	Efavirenz	Image 9	HIV	315.68	69.73
10	Famciclovir	Image 10	VZV	321.33	75.77
11	GRL-0617	Image 11	SARS	304.39	67.70
12	Isobavachalcone	Image 12	SARS	324.37	52.81
13	Mercaptopurine	Image 13	SARS	152.18	60.27
14	Neobavaisoflavone	Image 14	SARS	322.35	110.53
15	Oseltamivir	Image 15	AIV, H1N1, IBV	312.41	94.09
16	Papyriflavonol A	Image 16	Antifungal	438.47	56.56
17	Psoralidin	Image 17	Anticancer	336.34	120.99
18	Ribavirin	Image 18	HCV	244.21	71.04
19	Sofosbuvir	Image 19	HCV	529.45	124.29
20	Tioguanine	Image 20	Anticancer	167.19	37.43
21	4′-O-methylbavachalcone	Image 21	SARS	352.42	82.85

Structure and preparation of PLpro for docking studies

SARS-CoV-2 papain-like protease (PDB: 6W9C, chains A, B, C) was used as the template for molecular docking. The 6W9C is solved by cryo-electron microscopy with 2.70 Å resolution, which is a trimer of three homologous amino acid chains (A, B, and C). Water molecules and ligands were removed from the PBD using Pymol.

Molecular docking analyses and visualization

Blind dockings were conducted based on docking ligands into a protein target using Discovery Studio (BIOVIA, 2016). In the key steps of the dockings, the protein remained rigid, while the ligand was fully flexible. Energy minimization of the protein was performed using the DS with CHARMM force field to relax the structure and remove the steric overlap. Then, polar hydrogen atoms were added. The centroid of the site model was used to define the active site box required by LibDock. The SARS-CoV-2 PLpro-binding site in the X-ray crystal structure of ligand complex was determined by co-crystallization (X: -46.0103, Y: 14.2997, and Z: 29.9464), with a radius of 13.8 Å. The ligand conformation generation method was selected as “best”, the number of binding site hotspots was set to 100, and other parameters were set by default. After each compound was docked, the best conformations and LibDock score were obtained. Considering their interaction energy and binding free energy, the compounds were screened. The binding modes (3D or 2D ligand-receptor interaction simulation map) of the selected compounds were analyzed.

Results and disscussion

Docking and scoring software is used widely to enhance the drug design and predict the interaction between drugs and macromolecules in the pharmaceutical products. The docking method involves searching for binding sites on the whole macromolecular surface. Therefore, the papain-like protease of COVID-19, SARS-CoV-2 PLpro was used for blind docking analysis of some known drugs and bioactive substances. Fig. 1 A shows the overall workflow of Discovery Studio (DS) LibDock algorithm. Then SARS-CoV-2 PLpro structure was submitted to the SARS-CoV-2 PLpro structure hotspots methods to predict the protein-ligand binding sites.

Fig. 1

Workflow of molecular docking, architecture of SARS-CoV-2 PL, and sequences alignment of PLfrom coronaviruses. (A) A schematic workflow of protein-ligand complex formation by Discovery Studio. (B) Solvent-accessible surface representation of SARS-CoV PLpro is shown in blue. The two ubiquitin, ubiquitin-like binding subsites of SARS-CoV-2 PLpro and catalytic domain are shown in the solvent accessible surface representation with SUb1 shaded cyans, SUb2 shaded golden and red circle, respectively. (C) The PLpro from SARS-CoV-2 (7JIW), SARS (4OW0), MERS (4RF1), Rhinolophus affinis (KF569996.1) and MHV (5WFI). The secondary structure shown is the predicted by DSSP for SARS-CoV PLpro (5E6J). Similarity and alignment calculations were performed using ClustalW. Residue positions form the catalytic triads are marked with red triangle, while key residues forming ubiquitination and ISGylation binding sites are marked with blue stars. T, turn. The blocking loop2 (BL2) is boxed in gold. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

In addition to the recommend drugs for COVID-19, we have also conducted molecular docking studies with SARS-CoV-2 PLpro for various natural compounds, antiviral drugs, antifungal drugs and anticancer drug (Table 1). Among SARS-CoV-2 PLpro, SARS PLpro, MERS PLpro, Rhinolophus affinis PLpro, and MHV PLpro, classic catalytic triad, two different ubiquitin/ubiquitin-like conjugates and ISG15 modified protein binding sites are quite conservative (Fig. 1B). And some of them may be related to protein stability, kinetics, ligand binding and catalytic performance [9,[22], [23], [24]]. Amino acid sequence alignment showed that SARS-CoV-2 and SARS PLpro (82.86%) as well as Rhinolophus affinis PLpro (83.49%) were highly similar, and there were significant differences between SARS-CoV-2 and MERS (29.84%) and MHV PLpro (29.52%) (Fig. 1C). Therefore, SARS-CoV-2 PLpro is more closely related to SARS and Rhinolophus affinis than MERS and MHV. Examples include the active conservative catalytic trinity of Cys111, His272 and Asp286. The catalytic Trinity accepts Cys111 in the form of –SH as nucleophilic agent, His272 as general acid and alkali, and Asp286 promotes the deprotonation of Cys111 [25,26].

Aminoethyl, corylifol A, chromen, darunavir, ebselen, GRL-0617, isobavachalcone, papyriflavonol A, sofosbuvir and 4′-O-methylbavachalcone shows the formed interactions among the drug ligands near the catalytic triad. Aminoethyl and darunavir form conventional hydrogen bonds with Asp286 (Fig. 2 A and B). In addition, aminoethyl constructs ten van der waals with Leu162, Gly163, Val165, Asn267, Tyr268, His272, Gly287, Ala288, Leu289 and Thr301, two π-alkyl interactions with Lys105 and Pro248, two π-sigma with Asp164 and Tyr264, and two unfavorable donor-donor interactions of Trp106 and Tyr273. Darunavir constructs two carbon hydrogen bond with Asn267 and Tyr268, eighteen van der waals with Ala107, Asp108, Asn109, Lys157, Gly163, Asp164, Arg166, Glu167, Met208, Ala246, Pro247, Tyr264, Gly266, Gln269, Tyr273, Gly287, Leu289 and Thr301, four π-alkyl interactions with Lys105, Trp106, Leu162 and Pro248. However, the donor-donor interaction between Darunavir and Ala288 is unfavorable. Darunavir may block the binding of PLpro with target protein and exert its enzymatic activity. Chromen and corylifol A has almost the same site binding pattern, forming two conventional hydrogen bonds with Lys105 and Ala288 (Fig. 2C and D), a carbon hydrogen bond with Tyr268, ten van der waals with Lys92, Lys94, Leu162, Gly163, Asp164, Arg166, Glu167, Tyr264, Asn267 and Gly287, three π-alkyl interactions with Trp93, Ala107 and Leu289, a π-sigma with Trp106, and a unfavorable acceptor-acceptor sites with Asp286. Ebselen is an organic selenium compound that mimics the activity of glutathione. This compound forms a covalent bond between the selenium atom and thiols in the cysteine residue by forming a selenylsulfide (-Se-S-) linkage to inhibit HIV-1 capsid dimerization [27]. In the paper, we found that ebselen constructed a conventional hydrogen bonds with Lys105 (Fig. 2E), ten van der waals with Trp106, Ala107, Asp164, Glu167, Asn267, Tyr268, Tyr273, Asp286, Ala288 and Leu289. Although ebselen could not form -Se-S- with Cys111, it interacts with the PLpro active site. It GRL-0617 interacts with the residues in the active site through hydrogen bonding with Lys105 and Trp106, π-alkyl interaction with Pro247, Pro248, Ala288 and Leu289, van der waals forces of attraction with other residues (Fig. 2F). GRL-0617 is contained in the substrate cleft formed between the β11-β12 strands and the loop connecting α3 and α4 (α3-to-α4 loop), which occupies the S3–S4 pockets [4,6]. It leads to the loss of PLpro activity. Isobavachalcone and papyriflavonol A form van der waals forces with Asp286, π-alkyl hydrophobic interaction with Cys111, His272 and other sites (Fig. 2G and H). Sofosbuvir and 4′-O-methylbavachalcone construct four and two conventional hydrogen bonds, respectively (Fig. 2I and J). They also form van der waals with Asp286 and with other sites.

Fig. 2

The minimum docked poses of the compounds along with their corresponding catalytic triad interaction plots within the active site of SARS-CoV-2 PLpro.

The PLpro can also inhibit host innate immune response by reversing host ubiquitination and ISGylation events [8]. In addition to blocking catalytic sites, PLpro binding, processing active sites of ubiquitin and ISG15 protein sites are highly desirable [22]. Azvudine, bavachinin, disulfiram, efavirenz, famciclovir, mercaptopurine, neobavaisoflavone, oseltamivir, psoralidin and ribavirin show interactions proximity to the sites, locating around Trp106, Asn109, Met208, Lys232, Pro247, Tyr268, Gln269, Thr301 (Fig. 3 A–J). Tioguanine (6-thioguanine, 6 TG) has been proven to be effective in inhibition the enzyme activity of SARS PLpro by forming of a hydrogen bond with Cys112 [19,28]. Here, it is found that 6 TG forms two conventional hydrogen bonds with Tyr35 and Thr54, a carbon hydrogen bond with Asn146, six van der waals with Asn13, Leu36, Gly38, Phe55, Tyr83 and Met84, π-alkyl with Leu87, π-π T-shaped interaction forces with Tyr56 and other forces as depicted in the Fig. 3K. It may affect the host ubiquitination and ISGylation events through the bond interaction with the N-terminal UBL domain.

Fig. 3

The minimum docked poses of the compounds along with their corresponding ubiquitination and ISGylation interaction plots within the active site of SARS-CoV-2 PLpro.

The results of blind molecular docking of 21 different classes of compounds showed that they might be potential anti-SARS-CoV-2 drugs because they bind to the active site of the papain-like protease of SARS-CoV-2 PLpro. Inhibition of this enzyme/protease is the key to blocking virus replication and reproduction. The inhibitor binding sites in SARS-CoV-2 protease are composed of key residues Cys111, His272 and Asn286 as evident from the recent study on the GRL-0617 and ebselen inhibitors for SARS-CoV-2 PLpro [4,29]. Through the blind docking study of all 21 molecules with SARS-CoV-2 protease, we found that these molecules are usually surrounded by the above mentioned residues, ubiquitination and ISGylation sites, which clearly indicates that those molecules can prevent the replication of SARS-CoV-2 virus and evade the host immune response. The changes areas of the important interacting residues in the ligands and protease active sites were listed in Supplementary Table 1. Although blind docking studies have shown that all molecules can be used as potential drugs for the treatment of COVID-19, based on the estimated binding energy (ΔG) value, we can infer that among all the compounds studied, neobavaisoflavone (with the highest negative minimum ΔG value (−459 kcal/mol) may be the best potential inhibitor of SARS-CoV-2, followed by oseltamivir (−121.55 kcal/mol) ˃ sofosbuvir (−119.44 kcal/mol) ˃ famciclovir (−85.61 kcal/mol) ˃ isobavachalcone (−84.75 kcal/mol) ˃tioguanine (−78.64 kcal/mol) ˃ chromen (−67.92 kcal/mol) ˃ efavirenz (−66.98 kcal/mol) ˃ papyriflavonol A (−51.99 kcal/mol) ˃ ebselen (−50.99 kcal/mol) > corylifol A (−46.78 kcal/mol) ˃ mercaptopurine (−46.61 kcal/mol) ˃ 4′-O-methylbavachalcone (−42.64 kcal/mol) ˃ azvudine (−33.1 kcal/mol) ˃ psoralidin (−29.89 kcal/mol) ˃ ribavirin (−26.49 kcal/mol) ˃ bavachinin (−25.59 kcal/mol) ˃ disulfiram (−24.84 kcal/mol) ˃ GRL-0617 (−24.62 kcal/mol) ˃ aminoethyl (−20.87 kcal/mol) ˃ darunavir (−8.74 kcal/mol) (Fig. 4 ).

Fig. 4

The binding energy of the interaction between SARS-CoV-2 PLpro and 21 compounds.

In this study, the free energy of binding energy of compounds and PLpro was ranged from −459 kcal/mol to −8.74 kcal/mol, indicating that those compounds could interact with the SARS-CoV-2 PLpro. Our finding also suggests the potential binding modes with known therapeutic drugs for COVID-19. Since this study was performed using computational methods, special attention should be paid to these results and further experiments in vivo and in vitro are needed.

Funding

This work was supported by grants from National Key Plan for Research and Development of China [2016YFD0500300], [82072270 and 81871663], National Major S & T Project for the Prevention and Treatment of Major Infectious Diseases in China [2017ZX10004206-007], the Innovation Project of Shandong Academy of Medical Sciences, and Academic Promotion Programme of Shandong First Medical University [2019LJ001].

Declaration of competing interest

We declare no competing interests.