In silico identification of widely used and well-tolerated drugs as potential SARS-CoV-2 3C-like protease and viral RNA-dependent RNA polymerase inhibitors for direct use in clinical trials.

Despite strict measures taken by many countries, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) continues to be an issue of global concern. Currently, there are no clinically proven pharmacotherapies for coronavirus disease 2019, despite promising initial results obtained from drugs such as azithromycin and hydroxychloroquine. Therefore, the repurposing of clinically approved drugs for use against SARS-CoV-2 has become a viable strategy. Here, we searched for drugs that target SARS-CoV-2 3C-like protease (3CLpro) and viral RNA-dependent RNA polymerase (RdRp) by in silico screening of the U.S. Food and Drug Administration approved drug library. Well-tolerated and widely used drugs were selected for molecular dynamics (MD) simulations to evaluate drug-protein interactions and their persistence under physiological conditions. Tetracycline, dihydroergotamine, ergotamine, dutasteride, nelfinavir, and paliperidone formed stable interactions with 3CLpro based on MD simulation results. Similar analysis with RdRp showed that eltrombopag, tipranavir, ergotamine, and conivaptan bound to the enzyme with high binding free energies. Docking results suggest that ergotamine, dihydroergotamine, bromocriptine, dutasteride, conivaptan, paliperidone, and tipranavir can bind to both enzymes with high affinity. As these drugs are well tolerated, cost-effective, and widely used, our study suggests that they could potentially to be used in clinical trials for the treatment of SARS-CoV-2-infected patients. Communicated by Ramaswamy H. Sarma.

Introduction

Coronaviruses (CoVs), a large family of single-stranded RNA viruses, can be transmitted from animals to humans, potentially resulting in serious upper respiratory tract disease (Su et al., 2016). The first CoV-related outbreak, caused by the severe acute respiratory syndrome coronavirus (SARS-CoV), was observed in China between 2002 and 2003 (Drosten et al., 2003; Ksiazek et al., 2003). A decade later, Middle East Respiratory Syndrome coronavirus (MERS-CoV) was isolated in patients displaying symptoms of pneumonia in Middle Eastern countries (Zaki et al., 2012). More recently, a novel coronavirus SARS-CoV-2, the causative agent of coronavirus disease 2019 (COVID-19), has emerged in Wuhan, China, becoming a global health issue (WHO, 2020).

Sequence analysis of the full-length genome of SARS-CoV-2 shows 96% identity to bat coronavirus (Wu, Zhao, et al., 2020; Zhou et al., 2020). SARS-CoV-2 contains a linear (+)sense ssRNA genome that encodes proteins commonly produced coronaviruses (Wu, Liu, et al., 2020; Wu, Zhao, et al., 2020) and analyses of its genome organization suggests conserved RNA synthesis machinery and viral life cycle. All coronaviruses consist of genes encoding the following components: replicase (Rep 1a and 1b), spike (S), envelope (E), membrane (M), and nucleoprotein (N) (Graham et al., 2008; Habibzadeh & Stoneman, 2020; Ziebuhr, 2005). They also possess genes encoding non-structural proteins (nsps) placed between these genes (Gorbalenya et al., 2006). The replication of coronavirus in the cell is initiated by the translation of two overlapping open reading frames (ORF1a and ORF1b) (Scandella et al., 2006; St John et al., 2015; Su et al., 2016; Ziebuhr, 2005) to produce pp1a and pp1ab, generated by a (−)1 ribosomal frameshift mechanism (Firth & Brierley, 2012). These polyproteins are then proteolytically processed by two essential viral cysteine proteases, papain-like protease (PLpro, or nsp3) and 3C-like protease (3CLpro, also known as the main protease, Mpro, or nsp5) (Thiel et al., 2003; Ziebuhr, 2004). Cleavage of pp1a and pp1ab by these proteases results in the production of 16 nsps, with PLpro and 3CLpro responsible for cleavage at three and 11 sites, respectively (Prajapat et al., 2020). Nsp12 is the first protein encoded by ORF1b and has viral RNA-dependent RNA polymerase (RdRp) activity (Xu et al., 2003), requiring nsp7 and nsp8 co-factors for maximal activity (Kirchdoerfer & Ward, 2019). The nsp12-nsp7-nsp8 complex is responsible for transcribing both positive and negative RNA strands (Thiel et al., 2003).

The viral life cycle has multiple stages, including entry, replication of genetic material, protein translation, assembly, and release. Strategies for drug development, which includes computer-aided drug design, aim at targeting viral proteins and host receptors in order to interfere with different stages of the coronavirus life cycle (Boopathi et al., 2020; Das et al., 2020; Khan et al., 2020; Li et al., 2020; Sarma et al., 2020). Given the crucial role of RdRp and 3CLpro in the virus life cycle, these proteins are seen as promising drug targets. Similar approaches have yielded promising drugs against viruses such as human immunodeficiency virus (HIV) and hepatitis C virus (HCV) (De Clercq, 2007). Here, we performed in silico drug repurposing to find candidate drugs that may inhibit the enzymatic activity of 3CLpro and RdRp by using the U.S. Food and Drug Administration (FDA) approved drug library. Our study revealed 360 and 850 drugs with high affinity (< −7.5 kcal/mol) to 3CLpro and RdRp, respectively. After the identification of the drugs which may have potential to bind both viral enzymes, we performed detailed literature analysis to determine well-tolerated drugs that have demonstrated in vitro antiviral activity against other viruses. We found several drugs that may be directly used in clinical trials: antibiotic tetracycline and its derivatives, migraine drugs ergotamine and dihydroergotamine, vasopressin receptor antagonists conivaptan and tolvaptan, antiandrogen dutasteride, and trombopoetin receptor agonist eltrombopag. Our results suggest that these drugs may be used alone or in combination as adjuvants with current protocols for the treatment of SARS-CoV-2.

