Evaluation of potential anti-RNA-dependent RNA polymerase (RdRP) drugs against the newly emerged model of COVID-19 RdRP using computational methods.

INTRODUCTION: Despite all the efforts to treat COVID-19, no particular cure has been found for this virus. Since developing antiviral drugs is a time-consuming process, the most effective approach is to evaluate the approved and under investigation drugs using in silico methods. Among the different targets within the virus structure, as a vital component in the life cycle of coronaviruses, RNA-dependent RNA polymerase (RdRP) can be a critical target for antiviral drugs. The impact of the existence of RNA in the enzyme structure on the binding affinity of anti-RdRP drugs has not been investigated so far.
METHODS: In this study, the potential anti-RdRP effects of a variety of drugs from two databases (Zinc database and DrugBank) were evaluated using molecular docking. For this purpose, the newly emerged model of COVID-19 (RdRP) post-translocated catalytic complex (PDB ID: 7BZF) that consists of RNA was chosen as the target.
RESULTS: The results indicated that idarubicin (IDR), a member of the anthracycline antibiotic family, and fenoterol (FNT), a known beta-2 adrenergic agonist drug, tightly bind to the target enzyme and could be used as potential anti-RdRP inhibitors of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). These outcomes revealed that due to the ligand-protein interactions, the presence of RNA in this structure could remarkably affect the binding affinity of inhibitor compounds.
CONCLUSION: In silico approaches, such as molecular docking, could effectively address the problem of finding appropriate treatment for COVID-19. Our results showed that IDR and FNT have a significant affinity to the RdRP of SARS-CoV-2; therefore, these drugs are remarkable inhibitors of coronaviruses.

Introduction

COVID-19 has caused significant fatalities all around the world, and to date, no particular medicinal drug has been discovered for the treatment of this virus. There is no suitable vaccine or therapy for this deadly virus so far, and this could be attributed to the transient nature of coronaviruses. It has also been reported that clinical results differ with various groups and ages [1]. Coronavirus is a member of the Coronaviridae family, which is named for the crown-like spikes on its surface [2]. According to the nucleotide sequences of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), it is a member of Betacoronaviruses, such as SARS and MERS human coronaviruses (HCoVs) [3,4]. Seven different strains of HCoVs have been reported so far, including the 229E and NL63 strains of HCoVs (Alphacoronaviruses), and the OC43, HKU1, SARS, MERS, and SARS-CoV-2 HCoVs (Betacoronaviruses) [4,5]. Two groups of proteins characterize HCoVs: structural proteins like spikes marking all coronaviruses; and non-structural proteins, such as RNA-dependent RNA polymerase (RdRP) [6]. As a viral enzyme, RdRP is a vital component in the life cycle of RNA viruses, which has been targeted in various viral infections, including the Zika virus (ZIKV), hepatitis C virus (HCV), and coronaviruses (CoVs) [[7], [8], [9], [10], [11]]. Nucleotide and nucleoside analogs with polymerase inhibition activity are considered as an essential group of antiviral agents [12]. The active site of RdRP is thoroughly protected with two surface-accessible and successive aspartates in a beta-turn structure [[13], [14], [15]]. The structure of nsp12, the RdRP of SARS-CoV-2, which is bound to nsp7 and nsp8 with an overall view on the complex architecture has been previously specified [16], but the mechanism of RNA recognition and synthesis has yet to be determined. One of the most significant barriers to finding nucleotide-based drugs and their reliable evaluation is the absence of molecular detail on the substrate identification through the replication process. To address this problem, Wang et al. have specified the near-atomic-resolution arrangement of SARS-CoV-2 polymerase in complex with RNA when it is in the catalytic condition [17], which provides a considerable opportunity for the optimization of drug design and finding dependable drug candidates.

