Computational drug repurposing study of antiviral drugs against main protease, RNA polymerase, and spike proteins of SARS-CoV-2 using molecular docking method

Introduction

In the concluding weeks of 2019, a newly emerged human coronavirus (SARS-CoV-2) which is the etiological agent responsible for the 2019–2020 viral pneumonia outbreak reported in Wuhan, China [1, 2]. The novel coronavirus has been recognized as a new member of the betacoronavirus genus which is closely associated with bat coronavirus and severe acute respiratory syndrome coronavirus (SARS-CoV) [3]. According to the recent situation report from the World Health Organization on September 9, 2020, over 27 million cases confirmed positive for coronavirus disease of 2019 (COVID-19) and 894,983 deaths are registered [4]. Unfortunately, there are no approved drugs or vaccines for this disease and effective treatment options remain very limited and several clinical trials are in progress [5]. Novel coronavirus (2019-nCoV) is an enveloped, positive-sense, single-stranded RNA beta-coronavirus [6]. Similar to severe acute respiratory syndrome (SARS) and Middle East respiratory syndrome (MERS), the 2019-nCoV genome encodes two groups of protein, structural proteins, including spike (S), matrix (M), envelope (E), and nucleocapsid proteins (N), and nonstructural proteins such as RNA-dependent RNA polymerase (RdRp) and the main protease (Mpro) [7, 8]. Four structural proteins are essential for virion assembly and infection of SARS-CoV-2. This virus utilizes a densely glycosylated spike (S) protein composed of two subunits, S1 and S2, to gain entry into the host cell. The S protein is a trimeric class I fusion protein that exists in a metastable prefusion conformation that faces a structural rearrangement to fuse the viral membrane with the host cell membrane [9, 10]. This process is provoked when the S1 subunit binds to a host-cell receptor. SARS-CoV-2 uses angiotensin-converting enzyme 2 (ACE2) as its receptor [11]. After binding of the S1 subunit to the ACE2 receptor [12, 13], a conformational change in the S protein results in the exposure of a hydrophobic fusion peptide. The fusion peptide penetrates to the target cell membrane, initiating the virion-cell membrane fusion process. Fusion peptide insertion into the target cell membrane, followed by joining of two heptad repeats in S2 forming an antiparallel six-helix coiled-coil bundle. The formation of this bundle, also known as the trimer of hairpins, facilitates the fusion of viral and cellular membranes, leading to the ultimate release of the viral core into the cytoplasm of the host cell [10, 14]. At the next step, the genomic (+) single-stranded RNA (ssRNA) is not only used as a template for viral genome replication by RdRp but also translated directly to polyproteins. These two polyproteins encode nonstructural proteins to form the replication–transcription complex. The produced polyproteins are cleaved by virally encoded proteases, such as the Mpro, at specific sites to release the proteins used in SARS-CoV-2 replication [15-17].

Computational pharmacology approaches can be utilized for the design and discovery of drugs through drug repurposing methods [18, 19], which saves cost and time in the development of effective antiviral agents for clinical use [3, 20]. Drug repositioning of the approved pharmaceutical drugs offers an alternative strategy for the rapid identification of the potential drugs challenging novel diseases and manages emerging viral infections [21, 22]. Approved drugs for some diseases are safe for human use [5], and solely their effectiveness against the disease of interest needs to be determined [23]. However, clinical trials are necessary to ensure the results [24].

Various studies on the inhibition of SARS-CoV-2, through targeting different virus proteins with antiviral drugs, which are approved and used for other viral infections, have been performed [17, 19, 25], [26], [27]. In a study, the inhibitory effect of remdesivir, which is used to control the Ebola and SARS viruses, on the SARS-CoV-2 RdRp was investigated. Remdesivir bound to the active site of the RdRp enzyme with an estimated binding energy of −7.6 kcal/mol and contradicted the polymerase function [28]. In an in silico study conducted in 2020 by Kumar et al. the inhibitory potential of the drug Lopinavir–Ritonavir, which is mainly used to control human immunodeficiency virus (HIV) infection, on the Mpro of the SARS-CoV-2 virus was investigated. This compound was able to inhibit the active site of the protease enzyme with an estimated binding energy of −10.6 kcal/mol [29]. In another report, the Tegobuvir drug, a non-nucleoside inhibitor for the hepatitis C virus, has been shown to have a potential inhibitory effect on the SARS-CoV-2 spike protein, with an estimated binding energy of −8.1 kcal/mol [30]. Due to the key roles of the spike, Mpro, and RdRp proteins in SARS-CoV-2 infection, in the current study, these proteins were selected as a target for molecular docking-based virtual screening of more than 50 antiviral drugs. Hence, the study aims to evaluate the effect of antiviral drugs on blocking cleavage sites of spike protein, suppression of template binding, polymerization, and nucleoside triphosphate binding domains in RdRp, and also, inhibition of protease activity of the Mpro. Also, the drugs were assessed for their simultaneous inhibitory effect on all three target proteins. The predictions of this study will provide information that can be utilized for the selection of candidate drugs for in vitro or in vivo studies and clinical trials.

