Repurposing Drugs as Potential Therapeutics for the SARS-Cov-2 Viral Infection: Automatizing a Blind Molecular Docking High-throughput Pipeline

Abstract

In the context of the COVID-19 pandemic, scientists worldwide have been looking for ways to stop it using different approaches. One strategy is to look among drugs that have already proved safe for use in humans and tested for other illnesses. Several components from the virus and the infected cell are the potential therapeutic targets from a molecular perspective. We explain how we implemented a cavity-guided blind molecular docking algorithm into a high-throughput computational pipeline to automatically screen and analyze a large set of drugs over a group of SARS-CoV-2 and cell proteins involved in the infection process. We discuss the need to significantly extend the conformational space sampling to find an accurate target-ligand complex. Our results identify nine drugs with potential multi-target activity against COVID-19 at different stages of the infection and immune system evasion. These results are relevant in understanding the SARS-CoV-2 drug’s molecular mechanisms and further clinical treatment development. The code developed is available on GitHub [https://github.com/tripplab/HTVS].

1. Introduction

The coronavirus disease-2019 (COVID-19) is the third documented viral outbreak caused by a member of the Coronaviridae family. From 2002 to 2004, the severe acute respiratory syndrome coronavirus (SARS-CoV) spread to 29 countries, causing 8422 confirmed cases and 916 deaths, and is considered the first emerging epidemic of the twenty-first century [1, 2]. Later in 2012, the middle-east respiratory syndrome coronavirus (MERS-CoV) caused 2585 confirmed cases and 890 deaths to date [3]. In less than two decades since the appearance of SARS-CoV, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) emerged in late 2019 and has spread worldwide ever since by human-to-human transmission. As of April 29, 2022, there are more than 510 million confirmed cases and 6.2 million deaths related to SARS-CoV-2 infection, and it continues to increase at present [4]. The Coronaviridae family comprises a group of enveloped crown-shaped single-stranded positive-sensed RNA viruses (ssRNA+) with multiple domestic and wild animal reservoirs [5]. Lessons from previous and current outbreaks have shown the severity of cross-species transmission, which has led to concerns about health emergencies, such as COVID-19. The transmission of this disease occurs through an infected person’s respiratory droplets carrying the SARS-CoV-2, and the severity ranges from asymptomatic cases, mild and moderate flu-like symptoms, to critical illness requiring intensive care with mechanical ventilation, and death [6]. Global contributions and efforts following the COVID-19 outbreak have unraveled a considerable amount of information about viral infection, transmission, infection cycles, and immune evasion. Currently, the three-dimensional proteome structures of the SARS-CoV-2 are available on the RCSB protein data bank [7]. Therefore, it is feasible to evaluate drug-like small molecules against relevant targets in the viral infection cycle through a structure-based molecular docking approach. Blind molecular docking, unlike traditional molecular docking, does not require prior knowledge of target binding sites, which simplifies the automatizing of the process since it only needs the structural information of the target. In the past, this process was considered less accurate than the traditional. However, methods, such as CB-dock, have overcome this limitation by reducing the nonrelevant conformation sampling by directing the molecular docking on putative sites instead of the whole protein structure [8]. The integration of this tool into our customized high-throughput virtual screening pipeline allows the screening of N sorted-by-size cavities. The cavity-based search is an exciting scenario because protein-ligand interactions usually occur in large protein cavities or pockets that frequently contain the active site [9]. Moreover, the exploration of cavities in the vicinity of protein-protein interfaces (PPI) is also an attractive approach to searching for effective inhibitors since it plays an essential role in nearly all biological processes, including SARS-CoV-2 infection [10, 11]. In this context, screening already-known drugs with described pharmacology, dose, toxicity, formulation, and proven to be safe for use in humans represents a low-risk and cost-effective strategy to considerably shorten the time required for drug approval [12, 13]. We present an in-house customizable pipeline that integrates a cavity-guided blind molecular docking algorithm to extend the conformation space sampling on putative sites significantly. We also report the methodology to follow and results of the virtual screening of 47 drugs for potential repurposing against 16 structures of 10 viral and cell targets that are key in the SARS-CoV-2 infection cycle.

2. Overview of the SARS-CoV-2 infection cycle

The initial stage of the infection cycle starts with the recognition and anchoring of the SARS-CoV-2 spike protein complex into the host angiotensin-converting enzyme 2 (ACE2) through the receptor-binding domain (RBD) located at each one of the 3S proteins [14]. Then, the activation of the spike occurs at the surface or endosome level by transmembrane serine protease 2 (TMPRSS2) or cathepsin B/L proteases, respectively, to allow viral entry [15]. Once the virus membrane merges with the cell membrane, the genomic material enters the cell. The cell’s ribosomes then translate the viral RNA into pp1a/ab polyproteins, which will be later processed by cleavage through the enzymatic activity of the main protease (Mpro) and the papain-like protease (PLpro) [16]. This process will release 16 non-structural proteins (NSPs), including the RNA-dependent RNA polymerase (NSP12) and co-factors NSP7, and NSP8 of the RNA-replication machinery (Rdrp). After replication, expression of the structural proteins occurs, the genomic material is packaged, and the virion is assembled on a lipid membrane and matured for subsequent exocytosis. In addition, evidence suggests that the SARS-CoV-2 proteases and some of their cleavage products, besides their critical function for the proper infection process, interplay with the host’s innate immune response through different mechanisms [17]. In particular, PLpro-ISG15 interaction allows the virus to evade the innate immune response through deubiquitination and deISGylation activities of the protease [18, 19]. Interestingly, the process occurs at the same binding cavity as the PLpro known inhibitor, GRL0617 [20].

3. Methods

