Mechanistic insight on the remdesivir binding to RNA-Dependent RNA polymerase (RdRp) of SARS-cov-2.

The RNA-dependent RNA polymerase (RdRp) is a key enzyme which regulates the viral replication of SARS-CoV-2. Remdesivir (RDV) is clinically used drug which targets RdRp, however its mechanism of action remains elusive. This study aims to find out the binding dynamics of active Remdesivir-triphosphate (RDV-TP) to RdRp by means of molecular dynamics (MD) simulation. We built a homology model of RdRp along with RNA and manganese ion using RdRp hepatitis C virus and recent SARS-CoV-2 structures. We determined that the model was stable during the 500 ns MD simulations. We then employed the model to study the binding of RDV-TP to RdRp during three independent 500 ns MD simulations. It was revealed that the interactions of protein and template-primer RNA were dominated by salt bridge interactions with phosphate groups of RNA, while interactions with base pairs of template-primer RNA were minimal. The binding of RDV-TP showed that the position of phosphate groups was at the entry of the NTP channel and it was stabilized by the interactions with K551, R553, and K621, while the adenosine group on RDV-TP was pairing with U2 of the template strand. The manganese ion was located close to D618, D760, and D761, and helps in stabilization of the phosphate groups of RDV-TP. Further we identified three hits from the natural product database that pose similar to RDV-TP while having lower binding energies than that of RDV-TP, and that SN00359915 had binding free energy about three times lower than that of RDV-TP.

Introduction

Coronavirus 2019 (COVID-19) is a disease due to the SARS-CoV-2 virus which presently threatens human life across the globe. It is highly contagious and has affected the lives of millions of people with over one million deaths as of October 2020. There is therefore an urgent need to discover effective therapeutics for eliminating the disease considering its detrimental effects to human life.

SARS-CoV-2 is a positive-strand RNA virus whose genomic structure expresses high similarity with previous severe acute respiratory syndrome coronavirus (SARS-CoV). Like other coronaviruses, the proliferation of the SARS-CoV-2 requires non-structural protein 12 (nsp12) (RNA-dependent RNA polymerase, RdRp), a key enzyme which regulates the viral genome replication and transcription. It is thus the validated target for the development of COVID-19 disease therapeutics. Remdesivir (RDV) is a clinically-used drug which targets RdRp and its efficacy against Ebolavirus (EBOV) infection is well-understood [1,2], SARS-CoV-1 and Middle East respiratory syndrome coronavirus (MERS) [[3], [4], [5], [6]] has also been previously reported. Due to the high sequence similarity of RdRp structures of SARS-CoV-1 and SARS-CoV-2 (96% sequence identity), it has been therefore approved for emergency use of treating SARS-CoV-2 virus [7].

Remdesivir is a prodrug and nucleoside analogue which functions after biotransformation into its active triphosphate form (RDV-TP). It works as a substrate mimic of ATP, which is a substrate for RdRp for nucleotide addition on the primer RNA strand. However, RDV has been proposed to work differently during chain translocation. Several works proposed that the RDV induced delayed chain termination [8,9], however, its mechanism of action is not completely understood [10,11]. This study aims to elucidate the mechanism of action and dynamics of the RDV binding into RdRp by means of molecular dynamics (MD) simulation. The MD simulation approach has been widely accepted to get invaluable insight into drug binding in a specific target at the molecular level.

Previous computational studies on the interaction between Remdesivir and SARS-CoV-2 RdRp were reported. For instance, studies by Elfiky (2020a, 2020b) [12,13] indicated the strong binding of Remdesivir to RdRp. Kato et al. (2020) described intermolecular interactions between Remdesivir and SARS-CoV-2 RdRp by employing fragment molecular orbital calculations, which indicated that T687, N691, and D760 are critical residues [14]. Meanwhile, a study by Zhang and Zhou (2020) described the inhibitory mechanism of Remdesivir to SARS-CoV-2 RdRp [15]. However, the previous studies lacked of the RNA template-primer in their RdRp complex models, which is crucial in the nucleotide triphosphate binding. Here we reported the Remdesivir interactions with the RdRp complex along with RNA template-primer.

