Virtual screening of approved drugs as potential SARS-CoV-2 main protease inhibitors.

The global emergency caused by COVID-19 makes the discovery of drugs capable of inhibiting SARS-CoV-2 a priority, to reduce the mortality and morbidity of this disease. Repurposing approved drugs can provide therapeutic alternatives that promise rapid and ample coverage because they have a documented safety record, as well as infrastructure for large-scale production. The main protease of SARS-CoV-2 (Mpro) is an excellent therapeutic target because it is critical for viral replication; however, Mpro has a highly flexible active site that must be considered when performing computer-assisted drug discovery. In this work, potential inhibitors of the main protease (Mpro) of SARS-Cov-2 were identified through a docking-assisted virtual screening procedure. A total of 4384 drugs, all approved for human use, were screened against three conformers of Mpro. The ligands were further studied through molecular dynamics simulations and binding free energy analysis. A total of nine currently approved molecules are proposed as potential inhibitors of SARS-CoV-2. These molecules can be further tested to speed the development of therapeutics against COVID-19.

Introduction

In December 2019, a cluster of atypical pneumonia cases was reported in Wuhan, China (Wu et al., 2020; Zhu et al., 2020). Coronavirus disease (COVID-19) rapidly spread to other countries and by March 11th, the disease was declared a pandemic by the World Health Organization (WHO) (Cucinotta and Vanelli, 2020). According to the WHO, to date (Situation Report 83, April 15th, 2020) there are 1914916 COVID-19 cases and more than 120,000 COVID-19-associated deaths worldwide (World Health Organization, 2019). The causal agent of COVID-19 is a SARS-related coronavirus named SARS-CoV-2, which was likely passed to humans in live-animal markets (Jalava, 2020). There is strong evidence that SARS-CoV-2 originated in bats but was likely transmitted to humans through an intermediate animal, probably the pangolin (Lam et al., 2020).

SARS-CoV-2 is a betacoronavirus closely related to MERS-CoV and SARS-CoV, which have also caused outbreaks with pandemic potential (Chen et al., 2020a). These outbreaks have affected the global economy, causing high economic losses (Tripp and Tompkins, 2019). The SARS-CoV-2 genome is a single-stranded positive-sense RNA of about 30 kb in length and contains at least six open reading frames (ORFs) that code for a minimum of 16 non-structural proteins and 4 structural proteins (Chen et al., 2020a; Lai et al., 2020). ORF 1a/b is translated into a large protein that undergoes extensive proteolytic processing to yield the replicase complex, which mediates viral transcription and replication (Bartlam et al., 2005; Chen et al., 2020a). The protease responsible for the proteolytic processing is the main protease (Mpro) or 3C-like protease (3CLpro), which is matured by auto-cleavage into the dimeric active conformation (Xia and Kang, 2011). Given the relevance for the viral replication cycle, Mpro has been proposed as a target in the development of inhibitors against coronaviruses (Bartlam et al., 2005; Yang et al., 2006).

Drug repurposing, also known as drug repositioning, is the use of an active pharmaceutical ingredient to treat a novel medical condition different from the original intended condition and has arisen mainly by serendipity when beneficial off-target or secondary effects are noticed (Pacios et al., 2020; Pushpakom et al., 2019). The use of currently approved drugs to treat different diseases has the advantage of assuring medical safety because the drugs have already been tested in animal models and undergone clinical trials. Additionally, the infrastructure to manufacture at large-scale is already in place (Cha et al., 2018; Pushpakom et al., 2019). Drug repurposing is also a strategy that has been used to discover novel antibiotics or antiviral drugs (Dyall et al., 2018; Pacios et al., 2020). In the case of SARS-CoV-2, many drugs with repurposing potential are already being tested (Li and De Clercq, 2020). The attractiveness of repurposing has led to the evaluation of at least 35,000 drugs for more than one medical condition (Baker et al., 2018).

Another advantage of drug repurposing is a quick approval in emergencies such as the current COVID-19 pandemic. Taking this into consideration, we performed in silico evaluation of a set of approved drugs as potential inhibitors of Mpro from SARS-CoV-2; our findings show that several molecules warrant further analysis as treatment options against COVID-19.

Methods

