Prediction of Small Molecule Inhibitors Targeting the Severe Acute Respiratory Syndrome Coronavirus-2 RNA-dependent RNA Polymerase

The current COVID-19 outbreak warrants the design and development of novel anti-COVID therapeutics. Using a combination of bioinformatics and computational tools, we modelled the 3D structure of the RdRp (RNA-dependent RNA polymerase) of SARS-CoV2 (severe acute respiratory syndrome coronavirus-2) and predicted its probable GTP binding pocket in the active site. GTP is crucial for the formation of the initiation complex during RNA replication. This site was computationally targeted using a number of small molecule inhibitors of the hepatitis C RNA polymerase reported previously. Further optimizations suggested a lead molecule that may prove fruitful in the development of potent inhibitors against the RdRp of SARS-CoV2.


■ INTRODUCTION
The ongoing COVID-19 (coronavirus disease 2019) global crisis has resulted in an unprecedented level of international turbulence, both in terms of global health and economy. The first ever patient diagnosed to have fallen from an infection caused by the then strain of coronavirus, 2019-nCoV, is traced back to the 1st of December 2019. 1 Within a span of over three months, the diseaseCOVID 19has been declared a global pandemic and has infected 619533 individuals, resulting in 376,320 deaths as of 2nd June 2020. 2 Of the seven strains of coronavirus known to infect humans, SARS-CoV (severe acute respiratory syndrome coronavirus), MERS-CoV (Middle East respiratory syndrome coronavirus), and the current strain SARS-CoV2 (severe acute respiratory syndrome coronavirus-2) are known to cause severe lower respiratory tract dysfunction. 3 The primary symptoms of COVID-19 are fever, cough, and discomfort in the upper chest accompanied with labored breathing, similar to those of SARS-CoV and MERS-CoV infections. Severe cases are reported to show sepsis, secondary infections, and organ failure. Reports suggest a possible combination of natural selection and zoonotic transfer events that led to the emergence of this novel strain. 4 The ∼29.9 kbp (+) ss-RNA genome of SARS-CoV2 shares 96.2% identity with the bat CoV RaTG13 genome and 79.5% identity with the SARS-CoV genome. 5 The amino acid sequence encoded by the genome shares 82.3 and 77.2% identities with the amino acid sequences encoded by the genomes of bat CoV SL-CoVZC45 and SARS-CoV, respec-tively. 6 The first polyprotein coding DNA sequence of the SARS-CoV2 genome (isolate Wuhan-Hu-1, NC_045512) 7 is predicted to encode for 14 nonstructural proteins, of which " nonstructural protein 12 " (YP_009725307.1) is composed of two domains, RdRp (RNA-dependent RNA polymerase) and NiRAN (nidovirus RdRp-associated nucleotidyl transferase). The RdRp is a key enzyme in RNA viruses which mediates the replication of their RNA genomes. 8 The RdRp from nidoviruses in addition to the canonical finger, palm, and thumb subdomains, possesses two extra domainsthe interface and the NiRAN domain. 9,10 The interface domain primarily acts as a linker between the RdRp and NiRAN domains and is also predicted to have interactions with the accessory proteins of the replicase complexthe Nsp7 and the Nsp8. 11,12 The NiRAN domain has been observed to possess GTP (guanosine triphosphate)-and UTP (uridine triphosphate)-binding activities that aid in the replication initiation process. 9 RdRp enzymes in general are known to exhibit both primer-dependent and primerindependent modes of replication initiation, while some are exclusive to either mode. 8,13−16 Although SARS-CoV RdRp exhibits a primer-dependent activity in vitro, 17 the RdRp from its sister speciesEAV (equine arteritis virus)demonstrates a primer-independent activity. 18 The RdRp enzymes from SARS-CoV and EAV share a similar structural conservation, suggesting that the nidovirus RdRp may exhibit both modes of initiation in vitro. 16 Moreover, the nsp7−nsp8 complex has been observed to have a primase activity, possibly aimed at a primer-dependent initiation. 19 Interestingly, the primer-independent activity of RdRp is attributed to the binding of single NTP (nucleotide triphosphate), preferably GTP within or in the vicinity of the active site pocket. 20−22 In this report, using computer-aided homology modeling, docking, and molecular simulations, we have predicted the protein structure and probable small-molecule inhibitors against SARS-CoV2 RdRp (CoV2-RdRp). ■ PREDICTION OF THE GTP BINDING SITE AND THE   INITIATION COMPLEX FOR COV2-RDRP A BLASTp search with the amino acid sequence for CoV2-RdRp against the PDB database yielded greater than 93% sequence identity with the sequence for RdRp from SARS-CoV (CoV1-RdRp). A homology modelling using the cryo-EM structure of CoV1-RdRp (PDB ID: 6NUR 23 ) as the template yields a globular structure of CoV2-RdRp with two distinct structural architectures. The overall structure is constituted of an RdRp-like structure (∼570 amino acids) fused with a NiRAN domain (∼360 amino acids) at the N-terminus. Similar to the RdRps from viruses such as dengue, 24 hepatitis C, 25−27 polio, 28 and others, the 3D structure of CoV2-RdRp exhibits the canonical RNA polymerase fold comprising the finger, palm, and thumb subdomains ( Figure 1A). A superimposition of the CoV2-RdRp and CoV1-RdRp structures ( Figure S1A) highlights the conservation of the secondary structural elements, which is also evident from the low root mean square deviation (rmsd) value of 0.93 Å. While this study was in progress, an experimentally determined (cryo-EM) structure of CoV2-RdRp (Nsp12) in complex with Nsp7 and Nsp8 was reported (PDB IDs6M71 and 7BTF 11 ). The experimentally determined structure and the model generated in this study superimpose with each other with an rmsd of 0.43 Å, suggesting that the model structure is nearly identical to the experimentally determined structure ( Figure 1B). In light of this, the results of our study would remain unaltered largely should it be reproduced using the EM structure.

