![Free access]({"type":"image" , "src":"/pb-assets/3rdPartyLogos/FREE_ACCESS_ICON-1465051266437.png"}){kind=icon}
![ORCID Icon]({"type":"image","src":"/templates/jsp/images/orcid.png"} "ORCID Icon")
https://orcid.org/0000-0002-8530-8711![CrossMark Logo]({"type":"image" , "src":"/templates/jsp/images/CROSSMARK_Color_horizontal.svg"})
Abstract
?
![MathJax Logo]({"type":"image","src":"/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif"})
?Abstract
The novel corona virus (Covid-19) has become a great challenge worldwide since 2019, as no drug has been reported yet. Different clinical trials are still under way. Among them is Ivermectin (IVM), an FDA approved drug which was recently reported as a successful candidate to reduce SARS-CoV-2 viral load by inhibiting Importin-α1 (IMP-α1) protein which subsequently affects nuclear transport of viral proteins but its basic binding mode and inhibitory mechanism is unknown. Therefore, we aimed to explore the inhibitory mechanism and binding mode of IVM with IMP-α1 via different computational methods. Initially, comparative docking of IVM was performed against two different binding sites (Nuclear Localization Signal (NLS) major and minor sites) of IMP-α1 to predict the probable binding mode of IVM. Then, classical MD simulation was performed (IVM/NLS-Major site and IVM/NLS-Minor site), to predict its comparative stability dynamics and probable inhibitory mechanism. The stability dynamics and biophysical analysis of both sites highlighted the stable binding of IVM within NLS-Minor site by establishing and maintaining more hydrophobic contacts with crucial residues, required for IMP-α1 inhibition which were not observed in NLS-major site. Altogether, these results recommended the worth of IVM as a possible drug to limit the SARS-CoV-2 viral load and consequently reduces its progression.
Communicated by Ramaswamy H. Sarma
1. Introduction
In early December 2019, an outbreak of novel Coronavirus Disease-2019 (COVID-19) appeared in Wuhan, China(Huang et al., 2020; Zhu et al., 2020) that eventually becomes global ongoing pandemic owing to its contagious nature (Chen, Hu, et al., 2020). The Chinese center for disease control and prevention (CCDC) identified causative agent on January 7, 2020 from throat swab samples and named it Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). SARS-CoV-2 is a single-stranded RNA virus. The disease was named COVID-19 by the World Health Organization (WHO) (Sohrabi et al., 2020). The patients had developed symptoms of fever, dry cough, and sore throat (Chen, Mao, et al., 2020). Most of the cases were spontaneously recovered however; few had developed several fatal complications like organ failure, pulmonary edema, severe pneumonia, septic shock, and Acute Respiratory Distress Syndrome (ARDS). The SARS-CoV-2 outbreak was declared as a Public Health Emergency of International Concern (PHEIC) and a pandemic on January 30, 2020 and then on March 11, 2020 by the World Health Organization (WHO) (Vardhan). The total confirmed COVID-19 cases by 15 April 2020 reported to 1,914,916 across 210 countries with a total 123,0107 deaths (2020).
Various studies have been conducted and some are still in process to develop COVID-19 vaccine. According to recent research conducted by WHO 21 vaccine candidates against COVID-19 are evaluated clinically (www.who.int) (Calina et al., 2020). So far remdesivir is the only drug approved for conditional marketing in European Union. Although no antiviral drug has been clinically approved for SARS-CoV-2 (Chen, Mao, et al., 2020), however clinical trials on various vaccine candidates are under way (Chen, Hu, et al., 2020; Gao et al., 2020; Gautret et al., 2020). In this regard, the antiviral activity of an Food and Drug Administration (FDA) authority approved antiparasitic drug Ivermectin (IVM) (Buonfrate et al., 2019; Canga et al., 2008), was studied for novel causative virus (SARS-COV-2) (Caly et al., 2020). In vitro studies of IVM had revealed antiviral activity against various viruses such as dengue virus, influenza virus, Human Immunodeficiency Virus (HIV-1) (Götz et al., 2016; Lundberg et al., 2013; Tay et al., 2013; Wagstaff et al., 2012). The investigation of its inhibition mechanism revealed that IVM inhibit the association of importin (IMP) α/β 1 heterodimer which is required for the transport of different viral proteins into the nucleus. Consequently, IVM serves to impede nuclear import and viral replication(Wagstaff et al., 2012). IVM has also expressed several other activities, nevertheless it is originally known as inhibitor of nuclear import of host (Kosyna et al., 2015; van der Watt et al., 2016) and viral proteins, for instance simian virus SV40 large tumor antigen (T-ag) (Tay et al., 2013), and the nonstructural protein 5 (NS5) of Dengue virus (Wagstaff et al., 2011; 2012) . It has markedly reduce the infection of RNA Viruses like influenza (Chen, Hu, et al., 2020; Götz et al., 2016).Venezuelan equine encephalitis virus (VEEV), West Nile Virus (Yang et al., 2020), and DENV 1–4 (Tay et al., 2013). This effective broad range activity of IVM for many distinct RNA viruses is due to the fact that these all viruses depend upon IMP α/β 1 for infection and replication (Caly et al., 2012; Jans et al., 2019). Additionally, various in vitro and in vivo studies has also reported the potency of IVM drug against the DNA viruses i.e. pseudo rabies virus (PRV) (Lv et al., 2018). Likewise, in year 2014–2017 clinical trials of IVM against DENV were also carried out, which revealed dengue inhibition and safety both, however no clinical benefits and alteration in viremia was detected.
Considering all above information, Caly et al. (2020) investigated the therapeutic potential of IVM against SARS-COV-2. IVM showed a great inhibitory potential against IMPα/β protein and subsequently affect its nuclear transport. The IMPα/β is responsible for the transport of viral protein into and out of the nucleus and this movement is considered critical for several cellular processes such as differentiation and development which is also essential for disease states including oncogenesis and viral diseases. As mentioned earlier that the antiviral activity of IVM against both DENV and HIV-1 is strongly dependent upon importin α/β nuclear import of NS5 polymerase and HIV-1 integrase protein respectively. In case of COVID-19, causative agent SARS-CoV-2 a single stranded positive sense RNA virus displayed great resemblance to severe acute respiratory syndrome coronavirus (SAR-CoV). Investigation of SARS-CoV proteins indicated a vital role of IMPα/β1 during disease in signal-dependent nucleocytoplasmic shutting of the SARS-CoV nucleocapsid protein (Rowland et al., 2005; Timani et al., 2005; Wulan et al., 2015), that may affect host cell division (Hiscox et al., 2001; Wurm et al., 2001). Likewise, the accessory protein Open Reading Frames (ORF6) of SAR-CoV has shown to antagonize the antiviral activity of the STAT1 transcription factor by sequestering IMPα/β1 on the rough ER/Golgi membrane (Frieman et al., 2007; Vardhan). Above all, Caly et al. (2020) had also hypothesized that IVM dependent inhibition of nuclear transport could be a reasonable explanation of its inhibitory potential, but the basic inhibitory and binding mechanism of IVM with reference to IMP α/β1 is still unknown.
Therefore, in this study, we aimed to explore the inhibitory mechanism and binding mode of IVM in association with importin-α1 protein by utilizing various computational approaches. In the first phase, comparative docking of IVM was performed against two different binding sites (Nuclear Localization Signals (NLS) major and minor sites) of importin-α1 to predict the probable binding mode of IVM. While in the second phase, we performed classical molecular dynamics simulation of both complexes, i.e. (I) IVM/NLS-Major site and (II) IVM/NLS-Minor site.
2. Materials and methods
2.1. Protein preparation
The crystal structure of Importin-α1 in complex with a phosphomimetic peptide GM130 (PDB ID: 6K06) was accessed from RCSB Protein Data Bank with the resolution of 1.75 Å (Chang et al., 2019). As we need to explore the inhibitory mechanism of IVM within the binding cavity of target protein, the peptide and water molecules were removed and apo protein was selected and prepared. Furthermore, the protein structure was refined and optimized to relieve all the steric clashes via preparation wizard in MOE2019.01 (MOE, 2021). Afterwards, protein was charged and minimized by using AMBER99 forcefield to attain the lowest possible energy.
2.2. Molecule preparation
The 2-D (.smi format) Ivermectin (IVM) was extracted from PUBCHEM (PubChem CID 6321424) (Pubchem.) and converted to 3-D format by MOE Builder Suite of MOE2019.01 (MOE, 2021). The molecule was prepared, charged and minimized by the MM94x forcefield.
2.3. Molecular docking protocol
Since the binding mode of Ivermectin after binding with IMP-α1 has not been investigated yet. Therefore, first of all IVM molecule was docked into both NLS-major and minor site to investigate the comparative binding mode of Ivermectin within two different binding sites of IMP-α1 i.e. NLS-Major and NLS-Minor site. In this regard, the Dock module of MOE2019.01 (MOE, 2021) was utilized to analyze the suitable binding mode of Ivermectin. Before execution of docking, we benchmark different combination of scoring and placement methods of MOE Dock to find out the most suitable algorithm with respect to our target protein. In our case, induce fit docking protocol along with Triangle Matcher algorithm and LondonDG as initial scoring and re-scoring method gave the most reliable results with lowest binding energy, as compared to other combinations of algorithms. Afterwards, 30 conformations of IVM were generated and saved in .mdb file for each site. Finally, the best ligand conformation with highest scoring and lowest binding energy were visually inspected at molecular level to predict the probable binding mode of IVM within one cavity over other.
2.4. Molecular dynamic (MD) simulation protocol
To investigate the comparative inhibitory mechanism and stability dynamics of Ivermectin within two different binding pockets of importin-α1, Molecular Dynamic Simulation of three systems were carried out including (i) Apo-system having importin-α1 protein, (ii) IVM/NLS-Major site of IMP-α1, and (iii) IVM/NLS-Minor site of IMP-α1. The topology of Ivermectin for all the system was generated using Automated Topology Builder (ATB) web server (Malde et al., 2011). While the topology of target protein was generated by using pdb2gmx module in explicit solvent under periodic boundary condition via application of an GROMOS56a force field (Berendsen et al., 1995). Then the systems were solvated in a cubic box of SPCE water molecules with a distance cutoff of 1.0 nm between the solute and edge of water box. The perturbed and unperturbed charges of the system were neutralized by adding appropriate number of counter ions.
Energy minimization was performed at 10 KJ/mol with 50,000 steps of steepest descent algorithm by using Verlet cut off scheme (Páll & Hess, 2013) for each system to remove all the steric clashes and attain lowest possible energy. Then 1 ns equilibration was carried out under constant number of atoms, volume and temperature (NVT) at 300 K. The temperature was regulated by the velocity rescale algorithm. The LINCS holonomic constraints (Hess et al., 1997) were applied with the time step of 2 fs. The particle mesh Ewald (PME) method were used to treat the long-range electrostatic and van der Waals interactions (Hess et al., 1997). After that, during the second step, 1 ns equilibration following the, Parrinello–Rahman semi-isotropic pressure coupling (Parrinello & Rahman, 1981) was performed at constant pressure (1 atm) and temperature (NPT) at and 300 K. Finally, 80 ns long production MD was performed for each system to generate various sample snapshots for analyzing the possible inhibitory mechanism of IVM with respect to its conformational changes. All calculations were performed by GROMACS 5.1.2 (Abraham et al., 2015). The output trajectories of all systems were statistically analyzed to calculate various stability parameters involving root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (RoG) and hydrogen bonds.
2.5. Binding-free energy calculation
Finally, total binding free energy of each simulated system were calculated by the Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) method. According to the stability of each system, last 20 ns of the MD trajectory were used to calculate △Gbind by using the following equation: 
where GComplex illustrates the total energy of protein–ligand complex while Gprotein and Gligand highlights the individual total energy of protein and ligand in the presence of solvent.
3. Results and discussion
3.1. Molecular docking
Molecular docking simulation guided us to predict the binding mechanism of a given molecule with its suitable conformation. The r score is the representative of the binding affinity of molecule as a quantitative measure of quality. Considering this information, we docked IVM into the two different binding cavities of IMP-α1 to find the most suitable binding mode for this drug. Before docking, benchmarking was performed using several combinations of docking algorithms found in MOE (MOE, 2021). Among all the combination, the Triangle Matcher along with LondonDG scoring function was found to predict the highest binding affinity i.e. −13.371Kcal/mol for NLS-Major site, while −11.667Kcal/mol for the second NLS-Minor site, efficiently. The results of all the possible combinations of placement and scoring functions are summarized in (Table 1).
Table 1. Binding affinities of IVM within both NLS major and minor binding site by using different combination of scoring functions and placement methods in docking.
As the typical importin-α1 formed by two ends, N-terminal with importin-β binding domain (Ibb) and a C-terminal with NLS-binding sites, These NLS-binding sites were made up of ten repeating ARMs, from ARM 2–4 were reported for major site, while ARM 6–8 were associated with minor site (Holvey et al., 2015). Previous studies confirmed the role of major site in transportation of various viral proteins or peptides. The NLS-Minor has been used in past to develop small molecular inhibitors (antiviral drugs) that specifically bind to this site with sufficiently high affinity and compete with other nuclear proteins. Apart from it, the minor site was also reported in pathogenesis of malignancies. Both the sites have helical surfaces with hydrophobic interactions except the depth of both the binding cavities which is deeper for minor site as compared to major. As a result of molecular docking, the conformation with highest binding score of ivermectin (–13.37 Kcal/mol) in complex with IMP-α1 were chosen for further visual analysis at molecular level. Ivermectin in the NLS-Major site showed a number of electrostatic and hydrophobic contacts with crucial residues. It was observed that the amine group of Asn188 developed a strong hydrogen bond with oxygen atom of IVM at a distance of 2.83 Å. While the amine group of Arg238 established two hydrogen bonds with two oxygen atoms of IVM molecule at a distance of 2.94 and 2.92 Å, respectively. These significant residues (Ala148, Asp192 and Arg238) (Nakada et al., 2015) were reported as crucial amino acids that are responsible for the inhibition of IMP-α1 protein that eventually reduces the transport of viral protein within the nucleus. Apart from it, IVM also developed some hydrophobic contacts with Gly150, Asp192, Trp231, Lys240, Trp273, Asp280, Arg285and Arg315 of NLS-Major site which helped this molecule to enhance its stability within the binding cavity of target protein, attributable to the lowest binding energy score of IVM in complex with NLS-Major site as illustrated in Figure 1.
Published online:
05 May 2021Figure 1. Three-dimensional interaction diagram of IVM (green) within the binding cavity of NLS-Major site (pink) presenting three hydrogen bonds (red dotted lines) and multiple hydrophobic interactions.
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0001_c.jpg"})
Figure 1. Three-dimensional interaction diagram of IVM (green) within the binding cavity of NLS-Major site (pink) presenting three hydrogen bonds (red dotted lines) and multiple hydrophobic interactions.
Comparative to major site, NLS-Minor site (binding site II) showed low binding affinity (–11.66 Kcal/mol) toward IVM as shown in Figures 2. The detailed visualization of IVM within the binding pocket of NLS-Minor site highlighted less penetration of molecule deep into the cavity that might be because of the large size of IVM that is responsible for its less binding score within minor pocket of IMP-α1. Moreover, the visual inspection at molecular level showed very few hydrogen bonds and more hydrophobic contacts with crucial residue Asn361, Trp357, Glu396 and Trp399 (Nakada et al., 2015) involving single hydrogen bond between the amine group of Asn403 and oxygen atom of IVM at a distance of 2.89 Å along with a single hydrophobic contact between Ala364 and IVM. However, some more hydrophobic interactions were also observed between the aromatic rings of significant residue such as Trp399, Asp325, Arg366, and Glu396 and IVM that might help to stabilize the binding of IVM within this site.
Published online:
05 May 2021Figure 2. 3-D interaction pattern of IVM (green) within the NLS-Minor site (yellow) of IMP-α1 protein showing the number of hydrogen bonds (red dotted lines) and hydrophobic contacts.
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0002_c.jpg"})
Figure 2. 3-D interaction pattern of IVM (green) within the NLS-Minor site (yellow) of IMP-α1 protein showing the number of hydrogen bonds (red dotted lines) and hydrophobic contacts.
3.2. Molecular dynamic simulation
Since docking displayed a static picture of molecular mechanism, therefore we moved towards MD simulation to predict the comparative stability dynamics of IVM with both pockets with respect to time. In this study we calculated various stability parameters with reference to time that is usually difficult to attain by other experimental techniques. After exploring the binding mode of IVM with NLS-Major site via molecular docking, we moved toward molecular dynamics simulation to further understand the stability dynamics and conformational changes induced by IVM in complex with both binding sites. As a result, these analyses will help us in predicting the possible molecular mechanism adopted by IVM to inhibit virus replication and transmission including SARS-CoV-2. Unexpectedly, simulation brought forth surprising results as compared to docking studies, by displaying IVM/NLS-Minor site as most stable system than NLS-Major site. To investigate the key factors responsible for this stability shift of NLS-Minor site over NLS-Major site, the conformational and stability parameters of all three systems i.e. (i) IMP-α1 (apo), (ii) IVM/NLS-Major and (iii) IVM/NLS-Minor of IMP-α1 were analyzed by measuring their RMSD, RoG, RMSF and hydrogen bond calculations.
3.2.1. Root mean square deviation
The values of the positional differences against the backbone of apo-protein and protein–ligand complexes were calculated to investigate the comparative stability and deviation of NLS-Major and Minor sites after binding with IVM. Unexpectedly, RMSD plots (Figure 3) displayed that the system involving IVM bound with IMP-α1 NLS-Major site presented more deviations from its initial structure around 4.0 ± 0.2 Å as compared to the NLS-Minor site. Although this system showed a drastic increase in deviation during the initial simulation but after 50 ns the protein seemed to be stabilized with minimum deviations especially in the last 20 ns of simulation around 4.6 ± 0.2 Å. Comparatively, system having IVM within NLS-Minor site showed least deviations from its initial structure as compared to other two systems and remained stable around 3.8 ± 1.0 Å through-out 80 ns simulation highlighting more stability of IVM with this site. Apart from it, the apo structure showed the highest fluctuations because of the absence of IVM around 6.5 ± 0.5 Å. All these results highlighted that the association of IVM with NLS-Minor site induced more stability in IMP-α1 protein as compared to the major site.
Published online:
05 May 2021Figure 3. The comparative RMSD plots of IMP-α1 protein in apo (black) and complex form with IVM molecule within its both NLS-Major (red) and NLS-Minor site (green).
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0003_c.jpg"})
Figure 3. The comparative RMSD plots of IMP-α1 protein in apo (black) and complex form with IVM molecule within its both NLS-Major (red) and NLS-Minor site (green).
3.2.2. Radius of gyration
To analyze the comparative compactness and stability of IMP-α1 protein induced by binding of IVM within major and minor site, we measured radius of gyration (RoG) of target protein which calculates the mass-weighted positional distances of receptor atoms from the center of mass. RoG evaluate the level of compactness of the protein structure, the lower the RoG values, the more compact and stable is the protein, while higher values of RoG demonstrates high conformational entropy and disorderness in the protein(Naz et al., 2019). Figure 4 illustrated that the IVM/NLS-Major site showed higher structural stability during the initial 50 ns indicated by smaller Rg values around 28 ± 0.4 Å, but it immediately increased to a sharp peak around 30 ± 0.4 Å during the last 20 ns (from 60 to 80 ns). While on the other side, the IVM/NLS-Minor site showed higher Rg value of 29.5 ± 0.4 Å during the initial 30 ns but after that this system was observed to be stabilized around 29.0 Å in the later 50 ns which might be attributable to establish and maintain more electrostatic and hydrophobic contacts as compared to major site. While the apo protein due to the absence of IVM, undergo higher fluctuation with more disorderness in the system starting from 28.5 and reached to 30.0 ± 0.3 Å, while the later 40 ns displayed frequent Rg values around 29.5 ± 0.5 Å. It is quite visible from the graph (Figure 4) that association of IVM with NLS-Minor site decreased the conformational entropy of IMP-α1 protein and enhanced its compactness and stability during last 50 ns as compared to the other NLS-Major site which showed less compactness and higher entropy.
Published online:
05 May 2021Figure 4. Graphical representation of radius of gyration for three systems of IMPα1 i.e. apo (black), IVM/NLS-Major (red), and IVM/NLS-Minor (green) binding site with reference to time (ns).
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0004_c.jpg"})
Figure 4. Graphical representation of radius of gyration for three systems of IMPα1 i.e. apo (black), IVM/NLS-Major (red), and IVM/NLS-Minor (green) binding site with reference to time (ns).
3.2.3. Hydrogen bond analysis
To better understand the core factors responsible for the stability and compactness induced by IVM after binding with major and minor site of IMP-α1 protein, the total number of hydrogen bonds formed and deformed by IVM with the significant residues of major and minor sites were calculated through the entire simulation period of 80 ns. From the results (Figure 5), it was observed that IVM/NLS-Major complex established almost 32 different hydrogen bonds which were formed and deformed during simulation. Among all these, only seven bonds involving residues Ser149, Gly150, Thr155, Asp192, Asn239, Tyr277 and Asp280 showed stability with maximum occupancy of 13.1%. An important thing to note here was the conformational shift observed by the IVM molecule despite of establishing various hydrogen bonds. IVM displayed a shift in its dock conformation highlighting its less stability within NLS-Major site during whole simulation.
Published online:
05 May 2021Figure 5. The IVM molecule shift from its dock site to deeper into the NLS-Major binding site (left) which results in the formation of more stable and consistent hydrogen bonds as shown in occupancy graph (right).
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0005_c.jpg"})
Figure 5. The IVM molecule shift from its dock site to deeper into the NLS-Major binding site (left) which results in the formation of more stable and consistent hydrogen bonds as shown in occupancy graph (right).
Comparatively, the average number of hydrogen bonds observed between IVM and NLS-Minor site were low around 22 which further drops to five on the basis of their stability and occupancy calculation (Figure 6). The top five most stable bond involved residues Asp280, Gly281, Asn361, Ala364 and Ser406. Among these Asn361 being the most crucial residue responsible for the inhibition of this protein presented the most stable bond with maximum occupancy of 8.4% which is two times less than the occupancy observed in major binding site.
Published online:
05 May 2021Figure 6. IVM molecule remain at the same dock site of NLS-Minor binding site (left) during simulation with very weak and inconsistent hydrogen bond illustrated in occupancy graph (right).
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0006_c.jpg"})
Figure 6. IVM molecule remain at the same dock site of NLS-Minor binding site (left) during simulation with very weak and inconsistent hydrogen bond illustrated in occupancy graph (right).
The smaller number of hydrogen bonds between IVM and NLS-Minor along with more stable plot of RMSD and RoG encouraged us to investigate further the key factors responsible for the higher stability of IVM/NLS-Minor site over major site complex. Therefore, we visually inspected all the simulated trajectories at molecular level to find an answer to this mystery. After visualizing the trajectories, we came to know that despite of more hydrogen bonds established between IVM and NLS-Major site, it showed a shift in the conformation of IVM molecule from its dock site. This shifting of IVM removed its contact with one of the crucial residues, i.e. Ala148 and reduced its overall hydrophobic contacts. This shifting showed a drastic change in its interaction pattern observed in our docking study which might be a reasonable explanation for the less stability of IVM in this pocket during whole simulation. On the other hand, despite of few hydrogen bonds between IVM and NLS-Minor site, molecular visualization of its output trajectories displayed a more stable conformation of IVM molecule which maintained the same docked position by establishing more hydrophobic and van der Waals interactions. Unlike IVM/NLS-Major site complex, IVM molecule remained and maintained its dock conformation well within the NLS-Minor pocket by developing more hydrophobic contacts as compared to electrostatic interactions. Since IMP-α1 inhibition required more hydrophobic contacts than electrostatic interactions, therefore, IVM/NLS-Minor site presented more stability over major site that might be attributable to more hydrophobic contacts between IVM and NLS-Minor site. While the shift in the dock conformation of IVM molecule within NLS-Major site modified its interaction pattern, i.e. less hydrophobic contacts and more hydrogen bonds which was responsible for the overall instability IMP-α1 protein (highlighted by high RoG and RMSD values above) after association of IVM with NLS-Major site.
3.2.4. Molecular reasons for the stability of major over minor site
The ambiguous picture displayed by the stability plots and molecular interaction pattern; we further verify the factors responsible for the stability of IVM molecule within NLS-Minor site as compared to NLS-Major site. In this regard, we visualized the simulated trajectories of whole simulation at molecular level and saved the snapshot of IVM within both sites after every 20 ns time step. It helps to investigate the number of hydrogen bonds remain stabilized among IVM and major or minor site through entire simulation.
As depicted in Figure 7, with the passage of every nano second, the electrostatic interactions between IVM and NLS-Major site were kept on increasing but with insignificant amino acids which did not contributed toward the stability of this complex. In all the snapshots, during late 40 ns, IVM established a hydrogen bond with crucial residue Asp192 at a distance of 2.05 ± 0.5 Å while IVM shifting brought forth a salt bridge between IVM and Arg238 during late 20 ns. While the hydrophobic contacts between IVM and NLS-Major site residues were reduced, we observed only three hydrophobic contacts among IVM and Arg238, Lys240 and Tyr277 during the initial 30 ns. All these hydrophobic contacts were lost during the late 50 ns because of its shifting from its initial dock position.
Published online:
05 May 2021Figure 7. 3-dimensional interaction diagram of IVM with the crucial residues of NLS-major site after every 20 ns of entire simulation.
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0007_c.jpg"})
Figure 7. 3-dimensional interaction diagram of IVM with the crucial residues of NLS-major site after every 20 ns of entire simulation.
Comparatively, the second system comprising of IVM with NLS-Minor site showed more hydrophobic contacts with very few electrostatic interactions. As shown in Figure 8, IVM established one or two hydrogen bonds throughout whole simulation with minor site. It was observed that throughout 80 ns simulation, only one hydrogen bond was observed with crucial residue Asn361 during late 20 ns which is attributable to less occupancy of crucial hydrogen bonds between IVM and minor site as mentioned earlier. As compared to electrostatic interactions, the NLS-Minor site residues established more stable hydrophobic contacts with IVM molecule involving the crucial amino acids like Trp357, Ala364 and Trp399 which remain consistent through entire simulation and required for the inhibition of IMP-α1. From all these results in hand, we inferred that the hydrophobic contacts between both sites and IVM molecules were the key factors responsible for the overall stability of IMPα protein.
Published online:
05 May 2021Figure 8. Protein–ligand interaction diagram of IVM (green) with residues of NLS-minor site (yellow) after every 20 ns of entire simulation.
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0008_c.jpg"})
Figure 8. Protein–ligand interaction diagram of IVM (green) with residues of NLS-minor site (yellow) after every 20 ns of entire simulation.
3.2.5. Root mean square fluctuation
To probe the local minimal changes at the atomic or residual level, Root mean square fluctuation was calculated. In this regard, the RMSF of each backbone (C-α) atoms of the complex (IVM/NLS-major and IVM/NLS-minor) and apo structure were calculated (Figure 9). As shown in figure, the fluctuations are stronger (∼2.1 ± 0.6) in IMP-α1 (Apo) because of the absence of IVM. Whereas, the RMSF plot of IVM/NLS-Major site exhibited less stability with greater fluctuations around 2.1 ± 0.2 Å which might be attributable to the shifting of hydrophobic contacts to electrostatic interactions with IVM. As depicted in the plot, despite of more hydrogen bonds observed between IVM and NLS-Major site, it still showed some peak areas involving the residues of major binding site and also minor site. The fluctuation score clearly demonstrates the instability of IMP-α1 protein in complex with the ligand (IVM).
Published online:
05 May 2021Figure 9. Root mean square fluctuation plot of apo target protein (black) and in complex with IVM displayed the overall stability of NLS-Major site (red) over minor site (green).
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/tbsd20/0/tbsd20.ahead-of-print/07391102.2021.1906750/20210505/images/medium/tbsd_a_1906750_f0009_c.jpg"})
Figure 9. Root mean square fluctuation plot of apo target protein (black) and in complex with IVM displayed the overall stability of NLS-Major site (red) over minor site (green).
Comparative to major site, the association of IVM with minor site of IMP-α1 showed low fluctuations around 1.8 ± 0.2 Å. Our results showed that the binding of IVM molecule enhances the overall stability of IMP-α1 and NLS-Minor site by displaying least fluctuations in the binding residues (350–400). Therefore, despite of very few hydrogen bonds, IVM/NLS-Minor site showed higher stability with less conformational entropy after binding with IVM might be attributable to the formation of more hydrophobic contacts during whole 80 ns simulation.
3.2.6. Binding-free energy
Binding-free energy is the basic tool to analyze the change in ligand binding patterns by measuring its thermodynamic properties. Negative value of binding free energy (ΔGbinding) highlighted the stability of system while the positive values demonstrate instability (Table 2) summarized the binding energy components of MM-PBSA calculation which surprisingly showed opposite results from docking. IVM/NLS-Minor complex showed higher negative values of −151.5 ± 16.7 KJ/mol as compared to −104.2 ± 21.6 KJ/mol value of IVM/NLS-Major complex. All the other energetic components of these complexes were seemed to be in line with our docking and MD results highlighted by more negative value of electrostatic energy for IVM/NLS-Major site complex as compared to the minor one. The dislocation of IVM from its dock site because of the shallow and wide nature of major site brought forth a shift in its interaction pattern, i.e. less hydrophobic contacts and more hydrogen bonds that might be attributable to its less negative value of van der Waals energy and high negative value of electrostatic energy (Table 2). Because of this shift, IVM maintained more hydrogen bonds but lost its hydrophobic contacts with crucial residues required for IMP-α1 inhibition through-out 80 ns simulation. The greater number of hydrogen bonds between IVM and NLS-Major site was responsible for its higher negative value of electrostatic energy. While the minor pocket being deep and narrow, holds the compounds well in place by establishing more van der Waals and hydrophobic contacts instead of electrostatic interactions, highlighted by its low negative value for electrostatic energy and high negative value for van der Waals energy (Table 2). Since the inhibition of IMP-α1 requires more hydrophobic contacts with crucial of residues of binding pockets (major or minor). The similar pattern of interactions was observed in NLS-Minor site bound with IVM molecule. These results highlighted the possible inhibitory mechanism of IVM by binding with NLS-Minor binding site of IMP-α1 protein that eventually blocks the entry of SARS-CoV-2 protein within the nucleus and halts its transcription and replication.
Table 2. Summary of MM-PBSA binding energy components of both complex systems, i.e. IVM in complex with major and minor site.
3.2.7. Per-residue decomposition energy
Residue-based decomposition energy works on the principle of decomposition of whole binding free energy into smaller per-residue segments based on their electrostatic and hydrophobic terms which helped to signify specific hotspot residues that were responsible for increasing the binding affinity of ligand toward its protein. Since most of the significant amino acids involved in inhibition of IMP-α1 protein are of hydrophobic nature either they belong to NLS-Major or minor binding site. Therefore, we calculate the decomposition energy of distinct charged and hydrophobic residues involved in inhibition of IMP-α1 by establishing number of stable hydrogen bonds and hydrophobic contacts with IVM molecule. The residues that are stabilized by ligand showed negative values while the residues presenting positive energy values are not stabilized by ligand. In our study, results highlighted that the residues belonged to NLS-Major binding site displayed higher positive values of energy (63.736 to −0.552 Kcal/mol) Table 3 as they involved in establishing and maintaining more hydrogen bonds than hydrophobic contacts with IVM molecule. The residue Lys240 displayed the highest positive decomposition energy value as it did not maintain its hydrophobic contacts with IVM molecule because of shift in the molecule position within the major binding site. While the residues of NLS-Minor binding site presented less positive values and more negative values (0.731 to −0.056 Kcal/mol) (Table 3) for decomposition energy which is attributable to the establishment of more stable hydrophobic contacts as compared to electrostatic interactions which is attributable to the stable binding of IVM ligand within minor binding site. The higher negative values displayed by minor site with IVM is responsible for its overall stability and less conformational entropy with crucial residues of NLS-minor site. The per-residue decomposition energy displayed the higher affinity of IVM toward minor site as compared to major site because of the presence of strong and stable hydrophobic and van der Waals interactions between IVM and crucial residues of minor binding site (Asp280, Ile286, Trp399 and Gly407) through entire simulation. All these results highlighted the possible inhibitory mechanism of IVM via binding and establishing stable hydrophobic contacts with the crucial residues of NLS-Minor site than major site.
Table 3. Residue-based decomposition energy of the crucial residues of both NLS-major and NLS-minor binding site in complex with IVM molecule.
4. Conclusion
Ivermectin being an FDA approved drug, has well-established safety profile to be used in humans against various parasitic infections. Recent studies and reviews reported comparable safety of ivermectin at high doses as compared to its standard low dose except its safety profile in pregnancy is not evident yet. In recent study, Caly et al. reported antiviral activity of Ivermectin against SARS-CoV-2 (∼5000-fold reduction of viral RNA) with no toxicity. They hypothesized that this reduction in SARS-CoV-2 RNA is most likely through inhibition IMPα/β1 mediated nuclear import of viral proteins which needs further work up. Therefore, in our study, we explored the comparative binding mode and inhibitory mechanism of Ivermectin against two NLS-Major and NLS-Minor binding sites of IMP-α1. Our results discovered that the ivermectin reduces SARS-CoV-2 viral transport by inhibiting IMP-α1 protein after binding with its NLS-Minor site. The highest negative score in binding-free energy calculation brought forth the stable binding of IVM with NLS-Minor site as compared to NLS-Major site which was further confirmed during 80 ns simulation by investigating different stability parameters. During simulation, IVM stabilizes NLS-Minor site by establishing and maintaining more hydrophobic contacts than hydrogen bonds which is required for the inhibition of IMP-α1 protein. These hydrophobic contacts were not observed in case of NLS-Major site because of shifting from its stable conformation to unstable conformation with the pocket residues. Altogether, these results highlighted the worth of ivermectin along with its safety profile, to be considered as a possible drug to limit SARS-CoV-2 viral load and reduces the progression of this disease.
| S. no. | Docking method | Replacement methods | Scoring function | Re-scoring function | Dock score (kcal/mol) |
|---|---|---|---|---|---|
| Induce Fit | Triangle Matcher | London Dg | London Dg | –13.3721 (Major Site) –11.6671 (Minor Site) | |
| Triangle Matcher | London Dg | London Dg | –12.6937 | ||
| Triangle Matcher | Affinity Dg | London Dg | –11.4481 | ||
| Rigid | Triangle Matcher | .. | London Dg | –11.1412 | |
| Triangle Matcher | .. | Gbwsa | –8.2175 | ||
| Alpha Triangle | .. | .. | –7.6018 | ||
| Alpha PMI | Affinity Dg | Gbwsa | –7.1437 | ||
| Triangle Matcher | London Dg | Gbwsa | –6.1887 |
Except first combination, all other scores were reported for major site only.
| Energies | IVM/NLS-Major | IVM/NLS-Minor |
|---|---|---|
| Van der Waals energy | –180.9 ± 26.5 KJ/mol | –235.3 ± 18.8 KJ/mol |
| Electrostatic energy | –30.2 ± 8.7 KJ/mol | –17.5 ± 9.4 KJ/mol |
| Solvation energy | 124.5 ± 19.4 KJ/mol | 124.9 ± 20.4 KJ/mol |
| Total binding energy | –104.2 ± 21.6 KJ/mol | –151.5 ± 16.7 KJ/mol |
| IVM/NLS-Major site residues | IVM/NLS-Minor site residues | ||
|---|---|---|---|
| Residue no. | Energy Kcal/mol | Residue no. | Energy Kcal/mol |
| Ser149 | –0.003 | Asp280 | –0.001 |
| Asn188 | –0.005 | Gly281 | –0.003 |
| Asp192 | –0.019 | Ile286 | –0.056 |
| Trp231 | –0.103 | Asp325 | –0.008 |
| Arg238 | –0.131 | Trp357 | 0.731 |
| Asn239 | –0.552 | Asn361 | 0.005 |
| Lys240 | 63.736 | Glu396 | –0.012 |
| Tyr277 | –0.107 | Trp399 | –0.018 |
| Arg315 | –0.367 | Gly407 | –0.023 |
Acknowledgments
The authors also thank the Deanship of Scientific Research at King Saud University, Saudi Arabia for funding through the research group project no. RGP-1438-043.
Disclosure statement
The authors declare no competing financial interests.
Author contributions
SN and ZU hypothesized and designed the project. SM carried out the benchmarking and docking. UQ set up the MS system and carried out the simulations. SN, SM and UQ drafted the manuscript. ZU, M.N and SA facilitate the research work and reviewed the manuscript. All authors have read and approved the manuscript.