In silico identification of drug candidates against COVID-19.

The COVID-19 pandemic has caused unprecedented health and economic crisis throughout the world. However, there is no effective medication or therapeutic strategy for treatment of this disease currently. Here, to elucidate the inhibitory effects, we first tested binding affinities of 11 HIV-1 protease inhibitors or their pharmacoenhancers docked onto SARS-CoV-2 main protease (M pro ), and 12 nucleotide-analog inhibitors docked onto RNA dependent RNA polymerase (RdRp). To further obtain the effective drug candidates, we screened 728 approved drugs via virtual screening on SARS-CoV-2 M pro . Our results demonstrate that remdesivir shows the best binding energy on RdRp and saquinvir is the best inhibitor of M pro . Based on the binding energies, we also list 10 top-ranked approved drugs which can be potential inhibitors for M pro . Overall, our results do not only propose drug candidates for further experiments and clinical trials but also pave the way for future lead optimization and drug design.

Introduction

COVID-19 is an infectious disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which endangers the global health and economy. SARS-CoV-2 is a positive-sense, single-stranded RNA virus, which has 79.5% similarity of the genetic sequence to SARS-CoV [1]. Even though the fatality rate of COVID-19 is estimated to be lower than that of SARS (9.5%), its infectivity is higher than SARS-CoV [2]. To date, there have been 30 million confirmed cases and 961,400 deaths from the COVID-19 outbreak as of September 20, 2020 [3]. To fight SARS-CoV-2, a number of nonspecific antiviral drugs have been proposed and tested, such as remdesivir, favipiravir and lopinavir–ritonavir [4,5,6]. However, the debate over the drug's efficacy continues due to the conflicting experimental results. Currently, lack of effective medications or vaccines indirectly leads to an increase of infection numbers.

To discover effective drugs based on the therapeutic protein targets is a strategy to tackle viral threats. There are mainly five therapeutic protein targets, which are spike protein (S protein), angiotensin-converting enzyme 2 (ACE2), main protease (M), papain-like protease (PL), and RNA dependent RNA polymerase (RdRp). First, SARS-CoV-2 S protein interacts with the host cell receptor ACE2 to mediate SARS-CoV-2 entry into host cells [7]. Then, S protein is primed by human protease TMPRSS2. Camostat mesylate, which blocks the activity of TMPRSS2, can treat SARS-CoV-2-infected patients [8]. The lipopeptides EK1C4 was generated as an inhibitor against SARS-CoV-2 S protein-mediated membrane fusion, which was also applied to treat COVID-19 [9]. For ACE2, the uncertain effects of Renin–Angiotensin–Aldosterone System (RAAS) inhibitors remain due to the limited research [10]. After entry into host cells, the viral RNA genome is released and translated into polyproteins. Then, these polyproteins are cleaved into non-structural proteins such as RdRp by M and PL (Fig. 1 ). M and PL mediate viral replication by cleaving viral polyprotein precursors at certain sites [11]. There are 11 M cleavage sites and 3 PL cleavage sites [12]. Inhibition of M can block the synthesis of viral proteins. Hence, the pivotal role in the viral life cycle and the absence of closely related homologues in humans make M and PL attractive as drug targets. In additon, RdRp is also a critical therapeutic protein target catalyzing the synthesis of viral RNA. Following replication and translation of viral RNA, new viral components are assembled and released. So far, numerous drugs have been reported; however, there is still no significantly effective drug for COVID-19. Phytochemicals have been proposed as potential inhibitors against SARS-CoV-2 but are still under research [13]. The co-administered drugs lopinavir and ritonavir, which are HIV-1 protease inhibitors, have been tested for the patients with COVID-19; however, there was no significant benefit observed with lopinavir–ritonavir treatment [5]. This February, the pharmaceutical company Gilead used prodrug remdesivir to successfully cure a COVID-19 patient [14]. However, limited sample size cannot prove this drug is effective for all patients [15]. Even though remdesivir was regarded as one of most promising therapies by the World Health Organization (WHO) [16] and conditionally approved by some countries or regions [17], the details of inhibition need to be elucidated for future drug discovery. Furthermore, there are still many potential drug candidates without clear mechanisms. For example, an investigational drug ebselen has strong inhibitory effect with an IC50 value around 1 μM; however, the target protein and the inhibitory mechanism remain unclear [18,19]. Another prodrug favipiravir was proposed as a potential clinical intervention for COVID-19 but lack of inhibitory mechanism [4].

