key: cord-0963183-ig6n9cet
authors: Jiménez-Alberto, Alicia; Ribas-Aparicio, Rosa María; Aparicio-Ozores, Gerardo; Castelán-Vega, Juan A.
title: Virtual screening of approved drugs as potential SARS-CoV-2 main protease inhibitors
date: 2020-06-25
journal: Comput Biol Chem
DOI: 10.1016/j.compbiolchem.2020.107325
sha: e1db12b3919465c1c77de6f7863a9e9441e14400
doc_id: 963183
cord_uid: ig6n9cet

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.

In December 2019, a cluster of atypical pneumonia cases was reported in Wuhan, China (Wu et al., 2020; Zhu et al., 2020) . Coronavirus disease rapidly spread to other countries and by March 11 th , 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 15 th , 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 (Y. . 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 (Y. 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 . The protease J o u r n a l P r e -p r o o f 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 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 .

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.

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, J o u r n a l P r e -p r o o f 1WOF, 2AMQ, 2D2D, 3E91, and 3EA7 . 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) .

The predicted structural model was submitted to the CHARMM-GUI server to prepare the system (Brooks et al., 2009; Jo et al., 2014 Jo et al., , 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 (Abraham et al., , 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 (Danura 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.

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.

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 

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 97.7% conservation, which indicates that the active site of Mpro is highly conserved among the analyzed isolates ( Figure 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 (Figures 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 ( Figure 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). 

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 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. 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 ( Figure 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 Å. 

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. J o u r n a l P r e -p r o o f

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; . 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 (Y. W. 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 (Y. ; 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 Figure 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 J o u r n a l P r e -p r o o f 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) (A. 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 (N. ) and studies involving several complexes used shorter simulation times (Ge, 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 . 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) . 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; J. Wang, 2020) , there is little agreement between them on potential J o u r n a l P r e -p r o o f 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 (J. 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.

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.

Below you will find the CRediT author statement for each participant in the elaboration of this manuscript:

Alicia Jiménez-Alberto: Conceptualization, investigation, writing-review & editing, funding acquisition, visualization; Rosa M. 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, writingoriginal draft, visualization, funding acquisition.

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. 

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

The Gromacs development team

Indomethacin has a potent antiviral activity against SARS coronavirus

A bibliometric review of drug repurposing

Structural insights into SARS coronavirus proteins

Toward the estimation of the absolute quality of individual protein structure models

CHARMM: The biomolecular simulation program

Drug repurposing from the perspective of pharmaceutical companies

Emerging coronaviruses: Genome structure, replication, and pathogenesis

Prediction of the SARS-CoV-2 (2019-nCoV) 3C-like protease (3CL pro) structure: virtual screening reveals velpatasvir, ledipasvir, and other drug repurposing candidates

WHO Declares COVID-19 a Pandemic

Peptide Folding: When Simulation Meets Experiment

Identification of Combinations of Approved Drugs With Synergistic Activity Against Ebola Virus in Cell Cultures

Hydroxychloroquine and azithromycin as a J o u r n a l P r e -p r o o f treatment of COVID-19: results of an open-label non-randomized clinical trial

Molecular Dynamics-Based Virtual Screening: Accelerating the Drug Discovery Process by High-Performance Computing

Avidin

Critical Assessment of the Important Residues Involved in the Dimerization and Catalysis of MERS Coronavirus Main Protease

Antiviral Activity of Glycyrrhizic Acid Derivatives against SARS−Coronavirus

Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations

CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins

VMD: Visual molecular dynamics

First respiratory transmitted food borne outbreak?

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

Chapter Eight -CHARMM-GUI PDB Manipulator for Advanced Modeling and Simulations of Proteins Containing Nonstandard Residues

CHARMM-GUI: A web-based graphical user interface for CHARMM

Improvements to the APBS biomolecular solvation software suite

Identification of Chymotrypsin-like Protease Inhibitors of SARS-CoV-2 Via Integrated Computational Approach

E-pharmacophore modelling, virtual screening, molecular dynamics simulations and in-silico ADME analysis for identification of potential E6

Structure-based virtual screening, molecular dynamics simulation and MM-PBSA toward identifying the inhibitors for twocomponent regulatory system protein NarL of Mycobacterium Tuberculosis

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges

Identifying SARS-CoV-2 related coronaviruses in Malayan pangolins

CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field

Therapeutic options for the 2019 novel coronavirus (2019-nCoV)

idock: A multithreaded virtual screening tool for flexible ligand docking

The impact of interior dielectric constant and entropic change on HIV-1 complex binding free energy prediction

CaFE: a tool for binding affinity prediction using end-point free energy methods

Identification of potential molecules against COVID-19 main protease through structure-guided virtual screening approach

AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility

Unipro UGENE: a unified bioinformatics toolkit

Strategies to Combat Multidrug-Resistant and Persistent Infectious Diseases

UCSF Chimera--a visualization system for exploratory research and analysis

COVID-19 treatment by repurposing drugs until the vaccine is in sight

Scalable molecular dynamics with NAMD

Metabolic carbonyl reduction of anthracyclines -role in cardiotoxicity and cancer resistance. Reducing enzymes as putative targets for novel cardioprotective and chemosensitizing agents

Drug repurposing: progress, challenges and recommendations

Identification of COX-2 inhibitors via structure-based virtual screening and molecular dynamics simulation

GISAID: Global initiative on sharing all influenza data -from vision to reality

Clustal Omega, accurate alignment of very large numbers of sequences

ZINC 15 -Ligand Discovery for Everyone

Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set

Rapid Identification of Potential Inhibitors of SARS-CoV-2 Main Protease by Deep Docking of 1.3 Billion Compounds

Roles of Host Gene and Non-coding RNA Expression in Virus Infection

Fast Identification of Possible Drug Treatment of Coronavirus Disease -19 (COVID-19) Through Computational Drug Repurposing Study

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Discovery of the Novel Inhibitor Against New Delhi Metallo-β-Lactamase Based on Virtual Screening and Molecular Modelling

SWISS-MODEL: homology modelling of protein structures and complexes

Comparative Protein Structure Modeling Using MODELLER

Novel Coronavirus (2019-nCoV) situation reports [WWW Document

The SARS-CoV-2 outbreak: what we know

Activation and maturation of SARS-CoV main protease

Drug design targeting the main protease, the Achilles' heel of coronaviruses

Design of wide-spectrum inhibitors targeting coronavirus main proteases

Vgas: A Viral Genome Annotation System

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors

The novel coronavirus outbreak in Wuhan