In Silico Analysis of the Multi-Targeted Mode of Action of Ivermectin and Related Compounds

Abstract

Some clinical studies have indicated activity of ivermectin, a macrocyclic lactone, against COVID-19, but a biological mechanism initially proposed for this anti-viral effect is not applicable at physiological concentrations. This in silico investigation explores potential modes of action of ivermectin and 14 related compounds, by which the infectivity and morbidity of the SARS-CoV-2 virus may be limited. Binding affinity computations were performed for these agents on several docking sites each for models of (1) the spike glycoprotein of the virus, (2) the CD147 receptor, which has been identified as a secondary attachment point for the virus, and (3) the alpha-7 nicotinic acetylcholine receptor (α7nAChr), an indicated point of viral penetration of neuronal tissue as well as an activation site for the cholinergic anti-inflammatory pathway controlled by the vagus nerve. Binding affinities were calculated for these multiple docking sites and binding modes of each compound. Our results indicate the high affinity of ivermectin, and even higher affinities for some of the other compounds evaluated, for all three of these molecular targets. These results suggest biological mechanisms by which ivermectin may limit the infectivity and morbidity of the SARS-CoV-2 virus and stimulate an α7nAChr-mediated anti-inflammatory pathway that could limit cytokine production by immune cells.

1. Introduction

The spread of the SARS-CoV-2 pandemic around the world has spurred on a search for suitable drugs for therapeutic applications against this viral infection. Although vaccination is a proven strategy for containing coronavirus disease 2019 (COVID-19), a new challenge has developed due to the emergence of SARS-CoV-2 variants for which vaccines have offered lesser degrees of protection. Efforts have therefore been focused on the possibility repurposing existing, approved drugs, which do not require de novo design and lengthy testing.

Obstruction of the binding between SARS-CoV-2 spike protein and the ACE2 receptor on target human cells has been one focus of therapeutic intervention [1,2]. Although viral fusion to host cells and replication occur via ACE2, SARS-CoV-2 and several other viral strains, including other betacoronaviruses, initially attach to host cells via more abundantly-distributed glycoconjugate host binding sites [3,4]. Notable among these are sialic acid (SA) and the transmembrane glycoprotein receptor CD147, which are widely distributed in blood, endothelial and several other cells in the human body [4] and provide additional therapeutic targets. Nicotinic acetylcholine receptors, which are densely distributed in neuronal tissue and in cytokine-secreting inflammatory cells (for example macrophages and mast cells), may provide an additional attachment point for SARS-CoV-2 due to a “toxin-like” epitope on the viral spike protein and may thus provide druggable targets as well [5,6,7].

1.1. Binding of SARS-CoV-2 Spike Protein to Host Cell SA and CD147 Surface Molecules

The SARS-CoV-2 spike protein contains 22 N-linked glycosylation sites on each of its three monomers, with several of these associated glycans capped with terminal SA moieties in various forms [8,9,10,11]. Through those glycan bindings, the SARS-CoV-2 spike protein attaches to an SA-coated nanoparticle array, which is the basis for a sensitive viral detection technique [12]. Such a nanoarray is paralleled in densely distributed SA-tipped binding sites of glycophorin A molecules [13,14,15] and CD147 receptors [16,17,18] at the surface of red blood cells (RBCs), and SARS-CoV-2 was observed in the hemadsorption assay to clump with human RBCs [19]. Clinical confirmation of viral-RBC attachments was provided by the presence of SARS-CoV-2 spike protein punctae on 41% of RBCs from a series of hospitalized COVID-19 patients [20]. Binding between SARS-CoV-2 spike protein and the host cell receptor CD147 was likewise demonstrated by SPR, Co-IP and ELISA assays, and immuno-electron microscopy [21]. Meplazumab, a humanized anti-CD147 antibody, inhibited SARS-CoV-2 replication in vitro and reduced time to viral clearance from 13 to 3 days for COVID-19 patients in a small clinical study [21,22].

For viruses that bind to SA, including SARS-CoV-2, as noted above, that binding plays a key role in viral infectivity, as SA on host cells typically serves as the initial attachment point for viral spike protein [23,24,25,26,27,28,29]. But for such SA-binding viruses, the host limits viral attachment to host cell infectious targets through other entities in the body having SA-rich surfaces, including RBCs, platelets and leukocytes, along with mucins and plasma proteins [30,31]. Some viruses, in turn, dodge that defense through the expression of SA-cleaving enzymes that enable detachment from these blood cells and other snagging substances [28,31,32,33,34,35,36,37,38]. In particular, the SA-cleaving enzyme hemagglutinin esterase (HE) is expressed by the human betacoronaviruses that cause the common cold, OC43 and HKU1, but not by SARS-CoV, SARS-CoV-2 and MERS, the three deadly strains in that viral family [39,40,41]. It has been proposed that vascular occlusion, central to the morbidities of COVID-19 [42,43,44,45,46,47], is initially triggered by the clumping and snagging of SARS-CoV-2 with blood and endothelial cells, and that HE expressed by the common cold betacoronaviruses may limit these morbidities [4].

1.2. The Role of CD147 in the Inflammatory Response

In addition to its SARS-CoV-2 binding capability, CD147, a transmembrane glycoprotein receptor encoded in humans by the BSG gene [48], is of interest as a key mediator of inflammatory response, in particular, as related to vascular occlusion. In response to immunogenic stimuli, CD147 is upregulated in T cells [49,50], platelets [51,52] and endothelial cells [53], with upregulation of CD147 in endothelial cells occurring upon exposure to active or UV-deactivated betacoronavirus MHV-4 in vitro. CD147, in turn, has been observed to promote adhesion by RBCs [54,55,56], leukocytes [52,57,58,59] and platelets [57,58,60] to other blood cells and endothelial cells. Also, of particular interest are the indicated pro-infectious roles of CD147 and its binding partner cyclophilin A for SARS-CoV-2, SARS-CoV and other viruses [21,61,62,63]. In a broader clinical framework, the involvement of CD147 in the pathogenesis of a number of diseases, including lung inflammation, atherosclerosis, heart failure, ischemic myocardial injury and stroke [52,56,58,64,65], further suggests that CD147 antagonists or masking agents could mitigate a COVID-19 infection.

