Drug Repurposing to Identify Nilotinib as a Potential SARS-CoV-2 Main Protease Inhibitor: Insights from a Computational and In Vitro Study.

COVID-19, an acute viral pneumonia, has emerged as a devastating pandemic. Drug repurposing allows researchers to find different indications of FDA-approved or investigational drugs. In this current study, a sequence of pharmacophore and molecular modeling-based screening against COVID-19 Mpro (PDB: 6LU7) suggested a subset of drugs, from the Drug Bank database, which may have antiviral activity. A total of 44 out of 8823 of the most promising virtual hits from the Drug Bank were subjected to molecular dynamics simulation experiments to explore the strength of their interactions with the SARS-CoV-2 Mpro active site. MD findings point toward three drugs (DB04020, DB12411, and DB11779) with very low relative free energies for SARS-CoV-2 Mpro with interactions at His41 and Met49. MD simulations identified an additional interaction with Glu166, which enhanced the binding affinity significantly. Therefore, Glu166 could be an interesting target for structure-based drug design. Quantitative structural-activity relationship analysis was performed on the 44 most promising hits from molecular docking-based virtual screening. Partial least square regression accurately predicted the values of independent drug candidates' binding energy with impressively high accuracy. Finally, the EC50 and CC50 of 10 drug candidates were measured against SARS-CoV-2 in cell culture. Nilotinib and bemcentinib had EC50 values of 2.6 and 1.1 μM, respectively. In summary, the results of our computer-aided drug design provide a roadmap for rational drug design of Mpro inhibitors and the discovery of certified medications as COVID-19 antiviral therapeutics.

Introduction

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is accountable for the corona virus disease 2019 (COVID-19), a devastating pandemic with more than four million dead as of August 2021.[1−4] This novel RNA virus of the genus Betacoronavirus (BCoV) and family Coronaviridae(5) is an encircled positive-sense RNA virus with significantly bigger genomes with respect to the known RNA viruses.[5,6] SARS-CoV and MERS-CoV are two more members of this class.[5,6] Unlike the other BCoV members, SARS CoV-2 developed with a mutation that allowed the receptor binding domain of its spike glycoprotein to connect to the angiotensin-converting enzyme-2 with greater affinity.[7−10] A second mutation between the S1 and S2 glycoproteins created a novel polybasic furin cleavage site, amplifying the pathogen’s potency many times. Despite rigorous efforts of diagnostic testing, isolation, contact tracing, and immunization efforts, the number of infected and the rate of mortality are on a surge.[11−13] Rapid development of drugs and/or vaccines has, thus, been sought after, leading to intensive research to save the global economy from a near-collapse situation.[5,11,14−16]

There has been extensive research on drug repurposing to find FDA-approved drug candidates as COVID-19 antiviral therapeutics.[17] Drug repurposing has been emphasized as a priority given the prolonged time span required for the experimental validation of safety and efficacy with new potential drugs, thus saving a substantial amount of time and money.[18] Thus, the majority of the drugs under investigation are repurposed from FDA-approved therapeutics[19,20] Current repurposed drugs, that are either being prescribed or further developed for the COVID-19 treatment, include ritonavir, lopinavir, oseltamivir, remdesivir, ribavirin, and interferon-α.[21−25] Recent reports exhibited ritonavir and remdesivir to be efficacious against the SARS-CoV-2 in vitro and in vivo with the need of further experimental validations.[26,27] However, recent clinical trials have failed to produce very strong evidence for the efficacy of ritonavir and remdesivir when compared to placebos.[28,29] Ritonavir, if considered, should be administered well before the peak of the viral replication to observe an optimal effect.[30] A recent clinical trial with 1063 patients shows that remdesivir shortened the recovery time by four days with respect to the placebo.[28] Additionally, a significant difference was reported in the mortality rate between the drug (7.1%) and placebo (11.9%).[28] Thus, there has been constant efforts in identifying effective drug candidates with wider therapeutic indices against SARS-CoV-2. Recently, there are numerous reports on the computer-aided drug design (CADD) approach, aiming to detect repurposed drugs with potentially enhanced efficacy against various drug targets of the COVID-19 pathogen.[14,15,31−35]

Our goal was to find either FDA-certified or advanced clinical trial candidates, aiming at the SS7 main protease (Mpro). Similar to the previous coronaviruses, the duplication of SARS-CoV-2 relies on the production and proteolytic processing of two overlapping polyproteins, leading to multiple functional subunits.[36] Mpro or 3C-like protease (3CLpro) is responsible for site-specific hydrolysis of polyproteins. Thus, Mpro is one of the most crucial enzymes for virus replication and considered one of the most attractive drug targets among coronaviruses. One of the biggest advantages of the Mpro target is that human and SARS-CoV-2 proteases differ in cleavage specificity.[37] Hence, it is unlikely that an Mpro inhibitor will have side effects with a human protease. Recently, Jin and co-workers have reported the co-crystal structure of Mpro with a potent antagonist (6LU7, Figure ).[38] Two co-crystal structures of Mpro with sulfonamide-based antagonists were reported by Fearon et al. (5R70 and 5R7Y, Figure ).[39,40] There are numerous reports on drug repurposing virtual screening studies employing different methodologies against Mpro, resulting in several drug-like molecules with potential Mpro inhibitory activity.[15,16,35,41−46]

Figure 1

