In silico and in vitro evaluation of imatinib as an inhibitor for SARS-CoV-2

Abstract The rapid geographic expansion of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the infectious agent of Coronavirus Disease 2019 (COVID-19) pandemic, poses an immediate need for potent drugs. Enveloped viruses infect the host cell by cellular membrane fusion, a crucial mechanism required for virus replication. The SARS-CoV-2 spike glycoprotein, due to its primary interaction with the human angiotensin-converting enzyme 2 (ACE2) cell-surface receptor, is considered a potential target for drug development. In this study, around 5,800 molecules were virtually screened using molecular docking. Five molecules were selected for in vitro experiments from those that reported docking scores lower than −6 kcal/mol. Imatinib, a Bcr-Abl tyrosine kinase inhibitor, showed maximum antiviral activity in Vero cells. We further investigated the interaction of imatinib, a compound under clinical trials for the treatment of COVID-19, with SARS-CoV-2 RBD, using in silico methods. Molecular dynamics simulations verified that imatinib interacts with RBD residues that are critical for ACE2 binding. This study also provides significant molecular insights on potential repurposable small-molecule drugs and chemical scaffolds for the development of novel drugs targeting the SARS-CoV-2 spike RBD. Communicated by Ramaswamy H. Sarma


Introduction
In early December 2019, the Chinese health authorities reported several cases of pneumonia of unknown cause that had originated in Wuhan, a city in the Hubei province of China. The causative agent of this outbreak was identified to be a virus that belonged to the Sarbecovirus subgenus, Orthocoronavirinae subfamily which was previously referred to by its interim name 2019 novel coronavirus (2019-nCoV) (Organization, 2020;Zhu et al., 2020) and was later named as SARS-CoV-2 (of the International, 2020). Due to the rapid spread of COVID-19, the World Health Organization (WHO) declared it a global pandemic in March 2020 (Cucinotta & Vanelli, 2020). As of February 2022, over 425 million cases have been confirmed around the world, resulting in more than 5.8 million deaths (Dong et al., 2020). Although preventive vaccines are now available and several others are in the pipeline, only a few antiviral treatments are available ((CDC), 2021; (NIH), 2021; (WHO), 2020). Since supportive care is the primary recommended treatment, it is imperative to identify repurposable lead compounds to treat COVID-19 patients until a SARS-CoV-2-specific drug is developed.
Although the coronavirus genome consists of numerous conserved druggable enzymes, including papain-like protease (PLpro), 3 C-like protease (3CLpro), non-structural proteins RNA-dependent RNA polymerase (RdRp), and helicase, the development of clinically approved antiviral therapies has proven to be a difficult task (Zumla et al., 2016). The surface structural spike glycoprotein (S), a key immunogenic CoV antigen essential for virus and host cell-receptor interactions, is an important target for therapeutic development. The spike protein consists of an N-terminal S1 subunit (receptor binding) and a C-terminal S2 subunit (membrane fusion). The S1 subunit contains the receptor-binding domain (RBD) which attaches to the host membrane, thus playing an important role in viral entry. SARS-CoV-2 utilizes the ACE2 receptor for entry and the transmembrane protease, serine 2 (TMPRSS2) for spike protein priming (Hoffmann et al., 2020). Crystallographic studies have shown that SARS-CoV-2 binds to the ACE2 receptor, with a binding mode nearly identical to that of SARS-CoV (Lan et al., 2020;Letko et al., 2020;Wrapp et al., 2020;Zhou et al., 2020a). The binding affinity of the ACE2 receptor to the RBD of the SARS-CoV-2 spike protein is reported to be significantly higher as compared to SARS-CoV (Lan et al., 2020;Wrapp et al., 2020).
Based on the importance of virus membrane fusion events in the viral life cycle and its infectivity, the spike protein of SARS-CoV-2 was targeted for drug screening. This study utilizes in silico methodology followed by in vitro experimental validation to screen existing FDA-approved small molecule drugs specific to the RBD of the spike protein of SARS-CoV-2 to identify repurposable drugs targeting further clinical validation.

