Comparative study of the interaction of ivermectin with proteins of interest associated with SARS-CoV-2: A computational and biophysical approach
1. Introduction
SARS-CoV-2 is a novel virus belonging to the β-Coronavirus genus of the 2B group of the Coronaviridae family. This interesting virus contains only 29 proteins, 26 of which have been successfully expressed for in vitro studies to determine targets of interest for drug discovery [1]. For example, the conserved cysteine protease Mopar has been highlighted as an exciting target as it mediates the maturation cleavage of polyproteins during virus replication [2,3].
Despite great successes in the production and roll-out of vaccines against SARS-CoV-2, new variants are on the rise and there is still no globally accepted treatment for COVID-19 (https://www.who.int/publications/i/item/WHO-2019-nCoV-clinical-2021-1). During the last year, many laboratories focused on screening FDA-approved drugs for quick implementation in clinical settings [[4], [5], [6]]. A compound of interest from such work is the racemic mixture ivermectin, typically used to treat helminth infections. Most of the studies on the macrocyclic lactone ivermectin only consider the major constituent B1a in their dockings [7]. However, ivermectin is an approximately 80:20 mixture of two homolog derivatives of the compound avermectin B1, called 22,23-dihydro-avermectin B1a (HB1a) and B1b (HB1b) correspondingly, which differ in the presence of a sec-butyl and isopropyl group, at the C25 position, respectively (Fig. 1). Interestingly, this mixture has demonstrated in vitro antiviral activity against several single-stranded RNA viruses, such as Zika virus, dengue virus, Chikungunya virus, avian influenza A virus, Porcine Reproductive and Respiratory Syndrome virus, human immunodeficiency virus type 1, among others, including SARS-CoV-2 [8].

2. Materials and methods
2.1. Databases and structure selection
As the nuclear import for macromolecules is facilitated by importins, the structures of importin α1 subunit (PDB: 5KLR) from Mus musculus and importin β1 subunit (PDB: 2P8Q) from Homo sapiens were used as a model for the members of the nuclear import superfamily. The host nuclear import system can be bound and sequestered by pathogens such as SARS-CoV-2 allowing the transportation of viral proteins to the host nucleus leading to increased viral replication [[12], [13], [14], [15]]. Additionally, we also consider the multi-functional helicase (nsp13) of SARS-CoV-2 responsible for viral replication (PDB: 6ZSL) [[16], [17], [18]], and the main protease (Mopar) of SARS-CoV-2 (PDB: 6LU7) as it is a key enzyme of coronaviruses and has a fundamental role in mediating viral replication and transcription, making it an attractive target for drugs [11,[19], [20], [21], [22]]. All structures were obtained in PDB format from the RCSB protein database (https://www.rcsb.org/). The homologs structures of HB1a (CID_6321424) and HB1b (CID_6321425) that make up ivermectin was obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) in SDF format, and the SMILES online converter was used (https://cactus.nci.nih.gov/translate/) to get a PDB format. To avoid confusion throughout the study, we only refer to ivermectin when mentioning the mixture of homologs (a mixture of approximately 80:20 of the two homologs derivatives of the compound avermectin B1, called 22,23-dihydro-avermectin B1a and B1b), while for the individual study of each homolog, the type of avermectin analyzed (HB1a or HB1b) was always indicated. For a comparison between the ADME profiles of the avermectin homologs, the PubChem databases (https://pubchem.ncbi.nlm.nih.gov/) were used; SwissADME (www.swissadme.ch/) and the Molinspiration Property Engine v2018.10
2.2. Comparative molecular docking
We performed a comparative analysis using five popular molecular docking models for a rigorous prediction of the standard free energy (ΔG) of binding of ligand-protein complexes. The complexes were built in the programs MTiAutoDock (https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::MTiAutoDock), webinar (https://durrantlab.pitt.edu/webina/), DINC 2.0 (http://dinc.kavrakilab.org/), COACH-D (https://yanglab.nankai.edu.cn/COACH-D/) and doctor (https://dockthor.lncc.br/v2/) were used. They all represent some of the basic, improved, and more advanced versions of molecular docking associated with the efficient AutoDock algorithm either in the sampling stage, blind docking, or in the scoring function. To increase accuracy, a minimum of 10 runs per program was performed which implied approximately 106 evaluations per run. The rest of the parameters were considered by default. Additionally, to validate the docking results, the Pose&Rank server (https://modbase.compbio.ucsf.edu/poseandrank/) was used to score the protein-ligand complexes, using the statistical scoring function dependent on the atomic distance RankScore. As usual, all the water molecules were removed and the PDB files were separated into two different files, one containing the protein and the other containing the ligand structure. Only the three runs with the most favorable berth were considered in the sampling of the probabilistically most feasible and thermodynamically most favorable positions in the complexes. This criterion was used to discriminate the complexes that would be subjected to further analysis, including potential theoretical inhibition and molecular dynamics.
2.3. Comparative analysis of the hydrophobic characteristics of the binding sites
We use PockDrug-Server that predicts the drug delivery capacity in the pockets considering the Kyte-Doolittle Pocket Hydrophobicity Scale [23]. The hydrophobic characteristics of the binding sites were also analyzed using Biovia Discovery Studio 2021 [24].
