Design and Repurposing of Some Nucleoside Analogs Targeting Some Key Proteins of the Avian H5N1 Clade 2.3.4.4b Using Molecular Docking and Simulation

1. Introduction

Influenza virus Type A (IVA) is an RNA, negative sense, and has a segmented genome. The viral genome consists of 7-8 segments [1]. Each viral segment encodes at least one protein. These viral proteins include the hemagglutinin (HA), neuraminidase (NA), the matrix protein (M), the viral nucleoprotein (NP), the viral polymerase called (RNA dependent RNA polymerase (RdRp), which consists of three subunits (PA, PB1, and PB2), as well as the nonstructural protein-1 (NSP-1) [1]. Influenza viruses have a unique feature during their transcription. The viral genome is translocated into the host cell nucleus to acquire the cap from the host cell DNA polymerase to prime the viral transcription [2]. The avian influenza virus (AIV) belongs to the Influenza virus type-A. They have been classified based on their virulence into the low pathogenic (LPAIV) and the high pathogenic (HPAIV) viruses [3]. The influenza viruses are experiencing dynamic changes at the genomic level. There are several reasons for the emergence/reemergence of new viruses, strains, clades, or subclades of avian influenza viruses, including (1) Some point mutations across different proteins called antigenic drift, the exchange of some viral genome segments in case of double infection with more than one type of viruses called antigenic shift, and (3) the poor proofreading capabilities of the viral RdRp [4,5]. The IVA is classified based on their key proteins into H1-H16 and N1-N9 [6]. This pattern of frequent changes in the IVA on the genomic levels. These frequent changes of the viruses on the genomic levels may generate new viral isolates that are resistant to the currently used drugs and vaccines. This is hammering the presently used vaccine and antiviral therapy that is prepared against the previously seeded viruses. Thus, there is an urgent need to monitor these viruses on the genomic levels through active, vigilant molecular surveillance as well as frequent upgrading of the vaccine and drugs that match the most recently circulating strains of the virus. The highly pathogenic avian influenza virus continues to pose a great risk to human health and the poultry industry worldwide. During 2021, the H5N1 clade 2.3.4.4b started to cause several outbreaks in poultry farms across the globe. It spread rapidly to different parts of the world [7]. In early 2024, this clade was reported in cattle and in the milk in many cattle herds across the USA [8]. Human infections due to the same clade were reported in many countries across the world [9]. Unfortunately, there is no current vaccine or even antiviral drugs specific to prevent/treat the infected hosts with this newly emerged clade of the HPAIV. Several antiviral drugs were used in the past to treat many IVA in humans. There are several strategies to use for these antiviral drugs by targeting various key IAV proteins [10,11,12,13]. Most of these strategies are mainly targeting some key viral proteins involved in the viral replication and maturation in addition to some other host proteins that play important roles during the viral replication in the host cells [14]. Here, are focusing on some antiviral strategies that involve the AIV proteins. These strategies mainly focus on the inhibition of the viral entry, release, or the processing and maturation of the viral virions [14]. Some of these strategies mainly depend on the generation of some monoclonal antibodies that mainly neutralize the activities of the viral expressed proteins on the surface of the viral envelope,, particularly the HA and NA [15,16]. In this study, we are mainly focusing on the synthetic antiviral compounds that target various key proteins of the virus. These strategies include (1) inhibitors of the HA maturation, (2) inhibitors of the NA protein, (3) inhibitors of the RdRp, (4) the ion channel (M2) blockers, and (5) inhibitors of the NSP1 [14]. There are some compounds that inhibit the maturation of the HA protein of the IAV, such as Nitazoxanide [17]. The NA protein plays some crucial roles in AIV replication, including the degradation of the mucous layers in the respiratory tracts, thus granting access to the virus to be in touch with the epithelial cells of the respiratory tracts. The transmembrane protein NA also plays a critical role in the removal of the sialic acid from the HA protein. Thus, targeting the NA protein of the AIV is one of the most promising antiviral strategies. Several drugs are currently FDA-approved targeting the NA protein of the AIV, such as (Zanamivir, Oseltamivir, Peramivir, and Laninamivir) [18,19,20]. The viral polymerase RdRp is an attractive target for some antiviral drugs. Several antiviral compounds were shown to inhibit various subunits of this viral polymerase [21,22]. Several compounds targeting the cap-binding domain of the PB2 subdomain (PB2/CBD) have shown some inhibitory effects on the viral replication, such as the Pimodivir [23,24]. The nucleoside analogue antiviral drugs such as Ribavirin showed a broad-spectrum antiviral drug against many RNA viruses, including IAV [25]. The GS-441524, as a nucleoside analog, has been recently used as a powerful inhibitor of the feline infectious peritonitis virus main protease (Mpro) but was never tested as an anti-influenza virus potential drug [26,27]. In the current study, we are comparing the potential impacts of the currently known anti-AIV drugs with some new potential drugs using the molecular docking and simulation approach. We believe these AI tools will pave the way for screening a large number of synthetic and natural compounds to be used as new drugs against many viral diseases, including the currently circulating H5N1 clade 2.3.4.4b. However, further laboratory validations are required to confirm this prediction and to identify the mechanism of actions of these potential drugs.

2. Results

2.1. Homology Modeling of Some Key Proteins of the Avian H5N1 clade 2.3.4.4b (HA, NA, PB2/CBD, and M2)

To perform the homology modeling, we searched for experimentally determined template reference sequences closely related to the query sequences and shared at least 35-65% sequence identity and similarity. The multiple sequence alignment of HA, NA, and PB2/CBD with the templates is shown in (Figure S1-3). The sequences of the target proteins from H5N1 as receptors are HA (804 AA.), NA (797 AA), PB2/CBD (173 AA) of viral RdRp, and four chain M2 (24 AA each chain) channels.

Results of the homology models of the target proteins of the H5N1 clade 2.3.4.4b (HA, NA, PB2, and M2) based on their DOPE scores and energy scores are listed in (Figure S4). We selected the top-ranked models per protein based on their PDF total energy and the lowest DOPE scores (Table 1). Details of the ranked model per protein were selected to predict the 3D model per protein are listed in (Table 1).

The Ramachandran plot results indicate that the majority of the amino acid residues per each protein were located in the most favored regions. The remaining amino acid residues were located in the proper locations. Our results showed high accuracy of all the predicted homology models for the target proteins (HA, NA, PB2, and M2) of the H5N1 clade 2.3.4.4b (Figure 1). Results of both the Ramachandran plot and Verify-3D scores suggest the generated models per each protein are very stable and have the highest scores of the top-quality homology models (Table 1).

2.2. Results of the Binding Affinity of the Antiviral Compounds/Ligands with the Avian H5N1 Clade 2.2.4.4b HA Protein Using the Molecular Docking Approach

Our molecular docking studies revealed a high binding affinity of nearly all the selected compounds. The CDOCKER interaction energy scores suggested binding affinity score (-CDOCK scores) as listed in (Table 2). The GS445214 compound interacted with the active site residues of the same binding pocket as reference standard F0045(S) (Figure 2). The 2D interaction patterns of top hits of F0045(S) (standard) and GS441524 are illustrated in Figure S5.

The highest -CDOCKER interaction energy score (positive value) referred to the most favorable binding of ligand to the protein. The HA protein was docked with GS445214 and standard compound F0045(S), with the -CDOCKER Score of HA- GS441524 complex is 26.62 is higher as compared to the HA complex with a reference standard F0045(S), which is 22.92. The sofosbuvir was not docked in the same binding pocket of HA where the Standard F0045(S) and GS441524 were docked. The standard ligand F0045(S) does not make any hydrogen bond but makes hydrophobic interactions with interacting residues HIS24, TRP366, ILE390, and VAL393 of HA protein. However, the best hit of GS441524 with HA established four hydrogen bonds with residues ALA45, THR331, ILE390, and THR394 in the active site of the target HA protein. The docking results of NA also showed that the CDOCK Scores of standard reference compounds Zanamivir and Oseltamivir with NA are 49.26 and -65.74. However, The -CDOCK scores of sofosbuvir and GS445214 with NA are 53.43 and 32.34 (Table 2). The binding affinity of reference standard Oseltamivir-NA complex is exhibited good as its -CDOCK score is in negative values as compared to sofosbuvir-NA, GS445214, and Reference standard-Zanamivir which are showing positive -CDOCK score (Table 2). However, the binding affinity score (-CDOCK score) of NA-Sofosbuvir is higher than that of the NA-Zanamivir complex. The NA- GS441524 also showed good binding affinity as compared to the Standard Zanamivir-NA docking complex. Due to the interaction of NA-Zanamivir, seven hydrogen bonds were established, and zanamivir interacted with NA binding site residues TYR344, ARG368 (four hydrogen bonds), ASP151, and SER367. However, the top hit of ligand sofosbuvir established eight hydrogen bonds through NA protein interacting residues TYR344, ARG368, TYR402 (two hydrogen bonds), ARG225, LYS432, ASP151, and GLU278 in the active site. On the other hand, GS441524 established 10 hydrogen bonds with interacting residues ARG368 (four hydrogen bonds), VAL149, ASP151 (three hydrogen bonds), SER367, PRO431 (three hydrogen bonds) and LYS432. The GS445214 and Sofosbuvir interacted with active site residues of the same binding pocket as reference standard Zanamivir (Figure 3). The 2D interaction patterns of top hits of Zanamivir (standard), sofosbuvir and GS441524 are illustrated in Figure S6.

Further, the docking of compound GS445214 and Sofosbuvir with PB2/CBD protein of H5N1 also showed promising results as compared to reference standard compound PB2-39 and MGT (CAP) for protein PB2-CBD. The binding affinity scores of sofosbuvir-PB2/CBD and GS441524 docked with PB2/CBD exhibited by -CDOCK score are 50.88 and 34.23. However, the -CDOCK score of reference standard PB2-39 compound and MGT (CAP) for PB2 is 54.55 and 80.67, slightly higher than the binding affinity score of PB2-Sofosbuvir and PB2/CBD-GS441524 docking complex (Table 2). The standard ligand PB2-39 makes only two hydrogen bonds with interacting residues ARG355 and ASN429 of PB2/CBD protein, while MGT (CAP) interacted with residues LYS339, ARG335 (three hydrogen bonds), PHE404 (two hydrogen bonds), GLU361, HIS357 (two hydrogen bonds) and PHE323 established 10 hydrogen bonds with PB2 protein. The interaction of sofosbuvir with PB2/CBD established three hydrogen bonds with residues LYS339, LYS376, and HIS357 in the active site. However, the top hit of GS441524 with PB2/CBD suggested that six hydrogen bonds were established due to interaction of GS441524 with residues ARG355, LYS376, PHE404, GLN406, GLU361 and SER324 in the active site of PB2/CBD protein (Figure 4). The 2D interaction patterns of top hits of MGT and PB2-39 (standard), Sofosbuvir and GS441524 are illustrated in Figure S7. The presence of several hydrogen bonds, along with their favorable interaction score inside the binding cavity of the HA, NA, and PB2/CBD, predicts the validity of the GS441524 as a potential potent inhibitor.

The transmembrane domain of the M2 protein is a tetrameric channel, the target of the anti-influenza drugs Amantadine and its methyl derivative Rimantadine. The M2 protein was docked with Rimantadine and Amantadine as the standard compounds. The CDOCK Score of the M2- Rimantadine complex is 33.02, and that of the M2-amantadine complex is 27.55 (Table 2).

The sofosbuvir showed a good binding affinity with the M2 at the same binding pocket with the -CDOCKER Score of 62.05, which is higher as compared to the M2 complex with the Amantadine and Rimantadine compounds (Table 2).

The GS441524 compound docked with M2 exhibited good binding affinity with the CDOCK Score of 42.59, which is comparatively higher than the binding affinity with the known drugs (amantadien and Ramantadien). However, this binding affinity score is slightly lower than that of sofosbuvir with M2. The M2 interacting with standard ligands and the test ligands established hydrogen bond and hydrophobic bond interactions are listed in (Table 2). The top hit of the sofosbuvir and the GS441524 compounds with the M2 created six and four hydrogen bonds, respectively (Figure 5). The 2D interaction patterns of the top hits of the Rimantadine, the Amantadine, sofosbuvir, and the GS441524 compounds are presented in (Figure S8).

2.3. Results of the Calculation of the Ligands- H5N1 Clade 2.3.4.4b Proteins Binding Energy – (MM-GBSA)

The MM-GBSA calculations for the compounds GS441524 and Sofosbuvir were performed using data from the top 10 docking hits, along with the respective reference standard compounds: F0045(s) for HA, PB2-39 for PB2/CBD, and Zanamivir for the NA protein.

Based in the interaction with the active site of the HA protein, GS441524 exhibited a binding energy of −61.43 kcal/mol, compared to −48.86 kcal/mol for the reference standard. In the case of NA protein interactions, the top docking hits of Sofosbuvir and GS441524 showed binding free energies of −138.23 kcal/mol and −83.10 kcal/mol, respectively. These values were compared to the −201 kcal/mol binding energy of the top hit for the reference compound Zanamivir, indicating the complexes’ relative stability and supporting the docking results’ validity.

The interaction of PB2/CBD with the standard compounds PB2-39 and MGT (CAP) resulted in binding free energies of −236.79 kcal/mol and −236.00 kcal/mol, respectively, which are significantly more negative compared to the binding energies observed for PB2/CBD-Sofosbuvir (−90.00 kcal/mol) and PB2/CBD-GS441524 (−99.00 kcal/mol). These findings suggest greater complex stability for the standard compounds and support the reliability of the docking results (Table 2). For the M2 ion channel protein, the reference compounds Rimantadine and Amantadine showed binding energies of −112.23 kcal/mol and −104.46 kcal/mol, respectively. In contrast, the test compounds Sofosbuvir and GS441524 demonstrated more favorable (i.e., more negative) binding energies of −116.11 kcal/mol and −120.65 kcal/mol, indicating potentially stronger interactions with the M2 protein than the reference compounds.

2.4. Molecular Dynamics Simulation Results for Avian H5N1 Clade 2.2.4.4b (HA, NA, PB2, and M2) Proteins and Some Selected Compounds/Drugs

Proteins play essential structural and functional roles in biological systems, including microbial pathogenesis, by mediating receptor-based internalization and replication. While molecular docking offers insights into the stability of protein-ligand complexes, it lacks the ability to capture their dynamic behavior over time. To address this limitation, molecular dynamics (MD) simulations were performed for 50 nanoseconds (ns) to evaluate the interactions of the GS441524 compound with the HA, NA, and PB2/CBD proteins of avian H5N1 clade 2.2.4.4b. Reference standard complexes—HA-F0045(S), NA-Zanamivir, and PB2/CBD-PB2-39 were also simulated under the same conditions for comparative analysis.

The MD simulations provided detailed insights into the molecular behavior of each complex by capturing time-dependent structural deviations and flexibility. Root Mean Square Deviation (RMSD) was used to assess conformational stability, while Root Mean Square Fluctuation (RMSF) and Radius of Gyration were calculated to evaluate residue-level mobility and compactness of the complexes, respectively. These parameters are critical for understanding the ligand-bound protein systems’ overall dynamics, structural stability, and flexibility. Additionally, hydrogen bond monitoring and hydrogen bond distance analyses were conducted to further explore each complex’s dynamic interactions and stability. The comparative trajectory analysis of GS441524-bound complexes (HA, NA, PB2/CBD) with their respective reference standards is presented and discussed in this section.

2.4.1. Analysis of the RMSD Calculations Conferring Stability of Ligand-Protein Complexes

The RMSD plot analysis showed that all complexes, including the reference compounds, remained relatively stable when RMSD values were within the 3 Å threshold, as lower RMSD values typically indicate greater structural stability of protein-ligand complexes. For the HA-GS441524 complex, however, the RMSD values exhibited significant fluctuations compared to the HA-F0045(S) complex, suggesting lower dynamic stability of the GS441524-bound HA complex during the simulation (Figure 6A). The RMSD of the backbone atoms for HA-GS441524 and HA-F0045(S) ranged from approximately 1.8 Å to 5.5 Å. The average RMSD values were 3.2 Å for HA-GS441524 and 4.0 Å for HA-F0045(S), both exceeding the 3 Å stability threshold (Figure 6A). These results indicate that neither GS441524 nor the reference compound F0045(S) maintained consistent structural stability over the 50 ns simulation period.

The RMSD analysis of the neuraminidase (NA) complex with its standard inhibitor Zanamivir (NA–Zanamivir) revealed minimal structural fluctuations, with an average backbone RMSD of 3.13 Å—slightly exceeding the commonly accepted 3 Å stability threshold. In comparison, the NA–GS441524 complex showed a lower average RMSD of 2.5 Å (Figure 6B), suggesting enhanced structural stability. Thus, over the 50 ns simulation period, GS441524 exhibited superior dynamic stability when bound to the NA relative to the NA–Zanamivir complex.

RMSD analysis of the PB2/CBD complexes with PB2-39 and GS441524 indicated that the PB2/CBD–GS441524 complex exhibited the highest stability, maintaining RMSD values below 3 Å throughout the simulation. A slight increase to approximately 3 Å was noted after 30 ns; however, the RMSD remained steady between 2.7 Å and 3 Å from 30 to 50 ns without significant fluctuations (Figure 6C). In contrast, the reference complex PB2/CBD–F0045(S) displayed an average RMSD around 3 Å but underwent noticeable fluctuations between 17–21 ns and 32–35 ns. These findings suggest that GS441524 forms a more stable complex with PB2/CBD than the reference compound. Overall, the PB2 protein maintained a stable dynamic profile upon ligand binding, supporting the potential of GS441524 as an effective PB2 inhibitor.

2.4.2. Analysis of the Stability and Mobility of the Interaction Between the H5N1 Clade 2.3.4.4b Proteins and the Selected Drugs/Compounds -RMSF

The RMSF (Root Mean Square Fluctuation) plots for the C-α atoms of all complexes were generated from the MD simulation trajectories. In both the HA–F0045(S) and HA–GS441524 complexes, residues within the 337–393 region exhibited the highest flexibility. The HA–F0045(S) complex showed a peak fluctuation of approximately 6 Å, whereas the HA–GS441524 complex displayed slightly reduced fluctuations, peaking around 4.5 Å. Outside this region, most residues in both complexes fluctuated within a 0.5–2.5 Å range throughout the simulation. The overall RMSF profiles of the HA–GS441524 complex were comparable to those of the reference HA–F0045(S) complex, suggesting similar residue mobility and dynamic behavior (Figure 7A).

RMSF analysis of the NA–Zanamivir and NA–GS441524 complexes revealed that residues spanning positions 169 to 270 exhibited the most prominent fluctuations during the simulation. In this region, both complexes reached a maximum fluctuation of approximately 2.5 Å. Outside this range, most residues fluctuated between 0.5 and 1.5 Å over the 50 ns simulation period. The overall RMSF values for the NA–GS441524 complex were slightly lower but remained comparable to those of the NA–Zanamivir complex, suggesting similar residue flexibility and dynamic behavior in both complexes (Figure 7B).

The RMSF analysis of the PB2/CBD-PB2-39 and PB2/CBD-GS441524 complexes revealed that residues between positions 85 and 110 underwent the highest fluctuations during the simulation. The reference complex PB2/CBD-PB2-39 showed a peak fluctuation of approximately 6 Å, whereas the PB2/CBD-GS441524 complex exhibited a slightly lower peak of around 5 Å. Outside this region, residue fluctuations remained within 0.5–2.5 Å over the 50 ns trajectory. Notably, the PB2/CBD-GS441524 complex demonstrated lower overall RMSF values compared to the reference, suggesting reduced residue mobility and improved structural stability. Both complexes exhibited comparable fluctuation patterns around 3 Å within the same residue range, indicating that the PB2 protein preserved structural integrity upon ligand binding. These findings, reinforced by consistently low Radius of Gyration (Rg) values, highlight the compactness and stability of both the reference and test complexes (Figure 7C).

2.4.3. Assessment of the Compactness and Structural Integrity of the H5N1 Clade 2.3.4.4b Proteins-Ligands Complex over Time-Based on the Radius of Gyration (Rg) Using the Molecular Dynamics (MD) Simulations Approach

The average Rg value for the HA-F0045(S) complex was approximately 38 Å, while the HA-GS441524 complex exhibited a slightly lower average Rg of about 37.5 Å (Figure 8A). These similar values indicate that HA-GS441524 maintains a comparable level of compactness to the reference complex. Throughout the simulation, the Rg of HA-GS441524 remained steady around or above 37.5 Å, with fewer fluctuations than the HA-F0045(S) complex, suggesting a relatively stable structural conformation. Although both complexes showed minor fluctuations, the HA-GS441524 complex appeared slightly less compact. Overall, all analyzed complexes demonstrated acceptable levels of compactness, with some conformational changes attributed to ligand binding.

In the case of the NA protein, the average Radius of Gyration (Rg) values for the NA-Zanamivir and NA-GS441524 complexes were 20.1 Å and 20.2 Å, respectively (Figure 8B). These closely aligned values suggest that both complexes preserved structural integrity throughout the simulation, exhibiting only minor fluctuations. Similarly, for the PB2/CBD protein, the PB2/CBD-PB2-39 complex had an average Rg of approximately 16.4 Å, whereas the PB2/CBD-GS441524 complex showed a slightly lower average Rg of 16.2 Å (Figure 8C). The lower Rg value and minimal variation in the PB2/CBD-GS441524 complex indicate a more compact and structurally stable configuration compared to the reference. Overall, the Rg analysis supports that all protein-ligand complexes presented in this study maintained satisfactory compactness during the 50 ns simulation, with no significant conformational changes induced by ligand binding.

2.4.4. Assessment of H5N1 Clade 2.3.4.4b Proteins–Ligands Interactions Based on the Hydrogen Bonds Analysis of the Molecular Dynamics Simulations

The hydrogen bond (HB) formation in the complex during the simulation showed the strength of the complex and dynamic behavior with each bond formed during the time frame. To investigate the stability of the ligand interaction with protein, the measurement of intermolecular hydrogen bond development and hydrogen bond distance between ligand and interacting residues of proteins were evaluated.

The hydrogen bond monitor of complex HA-F0045(S) suggested that around two hydrogen bonds formed throughout the simulation time frame. However, two other hydrogen bonds were also established in various frames during the simulation (Figure 9A). The complex HA-GS441524 formed four intermolecular hydrogen bonds throughout the simulation frames and two hydrogen bonds in various frames during the complete 50 ns simulation time (Figure 9B).

In the case of NA protein hydrogen bond monitoring, the complex NA-zanamivir (standard) maintained four hydrogen bonds throughout the entire simulation time (Figure 9C). However, the complex NA-GS441524 formed three stable hydrogen bonds. (Figure 9D).

The hydrogen bond analysis of the PB2/CBD complex with the standard ligand PB2-39 (PB2/CBD-PB2-39) revealed the formation of two stable hydrogen bonds throughout the 50 ns simulation (Figure 9E), with two additional hydrogen bonds forming intermittently. In contrast, the PB2/CBD-GS441524 complex consistently maintained five stable hydrogen bonds, along with two more that appeared intermittently during the simulation (Figure 9F). These results indicate that the PB2/CBD-GS441524 complex exhibited greater stability compared to the reference PB2/CBD-PB2-39 complex.

2.4.5. Results of the Estimation of the Hydrogen Bond Distances as a Measure of the Complex Stability Between the H5N1 Clade 2.3.4.4b with the Selected Legends

To evaluate the stability of hydrogen bonds, changes in the distances between the compounds and key amino acid residues involved in hydrogen bonding were analyzed throughout the MD simulations. This was done by monitoring the variation in distance between the Cα atoms of the binding pocket and the ligands. For stable hydrogen bonds, the average distance should ideally remain below 3 Å.

In the HA-F0045(S) complex, the standard ligand F0045(S) formed three hydrogen bonds with residues GLN46, THR331, and HIS44 of the HA protein. The average bond distances were approximately 3.0 Å, 3.1 Å, and 5.0 Å, respectively (Figure 10A). However, these interactions showed fluctuations ranging from 2.2 Å to 7.0 Å over the course of the simulation, indicating relatively unstable bonding, particularly in the case of HIS44.

In contrast, the HA-GS441524 complex displayed interactions between GS441524 and four HA residues: THR331, HIS44, GLN46, and HIS465. The average hydrogen bond distances were about 3.2 Å, 3.3 Å, 4.1 Å, and 4.2 Å, respectively, with fluctuations ranging from 2.9 Å to 12.0 Å (Figure 10B). These larger variations further suggest that the interactions between GS441524 and the HA protein were also relatively unstable during the simulation.

For the NA-zanamivir complex (standard ligand), three hydrogen bonds were observed with residues TYR344, ARG368, and LYS432. The average hydrogen bond distances were 2.0 Å, 2.1 Å, and 2.7 Å, respectively, indicating stable interactions between zanamivir and the NA protein (Figure 10C).

In comparison, the NA-GS441524 complex formed hydrogen bonds with four residues: LYS432, ARG368, ARG430, and ARG18. The average bond distances were 2.7 Å, 3.0 Å, 3.2 Å, and 3.5 Å, respectively. These values suggest less stable interactions compared to the standard zanamivir complex (Figure 10D).

For the PB2/CBD-PB2-39 complex (standard ligand), three hydrogen bonds were established with residues ARG355, ASN429, and GLN406, with average bond distances of 3.0 Å, 3.1 Å and 3.3 Å, respectively (Figure 10E).

In contrast, the PB2/CBD-GS441524 complex formed five hydrogen bonds with GLU361, PHE404, LYS367, ARG355, and GLN406. The average hydrogen bond distances were 2.0 Å, 2.2 Å, 2.3 Å, 2.5 Å, and 3.0 Å, respectively, indicating more stable interactions than those seen in the PB2/CBD-PB2-39 complex (Figure 10F).

3. Discussion

This study utilized an integrated computational strategy—including homology modeling, molecular docking, molecular dynamics (MD) simulations, and MM-GBSA binding free energy calculations—to evaluate the structural stability and inhibitory potential of GS441524 and sofosbuvir against key viral proteins of the currently circulating H5N1 clade 2.3.4.4b, namely hemagglutinin (HA), neuraminidase (NA), polymerase basic protein 2 C-terminal domain (PB2/CBD), and matrix protein 2 (M2). These proteins are critical for viral entry, replication, and host adaptation, making them promising therapeutic targets for combating avian influenza [28,29,30,31].

The high-quality homology models of HA, NA, and PB2/CBD were validated using multiple structural assessment parameters, including favorable DOPE scores, low PDF total energy values, and Ramachandran plot analyses. The majority of amino acid residues were located within favored regions of the Ramachandran plots, confirming appropriate backbone dihedral angles and overall model stability. Additionally, the Verify-3D scores for all models exceeded the minimum acceptable thresholds and approached or surpassed expected high-quality benchmarks, further supporting the structural reliability of the predicted protein models [32]. These findings confirmed the reliability of the models for use in subsequent docking and simulation analyses.

Molecular docking results revealed that GS441524 and sofosbuvir exhibit strong binding affinities toward the HA, NA, and PB2/CBD proteins, often matching or surpassing those of established reference ligands (Table 2). Notably, GS441524 showed higher -CDOCKER interaction energy scores than F0045(S) when docked with HA and formed multiple hydrogen bonds with critical residues like ILE390 and THR394, which play key roles in the viral fusion process [33]. Similarly, both sofosbuvir and GS441524 exhibited strong interactions with the NA protein, forming up to 10 hydrogen bonds with key catalytic residues such as ARG368 and ASP151, which are essential for neuraminidase enzymatic function [34].

The PB2-39 is among the identified PB2/CBD inhibitors and has been recognized as a potent agent against the replication of various influenza A virus subtypes. When used in combination with the viral release inhibitor zanamivir, PB2-39 demonstrated a synergistic antiviral effect, enhancing the overall efficacy of the treatment [30].

In the PB2/CBD docking analysis, GS441524 showed extensive interactions with catalytically and structurally significant residues, such as ARG355 and PHE404, forming multiple hydrogen bonds. Its binding energy was comparable to, or slightly lower than, that of the reference compound PB2-39. MM-GBSA calculations supported these findings, revealing favorable binding energies for GS441524 with HA (−61.43 kcal/mol), NA (−83.10 kcal/mol), and PB2/CBD (−99.08 kcal/mol) (Table 2). Notably, the binding energy of GS441524 with NA exceeded that of Zanamivir, a clinically approved neuraminidase inhibitor. [35]. These results highlight GS441524’s potential as a viable anti-influenza agent. Although its binding energies with PB2/CBD were slightly less favorable compared to MGT or PB2-39, the overall interaction profile and energy scores reinforce its promise as a strong candidate for inhibiting influenza virus replication.

The M2 tetrameric channel protein of the influenza A virus serves as a crucial membrane-bound proton channel that acidifies virions enclosed within endosomes, thereby mediating the viral uncoating process and enabling the release of viral RNA into the host cell cytoplasm [36]. In our molecular docking analysis, the binding affinities of standard antiviral agents (amantadine and rimantadine) and test compounds (sofosbuvir and GS441524) were assessed against the M2 transmembrane domain of influenza A virus. Rimantadine and amantadine exhibited moderate binding affinities, with CDOCKER interaction scores of 33.02 and 27.55, respectively, and corresponding binding energies of −112.23 kcal/mol and −104.46 kcal/mol. However, known mutations within the M2 membrane-spanning domain that confer resistance to amantadine have been shown to generate ion channel activity that is no longer inhibited by the drug [37].

The binding sites for amantadine and rimantadine within the M2 ion channel are located near the histidine gate (His37) and pore-lining residues (Leu26, Ala30, and Ser31), where these inhibitors function as channel blockers. Binding interaction analysis showed that both standard and test compounds formed key hydrogen bonds and hydrophobic interactions with conserved residues of the H5N1 clade 2.3.4.4b M2 channel (Table 2). Notably, the test compound sofosbuvir demonstrated a substantially higher CDOCKER score (62.05) and more favorable binding energy (−116.11 kcal/mol) than the standard inhibitors, suggesting a stronger and more stable interaction with the M2 channel. Likewise, GS441524—the parent nucleoside of remdesivir, resulted in a solid CDOCKER score of 42.59 and the most negative binding energy (−120.65 kcal/mol) among all tested compounds. This indicates a potentially higher binding affinity and robust interaction with the M2 protein, particularly through engagement with key residues such as VAL27, ALA30, and HIS37. These findings, closely aligned with known interaction patterns of standard drugs, support the potential of GS441524 as a promising M2 channel inhibitor and affirm the reliability of the docking analysis [38,39,40].

The molecular dynamics simulations conducted over a 50 ns period revealed varying levels of stability across the protein-ligand complexes. In the case of HA, both GS441524 and the reference ligand F0045(S) showed elevated RMSD values exceeding 3 Å, indicating moderate conformational fluctuations and reduced stability. In contrast, the GS441524 complexes with NA and PB2/CBD displayed consistent RMSD values at or below 3 Å throughout the simulation, reflecting strong dynamic stability and limited structural deviations [41,42]. The RMSF analysis revealed minimal fluctuations at the residue level across all complexes, with PB2/CBD-GS441524 exhibiting slightly greater residue stability compared to the reference PB2-39, indicating tighter binding and decreased flexibility. Additionally, the radius of gyration (Rg) measurements confirmed the overall structural compactness of the ligand-protein complexes. Notably, PB2/CBD-GS441524 displayed a more compact and stable conformation than the standard reference, further supporting its enhanced binding affinity and stability [43,44].

The hydrogen bonding analysis of GS441524 with NA demonstrated favorable interaction patterns, although these were slightly less stable compared to those observed with PB2/CBD. GS441524 formed more hydrogen bonds with HA than the reference ligand F0045(S), but the dynamic stability of this complex was comparatively weaker, indicating that further optimization may be needed to improve binding efficacy with both HA and NA. Conversely, hydrogen bond analysis strongly supported the stability of GS441524’s interaction with PB2/CBD. In the PB2/CBD-GS441524 complex, up to five persistent hydrogen bonds were maintained, with bond lengths consistently below 3 Å, reflecting stable and specific binding. This suggests GS441524 is a promising inhibitor of PB2/CBD, potentially outperforming the standard inhibitor PB2-39. These findings are especially important given the critical role of the PB2 protein in viral polymerase function and host adaptation, where its inhibition could significantly disrupt viral replication [45,46].

These computational results are consistent with previous studies that have demonstrated the broad-spectrum antiviral efficacy of GS441524 against various RNA viruses, such as SARS-CoV-2 and feline infectious peritonitis virus (FIPV) [27,47,48,49]. Its ability to target highly conserved functional domains in influenza virus proteins indicates that GS441524 holds promise as a novel therapeutic candidate against avian influenza strains like H5N1, which remain significant zoonotic and pandemic threats. In summary, this study provides important insights into the drug repurposing potential of GS441524 and highlights the effectiveness of computational methods in accelerating the discovery of broad-spectrum antiviral agents against emerging and re-emerging viral pathogens.

4. Materials and Methods

4.1. The H5N1 Clade 2.3.4.4b Sequences Retrieval from the GenBank

In this study, we retrieved the full-length genome sequences of 279 isolates representing the H5N1 clade 2.3.4.4b from the GenBank (NCBI) (https://www.ncbi.nlm.nih.gov/protein) database. The sequence of each protein was downloaded in a fasta format. The sequences of these isolates included (115 from chickens, 40 from ducks, 30 from migratory birds, and 51 from the Canadian geese). We extracted the nucleotides and protein sequences of the (HA, NA, PB2, and M2) proteins out of these sequences.

4.2. Multiple Sequence Alignment (MSA) and Generation of the Consensus Sequences

We performed the multiple protein sequence alignments using Geneiuos software ((https://www.geneious.com/). The consensus sequences per protein were generated based on the individual’s MSA, which was one per protein. The generated consensus sequences will be used for each protein’s downstream prediction and homology modeling.

4.3. The Homology Modeling of the H5N1 Clade 2.3.4.4b (HA, NA, PB2 abd M2) Proteins

We used the MODELLER tool of the Biovia Discovery Studio (v24.1.0.321712) and applied the comparative modeling approach to build the homology model per each target protein. This approach enabled the creation of the 3D model structure for the query protein sequences, as previously described [50]. Briefly, the sequences of the target viral proteins, HA, NA, and the cap-binding domain of PB2 (PB2/CBD) of the viral RNA polymerase complex RdRp were retrieved and downloaded in FASTA format from the NCBI database (https://www.ncbi.nlm.nih.gov/). The homology modeling procedure requires the alignment of the query sequences with template protein sequences from the BLAST search tool. The query protein sequence aligned with 100 most identical template sequences with a maximum similarity of 99% and a minimum of 30% in BLAST with BLOSSUM62 algorithm protocol. Five identical template sequences (similarity: maximum 99-minimum 35%) were chosen to align sequences to the templates. Further, the 3D-build homology model protocol was run to align template structures and predict query sequence 3D confirmation based on template proteins. Discrete Optimized Protein Energy (DOPE) score was used to assess the quality of the generated different 3D protein models to predict five different homology models per query sequence. The DOPE Score and the PDF total energy score are used to evaluate the accuracy of each protein model. The more negative values of the DOPE score and the more positive PDF total energy values, the more stable and accurate the predicted protein model. The best model having the lowest DOPE score was selected for further in silico computational studies as previously described [51]. All the predicted 3D structure models were verified through the Ramachandran plot and the Verify 3D model tool of Biovia Discovery Studio. The Ramachandran plot suggested the stability of structure based on amino acid residues located in the most favored, highly allowed, and allowed region or site of the plot. If the verified score results from the Verify 3D method of the modeled protein is higher than the verified expected low Score value, then the model is considered of acceptable quality. The closer the verified score result is to the verified expected high score value, the better the quality of the predicted protein structure model. The prediction of structure is further visualized using the Biovia Discovery Studio v24.1.0.321712.

4.4. Protein and Ligand Preparation

The prepared proteins, through homology modeling, are used as receptors for ligand interaction. Compounds that are used as ligands in this study are retrieved in 3D-SDF file format from PubChem database: https://pubchem.ncbi.nlm.nih.gov. The compound IDs of all the ligands are listed in (Table 2). For energy minimization of the modeled protein, the smart minimizer at 10000 steps was run to minimize the energy of the prepared protein through the CHARMm forcefield. Further, the preparation of the ligand was done following the energy minimization through fill minimization protocol through smart Minimizer at 10000 steps with CHARMm forcefield using the Biovia Discovery Studio v24.1.0.321712 as previously described [50].

4.5. Molecular Docking of Compounds with Target Proteins

The molecular docking approach can be used to model the interaction between HA, NA, and PB2/CBD protein with some selected and reference standard anti-influenza virus compounds. The interactions between these small molecules and proteins at the atomic level, molecular docking, help us understand how these molecules bind and function within biological processes [52,53]. In the CDOCKER tool of the Biovia Discovery Studio docking protocol, we typically use 10 docking poses per ligand. For each dock/ligand pose, the higher the positive values of the CDOCKER interaction energy score (-CDOCK score), the more strongly the binding affinity is expected [53,54,55].The computational tool of the Biovia Discovery Studio also enables the visualization of the ligand-target interaction (molecular docking) and the identification of the compounds that bind more efficiently with the target [53]. This analysis typically involves examining the docking scores, ligand-protein interactions, and the visualization of the docked complexes. To describe the defining binding site in protein, -CDOCK score through Biovia Discovery Studio v24.1.0.321712. The docking figures were also created with the BIOVIA Discovery Studio, and the docking scores were ranked, with the most positive scores for each complex chosen.

4.6. Calculations of the Molecular Mechanics-Generalize-Born Surface Area (MM-GBSA)

The binding free energy of the protein-ligand complexes in this study was measured using the MM-GBSA approach, which integrates the molecular mechanics (MM) force fields with a generalized Born and surface area continuum with a none-implicit solvation model. The MM-GBSA calculation incorporates the CHARMm force field, with partial charge estimation using Memory Rone. The Di-electric constant was set to be one, and the minimum non-bond higher and lower cutoff distance was set to be 12 and 10 Å. The MM-GBSA was determined in this study using the equation (ΔG = E_complex (minimized)– [E_ligand (minimized) + E_receptor (minimized)) using the Biovia Discovery Studio v24.1.0.321712. The default setting employed to compute MM-GBSA involved rendering all protein atoms rigid while the ligand atoms are relaxed through the Biovia Discovery Studio software.

4.7. Molecular Dynamics (MD) Simulation

Based on in vitro analysis, docking energies, and conformational pose analysis, the top four complexes were selected for MD simulation as previously described [56]. To further explore the dynamic behavior of ligand-protein complexes in this study, the best predicted top hits of compounds and reference standards with respect to HA, NA, and PB2/CBD shown in (Table 2) were selected to perform 50,000 picoseconds (ps) through standard dynamics cascade simulation. To generate the molecular topology files for the protein complex and to create the topology of ligands, the CHARMm36 force field was used. The simulation system consists of an explicit boundary solvent model, orthorhombic box with a minimal distance of 7 nm between the protein surface and the edge of the box neutralized with the inclusion of cation type sodium (Na) and anion type Chloride (Cl) counter ions.

For the energy minimization, the steepest descent (minimization 1) with RMS gradient one and conjugate gradient (minimization 2) with RMS gradient 0.1. Both minimization algorithms were used for 50,00 steps. The heating phase was performed using simulation time 10 ps with time step 2 femtosecond (fs); for immersion, the initial temperature is 50 and the target temperature 300 K with save results interval 2 ps.

The equilibration of complexes was carried out for a 10000 ps simulation run with a two fs (femtosecond) time step, and the save result interval is two ps. The Particle-Mesh-Ewald (PME) algorithm was used for long-range electrostatic interactions with fourth-order cubic interpolation and kappa 0.34 Å grid spacing. The advanced dynamic integrator used Leapfrog Verlet algorithm with applied shake constraint. The implicit solvent model was used with dielectric constant 1, Nonbond list radius cutoff 14 in which nonbond higher cutoff distance is 12 and nonbond lower cutoff distance is 10. The production step of a standard dynamic cascade of MD simulation was carried out for 50 ns (50000 picoseconds). Trajectory analysis was done to confirm the number of hydrogen bonds established per each confirmation, hydrogen bond distance variation of ligand and interacted residues per confirmation, Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and Radius of Gyration (RG) of each system [57]. The stability of the complex is indicated by the highest potential inhibitor from the stable complex protein-ligand through Biovia Discovery Studio software.

5. Conclusions

This comprehensive in silico study provides strong evidence that GS441524, the active metabolite of remdesivir, is a promising multi-target antiviral agent against H5N1 avian influenza virus clade 2.3.4.4v. Using structural modeling, molecular docking, MM-GBSA binding energy calculations, and 50 ns molecular dynamics simulations, GS441524 demonstrated stable and high-affinity binding to four key viral proteins: hemagglutinin (HA), neuraminidase (NA), and the cap-binding domain of PB2 (PB2/CBD). Among these, GS441524 showed especially notable dynamic stability and binding strength with NA and PB2/CBD, supported by favorable -CDOCKER scores, robust hydrogen bonding, consistent RMSD/RMSF values, lower radius of gyration, and shorter hydrogen bond distances. These results suggest GS441524’s potential dual inhibition of viral entry and replication. Consistent with prior reports of its broad-spectrum antiviral activity against multiple RNA viruses, GS441524’s targeting of conserved influenza protein domains positions it as a viable candidate for repurposing against H5N1, a significant zoonotic and pandemic threat. However, in silico findings require experimental validation. Future work should include in vitro antiviral testing, in vivo efficacy and toxicity studies, pharmacokinetics, and exploration of resistance mechanisms. Evaluating combination therapies with existing antivirals like PB2-39, Zanamivir, oseltamivir, or Baloxavir could further advance GS441524’s clinical development.