The code for the cavity-detection guided blind docking (CB-Dock) [8] stand-alone version is freely available at Yang Cao’s Lab webpage [http://clab.labshare.cn/cb-dock/php/manual.php#download].

The customized high-throughput virtual screening pipeline we developed can be accessed at GitHub [https://github.com/tripplab/HTVS].

3.1 Drug selection and modeling

We conducted an extensive scientific literature search for drugs reported as potentially able to prevent SARS-CoV-2 infection. The search included in silico, in vitro, and in vivo studies, covering different stages of the viral cycle. We grouped the reported ligands into five sets: the fusion and viral entry into the host cell (RPA), the polyprotein processing by viral proteases (RPB and RPD), the RNA replication machinery (RPC), and other drugs with alternative or unknown mechanisms (EXT).

We performed the molecular in silico modeling of each ligand’s configuration using the PubChem compound identifier (CID) or, in its absence, using UCSF chimera 1.15 from scratch [21]. The solvent, ions, and other small molecules were removed in all cases, while charges and hydrogens were fixed at neutral pH. Then, ligands were subjected to energy minimization by 10,000 steepest descent steps and 1000 conjugate gradient steps to ensure the proper molecular conformation, saving the final structure in MOL2 format. The next step was to generate the files in PDBqt format using AutoDock Tools, considering the torsional degrees of freedom [22]. We used the PDBqt and MOL2 files as input for the high-throughput virtual screening pipeline. A list of all the ligands studied in this work is shown in Table 1.

ID	Drug	CIDa	References
RPA01	Losartan	3961	[23, 24]
RPA02	Telmisartan	65,999	[23, 25]
RPA03	Arbidol	131,411	[26, 27]
RPA04	Camostat mesylate	5,284,360	[26, 28]
RPA05	Rimantadine	5071	[29]
RPA06	Chloroquine	2719	[26, 30]
RPA07	Hydroxychloroquine	3652	[26, 30]
RPA08	Baricitinib	44,205,240	[31, 32]
RPA09	Colchicine	6167	[33, 34]
RPA10	Disulfiram	3117	[35, 36]
RPA11	Ebselen	3194	[35, 36, 37]
RPA12	Hesperidin	10,621	[38, 39]
RPA13	Qingdainone	3,035,728	[40]
RPA14	Nafamostat	4413	[41, 42]
RPA15	Dipeptidyl nitrile-derivative	Compound 10	[43]
RPB01	Lopinavir	92,727	[26, 44]
RPB02	Ritonavir	392,622	[26, 44]
RPB03	Darunavir	213,039	[29, 45]
RPB04	Cobicistat	25,151,504	[45]
RPB05	Isatin-derivative	Compound 26	[46]
RPB06	Rupinatrivir	6,440,352	[47, 48]
RPB07	E-64	123,985	[49]
RPB08	N3 inhibitor	405,067,310	[50]
RPC01	Ribavirin	37,542	[51, 52]
RPC02	Sofosbuvir	45,375,808	[51, 52]
RPC03	Molnupiravir	145,996,610	[51, 52]
RPC04	Nilotinib	644,241	[53, 54, 55]
RPC05	Saquinavir	441,243	[29, 51, 52, 55]
RPC06	Tipranavir	54,682,461	[51, 52, 55]
RPC07	Lonafarnib	148,195	[55]
RPC08	Tegobuvir	23,649,154	[51, 52, 55]
RPC09	Simeprevir	24,873,435	[51, 52]
RPC10	Filibuvir	54,708,673	[51, 52, 55]
RPC11	Cepharanthine	10,206	[55]
RPC12	Redemsivir	121,304,016	[26, 51, 52]
RPC13	Favipiravir	492,405	[26, 51, 52]
RPD01	rac5c	76,853,649	[19]
EXT01	Ascorbic Acid	54,670,067	[56]
EXT02	Ergocalciferol	5,280,793	[57, 58]
EXT03	Cholecalciferol	5,280,795	[57, 58]
EXT04	Ivermectin	6,321,424	[59]
EXT05	Azithromycin	447,043	[60]
EXT06	Heparin	772	[61]
EXT07	Methylprednisolone	6741	[62]
EXT08	Carvacrol	10,364	[63]
EXT09	Ursolic acid	45,358,157	[64]
EXT10	Oleanolic acid	485,707	[64]

Table 1.

Ligand information, CID number, and reference of 47 drugs with potential activity against the SARS-CoV-2 viral cycle.

^a

In the absence of the CID, the reference to the original investigation and the compound number are provided.

3.2 Target selection and modeling

We included viral and cellular targets involved in the SARS-CoV-2 infection cycle, covering the entry, polyprotein processing, and replication. The targets’ three-dimensional structures were obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB) in PDB format [65]. The complete structure of the spike homotrimer complex (PDB: 6VXX-1-1-1) was retrieved from the CHARMM-GUI Archive-COVID-19 proteins library [66]. We took special consideration to the spike complex given its large size and quaternary structure. We focused on four independent spike-based structures to extend the cavity sampling: the full-length spike’s homotrimer complex, the homotrimer head (S1), 1 S protein monomer, and one isolated receptor-binding domain (RBD).

Water, ions, glycosylations, and co-crystallized ligands were removed from all targets. Charges and hydrogens were fixed at neutral pH using chimera 1.15, their structure optimized, and the final configuration saved in PDB format [21]. In total, 16 structures of 10 targets were curated, as summarized in Table 2.

ID	Target	PDB IDa	SARS-CoV-2 infection step	References
H00	ACE2	1R4L_A	Viral recognition	[67]
H01	ACE2 (B0AT1 closed complex)	6M18_B	Viral recognition	[68]
H02	ACE2 (B0AT1 open complex)	6M1D_B	Viral recognition	[68]
H03	TMPRSS2	7MEQ_A	Viral priming	[69]
H04	Cathepsin B	3AI8_B	Viral priming	[70]
H05	Cathepsin L	2NQD_B	Viral priming	[71]
V01	Spike homotrimer	6VXX1-1-1	Viral recognition	[72, 73]
V01H	Spike homotrimer head	6VXX1-1-1	Viral recognition	[72, 73]
V02	S protein	6VXX1-1-1	Viral recognition	[72, 73]
V02R	S protein’s RBD	6VXX1-1-1	Viral recognition	[72, 73]
V03	Mpro	6LU7_A	Polyprotein processing	[74]
V08	PLpro	7JRN_A	Polyprotein processing	[20]
V04	NSP12	7AAP_A	RNA replication	[75]
V05	NSP7	6M71_C	RNA replication	[76]
V06	NSP8	6NUR_B	RNA replication	[77]
V07	Rdrp-complex (NSP12-NSP7-NSP8)	6M71_ABC	RNA replication	[76]

Table 2.

Structural information and PDB entries of viral (V) and host (H) targets included.

^a

Underscore denotes the chain selected from PDB coordinates files.

3.3 Extended conformational space sampling maximizes the prediction accuracy of the target-ligand complex

Molecular docking is a computational method that allows to sample the conformational space and rank the ligand poses through an energy scoring function. It attempts to generate an optimized target-ligand complex conformation with the lowest binding free-energy change estimate, predicting the interaction of the two molecules in the energy minimum. This task is a cyclic process performed by systematic or stochastic search methods. However, the latter is the choice of preference since it increases the probability of finding an energetic global minimum conformation because the search initiates from different random points [78]. For this reason, the results of two or more molecular docking cycles are not necessarily the same due to the random nature of the conformational search method. Therefore, performing as many cycles as necessary to get as close as possible to the energetic global minimum conformation is essential.

Easy customization of this parameter in the developed high-throughput virtual screening code offers the user the possibility of an exhaustive sampling of the conformational space that maximizes the accuracy of target-ligand complex prediction.

3.4 High-throughput virtual screening pipeline

We have developed in-house bash scripts that integrate the CB-Dock’s cavity-guided blind molecular docking method, which automatically identifies binding sites by calculating putative cavities through a curvature-based detection approach. Molecular docking analysis is conducted in these putative cavities to sample and rank ligand poses and estimate the best target-ligand complex binding energy scores per cycle.

The pipeline has three phases comprised of nested loops, schematized in Figure 1 as a flowchart. First, each target T is subject to a cavity detection step based on a spatial geometry measure of curvature distribution on the protein surface [79]. Cavity identification is achieved by clustering the resulting surface points by density and curvature factor [80]. All cavities are then sorted by size, considering their solvent-accessible surface area. Second, the algorithm automatically configures a docking box for each cavity by defining its center and size, considering the cavity space location and the ligand L size. Finally, in the third step, the blind molecular docking is performed by the AutoDock VINA algorithm [8, 81] for the user-defined top N cavities for each ligand L and each target T. This protocol will be repeated for K-independent rounds.

Figure 1.

Flowchart of the customized high-throughput virtual screening pipeline implemented in this work. Four phases are involved, i) target and ligand molecular modeling (blue), ii) target cavity detection (green), iii) docking box optimization (orange), and iv) target-ligand docking (red).

