Detailed in silico Evaluation of WNV Proteins: Dynamic and Thermodynamic Insights into Doravirine as a Potential Multitarget Agent

Introduction

West Nile Virus (WNV) is one of the most widespread and emerging mosquito-borne viral (MBV) zoonoses in Europe where climatic factors, such as temperature and precipitation, play a crucial role in the distribution and transmission of WNV. As climate projections indicate more severe climate extremes, such as drier and warmer springs and summers, the risk of WNV outbreaks is expected to rise.¹ An estimated 70–80% of human WNV infections are subclinical or asymptomatic and symptomatic persons experience an acute systemic febrile illness that often includes headache, weakness, myalgia, or arthralgia. Less than 1% of infected persons develop neuroinvasive or severe disease (WNND), which typically manifests as meningitis, encephalitis, or acute flaccid myelitis. Rarely, cardiac dysrhythmias, myocarditis, rhabdomyolysis, optic neuritis, uveitis, chorioretinitis, orchitis, pancreatitis, and hepatitis have been described in patients with WNV disease.^2,3

WNV disease should be considered in any person with an acute febrile or neurological illness who has had recent exposure to mosquitoes, blood transfusion, or organ transplantation, especially during the summer months in areas where virus activity has been reported. In this context, no WNV vaccines are licensed for use in humans and there is no specific treatment for WNV disease, and the clinical management of severe cases is supportive.^4,5

WNV belongs to the Flaviviridae family and the Flavivirus genus. The spherical WNV particle, approximately 50 nm in diameter, comprises a host-derived lipid bilayer envelope surrounding a nucleocapsid core, which contains a single-stranded, positive-sense RNA genome of about 11,000 nucleotides.⁶ This RNA genome encodes a polyprotein precursor of approximately 3430 amino acids, which is cleaved into three structural proteins (capsid [C], precursor membrane/membrane [prM/M], and envelope [E]) and seven nonstructural proteins (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5).⁷ Several of these proteins are potential antiviral targets, and various inhibitors have been identified. However, except for the envelope protein and NS3 and NS5, no inhibitors have been identified for the remaining viral proteins.^8,9

Drug discovery is a slow, expensive, and often unsuccessful process. Consequently, drug repurposing (the strategy of finding new therapeutic indications for approved drugs) has gained increasing prominence in modern pharmaceutical research.¹⁰

It is estimated that the total average development cost of a new drug ranges between $2 and $3 billion, with a development timeline of at least 13–15 years. Moreover, only about 10% of the drugs entering Phase I clinical trials ultimately receive approval, with most failures attributable to toxicity or lack of efficacy. A key advantage of drug repurposing is that, for existing drugs, extensive preclinical and clinical data—such as pharmacokinetic, pharmacodynamic, and toxicity profiles—are already available, thereby significantly reducing both development time and risk. Drug repurposing can be achieved through experimental studies or computational methods (in silico drug repurposing), with in silico drug repurposing offering an efficient and cost-effective means of discovering new therapeutic options, particularly for diseases that currently lack effective treatments.¹¹

Given the absence of available treatment for curing WNV infection, it is crucial to identify drugs that can be effectively used against WNV. To this aim we have searched in literature for well-known and extensively used antiviral drugs targeting structural and nonstructural viral proteins of flaviviruses, including WNV. RNA-dependent RNA polymerase (RdRp) of Flaviviruses is the most studied protein and remdesivir, a nucleotide analogue, seems efficient in inhibiting flaviviral RdRps in vitro. Moreover, rilpivirine, an antiretroviral drug used for treating HIV infection belonging to the non-nucleoside reverse transcriptase inhibitors (NNRTIs), has also demonstrated to inhibit the RdRp activity of WNV.^12–14 Finally, doravirine, a recently approved NNRTI for clinical use against HIV, was included in this study due to its structural similarity to rilpivirine and its distinct chemical scaffold, which confers improved pharmacokinetic properties.¹⁵ Although no experimental data is currently available regarding doravirine’s interaction with flaviviral proteins, its involvement aimed at investigating whether other NNRTIs might exhibit comparable or enhanced inhibitory potential against WNV targets.

Therefore, this computational work was designed as a focused, mechanistically driven proof-of-concept analysis rather than a large-scale screening. Our goal was to repurpose three FDA-approved drugs - two used to treat HIV infection (rilpivirine and doravirine) and one approved for SARS-CoV-2 infection (remdesivir) (Figure 1) - to determine which compound most effectively interacts with pivotal sites of WNV proteins. To achieve this, we employed different computational tools in order to analyze and deepen our understanding of key sites described in the literature for different viral targets, examining their interactions with the three antiviral drugs. Once validated, the proposed computational pipeline can be readily extended to high-throughput virtual screening to identify additional candidate compounds.

Figure 1 Two-dimensional (2D) chemical structures of remdesivir, rilpivirine, and doravirine.

Materials and Methods

For structure visualization and image generation, we used Chimera UCSF 1.8 software¹⁶ along with the graphical interface implemented in Maestro v4.1.¹⁷

After selecting the viral protein targets - core protein (PDB ID: 1SFK),¹⁸ NS3 protease (PDB ID: 5IDK),¹⁹ NS5 polymerase (PDB ID: 2HFZ),²⁰ and NS5 methyltransferase (PDB ID: 2OY0)²¹ - a meticulous preparation of these structures was performed. Using the Protein Preparation Wizard v4.1 tool implemented in Maestro operating OPLS-2005 as force field,²² the workflow included removing crystallographic buffer components and water molecules. The preparation then involved adding hydrogen atoms, assigning side-chain protonation states at a physiological pH of 7.4 (± 0.2), and executing an energy minimization simulation.

