Screening of Clinically Approved and Investigation Drugs as Potential Inhibitors of SARS-CoV-2: A Combined in silico and in vitro Study.

In the current study, we used 7922 FDA approved small molecule drugs as well as compounds in clinical investigation from NIH's NPC database in our drug repurposing study. SARS-CoV-2 main protease as well as Spike protein/ACE2 targets were used in virtual screening and top-100 compounds from each docking simulations were considered initially in short molecular dynamics (MD) simulations and their average binding energies were calculated by MM/GBSA method. Promising hit compounds selected based on average MM/GBSA scores were then used in long MD simulations. Based on these numerical calculations following compounds were found as hit inhibitors for the SARS-CoV-2 main protease: Pinokalant, terlakiren, ritonavir, cefotiam, telinavir, rotigaptide, and cefpiramide. In addition, following 3 compounds were identified as inhibitors for Spike/ACE2: Denopamine, bometolol, and rotigaptide. In order to verify the predictions of in silico analyses, 4 compounds (ritonavir, rotigaptide, cefotiam, and cefpiramide) for the main protease and 2 compounds (rotigaptide and denopamine) for the Spike/ACE2 interactions were tested by in vitro experiments. While the concentration-dependent inhibition of the ritonavir, rotigaptide, and cefotiam was observed for the main protease; denopamine was effective at the inhibition of Spike/ACE2 binding.

Introduction

Coronaviruses (CoVs) are the family of viruses containing single‐stranded RNA (positive‐sense) which is encapsulated by a membrane envelope. They are classified in the Nidovirales order, Coronaviridae family, which is comprised of two sub‐families and about 40 known species. These species are divided and characterized into four genera (alpha, beta, gamma and delta), and only the alpha and beta‐ strains are identified to be pathogenic to human and other mammals.[ , ] Before 2019, six coronaviruses were known to cause respiratory and enteric diseases in humans, especially the two of them belonging to betaviruses cause severe illness: SARS (Severe Acute Respiratory Syndrome)‐CoV and MERS (Middle East Respiratory Syndrome)‐CoV. A novel coronavirus is discovered in Wuhan, China in late 2019, and officially named as SARS‐CoV‐2 (formerly 2019‐nCoV) due to its genomic similarity to SARS‐CoV.[ , , , ] The disease caused by this virus is officially named as Coronavirus Disease 2019 (COVID‐19) by World Health Organization (WHO). Like SARS‐ and MERS‐CoVs, SARS‐CoV‐2 mostly affects the lower respiratory tract to cause pneumonia, and may also affect the gastrointestinal system, kidney, heart and central nervous system, with the common symptoms including fever, cough and diarrhea. On 11th of March 2020, WHO declared the COVID19 as pandemic. The first emergence of the virus was witnessed at the penultimate days of 2019 as pneumonia concentrated in Wuhan, China. The outbreak in China was then spread very quickly to the other countries and as of June 2021, more than 177 million individuals have been infected by SARS‐CoV‐2 virus and it is expected that the number will increase in the following months. Thus, drugs and vaccines are highly in demand to control the outbreak.

The genomic sequence of SARS‐CoV‐2 is available (GenBank ID: MN908947.3) and the initial analyses indicate that different essential enzymes found in SARS‐CoV and SARS‐CoV‐2 have the sequence similarity of around 80 %. Furthermore, the catalytic sites of the four key enzymes that could be the antiviral targets are vastly conserved between the two coronaviruses. SARS‐CoV‐2 is also reported to utilize the same cell‐entry receptor for infection, ACE2 (Angiotensin‐Converting Enzyme 2), as SARS‐CoV.[ , ] The genome of SARS‐CoV‐2 encodes for various proteins and important ones are 3‐chymotrypsin‐like protease, 3CLpro aka Main Protease, papain‐like protease, helicase, and RNA‐dependent RNA polymerase which construct the non‐structural proteins and Spike glycoproteins which belong to structural proteins.[ , ] While the Spike protein is the key for the virus to enter the cell through the interaction with ACE2, enzymes such as Main Protease are crucial for the life cycle of the virus. These proteins are the most attractive targets for the development of new drugs against SARS‐CoV‐2 due to their pivotal role in entry into the host cell and replication and transcription of the virus. One of the advantages of targeting these proteins is that although the mutagenesis rate is high in viruses, not many mutations observed in these proteins, except the Spike RBD, since any mutation here can be lethal for the virus.

There are more than 100 3D protein structures for SARS‐CoV‐2, mostly for the Main Protease structure in apo‐ and holo‐ states, resolved via X‐ray diffraction or cryo‐electron microscopy, deposited and available in Protein Data Bank (RCSB PDB). In silico studies related to Main Protease are increasing in terms of protein‐inhibitor interactions and drug screening.[ , ] For a broad review of additional background, patents, developments and perspectives in COVID‐19 and other diseases related to coronavirus, the reader is referred to see Cao et al.