1.3. Competitive Binding of Ivermectin to SARS-CoV-2 Spike Protein Binding Sites

Given that the attachment of the SARS-CoV-2 spike protein to host cell targets, including ACE2, SA and CD147, is central to viral infectivity and morbidity, the capability for competitive binding to limit such attachments has been one focus in the search for repurposed COVID-19 therapeutics [3]. Four molecular modeling studies that collectively screened over 800 such molecules were conducted toward that goal [66,67,68,69]. The strongest or close to strongest binding affinity in each study was obtained for ivermectin, a macrocyclic lactone with multifaceted antiparasitic and antimicrobial activity which has been distributed in 3.7 billion doses worldwide since 1987 [70,71,72,73]. Additional molecular modeling studies of competitive binding to SARS-CoV-2 spike protein sites that focused on ivermectin in particular likewise found strong binding affinities for that agent [74,75,76,77,78,79].

These findings are of interest given clinical, animal and epidemiological studies, including most of the 20 randomized clinical trials (RCTs) conducted to date, indicating the efficacy of ivermectin against COVID-19 [70,80,81], although interpretations of which of these RCTs are most reliable have been controversial. Ivermectin is suitable for mass use on a global scale, having been the mainstay of two worldwide campaigns to eliminate two devastating scourges affecting millions, onchocerciasis and lymphatic filariasis [82]. It is safe even at much higher doses than the standard dose of 200 μg/kg [83,84], and its limited side effects were noted in the Nobel Committee’s 2015 award honoring its discovery and its record of improving the health and wellbeing of millions [85]. However, a biological mechanism initially proposed for ivermectin activity against SARS-CoV-2, entailing blockage of its transport into the host cell nucleus, was proposed in conjunction with in vitro studies conducted at much greater than physiological concentrations and has been questioned [86,87,88].

1.4. Nicotinic Acetylcholine Receptors: Anti-Inflammatory Modulation and Blockage of Viral Bindings

Another biological mechanism of activity that may underlie the observed clinical benefits of ivermectin treatment of COVID-19 is a potent anti-inflammatory and immune modulatory effect mediated by its action as a positive allosteric modulator of the alpha-7 nicotinic acetylcholine receptor (α7nAChr) [89]. The core receptor of the cholinergic anti-inflammatory pathway is α7nAChr, which is under the control of the vagus nerve [90] and plays a crucial role in balancing of the body’s response to inflammation and sepsis [90,91]. This anti-inflammatory pathway connects the involuntary parasympathetic nervous system innervating all major organs to cytokine-producing cells such as TNF, IL1 and IL6-secreting macrophages, lymphocytes and mast cells [90,91], which are reported to play a major role during the inflammatory phase of COVID-19 infection (i.e., the cytokine storm [92]). The ivermectin-induced enhancement of this pathway might rapidly lower pro-inflammatory cytokine levels and decrease expressions of chemokines as well as adhesion molecules at the inflammatory sites [90,91]. Importantly, the marked increase in Ca++ current evoked by acetylcholine (ACh) in the presence of micromolar concentrations of ivermectin (e.g., a 20-fold shift of the affinity of ACh [89]) may also potentially explain the reported clinical activity of ivermectin during the late (i.e., inflammatory), critical phase of severe COVID-19 cases [93].

Recent in silico docking studies have indicated a potential direct interaction between the SARS-CoV-2 spike glycoprotein and α7nAChr, due to a “toxin-like” epitope on the spike glycoprotein, with homology to a sequence of a snake venom toxin [5,6]. Of interest, the α7nAChr receptor, which is densely distributed on neuronal tissue, has previously been shown to serve as the port of entry in the human body for another RNA virus endowed with strong neurotropic action, the rabies virus [7]. The loss of smell (anosmia) and/or taste (ageusia) are considered hallmarks of COVID-19 infection and are likely consequences of the direct SARS-CoV-2 infection of the olfactory and gustatory nerve [94]. Ivermectin high affinity binding to α7nAChr may therefore interfere with the attachment and internalization of SARS-CoV-2 on the olfactory/gustatory nerves, as recently reported in both animal models [94] and human patients [95].

1.5. Subdomains of the SARS-CoV-2 Spike Protein, S1 Region

The SARS-CoV-2 spike protein pockets selected for study as potential ivermectin binding sites were governed by the arrangement of subdomains of interest. The SARS-CoV-2 spike protein is a heavily glycosylated, type I transmembrane protein with 1,273 amino acid residues, assembled into trimers and attached on the virion surface, giving the virus its distinctive “corona” or crown-like appearance. Each spike protein trimer contains a central helical stalk consisting of three joined S2 subunits, each capped with an S1 subunit head in a mushroom-like shape [96,97]. The ectodomain of the spike protein S1 attaches to a host cell membrane, after which the S2 stalk engages in fusion, enabling the internalization and replication of the virus [96,97,98,99].

That viral-host engagement proceeds, in particular, through fusion of the receptor binding domain (RBD) of the SARS-CoV-2 spike protein S1 to an ACE2 receptor on a host cell [97,98,99]. As shown in Figure 1, the S1 N-terminal domain (NTD) contains eight of the 22 N-linked glycans on the SARS-CoV-2 spike protein. These eight N-linked glycans on the NTD form initial attachments to host cells’ glycoconjugates, including CD147 and others which have SA terminal residues [8,18,100,101,102,103,104]. The RBDs of the virus, one on each spike protein monomer S1 subunit, switch constantly between open (“up”) and closed (“down”) configurations, with the former enabling both ACE2 binding and immune surveillance and the latter blocking both of those functions [96,105].

2. Materials and Methods

2.1. Ligand Database Preparation

A set of 15 test compounds, composed of ivermectin and 14 similar molecules, was collected from the PubChem database and used for the docking studies. Compounds structurally similar to ivermectin were adopted from DrugBank Online under the “similar structures” section of the drug ivermectin [107]. This option provides users the capability to search rapidly for structurally similar small molecules, without having to redraw the molecule and perform additional database searches through the ChemQuery interface. Before evaluating any interaction, the ligand set database was prepared through a “Wash” wizard of the Molecular Operating Environment (MOE) software package. At this stage, the 3D dominant protonation state of each molecule was generated at the physiological pH of 7, followed by a short MOE built-in energy minimization procedure.