Recently discovered SARS-CoV-2 Mpro co-crystal structures. (A) SARS-CoV-2 Mpro complex with native ligand N3 (PDB: 6LU7). (B) SARS-CoV-2 Mpro complex with native ligand Z18197050 (PDB: 5R80). (C) COVID-19 Mpro complex with native ligand Z45617795 (PDB: 5R7Y).

In this current study (Figure ), a structure-based (SB) pharmacophore model was developed to carry out virtual screening experiments on the Drug Bank library, consisting of 8823 active pharmaceuticals either approved or under investigation. Pharmacophore-centered screening led to 502 candidates with potential to bind to the active site of Mpro in a noncovalent manner. The candidates from the pharmacophore screening were subjected to additional screening employing molecular docking to identify potential noncovalent inhibitors. The molecular docking study led to 44 promising compounds with binding affinities of −8.0 kcal/mol or above. Selected poses of top 44 candidates underwent molecular dynamics simulations (MD simulation) to assess the stability of the protein ligand complexes. The MD simulation study led to six top hits with potential Mpro inhibitory activity. A quantitative structure–activity relationship (QSAR) study was conducted on the 44 most promising docking hits. Prescreening of thousands of potential drug candidates based on their pharmacological potency in rational drug design and discovery is labor-intensive and time-consuming in the pharmaceutical industry. The use of CADD, computational chemistry, and molecular docking in conjunction with QSARs and multivariate analysis such as principal component analysis (PCA) and partial least square (PLS) regression has proven critical and their applications in rational drug design.[47−49] Accordingly, multivariate regression analysis of the QSAR of drug candidates was used for pattern recognition and to model and predict the binding energies of individual agents. Finally, 10 promising commercially available compounds with very good binding affinity and interactions at His41 and Met49 for the SARS-CoV-2 Mpro active site were selected for an in vitro SARS-CoV-2 cytopathic effect experiment.

Figure 2

(A) Workflow to identify repurposed therapeutics as potential SARS-COV-2 Mpro inhibitors. (B) Top 10 drugs biologically evaluated for SARS-COV-2 replication inhibition.

Methods

Computational Details

We used a series of computational and bioinformatics tools to predict the best drug candidate. We started with SB pharmacophore hypothesis generation, and virtual screening was performed using LigandScout 4.4/Advance software. Molecular docking experiments were carried out employing PyRx/AutoDock Vina software. MD simulation was accomplished using the Amber 18 MD package. Visualization of molecular interactions was done using Discovery Studio Visualizer and VMD. QSAR was carried out applying Unscrambler Chemometrics software (CAMO Software, 9.8, Oslo, Norway).

Chemicals and Cells

The American Type Culture Collection provided Vero E6 cells, which had been cultured in minimal essential medium enriched with 100 μg/mL streptomycin, 100 U/mL penicillin, and 5% fetal bovine serum. Chemicals were received from the MedChem Express Library (MCE; Monmouth Junction, NJ, USA); and the Cayman Chemical Library (Ann Arbor, MI, USA). The purity of the chemicals ranged from 96 to 99 percent. Chemical dilutions and dosage response concentrations were taken care of manually.

Generation of the SB Merged Pharmacophore

Three co-crystal structures of SARS-CoV-2 Mpro with the respective noncovalent inhibitors, namely, N3 (PDB ID: 6LU7, resolution 2.16 Å), Z18197050 (PDB ID: 5R80, resolution 1.93 Å), and Z45617795 (PDB ID: 5R7Y, resolution 1.65 Å), were used to generate three individual energy-optimized SB pharmacophore hypotheses using LigandScout’s SB module. Three structural coordinates of the Mpro–ligand complexes were downloaded from the PDB. The water molecules within 5 Å of the native ligand were eliminated during refinement of protein–ligand complexes. Protein structure adjustments, including addition of missing hydrogens and loops, were performed followed by controlled energy minimization using the GROMOS96 43B1 force field.[50] These steps were performed using Swiss-pdb viewer 4.1.0, AutoDockTools-1.5.7, and Discovery Studio Visualizer software. Prepared protein–ligand complexes were imported to LigandScout’s SB module workspace to map default pharmacophore features, including (a) hydrophobicity, (b) hydrogen bond donor, and (c) hydrogen bond acceptor. The alignment module of LigandScout was then employed to align three SB pharmacophores with respect to the 6LU7 SB hypothesis as the reference point. The alignment module was further utilized to generate the final merged pharmacophore hypothesis in two steps, namely, (a) generate a shared feature pharmacophore and then (b) merge pharmacophores and interpolate overlapping features.

Virtual Screening using the Merged Pharmacophore Model

Chemical structures of 8823 drugs candidates/ligands from the DrugBank database underwent virtual screening using the developed merged pharmacophore. The energy of the ligands was minimized employing the MMFF94 force field using the SB module of LigandScout 4.4/Advance software before undergoing screening. Using LigandScout’s screening module, probable stereoisomers and tautomers were created, and protonation states were designated. Finally, virtual screening was accomplished using LigandScout’s screening module to identify the drug candidates with important pharmacophore features responsible for improved binding to Mpro. A pharmacophore fit score of 40% or above was considered to select best hits (Table S1).

Virtual Screening on Molecular Docking