In the battle against coronavirus, scientists have utilized three different approaches for developing drugs [18]. The first approach is to use the current wide range of antivirals [19]. Although the metabolic features, employed dosages, potential effectiveness, and adverse effects of these therapies are evident, they cannot destroy coronaviruses in a targeted way because of their broad spectrum; also their unfavorable effects cannot be neglected [20]. The second approach is to find molecules that could have therapeutic effects on coronaviruses by screening the molecules present in the existing molecular databases [21,22]. Due to the high-throughput screening of this method, many new operations for therapeutic agents could be determined. In the third approach, coronaviruses are categorized based on their genomic details and pathological features for developing new targeted drugs from the beginning [20]. Although the drugs that will be found through this strategy are more targeted for coronaviruses, the procedure may take several years or even more than ten years [23].

There are also some state-of-the-art developments in discovering drugs against COVID-19. These approaches could affect this area of research considerably. One of these approaches was utilized by Hu et al. They conducted a normal mode analysis (NMA) of spike protein to comprehend the connection between their nanomechanical features and virus lethality as well as infection rate. Their investigation revealed that epidemiological virus properties could be directly connected to pure nanomechanical features [24]. Panda et al. used a structure-based drug designing and immunoinformatics approach to evaluate the effects of different antiviral drugs against the spike glycoprotein, main protease, and receptor-binding domain (RBD) of SARS-CoV-2 [25]. Another crucial approach to manage the COVID-19 pandemic is preparing sensitive diagnostic systems. Electrochemical SARS-CoV-2 biosensing supported by artificial intelligence is one consequential part of these systems, which was adequately reviewed by Kaushik et al. [26]. Furthermore, the role of nanomaterials and antimicrobial nano-coatings against COVID-19 is indisputable and was thoroughly addressed by Rai et al. [27]. In another study, the CRISPR/Cas13 system as a potential approach for the treatment of COVID-19 was discussed, which could attract considerable attention in the near future [28].

The fastest way to find therapeutic agents against coronaviruses is to identify the potential therapeutics among the approved or developing drugs. For this purpose, computational studies such as molecular docking or molecular dynamics could be used as beneficial strategies [17,29,30].

In this study, we used molecular docking to evaluate the potential of different drugs to target the post-translocated catalytic complex of COVID-19 RdRP (PDB ID: 7BZF), which is the reaction center of RdRP. For this screening, five drugs were chosen from the literature [17,31], and 23 more drugs with potential RdRP inhabitation were selected from a Food and Drug Administration (FDA)-approved drug database (Zink drug database, ZDD) [20] and the DrugBank database.

Materials and methods

Sequence alignment and homology

To gain a better insight into the recently emerged COVID-19 RdRP post-translocated catalytic complex, the sequence of this practical model was aligned with the structure of SARS-CoV-2 RdRP in complex with cofactors in reduced conditions (PDB ID: 7BTF), which is the template of RdRP post-translocated catalytic complex [17].

Protein preparation

