In Silico Comparative Analysis of Ivermectin and Nirmatrelvir Inhibitors Interacting with the SARS-CoV-2 Main Protease

1. Introduction

Emerging in Wuhan, China, towards the end of 2019, a new acute respiratory illness caused by the novel coronavirus (SARS-CoV-2) was officially named COVID-19. The World Health Organization (WHO) declared it a global pandemic on March 11, 2020 [1]. This virus is known for its high transmissibility. It primarily spreads through human-to-human contact via contaminated respiratory secretions released during coughing, sneezing, or talking [2,3]. When in close proximity, aerosol particles can also transmit the infection, initiating the viral replication process in a new host [4].

Previously, smaller-scale epidemics, such as those in 2002, 2003, and 2012, were linked to the Coronavirus, impacting various hosts, including bats, rodents, mammals, and birds [5,6] and by 2024, just over four years after the initial case in Wuhan, SARS-CoV-2 infections affected nearly 704 million people, with approximately 7 million resulting in fatalities [7]. It is worth noting the well-known challenges in official reporting, which may lead to underestimating these values [8].

SARS-CoV-2 comprises three structural proteins and 16 non-structural proteins, among which the ${\mathrm{M}}^{\mathrm{pro}}$, or ${\mathrm{3CL}}^{\mathrm{pro}}$, holds a central role in the virus’s replication mechanism, acting as a cysteine protease [9,10] — positioned strategically at the active site of ${\mathrm{M}}^{\mathrm{pro}}$, a cysteine acts like molecular scissors, facilitating the cleavage of viral polyproteins. These lengthy protein chains, synthesized from the viral genome, are broken down into smaller, functionally independent units crucial for viral replication, assembly, and release. Given the vital importance of ${\mathrm{M}}^{\mathrm{pro}}$, various research groups are currently actively engaged in the development of targeted inhibitors, aiming for potential therapeutic applications for COVID-19 [11]. Several promising candidates are in advanced stages of clinical trials, raising optimism about their eventual incorporation into the antiviral arsenal against this persistent pathogen [12].

In this context, Alugubelli and colleagues explored the potential of Boceprevir, a known inhibitor of HCV’s NSP3, as a repurposed treatment for COVID-19 [13]. This compound exhibited promising inhibitory effects on the ${\mathrm{M}}^{\mathrm{pro}}$ of SARS-CoV-2, featuring an $\alpha$-ketoamide group as a warhead, along with specific structural elements including a P1 $\beta$-cyclobutylalanine segment, a P2 dimethyl cyclopropyl proline motif, a P3 test-butyl glycine unit, and a P4 N-terminal start-butyl carbamide. Through modifications at all four positions, 20 ${\mathrm{M}}^{\mathrm{pro}}$ inhibitors based on Boceprevir were synthesized [13]. One such derivative, PF-07321332, underwent characterization for its potency against ${\mathrm{M}}^{\mathrm{pro}}$ both in vitro and in cellulo, using test tubes and 293T cells, respectively. The study involved analyzing the crystal structures of ${\mathrm{M}}^{\mathrm{pro}}$ bound to 10 of these inhibitors, alongside assessing the cytotoxicity and antiviral efficacy of four selected inhibitors.

By substituting the P1 site with a $\beta$-(S-2-oxopyrrolidin-3-yl)-alanil (Opal) residue and altering the warhead to an aldehyde, a high in vitro potency was achieved [13]. Additionally, the original structural components in P2, P3, and the P4 N-terminal positions of Boceprevir demonstrated superior performance compared to alternative chemical groups tested for their in vitro potency. In crystal structures, all inhibitors formed a covalent adduct with the cysteine at the active site of ${\mathrm{M}}^{\mathrm{pro}}$[13]. Specifically, the Opal P1 residue, P2 dimethylcyclopropylproline, and P4 N-terminal tert-butylcarbamide established robust hydrophobic interactions with ${\mathrm{M}}^{\mathrm{pro}}$, elucidating the high in vitro potency associated with compounds containing these structural features.

Andi and colleagues conducted a study focused on developing novel antivirals and vaccines targeting ${\mathrm{M}}^{\mathrm{pro}}$[14]. Their findings, supported by crystallographic evidence and binding assay data, reveal that three drugs approved for hepatitis C virus treatment, along with two drug-like compounds, establish covalent bonds with the catalytic residue Cys145 within the active site of ${\mathrm{M}}^{\mathrm{pro}}$. Furthermore, molecular docking studies offered additional insights, aiding in the formulation of new antiviral inhibitors for SARS-CoV-2 by leveraging these drugs as lead compounds. It is worth mentioning that numerous other investigations into potential ${\mathrm{M}}^{\mathrm{pro}}$ inhibitors have been documented in the literature [15,16,17,18].