SARS-CoV-2 genome sequences retrieval and homology modeling of Mpro

A total of 111 SARS-CoV-2 genome sequences were retrieved from the GISAID platform (Shu and McCauley, 2017) and aligned with Clustal Omega through the UGENE platform (Okonechnikov et al., 2012; Sievers and Higgins, 2014). For homology modeling, the BetaCoV/Wuhan/WIV02/2019 genome was analyzed with VGAS (Zhang et al., 2019) to predict the Open Reading Frame (ORF) corresponding to ORF1a, which contains the Mpro sequence. This sequence was used to predict the structure of Mpro in its biologically active conformation (dimer) by using Modeller (Ho et al., 2015; Webb and Sali, 2016); the following structures were used as templates (PDB ids): 2AMD, 1WOF, 2AMQ, 2D2D, 3E91, and 3EA7 (Yang et al., 2005). A total of 20 models were generated and the DOPE (Discrete Optimized Protein Energy) score was used to select the best structural model. Global and local structural quality was evaluated with QMEAN, which is a scoring function that measures the global and local quality of protein models, estimating the degree of structural ‘nativeness’. QMEAN uses a linear combination of structural descriptors that include long-range interactions, torsion angles, and solvation potential. Scores calculated form the structural descriptors are transformed into Z-scores to compare them with high-resolution crystal structures. QMEAN is available in the SWISS-MODEL server (Benkert et al., 2011; Waterhouse et al., 2018). Sequence conservation analysis was done with Chimera (Pettersen et al., 2004).

Molecular dynamics simulation of Mpro

The predicted structural model was submitted to the CHARMM-GUI server to prepare the system (Brooks et al., 2009; Jo et al., 2014, 2008; Lee et al., 2016). The Solution Builder module was used to prepare the protein inside a water cube (TIP3P model) and potassium chloride (KCl) was used to neutralize the system charge and to adjust the salt concentration to 0.15 M. The CHARMM36m force field was used and input files for GROMACS were generated and downloaded (Huang et al., 2017). The molecular dynamics simulation was performed with GROMACS (Abraham et al., 2019, 2015) in three stages: first, a minimization stage (steepest descent) consisting of 5000 steps was performed to eliminate major atomic clashes in the system. Then, an equilibration stage was performed in which protein movement was constrained to allow the solvent and ions to contact the protein. Harmonic force constants of 400 kJ mol−1 nm-2 for protein backbone and 40 kJ mol−1 nm-2 for sidechains were used, with a total equilibration time of 250 ps and a time step of 1 fs at 310 K. Lastly, the production stage was performed without position restrains for a total of 100 ns with a timestep of 2 fs at 310 K. The resulting trajectory was visualized and analyzed with VMD (Humphrey et al., 1996) and GROMACS. Conformers from the resulting trajectory were clustered with the GROMOS method (Daura et al., 1999), which generates non-overlapping clusters. Briefly, the method starts with the RMSD calculation between all pairs of conformers in the trajectory. The first cluster is formed when the biggest set of conformers sharing less than the cut-off is detected. These conformers are eliminated from the pool of clusters. The process is repeated for remaining conformers. A cut-off of 0.18 nm was used to generate approximately 20 clusters from the 1000 conformers resulting from the molecular dynamics simulation.

Docking-mediated virtual screening

The strategy to discover molecules that potentially inhibit Mpro involved molecular docking of each molecule into the active site of Mpro. iDock (Li et al., 2012) was used as the molecular docking engine and the procedure was first validated by redocking ligand 9NI to the SARS-Mpro structure (PDB id: 2AMD). iDock is a molecular docking program that was developed with a focus on virtual screening, being capable of processing a massive number of ligands in a relatively short time. It is similar to AutoDock Vina (a fast and accurate evolution of AutoDock) with regard to the scoring function and to the optimization algorithm, with emphasis on a faster execution. The SARS-CoV-2 Mpro structure and two of its main conformers, extracted from the molecular dynamics simulation trajectory file, were processed with AutoDockTools (Morris et al., 2009). The small-molecule database consisted of 4384 molecules and was downloaded from ZINC15 (Sterling and Irwin, 2015). These molecules corresponded to the “world” subset; according to ZINC15, these molecules have been approved for human use in major jurisdictions, including the U.S. Food and Drug Administration. The database was processed with the prepare_ligand4.py script included in AutoDockTools to obtain the input files needed for virtual screening. Box coordinates for virtual screening were center (x, y, z) = 93.3, 10.5, 2.0 and size (x, y, z) = 22.3, 22.1, 24.2. The candidate molecules were selected according to the docking score predicted by iDock; the top 10 molecules from each conformer were further analyzed to select the potential inhibitors of Mpro.