Materials and methods

Protein and ligand preparation for docking simulations

We initially performed pairwise multiple sequence alignment to determine the degree of coene">nservatioene">n amoene">ng SARS-CoV-2 and other coronavirus family members for 3CLpro and RdRp using EMBL-EBI server (Madeira et al., 2019) and JalView (2.11.0). SARS-CoV RdRp (PDB ID: 6NUR) and 3CLpro (PDB ID: 6Y84) (Owen et al., 2020) structures were retrieved from The Protein Data Bank (www.rcsb.org) and used in docking studies. Proteins were prepared for docking simulations via UCSF Chimera, Dock Prep module as described in (Tardu et al., 2016). Excess water molecules and, in case of presence of alternate locations for residues, the ones with the lower occupancies were removed from structures. Any missing side chains were completed by using the Dunbrack 2010 rotamer library. Polar hydrogens were added. Next, non-polar hydrogens were merged with atoms that they are bound. Finally, Gasteiger charges were introduced to structures with Auto Dock Tools (ADT) suite. Structures of drugs were obtained from Zinc15 database catalogue of FDA-approved drugs that are imported from the compounds within the U.S. Environmental Protection Agency’s (EPA) Distributed Structure-Searchable Toxicity (DSSTox) database. The 3948 drugs and inhibitors of HKU4-CoV 3CLpro were prepared for docking by using ADT suite. Active site of 3CLpro containing the catalytic dyad (Cys145 and His41 in SARS-CoV-2) was used as grid center for docking simulations of the protease. The binding site of nsp8 protein was identified as the grid center for RdRp (centering Thr324 and Val330 residues). Grid box sizes for the 3CLpro and RdRp are determined as 4050 Å3 and 9600Å3, respectively. All docking studies were performed via AutoDock Vina (version 1.1.2) (http://vina.scripps.edu/) (Trott & Olson, 2010).

Molecular dynamics simulations and calculation of binding free energy

Protein structures were solvated in a rectangular box with size of 1.04 × 106 Å3 and 3.97 × 105 Å3 using TIP3P water molecules for RdRp and 3CLpro, respectively. Systems were neutralized with sodium and chloride ions and then additional 150 mM Na+ Cl- ions were added to system to simulate physiological conditions (Humphrey et al., 1996). Energy minimization was performed for 20,000 steps by conjugate gradient method. Systems was gradually heated and equilibrated in NPT ensemble for 1.4 ns. Throughout equilibration, constrains on the proteins were gradually removed starting from 2 kcal/mol/Å2. At each stage constrains were reduced by 0.5 kcal/mol/Å2 and system were equilibrated for 0.4 ns. After equilibration, a production run of 20 ns for each complex was performed. Time step was set to 2fs. van der Waals (12 Å cut-off) and long-range electrostatic interactions (via particle-mesh Ewald) were included to calculate the force acting on the system. Production runs were performed at 310 K and under 1 atm pressure using a Langevin thermostat and a Langevin barostat for 20 ns, respectively. All simulations were performed by using NAMD (Phillips et al., 2005) software and CHARMM36m forcefield (Huang et al., 2017). The numerical calculations reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources). All analysis of MD simulations were performed using VMD, and protein representations were prepared by using Pymol (DeLano, 2009). Parameters of drugs for molecular simulations were generated through CHARMM-GUI (Jo et al., 2008; Kim et al., 2017). Structural analysis of proteins and protein-drug interactions were reported by i) calculating Root Mean Square Deviation (RMSD) of proteins, ii) analyzing the contact frequencies between residues and drugs, iii) calculating distances between centroid (geometric center) of drugs and Cα atoms of target residues, and iv) calculating binding free energies of drugs to target proteins as functions of time. RMSD of simulations and contact frequencies were calculated by using RMSD trajectory tool and timeline functions in VMD, respectively.

Binding free energy (BFE) of drugs (8 for 3CLpro and 4 for RdRp) were calculated by molecular mechanics generalized born surface area (MM/GBSA) method. MMPBSA.py script of AmberTools18 was used for MM/GBSA calculations (Case et al., 2018). For each calculation, 1000 frames obtained from 20 ns simulations were used. BFE between drug and protein is calculated based on Equation 1: where Gcomplex represents the energy of protein-drug complex, Greceptor and Gligand represent the energies of protein and drug in their unbound states, respectively. Each energy term consists of mechanical energy term (EMM), solvation energy term (Gsolv), and entropic term (TS) (Equation 2). In detail, EMM includes electrostatic, van der Waals and bounded energy terms; Gsolv includes polar and nonpolar solvation terms, and TS includes vibrational, rotational and translational energy terms.

Results and discussion

Selection of target proteins

Various strategies can be used to develop antiviral drugs, such as blocking viral binding to human cell receptors and inhibiting viral replicatioene">n and self-assembly (Boopathi et al., 2020; Islam et al., 2020; Khan et al., 2020; Li et al., 2020; Sk et al., 2020; Wu, Liu, et al., 2020). Essential SARS-CoV-2 proteins 3CLpro and RdRp (nsp12) are prime targets due to their function in viral polyprotein processing and genome replication, respectively (Gorbalenya et al., 2006; Wu, Zhao, et al., 2020). The FDA-approved drug library was used to screen for the identification of molecules with high affinity to the active site of 3CLpro and nsp8 binding site of RdRp (Figure 1).

Figure 1.

Outline of the study.