The structure of COVID-19 RdRP post-translocated catalytic complex (PDB ID: 7BZF) with an electron microscopy method and the resolution of 3.26 Å was downloaded from the PDB database (http://www.rcsb.org/pdb). Maestro molecular modeling platform (version 12.2) by Schrödinger program was used for these calculations [32]. To prepare the structure, the protein preparation module [33] from the Schrödinger Suite was used to add the hydrogen atoms, remove the waters beyond 5 Å from the binding sites, and optimize the structure for creating an H-bond network. Finally, the energy was minimized using optimized potentials for liquid simulations (OPLS3e) force field with a default setting of 0.30 Å root mean square deviation (RMSD). The overall stability and stereochemical quality of the protein in the 3D structure were determined using the Ramachandran plot [34].

Prediction of active sites

Before performing docking, the identification and characterization of the binding site should be carried out for a more reliable and accurate molecular docking. The active sites for target-ligand binding interaction were predicted using the SiteMap module of Schrödinger software (SiteMap, Maestro molecular modeling platform (version 12.2) by Schrödinger program). SiteMap can determine the binding sites in a large-scale validation with best results for the sites that bind ligands with the subnanomolar association, and in addition, it can propose an adjusted version of the score that accurately classifies the druggability of proteins [35].

Generation of the receptor grid

The receptor grid was developed using Glide (Glide, Schrödinger, generated LLC, New York, NY, 2020) by preferring the predicted active sites. The box was centroid of the active site of the protein receptor with the default parameters of van der Waals scaling factor of 1 Å, partial charge cut-off at 0.25, and docked ligand length of 20 Å.

Ligand preparation

Twenty-eight compounds from different databases and literature were identified as test ligand molecules against the COVID-19 RdRP post-translocated catalytic complex. The preparation of ligands was carried out utilizing the LigPrep module of Schrödinger [36], and the OPLS3e force field was used for generating the ionization and tautomeric states at pH 7.0 ± 2.0 using the Epik module of LigPrep. The name and structures of the chosen drugs are presented in Table 1 .

Table 1

Potential anti-RdRP drugs from the Zinc database, DrugBank database, and literature.

Remdesivir [17] <image>	Favipiravir [17] <image>	Ribavirin [6] <image>
Galidesivir [6] <image>	Tenofovir [6] <image>	Setrobuvir [6] <image>
CHEMBL180948 (DugBank) <image>	CHEMBL221046 (DugBank) <image>	R1626 (DugBank) <image>
Sofosbuvir (DugBank) <image>	Valganciclovir (Zinc) <image>	Diphenoxylate (Zinc) <image>
Ceftibuten (Zinc) <image>	Fenoterol (Zinc) <image>	Fludarabine (Zinc) <image>
Itraconazole (Zinc) <image>	Cefuroxime (Zinc) <image>	Atovaquone (Zinc) <image>
Chenodeoxycholic acid (Zinc) <image>	Cromolyn (Zinc) <image>	Pancuronium bromide (Zinc) <image>
Cortisone (Zinc) <image>	Tibolone (Zinc) <image>	Novobiocin (Zinc) <image>
Silybin (Zinc) <image>	Idarubicin (Zinc) <image>	Bromocriptine (Zinc) <image>
Benzylpenicilloyl G (Zinc) <image>

Virtual screening and molecular docking

Rigid docking

High-throughput virtual screening and molecular docking of the chosen ligands against the target protein was done using the Glide module of the Schrödinger Suite. The prepared output files of LigPreb were used as input files for screening and combined with redistribution for sub-jobs. The unique compounds were identified by generating unique properties for each input compound, and the Epik state penalties were used for docking. Standard-precision (SP) and extra-precision (XP) docking calculations were performed. The flexible method was utilized for docking, and 10% of the best compounds were kept after each step with retaining all good scoring states. The values of the scaling factor and partial charge cut-off were set at 0.80 and 0.15, respectively. The ligands that had a sufficient association with the RdRP were authorized. In the final step, the pose viewer was used for analyzing the interactions of the selected ligands and the docked protein complex.

Induced-fit docking

Considering the fact that rigid docking has some limitations as the protein is kept rigid [37,38], the accuracy of its results has to be further evaluated. The induced-fit docking module of Schrödinger was utilized to accomplish induced-fit docking [39,40]. In brief, the ligand is docked into the active site of the protein utilizing the Glide protocol (SP) [41], while the active site residues of the protein are maintained rigid. Next, the protein side chains or backbone are refined using the prime refinement module [42,43]. Finally, the ligand is redocked into the refined protein conformation, and the score is calculated to rank the protein-ligand complex. The centroid of the residues was selected for the box center option, and the residues of active sites were selected. For the conformational sampling option, the sample ring conformation was chosen with the energy window of 2.5 kcal/mol. The receptor van der Waals scaling and the ligand van der Waals scaling were set to 0.7 and 0.5, respectively. The residues were refined within 5 Å of the ligand poses, and Glide redocking was performed on the structures within 30 kcal/mol of the best structure.

Results

Sequence alignment

Not surprisingly, as shown in Fig. 1 , the COVID-19 RdRP post-translocated catalytic complex was highly conserved between SARS-CoV-2 RdRP in complex with cofactors in reduced conditions. The homology of these two structures was as high as 97%.

Fig. 1

Sequence alignment of a post-translocated catalytic complex of COVID-19 RdRP among SARS-CoV-2 RdRP in complex with cofactors.

Reliability of protein structure and prediction of active site

The quality and the reliability of the target protein in the 3D form are of importance for drug design. The Ramachandran plot was obtained using the Schröodinger Suite to display the allowed and disallowed regions. As can be seen in Fig. 2 , more than 98.9% of the residues were in the allowed regions (95% in the favored regions and 3.9% in the allowed regions), and only 1.1% of the residues were present in the disallowed regions implying the reliability of the structure.

Fig. 2

Ramachandran plot of 7BZF. The orange regions are the most favored regions, the yellow regions are the allowed regions, and the white regions are the disallowed regions. Glycine is plotted as triangles, proline is plotted as squares, and all the other residues are plotted as circles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The active sites of the protein are presented in Fig. 3 . Since the site score of 1 and above confirms the high druggability of a site [44], the first two sites with the site scores of 1.040 and 1.019 were chosen as the target regions. The details of the five active sites are given in Table 2 . In Fig. 3B and C the second active site with the size of 551 that indicates different parts of the second site is presented. The green regions show the hydrophilic domains of the site, while the blue and red regions are hydrogen-bond donors and hydrogen-bond acceptors, respectively. The residues that are included in the active sites are presented in Table 3 .

Fig. 3

A: The top-ranked active sites of the protein (white pocket). B: Magnification of the second site that shows different domains of this site. It presents the regions of hydrophilic domains (green regions) and hydrophobic ones (yellow regions). C: Hydrogen-bond donors (blue regions) and hydrogen-bond acceptors (red regions). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2

Description of different active sites.

Site No.	Site score	size	D score
Site 1	1.040	130	1.064
Site 2	1.019	551	0.942
Site 3	1	117	1.019
Site 4	0.995	199	1.031
Site 5	0.985	136	0.990

Table 3

Residues in the first and second active sites.

Site No.	Residues
Site 1	Arg249, Leu251, Thr252, Ser255, Tyr265, Ile266, Trp268, Leu270, Val315, Leu316, Ser318, Thr319, Val320, Phe321, Pro322, Pro323, Thr324, Ser325, Phe326, Arg349, Glu350, Thr394, Cys395, Phe396, Tyr456, Arg457, Asn459, Leu460, Pro461, Thr462, Pro627, Asn628, Met629, Ser664, Val675, Pro677
Site 2	Val410, Lys411, His439, Phe440, Phe441, Phe442, Gln444, Asp452, Tyr455, Tyr456, Ala490, Asn491, Val493, Ile494, Val495, Asn496, Asn497, Leu498, Asp499, Lys500, Lys511, Arg513, Tyr516, Tyr521, Gln524, Met542, Leu544, Lys545, Tyr546,Ala547, Ile548, Ser549, Ala550, Lys551, Arg553, Ala554, Arg555, Thr556, Val557, Ala558, Lys577, Ala580, Ala581, Arg583, Gly590, Gly616, Trp617, Asp618, Tyr619, Pro620, Lys621, Cys622, Asp623, Arg624, Thr680, Ser681, Ser682, Thr687, Ala688, Asn691, Ser759, Asp760, Asp761, Ala762, Lys798, Trp800, His810, Glu811, Phe812, Cys813, Ser814, Gln815, Pro832, Asp833, Arg836, Ala840, Val844, Asp845, Lys849, Arg858

Drug binding to COVID-19 RdRP post-translocated catalytic complex

Two active sites of the protein were selected to be used as the target regions. The first active site was chosen because it had the highest activity (site score: 1.040), and the next one was chosen due to its size (size: 551). However, the site score of the second site was heightened too. The results of molecular docking are presented in Table 4 . The top three docking scores of rigid docking in the first site belonged to idarubicin (IDR) (−6.424), fludarabine (−6.401), and cromolyn (−6.110). The top three docking scores in the second site belonged to IDR (−8.611), fenoterol (FNT) (−8.562), and silybin (−7.786). As can be seen, IDR had the most binding affinity to both sites. The ligand among the residues of the active site is shown in Fig. 4 . The ligand interaction of IDR is presented in Fig. 5 . As can be seen, the ligand had direct contact with Arg 555 and Lys 500. The docking scores of induced-fit docking revealed that FNT is comparable to IDR because the docking score of FNT is higher than that of IDR in the second site. In Fig. 6 FNT is presented among the residues and in Fig. 7 the interactions of FNT with RdRP are described.

Table 4

Docking scores of anti-RdRP drugs against two active sites with the highest site score and largest size.

No.	Drug name	Docking score of rigid docking. First site	Docking score of rigid docking. Second site	Docking score of induced-fit docking. First site	Docking score of induced-fit docking. Second site
1	Remdesivir	−5.591	−4.364	−5.413	−6.718
2	Favipiravir	−3.827	−5.179	−6.323	−4.822
3	Ribavirin	−5.600	−6.729	−6.490	−6.953
4	Galidesivir	−3.214	−7.580	−6.187	−7.154
5	Tenofovir	−4.397	−4.994	−5.241	−5.795
6	Setrobuvir	−5.052	−0.705	−4.601	−6.088
7	CHEMBL 180948	−4.073	−3.622	−3.598	−3.870
8	CHEMBLE 221046	−2.320	−3.273	−3.929	−3.561
9	R1626	−4.757	−3.286	−3.506	−4.781
10	Sofosbuvir	−5.376	−4.343	−6.918	−6.529
11	Valganciclovir	−5.884	−6.822	−6.485	−5.085
12	Ceftibuten	−4.718	−5.896	−6.302	−5.105
13	Fenoterol	−4.934	−8.562	−8.521	−7.561
14	Fludarabine	−6.401	−5.299	−5.937	−7.042
15	Itraconazole	−4.023	−1.070	−6.203	−5.025
16	Cefuroxime	−4.672	−5.711	−4.961	−6.135
17	Atovaquone	−3.926	−3.454	−4.882	−4.103
18	Chenodeoxycholic acid	−3.189	−4.807	−6.113	−4.951
19	Cromolyn	−6.110	−4.420	−5.917	−4.559
20	Dabigatran etexilate	−5.291	−3.975	−6.375	−6.367
21	Cortisone	−5.298	−4.088	−5.110	−5.820
22	Tibolone	−3.354	−3.914	−5.109	−5.109
23	Novobiocin	−5.203	−1.010	−5.492	−3.407
24	Silybin	−4.922	−7.786	−5.625	−6.997
25	Idarubicin	−6.424	−8.611	−7.806	−7.993
26	Bromocriptine	−3.763	−4.865	−6.708	−6.661
27	Diphenoxylate	−2.878	−5.486	−5.273	−4.528
28	Benzylpenicilloyl G	−5.639	−4.202	−5.179	−6.922

Fig. 4

The 3D structure of idarubicin between the residues of the 7BZF. Dashed yellow lines present the hydrogen bonds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

The interactions between idarubicin and the protein. Dashed lines in Ligplot indicate H-bonds. Carbons are in black, nitrogens are in blue, and oxygens are in red. Semicircles are hydrophobic contacts. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

The 3D structure of FNT between the residues of the 7BZF. Dashed yellow lines present the hydrogen bonds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7

The interactions between idarubicin and the protein. Dashed lines in Ligplot indicate H-bonds. Carbons are in black, nitrogens are in blue, and oxygens are in red. Semicircles are hydrophobic contacts. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

Despite all the efforts since the emerging of new COVID-19, no particular drug with a direct effect against this coronavirus has been discovered or created so far. The transient nature of the cause of this global disaster could be considered as a practical reason for this setback. Since developing antiviral treatment is time-consuming, the screening of approved drugs could be a practical approach for finding a probable treatment for coronavirus. In recent years, in silico screening of bioactive components based on molecular docking has become a common strategy in the pharmaceutical industry [45]. Moreover, the high-speed development of computer technology and the rapid increase of biological data on therapeutic targets, make it understandable that the use of computational methods and in silico approaches has increased remarkably [[46], [47], [48], [49]]. In fact, in silico approaches enable the possibility of virtual screening among millions of compounds, which leads to an affordable cost of discovering the desired drug candidates [50]. In this study, potential anti-RdRP drugs and some antiviral drugs were chosen to be evaluated against the newly emerged COVID-19 RdRP post-translocated catalytic complex. The active sites of the structure were distinguished to find more reliable results, which could indicate whether a candidate site could make a tight ligand binding or not. The SiteMap module of Schrödinger software could provide this opportunity. In addition, the SiteMap results could categorize the sites into the druggable and undruggable sites [51], which gives a better insight into the physical basis of classification. The sequence alignment of 7BZF and its template, 7BTF, showed that the gaps between these two structures were less than 3%. The significant gap was from glycine 897 to aspartate 910, which consisted of Gly897, His898, Met899, Leu900, Asp901, Met902, Tyr903, Ser904, Val905, Met906, Leu907, Thr908, Asn909, and Asp910. Theses residues were the same in SARS-COV-2, Bat-coronavirus-RaTG13, SARS-COV, Pangolin-coronavirus, and BtRs-BetaCoV/YN2018D, but the residues 908 and 909 in MERS-CoV were cysteine and glycine, respectively [17]. For the molecular docking, the Glide module of the Schrödinger Suite was utilized, which provides more reliable results than a variety of different docking programs such as AutoDock, GOLD, and FlexX [52]. After a high-throughput screening and considering the docking score of different drugs, IDR exhibited a tight binding affinity to the host protein. IDR is a member of the anthracycline antibiotic family, which is usually used for the treatment of acute myeloid leukemia [53]. It has been reported that IDR can interrupt the DNA topoisomerase IIA subunits [54] and effectively block the viral protein and RNA synthesis [55]. Charak and Mehrotra, using UV–visible spectroscopy, observed that IDR gets intercalated between the DNA bases [56]. The docking results against two different sites of 7BZF showed that IDR could tightly bind to the RdRP, especially at the second active site. There were different interactions between IDR and 7BZF. IDR created hydrogen bonds with two residues Lys500 and Arg555, with the bond lengths of 2.80 and 3.11 Å, respectively. The only protein residue that had a hydrophobic interaction with IDR was Arg836. The most important thing about this site was that it contained RNA, and IDR could form multiple bonds with RNA. IDR and RNA created four hydrogen bonds and a hydrophobic interaction. Bases A0 and C1 formed hydrogen bonds with the lengths of 2.84 and 2.80 Å, respectively. Base A4 formed two hydrogen bonds with the lengths of 3.01 and 3.06 Å, respectively. The only base with hydrophobic interaction was C3. It is worth noting that the hydrogen bond of A4, hydrophobic interaction of C3, and hydrogen bond of Arg555 were created with oxygen number 9 of IDR. The position of this oxygen is presented in Fig. 5. FNT formed six hydrogen bonds, and only one of these hydrogen bonds was with residues. The rest of the hydrogen bonds were formed between the drug and the RNA bases. These interactions indicate the importance and impact of RNA in the structure. The high affinity of IDR and FNT to the second active site and RNA gives us proof of the high potential anti-RdRP activity of these FDA drugs for the treatment of coronaviruses.

Conclusions

Due to the time-consuming development of antiviral drugs and the existence of high-speed computational approaches, utilizing in silico methods is a versatile approach to find reliable therapies for coronaviruses. In this study, the potential anti-RdRP activity of IDR and FNT, due to the tight binding of these drugs to the RdRP, revealed that these drugs could be considered as an appropriate treatment for COVID-19.

Declaration of Competing Interest

The authors declare that the investigators in this study have no conflict of interest.