Molecular dynamics simulations of Mpro-ligand complexes

The 30 complexes were inspected visually with Pymol (https://pymol.org) and Discovery Studio Visualizer (https://www.3dsbiovia.com); those molecules that made contact with the catalytic residues of Mpro (His41 and Cys145) were subjected to molecular dynamics simulations. The input files were prepared with CHARMM-GUI as mentioned above, with simulation times of 20 ns. The resulting trajectories were visualized and analyzed with VMD and GROMACS. The stability of the ligands in the active site of Mpro was estimated by calculating the RMSD (Root-Mean-Square Deviation) and the binding energy (ΔG) was calculated with the CaFE plugin (Liu and Hou, 2016) in VMD. The CaFE plugin consists of a set of Tcl scripts that allow the calculation of binding free energies through the MM/PBSA and LIE methods and is powered by VMD, NAMD (Phillips et al., 2005) and APBS (Jurrus et al., 2018) programs. Binding free calculations were performed with the MM/PBSA method. Briefly, receptor and ligand conformations are extracted from the trajectories, then three energetic components are calculated: The gas-phase energy difference between the complex, the receptor, and the ligand is obtained by calling NAMD. Then the polar solvation free energy is calculated with APBS. Subsequently, the difference of solvent accessible surface area (SASA) is measured and the nonpolar solvation free energy is estimated by its approximate linear relation with SASA. Finally, the binding free energy is summed and averaged throughout the selected conformations.

Results

Homology modeling of Mpro and molecular dynamics simulations

There are several experimental structures reported for the SARS-CoV-1 Mpro, which vary in amino acid sequence, resolution, and presence of additional molecules, such as inhibitors or solvents. Consequently, we decided to include structures 2AMD, 1WOF, 2AMQ, 2D2D, 3E91, and 3EA7 to model the structure of the SARS-CoV-2 main protease. Sequence identity of SARS-CoV-2 Mpro with SARS-CoV-1 sequences ranged from 94.4 % (sequence corresponding to PDB code 2AMQ) to 95.1 % (PDB codes 2AMD and 1WOF). The resulting model is in the classical heart-shaped conformation (Fig. 1 ), with the binding site near the dimer interface. The selected model had an RMSD=0.266 Å when compared with template 2AMD.

Fig. 1

Modeled structure of the SARS-CoV-2 Mpro dimeric protease. One monomer is shown in magenta and the other in cyan. Inhibitor 9IN, from structure 2AMD, is shown as spheres. Image generated with Pymol.

The final Mpro model was analyzed with the SWISS-MODEL server, resulting in a global QMEAN score of -0.41 with only minor deviations in the N- and C-termini (Fig. S1). These values indicate a good structural model, suitable for use in the following analyses. To observe the amino acid conservation of Mpro, amino acid sequences from 111 reported isolates were aligned and the degree of conservation was mapped into the Mpro structure. The most variable position had a 97.7 % conservation, which indicates that the active site of Mpro is highly conserved among the analyzed isolates (Fig. S2).

Next, solvent-explicit molecular dynamics simulations were performed on Mpro; the resulting trajectory showed that the protein has a highly flexible active site as the amino acids surrounding the binding site had high RMSF (Root-Mean-Square Fluctuation) values. The regions surrounding the binding site were the most mobile during the simulation (Figs. 2 and S3). This flexibility makes drug discovery a challenge because a ligand that could bind to the closed site could not bind to the open state and vice versa; therefore, we proceeded to analyze the MD trajectory, and from 19 clusters that encompass the trajectory, the two most representative conformers were selected (Fig. S4). Thus, the virtual screening was done with the starting structure, as well as two conformers: Mpro_0 (initial structure), Mpro_412 (conformer 412, which is representative of cluster 1), and Mpro_837 (conformer 837, which is representative of cluster 2).

Fig. 2

Flexibility of Mpro as determined from molecular dynamics simulations. The image on the left shows Mpro rendered as a transparent surface and the backbone as “putty”. Regions with high RMSF values are shown as thicker; colors range from blue (most stable) to red (most flexible). The highlighted area is zoomed in on the right. Ligand IN9 is displayed as “sticks” and is shown to identify the active site.

Virtual screening of approved drugs

Virtual screening was performed with the Mpro conformers and the 4384-molecule bank. The top 10 binders for each conformer were selected for further analysis (Table 1 ). Docking score values ranged from −9.17 kcal/mol (dihydroergotamine) to −10.25 kcal/mol (bisoctrizole). Noticeably, there are several drug metabolites (glucuronides) and drugs used in chemotherapy (sorafenib beta-d-glucoronide, N-Trifluoroacetyladriamycin, amrubicin, daunorubicin, idarubicin, midostaurin, maraviroc, and celsentri). Ligands that contacted catalytic amino acids His41 or Cys145 were the most promising candidates and thus were selected for further molecular dynamics simulations.

Table 1

Top 10 molecules with the best docking score.

Ligand Name	ZINC ID	iDock score (kcal/mol)
Mpro_0
Bisoctrizole	ZINC000011677911	−10.25
Ironophore	ZINC000017545546	−10.08
Sorafenib Beta-d-Glucuronidea	ZINC000077313075	−10.08
Oxazepam Glucuronide	ZINC000031460595	−10.04
Ergoloidb	ZINC000003995616	−9.95
Eltrombopag	ZINC000011679756	−9.85
Telmisartan	ZINC000001530886	−9.74
Bromocriptine	ZINC000053683151	−9.73
Dutasteride	ZINC000003932831	−9.72
Ergotamine	ZINC000052955754	−9.67
Mpro_412
N-Trifluoroacetyladriamycinc	ZINC000095618916	−9.96
Amrubicin	ZINC000003780800	−9.77
Daunorubicin	ZINC000003917708	−9.75
5a Reductase inhibitor	ZINC000014880001	−9.43
Carindacillin	ZINC000003871889	−9.42
Simeprevir	ZINC000253632968	−9.41
Idarubicin	ZINC000003920266	−9.38
Midostaurin	ZINC000100013130	−9.38
Lorazepam Glucuronide	ZINC000031290884	−9.33
Meclocycline	ZINC000084402690	−9.30
Mpro_837
Maraviroc	ZINC000101160855	−9.38
Glycyrrhizinate Dipotassium	ZINC000096015174	−9.33
Teniposide	ZINC000004099009	−9.32
Indomethacin Glucuronide	ZINC000084386263	−9.32
Celsentri	ZINC000003817234	−9.31
Lorazepam Glucuronide	ZINC000031290884	−9.26
Estrone Glucuronide	ZINC000004099004	−9.23
5a Reductase inhibitor	ZINC000014880001	−9.22
Ketotifen-N-Glucuronide	ZINC000095618575	−9.20
Dihydroergotamine	ZINC000003978005	−9.17

Glucuronides are products of metabolization mainly by the liver.

Names in bold correspond to molecules that contacted either catalytic amino acid (His41 or Cys145).

A metabolite of Valrubicin.

MD simulations for the selected Mpro-ligand complexes were carried out and the relative mobility of the ligand within the active site of Mpro was monitored through evaluation of ligand RMSD (Fig. 3 ). In the case of Mpro_0, the most mobile ligand was ergoloid, with an average RMSD of 7.6 Å, and the most stable ligand was ergotamine, with an average RMSD of 3.6 Å. In the case of Mpro_412, the most mobile ligand was daunorubicin with an average RMSD of 7.2 Å, and the most stable ligand was lorazepam glucoronide, with an average RMSD of 3.4 Å. Lastly, for Mpro_837 the most mobile ligand was lorazepam glucoronide with an average RMSD of 7.8 Å, and the most stable ligand was ketotifen-N-glucoronide, with an average RMSD of 3.5 Å.

Fig. 3

Ligand RMSD within the binding site of Mpro. a) Ligands complexed with Mpro_0; b) ligands complexed with Mpro_412; and c) ligands complexed with Mpro_837. Vertical axis units are Å and horizontal axis units are picoseconds.