Multiple sequence alignments of 3CLpro amino acid sequences from SARS-CoV, MERS-CoV, HKU4, and SARS-CoV-2 revealed the conservation of a catalytic dyad (Cys148 and His41) (Figure 2(A)). The crystal structure of SARS-CoV-2 3CLpro has recently been released (PDB ID: 6y84; (Owen et al., 2020). Comparison of this structure with SARS-CoV (Anand et al., 2003) and MERS-CoV 3CLpro (Tomar et al., 2015) showed high structural identity (Figure 2(B)). Experimental verification of the active site was confirmed by developing inhibitory molecules against HKU4 3CLpro (St John et al., 2015). Based on this information, we targeted the catalytic dyad pocket of 3CLpro for in silico screening (Figure 2(C)).

Figure 2.

The structure of 3 chymotrypsin like protease (3CLpro) A) Sequence alignment of 3CLpro from SARS-CoV, MERS, HKU and SARS-CoV2 viruses. Conserved residues are highlighted in red. Arrows show the conserved residues (His41 and Cys145) which are critical for the protease activity. B) Superimposing the crystal structures of 3CLPro from SARS shown in cyan (PDB ID: 3VB3), MERS shown in green (PDB ID: 4ylu), HKU4 shown in navy (PDB ID: 4yo9) and SARS-CoV-2 shown in white (PDB ID: 6y84) cartoon (left) or surface representations (right). Conserved catalytic residues (His41 and Cys145) are shown in red stick representation. C) Catalytic pocket of SARS-CoV-2 of 3CLpro is shown in red. Protein is shown in white cartoon (left) or surface (right) (PDB ID: 6y84). Catalytic residues His41 and Cys145 are shown with sticks (red) representation.

RdRp requires interactioene">n with nsp7 and nsp8 for maximal activity. The crystal structure of the SARS-CoV RdRp-nsp7-nsp8 complex has been solved (Kirchdoerfer & Ward, 2019). This enabled us to identify the critical region of RdRp that interacts with co-factor nsp8 (Figure 3(A)). Alignment of RdRp amino acid sequences from SARS-CoV and SARS-CoV-2 exhibited 95% identity, while RdRp residues that interact with nsp8 are 100% conserved (shown in rectangular boxes in Figure 3(B)). Thus, we targeted the nsp8 binding site of RdRp for in silico screening with a goal to find molecules that inhibit the interaction between RdRp and nsp8 and, in turn, may interfere with the replisome activity.

Figure 3.

The structure of RNA-dependent RNA polymerase (RdRp) A) RdRp is shown in blue cartoon (left) or surface representations (right), binding partners nsp8 are nsp7 are shown in orange and cyan cartoon, respectively (PDB ID:6nur). B) Sequence. alignment of RdRp from SARS and SARS-CoV2. Conserved residues are highlighted in red. Residues interacting with nsp8 are indicated with rectangle boxes.

Strategy of identification of candidate drugs