■ RESULTS AND DISCUSSION
In order to derive insights into the structural similarities and differences between the structures of CoV2-RdRp and HCV-RdRp (RdRp from the hepatitis C virus), we superimposed them. The overlay resulted in an rmsd value of about 5 Å and showed an overall conservation of the RdRp architecture, with exceptions to a few secondary structural elements and flexible loops. Removal of the nonaligning regions resulted in better superimposition of the three RdRp subdomainsfingers, palm, and thumb ( Figure 1C), with an rmsd value of 3.2 Å. Similar to those in the HCV-RdRp polymerase, two extended structures [Δ1 (400−447 aa) and Δ2 (544−556 aa) motifs] protruding from the finger subdomain make interactions with the thumb subdomain. These two loop-helix-loop motifs bridge the finger and thumb subdomains. Similar to other known polymerases, the CoV2-RdRp active site is situated in the palm subdomain with two conserved aspartic acid residues, as evident from the pairwise sequence alignment of the superimposing sequences (Asp 760 and Asp 761 ) ( Figure 1D). Regardless of the variations in the detailed structures, their general architectures are largely similar. For further analysis, only the RdRp region of the model (amino acids 360−932) was chosen.
Previous reports suggest the binding of a GTP molecule adjacent to the first NTP near the active-site pocket for the de novo initiation of RNA replication by RdRp enzymes. 20,29,30 This GTP is added by a nucleotidyl transfer reaction. 20 Figure S1B). Further stabilization of the model GTP complex yielded the minimum energy conformations model for GTP binding at the predicted site (Figures 2A and S1C). The α-, β-, and γ-phosphate groups of the GTP molecule form H-bond interactions with Lys798, Pro620, and Arg553, respectively ( Figure S1D). The guanosine group makes potential H-bonding interactions with Asp761, Ala762, Glu811, and Ser814 residues of the enzyme. S3A,B), suggesting that GTP binding is an essential aspect of all RdRp enzymes with possible implications toward the initiation of the RNA replication process.
Furthermore, we attempted to generate a probable model of the initiation complex of the CoV2-RdRp. A comparison of the reported initiation models for the Pseudomonas phage φ6-RdRp (PDB ID: 1HI0 36 ) ( Figures 2B and S4A) and the mammalian orthoreovirus 3-Dearing-RdRp (PDB ID: 1N1H 33 ) ( Figure S4B) with the GTP-bound model of the CoV2-RdRp revealed that the GTP in CoV2-RdRp sits in the vicinity of the template RNA, similar to the GTP of Pseudomonas phage φ6-RdRp and the enterobacterial phage λ-RdRp. This suggests that the predicted GTP binding pocket of CoV2-RdRp possibly plays an important role in the initiation process. A recently published study demonstrates the CoV2-RdRp bound to the template and primer RNA molecules (PDB ID: 7BV2). 37 A comparison of the GTPbound CoV2-RdRp model shows interactions between the primer RNA and the bound GTP. This further suggests that the bound GTP may play a role in primer-independent replication or the GTP binding site may act as the receptor site for incoming nucleotides prior to its addition to the growing RNA chain. The superimposition of the GTP-bound model structure and the RNA-and remdesivir (monophosphate)bound cryoEM structure of CoV2-RdRp 37 is presented in Figure 2C. The GTP in the model structure and the remdesivir molecule in the cryo-EM structure occupy overlapping regions in the active-site region of the CoV2-RdRp; however, the GTP is stabilized via interactions with amino acid residues, remdesivir shows interactions with amino acid residues and the nucleotides from the bound RNA molecule ( Figure S5). Taking together the aforementioned interaction and comparison of the model and experimentally determined structures, we propose the probable initiation complex of CoV2-RdRp bound to RNA and GTP molecules in Figure 2D.
Exploration of Inhibitors Targeting the GTP Binding Site of CoV2-RdRp. Over the past several years, a targetbased approach has led to the design of a number of small molecule inhibitors against HCV-RdRp. 26,27,38−43 Importantly, the current drug (sofosbuvir) 44 for HCV treatment targets the HCV-RdRp. The binding of GTP is a key step in the RNA replication process, suggesting the therapeutic importance of the GTP binding site pocket. 20,29 Previous reports 26,27 involving studies with thiophene scaffold inhibitors against HCV-RdRp suggest that these inhibitors target the thumb site II, nearly 35 Å from the polymerase active site. The inhibitor binding to this site may interfere with the GTP binding function of the RdRp. 45 We sought to explore whether these inhibitors and the current anti-HCV drug (sofosbuvir) and a broad spectrum antiviral drug (remdesivir) exhibit any binding affinity with the CoV2 polymerase, specifically targeting the GTP binding site. The 2D structures of six small molecules are presented in Figure 3A−F (left panels). Directed docking of a few proposed HCV-RdRp inhibitors along with the triphosphate/metabolite forms of sofosbuvir and remdesivir (Figure 3, right panels) exhibited binding interactions within the GTP binding pocket ( Figure S6). The docking properties for these candidate molecules are listed in Table 1 (rows 1−7). The detailed view of the molecular interactions within the GTP binding site for the three molecules with the most favorable docking properties is presented in Figure S7. A comparison of the remdesivir (triphosphate)-bound model structure and the remdesivir (monophosphate)-bound cryo-EM structure of the RdRp revealed that the molecules occupy similar locations within the active-site region of the molecule ( Figure S8A). A comparison of the atomic interactions of these molecules with  the amino acids of the RdRp shows common interactions with Lys545, Arg555, Val557, and Cys622 ( Figure S8B). Based on the molecular interactions and docking scores, we attempted to optimize the top three candidate molecules in order to achieve a stronger binding affinity. On computationally mutating the small molecule JPC ("3-[(2,4-dichlorobenzoyl) (isopropyl)amino]-5-phenylthiophene-2-carboxylic acid") to "3-[(2-bromo,4-chlorobenzoyl)(isopropyl)amino]-5phenylthiophene-2-carboxylic acid" (2D structure presented in Figure 4A, left panel), a significant rise in the docking affinity was observed ( Figure 4A, right panel). The interactions and the docking properties for this lead optimized molecule are presented in Figure 4B and Table 1 (Row 8), respectively. The detailed view of the interactions for this molecule is presented in Figure 4C.
The ADME/T (absorption, distribution, metabolism, elimination, or toxicity) property values of all the eight molecules satisfied the acceptable ranges in terms of Lipinski's rule of five 46 and other pharmacokinetic parameters such as probable H-bonding atoms, human oral absorption (<25%: poor; >above 80%: good), and the IC 50 values for blockage of hERG (human ether-a-go-go-related gene) K+ channels (<−5: satisfactory). Interestingly, the lead optimized molecule is predicted to have 94.7% of oral absorption and its hERG value evidenced the degree of drug likeliness. The detailed results of this analysis are listed in Table 2.
Molecular Dynamics simulation studies of the native and ligand-bound complexes of CoV2-RdRp. MD (Molecular dynamics) simulations were performed for the modelled structure of the RdRp protein and docked complexes for the GTP, lead optimized, and FIH compounds for a 50 ns time period. Here, the GTP molecule was considered for the simulation studies because it plays a vital role in the RdRp functional mechanism; thereby, the structural stability and potential interactions of the GTP in the active-site pocket of the RdRp protein can be understood. Additionally, the complexes of the top scored ligand molecules (lead optimized and FIH) were considered for understanding the structural  stability and potential interacting residues for the inhibitory mechanism of the complexes. In MD simulation studies, rmsd and RMSF (root mean square fluctuation) are the indicators for analyzing the stability of the complexes. The rmsd, RMSF, SASA (solvent accessible surface area), and H-bonding analysis results are presented in Figure 5 and Table 3. A comparison of the simulation statistics for the native and GTP-bound RdRp suggests that while the GTP-bound complex initially undergoes higher deviations and fluctuation, the complex gets stabilized over the course of the simulation. Additionally, there is a significant rise in the accessible surface area of the GTPbound complex, suggesting that the RdRp assumes an open conformation following the binding of GTP, possible assuming a preinitiation state in order to accommodate the incoming template and/or the primer RNA molecule. Interestingly, all the phosphate groups of the GTP molecule displayed conformational changes during the simulation and the first 50 residues at the N-terminal from the palm subdomain were found relaxed. Also, in the first 5 ns, the GTP exhibited approximately 15 H-bond interactions within the binding site, which got stabilized to 8 H-bond interactions throughout the simulation period. The complexes of RdRp with both FIH and the lead optimized compound exhibited significantly lesser deviations and fluctuations along with a relative decrease in the solvent accessibility as compared to the GTP complex. This suggests stronger competitive binding of the inhibitor molecules along with the generation of a possible half-closed RdRp conformation, indicating the inhibitory potential of these molecules. Interestingly, the FIH complex showed relatively lesser fluctuations and deviations as those of the lead molecule complex, suggesting that in a physiological state, the FIH molecule may exhibit stronger inhibition. However, the lead molecule exhibited higher H-bond interaction as compared to FIH. Of note, in the rmsd plot of the lead moleculeRdRp complex, an abrupt deviation is observed between the 15 and 25 ns time period and the remaining periods are highly stable. During this sudden deviation, the N-terminal region (50 residues) was found relaxed from the palm subdomain at 15 ns, after which the region again tends toward to the subdomain and stabilizes the complex during the rest of the cycle. These suggest significant conformational changes following the binding of the molecule, possibly causing the RdRp molecule to assume a closed state.
Density Functional Theory Studies of the FIH and Lead Optimized Compound. The DFT (density functional theory) analysis investigates the electronic features of the atoms in a charge transfer reaction, using the HOMO (highest occupied molecular orbital), LUMO (lowest unoccupied molecular orbital), and MESP (molecular electrostatic surface potential) analyses. This provides key information about the global and local indices of the compounds that help to understand its chemical reactivity. The electron donor/ acceptor properties of the molecules are indicated by the distribution of frontier molecular orbitals. The DFT study was performed for the top two docking scored molecules (FIH and lead optimized compound). These two molecules show low HOMO−LUMO energy gaps, signifying higher stability and chemical reactivity of the molecules. The calculated stereoelectronic properties and profiles of these compounds are presented in Table 4 and the HOMO, LUMO, and MESP contour maps are shown in Figure S9. In FIH, the HOMO regions are localized in the thiophene-2-carboxaldehyde group (energy value −0.22 eV), while the LUMO regions are localized at the sulfonyl amino group (energy value −0.02 eV). Together, the higher stability of the energy value of the HOMO as compared to the LUMO suggests that the compound may favor a nucleophilic reaction. The calculated HOMO−LUMO energy of the gap of the compound is −0.20 eV. The electrostatic charge surface potential analysis revealed that highly negatively charged surfaces are present around the molecule, while positively charged surfaces are found in the sulfonyl and −OH group. For the lead molecule, the HOMO regions are present in the sulfur and hydrogen atom of thiophene group and the LUMO regions localize in the hydroxyl group of the benzoyl ring. The calculated HOMO and LUMO energies of the compound are −0.24 and −0.002 eV, respectively, suggesting that the compound may favor the nucleophilic reaction. The charged surface density of the lead molecule is extremely low with only the negatively charged surface observed in the molecule. The calculated energy gap of the lead molecule (−0.22 eV) is higher than that of the FIH molecule (−0.20 eV), suggesting that the FIH molecule may demonstrate a higher inhibitory reaction potential under physiological conditions.