Binding free energy was calculated for each complex (the initial 5 ns were not considered for calculation) and the results are indicated in Table 2 . It is noteworthy that complexes with conformer Mpro_0 and Mpro_412 had the lowest binding free energy values. After considering all the bioinformatics analyses, we posit that ligands daunorubicin, ergotamine, bromocriptine, meclocycline, amrubicin, ergoloid, ketotifen-N-glucoronide, N-trifluoroacetyladriamycin, 5a reductase inhibitor are good candidates as Mpro inhibitors and therefore warrant further evaluation as potential treatment against SARS-CoV-2 infection.

Table 2

Binding free energies of Mpro-ligand complexes.

Conformer	ZINC code	Drug	ΔG (kcal/mol)
Mpro_412	ZINC000003917708	Daunorubicin	−138.8
Mpro_0	ZINC000052955754	Ergotamine	−119.2
Mpro_412	ZINC000003780800	Amrubicin	−117.5
Mpro_0	ZINC000053683151	Bromocriptine	−116.7
Mpro_412	ZINC000084402690	Meclocycline	−115.1
Mpro_0	ZINC000003995616	Ergoloid	−109.1
Mpro_837	ZINC000095618575	Ketotifen-N-Glucuronide	−89.0
Mpro_412	ZINC000095618916	N-Trifluoroacetyladriamycin	−46.2
Mpro_412	ZINC000014880001	5α Reductase Inhibitor	−31.0
Mpro_412	ZINC000031290884	Lorazepam Glucuronide	−1.1
Mpro_837	ZINC000031290884	Lorazepam Glucuronide	−0.5
Mpro_837	ZINC000084386263	Indomethacin Glucuronide	11.9
Mpro_837	ZINC000096015174	Glycyrrhizinate Dipotassium	36.1