Initial screening using pharmacophores yielded 502 best hits with probability to interact strongly with the Mpro catalytic site. This subset of drug candidates was then subjected to molecular docking against the 6LU7 Mpro–ligand complex to screen the most promising candidates. Molecular docking-based virtual screening was performed employing PyRX/AutoDock Vina software. A grid box was created by considering the native ligand (N3) at the center of it. It was made sure that all nine hotspot residues are inside the grid box. The dimensions of the box were fixed at 23 × 28 × 24 Å. Docking score was applied to rank the candidates based on their binding affinities. Candidates with a docking score of −8 kcal/mol or above were considered as hits, leading to 44 most promising drug candidates. All cysteine residues kept their natural (protonated) form.

Validation of the Docking Process

A control experiment was set up to validate the docking process. The native ligands in the crystal structures 6LU7 (N3), 6LZE (FHR), 6M0K (FJC), 6WTK (UED), and 6XHL (V2M) underwent docking against the refined protein while maintaining the uniform grid box. The docking score of this control experiment was set as a benchmark for comparing scores of other drugs. This experiment was accomplished using AutoDock Vina.

MD Simulations

System Preparation and Setup

The initial coordinate of the protein was derived from the recently crystallized COVID-19 Mpro complex consisting of an antagonist N3 (PDB code: 6LU7, solved at a resolution of 2.16 Å). The drug was docked into the 6LU7 receptor using AutoDock Vina. Lack of the hydrogen atoms in the protein was fulfilled by the LEAP module of the Amber18 package. Atomic partial charges and the absent parameter for the drug molecule were acquired by a semiempirical method (AM1-BCC) using the antechamber module of Amber18. The generalized Amber force field was employed for all the drug molecules. After that, the subsequent system undergoes solvation in a rectangular box of transferable intermolecular potential 3-point water model stretching up to a distance of 8 Å from the surface of the protein. After checking the total charge of the system, to neutralize it, few Na+ ions had been introduced into the surface of the protein. Each MD simulation utilized the Amber ff14SB force field for the protein residues.

MD Simulations

After proper parametrization, minimization was performed on the system’s geometries to eliminate all the nasty contacts and to relax the system (5000 steps by the steepest descent and 5000 steps by the conjugate gradient approach). Afterward, to adjust the temperature, the system was gently annealed for 50 ps, under the NVT ensemble, where a weak restraint was applied on protein residues. Subsequently, the system undergoes density equilibration for 1 ns under the NPT ensemble (constant pressure of 1.0 atm and constant temperature of 300 K). To maintain the constant temperature and pressure, a Langevin thermostat with 2 ps collision frequency and Barendsen barostat with 1 ps pressure relaxation time were used, respectively. This density equilibration over 1 ns is a MD simulation performed with weak restraints under periodic boundary conditions until the system obtained a uniform density. After 1 ns of density equilibration, all the restraints applied before (during heating and density equilibration) had been eliminated, and the system was further equilibrated for 3 ns followed by 100 ns of MD production run. Using the SHAKE algorithm, all the covalent bonds consisting of hydrogen had been kept constrained during all MD simulations. Particle mesh Ewald was utilized to manage long-range electrostatic interactions. The GPU version of the Amber18 package had been used for MD simulation of each complex, and the CPPTRAJ module was used for the analyses of the trajectories.

MMPBSA Calculation

The molecular mechanics generalized Boltzmann surface area (MMGBSA) calculations had been conducted over the 2000 frames selected from the trajectory, in which root mean square deviation (rmsd) of the complex was found to be well converged. Water molecules and Na+ ions present in the system had been removed before the MMGBSA analysis. For solute, the dielectric constant was 1, and for the surrounding solvent molecules, it had been kept at 80. MMGBSA calculation was carried out using the MMPBSA.py module of Amber18. We used the Onufriev and David Case GB model for MMGBSA calculations.

QSAR Analysis

The QSAR dataset used for the PCA and PLS regression is made up of 44 drug candidates (Table ). Drug candidates in Table were divided into two datasets, training set and validation dataset. The training set comprised 23 drug candidates. However, 19 drug candidates were used as validation or independent samples (Table S2). The drug candidate’s QSAR including the molecular weight, number of each atom present in the molecule, presence of single, double, and triple bonds, and benzene ring (x-variables) was modeled with subjecting to PCA and modeled by PLS regression for the prediction of drug candidate binding energies (x-variables).[51,52] Unscrambler Chemometrics software (CAMO Software, 9.8, Oslo, Norway) was used to conduct the drug candidates’ QSAR regression analysis.

Table 1

Drug Candidates with Highest Docking Scores (−kcal/mol) toward Mproa

