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

: 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 afﬁnity 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 identiﬁed 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-inﬂammatory pathway controlled by the vagus nerve. Binding afﬁnities were calculated for these multiple docking sites and binding modes of each compound. Our results indicate the high afﬁnity of ivermectin, and even higher afﬁnities 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-inﬂammatory pathway that could limit cytokine production by immune cells.


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

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 antiinflammatory 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].

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 1273 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].
Computation 2022, 10, x 4 of 25 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].

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].

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.

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.

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. 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]. 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).  [11]. Fantini et al. [101], meanwhile, proposed a new type of gangliosidebinding 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-parisdiderot.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].

Binding
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.

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 ligandprotein 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.

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.

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.

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 Tables 2 and 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.

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 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 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.

Molecular Dynamics Simulations and RMSD Analysis
A 100ns-long MD simulation was performed to check the stability of each proteininhibitor complex and to discriminate between stable and unstable docked poses. The topdocked 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 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. 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 ligandprotein interactions, it was confirmed that concanamycin underwent major binding adjustments with regard to its initial docked conformation and ended up leaving the binding site.

Molecular Dynamics Simulations and RMSD Analysis
A 100ns-long MD simulation was performed to check the stability of each proteininhibitor 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 site 2 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.

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). 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.

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].

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 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.

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].

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 (bioinfo3 d.cs.tau.ac.il/PatchDock/ accessed on 2 March 2022) 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. [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. Figure 6. Spike-α7nAChr complex model. Spike protein trimer is colored in dark blue (chain A), red (chain B) and green (chain C). α7nAChr pentamer is colored in cyan (chain A), pink (chain B), yellow (chain C), brown (chain D) and gray (chain E). The yellow and green parts of α7nAChr are interacting with the dark blue and gray monomers from spike protein. Chain B (red) and chain C (green) of α7nAChr are interacting with the chain A(cyan) and chain E(gray) of spike protein.

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 Figure 6. Spike-α7nAChr complex model. Spike protein trimer is colored in dark blue (chain A), red (chain B) and green (chain C). α7nAChr pentamer is colored in cyan (chain A), pink (chain B), yellow (chain C), brown (chain D) and gray (chain E). The yellow and green parts of α7nAChr are interacting with the dark blue and gray monomers from spike protein. Chain B (red) and chain C (green) of α7nAChr are interacting with the chain A(cyan) and chain E(gray) of spike protein.

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 hemagglutinationrelated 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.

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-2induced 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.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/computation10040051/s1, Figure S1: Binding sites of CD147 protein obtained from MOE (Site Finder); Figure S2: Binding sites of (A) desensitized, (B) activated and (C) resting conformations of α7nAChr protein obtained from MOE (Site Finder); Figure S3: Positive control docking of epibatidine of (PDB:7KOX) of the α7nAChr receptor; Figure S4: Decoy compounds; Figure S5. Time-evolution of the RMSD of top-ranked inhibitors with respect to spike; Figure S6. Time-evolution of the RMSD of top-ranked inhibitors with respect to the CD147 receptor; Figure S7. Time-evolution of the RMSD of top-ranked inhibitors with respect to the α7nAChr receptor; Figure S8. Ligand interaction plots of compounds selected for spike inhibition; Figure S9. Ligand interaction plots of compounds selected for CD147 inhibition; Figure S10. Ligand interaction plots of compounds selected for α7nAChr inhibition; Figure S11: Protein-protein interaction between (A) chain E (gray) α7nAChr and chain C (green) of spike protein, (B) chain A (cyan) α7nAChr and chain C (green) of spike protein, (C) chain E (gray) α7nAChr and chain B (red) of spike protein and (D) chain A (cyan) α7nAChr and chain B (red) of spike protein; Table S1: Binding sites of CD147 obtained from MOE (Site Finder); Table S2: Binding sites of α7nAChr obtained from MOE (Site Finder); Table  S3: Results of the docking analysis of ivermectin for spike protein S1 on all binding sites on NTD and RBD in open and closed positions; Table S4: Results of the docking analysis of ivermectin on all sites of CD147; Table S5: Results of the docking analysis of ivermectin on all sites of α7nAChr; Table S6: Amino acid mutations of SARS-CoV-2 Alpha, Beta, Gamma and Delta variants, with a focus on spike protein.

Data Availability Statement:
The datasets used and analyzed in the current study are available from the corresponding author upon request.

Conflicts of Interest:
A.D.S. reports grants from PUMA, grants from IMMUNOMEDICS, grants from GILEAD, grants from SYNTHON, grants and personal fees from MERCK, grants from BOEHINGER-INGELHEIM, grants from GENENTECH, grants and personal fees from TESARO and grants and personal fees from EISAI. The other authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: α7nAChr alpha-7 nicotinic acetylcholine receptor ACE2 angiotensin