Unveiling the Role of ZIKV NS4A Mutants F4L and E8D through Molecular Docking and Dynamics Simulation: Implications for MAVS-Mediated Immune Evasion

1. Introduction

Zika virus (ZIKV) has recently become an important pathogen in humans. The first isolation of the Zika virus was made by scientists of the Yellow Fever Research Organization, who removed it from a rhesus monkey in the Zika forests of Uganda (Dick, Kitchen, & Haddow, 1952) Zika virus had a major epidemiological impact in tropical and subtropical countries when the Aedes spp. a “cosmopolitan” vector, spreading widely in tropical areas. History by December 2016, Brazil had documented 440,000–1,300,000 suspected cases and 2,975 cases of ZIK-associated microcephaly. On its own, it had already recorded over 200,000 confirmed cases, leading the WHO to declare PHEIC in February 2016. This outbreak resulted in intense public health responses, including travel warnings and expedited funding of research (Musso & Gubler, 2016) .It is indicative that it is vector-borne and it is a part of the Flavivirus genus within the Potyviridae family, with about 70 viruses including those that are transmitted by mosquitoes or ticks. This family of viruses also includes some of the most threatening human emerging viruses, including DENV, yellow fever virus (YFV), WNV, JEV, and tick-borne encephalitis virus. Zika virus, an Aedes mosquito-borne virus (Boyer, Calvez, Chouin-Carneiro, Diallo, & Failloux, 2018) that can be transmitted through congenital-maternal, blood transfusion, and organ transplant. It replicates in the midgut, salivary glands, and human skin cells (Elong Ngono & Shresta, 2018) and causing relatively mild complaints, such as rash, fever and joint pain. While such symptoms are usually mild, it represents a significant threat to pregnant women, because of triggering microcephaly (Hoen et al., 2018). Zika infection was formerly only diagnosed in research labs, but the use of RT-PCR for detection of viral RNA allows for quite specific identification, albeit limited to the first few days of acute illness (Lanciotti et al., 2008). ZIKV is an enveloped virus with a 10.8-kb positive-sense, single-stranded RNA genome (Russo, Jungmann, & Beltrão-Braga, 2017) (Russo et al., 2017). Its genome contains a single open reading frame coding for a polyprotein of 3,423 amino acids that is cleaved into three structural proteins, capsid (C), membrane (M, generated from the premembrane prM), and envelope (E), and seven nonstructural proteins (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5). The virus particle is composed of the structural proteins as the name suggests. The non-structural proteins support genome replication and packaging, manipulate host pathways to the advantage of the virus (Sirohi & Kuhn, 2017) . Face-length of non-structural protein of the Zika virus (ZIKV) NS4A is 127 AAs (~16 kDa). Localization of the replication complex and the processing of polyprotein in the host cell involve NS4A (Panwar and Singh, 2018). Recent investigations have revealed that of 175 ZIKV isolates studied, 17 isolates harbored six distinct AA replacements in the protein of NS4A. Two simultaneous mutations in NS4A protein (F4L and E8D) existed in both the substitutions. It is also noteworthy that in this isolate no deletion/insertion mutations in the NS4A protein was located in the entire study. The Mitochondrial Antiviral Signaling Protein (MAVS) is a human innate immune effector. MAVS helps in identification of the viral RNA in the infected cells which causes an immunological response against the said virus (Hamel et al., 2015). It is found on the surface of mitochondria and is the essential adaptor protein in the RIG-I/MDA5 dependent signaling induced by viral RNA sensing (Loo & Gale, 2011). When infected with the virus, it is RIG-I and MDA5 that attach to viral RNA in the cytoplasm that causes the proteins to undergo a conformational change. RIG-I and MDA5 activated signaling cascades all meet at MAVS. Their interaction results in MAVS and causes it to produce a filament on the surface of mitochondria. This triggered MAVS therefore serves as a scaffold in further signaling responses leading to activation of transcription factors, such as NF-kB and IRF3 (Nuclear Factor kappa B and Interferon Regulatory Factor 3, respectively) (Otsuka et al., 2005). This is followed by the subsequent NF-kB and IRF3 translocation to the nucleus followed by the activation of the transcription of the genes that are involved in interferon production (Shimizu et al., 2012). Once synthesized, interferon dissociates and binds to other cells in the surrounding that causes antiviral reaction. (Pettersson et al., 2016). This pathway is specifically inhibited by NS4A, which interacts with MAVS (interaction between NS4A and MAVS is impaired) to inhibit downstream transcription factor response. This would lead to the suppression of the normal development of the interferons (Muñoz-Jordán, Sánchez-Burgos, Laurent-Rolle, & García-Sastre, 2003) , so the organism is left unnoticed and the Zika virus stands a greater chance to reproduce and live (IC50). Computational structural biology and molecular modeling Molecular modeling and computational structural biology is an interdisciplinary area of biological research into biomolecular structure using computational methods of bioinformatics, computational chemistry and biophysics to identify the hierarchy in hummock relationship of biomolecules. Homology modeling, molecular docking, MD simulations, and so on, enable us to investigate the dynamics and functional features of biomolecules on an atomic scale (Hollingsworth & Dror, 2018). Molecular modeling in the drug design process speeds up virtual screening, predicting binding affinity and optimization of a compound, resulting in fast drug candidate discovery. (Senior et al., 2020). Here we explore the Zika virus immune evasion tactics with respect to the NS4A protein. Structural studies reveal that the F4L and E8D mutations augment the binding of NS4A to MAVS, interfering with signaling and Interferon production. Molecular dynamics simulations present us with insights into atomic-level changes in these associations, highlighting the importance of computational modeling in understanding virus-host dynamics. Our study findings could be utilized in targeted therapies to interfere with these interactions, preventing the survival of Zika virus.