2.2. Protein Preparation

The crystal structure of CD147 was obtained from its Protein Data Bank (PDB: 3B5H). Chain A has the strongest electron density and thus was used for analysis. The CHARMM-GUI Archive of COVID-19 Proteins Library [108] was used to collect the structures of the two spike protein conformations, i.e., the closed (PDB: 6VXX) and open (PDB: 6VSB) states [9,109,110]. The NTD (aa 18–292) and RBD (aa 318–513) of one monomer were considered separately in the following analysis. The atomic coordinates of three possible conformations of α7nAChr, namely resting (PDB: 7KOO), desensitized (PDB: 7KOQ) and activated (PDB: 7KOX), were obtained from the PDB [111]. Only the extracellular region of the protein (aa 1–207) was considered in the docking analysis. All proteins were prepared in MOE, adjusting their protonation state according to a physiological pH of 7 and minimizing the potential energy.

2.3. Binding Sites

We employed the Site Finder module in MOE [112] to detect the possible binding sites in the NTD and RBD domains of the spike protein. All the sites we identified using the MOE software had already been reported in the literature, as we summarized them below, so we performed our molecular modeling calculations using the sites specified in

Table 1. Binding sites of spike NTD and RBD obtained from the literature (as reproduced from Table 3 in Aminpour et al. [3], with reference citation numbers adjusted to this reference list).

Binding Site	Reference	Binding Site Type	Residues	NTD/RBD
Site 1	Milanetti et al. [103]	sialoside	L18-Q23, H66-T78 and G252-S254	NTD
Site 2	Behloul et al. [113]	sialoside	E154, F157, Y160 and the so-called stabilizing loop (N122-N125)	NTD
Site 3	Baker et al. [12]	sialoside	(R21, Q23, L24, H69, F79, P82 and R246)	NTD
Site 4 (P1)	Di Gaetano et al. [114]	sialoside	R21, T22, Q23, L24, P26, R78, P82, V83, L110, F135, C136, N137 and R237	NTD
Site 5 (P2)	Di Gaetano et al. [114]	sialoside	F92, S94, E96, K97, S98, R102, N121, V126, I128, M177, D178, K182, N188, R190, F192, I203, L226, V227 and L229.	NTD
Site 6–14	Watanabe et al. [11]	glycosylation	N122, N149, N165, N17, N61, N74, N234, N282	NTD
Site 15	Fantini et al. [101]	ganglioside	Domain (111–158)- core Q-134 to D-138	NTD
Site 16	Carino et al. [117]	-	F342 N343 A343 T345 R346–W436 N437 S438–L441 D442 S443–G446–N448–Y451 L452	RBD
Site 17	Carino et al. [117]	-	S375–G404 D405–V502 G503–Q506–Y508	RBD
Site 18	Carino et al. [117]	-	E340 V341–F347 A348–N354 R355 K356–S399 F400 V401–V512	RBD
Site 19	Carino et al. [117]	-	F374–N388–Y495 G496 F497	RBD
Site 20	Carino et al. [117]	-	T376 F377 K378 C379 Y380–V407 R408–I410–V433 I444 A445	RBD
Site 21–22	Watanabe et al. [11]	glycosylation	N331–N3443	RBD

. We manually calculated the center of the binding sites from the residues involved.

Several sialoside-, glycosylation- and ganglioside- binding sites have been reported in the literature. Milanetti et al. [103] proposed a potential sialoside binding site containing three divergent loop regions (site 1). They supported the hypothesis of a structural resemblance between MERS-CoV and SARS-CoV-2 using iso-electron density mapping. Behloul et al. [113] compared the structural features of the SARS-CoV-2 spike protein S1-NTD with BCoV and consequently characterized a binding pocket that has the capability to bind SA species such as Neu5,9Ac2 (site 2). Baker et al. [12] aligned the sequences of the SARS-CoV-2 spike protein, mainly focusing on human coronavirus OC43 as the SA-binding protein. They identified a potential SA binding site, associating its glycan-binding characteristic utilizing glyconanoparticles for the detection (site 3). Gaetano et al. [114] calculated the druggability of all available ligand-binding pockets within the NTD segment of the spike protein S1 using SiteMap of Schrodinger software [115]. As a result, among all of the three hypothesized sialoside-binding pockets in the literature, site 3 by Baker et al. [12] is part of a cavity with a druggable property identified by Gaetano et al. (site 4 in Table 1, or site P1 as Gaetano et al. referenced in their paper [114]). Gaetano et al. also identified an unexpected binding pocket (site 5 in Table 1, or P2 as they referenced in their paper) within S1-NTD. Site 5 (P2) aligns with the recent experimental findings by Bangaru et al. [116].