Remdesivir, rilpivirine and doravirine were downloaded from PubChem²³ (https://pubchem.ncbi.nlm.nih.gov, accessed on 30 May 2025) and prepared through the LigPrep v4.5 Tool.²⁴ Thus, hydrogens were added, and ionization states were calculated using Epik at pH 7.4 (± 0.2). Each structure was submitted to the default energy minimization steps using OPLS_2005 as a force field.²⁵

Binding Sites Evaluation

To assess and analyze the druggability of the binding sites reported in the literature, particularly for targets without co-crystallized ligands (core protein, NS5 polymerase and NS5 methyltransferase in KDKE Motif), we employed the SiteMap program in Maestro.²⁶ This algorithm calculates protein surface regions suitable for ligand binding, providing insights into hydrophobic and hydrophilic areas. Specifically, the program calculates a SiteScore value for each pocket, combining geometric descriptors with physicochemical parameters to enable reliable identification and comparison of potential binding sites. These include: (i) pocket size and enclosure (number of site points/volume), (ii) hydrophobicity, reflecting the balance between polar and nonpolar regions, and (iii) hydrogen-bonding capability, indicating the potential to form favorable donor and acceptor interactions. Thus, this numerical descriptor was crucial for confirming the targetability of the identified sites, as only pockets with a SiteScore greater than 0.8 were considered targetable, thereby validating their suitability for subsequent molecular docking simulations. As reported by Thomas A. Halgren, setting a SiteScore threshold of 0.80 allows SiteMap to effectively discriminate ligand-binding sites from regions not known to bind ligands.²⁷

Molecular Docking Simulations

For each analyzed site of the selected proteins, a rigid receptor grid was generated by centering a 10 × 10×10 Å inner box as follows: for the C-terminal tunnel and the N-terminal hydrophobic pocket of the core protein, the grid was centered on the binding regions identified by the SiteMap tool; for the proteolytic site of NS3, on the co-crystallized ligand; for the Motif A/C region of NS5 polymerase, on three residues (Asp536 in Motif A and Asp668–Asp669 in Motif C); for the SAH-binding site of NS5 methyltransferase, it was centered on the cofactor of the crystal structure, while for the KDKE Motif, on the residues Lys61, Asp146, Lys182, and Glu218. Docking simulations were performed using Glide software v7.8, employing the extra precision (XP) algorithm and generating 10 poses per ligand for each identified site.²⁸

Molecular Dynamics Simulations

All complexes between the viral targets and the three ligands generated in the prior docking simulations underwent 500 ns of MDs using Desmond v4.4,²⁹ with OPLS_2005 as force field.²⁵ Each system was placed in an orthorhombic simulation box in an explicit solvent with TIP4P water model parameters,³⁰ maintaining a 10 Å buffer around the solutes. Counterions were added to achieve charge neutrality for all systems. Thus, following solvation, the systems were subjected to energy minimization and equilibrated using the Martyna–Tobias–Klein (MTK_NPT) isobaric–isothermal ensemble. Equilibration was first performed under NVT conditions, followed by NPT equilibration at a pressure of 1 atm and reaching a temperature of 300 K with the Berendsen thermostat–barostat. The same NPT ensemble was maintained for the subsequent production run. The Simulation Interaction Diagram (SID) tool in Desmond²⁹ was used for trajectory analysis, with trajectory frames recorded every 500 ps intervals to evaluate ligand and protein behavior, analyzing the Root-Mean-Square Deviation (RMSD), calculated on the heavy atoms and backbone, respectively, and the Root-Mean-Square Fluctuation (RMSF). Additionally, the radius of gyration (Rg) was calculated for both apo proteins and doravirine-bound complexes at 500 ps intervals to assess the global conformational stability and overall compactness of each target. The Rg values were obtained over the entire 500 ns trajectory using the Simulation Event Analysis tool in Desmond, defined as the mass-weighted root-mean-square distance of all protein atoms from the molecular center of mass.

MM-GBSA Analysis

Complexes in which the ligand showed good stability during the MDs were subjected to Molecular Mechanics Generalized Born/Surface Area (MM-GBSA) free energy calculations using Prime v6.2.³¹ The calculations were performed at a frequency of one every 10 frames, for a gtotal of 100 frames, employing molecular mechanics and continuum solvation models to determine the binding free energies based on the following equation:

ΔG_bind = G_comp − G_pro −G_lig = ΔE_ele + ΔE_vdw + ΔE_int + ΔE_GB + ΔE_surf

where G_comp, G_pro, and G_lig denote the free energy of the complex, protein, and ligand; by splitting the energy contribution, it refers to ΔE_ele, ΔE_vdw, and ΔE_int as the gas-phase interaction energy between protein and ligand, thus including the electrostatic energy term, the van der Waals energy term, and the bond, angle, and dihedral terms, respectively. On the other hand, ΔE_GB and ΔE_surf indicate the polar and nonpolar desolvation free energy, respectively. The implicit solvation was calculated using the GB model,³² and the non-polar solvation energy was calculated using the solvent accessible surface area algorithm. The ΔG_bind reported in this study omitted the entropy contribution due to its relatively high computational demand and the lack of information on conformational entropy, which could lead to the introduction of additional error into the results.³³

Results

Identification of WNV Protein Targets and Computational Analysis Pipeline

To analyze the behavior of remdesivir, rilpivirine and doravirine towards the structural and non-structural proteins of the WNV, we first conducted a search to identify target proteins with a three-dimensional structure resolved and available in the Protein Data Bank (PDB, https://www.rcsb.org/, accessed on 30 May 2025).³⁴ Thus, five viral proteins were identified: two structural proteins, E (envelope glycoprotein) and C (core protein), and three non-structural proteins, NS1 (immuno-evasive functions), NS3 (protease), and NS5 (RNA-dependent RNA polymerase and methyltransferase). The E protein was excluded from further analysis, since only an open post-fusion monomer conformation is available, without detectable pockets as potential binding sites. Moreover, the non-structural protein NS1 was also excluded from further analyses, as no defined site was identified that could be targeted by small molecules. Based on structural resolution and other crucial features, such as the presence of an inhibitor in complex with NS3, a magnesium atom essential for NS5 polymerase activity,²⁰ and the SAH cofactor for the methyltransferase domain of NS5, the following proteins were selected for further study: core protein (PDB ID 1SFK),¹⁸ NS3 (PDB ID 5IDK) in complex with a capped dipeptide boronate inhibitor,¹⁹ the RNA-dependent RNA polymerase (RdRp) domain of NS5 (PDB ID 2HFZ) with a magnesium ion,²⁰ and the methyltransferase domain of NS5 (PDB ID 2OY0) containing the SAH cofactor (Figure 2).²¹ Thus, for each of the selected targets, we applied the workflow reported in Figure 3. Specifically, after preparing the protein, we energy optimized it, and we identified its plausible binding pockets. Afterwards, we performed molecular recognition studies on the three antiviral drugs, and the obtained complexes were geometrically and thermodynamically analyzed by means of MDs.

Figure 2 X-ray crystallographic structures of the WNV proteins, available and deposited in the Protein Data Bank (PDB), used for computational analysis. The visualization of the structures and the generation of images were performed using the Chimera software.14

Figure 3 Computational workflow applied for each of the identified targets: 1) Protein preparation and minimization; 2) Binding sites identification; 3) Binding site evaluation; 4) Molecular docking simulations; 5) Molecular dynamics simulations; 6) MM-GBSA analysis.