In our study, we found that the optimal number of independent rounds is $K=30$ since it is at this point that the conformational search converges to the lowest energy binding pose; that is, more rounds do not improve the prediction. The calculations were performed on the top $N=10$ cavities for each $T-L$ pair to significantly extend the cavity and conformational space sampling. The value of these parameters is easily customizable at the top section of the bash script.

3.5 Target-ligand co-crystallization complex prediction

The method we used for automatizing the virtual high-throughput screening process is blind; that is, it does not require any information on the binding site. Hence, we validated its predictions by reproducing the enzymatic targets’ experimental binding complexes. We gathered a set of ligands with available complex co-crystallized data. Eight known enzymatic inhibitors were modeled, optimized, and evaluated under the same methodology conditions as the rest of the ligands included in this study. The ligands in the control set are listed in Table 3.

Target ID	Target name	Ligand	PDB ID	Reference
1. Co-crystallized reproducibility
INH01	ACE2	MLN-4760	1R4L	[67]
INH02	TMPRSS2	4-Guanidinobenzoic acid	7MEQ	[69]
INH03	Cathepsin B	Nitroxoline	3AI8	[70]
INH04	Cathepsin L	4-Bipheylacetyl-cys-(D)-ARG-TYR -N-(2-Phenylethyl) Amide	2NQD	[71]
INV01	Mpro	Narlaprevir	7JYC	[20]
INV02	NSP12	Remdesivir	7BV2	[82]
INV03	NSP12	Favipiravir	7AAP	[75]
INV04	PLpro	GRL0617	7JRN	[20]
2. Negative controla
INV05	Spike	Amantadine	NA	[83]

Table 3.

Modeled ligands to validate that the method is capable of reproducing the co-crystallized complex conformations and previous in vitro findings (negative control).

^a

Does not prevent ACE2-Spike interaction despite inhibiting in vitro SARS-CoV-2 replication.

Furthermore, at this time, a small drug-like co-crystallized molecule in complex with the spike homotrimer does not yet exist. We included amantadine (INV05) in our set as a negative control since it inhibits the SARS-CoV-2 infection but does not prevent spike-ACE2 interaction [83].

3.6 Data analysis and selection criteria

We inspected the top 10 size-ranked putative cavity sites screened for each target. We selected those that either had the active site (targets ACE2, TMPRSS2, cathepsin B/L, Mpro, and NSP12), or were inside a quaternary interface (targets spike, PLpro, and Rdrp). We selected the $T-L$ complex conformation with the best affinity estimation, that is, the conformation with the lowest energy scores after $K=30$ independent rounds for each target-ligand pair. We organized the data in matrix form and analyzed it with the statistical R package function heatmap.2. Rows (ligands) or columns (targets) were scaled to have average = 0 and standard deviation = 1 and generated a Z-score heatmap representation. Finally, we identified potential drugs for repurposing as those ligands with the best energy score estimate at least one standard deviation away from the mean toward more negative values. The data matrix of the VINA scores of the conformation with the lowest scores after $K=30$ independent cycles of each target-ligand pair for known inhibitors and the set of ligands evaluated are provided in the appendix section as Tables A-1 and A-2.