We firstly built the RdRp complex along with the primer RNA, template RNA, and manganese ions prior to docking RDV-TP and ATP to the RdRp model. We then employed the RDV-TP bound to the RdRp model to identify molecular hits of potential RdRp inhibitors from the natural product database through pharmacophore screening. Natural products are continuously showing their importance in drug discovery and development as up to 80% of drugs including those that inhibit viral replication approved by the Food and Drug Administration between 1981 and 2014 are linked to natural products [16]. For instance, catechin and epigallocatechin gallate (EGCG) were reported to have antiviral activity against SARS-CoV and MERS-CoV, respectively [17]. We validated the binding of our natural product hits using 200 ns MD simulation for each compound and performed MM-PBSA energy calculations to evaluate their affinities to RdRp. Interestingly, the three molecular hits had lower binding affinities than that of RDV-TP.

Experimental section

Homology modeling and ligand preparation

The RdRp model was constructed using Maestro's homology modeling tools [18]. The electron microscopy solved structure of SARS-CoV-2 non-structural protein 12 (NSP 12, PDB ID: 6M71) [19] and X-ray solved structure of RNA primer template hepatitis C virus NS5B (4WTG) [20] was employed. The charge state of the RdRp complex was optimized at pH = 7. The complex was further preprocessed, optimized, and minimized with Maestro's protein preparation wizard tool [18].

Using the prepared RdRp structure, a receptor grid file was generated around the RDV-TP binding site and the ionization/tautomeric state of RDV-TP was generated at pH = 7 using Maestro's Epik tool [21,22]. The ligand was then docked to the receptor using the Glide XP docking [23,24] scoring function. The MD simulation was then performed to check the binding dynamics of RDV-TP.

For comparison and validation of both our constructed RdRp complex and modeling methods for accuracy and consistency, we utilized the recently deposited crystal complex of RNA polymerase nsp12-nsp7-nsp8 with bound template-primer RNA and RDV-TP (PDB ID: 7BV2) [25]. After docking, we performed a sequence alignment of our newly bound polymerase-RDV-TP complex onto 7BV2, including the RNA template-primer and RDV-TP models and corresponding metal ions for best comparison. Figure S1 shows that our model is consistent with the geometry of 7BV2 in regards to partial protein conformations, RNA template-primer shape and RDV-TP binding poses. This information encouraged us to further continue with this model in the future simulations.

Pharmacophore screening was performed using active Remdesivir (Figure S2) against 325,187 natural product molecules [26] by means of Ligand Scout 4.3 [27]. In our pharmacophore model, we defined seven pharmacophore features including two negative ionizable, one aromatic ring and four H-bond acceptors and their corresponding positions. The natural origin compounds were our greatest priority in the present study since they exhibited inhibition to the RdRp of other RNA viruses, high bioavailability, low toxicity, and received greater importance for their antiviral applications [28,29]. Before docking, the ionization/tautomeric states of each compound were generated at pH = 7 using Maestro's Epik tool [21,22]. The lowest ionization/tautomeric states were selected, and then their geometries were minimized to their most energetically favorable states. To check if our docked ligand poses were reasonable, we compared the docked poses of the 65 ligands with the RDV-TP pose. Clearly, the 65 ligands bind similarly to the binding pocket of RDV-TP but with some differences. This provided a reasonable starting pose for future MD simulations which could further refine the binding pose given the full conformation flexibility from the simulations.

MD system setup