Core Protein

Two pivotal sites have been identified for the structural C protein: a central tunnel, formed by the four C-terminal α-helices of each monomer in the tetramer, containing positively charged residues responsible for the interaction with viral RNA, and a hydrophobic pocket, proposed to anchor the C protein to the membrane, present in each of the two dimers.^18,35 Following preparation and energy minimization of the protein, by using SiteMap,²⁶ we investigated both binding sites and we obtained a SiteScore value ≥ 0.8, indicating that both sites are druggable. Additionally, for the hydrophobic pocket, a Phobic Score over twenty times greater than the Philic Score was measured, confirming its hydrophobic nature, which is essential for interacting with the membrane (Figure 4A) (Table S1).

Figure 4 (A) SiteMap analysis of the C-terminal tunnel and the N-terminal hydrophobic pocket of the core protein, providing quantitative and graphical information to enhance the assessment of the binding sites. For each analyzed site, the SiteScore was reported, with values ≥0.8 indicating a druggable site. Moreover, the hydrophobic and hydrophilic contributions were quantified using the Phobic Score and Philic Score, respectively. Graphically, the protein is shown as a grey cartoon, while hydrophobic, hydrophilic hydrogen bond acceptor and hydrophilic hydrogen bond donor regions are depicted in yellow, red and blue, respectively. (B) 3D structure of the best binding pose generated by molecular docking simulations of the three antiviral drugs with C-terminal tunnel and N-terminal hydrophobic pocket of the structural core protein. Remdesivir is shown in red, rilpivirine in green, and doravirine in purple. The dashed yellow and green lines highlight hydrogen bonds and cation-π interactions, respectively. (C) MDs analysis for the three antiviral drugs in complex with C-terminal tunnel and N-terminal hydrophobic pocket of the core protein. In the top panel, the ligand RMSD plots, calculated on the heavy atoms of the ligand and expressed in Å, are reported. In the bottom panel, the average RMSD values are shown.

Thus, for the two identified sites, we performed molecular docking simulations of the three antiviral drugs, selecting the best binding pose for each compound based on the XP GScore value (Table S2).²⁸ Regarding the C-terminal tunnel, from the best pose generated, only doravirine perfectly arranges at the entrance of the tunnel among the four α-helices, forming two hydrogen bonds with Gln74 and Asn95, respectively, and a cation-π interaction with Lys73. Rilpivirine interacts only with one α-helix, forming a hydrogen bond with Lys73. In contrast, remdesivir interacts with two α-helices, establishing a cation-π interaction with Lys73 and two hydrogen bonds with Val71 and Asn72, respectively. These interactions — particularly the cation-π contacts formed by all three antivirals with Lys73 — may contribute to anchoring the compounds at the entrance of the C-terminal tunnel, potentially hindering access of the viral nucleic acid. On the other hand, within the N-terminal hydrophobic pocket, doravirine and remdesivir form hydrogen bonds with Ala32 and Lys30, respectively, whereas rilpivirine does not show any interaction (Figure 4B). These hydrogen bonds may contribute to maintaining the antivirals within the pocket, thereby hindering direct interactions between the pocket and the membrane.

To further investigate the behaviour of the three antivirals in relation to the two sites of the core protein, the obtained complexes were subjected to 500 ns MDs with explicit water solvent. During the MDs, the stability of the ligand binding mode was assessed by analyzing the Root-Mean-Square Deviation (RMSD) calculated on the heavy atoms of the ligand (all non-hydrogen atoms), with the complex initially aligned to the protein backbone of the first MD frame. At the entrance of the C-terminal tunnel, only doravirine maintains interactions with the four α-helices throughout the whole simulation, remaining at the center of the tunnel’s entrance. However, the RMSD trend exhibits partial instability, despite an average value of 5.24 Å. Remdesivir shows significantly higher RMSD fluctuations, moving out of the tunnel entrance with values exceeding 15 Å, while rilpivirine completely dissociates from the site, reaching RMSD values greater than 30 Å. In contrast, within the N-terminal hydrophobic pocket, Remdesivir displays partial instability but remains within the pocket. The other two drugs, particularly doravirine, demonstrate high stability with minimal fluctuations (Figure 4C) (Figure S1). Moreover, the comparison between the docking poses and the most representative conformations obtained from MD clustering (Figure S2A and B) demonstrated overall conservation of the binding modes.