Methods

Results

The PDB file of Mpro (PDB ID: 6LU7), spike (PDB ID: 6VSB), and RdRp (PDB ID: 6M71) proteins were obtained from PDB. The structure of the protease protein active site area, including T24, T25, T26, F140, G143, S144, C145, H163, H164, E166, D187, R188, Q189, and T190 residues was demonstrated in (Figure 1A). Also, (Figure 1B) illustrates the spike protein accompanied by the cleavage site of fusion peptide which includes K811, R815 residues. The three regions of the RdRp protein active site including, the polymerization site with the FSMMILSDDAVVCFN motif sequence, the template strand binding site with the DKSAGFPFNKWGKAR motif sequence, and the nucleotide triphosphate (NTP) binding site with the LKYAISAKNRAR motif sequence, are marked on the protein structure in (Figure 1C).

Figure 1:

Three-dimensional structure of severe acute respiratory syndrome coronavirus target proteins (cartoon format). (A) Main protease; specified area depicts the protease active site, (B) spike protein; the marked area shows the cleavage site of the fusion peptide, and (C) RNA-dependent RNA polymerase protein; yellow (template strand binding site), red (nucleotide-binding site), and blue (polymerization site). Figures were produced using UCSF Chimera [36].

During the modeling process, some bad contacts, unsuitable torsion angles and, etc. will produce inevitably. According to the observations that native protein structures correspond to a system at thermodynamic equilibrium with minimum free energy, hence the most straightforward way to increase the accuracy of a partially incorrect protein structure model is to search for the lowest-energy structure, also known as native structure [37], [38], [39]. Therefore, the energy levels of all three target proteins were minimized by 10,000 times by applying the CHARMM force field [40]. The energy level for spike, Mpro, and RdRp reached from −2,922.26 kcal/mol, −6,247 kcal/mol, and −4,794.346 kcal/mol to −57,077.606 kcal/mol, −22,577.139 kcal/mol, and −73,863.421 kcal/mol, respectively. The results show that all three proteins are reached their most stable energy levels and are ready to enter the next step.

Protein–ligand docking is a widely used computational tool in the drug designing and discovery process as it allows the prediction of the conformation and orientation of the most favorable structure of the complex formed between a given protein (target) and ligands within the binding site of a target molecule. Indeed, molecular docking studies have been employed to calculate structural parameters such as binding free energy (ΔG) of the drug molecule, their estimated inhibition constant (K _i), etc. [41]. The AutoDock applies restricted flexibility in the protein as a receptor and computes the interaction energy between ligand and receptor using a Lamarckian genetic algorithm [42, 43]. Hence, it requires grid maps that describe the field of interaction energies around the interested target macromolecule. At the end of each docking run, the docked conformation of the ligand, which represents a specific state of the ligand, is characterized by an estimated free binding energy. Also, AutoDock reports van der Waals energy and electrostatic energy for each atom of the ligand. Estimated free binding energy (ΔG_binding) is calculated from the following formulas [44], [45], [46]:

As mentioned before, Autodock 4.2 software was used for the virtual screening of antiviral drugs against spike, Mpro, and RdRp proteins. More than 50 drugs belonging to the category of antiviral drugs were downloaded from the Drug Bank database as PDB files. This library (Table 1) was used to screen the target proteins.

The comparison of the docked cocrystal ligand position to Mpro and complex protein (6LU7) confirmed the validation of docking operations. Since the conformations and locations of both forms were perfectly superimposed, ligands bound to their active site are shown in (Figure 2). Also, three compounds including, phenol, acetate, and benzene were docked to the target proteins as negative controls.