FDA approved drugs were retrieved from the ZINC database (http://zinc15.docking.org/catalogs/fda) and prepared for docking simulatioene">ns (Irwin & Shoichet, 2005). Virtual screening was performed as described previously (Armutlu et al., 2009; Çakir et al., 2012; Doruk et al., 2020). AutoDock Vina was used to predict the optimal conformations of the receptor-ligand complex. Binding affinity energies were calculated by assuming a rigid receptor (protein) and a flexible ligand. To determine a cut-off value that could be used for the identification of drugs with high docking binding energies to 3CLpro, we calculated binding energies between previously identified HKU4 3CLpro inhibitors (4f4, rfm, and 4f5) and 3CLpro (St John et al., 2015). As the binding energies of the inhibitors were between −7.5 kcal/mol and −7.0 kcal/mol (Table 1), we set −7.5 kcal/mol as the cut-off value during the selection of drugs from docking simulations. A histogram of Vina binding affinities of all investigated drugs is shown in Figure S1A. We eliminated the possibility of transient interactions between targets and molecules by performing molecular dynamics (MD) simulation for 20 ns. Root-mean-square deviation (RMSD) calculation of Cα atoms for each protein-drug simulation indicated when MD simulations had reached equilibrium. After initial increases in RMSD values during the first 2 ns of simulations, each equilibrated at approximately 1.5–3 Å (Figure S2). The stability of drug-protein interactions was also verified by visual inspection of MD trajectories. Critical amino acid residues of the target proteins that interact with bound drugs were then identified. Additionally, the distance between the geometric center of individual drugs and interacting amino acid residues were measured. Finally, the binding free energy (BFE) and amino acids that contribute to the BFE of target proteins were calculated from trajectories with the molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MM-GBSA) method.

Table 1.

List of top 15 drugs having the best binding affinity to 3CLpro and 3CLpro inhibitors from MERS-CoV.

ZINC ID and/or Drug Name	2D Structure	Vina Binding Affinity (kcal/mol)	Pharmacology / Indication
ZINC000084480349 / Metacycline		−9.1	For the treatment of acute bacterial exacerbations of chronic bronchitis
ZINC000001612996 / Irinotecan		−9.1	For the treatment of metastatic colorectal cancer
ZINC000003932831 / Dutasteride		−8.9	Indicated for the treatment of symptomatic benign prostatic hyperplasia (BPH) in men with an enlarged prostate gland
ZINC000004099009 / Teniposide		−8.9	Teniposide is used for the treatment of refractory acute lymphoblastic leukaemia
ZINC000003925861 / Vorapaxar		−8.9	Vorapaxar is indicated for the reduction of thrombotic cardiovascular events in patients with a history of myocardial infarction or peripheral arterial disease
ZINC000003938684 / Etoposide		−8.9	For use in combination with other chemotherapeutic agents in the treatment of refractory testicular tumors
ZINC000004214700 / Paliperidone		−8.8	Used in the treatment of schizophrenia
ZINC000003831508 / Teniposide		−8.8	Teniposide is used for the treatment of refractory acute lymphoblastic leukaemia
ZINC000012503187 / Conivaptan		−8.7	For the treatment of euvolemic or hypervolemic hyponatremia in hospitalized patients
ZINC000003977777 / Amcinonid		−8.7	Shows anti-inflammatory activity
ZINC000059364574 / Bromocriptine		−8.6	For the treatment of galactorrhea due to hyperprolactinemia, prolactin-dependent menstrual disorders and infertility, prolactin-secreting adenomas, prolactin-dependent male hypogonadism
ZINC000003978005 / Dihydroergotamine		−8.6	For the acute treatment of migraine headaches with or without aura and the acute treatment of cluster headache episodes.
ZINC000003830729/ Doxorubicin		−8.6	Doxorubicin is used to produce regression in disseminated neoplastic conditions like acute lymphoblastic leukemia, acute myeloblastic leukemia, Wilms’ tumor, neuroblastoma, soft tissue and bone sarcomas, breast carcinoma, ovarian carcinoma etc.
ZINC000052955754 / Ergotamine		−8.6	For use as therapy to abort or prevent vascular headache
ZINC000003833846/ Nelfinavir		−8.6	Used in combination with other antiviral drugs in the treatment of HIV in both adults and children
4F4		−7.5	Inhibitor of 3CL in MERS
RFM		−7.5
4F5		−7.0

Docking and MD simulations of selected drugs against SARS-CoV-2 3CLpro

Docking results of 3CLpro revealed 360 molecules with docking binding energies less than −7.5 kcal/mol (top 100 molecules are listed in Table S1). The top 15 molecules, with binding energies < −8.6 kcal/mol, are shown in Table 1. A literature survey of the top hits revealed that these drugs have antiviral, antibacterial, anticancer, and analgesic effects. The following drugs were determined as widely used and well tolerated: tetracycline and its derivatives, dutasteride, ergotamine, dihydroergotamine, bromocriptine, nelfinavir, and paliperidone.

In silico analyses revealed that tetracycline and its derivatives exhibit high affinity to 3CLpro. Docking binding energies of tetracycline, metacycline, doxycycline, oxytetracycline, rolitetracycline, chlortetracycline, and tigecycline were −8.1, −9.1, −8.2, −8.5, −8.2, −8.2, and −6.9 kcal/mol, respectively. In our docking model, tetracycline occupies the catalytic site of 3CLpro and interacts with catalytically active His41 and Cys145 residues (Figure 4(A)). In addition, tetracycline interacts with Ser46, Asn142, and Glu166 residues, forming hydrogen bonds with Ser46 and Glu166 (Table 2). We assessed the persistency of interactions between tetracycline and 3CLpro under physiological conditions by analyzing MD trajectories. We observed that hydrogen bond interactions between tetracycline and Ser46/Glu166 were maintained and additional interactions were formed with Met49 and Gln189 during the simulation (Figure 4(A)). Measuring the distance between critical residues (His41, Met49, and Glu189) and the geometric center of tetracycline further confirmed that the drug kept its initial docking position (Figure 5(A)). We then determined the BFE between tetracycline and 3CLpro from MD simulations using the MM/GBSA method. The BFE of tetracycline was calculated as −14.6 kcal/mol, which is comparable to the BFEs of control molecules 4f4 and 4f5 (−16.71 kcal/mol and −24.96 kcal/mol, respectively) (Figures S3 and S4). Detailed energy decomposition analysis showed that tetracycline had the highest binding affinity to His41 and Met49 located in the active site of the protease (Figure 6), which was similar to the control molecules (Figure S5).

Figure 4.

Binding mode of top 4 candidate drugs to SARS-CoV-2 3CLpro. A) tetracycline B) dihydroergotamine C) dutasteride, D) nelfinavir, and E) paliperidone. Docking position of each drug is shownat left panels. Proteins are shown in ribbon and drugs are shown with stick representations.20 ns MD simulations were performed and the frequency of interacting amino acid residues are shown in middle panels. Right panel shows the amino acid resides of the SARS-CoV-2 3CLpro that interacts with drugs after molecular simulations. The following coloring scheme is used in figures: carbons are white, nitrogens are blue, oxygens are red, sulphurs are yellow, however, carbons in ligands are shown in cyan.

Table 2.

Interacting residues of selected drugs with 3CLPro.

Drug	Interacting residues
Tetracycline	Thr26, His41, Cys45, Ser46, Met49, Asn142, Met165, Glu166, Gln189
Dihydroergotamine	Thr25, Thr26, His41, Thr45, Ser46, Met49, Gly143, Cys145
Dutasteride	Thr25, His41, Ser46, Asn142, Cys145, Glu166, Pro168, Gln189
Nelfinavir	Thr25, His41, Ser46, Leu141, Asn142, Gly143, Ser144, Cys145, Glu166, Pro168, Gln189
Paliperidone	Thr25, His41, Ser46, Leu141, Gly143, Cys145, Glu166, Gln189

Figure 5.

Distance between geometric center of drugs A) tetracycline, B) dihydroergotamine, C) dutasteride, D) nelfinavir, E) paliperidone and Cα atoms of critical residues (His41, Gln189, and Met49) on 3CLpro.

Figure 6.

Contribution of each residue to binding free energy in drug-3CLpro simulations.

Contribution of each residue to binding free energy in drug-3CLpro simulatioene">ns.