Since only doravirine demonstrated the ability to maintain its binding at the entry site of the C-terminal tunnel, it was the only compound submitted to the thermodynamic analysis. During the 500 ns MDs, the average ΔG_bind value was −45.02 kcal/mol. For the N-terminal hydrophobic pocket, the thermodynamic contributions of all three antivirals were evaluated. For remdesivir, rilpivirine, and doravirine the average ΔG_bind values during the simulation were −57.55 kcal/mol, −82.92 kcal/mol, and −62.72 kcal/mol, respectively (Table 1).

Table 1 Molecular Mechanics Generalized Born/Surface Area (MM-GBSA) Free Energy Calculations Performed Every 10 Frames, Utilizing Molecular Mechanics and Continuum SolvationModels to Estimate the Binding Free Energies (Δgbind) of the Core Protein Complexes in Which the Antiviral Exhibited Good Stability Throughout the MD Simulations

To investigate the core protein’s dynamic behavior in response to interactions with these specific antivirals, we assessed its stability by analyzing the Root-Mean-Square Deviation (RMSD) calculated on the protein backbone. The complexes were initially aligned to the protein backbone of the first MD frame, allowing for a comparative analysis of RMSD trends between the apo-protein and the complexes throughout the simulation. Regarding the C-terminal tunnel, the apo-protein was compared exclusively with the doravirine-bound complex, as doravirine was the only drug that remained stably bound to this site. The RMSD trends and average values were highly similar, with values of 2.72 Å for the apo form and 2.48 Å for the doravirine-bound complex. The protein complexed with doravirine exhibited slightly enhanced stability, particularly during the mid-phase of the simulation (Figure 5A). This similarity in the behavior was further confirmed by structural superimposition, which revealed no significant conformational differences. Conversely, for the N-terminal hydrophobic pocket, a comparative analysis was conducted across all three complexes. In this case as well, the RMSD trends and average values were comparable, with remdesivir, rilpivirine, and doravirine exhibiting mean RMSD values of 3.00 Å, 2.76 Å, and 3.29 Å, respectively. However, structural superimposition at the end of the MD simulations revealed a key conformational change: all three drugs, upon insertion into the hydrophobic pocket between the two α-helices, induced a widening of this region, rendering the hydrophobic pocket more exposed. In contrast, in the apo-protein, the absence of a ligand and the direct interaction with water molecules caused the two helices to shift closer together, reducing the exposure of the membrane-anchoring pocket (Figure 5B).

Figure 5 Molecular Dynamics analysis assessing protein backbone stability over the course of the simulation, based on Root Mean Square Deviation (RMSD) calculations for the core protein. Stability was evaluated by comparing RMSD profiles of the apo form with those of the protein bound to the analyzed compounds within the C-terminal tunnel (A) and the N-terminal hydrophobic pocket (B). At the bottom of the panel, the structural superimposition of the final simulation frames shows the apo core protein (gray) overlaid with the structures bound to remdesivir (red), rilpivirine (green), and doravirine (purple) in its N-terminal hydrophobic pocket (outlined in dashed circle).

NS3 Protease

The NS3 non-structural protein, in conjunction with its essential cofactor NS2B (NS2B-NS3pro), functions as a serine protease and is crucial for the post-translational cleavage of the viral polyprotein precursor.³⁶ Thus, the site selected for targeting was located near the catalytic serine. Also, the deposited and available models are preferentially in complexes with inhibitors covalently bound to the catalytic serine. In our case, the model used features a capped dipeptide boronate inhibitor. Since the site is well-defined and a co-crystallized inhibitor is present, no pocket analysis was carried out using SiteMap. Consequently, after breaking the covalent bond between the inhibitor and the serine, and following target preparation and minimization, a grid was generated on the inhibitor. From the subsequent molecular docking simulation (Table S3), all three antivirals showed the ability to form a hydrogen bond with the catalytic serine (Ser135). Specifically, remdesivir also can interact via hydrogen bonds with Tyr130, Tyr161, Gly151 and Gly153. Additionally, it forms a cation–π interaction with His51. Rilpivirine, in addition to binding Ser135, establishes hydrogen bonds with His51, Gly151, and Ile155. Finally, doravirine can interact with Ser135, His51, Gly153, and Ile155 through hydrogen bonds. Furthermore, it forms a hydrogen bond with Leu87 of NS2B (Figure 6A).

Figure 6 (A) 3D structure of the best binding pose generated by molecular docking simulation for the three antivirals with proteolytic serine site of non-structural protein NS3. Remdesivir is shown in red, rilpivirine in green, and doravirine in purple. The dashed yellow and green lines highlight hydrogen bonds and cation-π interactions, respectively. (B) MDs analysis for the three antiviral drugs in complex with proteolytic site of NS3. In the top panel, the ligand RMSD plots, calculated on the heavy atoms of the ligand and expressed in Å, are reported. In the bottom panel, the average RMSD values are shown. (C) Distance evaluation between the cyano groups of remdesivir, rilpivirine and doravirine and the serine hydroxyl group in the 500ns of simulation. In the top panel, the cyano group – serine distance plots, expressed in Å, are reported. In the bottom panel, the average distance values are shown. For rilpivirine, which has two cyano groups, these are shown in dark-green and light green, respectively.