The substrate binding site of ${\mathrm{M}}^{\mathrm{pro}}$, a region known for its remarkable conservation across various coronavirus genera, was investigated by Li and colleagues [19]. They specifically examined Nirmatrelvir (PF-07321332), an oral inhibitor developed by Pfizer that targets the ${\mathrm{M}}^{\mathrm{pro}}$ of SARS-CoV-2, with demonstrated effectiveness against other coronaviruses as well. Through their analysis of these structures, the researchers unveiled a conserved binding site shared among coronaviruses, offering valuable insights into the mechanism underlying the inhibition of viral replication [19].

Viewing it from this angle, Nirmatrelvir specifically targets ${\mathrm{M}}^{\mathrm{pro}}$ [20]. By inhibiting this protease at a critical stage of proteolysis before viral RNA replication, Nirmatrelvir effectively curbs the proliferation of SARS-CoV-2. It has exhibited potent antiviral activity against various human coronaviruses and notable off-target selectivity [20]. Additionally, animal experiments have confirmed its favorable safety profile and demonstrated its capacity to reduce lung viral load in a mouse model engineered to simulate SARS-CoV-2 infection [20]. Preclinical investigations have found no evidence of mutagenic interactions with DNA [21]. Moreover, a phase I clinical trial involving healthy volunteers showed plasma concentrations surpassing the threshold of in vitro cellular antiviral potency [21].

Conversely, Ivermectin is a broad-spectrum antiparasitic medication that initially garnered attention due to its in vitro ability to inhibit SARS-CoV-2 replication, sparking interest in its potential as a COVID-19 treatment [22,23,24]. However, its off-label use for this purpose is advised cautiously. Available in both human and veterinary formulations, Ivermectin is widely used clinically to combat parasites like Onchocerca volvulus (which causes river blindness) and Strongyloides stercoralis. Nonetheless, at elevated doses, it can lead to potentially severe adverse effects such as nausea, vomiting, diarrhea, dizziness, and even fatalities. Large-scale clinical trials investigating its efficacy against COVID-19 have yielded mixed results. While some studies fail to demonstrate significant benefits in reducing symptoms or duration, others suggest potential advantages [25,26]. Consequently, a consensus based on consistent, high-quality evidence regarding Ivermectin’s role in treating COVID-19 has yet to emerge in the literature.

Driven by the potential therapeutic implications, we conducted an in silico study to investigate the interactions between Ivermectin and Nirmatrelvir with ${\mathrm{M}}^{\mathrm{pro}}$. Employing molecular docking, we first examined how Ivermectin binds to ${\mathrm{M}}^{\mathrm{pro}}$, specifically focusing on its interaction with protomer A, denoted as ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$. Subsequently, we conducted molecular dynamics simulations for both Ivermectin and Nirmatrelvir bound to ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, facilitating a comprehensive analysis of their interaction profiles at a molecular level.

2. Metodology

2.1. Electronic Properties

We investigated the electronic properties of Ivermectin and Nirmatrelvir using density functional theory (DFT) and frontier molecular orbital (FMO) analyses within the Gaussian 16 software suite [27]. Molecular geometries were optimized through the semi-empirical Parameterization Method 6 (PM6) coupled with a steepest descent algorithm [28]. Following optimization, single-point energy calculations utilizing restricted DFT [29] were conducted to determine the energies of the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), and their energy gap [30]. Additionally, Mulliken population analysis and electric charge distribution were performed. For a comprehensive DFT analysis, Becke’s three-parameter Lee-Yang-Parr (B3LYP) exchange-correlation functional [31] and the def2-TZVP (valence triple-zeta polarization) basis set [32] were employed.

2.2. Molecular Docking

Ligands and proteins were prepared using AutoDock Tools 1.5.6 (ADT) [33]. Molecular docking simulations were employed to examine the non-covalent binding interactions between the ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$[20] (PDB: 7VH8) complex and Ivermectin. The crystallographic structure from the PDB contained the complex of ${\mathrm{M}}^{\mathrm{pro}}$ and the inhibitor Nirmatrelvir. Autodock Vina, integrated into PyRx, was integrated with the genetic algorithm set as the default configuration for the docking simulations. The characterization of interactions between Ivermectin and ${\mathrm{M}}^{\mathrm{pro}}$ was facilitated using ChimeraX and Discovery Studio programs. The simulation box encompassing the entire protein had dimensions of $60\times 40\times 65$ ${\AA }^3$, with a center point at (-1.74, 18.92, -26.05) Å. Ligand screening was performed with an accuracy of ±2 Å for ligand positions and ±0.01 Kcal/mol for binding affinities.

2.3. Molecular Dynamics