Notably, several in vitro screening studies have reported that doxycycline exhibits antiviral activity. Doxycycline is well tolerated, has great distribution to tissues, and has a bioavailability of over 80% when taken orally (Agwuh & MacGowan, 2006; Smilack, 1999). Consistent with our tetracycline results, doxycycline and other derivatives oxytetracycline and tigecycline were shown to exhibit high binding affinity to 3CLpro in a recent in silico study (Wu, Liu, et al., 2020). Doxycycline significantly inhibited viral propagation of vesicular stomatitis virus (VSV), porcine reproductive and respiratory syndrome virus (PRRSV), and replication of dengue viruses (DENV1-4) (Li et al., 2017; Rothan et al., 2014; Wu et al., 2015). These studies suggest that doxycycline may have a broader activity as an antiviral agent, in addition to its well-defined bacteriostatic activity.

Ergotamine and its derivative dihydroergotamine exhibited high affinity to 3CLpro with docking binding energy of −8.6 kcal/mol. In our docking models, dihydroergotamine bound to 3CLpro in a U-shaped form by means of His41, Cys145, Thr25, Thr26, Ser46, and Glu189 residues (Figure 4(B)), with hydrogen bond formation observed with Thr25/Glu189 (Table 2). MD simulations confirmed that interactions between the molecule and Thr25/His41/Met49/Ser46 were strong and stable (Figure 4(B)). Additionally, initial distances were maintained within 1–2 Å fluctuations, with a BFE of −16.22 kcal/mol, as calculated by MM/GBSA analysis (Figure 5(B)). Energy decomposition analyses further showed that Met49/Glu189/Thr25 had the highest contributions to the BFE (Figure 6). Ergotamine is an α-1 selective adrenergic agonist and is commonly used in the treatment of migraine disorders (Silberstein, 1997). Double-blind, randomized, placebo-controlled, parallel-group, multi-center studies indicate that ergotamine and dihydroergotamine are well tolerated and adverse events are generally mild or moderate and transient (Christie et al., 2003; Diener et al., 2002). However, combining ergotamine with protease inhibitors, such as nelfinavir and indinavir, is not advised due to potential cross-interaction through CYP3A4 (Mortier et al., 2001; Rosenthal et al., 1999). Bromocriptine is another semisynthetic ergot alkaloid derivative and functions as a dopamine receptor agonist and serotonin modulator (Kvernmo et al., 2006). Clinically, bromocriptine is used to treat hyperprolactinemia and Parkinson’s disease (Kvernmo et al., 2006; Via et al., 2010). Here, bromocriptine showed high binding affinity against both RdRp and 3CLpro, with −9.2 kcal/mol and −8.6 kcal/mol docking binding energies, respectively. Several in vitro and in vivo studies have shown that bromocriptine inhibits Zika virus NS2B-NS3 protease activity and interferes with dengue virus post-translation and/or RNA synthesis by interacting with the NS3 helicase (Chan et al., 2017; Kato et al., 2016; Yuan et al., 2017). However, a recent study tested bromocriptine in the Vero E6 cellular infection model and found no suppression of SARS-CoV-2 viral replication, even at doses of 100 µM (Lo et al., 2020).

Dutasteride is an oral syene">nthetic 4-azasteroid and displays antiandrogenic activity by blocking type I and II 5α-reductase, enzymes responsible for converting testosterone into dihydrotestosterone (DHT), a primary hormonal mediator that plays a role in the development and enlargement of the prostate gland (Andriole & Kirby, 2003; Clark et al., 2004). Dutasteride was approved by the FDA in 2001 for the treatment of symptomatic benign prostatic hyperplasia (BPH). Large open-label extension studies revealed that dutasteride is well tolerated and drug-related adverse events are reported in less than 2% of patients (Marberger, 2006; Roehrborn et al., 2002; Thomson, 2005). Docking analyses revealed that dutasteride has a high docking binding energy of −8.9 kcal/mol against 3CLpro (Figure 4(C)).

The BFE of nelfinavir, an HIV-1 protease inhibitor drug, to the target was calculated as −8.6 kcal/mol (Figure 4(D)). Interestingly, nelfinavir strongly inhibits replication of SARS-CoV (Yamamoto et al., 2004). Furthermore, a recent study showed that nelfinavir mesylate drastically inhibits cell-to-cell fusion in transiently transfected Vero cells with plasmids expressing SARS-CoV-2 spike protein (Musarrat et al., 2020).

Binding poses of dutasteride and nelfinavir in docking simulations displayed similarities. Both interact with the aforementioned tetracycline and dihydroergotamine binding polar residues via the 1,4-bis(trifluoromethyl) benzene and benzyl moieties of dutasteride and nelfinavir, respectively. While dutasteride forms hydrogen bond interactions with Thr25 and Glu166, nelfinavir can only interact with Glu166 of 3CLpro (Figure 4(C, D)) (Table 2). Contact analysis from MD simulations showed that dutasteride forms stable interactions with Thr25, Met49 and Glu189, and His41 and nelfinavir interacts with Met165 and Glu189 (Figure 4(C, D)). Despite deviations of dutasteride from its initial docking pose, its localization in the active site was verified by visual inspection (Figure 5(C)). Analyzing the distance between critical amino acid residues and nelfinavir indicated that its initial docking orientation was maintained throughout the MD simulation (Figure 5(D)). Met49, Cys145, Met165, and Glu189 are common residues that interact most favorably with dutasteride and nelfinavir, resulting in BFEs of −20.47 and −26.28 kcal/mol, respectively (Figure 6) (Table 3). These observations suggest that both drugs interact strongly with catalytically active Cys145 and other residues in the active site of 3CLpro, and may inhibit protease activity.

Table 3.

MM/GBSA BFE of protease inhibitors and promising drugs to interfere with SARS-CoV-2 3CLpro and RdRp activities.