To examine the binding features of remdesivir, rilpivirine and doravirine, the obtained complexes were subjected to 500 ns MDs in an explicit water solvent. The stability of the ligand binding was assessed during the simulations by analyzing the RMSD of the heavy atoms in the ligand after aligning the complex to the protein backbone of the initial MD frame. Notably, in this case, all three antivirals lift from the catalytic site and exhibited significant geometric instability, reflected in high mean RMSD values: 9.86 Å for remdesivir, 13.24 Å for doravirine, and over 30 Å for rilpivirine (Figure 6B).

All three drugs contain at least one cyano group, capable of undergoing nucleophilic addition to a triple bond by the serine hydroxyl group, potentially resulting in covalent docking. Thus, an analysis of the distance between the cyano group and the serine hydroxyl group over 500 ns of the simulation was performed, revealing that the cyano groups in all three drugs are positioned too far from the reactive serine residue, effectively precluding the possibility of a covalent docking (Figure 6C).

Due to the high structural instability and low binding affinity to the active site, the complexes were not subjected to thermodynamic evaluations.

NS5 RdRp

For the RdRp portion of NS5, a very intriguing site worthy of further analysis has been identified at the interface between Motif A and Motif C. This site, composed of three aspartate residues (Asp536 in Motif A and Asp668–Asp669 in Motif C) coordinating a magnesium ion, has been demonstrated to play a pivotal role in the catalytic mechanism of RNA polymerase.²⁰ Thus, after preparing the protein and performing energy minimization, we analyzed the identified site by means of SiteMap. Our investigation yielded a SiteScore value ≥ 0.8, thus suggesting the site is druggable (Figure 7A) (Table S4).

Figure 7 (A) SiteMap analysis of the Motif A/C of the NS5 RdRp portion, providing quantitative and graphical information to enhance the assessment of the binding sites. For each analyzed site, the SiteScore was reported, with values ≥0.8 indicating a druggable site. Moreover, the hydrophobic and hydrophilic contributions were quantified using the Phobic Score and Philic Score, respectively. Graphically, the protein is shown as green cartoon, while hydrophobic, hydrophilic hydrogen bond acceptor, hydrophilic hydrogen bond donor, and metal coordination regions are depicted in yellow, red, blue and dark pink, respectively. (B) 3D structure of the best binding pose generated by molecular docking simulation for the three antivirals with Motif A/C of the NS5 RdRp. Remdesivir is shown in red, rilpivirine in green, and doravirine in purple. The dashed yellow and light-blue lines highlight hydrogen bonds and π-π interactions, respectively. (C) MDs analysis for the three antiviral drugs in complex with Motif A/C of RdRp NS5. In the top panel, the ligand RMSD plots, calculated on the heavy atoms of the ligand and expressed in Å, are reported. In the bottom panel, the average RMSD values are shown.

We then performed molecular docking simulations of the three antiviral drugs at the identified site, analyzing the best binding poses (Table S5). Remdesivir showed interactions through three hydrogen bonds with Asp668 and Ser715. Rilpivirine formed two hydrogen bonds with Gln606 and Ser801, respectively, and a π-π stacking interaction with Tyr610. Finally, doravirine established only a single hydrogen bond with Lys694 (Figure 7B).

Also in this case, the obtained complexes between the three antivirals and the RdRp were subjected to 500 ns MDs in explicit water solvent, analyzing the stability of the ligand binding during the simulations by evaluating the RMSD of the ligand’s heavy atoms, with the complex aligned to the protein backbone of the initial MD frame. The results revealed that only remdesivir remained stable at the analyzed site, maintaining a mean RMSD value of 3.31 Å and consistently interacting with the divalent cation throughout the simulation. The persistent interaction with Mg²⁺, along with stable hydrogen bonds with Asp668, Asp669, and Ser715 and π–π interaction with Phe713, may underlie the high stability of the compound within the binding site, potentially altering the protein’s processivity. On the other hand, rilpivirine and doravirine quickly dissociated from Motif A/C, exhibiting no interaction with magnesium (Figure 7C) (Figure S3). As shown in Figure S2C, the alignment between the initial docking pose of remdesivir and the most representative MD clusters indicates that the ligand binding orientations remained largely consistent throughout the simulations.

Considering that Remdesivir was the only compound capable of maintaining its binding at Motif A/C, interacting with divalent cation, it was the only compound submitted to the thermodynamic analysis and, during the 500 ns MDs, the average ΔG_bind value was −17.09 kcal/mol (Table 2).

Table 2 Molecular Mechanics Generalized Born/Surface Area (MM-GBSA) Free Energy Calculations Performed Every 10 Frames, Utilizing Molecular Mechanics and Continuum SolvationModels to Estimate the Binding Free Energies (Δgbind) of the NS5 RdRp Complexes in Which the Antiviral Exhibited Good Stability Throughout the MD Simulations

To examine the dynamic behavior of the RdRp domain within NS5 upon interaction with remdesivir, the only compound demonstrating both stable and high-affinity binding to the targeted pocket, we evaluated the complex’s structural stability using RMSD analyses of the protein backbone, following the same protocol previously applied to the core protein. As shown in Figure 8, both systems exhibited highly similar RMSD patterns, with average deviations of 2.96 Å for the complex and 2.95 Å for the apo state. These findings were confirmed by structural superimposition, which revealed no significant conformational deviations between the two states.

Figure 8 Molecular Dynamics analysis assessing protein backbone stability over the course of the simulation, based on Root Mean Square Deviation (RMSD) calculations for the RdRp NS5 protein. Stability was evaluated by comparing RMSD profiles of the apo form with those of the protein bound to the analyzed compound within the RdRp domain of NS5.

NS5 Methyltransferase