Several MD simulations were performed: one 500 ns simulation for the apo-form of RdRp, three 500 ns simulations for the RdRp-RDV-TP complex, three 500 ns simulations for the ATP-RdRp complex, and a 200 ns simulation for each complex of SN00166900, SN00303170, and SN00359915. MD simulations were conducted by using AMBER16 software [30]. The AMBER ff14SB [31], RNA.OL3 [32,33], and GAFF2 [34] forcefields were applied for protein, RNA and ligands, respectively. A salt concentration of 0.15 M was applied to each system and were placed in a truncated octahedron TIP3P water box with a 10 Å radius. Each system underwent a 10000 step minimization three times, which consisted of 500 steps of steepest descents and 9500 conjugate gradients each. Then, minimizations were applied to each system to relax ion and water molecular positioning, with second and third minimizations being conducted with and without protein backbone atoms restraint, respectively. Each system was heated from 0 to 100, 100–200, and 200–300K, with each undergoing 50 ps under NVT ensemble with a time step of 0.0005 ps and a 5 kcal mol−1 A−2 restraint. Each system was then allowed three steps of relaxation under NPT ensemble, each with and without 3 kcal mol−1 A−2, 1 kcal mol−1 A−2 restraints, respectively.

The Langevin thermostat with a collision rate of 1.0 ps−1 was applied to all covalent bonds involving hydrogen atoms and were restrained using SHAKE algorithm [35]. The PME method was employed to treat long-range electrostatic interactions with an integration step of 2 fs [36]. The long-range non-bonded interactions were calculated with a cutoff distance of 9.0 Å by applying periodic boundary conditions.

To check the convergence of MD simulations, we calculated the RMSD values of protein Cα, RNA backbone, and ligand main atoms for each trajectory by employing the CPPTRAJ module [37]. In addition, the root mean square fluctuations (RMSF) of each individual residue of the protein, RNA and individual atoms of the ligands were calculated to characterize the local structural fluctuations in each complex.

Clustering analysis was employed to group similar RdRp-ligand structures with ions and solvent molecules being omitted. The backbone RMSD matrix was used as the metric for similarity, and the clustering method was accessed by DBSCAN [38] via CPPTRAJ module. The distance cutoff (epsilon) was 2.5 and minpoints 15, and the most abundant clusters were used for analysis. For the RDV-TP complex, the last 200 ns of simulation data were used, while for complex of SN00166900, SN00303170, and SN00359915, the last 100 ns of simulation data were used for clustering analysis.

The binding energy of ligands were calculated using the Molecular Mechanics-Poisson Boltzmann solvent accessible surface area (MM-PBSA) method [39] via the MMPBSA.py module. Trajectories of 200 snapshots of the last 100 ns of each complex was used for calculations.

Results and discussion

Homology model of RdRp complex

Before constructing the RdRp complex model, sequence alignment between SARS-CoV-2 non-structural protein 12 (NSP 12, PDB ID: 6M71) and RNA primer template hepatitis C virus NS5B (4WTG) was conducted. Superposition of the modelled and crystallographic structures is depicted (Figure S3-A), which shows that both structures align well with each other, particularly for amino acid residues 1–541. The interactions of protein and template-primer RNA mostly originated from contributions in the RNA phosphate group. Hydrogen-bonding interactions occurred with S501, R569, S592, Y595, and Y689, while salt bridge interactions occurred with K500, K577, and R583 (Figure S3–B). The template-primer RNA base pair interactions with the protein occurred between Uracil and S682 (H-bond) and K500 (pi-cation) (Figure S3, B and C). This indicated that the elongation step of RdRp function is not sequence specific which was also implied in the previous study [25].

During the MD simulations, we monitored the H-bond occupancy between RNA and protein during the 500 ns timescale (Table S2). It showed that the high H-bond occupancies occurred between U5 (U554) primer strand and U2 (U543) of the template strand with S814 (S445) and S682 (S313) occupancies of 96.87% and 82.53%, respectively. The modest occupancies occurred in many H-bonds interactions; for instance, those between A3 (A544) of the template strand with D684 (D315) (69.16%), A4 (A545) with R569 (R200) (66.03%), and U8 (U549) of template strand with Y595 (Y226) (58.40%). In addition, various H-bond interactions between protein and RNA and between RNA bases were observed to have low occupancies (Table S2).