3Clpro		RdRp
Drug Name	MM/ GBSA BFE (kcal/mol)	Drug Name	MM/ GBSA BFE (kcal/mol)
Tetracycline	−15.19 ± 2.76	Ergotamine	−24.65 ± 3.93
Dihydroergotamine	−16.22 ± 5.41	Eltrombopag	−35.33 ± 3.48
Dutasteride	−20.47 ± 6.29	Conivaptan	−23.83 ± 5.43
Nelfinavir	−26.28 ± 4.48	Tipranavir	−26.08 ± 3.86
Paliperidone	−18.23 ± 4.17
4F4	−16.71 ± 4.66
4F5	−24.96 ± 3.63

Paliperidone is an oral antipsychotic for the treatment of schizophrenia (Janicak & Winans, 2007). Three six-week trials in patients with acute schizophrenia reported that an extended release form of paliperidone is effective and well tolerated, compared with placebo (Davidson et al., 2006; Marder et al., 2007), undergoing limited hepatic metabolism thereby minimizing the risk of hepatic drug-drug interactions (Janicak & Winans, 2007). Paliperidone showed high binding affinity against both 3CLpro and RdRp, with docking binding energies of −8.8 kcal/mol and −8.6 kcal/mol, respectively. Similar to dihydroergotamine, paliperidone binds to 3CLpro in a U-shaped conformation and interacts with polar residues Thr25, His41, Met49, and Gln189 (Figure 4(E)). Consistent with interactions observed in the docking conformation, paliperidone interacted with these polar residues through almost the entirety of the MD simulation (Figure 4(E)). Distance analysis between critical residues and paliperidone suggested that the molecule adjusted its position as the simulation progressed. While paliperidone kept its initial distance to His41 and Met49, it got closer to Gln189 during the simulation (Figure 5(E)). Paliperidone exhibited high affinity with a BFE of −18.23 kcal/mol as assessed by MM/GBSA analysis. Additional energy decomposition analysis confirmed that Glu189, Met49, and His41 contributed favorably to the binding energy (Figure 6).

We validated our results by docking other antiviral drugs lopinavir and ritonavir, which have docking binding energies of −7.7 kcal/mol and −6.8 kcal/mol to 3CLpro, respectively. Lopinavir and ritonavir are antiretroviral protease inhibitors used in combination with other antiretrovirals in the treatment of HIV-1 infections. In vitro studies have shown that a combination of lopinavir and ritonavir can inhibit the replication of MERS-CoV and SARS-CoV to exert antiviral effects (Dyall et al., 2017). There are also several reports on the use of Kaletra (lopinavir/ritonavir) for the treatment of COVID-19 infections (He et al., 2020; Lim et al., 2020).

MD simulations of selected drugs showed that interactions between these drugs and 3CLpro are not transient but stable. Furthermore, comparable or better BFE values of selected drugs with respect to kene">nown inhibitors indicate that these drugs have high potential to interfere with 3CLpro activity. We also determined critical residues responsible for the high binding affinity of drugs to the protease, which may help to develop novel inhibitor molecules through rational drug design and quantitative structure-activity relationship (QSAR) studies.

Interestingly, seven of the top 20 drugs identified in docking studies have cytostatic activity and are used in cancer therapy: irinotecan, etoposid, teniposide, doxorubicin, epirubicin, and proscillaridin. It is well known that HIV protease inhibitors exert anticancer activity via two main mechanisms; endoplasmic reticulum stress-unfolded protein response pathway and Akt inhibition (Gantt et al., 2013; Maksimovic-Ivanic et al., 2017). Nelfinavir has also been shown to have in vitro anticancer activity (Gantt et al., 2013). HIV protease inhibitors, including nelfinavir and lopinavir, are peptidomimetics with similar structures and are expected to bind to the nelfinavir docking site in 3CLpro. Although several anticancer drugs have displayed in silico interactions with 3CLpro (≤ −8.6 kcal/mol), the use of these types of molecules as antiviral therapeutics may not be recommended due to their serious adverse or toxic effects.

Dockings and MD simulations of tested drugs against RdRp

Docking analyses revealed 850 drugs with docking binding energies lower than −7.5 kcal/mol against the nsp8 binding sites of RdRp (top 100 drugs are listed in Table S2). The top 15 compounds with docking binding energies lower than −9.2 kcal/mol are shown in Table 4, and include ergotamine, dihydroergotamine, bromocriptine, tipranavir, conivaptan, and eltrombopag, which are all well tolerated and widely used.

Table 4.

List of top 15 drugs having the best binding affinity to RdRp.

ZINC ID and/or Drug Name	2D Structure	Vina Binding Affinity (kcal/mol)	Clinical Trial
ZINC000003932831 / Dutasteride		−9.9	Indicated for the treatment of symptomatic benign prostatic hyperplasia (BPH) in men with an enlarged prostate gland
ZINC000052955754 / Ergotamine		−9.8	For use as therapy to abort or prevent vascular headache
ZINC000004198852 / d-Tubocurarine		−9.5	Used as a diagnosis agent for myasthenia gravis, and also to facilitate the intubation after induction of anesthesia in surgical procedure
ZINC000003784182 / Differin		−9.5	For the topical treatment of comedo, papular and pustular acne (acne vulgaris) of the face, chest or back.
ZINC000011679756 / Eltrombopag		−9.5	Thrombopoietin receptor agonists are pharmaceutical agents that stimulate platelet production in the bone marrow
ZINC000012503187 / Conivaptan		−9.4	For the treatment of euvolemic or hypervolemic hyponatremia in hospitalized patients
ZINC000100016058 / Tipranavir		−9.3	For combination antiretroviral treatment of HIV-1 infected adult patients with evidence of viral replication, who are highly treatment-experienced or have HIV-1 strains resistant to multiple protease inhibitors
ZINC000003978005 / Dihydroergotamine		−9.3	For the acute treatment of migraine headaches with or without aura and the acute treatment of cluster headache episodes
ZINC000064033452 / Lumacaftor		−9.3	lumacaftor is indicated for the treatment of cystic fibrosis (CF) in patients age 6 years and older who are homozygous for the F508del mutation in the CFTR gene.
ZINC000011681563 / Netupitant		−9.2	for use in combination with palonosetron for the prevention of acute and delayed vomiting and nausea associated with cancer chemotherapy including highly emetogenic chemotherapy
ZINC000253632968 / Simeprevir		−9.2	Indicated for the treatment of adults with chronic hepatitis C virus (HCV) infection:
ZINC000000538658 / Tolvaptan		−9.2	Treatment of symptomatic and resistant to fluid restriction euvolemic or hypervolemic hyponatremia associated with congestive heart failure, and cirrhosis.
ZINC000001757652 / Silibinin		−9.2	Currently being tested as a treatment of severe intoxications with hepatotoxic substances
ZINC000001550477 / Lapatinib		−9.2	Indicated in combination with capecitabine for the treatment of patients with advanced or metastatic breast cancer
ZINC000001542113 / Vilazodone		−9.2	Vilazodone is approved for treatment of major depressive disorder