For the NS5 methyltransferase, two crucial sites have been identified: the S-adenosyl-L-homocysteine (SAH) binding site, a product of the cofactor S-adenosylmethionine (SAM) after donating its methyl group, and the KDKE motif (Lys61-Asp146-Lys182-Glu218), which exhibits N-7 and 2′-O methylation activity.²¹ After preparing and energy-minimizing the protein structure, we utilized SiteMap to analyze the KDKE motif. This analysis yielded a SiteScore ≥ 0.8, indicating the presence of a druggable pocket. Notably, similar to the site identified in the RdRp portion, the KDKE motif is predominantly hydrophilic (Figure 9A) (Table S6).

After generating the docking grids for the SAH site and the KDKE motif, we conducted molecular docking simulations with three antiviral drugs (Table S7). At the SAH site, remdesivir demonstrated the ability to form hydrogen bonds with Lys105, Glu111, Asp131, Val132, and Asp146. Additionally, it interacted via π-π stacking with Trp87 and His110. Rilpivirine, on the other hand, formed only two hydrogen bonds, with Gly81 and Arg163. Doravirine was able to establish hydrogen bonds with Ser56, Gly81, and His110 while also forming an π-π stacking interaction with His110.

Figure 9 (A) SiteMap analysis of the KDKE Motif of the NS5 Methyltransferase portion, providing quantitative and graphical information to enhance the assessment of the binding sites. For each analyzed site, the SiteScore was reported, with values ≥0.8 indicating a druggable site. Moreover, the hydrophobic and hydrophilic contributions were quantified using the Phobic Score and Philic Score, respectively. Graphically, the protein is shown as a yellow cartoon, while hydrophobic, hydrophilic hydrogen bond acceptor and hydrophilic hydrogen bond donor regions are depicted in yellow, red and blue, respectively. (B) 3D structure of the best binding pose generated by molecular docking simulations of the three antiviral drugs with SAH site and KDKE Motif of NS5 methyltransferase portion. Remdesivir is shown in red, rilpivirine in green, and doravirine in purple. The dashed yellow and green lines highlight hydrogen bonds and cation-π interactions, respectively. (C) MDs analysis for the three antiviral drugs in complex with SAH site and KDKE Motif of NS5 methyltransferase. In the top panel, the ligand RMSD plots, calculated on the heavy atoms of the ligand and expressed in Å, are reported. In the bottom panel, the average RMSD values are shown.

In contrast, at the KDKE motif, remdesivir and doravirine began only hydrogen bonds with Ser150 and Arg213, and with Lys29, respectively. Instead, rilpivirine can interact via hydrogen bonds with Val35, Arg57, Lys61, and Glu218 and via cation-π interaction with Arg213 (Figure 9B).

In both analyzed binding sites, the hydrophilic nature of the pockets together with the formation of multiple hydrogen bonds could contribute to the stabilization of all three ligands within the MTase region, enabling them to displace the cofactor at the SAH site or the substrate, potentially reducing the enzymatic activity of the domain.

In order to further analyze the behaviour of remdesivir, rilpivirine, and doravirine at the two binding sites of the non-structural methyltransferase, we performed 500 ns MDs of the complexes in explicit water solvent. During the simulations, the stability of the ligand binding mode was evaluated by analyzing the RMSD of the ligand’s heavy atoms. For consistency, the complexes were aligned to the protein backbone of the initial MD frame. For both sites, doravirine exhibited the highest stability, maintaining strong interactions within the binding pockets, resulting in mean RMSD values of 5.26 Å for the SAH site and 4.22 Å for the KDKE Motif. Also, at the SAH site, remdesivir quickly left the cofactor binding site, while rilpivirine demonstrated a promising binding mode during the first 200 ns. Conversely, at the KDKE motif, both remdesivir and rilpivirine dissociated rapidly, showing high mean RMSD values (Figure 9C) (Figure S4). To further assess the consistency between the initial docking orientations and the conformations sampled during MD simulations, a clustering analysis was performed. The alignment between the docking poses and the most populated MD clusters (Figure S2D and E) revealed a general preservation of the starting binding modes.

Concerning the SAH site, both doravirine and rilpivirine were subjected to thermodynamic analysis using the MM-GBSA approach, yielding average ΔG_bind values of −66.48 kcal/mol and −54.16 kcal/mol, respectively. As observed in the RMSD trend, Rilpivirine exhibited a partial rearrangement of its binding mode after approximately 200 ns, while remaining in proximity to the SAH site. Therefore, in addition to the calculation performed over the entire 500 ns trajectory, we computed the average ΔG_bind values separately for the first 200 ns and the last 300 ns of the simulation. The resulting ΔG_bind values were –51.21 kcal/mol and –56.10 kcal/mol, respectively. Instead, regarding the KDKE motif, since only doravirine remains in the pocket, it was the only compound for which its binding mode was thermodynamically investigated, revealing an average ΔG_bind value of −38.88 kcal/mol during the simulation (Table 3).

Table 3 Molecular Mechanics Generalized Born/Surface Area (MM-GBSA) Free Energy Calculations Performed Every 10 Frames, Utilizing Molecular Mechanics and Continuum SolvationModels to Estimate the Binding Free Energies (Δgbind) of the NS5 MTase Complexes in Which the Antiviral Exhibited Good Stability Throughout the MD Simulations

Finally, to study the dynamic behaviour of the MTase domain of NS5 in response to ligand binding, we assessed the structural stability of the complexes through backbone RMSD analysis, using the same protocol. Two binding sites were examined: the SAH site and the KDKE motif. For the SAH site, comparisons were made between the apo-protein and complexes with doravirine and rilpivirine, as remdesivir dissociated from this pocket during the simulation. For the KDKE motif, only the doravirine-bound complex was analyzed, since doravirine was the only compound that remained stably bound at this site. As illustrated in Figure 10A and B, the protein in complex with doravirine showed an RMSD profile comparable to the apo form for both analyzed sites, further supported by structural superimposition, which revealed no significant conformational changes. In contrast, the rilpivirine-bound complex at the SAH site exhibited increased RMSD values, reaching approximately 4 Å by the end of the simulation. Structural superimposition at the final MD frame confirmed this observation, revealing a notable conformational shift, particularly near the SAH binding site and within the Asp36–Val49 loop region (Figure 10A).