The binding of remdesivir and ATP to RdRp complex

RDV-TP was docked at the i position against RdRp SARS-CoV-2 model which contains 6 bases in the primer strand, 8 bases in the template strand, and a manganese ion (Fig. 1 A). The conformation of RDV-TP was identical with the ATP pose (Fig. 1B), with both RDV-TP and ATP were located at the catalytic site of RdRp. The position of the phosphate group at the NTP entrance channel may reflect its role in inhibiting the entry of NTP to the active site.

Fig. 1

The 3D poses of RDV-TP (A) and ATP (B) in the RdRp SARS-CoV-2 model. Primer and template strands are colored red and cyan colors, respectively.

The detailed interactions of docked RDV-TP and ATP with RdRp are depicted (Fig. 2 ). Both phosphate groups of RDV-TP and ATP formed H-bond interactions with positively charged amino acid residues K621 and R555, while salt bridge interactions formed with K551 and R555. The phosphate groups of both ligands were additionally stabilized by manganese ions through metal coordination and salt bridge interactions. The adenosine motifs of both RDV-TP and ATP established H-bond interactions with U2 of template strands and pi-pi stacking interactions with U6 of primer strands. Additional pi-pi stacking interaction were formed with A3 of the template strand in the ATP conformation. Interactions with N691 were formed in both RDV-TP and ATP through their 3′-OH groups. In both complexes, D623 and D760 were located close to phosphate and manganese ions, while the nitrile group of RDV-TP was surrounded by T687 and A688. Clearly, the docked poses of both RDV-TP and ATP were almost identical.

Fig. 2

The detailed interaction of docked RDV-TP (A) and ATP (B) with RdRp.

During the 500 ns simulations, several H-bond interactions were conserved. The RDV-TP H-bond interactions occurring with U2 (U543) of template strand were retained at high occupancy (75.39%), while that occurring with D760 (D391) was maintained at modest occupancy (48.5%). In contrast, the H-bond with R555 (R186) was maintained at low occupancy (23.62%). Numerous H-bond interactions between RDV-TP and RdRp complex occurred with low occupancies (Table S3). In the case of ATP, the H-bond occupancies observed in the docked conformation were low such as those occurring between ATP and R555 (R186) (33.87%) and K551 (K182) (12.91%). Clearly, the RDV-TP was more capable at maintaining base pairing particularly with U2 (U543) of template strand. RDV-TP had more intense H-bond interactions during the MD simulation as compared with its ATP counterpart. The H-bond occupancies of ATP complex during 500 ns MD simulation are tabulated in Table S4.

The trajectory convergence during MD simulation was checked through the RMSD values averaged over three independent 500 ns MD simulation for backbone receptor Cα, RNA and ligand main atoms for apo-RdRp, ATP, and RDV-TP complexes (Fig. 3 ). It is clear in Fig. 3 that the receptor backbone quickly reached stability before 100 ns and remained stable throughout the simulation time for both RDV-TP and ATP systems. The RMSD of RNA main atoms showed greater fluctuations and reached stability at ~280 ns in both RDV-TP and ATP, while the average fluctuation during the whole simulation time was ~2 Å. The ligand RMSD of RDV-TP reached stability at ~50 ns, while those of the ATP system occurred at an earlier simulation time (~100 ns). The RMSD values of both RDV-TP and ATP systems were higher than that of the apo system, which indicated that the binding of RDV-TP and ATP induced greater fluctuations in RdRp complex during 500 ns simulation time. The RMSD values for the first, second and third 500 ns MD simulations of ATP and RDV-TP show similar patterns (Figure S4).

Fig. 3

RMSD values for the receptor Cα, RNA main atoms and ligand main atoms averaged over three 500 ns runs for apo RdRp and RDV-TP (A) and apo RdRp and ATP (B). Apo, RDV-TP, and ATP were colored red, green, and cyan colors, respectively.