The phenomenon known as drug repositioning or repurposing has gained attention as the development of new drug starting from the beginning becoming more costly in respect to both time and resources required.[ , , , ] Established favorable toxicological, pharmacokinetic and pharmacodynamic properties of approved/clinical drug molecules make them suitable to be used for new indications.[ , , , ] As the molecules considered in repurposing studies passed through several stages and have well‐defined profiles, they would not require prolonged pre‐clinical studies and hence, they would be excellent candidates in the cases of disease emergencies or outbreaks.[ , ] Drug repurposing studies have already been conducted for various kind of diseases (review articles for different diseases[ , , , , ]). Thanks to repositioning studies, different compounds have been found new usages than their original purposes[ , ] even though they have been failed in their original purpose and/or withdrawn from the market. Computational approaches such as virtual screening would decrease the time required for the identification of new targets for the existing drug molecules with the advantage of also being cost‐efficient as demonstrated in review and research articles.[ , , , ]

In our research group, virtual screening of different ligand databases including FDA approved drugs in recent years have been performed by an in‐house script and it is shown that the obtained results by this screening algorithm, which is a hybrid algorithm of ligand‐ and target‐driven based screening techniques, gave successful results.[ , , , , ] Thus, in the current study, this hybrid algorithm is applied for the identification of clinically approved or investigational compounds against SARS‐CoV‐2 essential target proteins (i. e., Main Protease and Spike receptor‐binding domain bound with ACE2).

Main protease has been studied by different research groups to find inhibitors capable of halting the activity of this enzyme and consequently the reproduction of the virus. After the SARS outbreak in 2003 many researches were conducted to target the Main Protease of the SARS‐CoV. Chen et al. utilized screening approaches using a 3D model of SARS‐CoV 3CLpro to screen the MDL‐CMC database that contains 8.000 compounds. Cinanserin was among the high‐ranked final compounds for which in vitro studies were performed. The proteolytic activity of the enzyme was shown to be inhibited by 70 to 90 % at 50 to 100 μM of Cinanserin. In a study performed by Liu et al. homology modeling was used to construct a model of main protease since SARS 3CLPro was not publicly available at the time of the work. Then, high throughput virtual screening was performed using different chemical libraries including the National Cancer Institute (NCI) Diversity Set (230.000 compounds total), ACD‐3D (Available Chemical Database, Release: ACD 3D 2002.2, 280.000 compounds in total), and MDDR‐3D (MDL Drug Database Report, Release: MDDR 3D 2002.2, 120.000 compounds in total). The final hits (40 compounds) were further tested in vitro to check the inhibition activity, and 3 of them were found to inhibit the protease activity up to 40 %. C3930 or calmidazolium which is the calmodulin antagonist was found as the best hit with the highest inhibition activity.

After the recent outbreak of SARS‐CoV‐2, many research groups have started to use screening methods to search for the inhibitors of Main Protease and Spike protein domain/ACE2 complex. In a recent paper, Li et al. have screened 8000 molecules including the approved or experimental compounds and small molecules derived from DrugBank. The protease protein with the PDB ID 5N5O was used as target. Compounds showing better docking score than −7.7 kcal/mol were selected as hits and experimentally unapproved ones, as well as those with strong side effects were removed from the list. The list was even shortened considering the marketability of the molecules. Prulifloxacin, Bictegravir, Nelfinavir and Tegobuvi are finally selected molecules. In a research conducted by Chen et al. apo‐enzyme structure of SARS‐CoV (PDB ID: 2DUC) was used to build a model for the main protease of SARS‐CoV‐2 and MTiOpenScreen web service was used to screen for purchasable drugs (Drugs‐lib). The library has 7173 compounds. Autodock Vina was used to screen the active site at chains A and B and finally 10 and 11 drugs were selected for these chains respectively based on the energy cut‐off. Another study which is conducted by Jin et al. targeted the Main Protease as well. This research group found some promising compounds by combining structure‐based drug design approaches with screening methods. In vitro cell‐based assays showed the high inhibitory effect of the chosen final compounds on the target enzyme and antiviral activities. Virtual screening studies were performed by Jin et al. using a model constructed based on the crystal structure of SARS‐CoV‐2 Main Protease in complex with N3 inhibitor (PDB ID: 6LU7). They used an in‐house library containing potential binding compounds. Cinanserin was found as the best binding affinity to the substrate‐binding pocket of the enzyme, and in vitro studies showed an IC50 value of 124.93 μM for Cinanserin. Cinanserin, a serotonin antagonist was previously found to also inhibit SARS‐CoV. Fluorescence Resonance Energy Transfer (FRET)‐based high‐throughput screening resulted in the finding of some FDA‐approved drugs (Disulfiram and Carmofur) and other compounds which are in preclinical/clinical trial (Ebselen, TDZD‐8, Shikonin, Tideglusib, and PX‐12). In their study, the compound Ebselen was found the most active compound with IC50 value of 0.48 μM.