■ CONCLUSIONS
The urgency arising from the current COVID-19 pandemic not only entails major research initiatives on the development of novel antiviral molecules but also drives the repurposing of the available drugs and inhibitors reported in previous studies, particularly on (+) ss-RNA viruses. This study reports the probable GTP binding site in the active site region between the finger and palm subdomains of the CoV2-RdRp. Comparison with the previously reported RdRp structures from RNA viruses shows the conservation of this GTP binding pocket and moreover suggests its close proximity to the RNA template and the RNA primer. The RdRp enzymes are known to function in both primer-dependent and -independent manner. This  possibly indicates the dual function of the predicted GTP binding site. In the case of a de novo initiation of the RNA replication process, this site could possibly incorporate the first and the second nucleotide, GTP being the second nucleotide, as experimentally demonstrated in a previous study. 20 Alternatively, in the case of a primer-assisted initiation process, this site can potentially function as a nucleotide-binding site for an incoming nucleotide prior to its addition to the growing RNA chain, as evident from the proximity of the GTP to the primer RNA. Thus, theoretically, the inhibition of this site by small molecules could effectively dampen both modes of initiation, thus rendering the RdRp enzyme dysfunctional. In addition, this study suggests that the HCV-RdRp inhibitors reported in previous studies are likely to bind the CoV2-RdRp enzyme at the predicted GTP binding site. Based on the predicted interactions, this study also reports a lead optimized molecule exhibiting a higher binding affinity and drug likeliness. Biochemical experiments need to be carried out to examine the in vitro inhibitory potential of these molecules against CoV2-RdRp.