4. Results

4.1 Blind docking correctly reproduces co-crystallized known-inhibitor binding

We found that the $T-L$ complex conformation with the lowest energy for the known co-crystallized inhibitors in our control set successfully reproduces the ligand binding at the active site with an RMSD below 1 Å in most cases, as shown in Figure 2. These findings strongly suggest that the implemented high-throughput blind docking cavity-guided protocol can accurately predict the binding mode of the $T-L$ data in the experimental set.

Figure 2.

Target-ligand complex superimposition of native co-crystallized inhibitors (yellow) and the best-predicted ligand conformation after $K=30$ independent blind docking pipeline rounds (green). The molecular targets (orange) are ACE2 labeled as H00, TMPRSS2 (H03), cathepsin B (H04), cathepsin L (H05), Mpro (V03), NSP12 (V04), and PLpro (V08), created with the visual molecular dynamics (VMD) [84].

4.2 Statistical analysis of the data: Sorting results by target

After doing all the blind docking calculations with an extended conformation sampling, we analyzed the most negative energy scores. We performed a Z-score transformation of the data for each independent column in the matrix (targets T). The graphical representation of the results is shown as a heatmap in Figure 3 using a six-color code based on the Z-score value.

Figure 3.

Target (columns) and ligand (rows) complex docking results. Heatmap of binding free-energy change estimates, using a color-based code according to the Z-score value through column analysis. Targets are grouped as host proteins (blue) and virus proteins (pink). Ligands are grouped by control set (green), potential repurposing drugs (orange), and others (brown). IDs correspond to those defined in Tables 1–3. Black shades represent ligands around the set’s mean. In green shades, ligands with at least one negative standard deviation from the mean. Red shades represent ligands with at least one positive standard deviation from the mean.

Since each column gathers the results for a different target, it is thus possible to identify which ligands had the best scores for each target (in green). It is worth noting that cathepsin L (H05) and PLpro (V08) co-crystallized inhibitors give a good binding free-energy estimate. Most of the co-crystallized inhibitors remained near the mean (in black, with respect to the experimental drug set), except for amantadine (INV05), which presents a positive Z-score value for the spike’s RBD (in red). The latter is concomitant to previous works, where amantadine fails to prevent the spike-ACE2 quaternary interaction [83].

4.3 Nine ligands showed potential for drug repurposing against targets involved in SARS-CoV-2 infection

Out of the 47 drugs screened, nine showed potential inhibition against viral or host targets of the SARS-CoV-2 infection cycle. Saquinavir, simeprevir, nilotinib, an isatin-derivative, telmisartan, tegobuvir, qingdainone, rac5c, and nafamostat achieved the selection criteria. Interestingly, all but rac5c and nafamostat showed the best scores against more than one target. The schematic representation of these results is summarized in Table 4.

ID	Drug	Target	VINA score	Z-scorea
RPC05	Saquinavir	ACE2 (H00)	−12.9	−1.70
RPB05	Isatin-derivative	ACE2 (H01)	−10.7	−2.06
RPC04	Nilotinib		−10.2	−1.72
RPD01	rac5c		−9.9	−1.51
RPC04	Nilotinib	ACE2 (H02)	−10.7	−1.85
RPC09	Simeprevir		−10.4	−1.62
RPC04	Nilotinib	TMPRSS2 (H03)	−8.9	−1.58
RPC05	Saquinavir		−8.8	−1.49
RPC04	Nilotinib	Cathepsin B (H04)	−9.6	−1.64
RPC09	Simeprevir		−10.3	−2.20
RPA02	Telmisartan	Spike (V02)	−9.4	−1.63
RPB05	Isatin-derivative	Spike (V01H)	−10.9	−1.67
RPA13	Qingdainone	Mpro (V03)	−9.4	−1.56
RPC04	Nilotinib		−9.7	−1.80
RPC09	Simeprevir	NSP12 (V04)	−9.1	−1.60
RPA13	Qingdainone	NSP7 (V05)	−7.4	−1.69
RPC08	Tegobuvir		−7.3	−1.59
RPC09	Simeprevir		−7.6	−1.89
RPA02	Telmisartan	NSP8 (V06)	−8.8	−1.65
RPC04	Nilotinib		−8.9	−1.73
RPC05	Saquinavir		−9	−1.80
RPC09	Simeprevir	Rdrp (V07)	−10	−1.96
RPA13	Qingdainone	PLpro (V08)	−9.9	−1.72
RPA14	Nafamostat		−10.1	−1.88
RPC04	Nilotinib		−9.8	−1.65
RPC08	Tegobuvir		−9.8	−1.65

Table 4.

Potential drugs for repurposing with the most negative free-energy change score found and their corresponding Z-score value grouped by the target.

^a

Z-scores were calculated by the target.

Li-Ta	H00	H01	H02	H03	H04	H05	V01	V01H	V02	V02R	V03	V04	V05	V06	V07	V08
INH01	−9.5	−7.3	−7.7
INH02				−6
INH03					−6.2
INH04						−9
INV01											−7.2
INV02												−7.5	−5.5	−6.7	−7.8
INV03												−7.6	−5.7	−6.6	−8.6
INV04																−9.8
INV05							−6	−6	−4.1	−3.4

Table A-1.

VINA scores of the conformation with the lowest scores after 30 independent cycles of each target–ligand pair for the enzymatic known inhibitors included.