Notably, ergotamine and its derivatives dihydroergotamine and bromocriptine exhibited high affinities against RdRp, with docking binding energies of −9.8 kcal/mol, −9.3 kcal/mol, and −9.2 kcal/mol, respectively. As explained above, these drugs also exhibited high binding affinities (< −8.6 kcal/mol) against 3CLpro. The binding mode of ergotamine showed that it forms hydrogen bonds with backbone oxygen atoms of Thr324, Phe326, and Pro328. In addition, ergotamine interacts with RdRp via hydrophobic interactions with Val330 and Val398 (Figure 7(A)). A complete list of amino acid residues that interact with ergotamine is given in Table 5. We investigated dynamic interactions between ergotamine and RdRp by performing 20 ns MD simulations. Contact analysis from MD showed that Thr324, Val330, and Val398 were critical for drug binding (Figure 7(A)). We also monitored distances between Val398/Thr324 and the geometric center of bound drugs. While ergotamine kept its position in the nsp8 binding site, it got closer to Val398 as the simulation progressed, with a BFE of −24.65 kcal/mol (Figure 8(A)) (Table 3). Although we did not have a reference molecule to gauge the BFEs of our selected molecules for RdRp, the values calculated for 3CLpro indicate that a BFE of approximately −20kcal/mol can be considered as high affinity. Detailed energy decomposition analysis indicated that Val398 and Thr324 interacted most strongly with ergotamine (Figure 9). Conservation of these interactions during MD simulations and their contribution to BFE suggests that ergotamine can stably interact with the nsp8 binding site of RdRp. These results indicate that ergotamine and dihydroergotamine are possible candidates for further in vitro testing and clinical evaluation as an anti-SARS-CoV-2 agents.

Figure 7.

Binding mode of top 4 candidates drugs to RdRp. A) ergotamine B) eltrombopag C) conivaptan and D) tipranavir. Docking position of each drug is shown at left panels for each drug. Proteins are shown in ribbon and drugs are shown with stick representations. 20 ns MD simulations were performed and the frequency of interacting amino acid residues are shown in middle panels. Right panel shows the amino acid resides of the SARS-CoV-2 3CLpro that interacts with drugs after molecular simulations. The following coloring scheme is used in figures: carbons are white, nitrogens are blue, oxygens are red, sulphurs are yellow, however, carbons in ligands are shown in cyan.

Table 5.

Interacting residues of selected drugs with RdRp.

Drug	Interacting residues
Ergotamine	Leu270, Thr324, Phe326, Val330, Ala379, Val398
Eltrombopag	Tyr273, Val329, Val330, Arg331, Val398
Conivaptan	Leu270, Leu271, Tyr273, Thr324, Phe326, Leu329, Val330, Val398, Val675
Tipranavir	Leu270, Leu271, Tyr273, Thr324, Leu329, Val330, Arg331, Ala379, Phe396, Val398, Met666

Figure 8.

Distance between geometric center of drugs A) ergotamine, B) eltrombopag, C) conivaptan, D) tipranavir and Cα atoms of critical residues (Val398 and Thr324) on RdRp.

Figure 9.

Contribution of each residue to binding free energy in drug-3CLpro simulations.

Contribution of each residue to binding free energy in drug-3CLpro simulatioene">ns.

Eltrombopag and conivaptan showed high binding affinities of −9.5 kcal/mol and −9.4 kcal/mol, respectively, to the nsp8 binding site of RdRp. While eltrombopag can form hydrogen bonds with Thr324, Pro328, and Leu329, and make hydrophobic interactions with Val330 and Val398, conivaptan interacts with only Leu270 and Val330 by hydrophobic interactions (Figure 7(B, C)) (Table 5).

Contact analysis indicated that all these polar and hydrophobic interactions were maintained throughout MD simulations for both drugs (Figure 7(B, C)). Measuring distances between the drugs and amino acid residues clearly showed that eltrombopag maintained its original docking positioene">n perfectly while conivaptan interacted with only Val330 (Figure 8(B, C)). BFEs of conivaptan and eltrombopag to RdRp were calculated as −23.83 kcal/mol and −35.33 kcal/mol by MM/GBSA analysis, respectively (Table 3). Analysis of the contributions of critical amino acid residues to BFEs indicated that Val398/Thr324 and Val329/Val330 had the highest impact on the binding of conivaptan and eltrombopag, respectively (Figure 9). All calculations and analysis suggest that eltrombopag can strongly bind to the nsp8 binding site of RdRp and may attenuate its polymerase activity. Interestingly, eltrombopag and conivaptan also showed high binding affinities of −8.2 kcal/mol and −8.7 kcal/mol, respectively, against 3CLpro.