Figure 2:

Representation of conformation and location of (A) the cocrystal molecule and (B) the docked cocrystal molecule to the protease (6LU7). The estimated binding energy of docked cocrystal molecules was −7.43 kcal/mol. Figures were produced using UCSF Chimera [36].

The results of molecular docking of antiviral drugs against protease as protease inhibitors are listed in (Table 1). Indinavir, lopinavir, and saquinavir, respectively, with the estimated binding energy of −12.14 kcal/mol, −12.11 kcal/mol, and −11.75 kcal/mol and with estimated K _i less than 3 nM, have shown to have the highest potential and affinity to protease inhibition. Since, the lower the estimated inhibition constant (K _i) of the inhibitors, the greater their tendency to inhibit proteins at lower concentrations, therefore, by considering these low values, it can be claimed that the mentioned inhibitors have a significant potential to inhibit the protein.

As demonstrated in (Figure 3A), indinavir formed a hydrogen bond with the Glu166 residue with the length of 1.876 Å and two hydrogen bonds with the Gln189 residue with the length of 2.610 Å and 2.254 Å. In terms of position and active site blocking, the indinavir is in a good condition. By forming two hydrogen bonds with Glu166 and Gln189 residues with lengths of 1.767 Å and 2.117 Å, respectively, lopinavir filled the active site space (Figure 3B). By establishing two hydrogen bonds with the length of 2.235 Å and 1.773 Å with Asn162 and Glu166 residues, respectively, saquinavir blocked the protease active site with proper estimated binding energy (Figure 3C). (Table 2) shows the energy of the docked compounds, including indinavir, lopinavir, and saquinavir to the SARS-CoV-2 Mpro.

Figure 3:

The antiviral drugs with the highest potential and affinity to protease inhibition. The binding pose of (A) the indinavir, (B) the lopinavir, and (C) the saquinavir with the active site of the protease. The red spiral refers to the hydrogen bond. The ligands are marked with a green margin. Figures were produced using UCSF Chimera [36].

Also, the results of molecular docking of antiviral drugs as spike inhibitors against viral spike are shown in (Table 1). In this protein, the fusion peptide is considered as a binding site after the proteolysis process. Three drugs, including remdesivir, saquinavir, and lopinavir, with the estimated binding energy of −8.15 kcal/mol, −7.93 kcal/mol, and −7.78 kcal/mol, respectively, are the best potential inhibitors of the fusion peptide in the spike protein. The estimated inhibition constant (K _i) of each drug is less than 3 μM, which indicates the high affinity of the inhibitors to the binding site.

Remdesivir binds to the fusion peptide region by forming one hydrogen bond with the Lys811 residue with the length of 2.066 Å, and two hydrogen bonds with the Asp808 residue with the lengths of 1.782 Å and 2.045 Å (Figure 4A). Although saquinavir has shown a low estimated binding energy for fusion peptides compared to other antiviral drugs, no hydrogen bond has been established between saquinavir and the fusion peptide region in the spike protein (Figure 4B). Lopinavir binds to the fusion peptide region by forming a hydrogen bond with the Lys811 residue with a length of 1.707 Å (Figure 4C). The data in Table 2 shows the energy of the three antiviral drugs, including remdesivir, saquinavir, and lopinavir, that docked to the SARS-CoV-2 spike protein.

Figure 4:

The antiviral drugs as the best potential inhibitors of the fusion peptide in the spike protein. The binding pose of (A) the remdesivir, (B) the saquinavir, and (C) the lopinavir to the fusion peptide region in the spike protein. The saquinavir, despite the low estimated binding energy for fusion peptide, has been established with no hydrogen bond. The red spiral refers to the hydrogen bond, and the ligands are marked with a green margin. Figures were produced using UCSF Chimera [36].