Li-Ta	H00	H01	H02	H03	H04	H05	V01	V01H	V02	V02R	V03	V04	V05	V06	V07	V08
RPA01	−9.8	−8.5	−7.9	−6.8	−8	−6.6	−9.7	−9.6	−8.9	−6.8	−7.7	−7.5	−5.6	−6.5	−7.1	−8.8
RPA02	−11.4	−8.8	−10.3	−8.2	−8.1	−7.5	−9.3	−10	−9.4	−7.7	−8.7	−8.2	−6.9	−8.8	−9.2	−9.2
RPA03	−8.6	−6.5	−6.9	−5.9	−6.3	−5.9	−7.6	−7.3	−5.9	−5.7	−6.6	−5.9	−4.7	−5.9	−5.9	−7.5
RPA04	−9.8	−8.7	−8.7	−7.6	−7.6	−7.8	−9.5	−9.4	−7.8	−6.6	−7.4	−7.3	−6	−7	−7.6	−8.4
RPA05	−6.1	−5	−6.6	−4.8	−5.3	−4.9	−6.8	−6.7	−4.6	−4.2	−4.6	−4.2	−4.3	−4.5	−4.3	−6
RPA06	−7.3	−7	−7.6	−5.7	−5.8	−5.9	−7.8	−7.8	−6.4	−5.2	−6.1	−5	−4.5	−5.8	−5.3	−7.3
RPA07	−7.3	−6.4	−7.7	−5.9	−6	−5.9	−7.9	−7.9	−6.4	−5.3	−6.3	−5.6	−4.7	−5.9	−5.7	−7.2
RPA08	−9.5	−8	−7.9	−7.3	−7.9	−6.7	−8.9	−8.9	−8.3	−5.9	−8.2	−6.9	−5.3	−6.2	−6.9	−8.1
RPA09	−8.9	−6.3	−7.1	−5.9	−7	−6.8	−7.4	−7.2	−6.3	−5.7	−7	−6.2	−5.7	−5.9	−6.9	−7.7
RPA10	−4.7	−5	−5.1	−3.9	−4.3	−4	−5.2	−5.2	−4.1	−3.4	−4.5	−3.6	−3	−3.6	−3.7	−4.7
RPA11	−7.8	−6.6	−7.5	−5.6	−6.7	−6	−7.8	−7.8	−6.2	−5.3	−6.6	−5.1	−5.1	−5.4	−5.6	−7.3
RPA12	−11.6	−9.1	−8.9	−8.7	−8.8	−8	−9.2	−9	−7.8	−7.5	−8.9	−8.9	−6.5	−7.1	−8.8	−8.3
RPA13	−12	−7.7	−9.5	−8	−8.5	−7.7	−8.8	−9.1	−7.9	−7	−9.4	−8	−7.4	−7.3	−8.6	−9.9
RPA14	−10.5	−9.2	−9.5	−8.3	−8.3	−7.8	−9.8	−9.9	−8.2	−7	−8.5	−8.4	−6.8	−7.7	−8.9	−10.1
RPA15	−10.9	−9.4	−9.3	−7.8	−7.8	−8.2	−10.1	−9.9	−8.7	−7.2	−8.7	−8	−6	−8	−7.9	−9
RPB01	−11.1	−8.7	−9.3	−7.6	−8.3	−7.7	−9.1	−9.8	−8.8	−7.5	−8.9	−7.3	−5.8	−7.8	−8.3	−8
RPB02	−11.3	−7.8	−8.8	−7.5	−7.7	−7.1	−8.4	−9.3	−8.2	−6.8	−7.7	−7.4	−5.7	−7.4	−7.9	−7.5
RPB03	−10.6	−7.8	−8.5	−7.5	−8.6	−7.6	−9.7	−9.6	−8.6	−6.9	−8	−7.3	−5.9	−6.8	−7.6	−7.8
RPB04	−11.8	−8	−9	−7.6	−8.1	−7.6	−8	−8.6	−8.2	−7.5	−7.9	−7.2	−6.2	−7.4	−8.2	−8
RPB05	−10.4	−10.7	−9	−8.2	−8.9	−8.4	−11	−10.9	−8	−6.9	−8.3	−7.6	−6.2	−8.1	−7.6	−9.7
RPB06	−10.7	−7.7	−9.4	−7.5	−8.8	−7.4	−9.1	−9.7	−7.9	−6.7	−7.8	−7.3	−6.1	−7.6	−8.2	−8
RPB07	−8	−6.3	−7.3	−6.3	−7.6	−6.4	−8.3	−8.1	−6.8	−5.8	−7	−6.6	−5	−5.5	−6.8	−7.6
RPB08	−11.6	−8.6	−9	−7.6	−8.2	−7.6	−8.6	−8.5	−8.4	−7.4	−7.6	−7.2	−6.1	−7.1	−7.9	−8.2
RPC01	−7.6	−6.1	−7.2	−7	−6.7	−5.9	−7.7	−7.7	−6.5	−5.4	−6.6	−5.9	−4.5	−5.1	−6.4	−6.4
RPC02	−10.3	−7.7	−8.3	−7.7	−8.2	−7.3	−9.9	−9.8	−7.4	−7.2	−8.1	−7.2	−6.1	−7	−7.3	−8
RPC03	−8.1	−6.7	−7.8	−7	−7.2	−5.9	−8.6	−8.6	−7.1	−5.8	−7.4	−6.8	−5.1	−5.8	−6.6	−6.9
RPC04	−12.3	−10.2	−10.7	−8.9	−9.6	−8.3	−10	−10.7	−8.8	−7.9	−9.7	−8.6	−7.2	−8.9	−9.1	−9.8
RPC05	−12.9	−8.7	−9.3	−8.8	−8.2	−8.4	−10	−10.4	−8.6	−7.7	−9.1	−8	−6.5	−9	−8.6	−8.4
RPC06	−11	−9.8	−9.8	−7.7	−8.6	−7.6	−9.9	−9.6	−7.6	−7.7	−8.1	−7.6	−6.2	−8	−8.1	−8.1
RPC07	−12.5	−8.9	−9.8	−8.6	−9	−7.8	−8.9	−9.5	−8.4	−7.2	−9.1	−8.7	−6.9	−8.1	−8.8	−8.9
RPC08	−11.5	−9.5	−10.2	−7.9	−9	−8.1	−9.7	−9.7	−8.6	−8.1	−8.6	−8.8	−7.3	−8.1	−8.7	−9.8
RPC09	−10.9	−9.3	−10.4	−8.6	−10.3	−8.6	−4.9	−8.7	−9.1	−8.1	−8.4	−9.1	−7.6	−8.3	−10	−9.3
RPC10	−11.2	−8.7	−9	−7.7	−8.4	−7.9	−9.2	−9.2	−8.6	−7.7	−8.6	−7.8	−6.4	−7.9	−8.1	−8.7
RPC11	−8.5	−7.4	−9.1	−7.7	−8.8	−7.4	−8	−7.9	−8.7	−6.5	−7.7	−8.8	−6.8	−6.6	−8.9	−7.2
RPC12	−10.3	−8.3	−8.9	−7.7	−8.1	−7.3	−8.9	−8.7	−8.1	−6.7	−8.2	−7.4	−5.8	−7.6	−7.5	−7.6
RPC13	−6.3	−5.5	−5.9	−6.2	−5.7	−5.4	−6.2	−6.2	−6.5	−4.4	−5.3	−4.7	−4.1	−4.5	−5.5	−5.5
RPC14	−6.4	−5.5	−6	−6.3	−5.8	−5.4	−6.2	−6.1	−6.6	−4.3	−5.1	−4.6	−3.9	−4.4	−5.4	−5.7
RPD01	−11	−9.9	−9	−8.1	−7.9	−8.3	−8.8	−8.6	−8.1	−7	−7.7	−7.4	−6.4	−7.8	−7.7	−9.1
EXT01	−6.3	−5.1	−5.5	−5.9	−5.8	−5.3	−6.2	−6.2	−6.2	−4.4	−5.1	−4.9	−3.9	−4.3	−5.3	−5.6
EXT02	−10	−8.4	−8.5	−6.5	−7.1	−7.2	−9.8	−9.6	−7.6	−6.3	−7.3	−6.6	−6	−7	−7.1	−7.2
EXT03	−9	−7.8	−8.3	−6.5	−6.9	−6.8	−10.1	−10	−7.1	−5.7	−6.8	−6.1	−5.7	−6.4	−6.8	−6.9
EXT04	−8	−9	−9.6	−8.4	−8.2	−7.9	0.9	−3.1	−8.2	−7.5	−8.2	−8.7	−7.2	−7.2	−9.3	−8.2
EXT05	−6.2	−6.1	−7.2	−6.5	−6.9	−5.6	−6.9	−6.7	−6.1	−5.4	−7.4	−7	−5.3	−5.7	−7.4	−5.2
EXT06	−9.6	−6.9	−7.4	−6.9	−7.6	−7.5	−5.1	−5.5	−6.5	−6	−7.7	−7.6	−5.4	−6	−8	−6.9
EXT07	−9.9	−7.3	−8.1	−6.7	−7.5	−6.3	−7.2	−7.8	−6.7	−6.3	−7.2	−6.6	−5.9	−6.5	−6.5	−7.1
EXT08	−6	−6.8	−6.1	−4.9	−5	−5	−6	−6	−5.3	−4.4	−4.8	−4	−4.4	−5	−4.5	−6.5
EXT09	−11.2	−7	−9.2	−6.7	−8.6	−6.6	−7.8	−7.8	−7.6	−6.4	−7	−7.4	−6.4	−7.2	−7	−6.1
EXT10	−9.7	−7.7	−9	−7.3	−8.3	−6.1	−7.1	−7.8	−6.8	−6.7	−7.7	−7.4	−6.3	−6.9	−7.6	−6

