Repurposing of Drugs for SARS-CoV-2 Using Inverse Docking Fingerprints.

Severe acute respiratory syndrome coronavirus 2 or SARS-CoV-2 is a virus that belongs to the Coronaviridae family. This group of viruses commonly causes colds but possesses a tremendous pathogenic potential. In humans, an outbreak of SARS caused by the SARS-CoV virus was first reported in 2003, followed by 2012 when the Middle East respiratory syndrome coronavirus (MERS-CoV) led to an outbreak of Middle East respiratory syndrome (MERS). Moreover, COVID-19 represents a serious socioeconomic and global health problem that has already claimed more than four million lives. To date, there are only a handful of therapeutic options to combat this disease, and only a single direct-acting antiviral, the conditionally approved remdesivir. Since there is an urgent need for active drugs against SARS-CoV-2, the strategy of drug repurposing represents one of the fastest ways to achieve this goal. An in silico drug repurposing study using two methods was conducted. A structure-based virtual screening of the FDA-approved drug database on SARS-CoV-2 main protease was performed, and the 11 highest-scoring compounds with known 3CLpro activity were identified while the methodology was used to report further 11 potential and completely novel 3CLpro inhibitors. Then, inverse molecular docking was performed on the entire viral protein database as well as on the Coronaviridae family protein subset to examine the hit compounds in detail. Instead of target fishing, inverse docking fingerprints were generated for each hit compound as well as for the five most frequently reported and direct-acting repurposed drugs that served as controls. In this way, the target-hitting space was examined and compared and we can support the further biological evaluation of all 11 newly reported hits on SARS-CoV-2 3CLpro as well as recommend further in-depth studies on antihelminthic class member compounds. The authors acknowledge the general usefulness of this approach for a full-fledged inverse docking fingerprint screening in the future.

Introduction

Severe acute respiratory syndrome coronavirus 2 or SARS-CoV-2 is a virus that belongs to the Coronaviridae family and is named after the crown serrations on its surface (Kahn, and McIntosh, 2005; Cui et al., 2019). It is a single-stranded positive sense RNA (+ssRNA) virus (Gorbalenya et al., 2020; Zhu N. et al., 2020). This group of viruses commonly causes colds but has tremendous pathogenic potential. In humans, an outbreak of SARS (severe acute respiratory syndrome) caused by the SARS-CoV virus was first reported in mainland China and Hong Kong in 2003, followed by 2012, when the Middle East respiratory syndrome coronavirus (MERS-CoV) led to an outbreak of Middle East respiratory syndrome (MERS) in Saudi Arabia, mainland China, United Arab Emirates, and the Republic of Korea (Hilgenfeld and Peiris, 2013; DeWit et al., 2016). More recently, swine acute diarrhea syndrome coronavirus (SADS-CoV) causing severe acute porcine diarrhea syndrome has also been described with high porcine pathogenicity on top of a variety of documented coronaviruses in other animals (Yang et al., 2020; Lin et al., 2021).

The emergence of the COVID-19 disease caused by the SARS-CoV-2 pathogen was reported in major media in December 2019 to have originated in Wuhan, Hubei, China (Wu et al., 2020), and spread worldwide in the first months of 2020, causing a pandemic of the COVID-19 disease (Li et al., 2020; Wang, et al., 2020). COVID-19 is a serious socioeconomic and global health problem that has claimed more than 4,294,225 lives at the time of writing this article (Nicola et al., 2020). Indeed, the majority of cases present with only mild symptoms, while a variable percentage (0.2%–5%) of patients progress to pneumonia and multiorgan failure, which can lead to death, especially without medical assistance at the secondary healthcare level (O’Driscoll et al., 2021; Malik et al., 2021). The medical and academic communities, as well as the pharmaceutical industry, have responded immediately with intensive research campaigns aimed primarily at uncovering pathogenicity mechanisms, researching new drugs and developing vaccines, accompanied by new social guidelines and the dissemination of information and good hygiene practices by the relevant authorities (Fry et al., 2020; Meier et al., 2020). Registered SARS-CoV-2 vaccines are available (Chen C. Z. et al., 2020; Amanat and Krammer, 2020), and they represent the forefront in battle against COVID-19, but the high viral mutation rate, which can lead to structural changes in key viral proteins, may render available vaccines ineffective (Naqvi et al., 2020). In late 2020, a novel SARS-CoV-2 alpha variant (B.1.1.7; Volz et al., 2021) and a beta variant (B.1.351; Tegally et al., 2020) were reported, followed by a gamma variant (P.1; Faria et al., 2021) and a new SARS-CoV-2 variant delta/delta+ in 2021 (B.1.617/AY.1), causing new infections and reinfections that are slowly spreading throughout the world (Moelling, 2021; Roy, et al., 2021). We can respond by developing novel vaccines, but even with novel technologies such as mRNA, the response time is substantial (Badgujar et al., 2020; Verbeke et al., 2021). Therefore, the development of other therapeutic options and novel drug approaches are essential for the future control of coronavirus infections (Kaddoura et al., 2020; Pooladanda et al., 2020; Sarkar et al., 2020).