Results and Discussion

In the present study, we used 7922 compounds from NIH Chemical Genomics Center (NCGC) Pharmaceutical Collection (NPC) database (https://tripod.nih.gov/npc/) and in order to eliminate the non‐specific binders, some criteria including molecular weight, between 100 to 1000 g/mol; number of rotatable bonds, <100; number of atoms, between 10 and 100; number of aliphatic and aromatic rings, <10; number of hydrogen‐bond acceptor and donors, <10 were set and as a result the total number of compounds was decreased to 6654. These 6654 compounds were then docked to the binding cavities of apo (PDB, 6M03) and holo (PDB, 6LU7) forms of SARS‐CoV‐2 main protease enzyme (Figure 1). Since the binding pocket conformations in apo and holo forms are slightly different, and both apo and holo forms are available at the RCSB database we decided to use both forms in our virtual screening protocol. The 6M0J PDB‐coded structure was used for the Spike protein/ACE‐2 target. In docking, standard precision (SP) protocol of Glide docking module of Schrodinger software was used. Tables S1 and S2 show the top‐100 docking scored compounds based on the docking scores at the main protease in holo and apo forms, respectively. Table S3 shows the corresponding docking scores for Spike/ACE‐2 complex.

Figure 1

Around 7000 FDA approved and drugs in clinical investigation from NPC database were screened at the apo (top‐left) and holo (bottom‐left) COVID‐19 Main Protease target. Hierarchical hybrid screening constructed by our group led to 6 hit compounds. Surface representation of one the identified hit compound Rotigaptide is shown in right.

Although recent studies have suggested that docking is a successful approach for selecting hits, since in the docking approach flexibility of both protein residues and docked ligand are not fully considered, hence the ranking and ordering of the compounds only by their corresponding docking scores may not potentially lead to the identification of correct compounds. Moreover, although molecular docking studies may give an initial insight into protein‐ligand interactions, it is always crucial to understand the maintenance of these interactions by performing dynamical studies such as molecular dynamics (MD) simulations. Therefore, we selected top‐100 compounds based on docking scores from each docking simulations and initially performed short (10‐ns) MD simulations for these complexes (in total 3‐μs MD simulations). An in‐house script was used for the preparation of simulation boxes as well as for the analysis of MD simulations. Desmond was used for MD simulations. Tables S4 and S5 represent average MM/GBSA scores using collected 1000 trajectory frames (strided by 10 throughout each simulations) of selected 100 top‐docking scored compounds from both holo‐ and apo‐based simulations, respectively. Table S6 shows average MM/GBSA scores of selected top‐100 docking scored compounds at Spike Protein/ACE‐2 interface. We also performed MD simulations for the co‐crystallized ligand‐bound structure (inhibitor N3 at PDB ID: 6LU7) using the same MD protocol for compound screening. Figure S1 shows protein‐ligand interaction diagram of inhibitor N3 at the Main Protease. The figure includes a timeline representation of the interactions and contacts (H‐bonds, hydrophobic, ionic, water bridges) representing which residues interact with the ligand in each trajectory frame. Interactions that occur more than 15 % of the simulation time in the selected trajectory (0 through 500 ns), are shown. The stacked bar charts are also normalized over the course of the trajectory (i. e., a value of 0.5 suggests that 50 % of the simulation time the specific interaction is maintained). Results showed that the following residues are crucial for ligand binding: Thr25, Thr26, His41, Ser46, Asn142, Met165, Glu166, and Gln189. Several water bridges and hydrogen bonding interactions dominate the interaction constructed from Glu166. The interactions between these residues and the screened compounds were also checked. Average MM/GBSA scores of co‐crystallized ligand N3 from long (500‐ns) MD simulations was found as −67.2±9.7 kcal/mol. Thus, we forwarded compounds that have better average MM/GBSA scores than a cutoff value (−70.0 kcal/mol) from short MD simulations to long MD simulations. We selected top‐21 compounds based on average MM/GBSA scores from short (10‐ns) MD simulations performed on holo form main protease target and conducted 100‐ns long MD simulations for each (Table S7). Other selected compounds from short MD simulations such as Arzoxifene and Truxicurium could not maintain their initial crucial interactions during the simulations. Based on average MM/GBSA scores from 100‐ns simulations, we also forwarded 13 compounds to 500‐ns MD simulations (Table S8).

Interestingly, while 21 compounds identified using the target retrieved from holo‐state of main protease short simulations fit the cutoff value, only 5 compounds were found from database screened at apo‐state. Long (100‐ns) MD simulations were also performed for compounds screened in apo form Main Protease enzyme and only two compounds were found as promising ligands (Table S9).

While Table 1 shows identified compounds targeting Main Protease enzyme, Table 2 represents determined compounds at the Spike protein domain/ACE2 interface. As can be seen from Table 1, for Main Protease enzyme average MM/GBSA scores show that following 11 compounds: Pinokalant, Bms181176‐14, Terlakiren, Bisnafide, Ritonavir, Cefotiam, Telinavir, Rotigaptide, Cefpiramide, Pimelautide, and Lopinavir have significant average MM/GBSA scores, and throughout the MD simulations the interaction between these compounds and the crucial residues of the target were maintained (Figures S2–S12). MM/GBSA scores of tested compounds at the main protease and at the Spike/ACE2 domain were also represented as box and Whisker plots at Figures 2 and 3, respectively. In the plots, while central lines in boxes represent corresponding median values of MM/GBSA scores of investigated compounds, bottom lines of the boxes show first quartile (25 %), and top line of the boxes show third quartile (75 %). The plots at Figure 2 show that within 11 identified hit compounds at the Main Protease, top lines of boxes of only two of them (Pimelautide and Lopinavir) have smaller (in absolute values) scores than cutoff value −70.0 kcal/mol. Although average MM/GBSA scores of Pimelautide and Lopinavir have similar values with positive control N3, both of them have smaller (in absolute values) average scores than −70.0 kcal/mol. The corresponding values of top lines of boxes of all of the other identified 9 hit compounds have better MM/GBSA scores than cutoff value. Figure 3 shows corresponding plot for the identified hits at Spike/ACE2 domain. While the top lines of two of the identified molecules (Benzquercin and Naminterol) have smaller scores (in absolute values) than −70.0 kcal/mol, other three hits (Bometolol, Denopamine, and Rotigaptide) have better scores. Figures S13–S16 represent ligand interactions diagrams for the identified compounds at Spike Protein/ACE2 interface. Following 3 compounds were identified as hit ligands at the Spike Protein/ACE2 domain: Bometolol, Denopamine and Rotigaptide.

Table 1

Selected hit compounds based on average MM/GBSA scores at Main Protease target. Long (100‐ns and 500‐ns) MD simulations are performed for these identified hits and average MM/GBSA scores were calculated using 1000‐trajectory frames (2000‐trajectory frames for 500‐ns) throughout the simulations. Table also shows the docking scores and corresponding mechanism of actions of the identified compounds.

Compounds	2D Structures	Docking Score (kcal/mol)	MM/GBSA (kcal/mol)	Mechanism of Action
Pinokalant		−8.2	−87.9±6.8	Broad‐spectrum cation channel blocker
Bms181176‐14 (Becatecarin)		−8.3	−87.6±7.7	Becatecarin is a small molecule, anticancer compound for the treatment of hepatobiliary duct tumors
Terlakiren		−8.1	−86.9±8.7	Antihypertensive, Renin inhibitor
Bisnafide		−8.3	−84.6±6.8	Bisnafide is a bis‐naphthalimide compound with anticancer activity.
Ritonavir		−8.2	−83.9±9.4	Ritonavir is an antiretroviral protease inhibitor
Cefotiam		−7.3	−82.2±7.7	Cefotiam is a parenteral second‐ generation cephalosporin antibiotic
Telinavir		−7.0	−81.1±8.7	Telinavir is an anti‐HIV aspartyl protease inhibitor
Rotigaptide		−9.0	−78.5±10.4	Rotigaptide (ZP‐123) is a drug under clinical investigation for the treatment of cardiac arrhythmias – specifically atrial fibrillation.
Cefpiramide		−7.2	−78.4±5.6	Cefpiramide is a third‐generation cephalosporin antibiotic.
Pimelautide		−9.0	−65.9±11.1	Immunostimulant. Pimelautide Built‐in Adjuvants Associated with an HIV‐1‐Derived Peptide
Lopinavir		−8.1	−63.1±12.1	Lopinavir is an antiretroviral of the protease inhibitor

Table 2

Selected hit compounds based on average MM/GBSA scores at ACE‐2/Spike Protein domain. Long (100‐ns) MD simulations are performed for these identified hits and average MM/GBSA scores were calculated using 1000‐trajectory frames throughout the simulations. Table also shows the docking scores and corresponding mechanism of actions of the identified compounds.

Compounds	2D Structures	Docking Score (kcal/mol)	MM/GBSA (kcal/mol)	Mechanism of Action
Bometolol		−7.4	−81.8±7.9	β‐adrenergic blocking agent
Denopamine		−6.9	−79.9±4.7	Denopamine (INN) is a cardiotonic drug which acts as a β1 adrenergic receptor agonist.
Rotigaptide		−7.0	−76.2±7.7	Rotigaptide (ZP‐123) is a drug under clinical investigation for the treatment of cardiac arrhythmias – specifically atrial fibrillation.
Benzquercin		−6.5	−71.8±7.4	A flavonoid compound
Naminterol		−6.7	−69.5±4.7	Naminterol is a β2 adrenoceptor agonist with bronchodilatory properties.

Figure 2

Box and Whisker plot representation of MM/GBSA scores of selected hit compounds at the SARS‐CoV‐2 Main Protease target. In the plots, central line in box represents median, bottom line of box is first quartile (25 %), top line of box is third quartile (75 %), bottom of whiskers is first quarter minus 1.5× interquartile range, top of whiskers is third quarter plus 1.5× interquartile range, and dots are outliers.

Figure 3

Box and Whisker plot representation of MM/GBSA scores of selected hit compounds at the SARS‐CoV‐2 Spike/ACE2 target. In the plots, central line in box represrents median, bottom line of box is first quartile (25 %), top line of box is third quartile (75 %), bottom of whiskers is first quarter minus 1.5× interquartile range, top of whiskers is third quarter plus 1.5× interquartile range, and dots are outliers.

Interestingly, Rotigaptide was identified as a promising compound both in Main Protease and Spike protein/ACE2 targets. Identified compound Rotigaptide is a drug under clinical investigation for the treatment of cardiac arrhythmias – specifically atrial fibrillation. Crucial residue interactions were formed by Thr26, His41, Glu166, Gln189 and Gln192 at the Main Protease binding pocket (Figure S9). While Cefotiam and Cefpiramide are antibiotics, Ritonavir, Lopinavir and Telinavir have been used in clinical studies for HIV infection. Terlakiren and Pinokalant are other identified compounds at the Main Protease which have anti‐hypertensive and broad‐spectrum cation channel blocker profile, respectively. Two of the identified molecules at the Main Protease (Bms181176‐14 and Bisnafide) were used in clinical studies as anti‐cancer compounds. Thus, usage of these anti‐cancer compounds in COVID‐19 may not be suggested.

We also docking studies for known Main Protease inhibitors40 from literature using same protocol. Then, top‐docking poses of these positive controls were used in 100‐ns MD simulations with the same MD protocol employed for the screened compounds and average MM/GBSA free energy scores are compared. Following average MM/GBSA energy scores were obtained for Ebselen, Disulfiram, Tideglusib, Carmofur, and Shikonin: −44.8, −44.5, −47.5, −46.2, −47.7 kcal/mol, respectively. Corresponding average MM/GBSA scores of identified compounds at the Main Protease were similar or even better than these compounds.

Two of the identified promising compounds (i. e., Bometolol and Denopamine) at the Spike/ACE2 domain are compounds targeting beta‐adrenergic receptors. As it is mentioned, Rotigaptide is another compound which strongly binds both Spike Protein/ACE2 and Main Protease targets.

Since several research papers highlight the importance of glycans that are playing a role not only in shielding, but in the functional dynamics of the system itself, we repeated simulations by using glycosylated SARS‐CoV‐2 Spike protein for the selected ligands at the Spike/ACE2 site.[ , ] For this aim, we considered two approaches: (i) Glycosylated SARS‐CoV‐2 Spike protein (PDB, 6VSB) was aligned with the ligand‐bound Spike/ACE2 complexes (PDB, 6M0J) from our docking simulations, and ligands poses were merged with the 6VSB‐PDB coded target protein that includes glycans (Positions of glycans were used from CHARMM‐GUI Archive – COVID‐19 Proteins Library). By this way, glycosylated SARS‐CoV‐2 Spike/ACE2 ligand‐bound ternary structures were prepared; (ii) Glycosylated SARS‐CoV‐2 Spike/ACE2 apo protein was used as target protein and proposed ligands at the Spike/ACE2 were docked using same previous docking protocols at the current study. In the simulations, one N‐glycan in RBD‐receptor‐binding‐domain (Asn343), and four N‐linked glycans in ACE2 (Asn53, Asn90, Asn322, and Asn546) were used. Obtained ligand‐bound Spike/ACE2 targets including glycans were used as input files at the MD simulations. Same MD simulations protocols and post‐processing MD analyses utilized for screening of identified compounds were performed. When we used approach‐I, following average MM/GBSA scores were obtained: Denopamine, −77.7 kcal/mol; Rotigaptide, −70.1 kcal/mol; and Bometolol, −65.8 kcal/mol. The results showed that glycans did not change the mean MM/GBSA score of Denopamine, while the mean MM/GBSA scores of Bometolol and Rotigaptide decreased slightly (at absolute values) (Figure S17). When we used approach‐II, we obtained results similar to approach‐I, but the decrease in the mean MM/GBSA scores of Rotigaptide and Bometolol was more remarkable (Figure S18). The binding interface of Spike/ACE2 region is very large thus compounds may form slightly different initial binding poses. Such as Figure S19 shows top‐docking poses of Rotigaptide at the binding pocket of Spike/ACE2 interface with and without glycans. As it is shown, ligand binds similar region at the interface and it has common binding residues in glycosylated and non‐glycosylated cases, however docking poses are slightly different (Figures S19 and S20). Although similar key residues form non‐glycosylated and glycosylated Spike /ACE2 interfaces with ligands tested, ligand‐protein interactions diagrams show that ligand‐protein contacts are formed and break faster in the latter case (Figure S21). In order to compare how the binding pocket volumes are changing throughout the simulations in case of glycosylated and non‐glycosylated cases, binding pocket volumes are calculated and compared throughout the simulations. Results showed that average binding pocket volumes, especially in the cases of Bometolol and Rotigaptide, are smaller when glycans were used in the simulations (Figure S22).

In order to verify the predicted screening results of integrated in silico analyses, selected hit compounds are tested using in vitro enzyme assays. Thus, four compounds (i. e., ritonavir, rotigaptide, cefotiam, and cefpiramide) for the main protease and two compounds (i. e., denopamine and rotigaptide) for the Spike/ACE2 interactions were ordered and tested by in vitro experiments. Two different inhibitor screening assay kits against Spike protein were used as ACE2:SARS‐CoV‐2 Spike S1 (#79945) and SARS‐CoV‐2 Spike:ACE2 (Cat. #79931), (BPS Bioscience). Both assays were designed for screening and profiling inhibitors of the Spike and ACE2 interactions, and both studied in 96‐well format. Denopamine was tested with SARS‐CoV‐2 Spike: ACE2 (Cat. #79931) while Rotigaptide was tested both inhibitor assay kits with different concentrations. The inhibition was measured via chemiluminescence, that can be measured using a chemiluminescence reader. 3CL Protease Assay Kit (#79955‐1 and #79955‐2, BPS Bioscience, San Diego CA, USA) was used to screen the activity of main protease.

Figure 4 represents the biochemical assay results of tested hit compounds ritonavir, cefotiam, and rotigaptide against the main protease. All three compounds showed their concentration‐dependent effects on the enzyme activity. Enzyme activity of the main protease was decreased around 26 % with 200 μM of ritonavir. When concentrations of these three ligands are increased, their effect on enzyme activity becomes clearer. However, cefpiramide did not significantly inhibited enzyme activity at any of the concentrations tested (Figure S23). Figure 5 shows the effect of Denopamine to the Spike/ACE2 binding. Denopamine at 200 μM concentration reduced the Spike/ACE2 binding effectively by 84 %. Figure S24 shows the corresponding results for the rotigaptide. While it decreases the Spike/ACE2 binding around 83 % at 100 μM concentrations, surprisingly higher concentrations (i. e., 200 μM) did not improve the decrease in binding percentage (Figure S24).

Figure 4

Change of main protease enzyme activity percentage with different inhibitor concentrations. All molecules were tested at least with four replicates. Mean and standard deviation results were provided at the figure.

Figure 5

Change in SARS‐CoV‐2 Spike/ACE2 binding with different inhibitor concentrations. The compound was tested at least with four replicates. Mean and standard deviation results were provided at the figure.

Experimental and theoretical studies overlap well with each other for the studied compounds. In both of in vitro and in silico studies Ritonavir was the most effective one within the studied four hits at the main protease. Moreover, binding free energy calculations showed that Denopamine is very effective for the inhibition of Spike/ACE2 binding and it is verified by the in vitro tests.

Conclusions

In this virtual drug repurposing study, we used 7922 FDA approved drugs and compounds in clinical investigation from NPC database. Both apo and holo states of SARS‐CoV‐2 Main Protease as well as Spike Protein/ACE2 were used for virtual screening. Initially, docking was performed for these compounds at target binding sites. The compounds were then sorted according to their docking scores which represent binding energies. The first 100 compounds from each docking simulations were initially subjected to short (10‐ns) MD simulations (in total 300 ligand‐bound complexes), and average binding energies during MD simulations were calculated using the MM/GBSA method. Then, we performed up to 500‐ns MD simulations for each system and for the free energy calculations of the screened ligands, we performed MM/GBSA. Similar length of MD simulations have been employed for the MM/GBSA calculations at the literature for the similar system sizes. However, in order to see the effect of the simulation time, we repeated all simulations with 1‐μs of MD simulations for each system. Results showed that all the studied compounds still show effective binding profiles at the target proteins (Figures S25 and S26). Here, a total of more than 25‐μs simulations were conducted at different lengths (short (i. e., 10‐ns) and long (100 and 500‐ns). Both docking and MD simulations followed by consecutive binding free energy calculations resulted that holo form of the target protein is more appropriate choice for virtual drug screening studies. As shown in Figures S27 and 28, high MM/GBSA scores for the selected hits are sustained throughout the simulations. These numerical calculations provided following 7 compounds as hit compounds for SARS‐CoV‐2 Main Protease: Pinokalant, Terlakiren, Ritonavir, Cefotiam, Telinavir, Rotigaptide, and Cefpiramide. In addition, following 3 compounds were identified as potential hits for SARS‐CoV‐2 ACE‐2/Spike protein: Bometolol, Denopamine, and Rotigaptide. The effect of the glycans to the binding profiles of the identified compounds at the Spike/ACE‐2 interface was also investigated. The results revealed that when using glycans in the simulations, the change in mean MM/GBSA score for Denopamine was very small, and that in Rotigaptide and Bometolol, although the binding energy values were still promising, they affected the MM/GBSA scores. Four compounds from the proposed compounds are tested by in vitro assays against main protease and two compounds were tested against Spike/ACE2 binding. Results showed that Ritonavir, Cefotiam, and Rotigaptide compounds are found effective at high concentrations at the main protease. Denopamine was effectively (84 %) blocked the interactions of Spike/ACE2 at 200 μM concentration. Overall, our results here showed that the computational simulations can effectively guide the drug repurposing studies, especially during the pandemic period.

Methods

7922 compounds were downloaded from NPC database (https://tripod.nih.gov/npc/). In order to eliminate the non‐specific binders, some criteria including molecular weight, between 100 to 1000 g/mol; number of rotatable bonds, <100; number of atoms, between 10 and 100; number of aliphatic and aromatic rings, <10; number of hydrogen‐bond acceptor and donors, <10 were set and as a result the total number of compounds was decreased to 6654. These ligands were prepared using LigPrep module of Maestro at neutral pH (LigPrep, Schrodinger v.2017). OPLS3 force field is used. In molecular docking, we used following protein structures: apo‐Main Protease, (PDB, 6M03); holo‐Main Protease, (PDB, 6LU7); and Spike Protein/ACE‐2, (PDB, 6M0J). Binding site of the main protease was defined by centering grids at the centroid of a set of three crucial residues in ligand binding, namely His41, Cys145, and Glu166. Ali and Vijayan stated a very strong and sustained salt bridge interactions between Lys417 of SARS‐CoV‐2 Spike RBD and Asp30 of ACE‐2. Thus, the corresponding residues at the Spike/ACE‐2 were used in grid generation. In the Spike/ACE2 target protein (PDB, 6M0J), glycans were not included at the simulations. However, in order to see the effect of glycans to the binding profiles of the proposed ligands at the Spike/ACE2 interface, simulations were repeated for the selected ligands at the glycosylated SARS‐CoV2 Spike protein (PDB, 6VSB). Open state conformation of glycosylated SARS‐CoV2 Spike protein was retrieved from CHARMM‐GUI Archive – COVID‐19 Proteins Library (http://www.charmm‐gui.org/?doc=archive&lib=covid19) as template and aligned with 6M0J‐coded structure, thus Spike/ACE2 target including glycans was generated (one N‐glycan in RBD (Asn343), and four N‐linked glycans in ACE2 (Asn53, Asn90, Asn322, and Asn546) were used). These proteins were prepared using Protein Preparation module of Maestro. PROPKA was used for determination of protonation states of amino acid residues. Restrained minimization was performed with OPLS3 force field for the protein using 0.3 Å heavy atom convergence. Docking was performed with Glide/SP using default settings. Protein‐ligand complexes were placed in the orthorhombic boxes with explicit TIP3P water models that have 10 Å thickness from the edges of target proteins. All systems were neutralized and 0.15 M NaCl salt solution added to the systems. The long‐range electrostatic interactions were calculated by the particle mesh Ewald method. A cut‐off radius of 9 Å was used for both van der Waals and Coulombic interactions. Simulations were performed at body temperature (310 K) and 1.01325 bar. Nose‐Hoover thermostat and Martyna‐Tobias‐Klein barostat was used at the simulations. The time step was 2 fs. The OPLS3 force field was used in simulations. Throughout the MD simulations, 1000 trajectory frames (for 10‐ns and 100‐ns simulations) and 2000 trajectory frames (for 500‐ns simulations) were recorded and Molecular Mechanics Generalized Born Surface Area (MM/GBSA) binding free energies of compounds were calculated. VSGB 2.0 solvation model at Prime module of Maestro was utilized during MM/GBSA calculations. In this method, the free energy of the binding ( ) of the small molecule (L) to a biological macromolecule (P) to form their complex (PL) can be represented as: . This equation can be decomposed as: , where ; ; ( . In the equations, and terms are the changes in the molecular mechanics (MM) energy in gas phase, solvation free energy and entropy upon ligand binding, respectively.48 Summation of electrostatic energies ( , and van der Waals energies ( ) gives changes in has polar ( and nonpolar ( contributions between the solute and the continuum (implicit) solvent, in which polar contribution can be calculated by generalized Born (GB) model, and nonpolar energy contribution can be estimated using solvent accessible surface area (SASA). It must be noted that the GB method gives an analytical expression for the polar solvation energy. Since the GB is much faster approach than the Poisson‐Boltzmann (PB) method, GBSA was considered in this study for the fast estimation of binding free energies of screened selected hit compounds. The MM/GBSA method is an effective procedure to improve the predictions of docking methods which was used for calculating the binding free energies of studied compounds because of its computational efficiency as well as its better performance for ranking molecules according to their binding affinities. The entropic contributions at MM/GBSA energies are not included in calculations. This is also one the reasons for highly negative values of free energies calculated by MM/GBSA algorithm. As mentioned in the review paper of Genheden and Ryde, entropy calculations could be costly in virtual screening of large databases and in most cases would not be accurate enough to improve the results. As such this term is generally omitted especially at the screening of large sets of compounds. SiteMap module is used to calculate binding pocket volume changes throughout MD simulations. SiteMap is a tool for the determination and evaluation of protein binding sites, yet this tool is also used for calculating the MD simulations trajectory binding site volumes. The atoms within 10 Å of the ligand are included for calculation.

Two different inhibitor screening assay kits against Spike protein were used as ACE2:SARS‐CoV‐2 Spike S1 (#79945) and SARS‐CoV‐2 Spike:ACE2 (Cat. #79931) Inhibitor Screening Assay (BPS Bioscience). Both assays were designed for screening and profiling inhibitors of the Spike and ACE2 interactions, and both studied in 96‐well format. Denopamine was tested with SARS‐CoV‐2 Spike:ACE2 (Cat. #79931) while Rotigaptide was tested both inhibitor assay kits at different concentrations. The inhibition was measured via chemiluminescence that can be measured using a chemiluminescence reader. SARS‐CoV‐2 Spike:ACE2 (Cat. #79931) aimed that the detection of His‐labeled ACE2 by HRP‐labeled Anti‐His. Protocol was followed as suggested by provider (BPS Bioscience). Briefly, 96‐well plate was coated with SARS‐CoV‐2 Spike protein and the plate was incubated overnight at 4 °C. After washing the plate with assay buffer and blocking with blocking buffer, the inhibitor molecules were added to Spike protein coated well plate with different concentration and then incubated at room temperature with slow shaking. Potential inhibitor molecules were tested in 500 μM, 200 μM and 100 μM. Next, ACE2‐His was incubated at room temperature with inhibitor molecules and SARS‐CoV‐2 proteins on the plate. As final step, chemiluminescence was produced by adding Anti His‐HRP and then HRP substrate and measured by using HIDEX Sense Microplate reader. Apart from previous assay kit, the ACE2: SARS‐CoV‐2 Spike S1 (#79945) kit aims the high sensitivity of detection of Spike S1‐Biotin protein by Streptavidin‐HRP. Protocol was followed as suggested by manufacturer (BPS Bioscience). As first step, ACE2 protein was added in order to attached the nickel‐coated 96‐well plate. After one‐hour incubation at room temperature, plate was washed with assay buffer; followed by adding inhibitor molecules to plate covered ACE2 protein with different concentration. Here, rotigaptide was used at 100 μM, 200 μM and 500 μM concentrations. Then, SARS‐CoV‐2 Spike S1‐Biotin was added to the plate with ACE2 and incubated at four hours at room temperature with slow shaking. Finally, Streptavidin‐HRP and HRP substrate were added respectively to produce chemiluminescence. The product was measured using HIDEX Sense Microplate Reader with 1 second of integration time without filter when measuring light emission as measuring condition.

3CL Protease Assay Kit (#79955‐1 and #79955‐2, BPS Bioscience, San Diego CA, USA) was used to screen the activity of main protease. 100 mM stock solutions of the compounds were prepared by 100 % DMSO. Different concentrations of each compound ranging between 100 to 500 μM were tested. Final DMSO concentrations were kept below 1 % for each condition. 5 ng of 3CL Protease enzyme was distributed to each well except “Blanks”. The kit comes with its own inhibitor (GC376) as an inhibitor control. 500 μM GC376 was added to the wells designated as inhibitor control. 5 μl of inhibitor compounds was added to their relative wells, and a 1X assay buffer/DMSO mixture was added to blanks and positive controls. 250 μM 3CL Protease substrate was added to each well to start the reaction, and its final concentration was 50 μM in 25 μl volume. After 4 hours incubation at room temperature, fluorescence was measured by a microtiter plate‐reader at a wavelength of 360 nm for excitation and 460 nm for emission. Blank values are subtracted from values of all other wells. The GraphPad Prism 8.0 software (GraphPad, San Diego CA, USA) was used to calculate activity percentage of each concentration of each compound, relative to no inhibitor control that was considered as 100 % inhibition. All molecules were tested at least four replicates. Mean and standard deviation results were provided.

Author Contributions

SD conceived, designed and directed the conducted study as principal investigator. SD, BA, and KS, and BD performed the simulations. SD, TA, BA, BD, KS, AS, and NBI analyzed and interpreted the data. MDO, SC, and TA performed and analyzed all the experiments. SD wrote the manuscript. All authors read the manuscript and provided comments and revisions to both scientific content and grammatical corrections. All authors read and approved the final form of the manuscript.

Conflict of interest

None declared.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.