Drug repurposing against the RNA-dependent RNA polymerase domain of dengue serotype 3 by virtual screening and molecular dynamics simulations

Abstract Dengue is an arboviral disease caused by the dengue flavivirus. The NS5 protein of flaviviruses is a potential therapeutic target, and comprises an RNA-dependent RNA polymerase (RDRP) domain that catalyses viral replication. The aim of this study was to repurpose FDA-approved drugs against the RDRP domain of dengue virus serotype 3 (DENV3) using structure-based virtual screening and molecular dynamics (MD) simulations. The FDA-approved drugs were screened against the RDRP domain of DENV3 using a two-step docking-based screening approach with Glide SP and Glide XP. For comparison, four reported DENV3 RDRP inhibitors were docked as standards. The hitlist was screened based on the docking score of the inhibitor with the lowest docking score (PubChem ID: 118797902; reported IC50 value: 0.34 µM). Five hits with docking scores and Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) energy lower than those of 118797902 were selected. The stability of the hit-receptor complexes was investigated using 100 ns MD simulations in an explicit solvent. The results of MD simulations demonstrated that polydatin and betiatide remained stably bound to the receptor, and formed stable interactions with the RDRP domain of DENV3. The hit-receptor interactions were comparable to those of 118797902. The average Prime MM-GBSA energy of polydatin and betiatide was lower than that of 118797902 during simulation, indicating that their binding affinity to DENV3 RDRP was higher than that of the standard. The results of this study may aid in the development of serotype-selective drugs against dengue in the future. Communicated by Ramaswamy H. Sarma


Introduction
Dengue virus is the causative agent of dengue haemorrhagic fever and severe dengue, which are associated with a higher risk of mortality. Studies have reported a marked global increase in the mortality rate due to dengue infections (Zeng et al., 2021). Using mathematical modelling, a study predicted that the number of individuals infected with dengue is expected to rise by 2.25 million by 2080 (Messina et al., 2019). The incidence of dengue during the COVID-19 pandemic in 2020 was reported to be higher than average in several countries, including Thailand, Singapore, Pakistan, Peru, and Ecuador (Brady & Wilder-Smith, 2021). Studies have also reported a marked rise in the number of dengue-SARS-CoV-2 coinfections, and the rising incidence of dengue during the pandemic crisis is likely to further cripple an already overburdened health system (Brady & Wilder-Smith, 2021). Therefore, there is a pressing need for developing drugs against dengue for combating outbreaks and epidemics, especially during the SARS-CoV-2 pandemic.
Previous studies have reported serotype-selective clinical differences in dengue viral infections, including haematological parameters, occurrence of internal haemorrhage, shock, disease severity, and the development of severe dengue (Perez et al., 2006;Yung et al., 2015). Additionally, numerous studies have demonstrated serotype-specific structural and functional differences in dengue viral proteins, including receptor usage for host entry, structural differences in the NS2B-NS3 protease complex, and intracellular localisation of the viral NS5 protein (Chandramouli et al., 2010;Hannemann et al., 2013;Thepparit & Smith, 2004). A previous study demonstrated that while the NS5 proteins of dengue serotypes 1 and 4 (DENV1 and DENV4) predominantly accumulate in the cytosol, those of dengue serotypes 2 and 3 (DENV2 and DENV3) largely localise in the nucleus (Hannemann et al., 2013). These serotype-specific structural and functional differences among dengue viral proteins led researchers to emphasise the importance of serotype-specificity during the identification of inhibitors. Metformin, an anti-hyperglycaemic drug, was repurposed as a pan inhibitor of dengue, and was prescribed as an adjunctive treatment strategy for diabetic individuals (Htun et al., 2018). However, a recent study demonstrated that metformin is ineffective as a pan inhibitor of dengue viral infections (Cheang et al., 2021). These findings suggest that serotype-selective inhibitors of dengue would have higher specificity and efficacy against dengue than pan inhibitors.
A recent study demonstrated that the mRNAs of DENV2 and DENV3 persist for longer durations in wastewater, over a wide temperature range, and are therefore more likely to cause wastewater-mediated outbreaks of dengue (Chandra et al., 2021). Stagnant wastewater is one of the major modes of dengue viral spread, and as DENV2 and DENV3 persist for longer durations in wastewater, these serotypes have a greater potential of causing outbreaks of dengue. Previous studies have reported an increasing spread of DENV3 worldwide, in addition to the emergence of novel subtypes of DENV3, in the past few decades (Franco et al., 2010;Messer et al., 2003). DENV3 has also been responsible for the occurrence of sudden outbreaks of dengue, severe dengue, and epidemics in several regions, including the Indian subcontinent, Africa, and Latin America, and is the predominant circulating serotype in certain regions (Faye et al., 2014;Messer et al., 2003;Tatura et al., 2021). Despite the rise in the number of sudden outbreaks of DENV3 infections, there are few studies on the development of inhibitors against DENV3. Although DENV2 is potentially more harmful than DENV3, the increasing global spread of DENV3, combined with the occurrence of sudden epidemics and the scarcity of studies, warrants the identification of DENV3 inhibitors.
Drug repurposing involves the use of existing approved drugs for treating a disease different from that for which it is originally intended. As the absorption, distribution, metabolism, excretion, and toxicity (ADMET) profiles of existing drugs are already known, this reduces the cost and time required for drug discovery, by bypassing phase I clinical trials (Nosengo, 2016). Drug repurposing can reduce the time and cost of drug discovery by more than half, and this strategy is currently being used or studied for treating SARS-CoV-2 infection (Pashaei, 2021;Singh et al., 2020). Virtual screening strategies are frequently employed for the discovery of novel inhibitors of the dengue virus. A recent study identified novel inhibitors of dengue viral envelope protein by combining structure-based pharmacophore modelling with MD simulation trajectory analyses (Hengphasatporn et al., 2019).
The NS5 protein of flaviviruses contains an RNA-dependent RNA polymerase (RDRP) domain and a methyltransferase domain. The RDRP activity is crucial for the replication of flaviviruses, including the dengue virus, and the NS5 protein of dengue is a potential therapeutic target (Lim et al., 2015;Malet et al., 2008). Inhibition of the RDRP activity of dengue inhibits viral replication, and this strategy has been employed for the identification of inhibitors of the dengue virus (Picarazzi et al., 2020). Using virtual screening combined with pharmacophore filtering and MD simulations, a recent study identified one lead against the RDRP domain of DENV3 from a library of 428 usnic acid derivatives (Roney et al., 2021). The aim of this study was to repurpose FDA-approved drugs against the RDRP domain of DENV3 by structure-based virtual screening, and adjudged with 100 ns MD simulations. All the FDA-approved drugs in the DrugBank database were screened against the RDRP domain of DENV3 using a twostep docking and scoring approach. The hitlist was screened based on the docking score of a reported inhibitor of DENV3 RDRP (PubChem ID: 118797902; reported IC 50 value: 0.34 mM) (Obi et al., 2021;Yokokawa et al., 2016). Top five scoring hits with docking scores and Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) energy lower than that of 118797902 were selected. The stability of the hit-protein complexes was investigated with 100 ns MD simulations, and the values of root mean square fluctuation (RMSF) and root mean square deviation (RMSD) of the hits were compared to those of the reported inhibitor, 118797902. The stability of the protein backbone following hit binding was compared to that of the protein bound to the standard inhibitor. The ligand-receptor interactions that were important for complex stability were identified from the occupancy of the interactions throughout the trajectory. The hits were finally screened on the basis of complex stability and the MM-GBSA energy averaged over the trajectory. The MM-GBSA energy of the hits was compared to that of 118797902, and hits with MM-GBSA energy lower than that of the standard inhibitor were selected. As DENV3 is a potential agent for sudden outbreaks of dengue, the results of this study may aid future drug development studies against DENV3.