drug name	Drug Bank ID	docking score (kcal/mol)	pharmacophore fit score	primary target
bemcentinib	DB12411	–10.4	45.7	AXL kinase inhibitor
nilotinib	DB04868	–9.2	55.8	tyrosine kinase (TK) inhibitor
SAR-125844	DB15382	–8.9	54.6	MET TK inhibitor
temafloxacin	DB01405	–8.8	65.7	fluoroquinolone antibiotic
enasidenib	DB13874	–8.7	55.7	IDH2 inhibitor
tegobuvir	DB11852	–8.6	55.5	HCV, chronic
NA	DB08553	–8.5	67.0	ser/thr-protein kinase B-raf
ziresovir	DB15145	–8.4	55.6	RSV
NA	DB08141	–8.4	65.1	cyclin dependent kinase-2
danoprevir	DB11779	–8.4	43.2	NS3/4A protease inhibitor
selinexor	DB11942	–8.4	42.8	refractory multiple myeloma
filibuvir	DB11878	–8.3	42.2	NS5B inhibitor for HCV
NA	DB02573	–8.3	55.3	ribonuclease pancreatic
A-443654	DB08073	–8.3	55.2	kinase inhibitor
danicopan	DB15401	–8.3	54.0	factor D inhibitor
AGG-523	DB15460	–8.3	52.3	aggrecanase-selective inhibitor
SB-649868	DB14822	–8.2	56.7	dual orexin receptor antagonist
NA	DB08166	–8.2	56.0	ser/thr-protein kinase pim-1
lufenuron	DB11424	–8.2	55.1	flea control
defibrotide sodium	DB04932	–8.2	42.2	antithrombotic, anti-inflammatory
nafamostat	DB12598	–8.2	42.0	synthetic protease inhibitor
indinavir	DB00224	–8.1	42.0	protease inhibitor for HIV
sisunatovir	DB15674	–8.1	66.0	RSV fusion inhibitor
montelukast	DB00471	–8.1	41.8	anti-inflammatory
NA	DB04020	–8.1	65.0	estrogen-β receptor
SB220025	DB04338	–8.1	57.1	mitogen-activated protein kinase 1/14
gemigliptin	DB12412	–8.1	56.5	DPP4
KW-7158	DB05498	–8.1	56.0	urinary incontinence
carfecillin	DB13506	–8.1	56.0	β-lactam antibiotic
NA	DB08445	–8.1	55.7	DPP4
brivanib	DB11865	–8.1	55.5	hepato-cellular carcinoma
lopinavir	DB01601	–8.0	41.6	antiretroviral protease inhibitor
NA	DB08588	–8.0	68.2	DPP4
verucerfont	DB12512	–8.0	56.4	CRF-1 antagonist
NA	DB06994	–8.0	56.3	DPP4
BMS-690514	DB11665	–8.0	56.3	inhibitor of EGFR and VEGFR
retagliptin	DB14898	–8.0	56.0	DPP4 inhibitor
CHEBI:40083	DB07145	–8.0	56.0	disintegrin and metalloproteinase domain-containing protein 17
salvianolic acid	DB15246	–8.0	55.5	oxidative stress
TAK733	DB12241	–8.0	54.8	inhibitor of MEK1
flufenoxuron	DB15006	–8.0	54.0	insecticide
NA	DB08590	–8.0	42.0	NA
NA	DB02510	–8.0	41.6	NA
NA	DB07568	–8.0	42.3	NA

NA: commercial name not available.

Biological Evaluation

In Vitro Screening

The in vitro screening experiment was performed following the protocol recently published by Jonsson and co-workers.[53]

Results and Discussion

Virtual Screening on the SB Pharmacophore