Protein structure prediction and validation
A SWISS-MODEL server (Waterhouse et al., 2018) was used to construct a homology model of the SARS-CoV-2 spike protein using the crystal structure of the SARS-CoV-2 spike protein (PDB:6VSB_chain A) as the template (Wrapp et al., 2020). The genome sequence Wuhan-Hu-1 (GenBank: MN908947.3) was used as a representative of the SARS-CoV-2. Spike protein sequence (GenBank: QHD43416.1) was used as the target sequence . The SWISS-MODEL Structure Assessment Tool was used to validate the quality of the predicted model.

Molecular docking
Around 5,800 compounds, including 3,700 nucleoside-like compounds from the Enamine Targeted Antiviral Library (enamine.net) and 2,100 Food and Drug Administration (FDA)-approved drugs from the ZINC15 database (Sterling & Irwin, 2015) were used for molecular docking. All molecules were prepared with obabel (O'Boyle et al., 2011) from .sdf or .mol2 format to .pdbqt format. The 3 D compound structures from the Enamine library were resolved by obabel -gen3d command. The docking file of the protein model was prepared with MGLTools v1.5.4 (Morris et al., 2009) and the molecules were docked at the RBD of the spike protein via Autodock Vina 1.1.2 (Trott & Olson, 2010). The grid box of 40 Â 60 Â 30 size with 1.0 Å spacing was fixed around the RBD (Thr323-Val511) of the spike protein. Each docking was done in three replicates, and the conformation with the highest binding score was recorded. The batch processing of docking and data collection was performed using an inhouse python script which is deposited in GitHub (https:// github.com/sfernando-BAEN/ELIXIR-A-Vina-Batch-Screening-Module). Data were analyzed statistically using R studio (Allaire, 2012) and graphs were constructed with ggplot in R (Warnes et al., 2016). The ligand-receptor interactions were visualized using Schr€ odinger Maestro (Release, 2017a), and molecules with high docking scores were selected from each screening library for further studies. The docking scores for the selected molecules were also validated using Schr€ odinger Glide SP mode with default protocol.

Molecular dynamics (MD) simulations
The recorded stable conformations of the selected molecules from docking studies were further analyzed through MD simulations which were performed using Schr€ odinger-Desmond (Release, 2017b). The protein-ligand complexes were prepared by Schr€ odinger's Protein Preparation Wizard (Wizard, 2014). Structures were refined by capping the C and N termini of the protein, adding missing hydrogens, optimizing hydrogen bonds, and were finally minimized using the OPLS3e force field (Roos et al., 2019). A molecular dynamics system was built for each of the processed protein-ligand complexes. The system was solvated in an orthorhombic box using the SPC solvent model with a buffer distance of 10 Å. The system charge was neutralized with Na þ or Clions using a 0.1 M concentration of NaCl. All simulations were performed using Desmond's default relaxation protocol under the OPLS3e force field. Heavy atoms in the system were initially minimized with restraints under 10 K, then further restrained by increasing the temperature to 300 K, and lastly relaxed under 300 K NPT ensemble to get the equilibrium status for the system. After relaxation, the simulations were performed under a 300 K NPT ensemble at 1.01325 bar pressure for 100 ns. Trajectories and related energy terms were recorded at an interval of 200 ps. Post-simulation analyses including complex root mean square deviation (RMSD), ligand/protein root mean square fluctuation (RMSF), and interactions analysis were performed via Schr€ odinger Simulation Interaction Diagrams. The python script (binding_sasa.py) from Schr€ odinger was used to calculate solvent accessible surface area (SASA) of ligand before and after receptor-ligand binding.

Binding free energy calculations
The Molecular Mechanism-Generalized Born Surface Area (MM-GBSA) binding free energies (Srinivasan et al., 1998) were calculated for the simulation trajectories. Five hundred frames (one every 200ps) were extracted over a duration of 100 ns. Overall, 250 frames from the last 50 ns (stabilized simulation trajectory) were used for MM-GBSA calculation for each complex. Prime MM-GBSA used the VSGB 2.1 solvation model  with the OPLS3e force field.
The binding free energy (DG) was calculated using the equation below: Further, the total DG ðbindÞ can be decomposed into the following contributing components: . Lower DG bind (more negative values) indicate a higher affinity of the ligand to the protein. However, the MM-GBSA DG bind estimates reported in kcal/mol are different from the true (experimental) binding free energy. The MM-GBSA method is based on end-point free energy calculations. It accounts for the molecular mechanics and solvation terms while omitting the conformational entropy terms. This approximation causes the calculated DG bind to be more negative than the actual affinity (Hou et al., 2011). For drug screening purposes, the MM-GBSA method is widely used since only relative binding free energies are needed .

Vero-TMPRSS2 cell line production
Vero-TMPRSS2 cells were produced by retroviral transduction.
To produce the retrovirus, 10 lg pQXCIH-TMPRRS2-HA was co-transfected with polyethylenimine (PEI) with 6.5 lg pBSgag-pol (Addgene #35614) and 5 lg pMD2.G in a 10 cm dish of 70% confluent HEK-293T cells in Opti-MEM I (1X) þ GlutaMAX. Retroviral particles were harvested at 72 hours post-transfection, cleared by centrifugation at 2000 x g, filtered through a 0.45 lm low protein-binding filter (Millipore), and used to transduce Vero cells. Polybrene (Sigma) was added at a concentration of 4 lg/ml to enhance transduction efficiency. Transduced cells were selected with hygromycin B (Invitrogen).

VSV Delta G rescue
The protocol for VSV-G pseudovirus rescue was adapted from Whelan and colleagues (1995). Briefly, a 70% confluent 10 cm dish of HEK-293T cells was transfected with 10 mg pVSV-eGFP-dG, 2 mg pCAG-VSV-N (nucleocapsid), 2 mg pCAG-VSV-L (polymerase), 2 mg pMD2.G (glycoprotein, VSV-G), 2 mg pCAG-VSV-P (phosphoprotein) and 2 mg pCAGGS-T7Opt (T7 RNA polymerase) using PEI at a ratio of 1:3 (DNA: PEI) in Opti-MEM I (1X) þ GlutaMAX. Forty-eight hours post-transfection the supernatant was transferred onto new plates transfected 24 hours prior with VSV-G. After a further 48 hours, these plates were retransfected with VSV-G. After 24 hours the resulting pseudoviruses were collected, cleared by centrifugation at 2000 x g for 5 minutes, and stored at À80 C. Subsequent VSV-G pseudovirus batches were produced by infecting VSV-G transfected HEK-293T cells with VSV-G pseudovirus at an MOI of 0.1. Titres were determined by preparing 10-fold serial dilutions in Opti-MEM I (1X) þ GlutaMAX. Aliquots of each dilution were added to monolayers of 2 Â 10 4 Vero cells in the same medium in a 96-well plate. Three replicates were performed per pseudovirus stock. Plates were incubated at 37 C overnight and then scanned using an Amersham Typhoon scanner (GE Healthcare). Individual infected cells were quantified using ImageQuant TL software (GE Healthcare). All pseudovirus work was performed in a Class II Biosafety Cabinet under BSL-2 conditions at Erasmus Medical Center.

Coronavirus S pseudovirus production
For the production of MERS-CoV, SARS-CoV, and SARS-CoV-2 S pseudovirus, HEK-293T cells were transfected with 10 mg S expression plasmids. Twenty-four hours post-transfection, the medium was replaced for in Opti-MEM I (1X) þ GlutaMAX, and cells were infected at an MOI of 1 with VSV-G pseudotyped virus. Two hours post-infection, cells were washed three times with OptiMEM and replaced with a medium containing anti-VSV-G neutralizing antibody (clone 8G5F11; Absolute Antibody) at a dilution of 1:50,000 to block remaining VSV-G pseudovirus. The supernatant was collected after 24 hours, cleared by centrifugation at 2000 x g for 5 minutes, and stored at 4 C until use within 7 days. Coronavirus S pseudovirus was titrated on VeroE6 cells as described above.

Pseudovirus assay
Transduction experiments were carried out by incubating pseudovirus with imatinib at concentrations ranging from 0-125nM in Opti-MEM I (1X) þ GlutaMAX for 1 hour at 37 C. Pseudovirus-imatinib mixes were added to monolayers of 2 Â 10 4 Vero or Vero-TMPRSS2 cells in a 96-well plate. Plates were incubated for 16 hours before quantifying GFP-positive cells using an Amersham Typhoon scanner and ImageQuant TL software.

In vitro toxicity of imatinib
To determine the toxicity profile of imatinib, we performed the MTT assay using a 1-hr and an 8-hr design. Briefly, a serial dilution of imatinib was prepared and incubated on Vero cells for 1 hr at 37 C. Subsequently, cells were washed, further cultured for eight hrs. In the 8-hr design, cells were incubated with a serial dilution of imatinib for eight hours without a washing step.

In vitro efficacy of imatinib
We tested serial dilutions of imatinib for its ability to neutralize SARS-CoV-2 (German isolate; GISAID ID EPI_ISL 406862; European Virus Archive Global #026 V-03883) using a plaque reduction neutralization test (PRNT) as previously described (Okba et al., 2020). Fifty lL of the virus suspension (200 spot forming units) was added to each well and incubated at 37 C for either 1 hr. Following incubation, the mixtures were added to Vero cells and incubated at 37 C for either 1 hr or 8 hrs. The cells incubated for 1 hr were then washed and further incubated in medium for 8 hrs. After the incubation, the cells were fixed and stained with a polyclonal rabbit anti-SARS-CoV antibody (Sino Biological; 1:500). Staining was developed using a rabbit anti-SARS-CoV serum and a secondary alexa-fluor-labeled conjugate (Dako). The number of infected cells per well was counted using the ImageQuant TL software.

Protein structure prediction and validation
A model for SARS-CoV-2 spike protein was constructed using the crystal structure (6VSB_chain A) to correct missing residues. The amino acid sequence identity between the target sequence (GenBank: QHD43416.1) and template (6VSB_chain A) was 99.58%. The SARS-CoV-2 model showed an RMSD of 0.4683 Å relative to the crystal structure (6VSB_chain A), superimposed structures shown in Figure 1A.
Structure assessment of the predicted model using the Ramachandran plot showed 90.04% residues in the most favored regions with 2.03% outliers ( Figure 1B). None of the outliers contained the residues present at the active site of the protein. This model was further used for in silico studies.

Virtual screening
In virtual screening, approximately 5,800 compounds from the Enamine Antiviral and ZINC15 FDA libraries were docked against the SARS-CoV-2 spike RBD protein. The distribution of the docking scores for the screened compounds amongst the respective sets is shown in Figure 2. The output was analyzed for common classes of drugs with the highest (most negative) docking scores (threshold: À6 kcal/mol) and manually reviewed to omit dye molecules such as trypan blue. Seven compounds were selected based on docking scores. Three compounds from Enamine Antiviral library, Antiviral825, Antiviral2038, and Antiviral2981 had docking scores À6.30, À6.20, and À6.00 kcal/mol; and the four compounds from ZINC15 FDA library, ponatinib, imatinib, ergotamine, and glecaprevir had docking scores of À7.63, À6.80, À7.70 and À7.20 kcal/mol respectively. The screened compounds had the highest scores within their respective sets and had one or more binding conformations at the RBD of the spike protein. The docking scores for these compounds were also validated using Glide (Supplementary data). Initially, when this study was performed, no glycan molecule was added to the RBD protein model. Later, a glycan molecule at Asn343 was added to the protein model and results were validated via Glide docking (Table S2 under Supplementary data). Abl tyrosine kinase inhibitors (ponatinib and imatinib) were a common class of drugs that received high docking scores. A detailed description of the screened drugs is given in Table S1 under Supplementary Data. The docking scores, type of interactions, and important interacting residues at the RBD are mentioned in Table 1. Individual interaction  diagrams for all the seven compounds with the spike protein are provided in Figure S1 under Supplementary Data.
The binding conformations and docking scores for the seven screened compounds at the RBD are shown in Figure  3A and B, respectively. The free energies of binding (DG) calculated over the last 50 ns of the MD simulations using the MM-GBSA method indicated that the compounds in the FDA-approved drug library had significantly more negative DG (i.e., higher affinities) values than the drugs in the enamine library, which are in agreement with the docking scores. Ponatinib, an Abl tyrosine kinase inhibitor used to treat leukemia (Rossari et al., 2018), displayed the most negative binding free energy À57.82 ± 6.14 kcal/mol, while imatinib, glecaprevir, and ergotamine also show binding free energies within a close range to ponatinib. The high affinity of the screened compounds is visible when compared with the negative control dimethyl sulfoxide (DMSO), which is ineffective against coronaviruses (Sisk et al., 2018).
Among the screened drugs, Abl tyrosine kinase inhibitor imatinib has been previously shown to inhibit SARS-CoV and MERS-CoV by blocking endosomal fusion at the cell-culture level (Basha, 2020;Coleman et al., 2016;Sisk et al., 2018;Zumla et al., 2016). As shown in Figure 3A (inset), extracted from the last frame of the MD simulation shows imatinib forming hydrogen bonds and pi-cation interactions with Tyr489 and pi-pi stacking interactions with Phe456, both of which are critical for interacting with ACE2 protein (Lan et al., 2020;Ortega et al., 2020).
The flexibility of protein-ligand complexes was evaluated by calculating the fluctuations of each residue. The RMSF diagram ( Figure 5)  Antiviral2038 had a higher SASA value initially but stabilized around the average SASA value similar to other compounds after 50 ns, as shown in Figure 6. All ligands in their receptorbound form reported a lower SASA value compared to their corresponding SASA values in unbound form. These results indicate that ligands remain bound to the protein throughout the simulation and the protein-ligand complexes are stable in solution.

In vitro efficacy of screened compounds
The viral plaque assay was performed to validate the in silico results and evaluate the in vitro activity of the screened compounds. (Due to a supply shortage, ponatinib and ergotamine were unavailable for purchase and hence could not be included in the initial viral plaque assays). Here, imatinib showed significant percent inhibition ($100%) as compared to the other compounds tested ($10-30%), as shown in Figure 7. Previous studies have shown imatinib to inhibit SARS-CoV and MERS-CoV by blocking endosomal fusion in cell-culture (Basha, 2020;Coleman et al., 2016;Dyall et al., 2014;Sisk et al., 2018;Zumla et al., 2016). It has been suggested that tyrosine-kinase inhibitors do not affect the cleavage of the spike protein but inhibit spike-mediated endosomal fusion (Coleman et al., 2016;Sisk et al., 2018;Zumla et al., 2016). The high affinity of tyrosine-kinase inhibitors towards the spike protein is deduced from the in silico results, where both imatinib and ponatinib have shown highly negative binding free energies. Based on promising in silico data, and initial viral plaque assay results, imatinib was chosen to be advanced for further experimental validation.

In vitro efficacy of imatinib
First, we evaluated the toxicity of imatinib when incubating the compound on Vero cells for one hour or eight hours. In the experiments where the compound remained on the cells for eight hours toxicity was measured at concentrations of 25 mM, 12.5 mM, 6.3 mM, and 3.2 mM. However, in the 1-hour design, no toxicity was observed. Next, we evaluated the ability of imatinib to inhibit replication and entry. At concentrations as low as 0.2 mM the compound was effective in suppressing 50% of plaque formation in the 8-hr design, and the  IC50 value determined using linear regression was 130 nM. Consistent with the toxicity data, toxicity was observed between 25 and 3.2 mM. The compound also showed efficacy in the 1-hour design, with higher IC50 values. These data indicate that imatinib inhibits virus replication in vitro as shown in Figure 8A and B.

Imatinib inhibits fusion
To evaluate if imatinib inhibits viral entry, we performed two fusion assays: endosomal (Vero) and plasma membrane (Vero-TMPRSS2) as shown in Figure 9A and B, respectively. Based on cytotoxicity, at concentrations below 15 nM, no toxicity was observed microscopically (red arrow in the graph). VSV-G control revealed 100% infectivity (cytopathic effect at every concentration below this, suggesting the inhibitor did not affect VSV-G entry. VSV-G particles cells do not carry spike proteins and thus, no significant entry inhibition occurred, suggesting that entry inhibition is likely mediated through the spike protein. However, the effect on Vero-TMPRSS2 cells was less clear for any of the coronaviruses used when compared to the VSV-G control. A similar level of toxicity was observed in these cells. It is worth noting that toxicity is probably the result of incubating cells with imatinib for 16 hours in the assay. Taken together, there is evidence that imatinib inhibits spike fusion and prevents viral entry, possibly by preventing endosomal entry.

Discussion
The interaction between human cell receptor ACE2 and SARS-CoV-2 spike RBD is the very first step of viral infection, thus, a critical target for drug development (Hoffmann et al., 2020). Studies have shown that cells expressing ACE2 receptors on their cell surface are susceptible to SARS-CoV-2 infection compared to those without ACE2 (Zhou et al., 2020b). Previously reported outbreak-causing coronavirus, SARS-CoV also uses the same cell receptor for entry (Hoffmann et al., 2020). However, the binding affinity of SARS-CoV-2 RBD to ACE2 is much stronger compared to that of SARS-CoV (Lan et al., 2020). Small molecules targeting the SARS-CoV-2 RBD Figure 7. Percent inhibition compared to the number of plaques on cells. Here, a 2-fold dilution of the compounds was done in duplo. Then 300 TCID50 of SARS-CoV-2 was added to each well, and plates were incubated at 37 C for 1 hour. Then, the mixes were added onto Vero cells and incubated for 8 hours at 37 C. Subsequently, cells were fixed for 15 min with 2% PFA, followed by another 15 min fixation with 70% ethanol. Fixed cells were stained with a monoclonal antibody, followed by AlexaFluor 488. Here imatinib shows significant percent inhibition as compared to other compounds tested.  could protect against the virus by blocking its entry into the host cell. Given the virus's pathogenicity and infectivity, there is an urgent need for anti-SARS-CoV-2 drugs. Health experts across the globe are trying to use existing clinically approved drugs to treat patients until a specific drug is developed. The present study, using a combination of computational techniques followed by in vitro experiments, identified imatinib, an FDA-approved anti-cancer drug as a potential treatment of SARS-CoV-2 infection. Our results indicate that imatinib inhibits SARS-CoV-2 in Vero cells at sub-micromolar levels while previously reported values were in the lower micromolar range (Coleman et al., 2016;Han et al., 2021;Sisk et al., 2018;Touret et al., 2021). This discrepancy may be due to the differences in experimental methods and assay conditions. Also, computational analyses suggest that imatinib forms strong interactions with the RBD residues critical for interaction with ACE2. Our results suggest that imatinib prevents viral replication and that, in analogy with other studies using coronaviruses, the most likely mode of action of imatinib is via inhibition of the spike protein. Nevertheless, we cannot completely exclude that inhibition of tyrosine kinase also resulted in significant inhibition of virus replication, as previously described (Coleman et al., 2016;Sisk et al., 2018). Therefore, it is likely that imatinib causes both inhibition of cellular kinase and virus fusion resulting in inhibition of virus replication. However, the potency of imatinib and its use as a promising repurposable candidate against SARS-CoV-2 remains a controversial topic due to inconsistent outcomes in preclinical evaluations (Bernal-Bello et al., 2021;Han et al., 2021;Lin et al., 2021;Morales-Ortega et al., 2021;Touret et al., 2021;Zhao et al., 2020).

Conclusion
In the present study, around 5,800 molecules were virtually screened using molecular modeling tools. Further, five molecules from those that reported docking scores lower than À6 kcal/mol were investigated for in vitro activity against SARS-CoV-2. Amongst these, imatinib showed the highest antiviral activity in Vero cells. MD simulations confirmed that imatinib interacts with RBD residues that are critical for ACE2 binding. Our results provide further evidence supporting the ongoing clinical trials (ClinicalTrials.gov Identifier: NCT04346147, NCT04357613, NCT04953052, NCT04794088, NCT05220280, NCT04394416, and NCT04422678) for COVID-19 patients with imatinib. This study also provides valuable insights on several other chemical scaffolds that could be optimized through medicinal chemistry to increase target specificity and achieve optimum potency against SARS-CoV-2.