Eltrombopag is a human thrombopoietin receptor agonist that is used to treat low blood platelet counts in adults with chronic immune (idiopathic) thrombocytopenia (ITP) (Erickson-Miller et al., 2009). The metabolism and disposition of eltrombopag were studied after a single oral administration of a solution dose of [(14)C] eltrombopag (75 mg, 100 μCi) and was shown to be a well-tolerated drug (Deng et al., 2011). In vitro analysis revealed that eltrombopag can inhibit severe fever with thrombocytopenia syndrome (SFTS) virus replication. This retrovirus causes SFTS, a newly recognized viral hemorrhagic fever that emerged in rural areas of China and spread to Japan and South Korea (Yu et al., 2011). In a recent in vitro screening study, eltrombopag was shown to have antiviral activity against SARS-CoV-2, with an IC50 value of <10 μM (Jeon et al., 2020).

Conivaptan inhibits both isotypes of the vasopressin receptor (V1a and V2) and is used for the treatment of euvolemic or hypervolemic hyponatremia (e.g. the syndrome of inappropriate secretion of antidiuretic hormone, or in the setting of hypothyroidism, adrenal insufficiency, pulmonary disorders, etc.) in hospitalized patients. Placebo-controlled, randomized, double-blind, multi-center studies have shown that oral conivaptan is well tolerated in cases with hyponatremia (Annane et al., 2009; Ghali et al., 2006). Conivaptan has antiviral activity against dengue (DENV) and Zika (ZIKV) viruses at high concentrations (IC50 >100µM) (Montes-Grajales et al., 2020). Another AVP receptor 2 antagonist, tolvaptan, also has a high binding affinity against RdRp of −9.2 kcal/mol. Like conivaptan, tolvaptan is also used to treat hyponatremia associated with various conditions.

Tipranavir is a sulfonamide-containing dihydropyrone and a nonpeptidic protease inhibitor that targets HIV protease. Interestingly, tipranavir shows high affinities −8.4 kcal/mol and −9.3 kcal/mol against both 3CLpro and RdRp, respectively. The binding mode of tipranavir to RdRp suggests that the drug forms a hydrogen bond interaction with backbone oxygen atom of Leu271 (Figure 7(D)). In addition, fluorine atoms on the trifluoromethylpyridine moiety of tipranavir are within hydrogen bond distance with Arg331. Val330, Phe396, and Val398 residues enable hydrophobic interactions between the drug and RdRp (Table 5). While hydrophobic interactions were maintained in MD simulations, interactions with Leu271 and Arg331 were not conserved, especially after the first 5 ns of the simulation. Instead, the drug started to interact more frequently with Thr324 and Leu389. Distance analysis between tipranavir and amino acid residues showed that tipranavir gets closer to Val398 and Thr324 (Figure 8(D)). BFE of tipranavir to RdRp was calculated as −26.08 kcal/mol by MM/GBSA analysis (Table 3). Decomposition of BFE further confirmed that tipranavir strongly binds to Val398, Leu387, Val330, and Thr324 (Figure 9).

Taken together, our docking and MD simulation analyses indicate that ergotamine, dihydroergotamine, eltrombopag, conivaptan, and tipranavir bind to the nsp8 binding site of RdRp and could be used to prevent nsp8 binding. These molecules are therefore strong candidates for use in therapies to attenuate the maximum activity of RdRp in SARS-CoV-2.

Conclusion

Our in silico analyses revealed several candidates as potential therapeutics against SARS-CoV-2. Tetracycline and its derivatives exhibited high affinity to 3CLpro. These antibiotics are cheap, well tolerated, and the patients may benefit from not only their antibacterial and anti-inflammatory effects, but also their antiviral activity. Dutasteride also exhibited high affinity to 3CLpro. Dutasteride has been shown to down regulate the expression of transmembrane serine protease 2 (TMPRSS2) (Lucas et al., 2008; Mostaghel et al., 2014). It was recently shown that SARS-CoV-2 uses the angiotensin-converting enzyme 2 (ACE2) for entry and TMPRSS2 for S protein priming. SARS-CoV-2 cell entry was blocked by TMPRSS2 protease inhibitor (Hoffmann et al., 2020). COVID-19 patients may benefit from dutasteride due to its effect on TMPRSS2 and ACE2 expression, and its possible antiviral activity. Additionally, ergotamine, dihydroergotamine, bromocriptine, conivaptan, eltrombopag, paliperidone, and tipranavir exhibited strong interaction potential with both 3CLpro and RdRp, and are also among the promising candidates (Table S3). As these drugs are cost-effective and widely used, they can be utilized alone or in combination as an adjuvant with currently used protocols in the treatment of SARS-CoV-2 infection. Additionally, each of the above drugs can be selected for different groups of patients based on their primary effects. For example, dutasteride could be selected for use exclusively in males because of its antiandrogenic effects and eltrombopag could be selected for COVID-19 patients displaying thrombocytopenia. As these drugs have high affinity to both target proteins, these compounds should be less sensitive to genetic variations in SARS-CoV-2 3CLpro and RdRp due to mutations arising during circulation of the virus amongst populations. Because of their target proteins, the above drugs should exert their antiviral effects at the early–mid stage of infection and impair the viral reproductive cycle.

As there are currently no approved therapeutics and vaccines for SARS-CoV-2 infection, in silico repurposing of approved drugs targeting viral proteins still represents a promising approach for discovering new antiviral molecules. Although the results of in silico analysis must be confirmed with in vitro studies and clinical testing, repurposing of commonly used, well-tolerated drugs can help to expedite the drug discovery process during urgent situations.

Click here for additional data file.