Molecular dynamics (MD) simulations were conducted using the Nanoscale Molecular Dynamics (NAMD) 2.9 computational package [34] along with the CHARMM36 force field [35], using a time step of 1.0 fs. The model protein/ligand complexes were maintained under constant pressure and temperature conditions within the NPT ensemble at 300 K and 1 atm, respectively. Equilibration occurred throughout 10 ns, followed by a 100 ns production run. To assess the temporal evolution of the system, Root-mean-square deviation (RMSD) and Fluctuation (RMSF) calculations were performed using Visual Molecular Dynamics (VMD) 1.9.1 [36]. Subsequently, re-docking was executed after assembling interactions derived from the initial docking poses.

3. Results

Initially, DFT calculations were performed to understand the behavior of molecules when isolated, particularly concerning the calculated frontier molecular orbitals (HOMO and LUMO). Figure 1 displays, on the left side, the results for HOMO and LUMO of the Ivermectin molecule, while, on the right side, it presents the same characteristics for Nirmatrelvir.

The figure analysis reveals favorable features for the interaction of Ivermectin with its target protein, indicating a potential enhancement of bioactivity. The high HOMO energy of -6.08 eV, coupled with a low LUMO energy of -0.96 eV, results in a significant energy gap of 5.12 eV for Ivermectin. In the case of Nirmatrelvir, the HOMO energy is -6.75 eV, and the LUMO energy is -1.65 eV. Consequently, the energy gap is 5.10 eV. This wide gap shows an efficient charge transfer within Ivermectin, potentially facilitating stable interactions with the target protein. Moreover, the high energy gap implies strong inhibition efficiency due to the minimum energy required for electron removal from the HOMO.

It is observed that ideal inhibitory molecules readily accept and donate electrons, exhibiting an electron-rich character associated with superior inhibition potential. Consistent with these principles, the previously discussed calculated LUMO energy for Ivermectin shows a good electron-accepting capacity, further supporting its potential as an effective inhibitor. Regarding the overall electronic characteristics of both Ivermectin and Nirmatrelvir, we find they are similar in both molecules. This trend implies that Ivermectin, like Nirmatrelvir, can potentially be used as an inhibitor ligand for ${\mathrm{M}}^{\mathrm{pro}}$, hindering the replication process of SARS-CoV-2.

In molecular docking, we prioritized ligand poses exhibiting more favorable binding affinity, as lower $\Delta$G free energy scores indicated, reflecting the most probable and stable binding modes [37]. The docking simulations revealed favorable binding affinities for Ivermectin, with the top-scoring pose presenting a $\Delta$G of -9.0 kcal/mol.