2. Materials and Methods

2.1. Retrieval of Protein Structures and Sequence

The amino acid sequence of NS4A protein was downloaded in FASTA format from Uniprot https://www.uniprot.org/ (Consortium, 2019). The 3D structure was predicted in Robetta server https://robetta.bakerlab.org) submitting the query sequence in a single-letter amino acid notation format. The server supported both comparative and de novo modeling. The input sequence was refrigeration and domestic hot water (DHW) demand, and it was further treated to predict a domain-specific (or full-length) structure (Bairoch et al., 2005)

2.2. Mutation Recognition

Specifically, ZIKV sequences at the time sets we used to collect ZIKV sequences were obtained from the NCBI database https://www.ncbi.nlm.nih.gov/ , and SNP data were obtained using the CovSurver tool from the GISAID database (https://www.gisaid.org/epiflu-applications/covsurver-mutations-app/ (accessed August 3, 2020)(Barrett et al., 2012). The necessary query sequences were formatted in FASTA and compared with the mutant strains. Also included in this analysis were additional residue positions within NS4A. (Kalia, Saberwal, & Sharma, 2021)

2.3. The Effect of Mutation on Structural Stability

2.3.1. Structure-Based Analysis

The effect of mutations on the stability of NS4A protein was evaluated using graph-based signatures and mCSM server ( http://biosig.unimelb.edu.au/mcsm/stability) Each mutation was calculated with RSA and ΔΔG. To further characterize the effects of these mutations on the protein dynamics and stability, we used an online server DynaMut2 (http://biosig.unimelb.edu.au/dynamut2) . This plug is based on the calculation of the Normal-Mode Analysis (NMA) to compute the i.e., a positive value is stabilizing, whereas a negative value is stimulating (Rodrigues et al., 2021). Combined, these analyses provided a general understanding of the effects of mutations on the structural stability and dynamics of NS4A. 2.3.2 Sequence-Based Analysis Protein stability effects of amino acid substitutions of NS4A were measured by winning pH and temperature configures at I-Mutant server (http://folding.biofold.org/i-mutant/i-mutant2.0.html) with default setting (ΔыG). The server requires input of mutant protein sequence and positions of wild-type (WT) residues. A positive value of ΔG 0 is stability, whereas a negative value of ΔG 0 is decreased stability (Calabrese et al., 2009).

2.4. Variant Modeling and Superimposition

To simulate the F4L and E8D mutations, the wild type (WT) NS4A structure was minimized by using the MOE Workstation. The in-silico mutants were compared with WT structure using PyMOL, where the mutants were superimposed with WT and the RMSD values were computed which would provide the estimate of differences regarding structures.

2.5. Protein–Protein Docking

Docking analysis of WT NS4A, its mutants and MAVS protein were performed by HDOCK server (Wu et al., 2017). Every docking run took ~ 30 minute and the top ten models were chosen considering empirical potential, docking score, Van der Waals interactions and ligand RMSD. A maximum of 100 models per complex were generated by the HDOCK server, and the lowest-energy conformation was used for each complex. Details of the interactions (nonbonded contacts, salt bridges, and hydrogen bonds) were calculated with PDBsum http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.The model that was produced by each docking run with the minimum energy score was taken to be the best (Muhammad et al., 2023).

2.6. Molecular Dynamics Simulation

The stability of the NS4A complexes was tested by performing molecular dynamics (MD) simulations with Amber20 using the Amber force field (Salomon-Ferrer, Gotz, Poole, Le Grand, & Walker, 2013). The systems were solvated with water TIP3P molecules and counter-ions were used to neutralize the systems. Simulations proceeded as follows in stages: energy minimization, heating, equilibration, and production. Energy refinement was performed for 9000 steps which consisted of 600 steps of the steepest decent (Watowich & Morimoto, 1988) followed by 300 steps of the conjugate gradient minimization. The system was then equilibrated at 300 K and 1 atm before 100 ns long production MD. Long-range electrostatics were handled by particle mesh Ewald (PME) (Salomon-Ferrer et al., 2013), and covalent bonds were constrained by the SHAKE algorithm (Kräutler, Van Gunsteren, & Hünenberger, 2001). Trajectory and MD analysis were performed using Amber20 CPPTRAJ package and PMEMD. CUDA, Tyche, and respectively (Roe & Cheatham III, 2013).

3. Results and Discussion

The Zika virus has mutated a number of times over the years to become more pathogenic and more contagious. The majority of these mutations have been found in both structural and non-structural proteins, impacting viral infectivity, disease severity, and patient outcomes (Elong Ngono & Shresta, 2018). Multiple works have shown that mutations located in the non-structural protein NS4A are the key determinants in escaping the human immune response. In particular, NS4A inhibits interferon signaling and generation by directly binding to MAVS, and hence inhibiting interferon response. Given that there are diverse mutations at NS4A identified, it is very important to clarify how each of the mutations affects the physical interaction of NS4A with MAVS.

Figure 1. The overall work flow, including step-wise approaches used in the research.

Figure 1. The overall work flow, including step-wise approaches used in the research.

Thus, in this study provides significant contributions to our understanding of the molecular interaction between NS4A and MAVS and the roles of the newly identified mutations contributed to immune escape. addition, these results may also provide a basis for designing new therapies targeting NS4A to promote the recovery of host immune functions. The general process and the stepwise operators used in this work is presented in Fig 1.

3.1. Identification of New Mutations in NS4A Protein

The GISAID online database was employed to search for the recently reported NS4A protein (AOA14215B9) sequence from UniProt database and for the examination of newly reported mutations. By contrasting the sequence to the wild type a few new mutations were found. The new strain has four amino acid substitutions in the NS4A region including L48M, E8D, K42E and F4L (Fig. 2 and Table 1).

3.2. Impact of Mutation on the Structural Stability of NS4A Protein.

Protein stability is a key factor that determines the structure, function, and regulation of proteins (Pace & Tanford, 1968).Hence, hereby various in silico tools like DynaMut2, mCSM and I- Mutant 2.0 were used to predict the effect of mutations on structural and functional stability of NS4A protein. I-Mutant 2.0 server calculations for four mutants generated ∆∆G values from 0.63 kcal/mol to-3.8 kcal/mol. Of these, the E8D mutation was identified to destabilize the protein (Table 2).

Figure 2. Schematic representation of mutation identified in NS4A protein.

Figure 2. Schematic representation of mutation identified in NS4A protein.

Additional analysis was performed with the DynaMut2 server in order to find destabilising variants. The ∆∆G values obtained in this analysis ranged from 0.29 kcal/mol to -0.19 kcal/mol. Among these four mutations, F4L and E8D showed destabilizing and K42E and L48M showed stabilizing effects on NS4A and hence are responsible for increasing protein stability in general (Table 2).

Using DynaMut2 and I-Mutant 2.0, dozens of destabilizing mutations were uncovered and partially confirmed by the mCSM server. In particular, substitutions including F4L (∆∆G = -0.289 kcal/mol), L48M (∆∆G = -0.604kcal/mol) and E8D (∆∆G = -0.462 kcal/mol) definitely destabilize the NS4A protein and influence its structural stability. The K42E mutated variant (∆∆G = 0.287 kcal/mol) was however categorized as stabilizing. For mapping the importance of highly destabilizing mutations (F4L and E8D) in immune evasion, an extra analysis was further conducted to study their influence on the network of NS4A–MAVS binding (Table 3).

3.4. Variants Modeling of NS4A Protein and Its Superimposition on WT NS4A

3D structure of a protein dictates its function and interaction with other molecules within the body458–460. In order to predict the 3D structure of NS4A protein, the amino acid sequence of NS4A was submitted to the Robetta server (https://robetta.bakerlab.org/) (Figure 3A). The server produced five models based on the input sequence. The model was further validated using ProSa-Web and PDBsum. The predicted structures were analyzed using Ramachandran plot and the model with the most residues in the favored regions and least outliers was chosen (Figure 3B). This ensured that the selected model possessed a biologically meaningful and structurally stable conformation. Furthermore, the quality of the best model was further evaluated, and structural discrepancies were revealed using the ProSa-Web tool. In order to examine the impact of destabilizing mutations (F4L and E8D) on NS4A–MAVS binding the mutations were installed in WT NS4A protein using Chimera program. The models of the mutant and WT proteins were then overlaid and the RMSD values were determined based on the structural differences between them (Fig. 4). Introducing these mutations induced changes in the protein secondary structure and conformation, and the effects of these mutations on the binding affinity between NS4A and MAVS

Figure 3. Structural modeling of NS4A showing mutant FL4 (A).And mutant E8D (B) while (C) show wild type.

Figure 3. Structural modeling of NS4A showing mutant FL4 (A).And mutant E8D (B) while (C) show wild type.

Figure 4. Superimposition of WT NS4A with mutants showing RMSD values for F4L (A, B) and E8D (C).

Figure 4. Superimposition of WT NS4A with mutants showing RMSD values for F4L (A, B) and E8D (C).

3.5. Bonding Network Analysis of WT NS4A-MAVS and Mutants NS4A-MAVS

The protein-protein interaction (PPI) is a problem where molecular docking methods have offered a helpful approach to the study of the structure and functionality of PPIs in relation to disease mechanisms. The mode of binding and protein conformation of known or putative associations of the proteins between the PPI is predicted by molecular docking, which plays a significant role in the molecular pathway of disease development. The critical role of NS4A in immune evasion is that it binds N-terminal caspase activation and recruitment domain (CARD) of MAVS. This event sequesters MAVS away from RLRs, prevents IRF3 activation, retains mRNA in the nucleus and inhibits interferon (IFN) production. Because NS4A is important for immune evasion, and MAVS plays a role in IFN regulation, the WT and mutants (F4L and E8D) were analyzed for their interactions. HDOCK was used for docking of MAVS–NS4A WT, and the details of interactions were analyzed with PDBsum. The FF12MC hydrogen bond analysis show the MAVS–NS4A WT complex has 195 non-bonded contacts and 2 hydrogen bonds was formed with Glu439–Thr63 and Tyr399–Leu28 having hydrogen bond, and salt bridge between Glu449–Arg33 in the MAVS–NS4A WT complex (Figure 5A). These comparisons indicate that the binding is even slightly better in the case of docking of MAVS–F4L with the HDOCK server (188 non-bonded contacts and 4 hydrogen bonds). In the MAVS–F4L complex, hydrogen bond were observed for residues Cys402–Trp93, Arg446–Gly1 (twice) and Tyr440–Met56, suggesting that containing mutation, F4L enhanced academy of NS4A by MAVS (Figure 5B)

Figure 5. Docking of NS4A WT and F4L mutant complexes with MAVS, where (A) shows key hydrogen bonding interactions of WT NS4A with MAVS (stick representation left, 2D interaction pattern right) and (B) shows key hydrogen bonding interactions of F4L mutant with MAVS (stick representation left, 2D interaction pattern right), with non-bonded contacts, salt bridges, and hydrogen bonds illustrated.

Figure 5. Docking of NS4A WT and F4L mutant complexes with MAVS, where (A) shows key hydrogen bonding interactions of WT NS4A with MAVS (stick representation left, 2D interaction pattern right) and (B) shows key hydrogen bonding interactions of F4L mutant with MAVS (stick representation left, 2D interaction pattern right), with non-bonded contacts, salt bridges, and hydrogen bonds illustrated.

Figure 6. Docking of E8D mutant complex with MAVS. Stick representation (left) shows key hydrogen bonding interactions, while the 2D diagram (right) highlights interaction patterns. Non-bonded contacts, salt bridges, and hydrogen bonds are indicated.

Figure 6. Docking of E8D mutant complex with MAVS. Stick representation (left) shows key hydrogen bonding interactions, while the 2D diagram (right) highlights interaction patterns. Non-bonded contacts, salt bridges, and hydrogen bonds are indicated.

Docking analysis using the MAVS–E8D complex showed 180 non-bonded contacts and 5 hydrogen bonds. The hydrogen bonds implicated residues Cys402–Trp93, Arg446–Gly1 (twice) and Tyr440–Met56 and Tyr399–Ile25. These findings indicate that the E8D mutation has strengthened the binding of NS4A to MAVS and can also be causing its augmented functioning to immune evasion (Fig 6). The outcome of molecular docking experiments showed an increase in binding affinity of NS4A mutants with MAVS as opposed to the WT and other mutations. In particular, F4L variant exhibited higher binding, which may enhance the NS4A capacity to avoid the host immune response. It is worth noting that the docking outcomes revealed the maximum binding affinity of E8D to MAVS especially in regard to the number of hydrogen bonds.

3.6. Dynamic Stability Analysis of Wild-Type and Mutant (F4L and E8D) NS4A–MAVS Complexes

Root Mean Square Deviation (RMSD) is a measurement of the average distance between atoms of aligned protein structures. It is also common in the field of molecular dynamics to track changes in conformational stability and deviations. A lower RMSD value implies stability of the structure and a higher value implies change or flexibility. It is necessary to study protein folding, the binding of ligands, and biomolecular dynamics in general calculated as the square root of mean squared deviations.(Grant, Rodrigues, ElSawy, McCammon, & Caves, 2006).The RMSD analysis of 100 ns molecular dynamics simulation revealed the different stability patterns among the NS4A-MAVS complexes. The wild-type complex (Figure a) was stable during the trajectory with RMSD values oscillating around 23A with no dramatic increase and indicating good conformational stability. By contrast, the F4L complex (Figure b) showed a slower variation in RMSD with time, starting off near 2A at the initial stage (0-20 ns) and gradually increasing to the range of 45A by 80-100 ns, which implies a moderate degree of structural drift and lower stability compared to the wild type. The E8D complex appeared by far the most unstable (Figure c), and the RMSD values increased rapidly within the first 20ns, between below 2A and approximately 5-6 A and finally, stabilized at this higher level throughout the remainder of the simulation, indicating significant conformational changes and loss of compactness. Collectively, these findings indicate that whereas the wild-type NS4A-MAVS complex remains stable in its respective structural dynamics during the duration of the 100 ns simulation, the F4L variant undergoes destabilization over time, and the E8D variant experiences an early and sustained destabilization, both of which could disrupt the integrity of interactions between NS4A-MAVS and potentially contribute to the signaling of host immunity by NS4A.

Figure 7. RMSD analysis of NS4A–MAVS complexes over 100 ns simulation. (a) RMSD profile of WT NS4A–MAVS complex, (b) RMSD profile of NS4A F4L–MAVS complex, and (c) RMSD profile of NS4A E8D–MAVS complex, each showing atomic deviations with respect to simulation time.

Figure 7. RMSD analysis of NS4A–MAVS complexes over 100 ns simulation. (a) RMSD profile of WT NS4A–MAVS complex, (b) RMSD profile of NS4A F4L–MAVS complex, and (c) RMSD profile of NS4A E8D–MAVS complex, each showing atomic deviations with respect to simulation time.

3.7. Residue-Level Fluctuation Analysis of Wild-Type and Mutant (F4L and E8D) NS4A–MAVS Complexes

The Root Mean Square Fluctuation (RMSF) provides residue-level information about the dynamic flexibility of proteins during molecular dynamics simulations. Unlike RMSD, which captures overall structural deviation, RMSF quantifies the average displacement of each residue from its mean position over the trajectory. Lower RMSF values are indicative of stable, less flexible regions such as structured domains, while higher values highlight flexible loops, terminal segments, or regions prone to conformational rearrangements. This makes RMSF particularly useful for identifying mutation-induced changes in local mobility that may influence protein stability and interaction with binding partners(Lindahl, Hess, & Van Der Spoel, 2001). When dealing with the NS4A–MAVS complexes, the wild-type protein had relatively low RMSF values in the majority of the individual residues, mostly in the 130-300 range, with only relatively small peaks at residues 40-60 and 170-190. This trend suggests that there is a lack of flexibility in the backbone and that the wild-type complex is in the compact and stable conformation throughout the simulation. F4L variant experienced a slightly higher fluctuation profile than the wild type, peaks went up to approximately 4A, especially around residues around ~30-50 and in the end (at around 200-220). These gains suggest local moderate flexibility, but not a global destabilization of the structure, suggesting only a partial compromise in stability. The highest peak was concentrated around residue ~100 in which the RMSF rose to a peak of about 8–9A, which was much larger than both the WT and F4L complex. Further increases were observed at the C-terminal end (~200-220) with values of about 5-6A. Such extreme fluctuations demonstrate that the wild-type NS4A-MAVS complex is characterized by stable residue-level dynamics and indicate that the F4L mutation causes localized increases in flexibility, whereas the E8D mutation causes overall instability, especially near residue ~100. These variations suggest that as E8D greatly destabilizes the NS4AMAVS interaction, F4L leads to only a minor perturbation that can affect the binding of MAVS, which inhibits the immune modulatory activity of NS4A in the pathogenesis of Zika virus.

Figure 8. Fluctuation analysis was performed for each residue in both the wild-type and mutant NS4A-MAVS complexes.

Figure 8. Fluctuation analysis was performed for each residue in both the wild-type and mutant NS4A-MAVS complexes.

3.8. Compactness Analysis of Wild-Type and Mutant (F4L and E8D) NS4A–MAVS Complexes

The Radius of Gyration (Rg) is a measure of structural compactness and provides insights into the overall folding and stability of protein complexes during molecular dynamics simulations. A stable and lower Rg value generally indicates that the protein remains tightly packed, while higher or fluctuating values suggest structural loosening, expansion, or unfolding tendencies. For the wild-type NS4A–MAVS complex (Figure a), Rg values were relatively stable, ranging between ~18.8 and 19.3 Å over the 100 ns trajectory. Only minor fluctuations were observed, indicating that the wild-type maintained a compact and well-folded structure throughout the simulation, reflecting its inherent structural stability(Lobanov, Bogatyreva, & Galzitskaya, 2008). In the case of the wild-type NS4A-MAVS complex (Figure a), the Rg value was relatively constant (between =18.8 and 19.3 A) throughout the 100ns trajectory, indicating that the Minor oscillations were revealed only, which indicated that the Rg profile of the F4L variant (Figure b) was more variable than the one of the wild type. The complex started forming around ~18.8 A yet with significant oscillations with a peak above 20.0 A especially in the initial stages of the simulation (020ns). It is associated with this behavior of higher and more unstable values of Rg, indicating reduced compactness and implying that the F4L mutation causes a moderate destabilization, presumably due to small conformational reorganizations.

Figure c showed that the E8D variant had an intermediate pattern between the wild type and F4L. Rg values were also consistent with the range of values measured between 18.6 and 19.2, however with a reduced range of values at the entire 100 ns trajectory. Although the system was not markedly expanded as in F4L, it also was not as tightly compact as the wild type was. Such sustained oscillations are indicative of reduced stability and inclination to structural looseness as compared to the WT.

Figure 9. Radius of gyration (Rg) analysis of NS4A–MAVS complexes over 100 ns simulation. (a) Rg profile of WT NS4A–MAVS, (b) Rg profile of NS4A F4L–MAVS, and (c) Rg profile of NS4A E8D–MAVS, showing changes in molecular compactness across the trajectory.

Figure 9. Radius of gyration (Rg) analysis of NS4A–MAVS complexes over 100 ns simulation. (a) Rg profile of WT NS4A–MAVS, (b) Rg profile of NS4A F4L–MAVS, and (c) Rg profile of NS4A E8D–MAVS, showing changes in molecular compactness across the trajectory.

3.8.1. Post-Simulation Hydrogen Bonding Analysis of Wild-Type and Mutant (F4L and E8D) NS4A–MAVS Complexes

Hydrogen bond are critical stabilizing forces that maintain protein–protein interactions during molecular dynamics simulations

(Kumar, Ma, Tsai, Sinha, & Nussinov, 2000). The H-bonds in the wild-type NS4A68-MAVS complex (Figure a) did not change but generally fluctuated between approximately 90 and 115 over the 100 ns trajectory, indicating that there is good intermolecular stability. A very similar profile was observed in the F4L variant (Figure b), which retained between 85 and 115 H-bonds with occasional values higher than 120, meaning that despite the moderate conformational drift seen in RMSD and Rg, the hydrogen-bonding network remained intact in the F4L variant. On the other hand, the number of H-bonds on the E8D variant (Figure c), was lower with a smaller deviation, which was more likely to be in the range, they were stronger and weaker, with weaker stabilizing interactions.

Figure 10. depicts the average hydrogen bonds in the wild-type NS4A-MAVS complex (a), F4L-MAVS complex (b), or E8D-MAVS complex (c). The average hydrogen bond analysis of the NS4A-MAVS complex.

Figure 10. depicts the average hydrogen bonds in the wild-type NS4A-MAVS complex (a), F4L-MAVS complex (b), or E8D-MAVS complex (c). The average hydrogen bond analysis of the NS4A-MAVS complex.

3.8.2. Binding Free Energy Calculations

The stability and favorability of the formation of the complexes between the NS4A and MAVS were estimated by the binding free energy of the complex using the MMPBSA method. The wild-type complex had an overall free energy (ΔG Total) of -61.73 kcal/mol indicating a stable association between NS4A and MAVS. The F4L mutant was relatively less likely to bind with an 8G Total of -54.05 kcal/mol, which would be in line with moderate destabilization. Interestingly, the E8D mutant showed a 0 Total of -56.25 kcal/mol, which though marginally stronger than F4L, nevertheless, showed reduced binding as compared to the wild type. Further energy decomposition showed that van der Waals (ΔEvdw) interactions dominated the other stabilizing factors in all the complexes with values of -108.30, -106.51 and -104.73 kcal/mol corresponding to the wild type, F4L and E8D, respectively. Electrostatic contributions (Eele) also played an important role, especially in the E8D variant (–112.31 kcal/mol), but were mostly compensated by the increased polar solvation penalty (EGB). Eventually, the binding profiles depended upon the balance between gas-phase energies (/170.00,/191.41,/217.04 kcal/mol) and solvation free energies (/108.27,/137.36,/160.79 kcal/mol). Overall, these findings indicate that the wild-type NS4A-MAVS complex is the one with the most preferable binding affinity, although both F4L and E8D forms are less stable, with F4L being the weakest binder. The destabilization of these mutants could affect the interaction of MAVS and consequently alter the immune modulation of NS4A in relation to Zika virus infection.

4. Conclusions

This study provides structural and energetic insights into how ZIKV NS4A mutations F4L and E8D enhance viral immune evasion. Although stability analyses revealed some destabilizing effects on NS4A structure, molecular docking and simulation consistently showed stronger affinities of both mutants toward MAVS compared to the wild type. The E8D variant emerged as the most impactful, forming a denser hydrogen bonding network and demonstrating a total binding free energy of –56.25 kcal/mol, surpassing the wild type (–61.73 kcal/mol). These results indicate that while F4L contributes to enhanced MAVS engagement, E8D may provide a greater evolutionary advantage by reinforcing NS4A’s ability to suppress interferon signaling. By integrating molecular dynamics and MM/PBSA analyses, this work highlights the role of NS4A mutations in shaping ZIKV immune evasion and underscores their potential as molecular targets for antiviral intervention.