According to the results of molecular docking of antiviral drugs as RdRp inhibitors against RdRp, which are shown in (Table 1), three regions including polymerase site, nucleotide-binding (NTP binding) site, and template strand binding site were targeted in the grid box. Three drugs, composed of, lopinavir, saquinavir, and nelfinavir, with the estimated binding energy of −9.75 kcal/mol, −9.73 kcal/mol, and −9.06 kcal/mol, respectively, are the best potential inhibitors of the RdRp. The estimated inhibition constant (K _i) of each of these drugs is 71.52 nM, 73.43 nM, and 227.62 nM, respectively. Despite the lowest estimated binding energy of the lopinavir–RdRp complex compared to other antiviral drugs, no hydrogen bond has been established between lopinavir and the polymerization site of RdRp protein (Figure 5A). Instead, due to the presence of phenylalanine residues around the lopinavir, the main interaction is van der Waals force. Saquinavir has been docked to the polymerization site of the RdRp by forming three hydrogen bonds with Ser861, Asp865, and Arg836 residues with lengths of 2.282 Å, 2.100 Å, and 1.853 Å, respectively (Figure 5B). Nelfinavir is bound to the polymerization site of RdRp through the formation of a hydrogen bond with the Lys593 residue with the length of 1.948 Å, two hydrogen bonds with Thr591 with the lengths of 3.647 Å and 2.605 Å, and two hydrogen bonds with Ser592 residue with the lengths of 1.854 Å and 2.783 Å (Figure 5C). The data in Table 2 reports the energy of three drugs, including lopinavir, saquinavir, and nelfinavir in binding to RdRp.

Figure 5:

The drugs with the lowest estimated binding energy against RNA-dependent RNA polymerase (RdRp). The binding pose of (A) the lopinavir, (B) the saquinavir, and (C) the nelfinavir in the polymerization site of RdRp. The red spiral refers to the hydrogen bond, and the ligands are marked with a green margin. Figures were produced using UCSF Chimera [36].

By examining the structure of the RdRp–drug complexes of the three mentioned drugs, it was found that the drugs have the highest affinity for binding to the NTP binding site of this protein rather than the two other sites. Also, to comparing the drugs with NTPs, docking was performed, and the results showed that the affinity of drugs to RdRp was significantly higher than NTPs. The results of molecular docking of antiviral drugs against RdRp are shown in (Table 1), as well as, (Table 3) reports the lowest estimated binding energy of the three mentioned drugs and four triphosphate nucleotides along with the number of times that the ligands docked to different regions of RdRp. The comparison of the inhibitory potential of RdRp protein by three drugs, including lopinavir, saquinavir, and nelfinavir, and relative to ATP, UTP, GTP, and CTP nucleotides in inhibition of polymerization binding sites, NTP binding site, and template binding site, indicated very satisfying results and the estimated inhibitory binding energy of all three mentioned drugs in the binding to all three sites of RdRp protein was several times higher than four nucleotides.

Discussion

A novel SARS-CoV-2 was identified from pneumonia patients in December 2019 over the last year, it has emerged as a serious threat to world public health [47]. However, there is no approved drug to effectively inhibit the virus so far. Similar to most other coronaviruses, the membrane spike glycoproteins, composed of S1 and S2 subunits, is the prime viral interacting protein with host cell targets (such as ACE2 and other cell adhesion factors) [10]. The S protein mediates subsequent fusion between the envelope and host cell membranes to aid viral entry into the host cell. The spike is a determinant for cell adhesion and virulence [48]. Once the viral genome enters the host cytoplasm, it is used as a template for viral genome replication by RdRp [16]. Also, it is translated to the polyproteins, which cleaved at specific sites by viral proteases, including the Mpro [31]. Therefore, considering the decisive role of these proteins in virus pathogenicity and its life cycle, these three proteins were selected as the target proteins in the present study.