To date, there are only a handful of therapeutic options to combat this disease, with only one direct-acting antiviral, remdesivir, conditionally approved in Taiwan, followed by a rapid succession of conditional approvals in the EU and Canada. Following these conditional approvals, an emergency approval for remdesivir (a prodrug of GS-441524) was granted in the US and Japan in May 2020 (Lamb, 2020). There is a tremendous research effort underway to develop novel drugs (Jin et al., 2020; Günther et al., 2021; Jukič et al., 2020), but given the immediate need for active compounds against SARS-CoV-2, the strategy of drug repurposing represents one of the fastest options toward this goal (Dotolo et al., 2021; Singh et al., 2020; Gatti and De Ponti, 2021). Most notably, the most commonly reported and direct-acting repurposed drugs include the antiviral agents favipiravir (Coomes and Haghbayan, 2020), lopinavir (dynamic)–ritonavir (kinetic) (Ye et al., 2020), ribavirin (Khalili et al., 2020), interferons (Zhou et al., 2020), the anthelmintic ivermectin (Schmith et al., 2020), and the antimalarials chloroquine (Cortegiani et al., 2020) or hydroxychloroquine (Meo et al., 2020), all shown in Figure 1. Multiple reviews on this subject beyond the scope of this article are available to the reader (Saha et al., 2020; Sourimant et al., 2021). Similar to remdesivir, favipiravir is a viral RdRp inhibitor, ribavirin inhibits IMPDH2, and lopinavir together with ivermectin inhibits viral 3CLpro. Chloroquine/hydroxychloroquine is thought to modulate viral endosome maturation and interact with sigma receptors (Abate et al., 2020; De et al., 2021). Other targeting approaches such as ACE2−RBD interaction have also been examined for drug repurposing (Hanson et al., 2020; Wei et al., 2020; de Oliveira et al., 2021) and compound sets evaluated on cell lines in vitro (Chen W. H. et al., 2020; Riva et al., 2020; Bakowski et al., 2021). Especially on examination of sigma receptor ligands, a key observation was made where phospholipidosis was a shared mechanism underlying the antiviral activity of many repurposed drugs (hydroxychloroquine, azithromycin, amiodarone). Mehanistically, this disrupts lysosomal lipid catabolism and trafficking and results in an in vitro correlation between drug-induced phospholipidosis and antiviral activity disrupting the clear mechanism-based design decisions (Tummino et al., 2021). This is especially evident in amphiphilic compounds and depends on the physicochemical properties (cLogP ≥ 3 and pKa ≥ 7.4) of drugs. Therefore, in this work, this was especially considered and compounds flagged in order to focus on molecules with therapeutic potential.

FIGURE 1

Commonly repurposed drugs against SARS-CoV-2 with proposed direct action (ritonavir is mainly a pharmacokinetic modifier for its partner in combination therapy).

In the context of anti-coronavirus therapeutics—all repurposing approaches with a large body of conflicting data—we sought to conduct a transparent in silico drug repurposing study using known FDA-approved drugs on a well-described SARS-CoV-2 target 3CLpro or Mpro (Zhu W. et al., 2020). This target is a first choice for repurposing campaigns due to its extensive experimental support (Gordon et al., 2020), available crystallographic data, and good biological evaluation data (Anand et al., 2003; Chiou et al., 2021; Osipiuk et al., 2021). The protease is an attractive target as it plays a central role in the viral life cycle by processing the viral polyproteins pp1a and pp1ab at multiple distinct cleavage sites and complementary reports on repurposing research are available, further contextualizing the work herein (Kuzikov et al., 2021). We proceeded to compute the inverse docking fingerprints to provide a focused outlook on a typical repurposing scenario and propose a consensus on the identified hit compounds. The work reveals a rarely studied in silico selectivity, and the authors have not overlooked the usefulness of this approach for future inverse docking fingerprint screening experiments to identify compounds that behave similarly across a large number of targets or explore the interactome (Sadegh, et al., 2020; Figure 2).

FIGURE 2

Virtual screening coupled to inverse docking protocol for drug repurposing. The first step is database selection; next is the VS campaign followed by a filtering step (filter off pains) with a final inverse docking fingerprint calculation to examine the target space where the compounds produce favorable binding poses (to this end, the ProBiS database was employed to conduct two experiments in parallel—one with all viral proteins and the other with focused Coronaviridae proteins). Hit selection was performed based on VS scores and prior available experimental data (ChEMBL) along with examination of inverse docking fingerprints.

The inverse docking procedure has already been used as a tool for drug repurposing before with positive results (Kharkar et al., 2014; Wang et al., 2019). As an example, it revealed the potential new targets for tanshinone IIA used in the treatment of acute promyelocytic leukemia (Chen, 2014). More recently, Ribone et al. (2021) applied inverse docking using multiple scoring functions for target proposal on SARS-CoV-2-repurposed drugs and discovered new potential targets for drugs with experimentally determined activity against SARS-CoV-2.

Materials and methods

Database preparation

The compounds were obtained from the Drug Bank Database as a subset of 2,511 drugs in sdf format (Wishart et al., 2018). The database was further expanded to 4,711 entries with calculation of tautomeric structures and ionization at pH 7.4 with further standard structure preparation steps such as enumeration of undefined chiral centers, removal of structural defects, and 3D structure minimization with optimization (using the OPLS3e force field) toward the final 3D conformation. For this work, the LigPrep tool from Schrödinger (Release Schrodinger 2020-4, Schrödinger, LLC, New York, NY, United States, 2020) was used. The exact parameters were as follows: ligprep -bff 16 -i 2 -ph 7.4 -pht 1.0 -s 8 -orig_file -orig_file_index 1 -isd input -osd output (Shelley et al., 2007; Greenwood et al., 2010).