To observe the fluctuations of amino acid residues, RNA residues, ATP and RDV-TP atoms during MD simulation, the root mean square fluctuation (RMSF) was recorded (Fig. 4 ). As shown in Fig. 4, the high peaks were observed at G16 (G385), V36 (V405), G63 (G432), N78 (N447), N127 (N496), W140 (W509), S149 (S518), G215 (G584), T274 (T644), R366 (R735), A428 (797), G454 (G823), T481 (T850), which corresponded to loop regions of the protein, while L526 (L895) and E541 (E919) corresponded to the C-terminals of the protein. Residues interacting with RDV-TP such as K551 (K182), R555 (R186), D623 (D254), N691 (N322), and D760 (D391) were observed to be stable.

Fig. 4

RMSF plot of protein, in which apo RdRp, RDV-TP, and ATP complexes are colored red, green, and cyan colors, respectively.

Meanwhile, RMSF of the RNA main atoms (Fig. 5 ) was highest at A1 (~2.5 Å) which is the first sequence of primer strand and did not interact directly with both RDV-TP and ATP. It is worth noting that the RMSF of U2 and A3 of template strands as well as U6 of primer stands in both RDV-TP and ATP systems were observed to be stable around 1 Å.

Fig. 5

The RMSF plot of RNA main atoms. Apo RdRp, RDV-TP, and ATP complexes are colored red, green, and cyan colors, respectively.

The RMSF of ligand atoms of both RDV and ATP system are depicted (Figure S5). The high peaks were recorded at O9, O10, and O11 (i.e. the oxygen atoms at γ phosphate group of RDV-TP and ATP). Atoms of the adenosine motif of RDV-TP fluctuated less than those of ATP system, which may indicate that base pairing of RDV-TP occurred with more intense interaction compared with the ATP system. Figures S6, S7, and S8 show the RMSF values for each three independent 500 ns MD simulation.

We then performed clustering analysis to generate structural representatives of the most abundant clusters for the 500 ns simulations. The most abundant clusters of RDV-TP and ATP complexes had 100% occupancy. The adenosine motif of RDV-TP interacted with U2 (U543) of template strand via H-bonding, which is absent in ATP interactions. The phosphate groups of RDV-TP and ATP interacted with positively charged residues such as K182 (K551), R184 (R553), and K252 (K621). The manganese ions were also located in the vicinity of the phosphate groups of RDV-TP and ATP, which indicated its role in ligand stabilization. Additionally, the manganese ions were also in close proximity with aspartate residues such as D249 (D618), D391 (D760), and D392 (D761). The nitrile group of RDV-TP was observed to be located in the empty space between C253 (C622), P251 (P620), and V188 (V557).

Our study is partly consistent with several previous studies. For instance, Hillen et al. (2020) implied that N691, S682 and D623 may recognize the 2′-OH group of NTP, thereby rendering the RdRp specific for the synthesis of RNA rather than DNA [40]. Yin et al. (2020) indicated that K545, R555, S682, N691, and D760 supported Remdesivir binding [25], while Zhang and Zhou (2020) implied that D618, S549, and R555 are major contributors to Remdesivir binding [15].

In addition, previous reports by Nudler (2009) and Burton et al. (2005) indicated that the substrate was pairing with the RNA template after entering through the main channel of the RdRp complex [41,42]. Its structural resemblance to natural ATP enables it to compete with ATP during viral synthesis [43]. We confirmed the interaction of U6 of primer strand at i site and RDV-TP at i+1 site and hence we speculated that the interaction would lead to delayed chain termination in i+3 or i+5 sites. Study by Gordon et al. (2020) and Tchesnokov et al. (2019) implied that the inhibition of the viral nucleic acid synthesis occurred through the delayed chain termination mechanism, in which following incorporation of RDV-TP, several more NTP molecules could still be added to RNA [6,8]. In brief, our study provides the details at an atomistic level on how Remdesivir behaves during its interaction with RdRp, which may be useful for finding more potent compounds for treating COVID-19 through RdRp targeting.