Drug library and receptor structure preparation
All the FDA-approved drugs in the DrugBank database, version 5.1.8, were retrieved in SDF format (Wishart et al., 2006). The small fragments and drugs that overlapped with the illicit and withdrawn drug categories were removed, and the remaining compounds were prepared with LigPrep, Schr€ odinger, using the OPLS4 force field Schr€ odinger Release 2021-2: LigPrep, Schr€ odinger, LLC, New York, 2021. The ionisation state and tautomeric state of the compounds were kept unaltered during ligand preparation, so as to ensure that the active structure of the drugs was not altered following ligand preparation. The specific chiralities were retained, and a maximum of 32 conformations were generated for each ligand. The lowest energy conformation of each drug was retained, and a library of 1623 drugs were finally obtained after ligand preparation.
The PDB structure 6IZZ, comprising the structure of the RDRP domain of DENV3 NS5 bound to the small molecule inhibitor, RK-0404678, was retrieved from the PDB (Berman et al., 2000;Shimizu et al., 2019). The resolution and R-value of the structure was 1.97 Å and 0.255, respectively. The PDB structure 6IZZ was selected in this study as the resolution of this structure was the highest among the other inhibitorbound wild-type structures of DENV3 RDRP available in the PDB, including 6J00 and 5IQ6. Although there are several structures of DENV3 RDRP in the PDB with a resolution higher than that of 6IZZ, including 5JJR, 5JJS, and 5I3Q, these structures contain one or more mutations, and were therefore not considered in this study. The PDB structures 5F3Z, 5F3T, 5F41, 5HMW, 5HMX, 5HMY, 5HMZ, 5HN0, 5I3P, 5JJR, 6H80, and 6H9R contain one or more mutations and have a resolution lower than that of 6IZZ, and were thus not considered in this study. Prior to docking, the structure was prepared using the Protein Preparation Wizard, Schr€ odinger (Madhavi Sastry et al., 2013). The missing side chains and loops were filled in using Prime, Schr€ odinger (Jacobson et al., 2002(Jacobson et al., , 2004. The hydroxyl groups and the protonation states of Asn, Gln, and His were optimised using ProtAssign. The water molecules that were 3.0 Å beyond the binding site of the standard inhibitor, RK-0404678, were removed. The hydrogen atoms were finally minimised using the OPLS4 force field, while keeping the heavy atoms restrained.

Grid preparation, validation of docking protocol
The grid was built around the binding site of the small molecule inhibitor, RK-0404678, in the crystal structure 6IZZ. The x, y, and z coordinates of the centre of the grid were À15.11, À20.02, and 39.17, respectively, and the dimensions of the grid were 23 Å Â 23 Å Â 23 Å along the x-, y-, and zaxes, respectively. The docking protocol was initially validated by redocking the small molecule inhibitor, RK-0404678, in the PDB structure, to its binding site using the same grid and docking parameters in Glide SP and Glide XP (Friesner et al., 2006;Halgren et al., 2004). The RMSD between the experimentally-determined and docked poses was calculated for determining the ability of the docking protocol in predicting the native binding mode of RK-0404678. The docking protocol was initially validated by analysing the interactions of RK-0404678 with the binding site residues and comparing with those reported in literature (Shimizu et al., 2019). The docking protocol was further validated by docking decoys of RK-0404678. A total of 10 decoys of RK-0404678 were selected from the Asinex and ZINC natural product databases based on molecular descriptors, using DecoyFinder 2.0. The minimum Tanimoto coefficient between the decoys and RK-0404678 was set to 0.90. The parameters used for decoy selection was provided in Supplementary Table S1. The decoys were docked to the binding site of DENV3 RDRP using the same docking and scoring parameters used for docking RK-0404678. The docking scores and ligand efficiency values of the decoys were compared to those of RK-0404678 for evaluating the discriminatory power of the scoring function.

Virtual screening
The prepared drug library was next screened against the RDRP domain of DENV3 using a two-step docking protocol. In the first step, the library was screened using Glide SP. In the second step, the top scoring 10% compounds from Glide SP were re-docked using Glide XP, and their binding energy was determined using Prime MM-GBSA version 3.0 (Prime Release 2021-2; Schr€ odinger, 2021). The VSGB continuum solvation model, which includes residue-dependent effects, was used for calculating the binging energy, using water as the solvent (Li et al., 2011). The OPLS4 force field was used for the MM-GBSA energy calculations. In addition, four experimentally reported standard small molecule inhibitors of DENV3 RDRP (PubChem IDs: 118797902, 118797901, 122172828, and 121232415 with IC 50 values of 0.34, 1.7, 1.5, and 0.048 mM, respectively), were docked to the binding site of DENV3 RDRP (Lim et al., 2015;Obi et al., 2021;Tarantino et al., 2016;Yokokawa et al., 2016). For hit selection, the hitlist was screened based on the Glide XP score of 118797902, which had the lowest docking score among the reported inhibitors. Top five scoring hits with Glide XP scores and MM-GBSA energy lower than those of 118797902 were selected and subjected to subsequent MD simulations for studying complex stability. The ADMET properties of the selected hits, including CNS activity, HERG channel blockage potential, oral absorption, intestinal permeability, skin permeability, binding potential to human serum albumin, were determined using QikProp, Schrodinger (Schr€ odinger Release 2021-2: QikProp, Schr€ odinger, LLC, New York, 2021). The number of violations of Lipinski's rule of five were also determined with QikProp.

MD simulations
The stability of the ligand-receptor complexes of the selected hits was determined using MD simulations in Desmond, Schr€ odinger, for 100 ns in explicit solvent using the OPLS 2005 force field (Bowers et al., 2006). For comparison, the ligand-receptor complex of the standard inhibitor, 118797902, was also subjected to 100 ns MD simulations. The complexes were first prepared in Desmond by solvating in an orthorhombic water box of TIP4P solvent with a buffer of 10 Å thickness on all sides. The solvated systems were then neutralised by the addition of the requisite number Clions. Finally, Na þ and Clions were added to the systems until the concentration of the solution was 0.15 M. The prepared systems were initially minimised using a multi-step minimisation protocol. Briefly, the systems were first minimised in the NVT ensemble for 100 ps, while keeping the heavy atoms of the solute restrained. The systems were then minimised in the NPT ensemble for 12 ps, while keeping the heavy atoms of the solute restrained. The temperature of the systems was slowly raised to 310 K, following which the systems were equilibrated for 24 ps in the NPT ensemble at a temperature of 310 K and pressure of 1 bar. Following equilibration, the systems were simulated for 100 ns in the NPT ensemble, under periodic boundary conditions. The temperature of the systems was maintained at 310 K using a Nos e-Hoover chain thermostat and a constant pressure of 1.01325 bar using a Martyna-Tobias-Klein barostat. The cutoff radius for short-range interactions was set to 9.0 Å, while the long-range interactions were calculated using the Particle Mesh Ewald method (Darden et al., 1993). The timestep was set to 2.0 fs.

Trajectory analyses and calculation of MM-GBSA binding energy from trajectory
The stability of the complexes was assessed by trajectory visualisation, the values of RMSD and RMSF. The values of RMSF and RMSD were compared with those of the standard inhibitor, 118797902 (Obi et al., 2021). The occupancy of the ligand-receptor interactions throughout the entire trajectory and in the last 20 ns were determined and compared to that of 118797902. The absolute binding energy of the hits and the standard inhibitor during simulation was determined using the Prime MM-GBSA approach. The MM-GBSA energy was calculated from the trajectories with the thermal_ mmgbsa.py script, using Prime MM-GBSA version 3.0, Schr€ odinger (Prime Release 2021-2; Schr€ odinger, 2021). The MM-GBSA binding energy of the hits was averaged over the entire trajectory and compared to that of the standard inhibitor, 118797902.

Validation of docking protocol
The docking protocol was initially validated by redocking the bound inhibitor, RK-0404678, in the crystal structure 6IZZ using the same grid and docking parameters used for screening the drug library. The RMSD between the pose of RK-0404678 in the crystal structure and the docked poses generated by Glide SP and Glide XP was 1.0951 Å and 0.2141 Å, respectively (Supplementary Figure S1b). This indicated that there was a good correlation between the binding conformation of RK-0404678 in the experimentally-determined crystal structure and the binding pose predicted by docking with Glide. This indicated that the docking algorithm was capable of accurately predicting the native binding pose of ligands in the receptor binding site.  Figure S1c), which corroborated with those reported in literature, thereby further emphasising the accuracy of the docking protocol (Shimizu et al., 2019). The Glide XP score and MM-GBSA energy of RK-0404678 was À2.687 and À36.757, respectively. The docking protocol was further validated by docking decoys of RK-0404678 with the same grid and docking parameters. The results demonstrated that the docking scores of the decoys were higher than that of RK-0404678 (Supplementary Table S2), indicating that the binding affinity of the decoys was lower than that of RK-0404678. Therefore, the results indicated that the docking protocol and scoring function was capable of differentiating between decoys and the inhibitor, RK-0404678.

Virtual screening and hit selection
The prepared library of FDA-approved drugs was next screened using Glide SP, and the 10% top scoring compounds from Glide SP were re-docked using Glide XP, and their MM-GBSA energy was determined with Prime. The Glide XP hitlist was screened on the basis of the Glide XP docking score (-6.698) and MM-GBSA energy (-31.622) of the standard inhibitor, 118797902. Finally, top 5 scoring hits, including tiaprofenic acid, ibandronate, polydatin, betiatide, and mupirocin, with Glide XP scores À6.698 and MM-GBSA binding energies À31.622 were selected. The docking scores, ligand efficiency values, and MM-GBSA binding energies of the hits are enlisted in Table 1. The Glide XP docking scores and MM-GBSA energy of the hits were markedly higher than those of the standard inhibitor. As the molecular weights of the hits and standard inhibitors were comparable, the lower ligand efficiency values of the hits indicated that the binding affinity of the hits normalised to molecular weight was higher than that of the standard inhibitors (Table 1).
The ADMET properties of the five selected hits, as determined using QikProp, are provided in Table 2, and indicates that tiaprofenic acid, ibandronate, and betiatide had no violations of Lipinski's rule of five, while polydatin and mupirocin had only one violation of Lipinski's rule of five. The ADMET properties were compared to the range of 95% drugs in QikProp. None of the hits were predicted to have CNS activity, with the exception of tiaprofenic acid, which was predicted to have very slight CNS activity ( Table 2). As the hits are not required to act in the CNS compartment, the fact that they did not possess CNS activity was not a great  concern. In QikProp, the HERG channel inhibition potential of compounds is calibrated in such a way that it becomes a concern below a value of À5. With the exception of polydatin, none of the hits had possible HERG blockade potential, and the value for polydatin was predicted to be À5.606, which was not very low compared to the threshold value of À5, indicating that the concern of HERG channel blockade by polydatin was not very high. Tiaprofenic acid had high skin permeability, while ibandronate had low skin permeability, as reflected in the values of log K p . On the other hand, polydatin, betiatide, and mupirocin had moderate skin permeability. All the hits were predicted to have low human serum albumin binding potential, as revealed by the values of QPlogKhsa. Tiaprofenic acid was predicted to have high oral absorption, while polydatin, betiatide, and mupirocin was predicted to have medium human oral absorption. On the other hand, the human oral absorption potential of ibandronate was predicted to be low (Table 2).

Hit-RDRP interactions
The ligand-receptor interactions of the five selected hits were determined using Maestro, Schr€ odinger, for analysing whether they interacted with the same residues in the RDRP domain of DENV3 NS5 as the standard inhibitor, 118797902 ( The results of ligand-receptor interaction analysis revealed that the hits interacted with the catalytic residues of RDRP, and some of the hit-receptor interactions were comparable to those of the standard inhibitor, 118797902 (Figure 1 and Supplementary Figure S1d). The stability of the interactions of the hits and 118797902 with the binding site of DENV3 RDRP was further assessed with 100 ns MD simulations, and the occupancy of the interactions was determined by trajectory analysis for identifying the interactions that were crucial for complex stability. Ligand-protein interactions of higher occupancy are considered to be more stable, and therefore more important for complex stability.

Trajectory analyses
The 5 hit-receptor complexes and the ligand-receptor complex of the standard inhibitor, 118797902, were simulated for 100 ns for assessing complex stability, determining the occupancy of the ligand-receptor interactions, and for comparing the values of RMSD and RMSF and ligand-receptor interactions of the hits and standards. The trajectories were visualised in VMD 1.9.3 and the values of RMSF and RMSD of the protein and ligands were separately determined with Desmond and VMD (Bowers et al., 2006;Humphrey et al., 1996). The stability of the ligand-receptor interactions was analysed by determining the occupancy of the interactions using the simulation interactions diagrams tool in Desmond.

RMSF and RMSD analyses, complex stability, ligand torsions
The local fluctuations in the backbone atoms of DENV3 were determined from the values of RMSF. The RMSF values of the backbone atoms of the unliganded DENV3 protein were compared with those of the protein bound to the hits and standard inhibitor, 118797902. The results demonstrated that the RMSF values of the DENV3 protein complexed with the hits and 118797902 were very similar to those of the unliganded protein (Figure 2a). This indicated that ligand binding did not alter the local residue fluctuations in the protein.
As depicted in Figure 2a, the RMSF values of the terminal residues were higher than those of the other residues of the protein. This could be attributed to the fact that these residues were terminally located, and therefore less restricted, and could undergo higher fluctuations than the other residues. The RMSF values of the regions with secondary structure elements, including alpha helices and beta sheets, were lower than those of the uncoiled regions. As beta sheets and alpha helices are more rigid than unstructured regions, they undergo less fluctuations than regions without regular secondary structure elements, which accounted for the lower values of RMSF in these regions.
The RMSD values of the protein backbone and ligands was determined throughout the trajectory, in Desmond, and are graphically represented in Figure 3. As depicted in Figure  3a, the values of RMSD of the backbone atoms of the unliganded protein became steady at $38 ns, and remained steady (RMSD range $ 0.4 Å) until the end of the production run. This implied that the unliganded protein stabilised after 38 ns of the production run and remained stable thereafter, indicating that the system had reached equilibrium after 38 ns. The RMSD values of the backbone atoms of DENV3 RDRP bound to the standard (118797902), tiaprofenic acid, ibandronate, polydatin, betiatide, and mupirocin, became steady after 12 ns, 13 ns, 60 ns, 22 ns, 25 ns, and 16 ns, respectively, indicating that the systems had reached equilibrium and remained stable after the indicated time periods (Figure 3a). The fluctuations in the RMSD values of the backbone atoms of DENV3 bound to the standard (118797902), tiaprofenic acid, ibandronate, polydatin, betiatide, and mupirocin following equilibration was in the range of 0.4-0.8 Å (Figure 3a).
The RMSD values of the standard (118797902) underwent a major fluctuation around 30 ns, as depicted in Figure 4b. Trajectory visualisation revealed that this fluctuation was attributed to a conformational change at around 30 ns, following which the conformation of 118797902 was preserved. However, the RMSD values of the standard underwent high fluctuations after 30 ns (Figure 3b). Trajectory visualisation and analysis of the torsion plot of 118797902 revealed that this was attributed to the marked rotational motion around the C10-C17 bond (Figure 4a), which resulted in marked orientational fluctuations in the 3-(thiophen-2-yl)prop-2-yn-1ol group during the production run. However, trajectory visualisation revealed that although the 3-(thiophen-2-yl)prop-2yn-1-ol moiety underwent orientational fluctuations, 118797902 remained stably bound to the binding site of DENV3 RDRP and did not dissociate from the complex. The  high values of per-atom RMSF of the 3-(thiophen-2-yl)prop-2yn-1-ol moiety was attributed to the conformational fluctuations of this moiety around the C10-C17 bond ( Figure 5).
The RMSD values of tiaprofenic acid became steady after 52 ns, and the changes in the RMSD values were in the range of 0.4-0.6 Å. Trajectory visualisation revealed that tiaprofenic acid remained stably bound to the catalytic site of DENV3 RDRP after 10 ns, until the end of the production run. However, there was a sharp fluctuation in the RMSD values between 64 and 72 ns (Figure 3b), which was attributed to conformational changes in the molecule, arising from rotations around the C5-C6, C5-C12, C7-C11, and C11-C13 bonds (Figure 4a). The rotational motion of these bonds contributed to the high values of RMSF of the terminal atoms of tiaprofenic acid ( Figure 5). As depicted in Figure 4a, the rotational motion of the C5-C6 and C5-C12 bonds were more pronounced than that of the C7-C11 and C11-C13 bonds. This could be attributed to the fact that the C7-C11 and C11-C13 bonds were bound by bulky aromatic groups and their rotational motion was consequently restricted by the presence of these bulky groups (Figure 4a). However, as the C5-C6 and C5-C12 bonds were terminally located, their motion was unhindered, which accounted for their higher torsional motion and high values of RMSF of the terminal O2 and O4 atoms (Figures 4a and 5b).
Trajectory visualisation revealed that ibandronate did not dissociate from the receptor throughout the production run, which was consistent with the values of RMSD ( Figure 3b). As depicted in Figure, the RMSD values of ibandronate became steady after 45 ns. However, there was a sharp fluctuation in the RMSD values of the ligand between 80 ns and 90 ns. Analysis of the torsion plot revealed that this was attributed to the rotational motion of the terminal hydroxyl groups (Figure 4b), which resulted in the high values of RMSF of these atoms (Figure 5b). Trajectory visualisation revealed that polydatin remained stably bound to the receptor throughout the simulation. However, polydatin underwent a translational shift around 43 ns, and was buried more deeply into the binding cleft of DENV3 RDRP, which corroborated with the sharp fluctuation in the RMSD values of polydatin around 43 ns, following which the RMSD values decreased and remained stable (Figure 3b). Visualisation of the trajectory revealed no marked changes in the binding site or conformation of polydatin beyond 40 ns, and that the compound remained deeply buried in the binding site thereafter, indicating that polydatin formed a stable complex with DENV3 RDRP.  Trajectory analyses revealed that betiatide remained stably bound to the binding pocket of DENV3 RDRP, and became more deeply buried towards the end of the trajectory, indicating the formation of a stable complex with DENV3 RDRP. This was consistent with the RMSD values of betiatide, which became steady after 20 ns (Figure 3b). The fluctuations in the RMSD values of betiatide were in the range of 0.4-0.8 Å following equilibration, which gradually lessened over time, indicating that the betiatide-DENV3 complex became steadier towards the end of the production run (Figure 3b). On the other hand, although mupirocin did not dissociate from the binding site, the compound did not form a very stable complex with the receptor, as observed from the values of RMSD. Trajectory visualisation revealed that the fragment containing the cyclic oxane and oxirane moieties (C1-C16) were buried in the binding site; however, the 9-hydroxynonanoic acid moiety of mupirocin (C16-O26) remained outside the binding site and underwent marked conformational fluctuations during simulation. This was attributed to the increasing torsional conformation of the rotatable bonds in the 9-hydroxynonanoic acid moiety, especially after 78 ns, and also to the rotational motion of the 9-hydroxynonanoic acid moiety around the C14-O15 bond ( Figure 4c). As depicted in Figure 3b, the RMSD values of mupirocin underwent sharp fluctuations, in the range of 8 Å, resulting in the sharp fluctuations in the values of RMSF (Figures 4c and 5b).

Analysis of ligand-protein interactions over the trajectory
The occupancy of the ligand-protein interactions was determined throughout the 100 ns trajectory and in the last 20 ns of the production run. For the calculation of hydrogen bonds, the cut-off distance between the donor and acceptor atoms was set to 2.5 Å. Analysis of the occupancy of the ligand-receptor interactions revealed that tiaprofenic acid formed 4 hydrogen bonds with Lys 756, Cys 780, Val 783, and Asp 808 in the last 20 ns, when the system had reached equilibrium (Table 3, Figure 6). Of these, the hydrogen bond with Lys 756 was stable and crucial for complex stability, as  (Obi et al., 2021;Yokokawa et al., 2016). indicated by the high occupancy (78%) of this interaction ( Figure 6, Table 3). The number of frames in which the interactions were observed are provided in Supplementary Table  S4. The aromatic phenyl ring of tiaprofenic formed a favourable pi-pi stacking interaction with the aromatic phenyl moiety of Tyr 882 in the binding pocket. However, this interaction was only observed for 30% of the trajectory, and the low occupancy of this interaction could be attributed to the torsional motion of the C11-C13 bond (Figure 4a). Tiaprofenic acid also formed hydrophobic interactions with Leu 880 and Tyr 882, of occupancies 30.8% and 88.6%, respectively ( Figure 6). The high occupancy of the latter hydrophobic interaction with Tyr 882 revealed that this interaction was highly stable and crucial for stabilising the complex. Although tiaprofenic acid remained bound to the receptor binding site throughout the production run, the occupancies of some of the ligand-receptor interactions were not very high owing to the rotational motion of the phenyl, thiophene, and carboxylic acid moieties around the C11-C13, C7-C11, C5-C6, C5-C12, and O2-C12 bonds (Figure 4a).
Although ibandronate formed eight hydrogen bonds in the last 20 ns following equilibration, the occupancy of these interactions was not high (Table 3). Of these eight hydrogen bonds, three of the hydrogen bonds with Lys756, Thr 805, and Met 809 were water-mediated ( Figure 6). Occupancy analysis revealed that the hydrogen bonds with Lys 756 and Glu 807 were the highest, being 69% and 56%, respectively, implying that none of the interactions between ibandronate and DENV3 were highly stable ( Figure  6, Table 3). Trajectory analysis revealed that the low occupancies of the hydrogen bonds were attributed to the rotational motion of the terminal hydroxyl groups, which prevented the formation of stable interactions with high occupancy (Figure 4b).
Polydatin formed two water-mediated hydrogen bonds with the side chain of Ser 776, of occupancy > 30% ( Figure  6, Table 3). In the last 20 ns, when the system had reached equilibrium, the occupancy of these hydrogen bonds increased to 92% and 64%, indicating that these interactions were highly stable and crucial for complex stability. The complex was also stabilised by hydrophobic interactions with Met 809, Thr 832, and Tyr 882, of which the hydrophobic interactions with Tyr 882 were critical for stabilising the complex, as indicated by the high frequency (111.9%) of these interactions (Table 3).
Betiatide formed 8 hydrogen bonds with the catalytic site of DENV3, of which 4 hydrogen bonds with Lys 756, Asn 777, Asp 808, Tyr 882, and Trp 833 were water-mediated (Table 3, Figure 6). The occupancy of the hydrogen bonds with Asn 777, Asp 808, Trp 833, Tyr 882 was 77%, 73%, 60%, and 82%, respectively, indicating that these interactions were highly stable and crucial for complex stability. The occupancy of the hydrogen bonds with Asn 777, Asp 808, Trp 833, Tyr 882 increased to 98%, 92%, 78%, and 98%, respectively, in the last 20 ns, which implied that these interactions became more stable after the systems had reached equilibrium ( Figure 6). Analysis of the hydrogen bond occupancy also revealed the formation of an intramolecular hydrogen bond between O2 and N9, of occupancy 69% ( Figure 6). The occupancy of this hydrogen bond was found to be 92% in the last 20 ns, indicating that this intramolecular interaction was stable. Analysis of the ligand-protein interactions of mupirocin revealed that the non-polar alkyl moiety of mupirocin did not fit well in the polar binding site. This was reflected in the occupancy of the hydrogen bonds between mupirocin and DENV3, became lower towards the end of the run, indicating that these interactions were not very stable towards the end of the simulation ( Figure 6, Table 3).

MM-GBSA analyses
The free energy of binding of the hits were determined throughout the trajectory, using the Prime MM-GBSA approach. The results demonstrated that the MM-GBSA binding energy of the hits determined from the trajectory were comparable to that determined from the docked hit-receptor complexes, with the exception of ibandronate (Tables 1 and  4). The highest contributor to the MM-GBSA energy of the standard inhibitor, 118797902, was the van der Waals energy term, followed by the lipophilic energy term and coulomb energy term (Table 4). The results of MM-GBSA analysis of tiaprofenic acid revealed that the van der Waals energy term was the most negative, implying that this energy term had the highest contribution to the binding energy of tiaprofenic acid with DENV3 RDRP (Table 4). The Coulomb energy and lipophilic energy terms also contributed to the free energy of binding of tiaprofenic acid, but their contribution was lower than the van der Waals energy term (Table 4). On the other hand, the coulomb energy and van der Waals terms had the highest contributions to the absolute binding energy of ibandronate (Table 4), while the lipophilic and hydrogen bond terms had very little contribution. Analysis of the MM-GBSA energy of polydatin revealed that the van der Waals term had the highest contribution to the binding free energy, followed by the lipophilic energy term. The major contributor to the MMGBSA energy of betiatide and mupirocin was the van der Waals energy term, followed by the coulomb and lipophilic energy terms (Table 4). The results of MM-GBSA analysis revealed that the free energy of binding of the hits were lower than that of the standard inhibitor, RK-0404678 (-36.757), indicating that the binding affinity of the hits to the receptor was higher than that of RK-0404678.

Discussion
Dengue continues to be classified as a neglected tropical disease, and the number of cases and deaths due to dengue have recently increased worldwide, especially during the SARS-CoV-2 pandemic (Brady & Wilder-Smith, 2021;Horstick et al., 2015;Zeng et al., 2021). DENV3 has a higher potential of causing wastewater-mediated dengue outbreaks, and there is an increasing global spread of novel DENV3 subtypes capable of causing sudden dengue outbreaks (Chandra et al., 2021;Franco et al., 2010;Messer et al., 2003). Therefore, DENV3 poses as a potential agent for unexpected future dengue epidemics. Despite this concern, there are no definitive treatments against DENV3 to date. Previous studies have indicated that serotype-selective inhibitors would have higher specificity than pan inhibitors of dengue. This study aimed to perform drug repurposing against the RDRP domain of DENV3 using computational approaches. The strategy involves the identification of novel targets of existing approved drugs, and can reduce the time and cost of drug development by more than half (Nosengo, 2016).
This study aimed to screen FDA-approved drugs against the RDRP activity of DENV3 by employing structure-based virtual screening and MD simulations. To this end, the FDAapproved drugs were retrieved from DrugBank and screened against the RDRP domain of DENV3 using Glide SP and Glide XP. The five top scoring hits, including tiaprofenic acid, ibandronate, polydatin, betiatide, and mupirocin, were selected by docking-based screening, by comparing their docking scores and MM-GBSA energy with those of a standard inhibitor, 118797902. The docking scores and MM-GBSA energy of these hits were lower than those of 118797902, implying that the binding affinity of these hits to the receptor was higher than that of 118797902. The stability of the hit-receptor complexes was determined by 100 ns MD simulations in explicit solvent and compared to that of 118797902. The hits were further screened on the basis of their MM-GBSA energy, and hits with MM-GBSA binding energy lower than that of 118797902 were selected, as a lower value of MM-GBSA energy indicates higher binding affinity.
Trajectory analyses revealed that all the five hits selected from structure-based virtual screening remained bound to the catalytic site of DENV3 RDRP throughout the simulation and did not dissociate from the complex. Although tiaprofenic acid formed stable interactions with Lys 756 and Tyr 882, reflected in the high occupancy of these interactions (Table 3), the MM-GBSA energy of tiaprofenic acid was higher than that of 118797902. This indicated that the binding affinity of tiaprofenic acid with DENV3 RDRP was lower than that of 118797902. Ibandronate, however, did not form stable interactions with DENV3 RDRP, as revealed by the low occupancy of the ligand-receptor interactions (Table 3). The low occupancy of the ligand-receptor interactions for ibandronate were attributed to the rotational motion of the terminal groups, as revealed by trajectory visualisation. In addition, the MM-GBSA energy of ibandronate was lower than that of 118797902 (Table 4), implying that the binding affinity of ibandronate to DENV3 RDRP was lower than that of 118797902.
Trajectory analysis revealed that both polydatin and betiatide remained stably bound to the receptor and formed highly stable interactions with DENV3 RDRP, as revealed by the high occupancy of the ligand-receptor interactions. Polydatin formed stable interactions with Ser 776 and Tyr 882, while betiatide formed stable interactions with Asn 777, Asp 808, Trp 833, and Tyr 882. The interacting residues of polydatin and betiatide, including Asn 777, Cys 780, and Trp 833, overlapped with those of the standard inhibitor, 118797902, implying that polydatin and betiatide interacted with the same residues of the receptor as 118797902 did. Additionally, the values of MM-GBSA energy of betiatide was much lower than that of 118797902, indicating that its binding affinity was much higher than that of the standard inhibitor, 118797902. However, the binding affinity of polydatin was slightly higher than that of 118797902, as reflected in the values of MM-GBSA energy. Although the MM-GBSA energy of mupirocin was much lower than that of 118797902, trajectory visualisation revealed that its 9-hydroxynonaneic acid moiety underwent abrupt fluctuations in the receptor binding site. Although mupirocin did not dissociate from the complex, the molecule underwent conformational fluctuations and did not form highly stable interactions with the binding site residues. Analysis of the ligand-receptor interactions revealed that mupirocin formed only one stable interaction with Glu 807. Binding free energy analysis revealed that the MM-GBSA energy of mupirocin was lower than that of the standard inhibitor, 118797902, implying higher binding affinity.
Altogether, the results demonstrated that polydatin and betiatide formed highly stable interactions with the catalytic site of DENV3, and the interacting residues were comparable with those of the standard inhibitor, 118797902. However, the occupancy of the interactions of polydatin and betiatide was higher than those of 118797902, implying that the stability of these hit-receptor complexes was higher than that of 118797902. This was further supported by the results of trajectory visualisation, which revealed that betiatide and polydatin remained stably bound, while 118797902 underwent conformational fluctuations during simulation. Polydatin is a naturally occurring stilbene glycoside and precursor of resveratrol that is found in Polyganum cuspidatum (Wang et al., 2013). Polydatin inhibits melanogenesis, and a topical formulation of polydatin has been approved for use as a whitening agent. In addition, polydatin is being presently evaluated in phase II clinical trials for treating septic or haemorrhagic shock (Jeong et al., 2010;Wishart et al., 2006). Recent studies have reported that polydatin possesses antiviral activity against SARS-CoV-2, in addition to its antiinflammatory, anti-oxidant, and anti-cancer properties (Bai et al., 2021;Wiwanitkit, 2021;Wu et al., 2020). Betiatide has been approved for intravenous use for preparing a radiopharmaceutical for renal imaging and assessing renal function (Vera et al., 2001;Wishart et al., 2006). Analysis of the binding free energy of the hits from the trajectory using the Prime MM-GBSA approach revealed that the MM-GBSA energy of polydatin and betiatide were lower than that of the standard inhibitor, 118797902 (Table 4). This indicated that the binding affinity of polydatin and betiatide to DENV3 RDRP was higher than that of the experimental standard, 118797902.
Trajectory analyses revealed that the hits selected herein might serve as potential inhibitors of DENV3 RDRP by occluding the catalytic site of the enzyme by interacting with the active site residues. However, further experimental studies are necessary for validating the inhibitory potential of polydatin and betiatide against the RDRP activity of DENV3.

Conclusion
The aim of this study was to repurpose FDA-approved drugs against the RDRP domain of DENV3 using docking-based virtual screening and adjudged through MD simulations. In this study, all FDA-approved drugs were screened against the RDRP domain of DENV3 using a two-step docking and screening protocol. Following hit selection, the hit-receptor complexes of the top five scoring hits were subjected to 100 ns MD simulations for determining complex stability. The ligand-receptor interactions that were crucial for complex stability were identified by trajectory analyses and compared to those of a reported inhibitor of DENV3 RDRP, 118797902, selected as standard. The results of this study demonstrated that polydatin and betiatide formed highly stable interactions with DENV3 RDRP. The MM-GBSA energy of the selected hits was lower than that of 118797902, implying that the predicted binding affinity of the hits was higher than that of the reported inhibitor. The results of this study may prove to be useful for the development of serotypeselective drugs against dengue in the future.

Disclosure statement
No potential conflict of interest was reported by the authors.