Ivermectin B1a/B1b and selamectin were modeled using the Avogadro chemical editor (Hanwell et al., 2012) and then optimized using Gaussian 16 (Frisch et al., 2016) in conjunction with the B3LYP method and the 6-31G(d) basis set, in the absence of a suitable 3D structure. The 3D structures of all other molecules used were obtained from the Drug Bank Database, ionized using LigPrep Ioniser software (Release Schrodinger 2020-4, Schrödinger, LLC, New York, NY, United States, 2020) at pH 7.4 ± 1 and their geometry optimized using the RDKit geometry optimization node in KNIME software.

Target preparation

We chose a well-described SARS-CoV-2 3CLpro complex with PDB ID: 6Y7M and a resolution of 1.9 Å (Zhang et al., 2020). The complex contains {tert}-butyl ∼{N}-[1-[(2∼{S})-3-cyclohexyl-1-[[(2∼{S},3∼{R})-4-(cyclopropylamino)-3-oxidanyl-4-oxidanylidene-1-[(3∼{R})-2-oxidanylidene-3,4-dihydropyrrol-3-yl]butan-2-yl]amino]-1-oxidanylidene-propan-2-yl]-2-oxidanylidene-pyridin-3-yl]carbamate (OEW), a peptide-like covalent inhibitor (MW = 585.69 g/mol). As described previously (Jukič et al., 2021), this complex comprises the peptidomimetic inhibitor OEW, which occupies all major pockets at the active site of the enzyme, leaving the S1 pocket accessible and the enzyme in the active conformation. After superposition with a reference structure (PDB ID:6LU7), the catalytic binding pocket was defined around Cys145 (Jin et al., 2020). The covalent OEW bond was cleaved, the small molecule removed, and the Cys145 amino-acid residue regenerated (open-source PyMOL, version 2.1; DeLano, 2002). The target was prepared using the protein preparation module of Schrödinger Small-Molecule Discovery Suite (Release Schrödinger 2020-4, Schrödinger, LLC, New York, NY, United States, 2020). Missing hydrogen atoms were added, the H-bond network was optimized using the PROPKA tool at pH 7.4, waters were removed, and restrained minimization was performed with convergence of heavy atoms toward 0.3 Å. Finally, a docking receptor (Morley and Afshar, 2004) was generated with the docking package CmDock (https://gitlab.com/Jukic/cmdock/) using the program cmcavity. The reference ligand method was used to calculate the cavity (receptor definition), where we used the OEW-cleaved regenerated ligand as a reference and a sphere of 7 Å around the ligand to calculate the docking volume. We calculated a total docking volume of 3,106.25 A3 and included the calculated cavity (Cavity #1) in the definition of the docking receptor. The parameters of cavity #1 were size of 24,850 points, min = (−33.5, −53.5, −8.5), max = (−13, −26.5, 12), center = (−24.4138, −38.9632, −0.179235), and extent = (20.5, 27, 20.5) Å (Figure4A).

FIGURE 4

(A): Prepared 3CLpro docking target in the green cartoon model with docking volume highlighted in blue mesh representation. (B): Calculated docking pose of darolutamide shown in the yellow stick model superposed on the OEW reference ligand in the blue stick model. Emphasized are 3CLpro individual binding pockets. (C): Calculated docking pose of nebivolol in the stick model colored magenta. (D): Calculated docking pose of amikacin in the orange-colored stick model. The 3CLpro protein is shown in a green-colored line model, with the active site surface in gray and catalytic Ser144 highlighted in a red-colored stick model.

Structure-based virtual screening

In the virtual screening experiment, we performed molecular docking using CmDock (CmD) software (https://gitlab.com/Jukic/cmdock/; Ruiz-Carmona et al., 2014). The inputs were the FDA-approved drug database precalculated by LigPrep and the prepared receptor (cavity #1) as described in the previous section. First, we performed a redocking experiment in which we successfully redocked a regenerated OEW reference ligand (PDB ID: 6Y7M) and obtained the binding conformation of the crystal complex with a root mean square deviation (RMSD) of 1.34 Å. The CmD parameters were the standard docking protocol (dock.prm) with 100 runs, no constraints, and no score filters. Using the same exhaustive docking protocol, we performed the virtual screening experiment with the 4,711 FDA-approved drugs and analyzed the docking results using KNIME software (Fillbrunn et al., 2017). The minimum docking score was −37.2, the maximum was 8.3, and the mean was -11.7 with a standard deviation of 5.8 after exhaustive docking with CmDock at the prepared receptor binding site (cavity #1).

Inverse docking

We applied an inverse molecular docking approach to 11,640 viral protein structures, including 836 protein structures from the taxonomic lineage of Coronaviridae. Small molecule-binding sites were identified and prepared for inverse molecular docking using the ProBiS-Dock system (Depolli et al., 2013; Konc et al., 2021). By reducing the size of the docking space and focusing on the binding sites, the time and complexity of inverse molecular docking are reduced. The creation of the ProBiS docking database, which served as a template for our viral database, is further described in Štular et al. (2016) and Konc et al. (2021) and has already been successfully used for mechanistic insights into the side effects of troglitazone and rosiglitazone (Kores et al., 2021).

The inverse docking CANDOCK algorithm (Fine et al., 2020) applies a hierarchical approach to small-molecule reconstruction from the atomic lattice using generalized statistical potential functions and graph theory. The docking scores represent an approximation to the relative free energies of binding and have arbitrary units. The algorithm works in several consecutive steps. First, a small molecule is taken and broken into fragments. Then, the fragments are docked into protein-binding sites from the database using knowledge-based scoring methods. Then, the best-docked fragments are selected and linked using a fast maximum clique algorithm (Konc et al., 2021). In the course of the reconstruction, iterative dynamics is used for better placement of the ligand in the binding cavity. In a second step, the conformation optimization procedure is performed (Štular et al., 2016; Fine et al., 2020; Kores et al., 2021). Validation of the algorithm has been extensively reported beforehand (Furlan et al., 2018; Kores et al., 2019; Fine et al., 2020; Kores et al., 2021). To sum up the validation procedures, we employed multiple methods. First is the redocking procedure, reported by Kores et al. (2019) and Kores et al. (2021). In this procedure, the known structures from the Protein Data Bank (PDB) with co-crystallized ligands are taken, and the ligand redocking is preformed to show that the used docking algorithm can produce similar poses with the highest docking score to the native (co-crystallized) ones. This was done for resveratrol (Kores et al., 2019) and also for troglitazone and rosiglitazone (Kores et al., 2021), where good agreement between predicted and experimentally determined poses was shown based on the RMSD of atomic positions. The second method combines the calculation of receiver operating characteristics (ROC), enrichment, and predictiveness curves (PC). The usage was reported by Furlan et al. (2018), Fine et al. (2020), and Kores et al. (2021). Here the experimentally confirmed protein targets from the ChEMBL database for ligands were used to determine if the CANDOCK protocol produced similar binding targets. This analysis was done for curcumin (Furlan et al., 2018) and troglitazone and rosiglitazone (Kores et al., 2021), where it showed that the CANDOCK protocol is expected to provide a good agreement with experiments. Top scoring hits identified by our methodology demonstrate also an experimental activity on 3CLpro (CHEMBL4495582, CHEMBL4495583). Therefore, our proposed protocol and the workflow are suitable for identification of novel potential actives on the studied target.

Inverse docking fingerprinting

We used five reference compounds, namely, chloroquine and hydroxychloroquine, ivermectin B1a and B1b, and selamectin, which are commonly reused drugs as control compounds (Figure 1). In addition, the 11 compounds that scored highest in the virtual screen (Table 1), all with previously reported activity on 3CLpro, and 11 compounds from our virtual screening experiment with the highest score that had not previously been reported to have activity on 3CLpro (CHEMBL4495582, CHEMBL4495583; Table 2) were examined. Using the CANDOCK inverse docking protocol, we elaborated their binding potential against two previously created databases of binding sites (viral protein target database—11,640 targets as well as 836 Coronaviridae protein targets) with a total of more than 325,000 individual docking experiments (Figure 3).

TABLE 1

Identified top-scoring drugs in the virtual screening repurposing experiment on the SARS-CoV-2 main protease 3CLpro with previously reported activity on 3CLpro (CHEMBL4495582, CHEMBL4495583).

No.	Structure	Mr (g/mol)	Name (INN) a	CmDock docking score b	Classification c	Phospholipidosis potential in silico d
1		585.6	Amikacin a	−37.2	Aminoglycoside antibacterial	No
2		451.5	Dibekacin a	−34.5	Aminoglycoside antibacterial	No
3		463.6	Micronomicin a	−32.3	Aminoglycoside antibacterial	No
4		398.9	Darolutamide	−30.2	Anti-androgen (androgen receptor antagonist)	No
5		417.8	Pexidartinib	−28.4	Antitumor agent (selective CSF1R inhibitor)	CAD (Slog P = 5.3, RDKit, pKa > 7.4); not a recorded phospholipidosis inducer
6		444.5	Canagliflozin	−27.9	Antidiabetic (SGLT2) inhibitor)	No
7		522.6	Ticagrelor	−27.7	Antithrombotic (P2Y12 platelet inhibitor)	No
8		405.4	Nebivolol	−27.5	Beta blocker	No
9		341.4	Propafenone	−27.4	Class 1C antiarrhythmic agent	No
10		425.5	Filgotinib	−27.3	Antirheumatic (JAK 1 selective inhibitor)	No
11		427.4	Cariprazine	−26.4	Atypical antipsychotic (D2 and 5-HT1A modulator)	No

Flagged compounds are antibacterials and should be treated with care.

Docking scores in kJ/mol.

Classification resting on ATC, codes.

As per Tummino et al.

TABLE 2

Identified novel drugs in the repurposing experiment with no prior reported 3CLpro activity data.

No.	Structure	Mr (g/mol)	Name (INN)	CmDock docking score a	Classification b	Phospholipidosis potential in silico c
12		386.5	Olodaterol	−25.8	Beta2-adrenergic agonist	No
13		311.3	Xanthinol	−24.5	Vasodilatator	No
14		291.4	Penbutolol	−24.4	Beta-adrenergic antagonist	No
15		418.4	Alloin	−24.1	Not classified, exp.	No
16		430.6	Piritramide	−23.9	Synthetic opioid	No
17		354.5	Sultopride	−23.8	Neuroleptic	No
18		516.7	Udenafil	−23.4	PDE5 inhibitor	No
19		425.7	Doravirine	−22.9	Non-nucleoside reverse transcriptase inhibitor	No
20		418.6	Valbenazine	−22.7	Monoamine transporter 2 inhibitor	No
21		475.6	Arzoxifene	−22.7	Not classified, exp.	CAD (SlogP = 5.5, RDKit, pKa> 7.4); not a recorded phospholipidosis inducer
22		401.5	Risdiplam	−22.4	mRNA splicing modifier	No

Docking scores in kJ/mol.

Classification resting on ATC, codes.

As per Tummino et al.

FIGURE 3

Inverse docking fingerprinting protocol and example of fingerprint comparison. Fingerprinting using heatmaps is useful for determining the similar binding patterns and find targets where compounds bind favorably. The darker the color in the fingerprint heatmaps, the less favorable the binding of the compound into target.

From the calculated conformations of all 27 compounds, we generated ranked lists of CANDOCK docking scores with all targets from the viral protein database or with targets from the Coronaviridae family of proteins (Figure 3). The fingerprint for a single molecule was then generated using the union of the 10 highest-scoring targets from all molecules in the study to obtain comparable fingerprints of equal length. The highest-scoring targets were carefully examined for organisms and protein families. The lists of CANDOCK scores obtained in this way can be used as fingerprints for each compound, and the results are presented as heatmaps for analysis. The authors acknowledge the limitations of the inverse docking method in accurately identifying targets (yet the method is validated as referenced) but postulate that this approach can be used to compare a range of compounds and infer on their structural properties based on in silico interactions with a large number of prepared targets or even pharmacophore models. This simple approach can also pave the way for the development of further in silico selectivity methods (Figure 3).

Results and discussion

Virtual screening of the FDA-approved drug library

To identify accessible hit compounds with potentially favorable physicochemical properties and suitable downstream properties for biological evaluation, rigorous post-docking filtering was performed to screen out pan-assay interference compounds (PAINS; Shoichet, 2006; Baell and Holloway, 2010; Saubern et al., 2011), aggregators (Irwin et al., 2015; to obtain compounds ready for biological evaluation in vitro), and structures with reactive functional groups (REOS; to prioritize toward non-covalent inhibitors; Walters et al., 1998; Zhu et al., 2013). The KNIME software with RDKit nodes was applied to compare all structures in the library with the selection of SMARTS-formatted PAINS, aggregator, and REOS libraries and remove flagged matches. The 10 highest-scoring compounds were selected with the Z-score cutoff of −2.5, clustered, and examined in detail (Table 1).

The identified top-scoring repurposing candidates belong to nine distinct therapeutic classes according to ATC and are detailed in Table 1, with the top three compounds belonging to a group of aminoglycoside antibacterials. Indeed, antibacterial compounds have been proposed as SARS-CoV-2-repurposing candidates, but we would like to highlight their value for the therapy of acquired COVID-19 bacterial coinfections rather than their role as direct-acting antiviral agents (Lai et al., 2020). Moreover, their high scaffold decoration allows for additional target contacts in a typical docking experiment, and we urge the reader to be aware of this fact and interpret the results with this in mind (Meyer-Almes, 2020; Gyselinck et al., 2021). We report that all compounds undergo non-covalent interactions and mainly occupy two different binding poses with classical P1-P2 pocket occupation (Figure 4B; dibekacin, micronomicin, darolutamide, propafenone) or P1′-P2 pocket occupation (Figure 4C), exemplified by pexidartinib, canagliflozin, nebivolol, and filgotinib. The first compound, amikacin, is slightly larger with a MW of 585.6 g/mol and occupies all three P1′-P1-P2 pockets of the active site of 3CLpro. The predicted bound conformations for the first 10 hit compounds are analogous with key contacts Thr25, Leu27, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, Asp187, Thr190, Gln189, and Gln192 at the active site of 3CLpro (Figure 4).

Our top-scoring compounds have also been previously reported and biologically evaluated against 3CLpro (Kumar et al., 2020). Results of extensive screening campaign can also be observed in the ChEMBL database under the comprehensive assay (SARS-CoV-2 3CL-Pro protease inhibition percentage at 20 µM by a FRET kind of response from peptide substrate) with more than 8,700 data points at the time of writing. We report these compounds as 3CLpro hit identification control and support their repurposing research via further in silico inverse docking fingerprint analysis. Nevertheless, we are compelled to report further 11 potential 3CLpro inhibitors (compounds 12–22) identified by our virtual screening experiment without any existing 3CLpro activity data (CHEMBL4495582, CHEMBL4495583) that can be of use to the SARS-CoV-2 repurposing research (Table 2).

New 11 reported hit compounds (12–22) belong to eight distinct therapeutic classes according to ATC with two compounds (15, 21) not yet classified as being experimental agents. Examining the ChEMBL database, no compound possesses activity on coronavirus targets as of yet. In this hit list, we also do not report antibacterials, as they are not favorable for drug repurposing causing additional selection pressure and antibacterial resistance problems (Andersson et al., 2016). The list represents small molecules suitable for further biological evaluation on SARS-CoV-2 3CLpro without phospholipidosis potential as identified by Tummino et al. (2021). Majority of compounds occupy the P1-P2 pocket at the 3CLpro active site near Cys145 (Figure 5A) analogous to the PDB ID: 6Y7M OEW reference ligand. Highest-ranking compound 12 (olodaterol) thus makes hydrophobic contacts with Met165 and Glu166; hydrogen bonds with Gly143, Ser144, Cys145, Glu166, and Gln192; and cation-π interactions toward His41 (Figure 5; additional details in Supplementary Material).

FIGURE 5

(A): Prepared 3CLpro docking target in the green cartoon model with labeled active site residues with calculated docking pose of olodaterol (12) shown in the light-blue stick model superposed on the OEW reference ligand in the white stick model. The binding site surface around the ligand is emphasized in transparent blue-gray color. (B): 2D projection of the calculated docking pose of olodaterol (12) indicating key residues in vicinity and ligand-binding site interactions.

Moody et al. also conducted a virtual screening experiment on 3CLpro and reported a list of 56 top scoring candidates on 3CLpro, selected according to the highest S-score (between −9.95 and −4.24). After biological evaluation, Moody et al. proposed six repurposing candidates with IC50 values ranging from 21.5 to 75.5 µM, namely, micafungin (4th; −9.60 S-score), an antifungal agent of the echinocandin class, ombitasvir (12th; −8.97 S-score), an HCV NS5A inhibitor, boceprevir (23rd; −8.42 S-score), a representative of an HCV protease inhibitor, ivermectin (32nd; −7.74), an antiparasitic that binds to glutamate-gated chloride channels found in invertebrate nerve and muscle cells, tipranavir (35th; −7.48 S-score), an HIV protease inhibitor, and paritaprevir (36th; −7.43 S-score), an HCV NS3-4A serine protease inhibitor (Mody et al., 2021). The observed differences in the hit list can be attributed to the screening software. Moody et al. used the MOE software suite with a different approach to scoring (S-score calculated using the London dG score for placement and the GBVI/WSA dG score for pose refinement). To illustrate the difference, micafungin, an antifungal drug with a molecular mass of 1,270.28 g/mol, was not included, while ombitasvir, boceprevir, ivermectin, tipranavir, and paritaprevir were scored by CmD as −7.81, −7.35, −9.29, −20.10, and −18.27, respectively. However, all reported compounds by Moody et al. are thus reported to have activity on 3CLpro, similarly as we demonstrated by our screening effort.

Inverse docking of repurposing candidates

To identify the protein targets to which our identified hit compounds express the highest binding potential, we performed rigorous post-docking filtering and expressed only top-scoring targets for reader benefit. We recorded the lists of proteins for each compound and considered only the highest-scoring targets (Table 3). Complete target protein lists are found in Supplementary Tables S1–S3. To generate inverse docking target interaction fingerprints, we combined the 10 highest-scoring viral protein targets of each compound into one list (example for amikacin: 2c86A: 0; 2gecB: −51.6166; 2q6fB: −73.2794; 3cl5A: −57.3933; 4f49A: −68.9848; 4h14A: 0; 4pt5A: −59.4445; 4rezA: −51.841; 4wurA: 0; 5c3nA: 0; 5gwzB: −75.8171; 5hyoA: 0; 5jilA: −67.1155; 5nfyA: −60.2988; 6jijA: 0; 6l5tA: −45.6625; 6nozA: −65.0647; 6qfyA: 0; 6u7hB: −39.7457; 6u7kB: −70.3975; 6y3yA: −35.5868; 6zgfC: 0). We removed duplicate proteins and included docking scores for all targets on the list for each compound. A target interaction space was then plotted as a heatmap as shown in Figure 6. The same procedure was used for the Coronaviridae protein subset experiment (Figure 7). Complete lists of successfully docked Coronaviridae family targets for each compound are found in Supplementary Tables S4–S19.

TABLE 3

Identified viral proteins with the highest scores in the inverse molecular docking experiment.

Ligand	Docking score (arbitrary units)	Organism	Protein name
Amikacin	−101.392	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−100.722	Feline immunodeficiency virus (isolate Petaluma)	Retropepsin
	−99.3221	Human rhinovirus type 5	Rhinovirus B5 VP4
	−98.5941	Human immunodeficiency virus type 1 (BRU ISOLATE)	HIV-1 protease
Canagliflozin	−71.9432	Paramecium bursaria Chlorella virus PBCV-1	Probable thymidylate synthase
	−64.0748	Human immunodeficiency virus type 2 (ISOLATE ROD)	Protease
	−62.3561	Human immunodeficiency virus type 1 (BRU ISOLATE)	HIV-1 protease
	−61.7061	Southampton virus (serotype 3)	Thiol protease P3C
Cariprazine	−77.0656	Human immunodeficiency virus type 2, HIV-2	HIV-2 protease
	−75.9674	Human immunodeficiency virus-1 HIV-1	Protease
	−74.7457	Human immunodeficiency virus type 1 (BRU ISOLATE)	Protease
	−74.6709	Human immunodeficiency virus type 1 (ARV2/SF2 ISOLATE)	Protease
Chloroquine	−69.6592	Paramecium bursaria Chlorella virus PBCV-1	Probable thymidylate synthase
	−62.3342	Human enterovirus EV68	Capsid protein VP1
	−61.6339	Human enterovirus EV68	Viral protein 1
	−61.5787	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
Darolutamide	−73.0195	Human rhinovirus type 5	Rhinovirus B5 VP4
	−71.54	Human immunodeficiency virus type 1 (BH10 ISOLATE)	Reverse transcriptase/ribonuclease H
	−70.4064	Dengue virus type 1 Singapore/S275/1990	Fusion protein of nonstructural protein 2B and nonstructural protein 3
	−69.5576	ZIKV	NS2B-NS3 protease
Dibekacin	−82.4989	Human immunodeficiency virus type 1 (BRU ISOLATE)	HIV-1 protease
	−82.4278	Human immunodeficiency virus-1 HIV-1	Protease
	−82.1087	Adeno-associated virus type 3B	Capsid protein VP1
	−81.8136	Influenza virus type A	Hemagglutinin HA2 chain
Filgotinib	−76.1691	Human rhinovirus type 5	Rhinovirus B5 VP4
	−71.0266	Human immunodeficiency virus type 1 (BH10 ISOLATE)	Reverse transcriptase/ribonuclease H
	−70.9686	Human immunodeficiency virus type 1 (NEW YORK-5 ISOLATE)	Capsid protein p24
	−69.6959	Human immunodeficiency virus-1 HIV-1	gp120
Hydroxychloquine	−75.4441	Human enterovirus CVA10	Capsid protein VP1
	−70.6169	Human enterovirus CVA10	Capsid protein VP1
	−70.4093	Human enterovirus CVA16	VP1
	−69.3791	Human immunodeficiency virus-1 HIV-1	Protease retropepsin
Ivermectin B1a	−111.313	Human immunodeficiency virus-1 HIV-1	Protease
	−108.816	Human immunodeficiency virus type 1 (Z2/CDC-Z34 ISOLATE)	Protease
	−104.269	Human immunodeficiency virus-1 HIV-1	Protease
	−104.208	Human immunodeficiency virus-1 HIV-1	Protease
Ivermectin B1b	−107.924	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−106.019	Human immunodeficiency virus-1 HIV-1	Protease retropepsin
	−105.357	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−104.804	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
Micronomicin	−88.8098	Human immunodeficiency virus-1 HIV-1	POL polyprotein
	−88.7453	Human immunodeficiency virus-1 HIV-1	Protease
	−86.5514	Adeno-associated virus type 3B	Capsid protein VP1
	−86.012	Human immunodeficiency virus-1 HIV-1	Protease
Nebivolol	−71.5908	Human immunodeficiency virus-1 HIV-1	HIV-1 capsid protein
	−71.0465	Human enterovirus EV68	Viral protein 1
	−70.2915	DG-75 Murine leukemia virus	Gag-pro-pol polyprotein
	−69.8426	Human immunodeficiency virus-1 HIV-1	Protease
Pexidartinib	−61.0336	Adeno-associated virus type 6	Capsid protein VP1
	−54.6906	Human immunodeficiency virus-1 HIV-1	Protease
	−54.0666	Bovine respiratory syncytial virus BRSV	Fusion glycoprotein F0, Fibritin
	−53.3484	Human immunodeficiency virus type 1 (NEW YORK-5 ISOLATE)	Capsid protein p24
Propafenone	−69.6736	Human enterovirus EV68	Capsid protein VP1
	−67.4871	Human enterovirus EV68	Viral protein 1
	−67.0906	Influenza virus type A	Hemagglutinin HA2 chain
	−66.5011	Human immunodeficiency virus type 1 (BRU ISOLATE)	Protease
Selamectin	−89.903	ZIKV	NS2B-NS3 protease
	−86.8758	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−85.6238	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−85.1895	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
Ticagrelor	−83.1521	Human immunodeficiency virus type 2, HIV-2	HIV-2 protease
	−79.2021	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−78.8308	Human immunodeficiency virus-1 HIV-1	HIV-1 protease
	−78.4946	Human echovirus E11	Echovirus 11 coat protein vp1

FIGURE 6

Heatmap representation of compound fingerprints for the complete viral target database (11,640 proteins) of all hit compounds. The total number of combined targets is 144, and the list of PDB IDs is found in Supplementary Table S20. The values shown in the heatmap were calculated and colored according to the normalized docking score for each compound. Interval 0 means that the compound did not interact with the target, while interval 5 means that the compound had a most favorable docking score. Protein targets are grouped according to their class and classes emphasized in blue color.

FIGURE 7

Heatmap representation of compound fingerprints for the Coronaviridae database subset (836 proteins) for all compounds. The total number of combined targets is 96, and the list of PDB IDs is found in Supplementary Table S21. The values shown in the heatmap were calculated and colored according to the normalized docking score for each compound. Interval 0 means that the compound did not interact with the target, while interval 5 means that the compound had a most favorable docking score. Protein targets are grouped according to their class and classes emphasized in blue color. Fingerprint profile similarities are colored green. Two “repurpose” subsets represent two different fingerprint profiles within the repurpose compounds group.

An examination of our inverse docking fingerprints reveals that they are unique to each compound, suggesting that subtle differences in compound conformational space and scaffold decoration (chloroquine and hydroxychloroquine, ivermectin B1a and B1b) have a profound effect on the scope of the target reach, making them very useful for the future development of models (Belyaeva et al., 2021; Kumar Das et al., 2021) and machine learning (D’Souza et al., 2020). Namely, the compound’s ability to conform to a particular set of protein targets is used as ligand structural information that can be of use also when comparing a set of diverse compounds (e.g., can two compounds, even structurally dissimilar conform to a larger set of protein targets and produce favorable binding modes?)

The most studied target of SARS-CoV-2 is the 3CLpro, and it is actually postulated as a target (Mody et al., 2021) for the previously reported repurposed drug ivermectin (Kaur et al., 2021; Kern et al., 2021; Nardelli et al., 2021). As is shown in Table 3, the major targets of ivermectin B1a and B1b, as well as the structurally similar selamectin, were identified as proteases. We can support previous reports of ivermectin targeting (Eweas et al., 2021), while with our approach we are able to also place the compounds in the context of other potential therapeutic targets (Zaidi and Dehgani-Mobaraki, 2021; Figure 6). The fingerprint profiles of ivermectins and selamectin are also similar (Figure 6). In comparison, the targets of chloroquine and hydroxychloroquine, calculated by inverse docking, are not proteases and their fingerprints differ significantly. Following the initial reports (Garcia-Cremades et al., 2020) on their efficacy and viral targeting (Adeoye et al., 2020), we are unable to report similar observations and we can treat these compounds (and their fingerprint profile)as negative control examples, as recently reported (Shah, 2021). The fingerprints described can thus be of use in an extensive inverse docking fingerprint screening of a large database of compounds to identify compounds with a similar fingerprint profile.

Upon inspection of the inverse docking fingerprints and top docking poses for all hit compounds 1-11 from the first virtual screening experiment with FDA-approved drugs, all drugs except filgotinib have high docking potential against proteases (Table 3). Beforehand, Ribone et al. (2021) identified protease (PLpro) as the main target of hydroxychloroquine, while further suggesting a human protein PIKfyve as a main target of pexidartinib. Therefore, we can support the potential of the identified FDA-approved drugs to inhibit SARS-CoV-2 3CLpro and also treat the compounds as controls (experimental 3CLpro activity; CHEMBL4495582, CHEMBL4495583). For example, we can postulate that the diverse set of compounds can produce favorable binding poses at the active site of HIV proteases so this set of compounds can access the conformational space that is similar or presents a similar set of transient pharmacophoric elements. This identified potential (similarity in accessing other targets) of a set of studied compounds can then be leveraged in further experimental and repurposing experiments to access novel chemical space.

We fingerprinted the compounds 1–22 using the Coronaviridae family subgroup in a similar manner. The proteins of compounds 1–11 controls (experimental 3CLpro activity) are used for fingerprint comparison with novel repurposing candidates without prior experimental information 12–22 (Figure 7). The highest-scoring Coronaviridae targets belong to the protease family (Supplementary Tables S4–S19), making all newly reported compounds except arzoxifene, a potential SARS-CoV-2 3CLpro inhibitor that should be evaluated further (Figure 7).

If we compare the fingerprint of propafenone to the other repurposing candidates, we can see that it has significantly more unfavorable targets than other compounds (Figure 7). We also observed that propafenone possesses a distinct fingerprint profile in comparison to other studied drugs. From this, one can postulate that propafenone could possess different mechanisms of action to the other compounds. Each fingerprint is unique; nevertheless, global trends emerge when examining multiple compounds in parallel (Figure 7). Namely, pexidartinib and arzoxifene show similar fingerprint profiles (Figure 7; bottom) which indicates less favorable binding to majority of targets, and both are incidentally also flagged as CAD compounds with phospholipidosis potential (Tables 2, 3). This information can be leveraged in further full-scale inverse docking screening, as stated beforehand.

To further elaborate on inverse docking fingerprints, two case studies were conducted. First is a case study of a fingerprinting target, papain-like protease (PLpro). This target was selected as a random representative of the most abundant protease class. Case compounds were also randomly selected, one from each fingerprint-like group: dibekacin, hydroxychloroquine, olodaterol, pexidartinib, and udenafil. As detailed in Supplementary Figure S1, all five compounds bind into the same binding pocket, validating the target database (DeLano, 2002, Adasme et al., 2021; Pettersen et al., 2004). More detailed analysis of the binding interactions shows that udenafil has hydrophobic interactions with Asp165 and Tyr269; hydrogen bonds with Leu163, Tyr265, Tyr269, and Tyr274; and salt bridges to Asp165 and Asp303, suggesting that it may have a more favorable binding potential as a repurposing candidate compared to pexidartinib, which shows hydrophobic interactions with Leu163, Gln270, and π-π stacking with Tyr269. This observation is also immediately apparent from the fingerprint profile in Figure 7. Second, a case study of the fingerprint compound olodaterol is presented. This compound was chosen because it is the highest-scoring repurposing candidate for which data on 3CLpro activity are not yet available. We investigated the binding positions of the repurposing candidate to different protein targets included in its inverse docking fingerprint. The representatives of the target proteins were randomly selected from each protein target class: 3C-like protease and papain-like protease, NSP16 transferase, and spike protein. A detailed analysis of the binding interactions (Supplementary Material) shows that the repurposing candidate has the best docking results and the most favorable binding conformations toward proteases compared to the other targets (3CLpro: hydrogen bonds with Leu141, Gly143, Ser144, Cys145, and Glu166; PLpro: hydrophobic interactions with Asp165, Tyr265, and Tyr269; and hydrogen bonds with Asp165, Arg167, Glu168, Tyr269, and Gln192). The observations thus support our findings and suggestions on the new SARS-CoV-2 repurposing candidates.

Conclusion

We conducted a thorough virtual screening experiment on SARS-CoV-2 3CLpro using CmDock software and a database of FDA-approved drugs and identified the highest-scoring drug repurposing candidates. All newly reported candidates (12–22) are readily available and show favorable non-covalent interactions at the active site of 3CLpro, and all newly identified hit compounds show a low propensity for phospholipidosis, with the exception of arzoxifene. The latter compound belongs to a typical cathionic amphiphilic drug scaffold but has not been identified as a known inducer of phospholipidosis. Furthermore, in reviewing the inverse docking fingerprints, we found that the majority of the identified FDA-approved repurposing candidates have favorable docking scores against Coronaviridae family proteases and viral proteases in general. In addition, we support current in silico studies on 11 top-scoring compounds (1–11) as well as on antihelminthic class member compounds. We postulate the application of this approach to future inverse docking fingerprint screening experiments to investigate the selectivity of compound interaction in silico and to identify similarly interacting compounds in large protein databases.