We then applied pharmacophoric screening against the natural database using pharmacophore features. Screening against natural product database using seven pharmacophoric features, i.e. two negative ionizable, one aromatic ring, and four H-bond acceptors (Figure S2), retrieved 56 molecular hits. We proceeded to select the three best hits for a 200 ns molecular dynamics simulation for each. The three best docked hits are subjected for 200 ns MD simulation, i.e. SN00166900, SN00303170, and SN00359915.

The best docked hit molecules show the similar poses with RDV-TP. The oxane group of SN00166900 was seen to sit after U6 of primer strand and its phenyl ring established pi-pi interactions with U6. More occurrences with N691 and R555 were seen through H-bond and salt bridge interactions. Additionally, the manganese ion interacted with OH group of ligand through salt bridge interactions. Fig. 6 displays the interaction of SN00166900 with RdRp.

Fig. 6

The interaction of SN00166900 with RdRp complex. Primer and template strands are colored red and cyan colors, respectively.

Meanwhile, the SN00303170 pose was also located in the U6 of primer strand and the pi-pi interaction was recorded between its phenyl ring and U6. The residues R555, N691, and manganese ion were also observed to interact with SN00303170 through H-bond and salt bridge interactions. Fig. 7 depicts the conformation of SN00303170 in the RdRp complex.

Fig. 7

The interaction of SN00303170 with RdRp complex. Primer and template strands are colored red and cyan colors, respectively.

The docked pose of SN00359915 was similar to other ligands, in which it occupied the space above the U6 primer strand. The pi-pi stacking interaction in SN00359915 binding occurred at A3 and A4 of the template strand. The H-bond interactions were recorded between the ligand with R555, N691, and C622. The manganese ion also formed a salt bridge interaction with SN00359915. Fig. 8 depicts the conformation of SN00359915 in the RdRp complex.

Fig. 8

The interaction of SN00359915 with RdRp complex. Primer and template strands are colored red and cyan colors, respectively.

MD simulation of molecular hits

To get the binding dynamics of each ligand, RMSD values for 200 ns were monitored. Fig. 9 shows the RMSD values for the receptor Cα, RNA and ligand main atoms of each complex for the 200 ns simulation.

Fig. 9

The RMSD values of RDV-TP and SN00166900 (A), RDV-TP and SN00303170 (B), RDV-TP and SN00359915 (C), and RMSF (D) values for each complex. The RDV-TP, SN00166900, SN00303170, and SN00359915 are colored green, cyan, purple, and pink colors, respectively.

The RMSD of protein Cα, RNA and ligand main atoms of each complex was stable during 200 ns. The stability of protein Cα was reached rapidly in the first 25 ns and for the rest of the simulation time for each SN00166900, SN00303170, and SN00359915 system. It was observed that the RMSD values of protein Cα of RDV-TP complex was lower than those of SN00166900, SN00303170, and SN00359915 complexes. Meanwhile, the RMSD values of RNA main atoms were slightly lower in SN00166900 and SN00359915 systems and was higher in SN00303170 than those of the RDV-TP complex. However, the ligand RMSD values were lower in SN00303170 and SN00359915 than the RDV-TP system and was higher in the SN00166900 system. Meanwhile, RMSF values of the protein were also monitored during 200 ns. High peaks were observed at V36 (V405), K61 (K430), T275 (T644), V368 (V737), T481 (T850), which corresponded to loop regions of the protein and V527 (V905) and E541 (E919), which was the carboxyl ends. While those residues interacting with ligands such as R555 (R186) and N691 (N322) were observed to be stable. Fig. 9D depicts the RMSF of protein residues for each complex.

Similar to the RMSF plot of RDV-TP, RMSF of RNA residues (Fig. 10 ) was high at A1 of the primer strand, which is the first sequence of the primer strand and was not involved in direct ligand interaction. The RMSF of U6 of primer strand and A3 and A4 of the template strand, which involved in the ligand binding, were observed to be stable.