Discussion

A race against time is currently underway to develop safe and effective antivirals against SARS-CoV-2. Drug repurposing is a straightforward strategy when time is a critical factor, such as in this COVID-19 pandemic (Phadke and Saunik, 2020; Pushpakom et al., 2019). Indeed, the U.S. Food and Drug Administration has encouraged the use of the App CURE ID for health care professionals to report novel uses of existing medicines for the treatment of infectious diseases (https://cure.ncats.io/), and there are several studies, some in clinical phases, that involve the repurposing of approved drugs (Pacios et al., 2020). On the other hand, bioinformatic studies are useful in detecting drugs that can be used in in vitro studies or clinical trials as treatment options against COVID-19. In this work, a set of nine currently approved drugs was identified as having the potential to inhibit the SARS-CoV-2 Mpro and therefore could be used as a COVID-19 treatment, either alone or in combination with other drugs.

Some of the most promising strategies involve the use of chloroquine alone or combined with azithromycin, as well as remdesivir (Gautret et al., 2020; Wang et al., 2020a). However, the efficacy of these therapies is as yet unknown. Regarding Mpro, currently, there is at least one clinical trial involving the use of protease inhibitors as a treatment for COVID-19 (clinical trial NCT04303299) and there are several reports of potential inhibitors of SARS-CoV-2 Mpro (Chen et al., 2020b; Khan et al., 2020; Li and De Clercq, 2020; Ton et al., 2020; Zhang et al., 2020). These reports use computational approaches to identify molecules with Mpro inhibitory potential; however, our work is the only one that considers the flexibility of the active site, which allowed us to propose a greater diversity of potential Mpro inhibitors. Mpro is an excellent target for computer-assisted drug discovery because it is critical in the early stages of viral replication (Chen et al., 2020a); however, we found that the substrate-binding site is highly flexible, therefore virtual-assisted drug design approaches should take flexibility into account. When this work started, there were no structures available for SARS-CoV-2 Mpro; however, two experimental structures have been reported to date (PDB id: 6LU7 and 6W63) (Jin et al., 2020). A rapid analysis showed that the b-factors of 6LU7 agree with our results (Supplementary Fig. S5), therefore the homology model constructed in this work is valid for our purposes.

MD simulations had a duration of 20 ns, and the first 5 ns were discarded for binding-free energy analysis to allow the system to equilibrate. There is no consensus about the duration of a molecular dynamics simulation, but it has been found that longer simulations do not necessarily improve binding free energy calculations (Hou et al., 2011). We chose 20 ns as a result of a literature search, in which simulation times ranged from 15 ns (one protein-ligand complex) to 200 ns (several protein-ligand complexes) (Kumar et al., 2019a; Mittal et al., 2020). To define the duration of a simulation run, the complexity of the phenomena under study and the size of the system should be considered, also depending directly on the available computer power. Overall, studies involving one complex normally used longer simulation times (Kumar et al., 2019b; Wang et al., 2020b) and studies involving several complexes used shorter simulation times (Ge et al., 2013; Razzaghi-Asl et al., 2018).

Binding free energies ranged from −138.8 kcal/mol to 36. 1 kcal/mol. At first sight, these values are noticeably higher than expected, considering that the experimental ΔG value for the avidin-biotin complex is −20.4 kcal/mol. The reason for this discrepancy is the solute dielectric variable used for our calculations; in fact, the use of different solute dielectric values can significatively shift the binding free energy calculation to a more negative or positive value (Li et al., 2018). In the case of the avidin-biotin complex, a solute dielectric variable of 1.0 agrees well with experimental values since calculations made with the CaFE plugin indicate a ΔG = −20.8 kcal/mol when a solute dielectric = 1.0 was used; this value is close to the experimental binding free energy (-20.4 kcal/mol) (Green, 1975). The choice of the solute dielectric value depends on the charges present at the binding interface and increases as more charges are present. Usually, values of 1.0, 2.0 and 4.0 are used depending on the nature of the binding site (Hou et al., 2011). More often, values ranging from 2.0 to 4.0 are used when studying binding sites with polar or charged amino acids are present. For proteases similar to Mpro, the solute dielectric value that correlates best with experimental results is 4.0 (Sun et al., 2014).

Three of the proposed drugs are used in chemotherapy (daunorubicin, amrubicin, and the valrubicin metabolite N-trifluoroacetyladriamycin) (Piska et al., 2017). These drugs as well as meclocycline, a highly toxic antibiotic, should be treated carefully because of side effects in the organism. Two of the drugs identified in this study, indomethacin and glycyrrhizinate, have also been tested in vitro and proposed as potential therapies against coronavirus (Amici et al., 2006; Hoever et al., 2005). However, the binding free energy analysis suggests that Mpro might not be the primary target of these molecules.

Despite the increasing number of publications related to virtual screening focused on Mpro (Jin et al., 2020; Ton et al., 2020; Wang, 2020), there is little agreement between them on potential candidates identified. Of all the reported drugs found in this work, only valrubicin also was identified in another study as a potential inhibitor of Mpro (Wang, 2020). This may be explained by the diversity of the selected databases and the algorithms used in docking and virtual screening. However, if the original purpose of the drug is considered, similarities arise between our findings and those of Jin et al. (2020) where drugs used in chemotherapy and antihistamines potentially inhibit Mpro. Given the emergency we are currently facing, this diversity of results may turn out to be beneficial because more options are available to repurpose these drugs as COVID-19 treatment options. Again, the drugs mentioned in this article should not be used as treatment against COVID-19 unless they have been tested in proper clinical trials.

Conclusions

The main protease of SARS-CoV-2 has a very flexible active site, an aspect that must be considered in rational drug design. Our approach involved the use of three conformers of Mpro and led to the identification of nine molecules that warrant in vivo testing or even in clinical trials so they can be repurposed as treatment of COVID-19.

CRediT authorship contribution statement

Alicia Jiménez-Alberto: Conceptualization, Investigation, Writing - review & editing, Funding acquisition, Visualization. Rosa María Ribas-Aparicio: Conceptualization, Resources, Writing - review & editing, Funding acquisition, Supervision, Formal analysis. Gerardo Aparicio-Ozores: Conceptualization, Investigation, Writing - review & editing, Funding acquisition, Data curation. Juan A. Castelán-Vega: Methodology, Software, Validation, Formal analysis, Data curation, Writing - original draft, Visualization, Funding acquisition.

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.