Table A-2.

VINA scores of the conformation with the lowest scores after 30 independent cycles of each target–ligand pair for the set of ligands evaluated.

4.4 Sorting results by ligand: The screened ligands showed a preference for ACE2, Spike, and PLpro targets

Also, we inspected the results by ligand, performing a Z-score analysis by row (ligands). Since the rows gather data from the $L-T$ complex, it is thus possible to identify which ligand had the best scores on any particular target, that is, which target T might be a potential pharmacological target for the ligand L. The results are presented in Figure 4 with a heatmap, using the same color-based code as previously described (see Section 4.2). The targets having most of the ligands in the experimental set with a negative Z-score are ACE2 H00, H01, and H02, and the spike protein structures V01, V01H, and V02 (see Table 2). The results suggest a greater acceptance of those two targets for the ligands as drug-like molecules, at least in the cavities evaluated. These findings are not a minor fact because those targets are directly involved in the first step of the viral infection.

Figure 4.

Target (columns) and ligand (rows) complex docking results. Heatmap of binding free-energy change estimates, using a color-based code according to the Z-score value through row analysis. Targets are grouped as host proteins (blue) and virus proteins (pink). Ligands are grouped by control set (green), potential repurposing drugs (orange), and others (brown). IDs correspond to those defined in Tables 1–3. Black shades represent ligands around the set’s mean. In green shades, ligands with at least one negative standard deviation from the mean. Red shades represent ligands with at least one positive standard deviation from the mean.