Fig. 10

Plot of RNA residue fluctuation was recorded as RMSF values as RDV-TP (green) SN00166900 (cyan), SN00303170 (purple), and SN00359915 (pink).

Clustering analysis shows the most dominant structure for each complex during the 200 ns simulations. The ligand SN00166900 maintained H-bond interactions with U553 (U4) and U554 (U5) of the primer strand in addition to D83 (D452), R184 (R553), and R186 (R555). The manganese ion was stabilized through interaction with the ligand OH group. Binding of SN00303170 was preserved through interactions with D249 (D618), T187 (T556), A393 (A762), and H441 (H810). While, SN00359915 interacted with D83 (D452), R184 (R553), R186 (R555), T187 (T556), C253 (C622), S312 (S681), D315 (D684), and N322 (N691).

Binding free energy calculation

The total binding free energy was calculated from electrostatic (ΔE ELE), van der Waals (ΔE VDW), non-polar term of solvation energy (ΔE PB), and polar term of solvation energy (ΔE PBSUR). The contribution of each energy term was depicted in Table 1 . The total binding energy (ΔE PBTOT) was −18.90 ± 7.59 kcal mol−1 for RDV-TP, −16.62 ± 10.34 for ATP, −35.33 ± 7.56 for SN00166900, −45.23 ± 5.53 for SN00303170, and −64.86 ± 4.66 for SN00359915. It is clear that the SN00359915 had lowest total binding free energy out of the four ligands, and it is around three times as negative than that of RDV-TP. The large contribution of SN00359915 binding to RdRp was originated from van der Waals energy term.

Table 1

Binding energy predicted by MM-PBSA protocol calculated for 200 snapshots at last 100 ns trajectory.

Ligand	ΔEELE (kcal/mol)	ΔEVDW (kcal/mol	ΔEPB (kcal/mol)	ΔEPBSUR (kcal/mol)	ΔEPBTOT (kcal/mol)
SN00166900	−61.92 ± 6.39	−36.41 ± 8.57	69.72 ± 5.39	−6.72 ± 0.18	−35.33 ± 7.56
SN00303170	−26.81 ± 4.47	−52.53 ± 6.28	40.36 ± 4.09	−6.26 ± 0.27	−45.23 ± 5.53
SN00359915	−9.02 ± 4.88	−71.15 ± 5.54	22.35 ± 3.78	−7.04 ± 0.21	−64.86 ± 4.66
RDV-TP	−52.48 ± 30.14	−18.19 ± 6.16	56.19 ± 26.08	−4.42 ± 0.23	−18.90 ± 7.59
ATP	−40.99 ± 32.16	−27.11 ± 6.32	56.00 ± 25.73	−4.52 ± 0.13	−16.62 ± 10.34

In addition, the three hits identified in this study showed higher predicted binding energies than RDV-TP. The superiority of SN00359915 binding, for instance, which were supported by more diverse interactions with key residues such as those with R555 and N691, while was additionally corroborated by favorable interactions with D452, R553, T556, C622, S681, and D684.

Conclusion

In this study, RdRp in complex with RNA and manganese ions were built using homology modeling. The RdRp model was stable during 500 ns MD simulation. The binding of RDV-TP to RdRp showed that the position of the phosphate group was at the entry of NTP channel, which may reflect its role in inhibiting the entry of NTP to the active site. Additionally, the phosphate group established interactions with K551, R553, and K621, while the adenosine motif formed base pairing with U2 of template strand. The metal ions stabilized the phosphate group of RDV-TP and were located close to D618, D760, and D761. Using pharmacophore models, the three natural product hits SN00166900, SN00303170 and SN00359915 were observed to have were identified which had lower binding energies than that of RDV-TP.

Author contributions

MA performed simulation, analyzed data, and wrote the manuscript. STW, DJB and NP performed simulation. CW designed the work, analysis data, and supervised the work. All authors approved the manuscript.