In this perspective, Figure 2 details the ligand-${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ interactions. At the top, we analyze Nirmatrelvir. This complex was obtained from Zhao and colleagues [20]. We observed that the molecule establishes interactions near the N-terminus of ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, encompassing methionines 49 and 165 (MET), histidine 41 (HIS), phenylalanine 140 (PHE), glycine 143 (GLY), threonine 190 (THR), cysteine 145 (CYS), glutamate 166 (GLU), glutamine 189 and 192 (GLN), and proline 168 (PRO). On the other hand, at the bottom of this figure, Ivermectin, in the reported best energy configuration, occupies a distinct pocket closer to the C-terminus of ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, interacting with alanine 154 (ALA), aspartic acid 197 (ASP), leucines 286 and 287 (LEU), methionine 176 (MET), and proline 168 (PRO). Remarkably, while Nirmatrelvir and Ivermectin bind to separate pockets near the N- and C-termini, respectively, they converge towards the central region, as evidenced by their shared interaction with proline 168 (PRO). The dominant interactions of Nirmatrelvir involve conventional H-bonds. In contrast, those of Ivermectin are dominated by alkyl and pi-alkyl interactions, showing that Nirmatrelvir interacts more strongly with ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$.

Building upon the molecular docking results for Ivermectin and the Nirmatrelvir/${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ complex discussed earlier, MD simulations of these systems were conducted to explore the interaction dynamics of both molecules with ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ under ambient temperature and pressure. This aimed to verify their stability in the presence of external factors and their persistence as inhibitors throughout the simulation.

Figure 3 illustrates the temporal evolution of RMSD (top left panel) and RMSF (top right panel) simultaneously for Nirmatrelvir (blue curves) and Ivermectin (red curves) over the 100 ns simulation time. The lower panels depict the number of h-bonds for Nirmatrelvir (left panel) and Ivermectin (right panel). The RMSD values for best pose score for ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ with Nirmatrelvir are 2.0 Å, and with Ivermectin, it is 1.87 Å. These values, hovering around 2 Å, suggest that the reference structures are highly similar to the initial structures, signifying that both ligands remain securely bound to the cavity throughout the simulation time. The RMSD curves exhibit notable similarity, characterized by low fluctuations, showing that Ivermectin has the potential to act as an effective ligand against ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$. Moreover, in the interval between 25 ns and 45 ns of simulation, our results indicate that Nirmatrelvir experiences more pronounced RMSD fluctuations than Ivermectin, followed by increased stability beyond this point. These results imply a more pronounced stability of ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ in complex with Ivermectin. RMSF, on the other hand, analysis was employed to assess the dynamic behavior of the ${\mathrm{M}}^{\mathrm{pro}}$ protein when bound to Nirmatrelvir and Ivermectin. Both docking complexes exhibited minimal amino acid fluctuations throughout the 100 ns simulations, as illustrated by their RMSF plots. Thus, both ligands significantly stabilize the ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ protein.

In terms of h-bond numbers, Figure 3 (lower panels) shows that the ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ complex with Ivermectin exhibited superior stabilization over 100 ns compared to the complex with Nirmatrelvir. The temporal progression of hydrogen bond formation between the protein and ligands was also scrutinized to glean insights into molecular interactions. The analysis uncovered that Nirmatrelvir initially formed four hydrogen bonds with the ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ protein, with two and three persisting and alternating throughout the simulation. In contrast, Ivermectin formed two to three initial hydrogen bonds, with only one remaining stable throughout the simulation.

From the above-discussed MD simulations, we emphasize two snapshots at 50 ns and 100 ns (see Figure 4). These visualizations unveil a notable difference in conformational flexibility, wherein both ligands maintain binding to ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ throughout the simulations. However, Nirmatrelvir exhibits a high conformational change compared to Ivermectin.

Finally, Figure 5 delineates the interaction profiles of Nirmatrelvir and Ivermectin with ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ at 50 ns and 100 ns, representing the molecular dynamics snapshots discussed earlier. At the 50 ns mark, Nirmatrelvir (right-upper panel) engages with amino acids, including methionine 49 and 165 (MET), histidine 41, 163, and 164 (HIS), phenylalanine 140 (PHE), cysteine 145 (CYS), glutamate 166 (GLU), and asparagine 142 (ASN). In the subsequent 100 ns interval (left-upper panel), the Nirmatrelvir/${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ complex maintains interactions with methionine 49 and 165 (MET), histidine 41, glutamate 166 (GLU), asparagine 142 (ASN), glutamine 192 (GLN), and leucine 167 (LEU). Conversely, Ivermectin’s interactions persist over 50 ns and 100 ns (lower panels) with amino acids tyrosine 237 and 239 (THR) and leucine 271, 272, and 287 (LEU). Notably, Ivermectin exhibits a more enduring interaction pattern with the same amino acids over time, while Nirmatrelvir interacts with approximately half of the amino acids at 50 ns compared to 100 ns. Consistent with the docking results, Nirmatrelvir predominantly relies on conventional hydrogen bonds, whereas Ivermectin primarily engages in alkyl and pi-alkyl interactions. Consequently, the more robust nature of hydrogen bonds contributes to Nirmatrelvir forming a higher number of interactions with ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$.

4. Summary and Conclusions

In conclusion, this study aims to investigate the distinct binding modes exhibited by Nirmatrelvir and Ivermectin with the SARS-CoV-2 ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, with the perspective of understanding the nature of these molecules as inhibitors in the viral replication process. In this regard, Nirmatrelvir demonstrates engagement with the ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, forming numerous interactions involving amino acids such as methionine 49 and 165 (MET), histidine 41, 163, and 164 (HIS), phenylalanine 140 (PHE), cysteine 145 (CYS), glutamate 166 (GLU), and asparagine 142 (ASN). Through molecular dynamics simulations, the RMSD analysis reveals a stable binding of Nirmatrelvir, with an RMSD value of around 2.0 Å, indicating its secure attachment to the binding pocket over the 100 ns simulation period.

On the other hand, Ivermectin exhibits good stability, maintaining prolonged interactions with specific amino acids, such as tyrosine 237 and 239 (THR) and leucine 271, 272, and 287 (LEU). The RMSD for the Ivermectin-${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$ complex is 1.87 Å, reflecting its consistent and secure binding throughout the simulation. It is worth noting that Ivermectin interacts with approximately the same set of amino acids at both 50 ns and 100 ns, demonstrating enduring interaction patterns over time.

The diversified exploration of ligand-protein interactions emphasizes the intricate nature of these binding profiles. Nirmatrelvir, relying on conventional hydrogen bonds, engages in a significant number of stronger interactions with ${\mathrm{M}}^{\mathrm{pro}\left(\mathrm{A}\right)}$, contributing to its strong binding. However, Ivermectin primarily utilizes alkyl and pi-alkyl interactions, indicating a different binding mode.