Figure 10 Molecular Dynamics analysis assessing protein backbone stability over the course of the simulation, based on Root Mean Square Deviation (RMSD) calculations for the MTase NS5 protein. Stability was evaluated by comparing RMSD profiles of the apo form with those of the protein bound to the analyzed compounds within the SAH site (A) and KDKE motif (B). In panel A (bottom), structural superimposition of the final MD frames shows the apo NS5 MTase (yellow) overlaid with the rilpivirine-bound structure (green), with the cofactor depicted as a purple ball-and-stick. The Asp36–Val49 loop region, exhibiting a conformational shift, is highlighted by a dashed circle.

Doravirine as Potential Multitarget Inhibitor

Thus, considering the promising behavior of doravirine, characterized by high structural stability, favorable thermodynamic results, and consistent interactions with the C-terminal tunnel and N-terminal hydrophobic pocket of the core protein, as well as with the KDKE motif and SAH site of the NS5 MTase domain, we further investigated its geometric stability by performing triplicate MD simulations for all complexes. As reported in Figure 11A, doravirine maintained overall stability and preserved key interactions throughout all simulations. Instead, higher instability was observed when doravirine was bound to the C-terminal tunnel of the core protein. Nevertheless, despite some fluctuations and an average RMSD of approximately 8 Å in the second simulation, the compound remained anchored at the tunnel entrance, maintaining persistent interactions with key residues, as shown in Figure 12.

Figure 11 Triplicate MD simulations of doravirine in complex with the C-terminal tunnel and N-terminal hydrophobic pocket of the core protein, and with the SAH site and KDKE motif of the NS5 methyltransferase. (A) For each analyzed site, the upper panels show the ligand RMSD plots, calculated on the heavy atoms of doravirine and expressed in Å. The lower panel’s report the corresponding average RMSD values. (B) RMSF analysis of the protein residues in complex with doravirine compared to their apo forms. Orange dotted circles highlight the regions exhibiting increased fluctuations corresponding to the specific binding sites.

Figure 12 Protein–ligand interactions of doravirine with the C-terminal tunnel and N-terminal hydrophobic pocket of the core protein (A), and with the SAH site and KDKE motif of the NS5 methyltransferase (B). Only interactions occurring for more than 30% of the 500 ns simulation time are shown. Hydrogen bonds are depicted as purple dashed lines and arrows, π–π interactions as green lines, and cation–π interactions as red lines.

To provide more in-depth computational insight, we also analyzed the root mean square fluctuation (RMSF) profiles of the proteins in complex with doravirine. For the core protein, the most variable regions—whether doravirine was bound to the C-terminal tunnel or the N-terminal hydrophobic pocket—were located within the N-terminal peptides of the hydrophobic pocket. These fluctuations were evident in all chains (E, F, G, and H) and were more pronounced, with higher RMSF values, when doravirine was bound to the N-terminal pocket compared with the apo protein. These fluctuations are highlighted by the orange dotted circles in Figure 11B and support the occurrence of direct interactions between doravirine and the N-terminal hydrophobic pocket.

Regarding the C-terminal tunnel, doravirine induced increased fluctuations in the regions corresponding to two α-helices of the tunnel entrance, highlighted by orange dotted circles in chains F and G. For the NS5 MTase domain, complexation with doravirine also resulted in an overall increase in flexibility compared to the apo form. A common region of increased fluctuations was observed within the Asp36–Val49 loop, consistent with the observed trends in backbone RMSD. In the SAH-binding site, doravirine caused higher fluctuations in three sequence segments involved in cofactor interactions (His110–Arg160), highlighted by orange dotted circles. Finally, when bound to the KDKE motif, doravirine induced a marked RMSF increase in a region located between the ligand and the SAH cofactor, encompassing Asp146 of the KDKE motif (Figure 11B).

These findings further support the potential multitarget and multisite interaction profile of doravirine with these essential WNV proteins. Moreover, doravirine maintained consistent and persistent interactions during the MD simulations in all the analyzed binding sites (Figure 12).

To further investigate the global conformational stability of the selected targets and to complement the local flexibility analysis obtained through RMSF, the radius of gyration (Rg) was calculated for both the apo forms and the doravirine-bound complexes. The Rg provides an estimate of the protein’s overall compactness, enabling the detection of large-scale conformational rearrangements that may occur upon ligand binding.

Across all the simulations, the NS5 MTase domain exhibited relatively stable Rg values in both the apo and doravirine-bound forms, confirming its intrinsic structural stability. Indeed, doravirine binding did not induce significant unfolding or compaction events. These results are consistent with the RMSF analysis, which—apart from the Asp36–Val49 loop—showed only moderate, localized increases in flexibility within regions directly involved in ligand accommodation. Specifically, the Rg values ranged from 18.57 Å for the apo protein up to 19.01 Å for the third simulation of doravirine bound at the SAH site, supporting a globally stable conformation.

In contrast, the core protein, characterized by its tetrameric organization, displayed greater intrinsic structural variability, reflected in a wider Rg range fluctuating by approximately 5 Å across all simulations. This behaviour suggests a naturally more dynamic assembly. Upon doravirine binding, the average Rg values remained comparable to those of the apo form, indicating that the ligand did not compromise the quaternary structure or global stability of the protein. However, higher average Rg values were observed in the second MDs for both the binding sites, consistent with the corresponding RMSD and RMSF trends, which revealed local flexibility. In conclusion, in agreement with the RMSF analysis, the broader Rg range observed for the core protein reflects its intrinsic conformational adaptability, particularly within the N-terminal peptides of the hydrophobic pocket (Figure 13).