■ COMPUTATIONAL METHODS
Protein Structure Prediction. The primary amino acid sequence of CoV2-RdRp was retrieved from the NCBI database (accession Y_009725307). A BLASTp search was performed against the PDB database to identify the structural homologues of CoV2-RdRp. The protein structure for RdRp-CoV1 (PDB ID: 6NUR) was considered as the template for homology modelling to building the CoV2-RdRp structure using the I-TASSER server. 47 The structure generated was used for further studies.
In Silico Protein and Ligand Preparation. The modeled structure of CoV2-RdRp was further improved using protein preparation wizard (OPLS-2015) incorporated in "Schrodinger Suite 2019" 48 to improve the quality of the structure by adding hydrogen atoms and hydrogen bonds and removal of the steric clashes. This minimized model of the CoV2-RdRp structure was used in further analysis. The four potential inhibitors of HCV-RdRp along with the known inhibitors of antiviral drugs remdesivir and sofosbuvir were selected for interaction studies with CoV2-RdRp. The HCV-RdRp inhibitors (JPC, CCT, SNH, IPC, and FIH) were retrieved from the reported crystal structures (PDBs ID: 1YVZ, 1YVX, 2D41, 1D3Z), while structures for the triphosphate forms of remdesivir (ID: DBMET03038) and sofosbuvir (ID: GS-461203) were retrieved from DrugBank. All the inhibitors were prepared by assigning appropriate geometries, bond orders, tautomers, and ionization states using the LigPrep module implemented in "Schrodinger Suite 2019".
GTP Binding Site Prediction and Docking Analysis. The minimized model was used to predict the GTP binding site using the SiteMap module in the "Schrodinger Suite 2019". In order to verify the accuracy of the site, the GTP binding site was regenerated using the Phyre server, 49 CASTp server, 50 and the results were compared. The results predicted the nearidentical sites in the core region of CoV2-RdRp. This site was then used to generate a grid in order to dock the inhibitors. The prepared inhibitors were docked within the GTP binding site pocket of the CoV2-RdRp protein using the Glide XP (extra Precision) module implemented in "Schrodinger Suite 2019".
Molecular Mechanics Generalized Born Surface Area Calculation. The binding free energies of the enzyme− inhibitor complexes were analyzed by molecular mechanics generalized Born surface area calculation. 51 The calculation was performed on the basis of the following equation where ΔE MM represents the minimized molecular energy changes in the gas phase, G SGB represents the surface calculation using the GB model, and G SA represents the accessible surface area of the complexes. Lead Optimization. Lead optimization is an important method in the drug development process to improve the potency and pharmacokinetic properties of the drug molecule. 52 Here, lead optimization was performed against the JPC compound using "Schrodinger Suite 2019". The JPC compound has two chlorine atoms in its benzoyl ring. In order to enhance the potency of the molecule, we modified the molecule with different combinations of fluorine (F), bromine (Br), and chlorine (Cl) atoms at the second and fourth position of its benzoyl ring. The modified lead molecules were analyzed for their docking and ADME/T properties. ADME/T Property Prediction. The seven inhibitors were analyzed for their drug likeliness and pharmacokinetic properties such as ADME/T. 53 The QikProp module implemented in the "Schrodinger Suite 2019" was used to estimate the drug likeliness of the inhibitors using Lipinski's rule of five 46 and other pharmacokinetic parameters such as probable H-bonding atoms, human oral absorption (<25%: poor; >above 80%: good) and the IC 50 values for blockage of human ether-a-go-go-related gene (hERG) K + channels (<−5: satisfactory). 54 MD Simulation Studies. The structural stability of native and inhibitor-bound complexes of SARS-CoV2-RdRp was analyzed using the GROMACS (Groningen Machine for Chemical Simulation) 55 software for a 50 ns time period. The Gromos56a force field with its default parameters was applied to the protein and inhibitor complex(es). The inhibitor topology was refined and generated using the PRODRG server 56 and incorporated into the structure for understanding the dynamics and behavior of the complex. A cubic box of 10 Å dimension was set with SPC216 (simple point charge) water molecules and the periodic boundary condition was applied in the all directions. Further, the system was neutralized by including 12 Na 2+ ions in the native and complex system(s). Energy minimization and equilibration are important steps in the generation of the MD. The steepest decent (SD) algorithm was performed for 50,000 steps with a tolerance of 1000 kJ/ mol to minimize the system. The thermostat and barostat algorithms such as NVT [canonical ensemble-substance (N), volume (V), and temperature (T)] and NPT [isothermal isobaric-substance (N), pressure (P), and temperature (T)] ensembles were used. The linear constraint solver (LINCS) and particle mesh Ewald (PME) method was used to constraint the covalent bonds and long-range electrostatic interactions of the systems, respectively. The temperature (310 K) and the pressure (1 atm) of the systems were equilibrated by the V-rescale coupling and Parrinello−Rahman (PR) method. The preprocessed systems were used to perform the MD simulations for a 50 ns time period.
DFT Calculation. In quantum mechanical calculation, DFT analysis determines the molecular electronic features such as electron density, frontier molecular orbital (HOMO and LUMO) density, and molecular electrostatic map to predict the chemical reactivity of proposed inhibitors. 57,58 The DFT study was carried out for the top two predicted inhibitors in order to understand the inhibitory reaction against the SARS-CoV2-RdRp. Geometry optimization of the compound was performed using a hybrid DFT approach at B3LYP (Becke's three-parameter exchange potential and the Lee−Yang−Parr correlation functional) with the 6-31G** basis set. 59,60 The Poisson−Boltzmann solver was used to calculate the energy under aqueous conditions, simulating physiological conditions, which provides the information about the global and local indices of the molecules to their chemical activity. The spatial distributions of electronic features in the charge transfer mechanism were understood from the HOMO and LUMO molecular orbitals. All DFT calculations were carried out using Jaguar, version 9.1. 61,62 ■ ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsomega.0c02096.