A pharmacophore represents an assembly of steric and electrostatic features of different molecules that are crucial to ensure best-possible supramolecular interactions with a given protein structure to stimulate or deactivate its pharmacological response.[26] Recently, pharmacophore-based screening has been employed extensively as a result of its capacity to identify potential hits from a large dataset within a short timeframe and with minimum computational powers.[54] In the present study, we have used the “SB” pharmacophore creation guideline available in LigandScout to generate pharmacophores indicative of nonbonding interactions of three individual Mpro protein–ligand complexes, namely, 6LU7, 5R7Y, and 5R80 (Figure A–C). The number of pharmacophore features varied from one for 5R7Y to two for 5R80 to a 12-feature model for 6LU7. We have then used the “alignment” protocol in LigandScout to align three SB pharmacophores by the reference point (6LU7 as a reference), generate a shared feature pharmacophore, and finally merge pharmacophores and interpolate overlapping features. The final pharmacophore model included 14 features (14FP, Figure D), including seven hydrogen bond donors, two hydrogen bond acceptors, and a positive ionizable feature, representing nonbonding interactions between the Mpro catalytic site and ligands. A group of BioVia scientists recently performed alanine scanning mutagenesis on the catalytic pocket residues of all known Mpro protein–ligand complexes, finding eight critical contacts (Gln192, Gln189, Pro168, Met165, Asn142, Met49, His163, and His41) in no less than one complex and three critical contact residues in no less than four complexes.[55] Six characteristics of the newly developed 14FP hypothesis correlate to interactions with one of the eight critical contact residues (Figure D). The 14FP hypothesis was applied as a 3D probe to filter 8823 candidates from the Drug Bank database (https://go.drugbank.com/) employing the “screening” module of LigandScout with the requirement that at least four features needed to match to identify a hit. During screening, the LigandScout screening module analyzed the alignment of the drug candidates with query features and ranked them on the basis of the pharmacophore fit score (Table S1). Pharmacophore-based screening resulted in 502 hits (Table S1) that were moved forward for the molecular docking-based screening against the 6LU7 binding site.

Figure 3

SB pharmacophore models created using LigandScout. (A) SB pharmacophore model for Mpro–N3 (PDB: 6LU7). (B) SB pharmacophore model for Mpro–Z45617795 (PDB: 5R7Y). (C) SB pharmacophore model for Mpro–Z18197050 (PDB: 5R80). (D) Total of 14 featured merged pharmacophores, generated by aligning 3 SB pharmacophore models by the reference point (6LU7 SB pharmacophore model), generating shared feature pharmacophores and then merging and interpolating overlapping features using LigandScout’s alignment module. The merged pharmacophore contains four hydrophobic features, eight hydrogen bond donor features, and two hydrogen bond acceptor features. Exclusion volumes are not shown. : hydrophobic interaction, : hydrogen bond donor, : hydrogen bond acceptor, and : positive ionizable area.

Top candidates from SB pharmacophore screening (n = 502 with a fit score of 40% and above) were further filtered through molecular docking-based screening. The docking experiment assesses nonbonding interactions of different poses of ligands with the binding pocket of Mpro to rank them based on the docking score (binding affinity: kcal/mol). The docking study also helps establish atomic-level interactions between the Mpro and candidates. Selected docking parameters were validated by redocking the 6LU7 native ligand N3, 6LZE native ligand FHR,[56]6M0K native ligand FJC,[56]6WTK native ligand UED,[57] and 6XHL native ligand V2M[58] in the respective binding site applying the exact procedure, yielding native ligands with a similar location and alignment-like co-crystal structure (Figures S1–S6), confirming the optimal choice of parameters. Seven hydrogen bonding interactions, including Gly143, His164, Phe140, Glu166, Thr190, and Gln189, stabilized the Mpro–N3 complex (Figure S2A). The docked structure demonstrated six of these hydrogen bonds, including the ones with Cys145 and His41 (Figure S2B). Other co-crystal structures, Mpro–FHR, Mpro–FJC, Mpro–UED, and Mpro–V2M complexes, have exhibited binding stabilization by four–six hydrogen bonding contacts with the binding site amino acids, including Phe140, Cys145, Gly143, His164, Glu166, and Gln189 (Figures S3–S6). The in silico-docked structures were in concert with these hydrogen bonding interactions.

The molecular docking-based virtual screening led to 44 promising potential inhibitors based on the affinity of binding of −8 kcal/mol or above (Table ). The docking score of the native ligand, N3, was found to be −8 kcal/mol (Table S3), providing the rationale behind setting it up as the threshold for elimination. The highest-scoring candidate bemcentinib (DB12411, −10.4 kcal/mol) is an experimental selective inhibitor of AXL kinase and currently being tested in phase II clinical trials for nonsmall-cell lung cancer and triple-negative breast cancer.[59−61] Interestingly, bemcentinib is currently under clinical trials against SARS-CoV-2 as the COVID-19-infected lung cells have been observed to upregulate AXL.[62] The AXL receptor signaling pathway is known to downmodulate the interferon-related host immune response and, thus, promote viral infection.[63] Nilotinib, Tasigna, (DB04868, −9.2 kcal/mol) is an approved antagonist of tyrosine kinase against chronic myelogenous leukemia.[64,65] Temafloxacin is a fluoroquinolone-based antibacterial agent approved in 1992 and marketed by Abbott.[66] However, it was withdrawn due to adverse side effects, including allergic reactions and hemolytic anemia. Enasidenib (DB13874, −8.7 kcal/mol) is an orally bioavailable approved drug against refractory acute myeloid leukemia in patients with specific mutation in isocitrate dehydrogenase 2 (IDH2).[67] Tegobuvir prevented viral RNA replication in vitro and was found to have superior efficacy in individuals with hepatitis C virus (HCV).[68] Last year, Ark Biosciences reported the successful completion of a phase II clinical trial of ziresovir (DB15145, −8.4 kcal/mol) in infants with respiratory syncytial virus (RSV) infection.[69]

Danoprevir is one of the most potent HCV protease inhibitors (DB15401, −8.3 kcal/mol) and is being studied for remedy of COVID-19 patients.[70] Selinexor (DB11942, −8.4 kcal/mol) is an anticancer medication that inhibits nuclear export selectively.[71] It has been approved by the FDA for the therapy of refractory multiple myeloma and diffuse large B-cell lymphoma.[72−74] Currently, selinexor is being investigated for the cure of SARS-CoV-2 for its potential to decrease the ability of the virus to replicate and reduce the release of a biochemical, IL-6, responsible for inflammation and organ damage.[75] Defibrotide sodium (DB04932, −8.2 kcal/mol) has been certified by the U.S. FDA in 2016 for the remedy of hepatic veno-occlusive disease, aka sinusoidal obstructive syndrome, with renal or pulmonary disfunction.[76] It is now being tested for the deterrence and therapy of respiratory distress and COVID-19 cytokine release syndrome.[77] Nafamostat (DB12598, −8.2 kcal/mol), a highly potent synthetic serine protease inhibitor, has been investigated for acute pancreatitis.[78] Nafamostat has been observed to significantly inhibit viral transmembrane serine protease 2 enzyme, leading to the inhibition of SARS-CoV-2 from infecting human lungs.[79] Nafamostat is being tested in clinical studies against COVID-19. Indinavir, a protease inhibitor, is utilized to reduce the amount of HIV in human body, so the body’s immune system can work better. It is known to keep the viral polyprotein precursor from cleaving into individual functional protein subunits found in the infectious HIV-1. There has been a number of recent reports recommending indinavir for COVID-19 treatment.[80−82] Montelukast, a leukotriene inhibitor, is an FDA-certified drug for the cure of asthma and allergies and is prescribed to control wheezing and shortness of breath. There have, recently, been many publications with potential benefits of montelukast in SARS-CoV-2 treatment.[83,84] Recent publications suggest that dipeptidyl peptidase 4 (DPP4) inhibitors, for example, gemigliptin and retagliptin, may help reduce virus entrance and reproduction in the respiratory system.[85,86] DDP4 inhibitors can also disrupt prolongation of cytokine storm and infection in the lungs of COVID-19 sufferers.[85,86] Recent findings suggest that COVID-19 treatment may stimulate cardiovascular impairment in people with no abnormalities.[87] The high morbidity and mortality number of COVID-19 are attributed to its adverse impacts on the lungs and heart.[88] The mitogen-activated kinase (MAPK) pathway is critically involved in the secretion of cytokines and linked with severe lung damage and myocardial dysfunction.[88] Thus, MAPK inhibitors, such as SB220025 (DB04338, −8.1 kcal/mol), are being looked into for their ability to significantly diminish cardiovascular damage and morbidity in COVID-19 patients.[87,88] Lopinavir, known for antagonizing antiretroviral protease, is applied against HIV infection. It is under investigation in conjunction with ritonavir for individuals with serious COVID-19 symptoms.[23] In summary, molecular docking-based virtual screening has led to many potential SARS-CoV-2 therapeutics, and many of these are under clinical trial to further validate our experiments. Figure represents conformations of selected bound drugs, with high binding affinities, in the enzyme binding site.

Figure 4

Drugs with high binding affinities in the catalytic site of Mpro (PDB: 6LU7). (A) Bemcentinib, (B) nilotinib, (C) enasidenib, (D) tegobuvir, (E) ziresovir, (F) selinexor, and (G) nafamostat.

As demonstrated in the earlier section, the docking and virtual screening revealed 44 drugs as good binders for the main protease. Next, we performed MD simulations for these 44 drugs to check their binding stabilities. During the MD simulation of all 44 drugs, we found that some of these drugs left the initial binding cleft and show very poor stabilities (see the rmsd plot in Figure ). The visual inspection of MD trajectories reveals that some of those drugs who maintained the same binding conformation as in molecular docking were also very flexible. Therefore, among the 44 drugs screened by molecular docking, only few drugs were showing a stable binding in MD simulations. To quantify these observations, we further performed thermodynamical calculations using MMGBSA for all 44 drugs. The results of the MMGBSA calculations along with the residue-wise interaction map are shown in Table S4. Here, we see that many drugs which have a good docking score have less binding energy, which is in accordance with the poor rmsd of these drugs. Considering the fact that the change in entropy of an FDA-approved drug is roughly around 25 ± 5 kcal/mol (calculated for the FDA-approved drug nilotinib, Table S5), we considered >33.0 kcal/mol as a screening cutoff for a good binder. Applying this screening cutoff, we found that DB04020, DB07145, DB11865, DB12411, DB11779, and DB00471 can be considered as better drugs for the main protease. Among these six drugs, three drugs, DB04020, DB12411, and DB11779, show excellent binding free energies (Table ). To scrutinize the key interactions of these three drugs, the total binding free energy had been decomposed into key residue interactions, as shown in Tables –5. A comparison of key interactions for all three drugs shows that His41 and Met49 are common residues for all three drugs, and they have significant contribution in the free binding stabilization energies. Moreover, DB04020 which shows the best binding affinity interacts very strongly with the Glu166 residue. Interestingly, the interaction of Glu166 was absent in all other drugs except for DB04020. The key interactions of these drugs are shown in Figure .

Figure 5

RMS deviation for all 44 drugs.

Table 2

Binding Free Energy of Top Three Drugs Selected from MD Simulations

drug ID	total ΔH (kcal/mol)
DB04020	–67.53 ± 5.6
DB12411	–46.51 ± 4.0
DB11779	–49.82 ± 4.1

Table 3

Residue-Wise Energy Decomposition for DB04020

drug (DB04020)	residue ID	total interaction energy (kcal/mol)
DB04020	HID 41	–4.43 ± 0.59
DB04020	MET 49	–4.20 ± 0.59
DB04020	PHE 140	–4.39 ± 0.74
DB04020	ASN 142	–3.47 ± 0.99
DB04020	HID 163	–3.68 ± 0.42
DB04020	GLU 166	–12.25 ± 1.38
DB04020	HID 172	–3.63 ± 1.39

Table 5

Residue-Wise Energy Decomposition for DB11779

drug (DB11779)	residue ID	total interaction energy (kcal/mol)
DB11779	THR 25	–8.82 ± 1.12
DB11779	THR 26	–3.13 ± 0.61
DB11779	HID 41	–4.428 ± 0.57
DB11779	SER 46	–6.645 ± 0.83
DB11779	MET 49	–5.483 ± 0.68
DB11779	ASN 142	–4.606 ± 1.91

Figure 6

Key interactions for the top three drugs during the MD simulations.

In Vitro Cytopathic Effect Experiment

We have selected top nine commercially available drug candidates (Table ) from the MD simulation and nilotinib for performing the in vitro SARS-CoV-2 cytopathic effect experiment. We have used remdesivir as a positive control and tested three native ligands, 6LU7, 5R80, and 5R7Y (Table S6) for direct comparison against the drug candidates. We have picked nilotinib because it had a binding free energy of above −25 kcal/mol with the Mpro protein, and it has been reported to have antiviral activity.[89] Moreover, nilotinib has shown strong predicted binding interactions with the SARS-CoV-2 papain-like protease (PLpro) active site in our preliminary study (Figure S7). Recently, nilotinib was demonstrated to be efficacious against SARS-CoV-2 in vitro, strengthening our computational calculation.[90] Thus, we wanted to make sure that it worked in our hands. Out of 10 compounds tested, nilotinib and bemcentinib were able to interfere with the repetition of SARS-CoV-2 in Vero-E6 cells in vitro with the respective EC50 values of 2.6 and 1.2 μM (Figure ), while the control remdesivir demonstrated an EC50 value of 2.8 μM. BMS-690514 (EC50 = 32.78 μM) and lopinavir (EC50 = 30.10 μM) had shown moderate activity, while montelukast was weakly active (Figure S8). Rest of the drug candidates were found to be inactive. Nilotinib and BMS-690514 demonstrated potent to moderate antiviral behavior with one of the highest SI values among the tested compounds, provided their low cytotoxicity (Table ). Thus, a higher concentration of nilotinib and BMS-690514 may be used in in vivo experiments without incurring undesired toxicity. These findings strengthen our computational screening campaign as both bemcentinib and nilotinib have come out to be top hits. Bemcentinib is being tested in clinical trials for SARS-CoV-2, and we are the second to experimentally evaluate the therapeutic potential of nilotinib against SARS-CoV-2. Although lopinavir has been in the clinical trial against severe hospitalized SARS-CoV-2 patients, we, for the first time, report BMS-690514[91] (an orally active EGFR and VEGFR inhibitor) to be a potential treatment for SARS-CoV-2. These findings suggest that nilotinib and BMS-690514 should be assessed in clinical trials for further experimental evaluations.

Table 6

Fifty-Percent Efficacy (EC50) against SARS-CoV-2 WA1 Strain and Vero E6 Cytotoxicity (CC50) Measured for the Drugs Identified. Remdesivir as a Positive Control Candidatea

95% confidence interval (CI); EC50: half maximum effective concentration; CC50: 50% cytotoxic concentration; SI: selectivity index; and IN: inactive. aIndeterminate.

Figure 7

Inhibitory activity of nilotinib, bemcentinib, and remdesivir (positive control) in Vero-E6 cells.

PCA and PLS Regression

The initial PCA and PLS regression modeling involved model optimization and model refinement.[51,52] Five (5) PCs were found to be sufficient for PCA and PLS regression analysis. Figure presents the score plots (PC1 vs PC2) of the PCA study. A total of 98% of the variability in the drug candidate QSAR and 6% of the variability in drug candidate binding energy can be explained by the first PC. PC2 accounted for 1% of the variability in the drug candidate’s QSAR and 26% of the variability in drug candidate binding energy. Groupings of drug candidates into four distinct and notable groups based on the similarity and differences in the drug candidate QSAR were observed in PCA analysis. For instance, S7, S9, S21, and S33 drug candidates are grouped on the far-right corner of the score plot. Interestingly, S7, S9, S21, and S33 have a strong relationship, as they have either imidazole or similar looking guanidine moieties. S7 (DB08553) is under investigation for the treatment of cancer and targeting serine/threonine protein kinase B-raf. S9 (DB08141) is under investigation for the treatment of cancer and targeting cyclin-dependent kinase 2. S21 (DB12598) is under investigation as a serine protease inhibitor for liver transplantation and post-reperfusion syndrome therapy. This is also approved for use in individuals who are on continuous renal replacement therapy owing to acute kidney damage.

Figure 8

Score plots of PCA showing groupings of drug candidates based on similarity and differences in the drug QSAR.

Drug candidates (S23, S29, S31, S35, S39, and S41) containing at least one crucial 5-membered heterocyclic or fused heterocyclic moiety in the molecule are notably identified in the center of the score plots. S23 (DB15674) is an under-investigation antiviral fusion inhibitor, with a benzimidazole moiety, for the treatment of RSV. S29 (DB13506) is a β-lactam antibacterial agent with a fused 5-membered thiazole moiety. S31 (DB11865) is under investigation for colorectal cancer. S35 (DB06994), five-membered oxadiazole moiety, is under investigation targeting DPP4 for T-cell activation and enhanced immune response. S39 (DB08590) is under investigation, with a bezimidazole moiety, targeting interleukin-1 receptor-associated kinase 4 (IRAK-4) known to be involved in signaling innate immune responses from toll-like receptors. Animals with knockout IRAK-4 are known to be more susceptible to the viruses and bacteria.

Drug candidates (S1, S3, S15, S17, and S25) with fused heterocyclic compounds with a single fluorine or chlorine atom are grouped at the left corner of the third quadrant of the score plots. S3 (DB15382) is under investigation for advanced malignant solid tumors. S15 (DB15401) is under investigation for patients with kidney malfunction, C3 glomerulopathy. S17 (DB14822) is under investigation as a dual orexin receptor antagonist for the treatment of insomnia. S25 (DB04020) is a heterocyclic compound, with a single chlorine atom, under investigation targeting the ER-β receptor.

Drug candidates (S5, S11, S19, S27, S37, and S43) with at least one trifluoromethyl group moiety and one aromatic six-membered ring (benzene, pyridine, pyrimidine, or triazine) are grouped in the fourth quadrant of the score plots. S5 (DB13874) is a FDA-certified candidate against AML. S11 (DB11942) is a recently approved drug for multiple myeloma, a type of cancer that originates from antibody-producing plasma cells. S27 (DB12412) is under investigation against type-2 diabetes. It has also been tested for cancer. S37 (DB14898) is under investigation for type-2 diabetes as well. S43 (DB15006) is under clinical trials for septic shock.

The obtained pattern recognition of drug candidates on the score plots is remarkable, which could enhance the ability to rapidly filter and estimate the efficacy and pharmacophore properties of future drug candidates with a similar QSAR in rational drug design.

Validation Study

The summary of the validation study conducted using 19 independent leads showing the calculated energy of binding stabilization by the PLS regression and the actual binding energy calculated from the pharmacophore and molecular docking is shown in Table S2. The obtained low % relative errors (% RE) of predictions revealed the potential of the PLS regression in precisely determining the energy of binding for individual candidates. Root-mean-square-percent-relative error (rms % RE) was applied to further evaluate the overall performance of the PLS regression for the determination of drug candidates’ binding energy. Overall, PLS regression determined the binding energy of individual drug candidates with a low rms % RE of 4.4% and impressively high accuracy of 95.6%.

Conclusions

Medication repurposing is one of the most effective techniques for combating the horrible medical emergency initiated by SARS-CoV-2 infections which are spreading at an exponential rate. Repurposing employs the techniques of screening and identifying known therapeutic agents aiming certain biological targets in SARS-CoV-2. The biggest advantage of drug reprofiling is the potential to bring them to the market in a shorter timeframe and at lower costs as opposed to the orthodox drug discovery. The Mpro of SARS-CoV-2 has earned considerable attention as a drug target because of its important significance in proteolytic processing of polyproteins, translated from the virus RNA, facilitating viral assembly and disease progression. In this study, a sequence of SB pharmacophore and molecular modeling-centered virtual screening allowed us to repurpose a subset of drugs, from the Drug Bank database, against COVID-19 Mpro. After multiple screening, 44 promising hits, out of 8823 agents from the Drug Bank, were recruited for MD simulation experiments to explore their strength of binding in the SARS-CoV-2 Mpro catalytic site. MD results demonstrated that three drugs (DB04020, DB12411, and DB11779) show very strong affinity for binding against the SARS-CoV-2 Mpro catalytic site. His41 and Met49 are common residues, which interact with all drugs. Met165 interacts with most of the candidates, demonstrating interactions with five of the top six hits from MD simulation, Table S4. Gln189 has consistently interacted with most of the candidates. Gly143 and His163 have been reported to be very attractive residues to form hydrogen bonding with potential inhibitors, as we have observed for DB11942 and DB08166.[38,56,92,93] The main chain of Leu141 and side chain of Ser144 have been reported to form hydrogen bonding with a number of ligands.[92] Two of the candidates with a low binding free energy, DB04020 (ΔH = −67.5326 kcal/mol) and DB08166 (ΔH = −28.2756 kcal/mol), have demonstrated interactions with Phe140, revealing its importance as reported previously.[92,94] MD simulations have also shown that an additional interaction with Glu166 is very crucial, and it enhances the binding affinity significantly (Table S4). Therefore, His41, Met49, Phe140, Met165, His163, Gly143, Glu189, and Gln166 could be interesting targets for SB drug design. We have then performed QSAR analysis on the 44 promising hits from molecular docking-centered screening. The use of PCA and QSAR promoted fast drug candidate pattern recognition. PLS regression accurately predicted the values of independent drug candidates’ binding energy with impressively high accuracy. Finally, we have selected 10 drug candidates, as shown in Table , for the in vitro SARS-CoV-2 cytopathic effect experiment. Four out of the 10 drug candidates, nilotinib (DB04868), bemcentinib (DB12411), lopinavir (DB01601), and BMS-690514 (DB11665), have demonstrated capability of interfering in SARS-CoV-2 duplication with EC50 values of 2.6, 1.1, 30.1, and 32.8 μM, respectively, in Vero-E6 cells. It is noteworthy that nilotinib and BMS-690514 have one of the highest SI values among the tested compounds, given their low cytotoxicity values. Currently, bemcentinib and lopinavir are in clinical trials against COVID-19. We are the second to report experimental validation of nilotinib against SARS-CoV-2. Also, we, for the first time, report BMS-690514 to be an experimentally validated potential drug for SARS-CoV-2. These findings demonstrate the strong antiviral potential of nilotinib, bemcentinib, lopinavir, and BMS-690514. While bemcentinib and lopinavir are currently in clinical trials for SARS-CoV-2 treatment, nilotinib and BMS-690514 hold strong potential for future experimental evaluation. Overall, the results of our study are significant and will facilitate rapid screening, assessment, and prediction of pharmacophore properties, behavior, and drug candidate potency. The results of this CADD are envisioned to lead to rational drug design and new discoveries against the COVID-19 pandemic.

Data and Software Availability

LigandScout software (4.4.5 version) used for the pharmacophore development and virtual screening can be downloaded for a one-month free trial at http://www.inteligand.com/. PyRx/AutoDock Vina software used for molecular docking can be downloaded free of charge at https://pyrx.sourceforge.io/. The Amber 18 MD package used for MD and simulation is available at https://ambermd.org/. BIOVIA Discovery Studio Visualizer used for 3D visualization is available free of charge at https://discover.3ds.com/discovery-studio-visualizer-download. VMD used for visualization is available free of charge at https://www.ks.uiuc.edu/Research/vmd/. The QSAR was investigated using Unscrambler Chemometrics software available at https://www.camo.com/unscrambler/.

Table 4

Residue-Wise Energy Decomposition for DB12411

drug (DB12411)	residue ID	total interaction energy (kcal/mol)
DB12411	HID41	–5.27 ± 0.89
DB12411	MET 49	–3.45 ± 0.52
DB12411	MET 165	–4.20 ± 0.59
DB12411	ASP 187	–4.35 ± 1.12
DB12411	GLN 189	–6.49 ± 1.37