The ligands such as arbidol, colchicine, qingdainone, nafamostat, and carvacrol exhibit a binding preference to PLpro (V08). It is important to highlight the essential function of the protease PLpro for processing the viral proteome and evading the host’s innate immune system. In the latter case, PLpro cleaves off post-translational modifications, such as ubiquitin and ubiquitin-like proteins from cell proteins, disrupting the inflammatory signaling pathway necessary for an appropriate immune response [18, 85]. Noteworthy, the potential PLpro inhibitors we have identified in the present work as repurposed drugs form a $T-L$ complex in the cavity where GRL0617 binds, located in the USP domain [20]. The inhibition of this site means blocking the interaction with the ubiquitin-like protein ISG15, evading the immune mechanisms and compromising its canonical enzymatic activity due to the proximity of the assessed site to the active site.

5. Scientific evidence to support our findings

5.1 Saquinavir and simeprevir targeting viral entry and Rdrp quaternary complex formation

Saquinavir is a peptide-mimetic HIV inhibitor. However, some reports suggest potential inhibitory activity against SARS-CoV-2 proteases [86, 87, 88] and other targets involved in the viral infection, such as the Rdrp replication complex [55, 89] and the spike-ACE2 PPI [90]. In our study, saquinavir showed the best energy scores against TMPRSS2, ACE2, and the NSP8-NSP12 interface of the Rdrp complex, as shown in Figure 5. The transmembrane serine protease 2 (TMPRRS2) is essential in several viral infections. Previous reports have shown that the inhibition of this target significantly reduces SARS-CoV-2 entry in lung cells at nM concentrations and therefore the viral infection [91]. Saquinavir also presented the best energy scores against the ACE2 active site, a critical host target needed to initiate entry through the formation of the spike-ACE2 quaternary complex. In this scenario, conformational changes upon ligand binding into the catalytic cavity may shift the relative positions of the receptor’s interface residues that bind to the spike protein and prevent the anchoring of the spike on host cells [92]. However, because saquinavir targets the catalytic site of ACE2, the main activity of this enzyme in the renin-angiotensin system requires further investigation of its biological effect as a competitive inhibitor [93]. In addition, our results show that this drug targets the Rdrp replication complex, which is consistent with the previous results reported in the literature [55, 94]. Interestingly, saquinavir appears to target two essential steps, compromising the entry and viral replication of the SARS-CoV-2.

Figure 5.

Target-ligand complex conformations of potential drugs for repurposing. Molecular docking against viral and host targets relevant in the SARS-CoV-2 infection cycle. A. Superposition of ACE2 target (H00, H01, and H02) docked with saquinavir (cyan), isatin-derivative (red), nilotinib (pink), rac5c (brown), and simeprevir (yellow). B. TMPRSS2 docked with saquinavir (cyan) and nilotinib (pink). C. Spike docked with telmisartan (purple) and isatin-derivative (red). D. Cathepsin B docked with nilotinib and simeprevir (yellow). E. Mpro docked with nilotinib (pink) and qingdainone (orange). F. PLpro docked with nafamostat (green), nilotinib (pink), qingdainone (orange), and tegobuvir (dark orange). G. Superposition of NSP12 and NSP7 and NSP8 cofactors docked with simeprevir (yellow), tegobuvir (dark orange), nilotinib (pink), telmisartan (purple), qingdainone (orange), and saquinavir (cyan), created with the visual molecular dynamics (VMD) [84].

On the other hand, simeprevir also showed the best energy scores on targets relevant to viral entry and replication, including the active cavities of ACE2, cathepsin B, NSP12, and the Rdrp complex interface. We show a molecular visualization of these results in panels A, C, and G of Figure 5. This drug is a protease inhibitor that has presented potent in vitro suppression of SARS-CoV-2 replication at μ M range in Vero E6 cell lines [95]. It is a macrocyclic drug that forms a non-covalent bond within the active site of the hepatitis C virus (HCV) NS3/4A protease, which has a similar three-dimensional arrangement to SARS-CoV-2 Mpro catalytic residues [96]. Simeprevir binds to ACE2 in a quaternary complex inhibition mechanism, analogous to saquinavir, concomitant with the reported disruption of spike-ACE2 PPI [90]. However, the binding does not occur directly at the active site as saquinavir but in the same but larger cavity. Additionally, this drug targets the peptidase activity of cathepsin B, which is a crucial step in spike activation and viral entry. ACE2 was previously proposed as a strategic target to limit viral infection by targeting the cathepsin-mediated entry pathway, decreasing the viral infection efficiency [15, 97]. Simeprevir also showed the best results for the Rdrp complex. Consistent with our results, biochemical assays show low Rdrp replication efficiency after treatment with this drug [95].

5.2 Nilotinib targeting viral entry, polyprotein processing, and Rdrp quaternary complex

Nilotinib is used to treat chronic myelogenous leukemia as a Bcr-Abl tyrosine kinase antagonist. Our results suggest the potential inhibition of six targets involved in the SARS-CoV-2 infection process, including the catalytic cavities of enzyme targets ACE2, TMPRSS2, cathepsin B, Mpro, and the PLpro-ISG15 and Rdrp’s NSP8-NSP12 interfaces. We show a molecular visualization of these results in Figure 5. Reports suggest that nilotinib can inhibit the SARS-CoV and SARS-CoV-2 infection processes, but not MERS-CoV. Interestingly, the latter does not use ACE2 as a cell receptor [98, 99]. This observation is particularly interesting since other reports suggest that nilotinib can destabilize the SARS-CoV-2 spike-ACE2 complex [54]. According to our results, nilotinib might prevent the spike priming and activation since it showed the best energy scores against TMPRSS2 and cathepsin B at the active site cavities. These findings represent a potential inhibition of two independent priming pathways. Moreover, nilotinib potentially inhibits the Mpro and Rdrp complex and is consistent with previous in vitro and in silico results [55, 100].

Interestingly, nilotinib also had the best energy scores against PLpro. In addition to PLpro’s essential protease activity in the processing of pp1a polyprotein, it is also implicated in host immune innate response evasion mechanisms as described in Section 4.4. The inhibition of PLpro decreases the exacerbated immune response, as described by other members of the Bcr-Abl inhibitors family, for example, ponatinib, which protects against cytokine storm in mouse models [101, 102].