Fig. 1

Schematic diagram of life cycle of SARS-CoV-2.

To accelerate the drug development of COVID-19, we attempted to virtually examine the efficacy of the proposed agents with computational approaches. Most of the proposed agents are FDA-approved drugs inhibiting HIV-1 protease or the nucleotide-analogs targeting to SARS-CoV-2 M or RdRp. Therefore, we docked these proposed drugs or prodrugs to their potential target proteins, M or RdRp of SARS-CoV-2. As SARS-CoV-2 is an RNA virus highly similar to the SARS virus and some prodrugs like remdesivir were proposed to target SARS originally, we docked the drug candidates to the target proteins from both viruses and compared the docking results. As a result, saquinavir and remdesivir are the best inhibitors for M and RdRp of SARS-CoV-2, respectively. We also performed virtual screening by docking 728 approved drugs from DrugBank (https://go.drugbank.com/) [20,21,[22], [23], [24]] onto SARS-CoV-2 M and proposed 10 top-ranked binders for further experiments and clinical trials (Fig. 2 ). Our results do not only propose drug candidates for further experiments and clinical trials but also pave the way for future lead optimization and drug design.

Fig. 2

Flowchart showing experimental design of this study.

Methods and materials

Protein structure preparation

The structures of HIV-1 protease (PDB ID: 2Q5K) [25], SARS-CoV M (PDB ID: 4MDS) [26], SARS-CoV-2 M (PDB ID: 6LU7) [18], and SARS-CoV-2 RdRp (PDB ID: 6M71) [27] were retrieved from RCSB's Protein Data Bank (https://www.rcsb.org/) [28]. HIV-1 protease and SARS-CoV M were selected as control groups. HIV-1 protease inhibitors were docked on SARS-CoV-2 M and inhibitory effects were compared between original target (HIV-1 protease) and new target (SARS-CoV-2 M). SARS-CoV M was selected as control because we intended to compare the inhibitory effects between two targets which share high genetic sequence similarity [1]. The RNA structure was extracted from PDB file 3H5Y [29]. Then, the structures of RdRp and RNA were combined using Maestro (Schrodinger, version 11.9). To dock the adenosine or other three nucleotides analogue inhibitors onto the active site on RdRps, the nucleotide GMP that would interact with the incoming nucleotide was mutated to UMP or the other corresponding nucleoside 5′-monophosphate.

All the protein structures were prepared by Protein Preparation wizard in Maestro [30]. The workflow of protein preparation contains three steps. The first step is Preprocess, which includes assigning bond orders, adding hydrogens, creating zero-order bonds to metals, creating disulfide, filling in missing side chains using Prime, deleting water molecules beyond 5.00 Å from het groups and generating het states using Epik () [31]. The second step is Optimization, which contains using PROPKA default setting () and performing optimization [32]. The third step is Minimization. This step was performed using the OPLS3e force field [33]. The converge heavy atoms to root-mean-square deviation (RMSD) was 0.30 Å (default setting).

Ligand preparation

The 3D molecular structures of HIV-1 protease inhibitors, the pharmacoenhancers of protease inhibitors, and nucleotide-analog inhibitors were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The 3D molecular structures of approved drugs for virtual screening were retrieved from DrugBank. The 3D molecular structures of Remdesivir-TP, Favipiravir-TP, and Galidesivir-TP were created in Maestro. All the compounds were prepared using Ligprep panel in Maestro. The force field was OPLS3e [33]. The preparation process consists of adding hydrogens, computing correct partial charges, and optimizing the structures.

Ligand-protein docking

To investigate the interaction between ligands and the target proteins and to estimate their binding energies, ligand-protein docking was conducted using Ligand Docking panel in Maestro. After the ligands and the target proteins were prepared using Ligprep and Protein Preparation, respectively, a receptor grid box was generated. For proteases, the receptor grids were generated according to the binding sites of existing inhibitors [18,25,26,34]. For SARS-CoV-2 RdRp and its mutant, the receptor grid was generated according to the position of incoming dNTP, because the inhibitors are supposed to compete for the binding site with incoming dNTP [27]. The size of the receptor grid box was set as default (20 Å). Then, the ligand-protein docking was performed with extra-precision (XP) mode.

MM-GBSA calculation

Prime MM-GBSA is a tool for calculating ligand binding energies within the MM-GBSA (molecular mechanics generalized Born surface area) continuum solvation model. The binding energy () between a protein and a ligand reflects how stably they bind to each other. Therefore, we examined whether an inhibitor tightly bound onto its target protein by calculating the MM-GBSA energies. Here, was estimated using the Prime MM-GBSA module in Maestro. The pose viewer files of the docked complex were uploaded into the MM-GBSA panel. The force field was OPLS3e [33].

Virtual screening

To screen more effective drugs for SARS-CoV-2 M, the structures of the approved drugs (2,635 drugs) in .sdf format were downloaded from DrugBank. Since the molecular weights of 11 HIV-1 protease inhibitors are in the range between 500 and 800, we picked out 728 drugs with molecular weights from 400 to 1,000. Then we docked the 728 drugs onto SARS-CoV-2 protease using XP mode. Then, according to docking scores of the HIV-1 protease inhibitors we examined in this study, we set a cut-off criterion for further screening. The docking scores of squinavir and cobicistat, whose binding affinities are better than the other HIV-1 protease inhibtors or their pharmacoenhancers, are −5.409 (saquinavir) and −6.655 (cobicistat). Hence, to identify potential inhibitors comparable to saquinavir and cobicistat, we set the docking score −5 as a cut-off criterion. There were 346 compounds with docking scores better than −5 in total. Next, we calculated the MM-GBSA energies of those compounds. Based on the MM-GBSA energies, the top 10 best compounds were picked out.

Molecular dynamics simulation

The structures of RNA and RdRp had been combined in Maestro because the complex structures of RNA bound RdRp was just published [35]. To obtain stable protein-RNA complex for the following docking experiments, we performed MD simulations to get the complex with low potential energies. The molecular dynamics (MD) simulations were performed using GROMACS version 2018.1 and CHARMM36 all-atom force field [36,[37], [38], [39],40]. The starting coordinates of the protein-RNA complex were obtained in Maestro. Then we used a cubic box as the unit cell and filled it with water. After adding ions, the complex was minimized for 50,000 steps of steepest descent minimization. Next, the complex was equilibrated using an NVT ensemble (constant Number of particles, Volume, and Temperature) and NPT ensemble (the Number of particles, Pressure, and Temperature). The target temperature for equilibration was 300 K. The last step consists of performing the simulations for 100 ns. Finally, the PDB files with low potential energies were outputted and prepared for docking in Maestro. From protein preparation to MM-GBSA calculation, OPLS3e was used as force field consistently.

Results

Inhibitory effects of HIV-1 protease inhibitors

The coronavirus’ M, which is responsible for proteolytic processing of viral proteins, is one of the significant anti-CoV drug targets. Here, to discover effective drugs for proteases and to demonstrate the inhibitory effects, the interactions between proteases and existing protease inhibitors were investigated via ligand-protein docking. We first selected 11 of HIV-1 protease inhibitors and their pharmacoenhancers [41] to dock onto SARS-CoV-2 M (PDB ID: 6LU7) [18] and SARS-CoV M (PDB ID: 4MDS) [26]. Meanwhile, HIV-1 protease (PDB ID: 2Q5K) [25] was used for docking as a positive control. Based on the docking results, MM-GBSA energies were calculated to describe the binding energies using Prime MM-GBSA. From Table 1 , we find that saquinavir and cobicistat show the best binding energies for SARS-CoV-2 M, −106.17 kcal/mol and −115.61 kcal/mol, respectively. However, cobicistat is a pharmacoenhancer whose major function mechanism is to inhibit human CYP3A, the liver enzymes metabolize medical compounds [42]. Like cobicistat, another HIV-1 protease inhibitor, ritonavir, inhibits liver enzymes and is often used in combination with other medications. Therefore, we propose that saquinaivir should be the best candidate among the 11 docked drugs. Cobicistat also shows the best MM-GBSA energy for SARS-CoV M (−99.54 kcal/mol), which is a little lower than that on SARS-CoV-2 M. For the HIV-1 protease, the MM-GBSA energy of lopinavir is −122.77 kcal/mol, which is the best one among 11 inhibitors. However, the MM-GBSA energies of lopinavir and ritonavir on SARS-CoV-2 M are only −78.10 and −92.36 kcal/mol, respectively, which indicates that lopinavir and ritonavir might not be effective enough for COVID-19 treatment. This result is consistent with the reported conclusion that no benefit was observed with lopinavir–ritonavir treatment on the patients with severe COVID-19 [5].

Table 1

The binding energies (kcal/mol) of 11 HIV-1 protease inhibitors and their pharmacoenhancers onto HIV-1 protease (2Q5K), SARS-CoV M (4MDS), and SARS-CoV-2 M (6LU7).

Inhibitors	Molecular Weight (g/mol)	HIV-1 protease (2Q5K)	SARS-CoV Mpro (4MDS)	SARS-CoV-2 Mpro (6LU7)
Lopinavir	628.8	−122.77	−87.97	−78.10
Indinavir	613.8	−122.53	−78.16	−76.80
Ritonavir	720.9	−119.08	−98.72	−92.36
Nelfinavir	567.8	−96.72	−79.02	−73.96
Darunavir	547.7	−88.62	−61.55	−55.69
Amprenavir	505.6	−76.45	−62.07	−69.57
Tipranavir	602.7	−101.62	−57.69	−67.62
Fosamprenavir	585.6	−72	−57.67	−60.67
Atazanavir	704.9	−116.65	−94.53	−85.98
Saquinavir	670.8	−112.07	−90.35	−106.17
Cobicistat	776	−118.25	−99.54	−115.61

To further understand the inhibitory effect of the top-ranked inhibitor saquinavir, the docking pose and 2D protein-ligand interaction are illustrated in Fig. 3 . There are three hydrogen bonds formed between saquinavir and ASN142, and three hydrogen bonds formed between saquinavir and GLU166. Saquinavir also has a hydrogen bond with both HIS164 and GLN189, respectively. Forming more hydrogen bonds might lead to better binding affinity of saquinavir.

Fig. 3

Docking pose and 2D ligand-protein interaction of saquinavir docked on SARS-CoV-2 M (PDB ID: 6LU7). The pink arrow indicates the hydrogen bond. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Inhibitory effects of nucleotide-analog inhibitors

Since SARS-CoV-2 is an RNA virus, blocking synthesis of viral RNA is another critical antiviral strategy. The processes of RNA replication and transcription are mediated by an RNA-dependent RNA polymerase (RdRp) complex. This complex consists of viral nonstructural proteins (nsps), which contains nsp7, nsp8, and nsp12. Here, we applied a cryo-EM structure of RdRp complex (PDB ID: 6M71) [27] as the receptor for docking experiments. There was no structure of the binding RNA determined in 6M71 when we conducted the docking experiment. However, the binding drugs must interact with RdRp and the elongating RNA molecule simultaneously. To simulate the inhibitory effects of nucleotide-analog inhibitors and to predict the docking pose correctly, we extracted the RNA structure from PDB file 3H5Y, the structure of Norovirus RNA polymerase binding with RNA and combined the RNA with 6M71.

Nucleotide-analog inhibitors are an important group of antiviral agents, which bind onto the active site of DNA or RNA polymerase to compete with nucleotide substrate. Therefore, we tested the inhibitory effects of nucleotide-analog inhibitors on SARS-CoV-2 RdRp. Herein, remdesivir, as a promising adenosine-analog prodrug, was used as main ligand to elucidate the mechanism of inhibitory effects. Meanwhile, we docked the other 11 nucleotide-analog inhibitors on as a control group SARS-CoV-2 RdRp (Table 2 ). Traditionally, molecular docking was mostly applied to simulate the binding between ligands and the target proteins. Usually, the other molecules such as DNAs and RNAs were not taken into account. The proposed binding on the empty active site of RNA polymerases does not reflect the actual binding mechanisms and binding propensities of the drug candidates; therefore, those high-scored drug candidates found to bind to the empty pocket of RdRp cannot be reliable. In our docking experiment on RdRp, the potential ligands are docked onto the 3′-terminal of the RNA chain and the docking results demonstrate that good inhibitors always interact with terminal nucleotides and the binding or catalytic residues on the protein simultaneously. During the docking process, the structures of prodrugs such as remdesivir were all changed into triphosphate (TP) form. As a result, the binding energy of remdesivir-TP is −58.94 kcal/mol, which is better than the other 11 inhibitors and is the only one inhibitor that has better binding energy than the relative substrate ATP (−52.12 kcal/mol). This result indicates that remdesivir-TP can compete with ATP for the active site to terminate the elongation of RNA chain. Fig. 4 shows that there are two hydrogen bonds formed between remdesivir-TP and RU2, which is similar to the interaction between AMP and UMP. Additionally, there are three hydrogen bonds and one Pi-Pi stacking formed between remdesivir-TP and RG8, which suggests that remdesivir-TP has stable interaction with RNA so that it can substitute for ATP to react with RNA and then terminate RNA synthesis.

Table 2

The binding energies (kcal/mol) of selected nucleotide-analog inhibitors docked onto SARS-CoV-2 RdRp (6M71).

Ligand	Nucleoside	Binding energy
ATPa		−52.12
Remdesivir-TP	A	−58.94
Favipiravir-TP	A,G	−45.72
Galidesivir-TP	A	−42.47
GS-4611203	U,T	−36.50
Aphidicolin	C	−41.45
Combivir	T	−38.98
Didanosine	A	−38.05
Zalcitabine	C	−37.62
Stavudine	T	−31.74
Lamivudine-TP	C	−47.54
Carbovir-TP	G	−49.63
(−)FTC-TP	C	−37.42

ATP is the substrate.

Fig. 4

Docking poses and 2D ligand-protein interaction of ATP and remdesivir-TP docked on SARS-CoV-2 RdRp (PDB ID: 6M71). (A) Docking pose of SARS-CoV-2 RdRp with ATP. (B) 2D ligand-protein interaction of ATP at SARS-CoV-2 RdRp binding site. (C) Docking pose of SARS-CoV-2 RdRp with remdesivir-TP. (D) 2D ligand-protein interaction of remdesivir-TP at SARS-CoV-2 RdRp binding site. The pink arrow indicates the hydrogen bond; the green line represents pi-pi stacking; the blue-red line indicates the salt bridge. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Furthermore, to test whether remdesivir-TP has the same binding ability to the mutated RdRp, we docked remdesivir-TP onto the RdRp of mutant P323L (Fig. 5 ). Meanwhile, we docked ATP, favipiravir-TP, and galidesivir-TP as the control group. Notably, even though the binding energy of remdesivir-TP decreases to −47.05 kcal/mol (Table 3 ), the binding energy of remdesivir-TP is still better than those of other two inhibitors. The binding energy of ATP also decreases from −52.12 kcal/mol to −40.43 kcal/mol, which is lower than remdesivir-TP's binding energy. We conclude that remdesivir-TP can still compete for the active site with ATP to show the inhibitory effect on SARS-CoV-2 mutant P323L.

Fig. 5

2D ligand-protein interaction of ATP (A) and remdesivir-TP (B) docked on SARS-CoV-2 RdRp mutant P323L. The pink arrow indicates the hydrogen bond; the blue-red line indicates the salt bridge. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 3

The binding energies (kcal/mol) of three adenosine-analog inhibitors docked onto SARS-CoV-2 mutant P323L RdRp.

Ligand	SARS-CoV-2 mutant P323L RdRp
ATPb	−40.43
Remdesivir-TP	−47.05
Favipiravir-TP	−39.07
Galidesivir-TP	−38.84

ATP is the substrate

Virtual screening of approved drugs

To discover more effective drugs for M, we downloaded a set of approved drugs from DrugBank to screen the potential inhibitors via virtual screening. Based on the calculated binding energies, we identified top 10 ranked compounds as potential drugs (Table 4 ). Acarbosa (DB00284) shows the best binding energy, −99.51 kcal/mol. However, acarbosa is used for treatment and management of type II diabetes [43], and it has never been applied as an antiviral drug. Rutin (DB01698), also named as quercetin-3-O-rutinoside, is the second best inhibitor among ten compounds (−92.44 kcal/mol). The antiviral effect of rutin has been demonstrated against avian influenza strain H5N1 [44]. Accordingly, rutin may be applied as a potential drug for treatment of COVID-19. The third-best drug is lumefantrine (DB06708), with binding energy of −88.91 kcal/mol. Notably, lumefantrine is an antimalarial agent, which also exhibits antiviral effect when used as a combination drug with arthemeter [45]. Dabigatran etexilate (DB06695), with binding energy of −88.60 kcal/mol, can be used to inhibit the formation of blood clots [46]. Dihydroergotamine (DB00320) shows the binding energy of −88.36 kcal/mol, which could be used as a vasoconstrictor [47]. Both dabigatran etexilate and dihydroergotamine have not been reported as antiviral agents [48,49]. Pentagastrin, a synthetic polypeptide to stimulate gastric acid secretion, has similar binding energy (−88.33 kcal/mol) with dabigatran etexilate and dihydroergotamine. Valrubicin (DB00385), cabazitaxel (DB06772), and paclitaxel (DB01229) are all anti-cancer drugs [50,51,52], whose binding energies are −85.96 kcal/mol, −85.57 kcal/mol, and −84.06 kcal/mol, respectively. Among these three anti-cancer drugs, valrubicin has been proved to inhibit SARS-CoV-2 M most recently [53]. The binding energy of plazomicin (DB12615) is −85.01 kcal/mol. Plazomicin is an aminoglycoside antibiotic, and there is no related report about antiviral effect of plazomicin so far [54].

Table 4

Results of virtual screening on SARS-CoV-2 M.

DrugBank Id	Name	ΔG_bind (kcal/mol)	Description	Reference
DB00284	Acarbosa	−99.51	Treatment and management of type II diabetes	[43]
DB01698	Rutin	−92.44	Decrease capillary fragility	[44]
DB06708	Lumefantrine	−88.91	Antimalarial agent	[45]
DB06695	Dabigatran etexilate	−88.60	An anticoagulant that prevents blood clots from forming	[46]
DB00320	Dihydroergotamine	−88.36	A vasoconstrictor	[47]
DB00183	Pentagastrin	−88.33	A synthetic polypeptide that stimulates gastric acid secretion	[55]
DB00385	Valrubicin	−85.96	Treatment of the bladder cancer	[50]
DB06772	Cabazitaxel	−85.57	Treatment of the prostate cancer	[51]
DB12615	Plazomicin	−85.01	Antibacterial activity	[54]
DB01229	Paclitaxel	−84.06	A chemotherapeutic agent	[52]

Among those 10 drugs, acarbosa and rutin showed better binding affinities than the other drugs. By comparing the 2D ligand-protein interactions of M-acarbosa and M-rutin (Fig. 6 ), we found that rutin also interacted with ASN142, GLU166, and GLN189 by forming hydrogen bonds, which was similar to saquinavir. However, acarbosa interacted with other amino acids in binding pocket such as THR24 and THR26. Hence, we speculate that rutin may exhibit similar inhibitory effect as saquinavir.

Fig. 6

2D ligand-protein interaction of acarbosa (A) and rutin (B) docked on SARS-CoV-2 M (PDB ID: 6LU7). The pink arrow indicates the hydrogen bond. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Discussion

We employed protein-ligand docking and virtual screening to search for potent drugs against COVID-19. No drugs have yet been proven to treat this disease effectively. Currently Gilead Sciences’ remdesivir, which has received emergency use authorization (EUA) from the US Food and Drug Administration (FDA) [56], is the top-choice medicine for the treatment of COVID-19. However, there have been conflicting outcomes of the efficacy of remdesivir [57]. Our computational simulations help us further understand the inhibitory roles and mechanisms of remdesivir and other potential drugs against COVID-19. First, as an adenosine-analog prodrug, remdesivir-TP can compete for the active site with substrate ATP, demonstrating that remdesivir-TP can block the replication of viral RNA. While we were preparing this manuscript, the cryo-EM structure of SARS-CoV-2 RdRp with remdesivir was published [35]. Compared with the cryo-EM structure, the nucleoside-analog moiety of remdesivir in the docked and cryo-EM structures are slightly different (RMSD2.14 Å) because the locations of the RNA molecules in two complex structures are slightly different; however, our docking results correctly predict the binding position and pose of remdesivir on SARS-CoV-2 RdRp (Figure S1). Second, compared with the other nucleotide analog inhibitors, remdesivir-TP shows the best binding energy to the wild-type SARS-CoV-2 RdRp. The binding energy of remdesivir-TP to the mutant P323L is not as good as that to wild-type SARS-CoV-2 RdRp, which explains that why remdesivir is not effective for all patients. Even though remdesivir has been proved to inhibit virus infection in a human cell line [58], remdesivir failed its clinical trial which was conducted in China [59]. Another study concluded that remdesivir is so far the most promising drug candidate and a similar drug candidate, favipinavir, has less strong supportive data to back its use [60]. It is consistent with our docking scores. A failure of clinical trial cannot conclude that the docking results are not reliable, as the failure of a clinical trial can be attributed to many factors including the severe side effects of the patients, non-specific binding to other proteins, mutations or drug resistance development. Our results propose the interactions between these compounds and the potential targets, which can help us rule out the non-binding scenario and focus on possible hypotheses. Similarly, lopinavir-ritonavir, which is a combination of HIV protease inhibitors, is still under debate. Despite no effect on COVID-19 patients, lopinavir-ritonavir is still regarded as a promising therapy by WHO [16]. Notably, in our results, lopinavir-ritonavir does not show the best binding energy on SARS-CoV-2 M, but saquinavir and cobicistat do. Although cobicistat is a pharmacoenhancer of protease inhibitor [41], it is not recommended to use with saquinavir [61]. Therefore, we speculate that either of them may be effective for treatment of COVID-19, instead of being used as a combination drug. In addition, we performed ligand-protein docking and MM-GBSA calculation on PL (PDB ID: 6W9C) with 11 HIV-1 protease inhibitors and their pharmacoenhancers. As a result, atazanavir and saquinavir showed better binding energies among 11 compounds (Table S1), which also proved that saquinavir was a potential drug candidate.

Two studies proposed that ebselen inhibits either PL or M based on both bioassays and docking experiments [18,19]. Compared to our data in Table 1, the MM-GBSA of ebselen is relatively low (−53.83 kcal/mol) [62], but the IC50 values suggest it is effective (IC50 is about 1 μM). One of the two studies proposed an inhibitory mechanism where ebselen forms a selenenyl-sulfide bond with Cys112, which is one of the three residues of the catalytic triad of PL. Most of the general docking scoring functions are designed based on non-covalent interactions between proteins and ligands; therefore, the calculated binding energy cannot reflect the irreversible binding stabilized by a covalent bond. Similarly, the nucleotide analogs such as remdesivir-TP also can covalently link to the 3’ end of the binding RNA and terminate its elongation. The docking pose overlapping the substrate binding site and the estimated binding energy comparable to the substrate demonstrate its ability to compete with the coming nucleotide and inhibit the viral replication.

Identifying effective drugs to treat COVID-19 is an urgent and important task. Generally, drug discovery is time-consuming and complicated [63,64]. Drug repurposing is an efficient strategy to obtain effective drugs with low risk, such as using remdesivir or HIV-1 protease inhibitors. However, due to the existing conflicting results of remdesivir and lopinavir-ritonavir, we still need to look for more effective drugs to fight COVID-19. Accordingly, to speed up the drug repurposing, we applied virtual screening with a dataset of 728 approved drugs. Even though the binding energies of top 10 drugs for SARS-CoV-2 M we list are not as good as those of saquinavir and cobicistat, acarbosa and rutin, the first two hits show better binding energies than lopinavir-ritonavir. Additionally, four of the top 10 drugs have antiviral properties, which are rutin, lumefantrine, pentagastrin, and valrubicin. Therefore, this virtual screening result can act as a starting point for further in vitro and in vivo testing. For RdRp, the virtual screening of approved drugs is ongoing.

Selecting drug candidates according to molecular weight is the first step of virtual screening. The molecular weights of 11 HIV-1 protease inhibitors and their pharmacoenhancers in this study are in the range between 500 and 800. To get the similar compounds from approved drugs, we selected 728 approved drugs with molecular weights from 400 to 1000. The limited number of approved drugs in the dataset may not cover all the potential drugs. In future studies, we should scale up the dataset and further group the drug candidates according to the structural features. Another limitation of this study is the limited tests on mutants. We used RdRp mutant P323L to test the inhibitory effects of selected inhibitors. Even the results show that remdesivir-TP can compete with ATP for active site, the inhibitory effects of remdesivir on other mutants are unknown. Therefore, more tests are needed on mutants to compare the inhibitory effects. All the results in this work were obtained by using in silico methods, without experimental confirmation; nevertheless, drug discovery is a process consisting of many steps including target identification, structure determination, docking or virtual screening, bench experiment validation, toxicity examination, clinical trial, and so on, which are difficult and not efficient to be finished in one research project of a single lab. The achievement and conclusion of every one or two steps can provide clues, evidence, and/or hints for the scientific community. Overall, our work demonstrates that saquinvir and remdesivir exhibit better binding energies compared with other potential inhibitors of SARS-CoV-2 M and RdRp. We also proposed 10 potential inhibitors for M via virtual screening, which can provide clues to discover more effective drugs for COVID-19 treatment.

Author contributions

Y.W. collected and analyzed the data. Y.W, L.G.E, and B.K.D carried out the docking experiments and L.L. performed all the MD simulations. Y.W. and L.G.E. drafted the manuscript. Y.W. coordinated the experiments. K.Y.C. and Z-R.X. conceived the project and revised the manuscript. Z-R.X. guided the research. All authors read and approved the final manuscript.

Funding

This research is supported by a start-up grant from the College of Engineering, .

Data availability statement

All datasets generated for this study are included in the article/Supplementary Material [LINK].

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.