Figure 13 Radius of gyration (Rg) profiles of the apo and doravirine-bound complexes at the C-terminal tunnel and N-terminal hydrophobic pocket of the core protein, and at the SAH site and KDKE motif of the NS5 MTase domain over 500 ns of molecular dynamics simulations. Triplicate simulations were analyzed for each target to assess global conformational stability and overall structural compactness.

Discussion

In this computational study, the antiviral drugs remdesivir, rilpivirine and doravirine were evaluated for their binding potential to several structural and non-structural proteins of the West Nile Virus. Four protein functions were selected for our analysis: the core protein (C), NS3 protease, and two domains of NS5 (RNA-dependent RNA polymerase and methyltransferase). An initial structural analysis identified druggable pockets on these proteins, followed by molecular docking and MDs to assess the binding stability and to investigate the thermodynamic properties of the drug-protein complexes. Our results revealed that doravirine exhibited the most promising binding mode, characterized by consistent interactions and high stability at the C-terminal tunnel and N-terminal hydrophobic pocket of core protein and SAH site and KDKE Motif of NS5 MTase. These interactions were maintained throughout extended MD simulations, further supported by RMSF analyses that confirmed local flexibility changes consistent with ligand binding. Remdesivir showed a marked stability with the RdRp catalytic Motif A/C, although the average ΔG_Bind value was not particularly favorable, and moderate stability within the hydrophobic pocket of the c protein. On the other hand, rilpivirine was associated to weaker or less stable interactions overall.

The observed multitarget profile of doravirine aligns with emerging research highlighting the potential of NNRTI-class antivirals to inhibit multiple flaviviral enzymes, including NS5.^12–14 Moreover, previous studies have reported the successful repurposing of clinically approved antivirals against related flaviviruses, supporting the rationale for evaluating doravirine as a novel multitarget candidate against WNV.^37,38

Doravirine, further to better binding profile than remdesivir or rilpivirine, also exhibits superior pharmacokinetics compared to both rilpivirine and remdesivir. Doravirine is rapidly absorbed after oral administration, with bioavailability unaffected by food intake,¹⁵ making it more practical than rilpivirine, which requires food for optimal absorption. Additionally, rilpivirine pharmacokinetics are influenced by gastric pH.³⁹ In contrast, doravirine pharmacokinetics remain stable across various patient factors, including age, gender, race/ethnicity, and both moderate hepatic impairment and mild to severe renal impairment.¹⁵ No dosage adjustments are required in these populations, enhancing its utility in diverse clinical settings. On the other hand, remdesivir, primarily used for COVID-19, is administered intravenously and available as a solution or lyophilized powder for infusion,⁴⁰ limiting its convenience for outpatient use. Thus, overall, doravirine’s promising pharmacokinetic profile, ease of administration, and minimal need for dosage adjustments provide clear advantages over rilpivirine and remdesivir.

The results presented here serve as a mechanistic foundation for subsequent experimental validation within the SPARROW Project (Seeking Preemptive Antiviral Responses and Rapid Diagnostic Tools for West Nile Virus Outbreaks in One Health Approach, Next Generation EU Funding, PNRR: M6/C2, Project Code: PNRR-POC-2023-12377826). Future in vitro assays (eg, viral titer reduction in Vero cells) and in vivo evaluations in murine models will be essential to confirm the inhibitory activity and pharmacodynamic potential of doravirine.

Despite the robustness of the computational approach, some limitations of this study should be acknowledged. First, our analyses are based exclusively on computational methods and therefore depend on the accuracy of the available crystallographic structures and force field parameters. Secondarily, while MD simulations provide valuable insights into binding stability and energetics, they cannot fully account for complex biological factors such as protein conformational heterogeneity, cellular uptake, or metabolic activation. Finally, the MM-GBSA estimations offer relative rather than absolute binding free energies and should be interpreted within a comparative rather than quantitative framework.

Overall, this study provides computational evidence supporting doravirine as a multitarget antiviral candidate with both structural and pharmacokinetic advantages, positioning it as a promising lead compound for future WNV therapeutic development.

Conclusion

WNV infection, already epidemic in several tropical regions, currently it is becoming endemic in new parts of the world including some European countries,⁴¹ and it will likely continue to disperse to naive areas. Although most of WNV infections in humans are asymptomatic, 20% of WN-infected persons show symptoms with fever being the most common clinical presentation. The elderly and immuno-compromised persons are at higher risk of developing West Nile severe clinical forms.⁴² Currently, no specific treatment or vaccine are approved for West Nile disease, therefore it is crucial to intensify research efforts to develop new therapeutic approach.

In this work, an integrated computational approach combining molecular docking, molecular dynamics simulations, and thermodynamic analyses was employed to explore the repurposing potential of three antiviral drugs—remdesivir, rilpivirine, and doravirine—against key WNV proteins. Among these, doravirine demonstrated the most promising multitarget activity, displaying stable and energetically favorable interactions with both structural core protein and non-structural NS5 methyltransferase viral proteins. Moreover, its favorable pharmacokinetics characteristic and the oral administration make doravirine as a potential candidate for early treating infected patients to accelerate viral clearance and to reduce disease severity. Other potential applications include prophylactic or early post-exposure prophylactic treatment for people living in WNV endemic areas, or for travelers visiting countries with ongoing outbreaks. Future efforts will include experimental validation of these computational predictions through in vitro antiviral assays and in vivo efficacy studies in relevant models. In addition, further computational analyses involving larger compound libraries or structure-based optimization of doravirine derivatives could help refine its antiviral potential and guide rational drug design against WNV and related flaviviruses.