5.3 An isatin-derivative and telmisartan targeting SARS-CoV-2 entry

Isatin-derivatives have shown potential antiviral properties, some of them with promising results against HCV, SARS-CoV [103, 104], and SARS-CoV-2 [46]. In particular, the compound 1-(naphthalen-2-ylmethyl)-2,3-dioxoindoline-5-carboxamide inhibits Mpro from SARS-CoV-2. Therefore, we decided to evaluate it against our whole set of targets. It presented the best score against the ACE2 active site, which might disrupt the spike-ACE2 interaction as discussed previously (see Section 5.1). Moreover, it also showed the best energy scores against the spike protein, precisely in the quaternary interface region of the homotrimer complex, and thus a plausible termination of the viral cycle at an early stage in the replication process.

Telmisartan is an anti-antihypertensive. There is evidence of a morbidity and mortality reduction in hospitalized patients infected with SARS-CoV-2 treated with this drug [105]. Telmisartan showed the best energy scores against a spike’s cavity in the homotrimer quaternary interface. Therefore, the isatin-derivative could inhibit two targets involved in the viral entry (spike and ACE2), while telmisartan might prevent the spike homotrimer formation and the Rdrp complex. We show the molecular visualization of these results in panels A, C, and G of Figure 5.

5.4 Tegobuvir, qingdainone, and nafamostat targeting quaternary interface regions

Tegobuvir is a non-nucleoside inhibitor of the NS5B polymerase of HCV. Our results suggest that this drug may prevent the formation of the Rdrp quaternary complex. Previously, in silico results reported tegobuvir as a potential inhibitor of Rdrp active site [106]. According to our data, tegobuvir did not achieve the selection criteria at the Rdrp active site. However, it shows a negative Z-score value at the NSP7-NSP12 interface region, which may compromise the RNA synthesis efficiency of the complex since its importance along with NSP8 for the Rdrp enzymatic activity [107]. Moreover, tegobuvir, nafamostat, and qingdainone presented the best binding free-energy change estimates on a cavity of PLpro in the vicinity of the interface of this target with ISG15, compromising an adequate immune response. In this manner, these drugs could avoid the formation of PLpro-ISG15 and the Rdrp quaternary complexes.

In addition, qingdainone also showed the potential inhibitory activity on Mpro active site, suggesting that this drug might completely disrupt the polyprotein processing stage by targeting both proteases, Mpro and PLpro. We show a molecular visualization of these results in Figure 5.

5.5 Nafamostat and rac5c as potential inhibitors of PLpro and ACE2

We included nafamostat and rac5c in our ligand sets due to evidence suggesting their inhibitory capacity against TMPRSS2 [108] and PLpro [19], respectively. Neither ligand achieved the selection criteria for their expected targets despite being on the borderline with scores of −8.3 and − 9.3 kcal/mol, which indicates the selection criteria’s exhaustiveness. However, nafamostat does achieve the best scores against PLpro’s USP domain, while rac5c presented the best score on the ACE2 active site. We show these results in panels A and F of Figure 5.

6. Conclusions

We have theoretically identified nine drugs or compounds for potential drug repurposing against SARS-CoV-2 through a cavity-based blind molecular docking protocol (Figure 6). Interestingly, seven of them present potential inhibitory activity on multiple targets at different stages of the viral infection cycle, including innate immune evasion. We have implemented an in-house high-throughput virtual screening pipeline that successfully reproduces experimental data and findings from previous works. After the target’s cavity detection and ranking by surface area, we used the pipeline to perform the numerous independent blind molecular docking rounds to achieve a sufficiently extensive conformational target-ligand complex search.

Figure 6.

Repurposing drugs (left) with corresponding potential inhibitory activity on multiple viral or host targets (right).

Experimental design is a critical step in every scientific study, for example, method validation by including a control group. Nonetheless, one has to be wary of the limitations of the methodology employed. In this case, molecular docking can be a good estimator for the most energetically favorable $T-L$ complex. However, the method does not explicitly consider solvent or thermodynamic parameters. Hence, molecular docking results should be taken as the input of other methodologies to further the study, for example, molecular dynamics.

We analyzed the molecular binding predictions through rigorous visualization and Z-score-based statistical algorithms to identify the potential drugs for repurposing. In this context, our findings suggest that:

• Saquinavir and simeprevir could target viral entry and Rdrp complex quaternary formation,

• Nilotinib could target viral entry, polyprotein processing, and Rdrp quaternary complex formation,

• An isatin-derivative and telmisartan could target SARS-CoV-2 entry into the host,

• Tegobuvir, qingdainone, and nafamostat could target quaternary interface Rdrp regions, and

• Nafamostat and rac5c could be potential inhibitors of PLpro and ACE2.

These results are relevant in understanding the SARS-CoV-2 drug’s molecular mechanisms and further clinical treatment development, either at a single or multi-target activity.

Acknowledgments

The authors acknowledge funding from the Consejo Nacional de Ciencia y Tecnología México [grant number 132376] and Fondo Sectorial de Investigación para la Educación [grant number A1-S-17041] to MCT. All computations and analyses reported here were performed at the bmdhpc computing resources of the Biomolecular Diversity Lab (tripplab.com) at CINVESTAV Unidad Monterrey, México.

Conflict of interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

ACE2
COVID-19	Coronavirus disease 19
HCV
HTVS
ISG15
MERS-CoV	Middle East respiratory syndrome coronavirus
Mpro
PLpro
PPI
RBD
Rdrp complex	RNA-dependent RNA polymerase complex
SARS-CoV	Severe acute respiratory syndrome coronavirus
SARS-CoV-2	Severe acute respiratory syndrome coronavirus 2
ssRNA
TMPRSS2	Transmembrane serine protease 2