Sites 6 to 14 are associated with the glycosylation binding sites proposed by Watanabe et al. [11]. Fantini et al. [101], meanwhile, proposed a new type of ganglioside-binding domain performing molecular dynamics (MD) calculations. The results of his simulations reveal a strong interaction between GM1 ganglioside and S1-NTD (site 15). Finally, Carino et al. [117] utilized the Fpocket server (https://bioserv.rpbs.univ-paris-diderot.fr/services/fpocket/, accessed 22 February 2022) and computationally identified sites 16 to 20 in the RBD fragment of the spike protein. They also studied the binding of several triterpenoids (e.g., glycyrrhetinic and oleanolic acids) and natural bile acids and demonstrated that their semisynthetic derivatives can reduce RBD adhesion to ACE2 in vitro. Sites 21 to 22 belong to the set of glycosylation binding sites proposed by Watanabe et al. [11].

Nine different binding sites in the CD147 dimer identified by MOE (Site Finder) are listed in Table S1 and are illustrated in Figure S1. The three potential N-glycosylation sites of CD147 are N44 (site 8), N152 (site 3) and N186 (site 3).

The putative binding pockets of the α7nAChr (desensitized, activated, resting) and CD147 structures were identified using the Site Finder tool in MOE, which computes the possible binding sites of a protein from its 3D structure using a geometrical approach. For each protein, only the sites characterized by a propensity for ligand binding (PLB) greater than or equal to 1 were considered in the docking experiments. The α7nAChr protein is a pentamer with a five-fold symmetry. To have a clear presentation, we summarized the common binding sites between monomers of the α7nAChr protein, excluding the common ones, in Table S2. In total, 37 binding sites with PLB > 1 were identified in all the three conformations of the α7nAChr protein. The binding sites identified by MOE (Site Finder) related to Table S2 are illustrated in Figure S2.

2.4. Molecular Docking Simulations

Docking was performed with a flexible ligand and a rigid receptor approach using the AutoDock Vina program [118] to predict the binding pose of the ligands. In the AutoDock Vina software, receptor–ligand binding affinities were predicted as negative Gibbs free energy (∆G) values (kcal/mol), which were calculated on the basis of the AutoDock Vina scoring function and classified on the basis of a numerical value referred to as the “Score”. The interactions of inhibitors with receptor proteins are predicted on the basis of the Score; the lower the Score (in negative value), the greater the interaction. The Vina scoring function incorporates two features from knowledge-based and empirical potentials. A cubic box with 30.0 Å size, required to delimit the docking area, was used on each binding pocket, centered at their center of geometry. The maximum number of poses to be generated for each docking calculation was set to 20. The minimum root mean square deviation (RMSD) to distinguish between two different poses was 1 Å. Every generated pose was energy-minimized in vacuo using Amber16 by keeping the protein fully rigid [119], with out of box poses then being discarded. Finally, the Vina Score function was used to re-score the poses after the minimization and the pose with the best Score was selected for each compound–receptor pair. The DockBox package was used to facilitate the preparation of docking inputs, the post-processing of the docking results and the rescoring procedure [120]. No constraints were applied in the docking studies. Although we minimized the ligand–protein structures after docking, we double checked the stability of compounds by running 100 ns MD calculations in explicit solvent on the unrestrained ligand–protein complex (see Section 2.5).

To the best of our knowledge, there are no effective therapeutics for COVID-19 which have biological mechanisms similar to those indicated for our 15 test compounds to compare with our docking results, which could be checked for competitive binding to the spike protein or the other host receptors. Therefore, we have the limitation of not being able to usefully check these results against known controls. In order to evaluate docking parameters for a given target prior to undertaking docking calculations on unknown ligands, however, it is always beneficial to perform control docking if the binding of known ligands is available in the crystal structure and if they have a non-covalent nature. Therefore, here, we also performed positive control docking calculation for the ligands that were experimentally available in the crystal structure of the proteins that we used in our study. It was not possible to use the NAG (N-acetyl-D-glucosamine) ligand (PDB: 6VSB) of the spike protein as a positive control since the nature of the binding was covalent. Also, we were not able to perform control docking on the CD147 protein (PDB: 3B5H) since there was no known ligand available in the crystal structure. The ligand Epibatidine (PDB:7KOX) of the alpha-7 nicotinic acetylcholine receptor was used as a control. We were able to successfully generate the same pose (RMSD = 1.2 Å) with a binding affinity of −8.73 kcal/mol (see Figure S3). We used decoys, which are molecules that are physically similar yet chemically dissimilar to the active ligands [121], as a negative control for the docking calculations. We used a state-of-the-art benchmark, the Directory of Useful Decoys (http://dude.docking.org, accessed 22 February 2022), to select decoys for ivermectin [122,123]. The structures of the decoy compounds are presented in the Supplementary Information (Figure S4). The binding affinities of the decoy compounds for the spike protein S1, CD147 and α7nAChr binding sites are in the range of (−3.345 to −5.496 kcal/mol), (−4.217 to −5.137 kcal/mol) and (−4.940 to −6.070 kcal/mol), respectively. The decoy compounds exhibited lower affinities than ivermectin and the most related compounds.

2.5. Molecular Dynamics (MD) Simulations

In order to establish the stability of each docked protein-inhibitor complex, MD simulations were run in explicit solvent using the Amber16 software. The computationally intensive all-atom MD simulations of the Open (6x) and Closed (3x) systems (each tallying ∼1.7 million atoms including explicit water, ions and membrane lipids) that was done by another research group on Texas Advanced Computing Center (TACC) achieved benchmarks of ∼60 ns/day on 256 GPU nodes [124]. To reduce the computation time, for NTD binding ligands, protein spike S1 was truncated from S698 to D1146, and from P322 to C590. For RBD binding ligands, protein spike S1 is truncated from M1 to E324 and from C590 to D1146. The hydrophobic part of the α7nAChr protein (T207 to L320) was removed in each monomer to prevent the exposure of the hydrophobic area in water. The breaks in all of the structures were capped with MOE’s Structure Preparation. All the residues of protein CD147 were kept. MD simulations were carried out on Compute Canada’s Graham cluster (V100 GPUs), as well as Cedar (P100 NVIDIA GPUs), depending on their respective availability. Each simulation was carried out on a single GPU. Using the AmberTools 16′ leap program, each complex was solvated in a cubic box with a side length of 12 Å using a three-points (TIP3P) water model. Na⁺ and Cl⁻ ions were added in such a way to adjust the salt concentration to the physiological value of 0.15 M and neutralize the system. The minimization of the complexes was achieved in two steps, using the steepest descent (5000 steps) and conjugate gradient (5000 steps) methods successively. At first, only solvent atoms were minimized, by restraining the protein–ligand complex. Next, the minimization was run with the same parameters without the restraint. After the minimization step, the MD simulations were conducted in three stages: heating, density equilibration and production. At first, each solvated system was heated to 298 K for 500 ps, with weak restraints on all backbone atoms. Next, density equilibration was carried out for 1 ns of constant pressure equilibration at 298 K, with weak restraints. Finally, MD production (one trajectory per complex) were performed without any restraints for all systems for 100 ns. The trajectory of the ligand–protein complex was visually investigated using the VMD package (the University of Illinois at Urbana-Champaign, Urbana, IL, USA). Time-evolutions of the RMSD of top-ranked inhibitors with respect to receptors (spike, CD147 and α7nAChr) were calculated using the CPPTRAJ module of the AMBER16 software. Clustering analysis was carried out on the protein-bound ligand poses where the trajectory reached a plateau using Amber’s CPPTRAJ program [67]. Consequently, the representative pose selected from the dominant cluster was considered as a predicted ligand pose.

2.6. Ligand Interaction Fingerprint

The Protein-Ligand Interaction Fingerprint application in the MOE software [112] was used to outline the interactions between ligands and proteins with a fingerprint scheme. Interactions such as hydrogen bonds, ionic interactions and surface contacts are classified in accordance with the residue of the origin and built into a fingerprint scheme which is representative of a given database of protein–ligand complexes.

3. Results

3.1. Molecular Docking Analysis

We docked the 15 test compounds on the SARS-CoV-2 spike protein (open and closed conformations), CD147 and α7nAChr (activated, desensitized and resting states). All docking scores and binding sites are reported in

Table 2. Results of the docking analysis for spike protein S1 binding sites on NTD and RBD in open and closed positions. Scores listed are maximum absolute values for the sites listed in Table 1 for NTD or RBD, open or closed, with the maximum for all four combinations shown in column 2. Compounds are sorted in descending order of that maximum |Score| (column 2).

Compound Name	Maximum |Score|		Open				Closed
			NTD		RDB		NTD		RBD
	Score (kcal/mol)	At Site	Score (kcal/mol)	Site	Score (kcal/mol)	Site	Score (kcal/mol)	Site	Score (kcal/mol)	Site
Ivermectin	−8.948	NTD-open site 10	−8.948	site 10	−8.256	site 17	−8.205	site 4	−7.735	site 22
Moxidectin	−8.902	NTD-open site 2	−8.902	site 2	−8.218	site 21	−7.659	site 2	−7.989	site 18
Doramectin	−8.885	NTD-open site 2	−8.885	site 2	−8.144	site 21	−8.867	site 9	−8.216	site 19
Oleandrin	−8.787	RBD-closed site 19	−7.787	site 10	−8.051	site 22	−8.083	site 14	−8.787	site 19
Selamectin	−8.774	NTD-closed site 10	−8.476	site 15	−7.432	site 19	−8.774	site 10	−8.142	site 16
Okadaic acid	−8.716	NTD-open site 10	−8.716	site 10	−8.067	site 21	−7.937	site 4	−8.25	site 18
Gitoformate	−8.514	NTD-open site 10	−8.514	site 10	−7.669	site 21	−7.88	site 10	−7.992	site 19
Amphotericin_B	−8.304	NTD-open site 15	−8.304	site 15	−7.516	site 21	−7.931	site 4	−7.332	site 21
P-57AS3	−8.045	NTD-open site 4	−8.045	site 4	−7.663	site 22	−7.704	site 5	−7.627	site 19
Eprinomectin	−7.646	NTD-open site 6	−7.646	site 6	−7.584	site 21	−7.088	site 6	−7.302	site 21
Concanamycin A	−7.564	NTD-open site 10	−7.564	site 10	−7.335	site 19	−7.347	site 3	−7.302	site 21
Natamycin	−7.529	RBD-open site 21	−7.388	site 13	−7.529	site 21	−7.359	site 4	−6.87	site 18
Nystatin	−7.333	RBD-open site 21	−7.226	site 6	−6.845	site 21	−6.867	site 14	−6.773	site 19
beta-Escin	−7.324	NTD-open site 10	−7.324	site 10	−7.333	site 21	−7.264	site 4	−7.296	site 19
Fusicoccin	−6.705	NTD-open site 2	−6.705	site 2	−6.123	site 22	−6.353	site 10	−6.381	site 18

and

Table 3. Results of the docking analysis on CD147 and α7nAChr. Compounds are sorted in descending order according to |Score| separately for CD147 and α7nAChr.

CD147			α7nAChr
Compound Name	Score (kcal/mol)	Site	Compound Name	Score (kcal/mol)	Site
Okadaic acid	−8.578	site 5	Ivermectin	−10.636	Activated site 2
Doramectin	−8.253	site 1	Doramectin	−10.243	Activated site 2
Selamectin	−8.082	site 5	Okadaic acid	−10.240	Activated site 2
P-57AS3	−8.010	site 1	Moxidectin	−10.142	Resting site 1
Concanamycin A	−7.847	site 9	Concanamycin A	−9.932	Activated site 2
Ivermectin	−7.527	site 5	P-57AS3	−9.799	Desensitized site 3
Amphotericin_B	−7.481	site 1	Gitoformate	−9.794	Resting site 1
Moxidectin	−7.469	site 1	beta-Escin	−9.711	Resting site 3
Oleandrin	−7.434	site 4	Natamycin	−9.611	Activated site 1
Gitoformate	−7.297	site 8	Oleandrin	−9.465	Activated site 2
Nystatin	−7.038	site 9	Selamectin	−9.397	Activated site 2
Eprinomectin	−6.827	site 9	Nystatin	−9.214	Resting site 3
beta-Escin	−6.755	site 1	Eprinomectin	−8.968	Resting site 3
Natamycin	−6.739	site 7	Fusicoccin	−8.814	Resting site 3
Fusicoccin	−5.872	site 1	Amphotericin_B	−8.811	Resting site 3

. In the following sections, we discuss the top five spike and α7nAChr and the top six CD147 protein inhibitors, as well as the common inhibitors within the top-ranked inhibitors for all the receptors.

All individual docking scores at all sites for ivermectin on the spike, CD147 and α7nAChr are presented in Table S3, Table S4 and Table S5, respectively.

3.2. Selection of the Most Promising Compounds

The top five inhibitors, in descending order of the absolute value of the Score, were found to be as follows: for the spike protein: ivermectin, moxidectin, doramectin, oleandrin and selamectin; for CD147: okadaic acid, doramectin, selamectin, P-57AS3, concanamycin A and ivermectin; and for α7nAChr: ivermectin, doramectin, okadaic acid, moxidectin and concanamycin A. The common inhibitors within the five top-ranked inhibitors were, for the spike and α7nAChr: ivermectin, doramectin, okadaic acid and moxidectin; for the spike and CD147: doramectin; and for CD147 and α7nAChr: okadaic acid, doramectin and concanamycin A. The only common top inhibitor for all the receptors was doramectin. The majority of compounds bound to the spike NTD, while the highest affinity could be observed towards the open conformation. As for α7nAChr, the activated and resting states were preferred with respect to the desensitized state. Moreover, compounds had higher affinities to the activated state of α7nAChr.

All the top inhibitors considered, with the exception of oleandrin, were found to bind to S1-NTD. Both ivermectin and selamectin bound to site 10 of S1-NTD (Figure 2A,C), which is a glycosylation binding site (N61). Moxidectin and doramectin bound to site 2 S1-NTD (Figure 2A), which is a sialoside binding site proposed by Behloul et al. [113]. Oleandrin bound to site 19 of S1-RBD proposed by Carino et al. [117] (Figure 2B). Site 19 includes the N388 glycosylation binding site.

The binding poses of all the compounds with high affinity for CD147 are shown in Figure 3. Okadaic acid, selamectin and ivermectin were found to bind to site 5, which is located in domain A of CD147 protein. Doramectin and P-57AS3 were found to bind to site 1 of CD147, which is in the interface of domain 1 and domain 2 of CD147. Concanamycin A was found to bind to site 9 of CD147.

The binding poses of all the compounds with high affinity for α7nAChr are shown in Figure 4. Ivermectin, doramectin, okadaic acid and concanamycin A were found to bind to site 1 of the activated conformation of α7nAChr (Figure 4A,C). Moxidectin was found to bind to site 1 of the resting conformation of α7nAChr (Figure 4B,D). In what follows, the interactions of the top inhibitors for the spike, CD147 and α7nAChr will be discussed.

3.3. Molecular Dynamics Simulations and RMSD Analysis

A 100ns-long MD simulation was performed to check the stability of each protein-inhibitor complex and to discriminate between stable and unstable docked poses. The top-docked pose (with the lowest docking Score) for each protein–ligand complex was used as an initial structure for the simulations. The binding stability was assessed by following the time evolution of the ligand RMSD in each trajectory, where we used the starting structure as a reference and RMSD alignment was carried out on protein atoms.

From the RMSD analysis and a visual inspection of MD trajectories, we found that, except for salemectin, all of the top five compounds in complex with the spike protein were relatively stable, reaching a RMSD plateau between 2 Å and 4 Å (Figure S5). Conversely, CD147 went through hinge movements during MD (Figure S6), which made it difficult to align the structures and caused fluctuations and higher RMSD values. Visual inspections and ligand–protein interaction analysis (Section 3.4) confirmed that all compounds, except for okadaic acid, maintained their binding to the same binding site during MD simulations (2 Å < RMSD < 6 Å) (Figure S6). Regarding α7nAChr, a common behaviour was observed for almost all of the compounds: before MD, binding to α7nAChr occurred through the interaction between the disaccharide group of each ligand and the activated site2 of α7nAChr inside the pore (except for moxidectin, which bound to the outer wall of α7nAChr). Benzofuran and spiroketal groups were pointed toward the center of the pore, with no apparent hydrogen bonds with any residue. After conducting MD simulations, the stable structure of compounds tended toward a conformation that maintained its binding with activated site 2, with extra binding through the benzofuran group, by getting close to the pore wall. Ivermectin, okadaic acid and moxidectin manifested a stable RMSD (1.5 < RMSD < 4) (Figure S7). An abrupt shift in the RMSD of doramectin was due to the detachment of the benzofuran group from one monomer and the attachment to another monomer due to the symmetry of the α7nAChr protein. The new conformation still bound, through disaccharide, with the same binding site, and it was as stable as the first conformation. During visual inspection and through ligand–protein interactions, it was confirmed that concanamycin underwent major binding adjustments with regard to its initial docked conformation and ended up leaving the binding site.

3.4. Analysis of the Protein–Ligand Interactions

In stable MD trajectories, the top representative pose of each compound was selected from the populationally dominant cluster using clustering analysis on all the trajectories for further ligand–protein interaction analysis. The Protein-Ligand Interaction Fingerprint module of MOE was used to summarize the interactions between ligands and proteins with a fingerprint scheme. N61, R415, F157 and D40 emerged as main residues of the spike protein due to their interaction with high-affinity compounds (Figure S8). As for CD147, the residues interacting with the selected compounds were L46, K87, R85 and H32 (Figure S9). In case of α7nAChr, four out of the five selected compounds bound to activated site 2 and interacted with P16, N106, W85 and N100, that are exposed on the interior surface of the protein channel. One compound, namely moxidectin, interacted with N110 of resting state α7nAChr, which is exposed on the outer surface of the protein (Figure S10).

The interaction mechanisms of ivermectin with the SARS-CoV-2 spike protein, CD147 and α7nAChr were analysed using MOE software. Binding energies were obtained through the GBVI/WSA forcefield-based scoring function, which uses the AMBER99 forcefield to compute electrostatic, solvation, van der Waals and surface area contributions to the free energy given the ligand pose. Two to four hydrogen bond acceptor interactions were characterized in the best pose of the compounds in all receptors.

Ivermectin remained in the same binding site for all the receptors during MD simulations (2 Å < RMSD < 4 Å). In case of the spike protein (RMSD~2.5 Å), N61 (the main residue of the glycosylation site 10) were involved, with a binding energy of −2.9 kcal/mol, with the benzofuran group of ivermectin and R415 were involved with the lactone group of ivermectin, with a binding energy of −2.9 kcal/mol (Figure 5A).

During MD simulations, CD147 went through hinge movements and gave rise to a relatively higher (RMSD < 6 Å) value for ivermectin. Ivermectin stayed stable after 60 ns and strongly bound to CD147 through its disaccharide group, featuring E43 and K54 residues with −2.7 kcal/mol and −4.1 kcal/mol of binding energies, respectively, and a lactone core group featuring L46 residue with −1.1 kcal/mol of binding energy (Figure 5B).

As for α7nAChr, strong hydrogen bond acceptor interactions were found with K86 (−6.8 kcal/mol of binding energy) and N106 (−3 kcal/mol of binding energy). Moreover, it was characterized by an additional hydrophobic interaction with H85 (−2.7 kcal/mol of binding energy) (Figure 5C). In addition to maintaining disaccharide group binding with α7nAChr through K86 and N106, the equilibrated structure formed an extra binding to α7nAChr through its benzofuran group with H85 compared to the initial docking pose. Ivermectin maintained its attachment to α7nAChr at the same binding site with (RMSD < 4 Å). The presence of the same type and number of interactions in the analyzed proteins may support the hypothesis of a multi-targeted action of ivermectin.

A summary of the amino acid mutations of the SARS-CoV-2 Alpha, Beta, Gamma and Delta variants with a focus on the spike protein is presented in Table S6. We do not expect that the given variant will have a significant effect on the binding of the selected compounds, considering that the mutations are not directly involved in the binding sites of ivermectin and related compounds. However, it is noteworthy to mention that allosteric interactions should be taken in consideration for a comprehensive and accurate evaluation.

3.5. Bioactivity of the Test Agents with Greatest Binding Strength

By Lipinski’s rule of five, agents with a molecular mass greater than 500 would tend to be suboptimally bioactive as oral agents. However, although among these test agents, ivermectin and doramectin, for example, have molecular masses of 875.1 and 899.1, respectively, both are well-absorbed with similar pharmacokinetics [125]. Ivermectin, in particular, is distributed throughout the human body within eight hours of oral administration [83,126,127], and its success in combatting diseases affecting hundreds of millions of people is well established [70].

3.6. Protein-Protein Interactions

The spike (PDB: 6VSB for open conformation) and α7nAChr (PDB: 7KOX) initial structures were obtained from the RCSB Protein Data Bank. PatchDock software (bioinfo3d.cs.tau.ac.il/PatchDock/) was used for protein–protein docking simulations [128,129]. PatchDock is a geometry-based molecular docking algorithm aimed at finding docking transformations that yield good molecular shape complementarity. Each candidate model is further evaluated by a scoring function that considers both the atomic desolvation energy and the geometric fit. The results obtained from PatchDock were further refined with the associated server FireDock, which delivers a further refinement of both the score function and of the complexes’ geometries. We present the highest scoring structure in Figure 6. The best docking pose indicates the interaction between two RBD segments of the spike trimer: from the RBD part of chain B (red) and chain C (green) to the outer surface of two subunits (chain A (cyan) and chain E (gray) of α7nAChr pentamer). We presented the highest scoring spike–α7nAChr complex in Figure S11. The Protein Contacts panel of the MOE software was used to study the interaction between the atoms of proteins. The interaction between the two proteins were evaluated using six types of contacts: Hydrogen bonds (Hbond), metal, ionic, arene, covalent and Van der Waals distance interactions (Distance). We identified the Van der Waals distance interactions between chain E (gray) of α7nAChr and chain C (green) of the spike protein (Figure S11A). There is a main interaction between the receptor-binding motif (aa 437–508) of the spike RBD and aa 186–192 of the extracellular domain of the nAChR 9 subunit. Previously, Farsalinos et al. reported aa 189–192 of the extracellular domain of α7nAChr as part of the a region which forms the core of the “toxin-binding site” of the nAChRs [130]. There is Van der Waals distance interactions between chain A (cyan) of α7nAChr and chain C (green) of the spike protein (Figure S11B). We also identified Hbond interactions between chain E (gray) of α7nAChr and chain B (red) of the spike protein (Figure S11C) and chain A (cyan) α7nAChr and chain B (red) of the spike protein (Figure S11D). Further details relating to the interaction between different chains are presented in Figure S11.

4. Discussion

Protein–ligand docking is a powerful and popular computational tool to simulate drug–target interactions. Several in silico studies [66,67,68,69,74,75,76,77,78,79] have explored whether competitive binding at subdomains of interest on the SARS-CoV-2 spike protein by ivermectin could explain its efficacy against COVID-19, as indicated in the several RCTs and animal studies related above. In one of these molecular docking studies, Lehrer and Rheinstein (2020) examined potential sites on the SARS-CoV-2 S1 RBD at which ivermectin might bind and competitively block attachment to ACE2, limiting viral replication [74]. They identified one such site at which ivermectin was predicted to dock with high binding energy.

The potential for competitive binding by ivermectin on the spike protein NTD, the subdomain with the highest concentration of glycan binding sites, however, is of interest, especially given the importance of the glycan bindings of SARS-CoV-2 for initial attachments to host cells and the possibilities for hemagglutination, as described above. In particular, the nanometer-scale spacing and the composition of terminal sugar molecules (including SA, galactose, mannose, fucose, N-acetylglucosamine (GlcNAc) and/or N-acetylgalactosamine (GaMAc) for the 22 N-glycosylation sites of the SARS-CoV-2 spike protein [131]) meshes closely with the spacing and terminal sugar composition of glycophorin A [132,133], a ubiquitous molecule on the RBC surface that has no known physiological purpose other than the clearance of viruses and other pathogens [30,31,134].

Here, the AutoDock Vina program was used to perform binding affinity computations for the 15 test compounds (ivermectin and 14 related molecules) for seven binding sites on RBD and 15 on NTD, as identified in the literature as potential SARS-CoV-2 spike protein binding sites of interest for druggability. The inclusion of NTD as well as RBD sites allowed for the consideration of potential competitive binding by ivermectin and related molecules to limit initial viral attachments to host cells and potential hemagglutination-related morbidities. We also examined potential bindings of ivermectin to CD147, an SA-tipped receptor that is densely distributed on RBCs, to provide some indication as to whether ivermectin might limit glycan bindings of SARS-CoV-2 at the host cell end as well. The potential for binding by ivermectin to α7nAChr to inhibit viral attachment to that receptor and activate the cholinergic anti-inflammatory pathway was also explored.

We calculated binding affinities using the Vina Score value, selecting for those bindings most likely to be physically realized. Several of the 15 test compounds, including ivermectin, had bindings of strong or moderate affinity to sites on the spike protein, CD147 and α7nAChr, as detailed here, but since ivermectin, a safe and widely available drug, has been the subject of closest study for COVID-19 treatment among these compounds, the discussion below focuses on the results for that agent and their significance.

As reported in Table S3, docking computations for ivermectin binding to the spike protein found the strongest binding (−8.948 kcal/mol) at site 10 of S1-NTD, which is a glycosylation binding site (N61), in the open position. A study of AutoDock binding energies calculated for a large set of HIV inhibitors and likely non-inhibitors against multiple ligands found that the selection of binding energy < −7.0 kcal/mol identified the inhibitors with 98% sensitivity and 95% specificity [135]. It is thus noteworthy that for the sites at the NTD and RBD in the open position and for NTD in the closed position, most of the binding energies were <−7.0 kcal/mol, and so, per the above, this indicates their capability to be physiologically active. Likewise, binding affinities < −7.0 kcal/mol for ivermectin at five of the 12 sites of CD147 and of 30 of the 37 sites of α7nAChr (for the desensitized, activated and resting states, total) indicate a capability for physiologically-manifested binding to these host receptors as well. Despite the above-cited indication of high sensitivity and specificity for physiological efficacy with binding energies < −7.0 kcal/mol, it is clear that a physiological relevance corresponding to this study’s results can only be clearly established through follow-up confirmation with in vitro and/or in vivo findings.

When considering the binding affinities of ivermectin to NTD sites, it is significant that glycan bindings from SARS-CoV-2 and other coronaviruses to host cells are generally weak when univalent but orders of magnitude stronger when multivalent [4,12,29]. Thus, ivermectin, the molecular dimensions of which span approximately 2 × 1 nm [136] (with the length of the spike protein being ~20 nm [137,138]), could block clinically relevant multivalent bindings from the spike protein to host cells by steric interference, even if its actual bindings to some glycan sites on the spike protein were somewhat weaker than predicted. It is relevant, here, that in a rough order of magnitude, considering an average of 0.41 spike protein punctae found attached per RBC in COVID-19 patients [20] and a peak serum concentration of 137.4 nM for ivermectin plus active metabolites after the ingestion of ivermectin with a fatty meal at a dose of 200–350 µg/kg, which is in the range of standard dosing, there would be about 126,000 molecules of ivermectin and active metabolites per spike protein molecule in blood [4].

Two especially significant consequences of these predicted multiple bindings of ivermectin along the spike protein with binding affinities mostly less than −7.0 kcal/mol are, first, that these could provide effective competitive binding for all variants of the virus and, second, since multivalent bindings govern spike protein attachments to host cells, it is noteworthy that competitive inhibitors of such multiple bindings having only moderate binding affinities at individual viral binding sites can have strong inhibitory effects on total attachment strength [139]. Binding affinities as computed < 7.0 kcal/mol at five of 12 sites of CD147 further indicate the potential for ivermectin to limit glycan bindings to meshing glycan binding sites on host cells, and also to limit inflammatory pathways mediated by that receptor.

Previous studies using both chick and human cells have demonstrated that a micromolar concentration of ivermectin strongly potentiates the ACh-evoked current of the α7nAChr receptor [89], as expressed on neuronal cells as well as on different airway cells, such as bronchial epithelial cells and type II alveolar epithelial cells, endothelial cells and cytokine secreting immune cells (i.e., macrophages, lymphocytes and mast cells) [140,141]. Importantly, α7nAChr is one of the main receptors under the control of the vagus nerve used by the parasympathetic nervous system to regulate multiple physiological homeostatic mechanisms commonly affected during SARS-CoV-2 infection, including the respiratory rate, heartbeat, blood pressure, vessel tone, hormone secretion, intestinal peristalsis, digestion and inflammation [142]. Our computational studies with ivermectin are consistent with in vitro experimental results obtained with chick and human cells [89] and confirm high affinity bindings to α7nAChr (i.e., Score of −10.636 kcal/mol, the highest of all the 15 test compounds, and Score < −7.0 kcal/mol at 30 of 37 sites for all three states total). Moreover, in agreement with previous reports [5,6], we were also able to demonstrate a potential direct binding of the SARS-CoV-2 spike-1 protein to α7nAChr, suggesting that this ubiquitous cholinergic receptor may represent an additional port of entry for SARS-CoV-2 into human cells. Taken all together, our computational results demonstrating the high-affinity binding of ivermectin to α7nAChr and a potential direct interaction of the cholinergic nicotinic receptor with SARS-CoV-2 spike 1 on neurons, cytokine secreting cells and endothelial cells, which might potentially explain multiple aspects of SARS-CoV-2 infection pathophysiology including but not limited to (a) the typical loss of smell and taste [94,95], (b) the triggering of the life threatening cytokine storm through inactivation of the cholinergic anti-inflammatory α7nAChr pathway on TNF/IL6/IL1 secreting macrophages [90,91,92] and (c) the impairment of the endothelium dependent acetylcholine-induced vasodilation caused by SARS-CoV-2 spike 1 infection of the lung vasculature [143]. Ivermectin high-affinity binding may therefore potentially shield from infection α7nAChr-expressing host cells while at the same time, through its allosteric agonistic function, potentiate the activation of the cholinergic pathway and attenuate SARS-CoV-2-induced parasympathetic dysregulation by restoring the function of these receptors.

5. Conclusions

In the present study, a computational investigation including molecular docking was conducted to explore the potential bindings of ivermectin and 14 similar compounds to three targets of interest (the spike, CD147 and α7nAChr) that are relevant for drug activity against COVID-19. Strong or moderate affinity bindings were found for ivermectin to multiple sites on the spike protein, CD147 and α7nAChr, which could provide effective competitive bindings to all variants of the virus. According to our calculations, ivermectin binds strongly to a glycosylation binding site (site 10: N61) of the spike protein S1-NTD in the open position and to several other sites on S1 NTD and RBD. We also examined the potential bindings of ivermectin to CD147. Ivermectin was found to bind to site 5, which is located in domain A of the CD147 protein, and to other sites on CD147, indicating that ivermectin might limit glycan bindings of SARS-CoV-2 at the host cell end as well.

Among all the targets, ivermectin has the highest affinity to the α7nAChr receptor. Protein–protein docking results reveal a potential direct binding of the SARS-CoV-2 spike-1 protein to α7nAChr, suggesting that this ubiquitous cholinergic receptor may mediate SARS-CoV-2 entry into cells, shedding light on multiple aspects of SARS-CoV-2 infection pathophysiology (i.e., the loss of smell and taste, the cytokine storm and impairment of the endothelium-dependent acetylcholine-induced vasodilation). In this context, the high affinity of ivermectin and related compounds to α7nAChr may both prevent viral entry and potentiate the activation of the cholinergic pathway and attenuate SARS-CoV-2-induced parasympathetic dysregulation by restoring the function of these receptors. Our preliminary results warrant further in vitro and in vivo testing of the 15 test compounds, in particular ivermectin, an available and safe drug, against SARS-CoV-2, with the hope of containing the virus and limiting its morbidity.