The timely development of effective antiviral agents for clinical use is extremely challenging because conventional drug development approaches normally take years of investigation and cost billions of dollars [20]. Therefore, repurposing approved pharmaceutical drugs provides an alternative approach to rapidly identify potential drugs, which have specific pharmacokinetic and toxicological properties [49]. Relying on this topic and repurposing concept, a virtual screening procedure employing docking of more than 50 FDA approved antiviral drugs was used to identify new potential small molecule inhibitors against three major proteins in SARS-CoV-2, including Mpro, RdRp, and spike. In addition to evaluating the inhibitory potential of each protein by different antiviral drugs, another goal of this study was to find a drug with the inhibitory potential of all three proteins that could be used to control the pathogenicity of the SARS-CoV-2 at different stages. Based on the results of molecular docking of antiviral drugs against all three target proteins, it is observed that five drugs including, indinavir, lopinavir, saquinavir, nelfinavir, and remdesivir have the highest inhibitory potential for all three proteins compared to other antiviral drugs. During the COVID-19 pandemic, remdesivir has been identified as a potential inhibitor of RdRp protein. By considering the consistent results of the present study with other reports [28, 50] on the remdesivir inhibitory effect on RdRp protein, the results indicate that remdesivir is more effective in the inhibition of the fusion peptide in the spike protein. Although remdesivir binds to the RdRp with an appropriate estimated binding energy of −8.67 kcal/mol, by sorting the ligands in terms of their binding energy, lopinavir, saquinavir, and nelfinavir were still in higher ranks. In this study, among the studied antiviral drugs, lopinavir and saquinavir were identified as the best inhibitor of all three proteins by considering their binding energies. Lopinavir is an approved drug for inhibiting the protease of HIV-1 infection, including its physicochemical properties, 628.8 g/mol molecular weight, Log P 5.9, 15 rotatable bonds, four hydrogen bond acceptors, and five hydrogen bond donors (Pubchem CID 92727). Also, saquinavir is another proven HIV-1 protease inhibitor with physicochemical properties including 670.8 g/mol molecular weight, Log P 4.2, 13 rotatable bonds, seven hydrogen bond acceptors, and five hydrogen bond donors (Pubchem CID 441243). According to the Drug Bank, lopinavir (DB01601) and saquinavir (DB01232) have a low oral bioavailability of ∼25 and 4%, respectively. To increase oral bioavailability, lopinavir and saquinavir co-administration with ritonavir as an enzyme inhibitor responsible for the first-pass metabolism of lopinavir and saquinavir. It is noteworthy that in this study, ritonavir also showed acceptable molecular docking results with binding energies of −9.72 kcal/mol, −7.69 kcal/mol, and −8.16 kcal/mol for Mpro, spike, and RdRp, respectively. Based on this, it can be predicted that this combination drug regimen could have an acceptable effect in inhibiting SARS-CoV-2 pathogenicity. Also, investigating the relationship between the structure of ligands and their inhibitory potential through analyzing hydrogen bonds, is as follows:

1. The amide group in the three drugs including lopinavir, saquinavir, and indinavir is most effective in forming hydrogen bonds with the Mpro active site.

2. The hydroxyl groups in lopinavir and remdesivir show more involvement in the formation of hydrogen bonds with the residues of the fusion peptide region in spike protein.

3. The hydroxyl and amide groups in nelfinavir and the amide group in saquinavir have the greatest effect on the formation of hydrogen bonds with residues of the RdRp active site.

Therefore, we hypothesize that the presence of amide groups in the ligand structure will contribute to the design or discovery of new drugs with antiprotease properties of the SARS-CoV-2. Also, if the goal is to disrupt the virus entry into the host cell, hydroxyl groups in the ligand structure will be effective in inhibiting the fusion peptide region. The three drugs including lopinavir, saquinavir, and nelfinavir have a high affinity for binding to the NTP binding site and the polymerization site, compared to the strand binding site in the RdRp protein. The pattern of hydrogen bonds of the amide group in saquinavir and nelfinavir, the hydroxyl group in nelfinavir, and the van der Waals interactions between the lopinavir aromatic rings and the benzene groups of the active site’s phenylalanine residues reflects the diversity of different groups in binding to different sites of RdRp. We suggest that the amide, hydroxyl, and aromatic rings can be present in the ligand structure to design or discover a drug against SARS-CoV-2 RdRp.

Conclusion

The present study successfully identified three antiviral drugs as potential inhibitors for each target protein of SARS-CoV-2, including Mpro, spike, and RdRp. According to estimated binding energy (kcal/mol), and estimated inhibition constant (K _i), and investigation of hydrogen bond of the mentioned five drugs in (Table 2), especially lopinavir and saquinavir, showed the greatest inhibitory potential of the target proteins in total. Since these compounds are FDA-approved and successfully passed various toxicology studies, therefore there is a hope that these drugs may be considered as potential inhibitors for spike protein, RdRp, and Mpro of SARS-CoV-2, which can be further explored to test against the disease in pre-clinical and clinical trial studies as a two-drug treatment regimen.