N1 neuraminidase of H5N1 avian influenza A virus complexed with sialic acid and zanamivir – A study by molecular docking and molecular dynamics simulation

Abstract Development of antiviral drugs is an urgent need to control and prevent the presently circulating H5N1 avian influenza virus which is affects the human respiratory tract. The complex crystal structure of N1–N–acetylneuranamic acid (sialic acid, SIA) is not available as complex and hence SIA and zanamivir (ZMR) are docked into the binding site of N1 neuraminidase. Based on the analysis, the initial complex structures have been simulated for 120 ns to get insight into the binding modes and interaction between protein–ligand complex systems. NAMD pair interaction energy and MM–PBSA binding free energy are calculated and show that there are two possible binding modes (BM1 and BM2) for N1–SIA and a single binding mode (BM1) for and N1–ZMR complex structures respectively. BM1 of N1–SIA is the most preferred binding mode. On contrary to the currently available drugs in which the chair conformation is distorted, in both the binding modes of N1–SIA, the binding pocket of N1 neuraminidase is able to accommodate SIA in 2C5 chair conformation which is the preferred conformation of SIA in solution state. In N1–ZMR complex, ZMR is bind in a distorted chair conformation. The neuraminidase binding pocket is also able to accommodate galactose of SIAα(2→3)GAL and SIAα(2→6)GAL. RMSD, RMSF and hydrogen bonding analyses have been carried out to identify the conformational flexibility and structural stability of each complex system. All the analyses show that SIA can be used as an inhibitor for N1 neuraminidase of H5N1 influenza viral infection. Communicated by Ramaswamy H. Sarma


Introduction
Carbohydrates are an important class of biomolecules, because of their indispensable role in molecular recognition as the interaction between carbohydrates and any other biological macromolecule (molecular recognition of carbohydrates) is the first step to start a biochemical reaction or biological process (Palomares et al., 2011). These interactions play a dominant role in biological processes such as cell adhesion, cell trafficking, signal transduction, egg-sperm interaction bacterial, cell-matrix interaction, fungi and parasite adhesion (Varki, 2017). Cell surface carbohydrates act as the receptor for the attachment of parasites, fungi, toxin, bacteria, and viruses (Doyle, 2011). Besides that, the oligosaccharides chains which terminate with Neu5Aca(2!3)Gal or Neu5Aca(2!6)Gal moiety act as the receptor for various types of viral infections such as influenza viruses, equine viruses and feline calci viruses (Priyadarzini et al., 2012a;Selvin et al., 2012;Stuart & Brown, 2007). The Influenza A virus pathogen is an eight segmented encoding 10-17 proteins, negative-sense, single-stranded RNA genome (Gallagher et al., 2017;Rajao et al., 2018;Ray & Banerjee, 2018). The viral particles are generally spherical in shape and their diameter is % 80-120 nm (Modena et al., 2007;Priyadarzini et al., 2012b). The cell membrane glycoproteins hemagglutinin (HA) and neuraminidase (NA) form the spikes on the surface virion with the spike ratio of 4:1 (Gaymard et al., 2016). The rod-shaped HA spikes are responsible for the attachment and penetration of viral particles into the host cell membrane containing sialic acid (SIA, Neu5Ac) receptor and the mushroom-shaped NA spikes are responsible for cleavage of SIA and release of viral particles into the host cell body (Das et al., 2009;Du et al., 2018).
In the process of development of antiviral drugs to influenza targeting to block the viral attachments to the host cell is an important step. Hemagglutinins and neuraminidases are important targets as they involve in the early stages of infection. Preventing the attachment of NA with cell surface receptors block the viral progeny and the further spread. The catalytic center of NA is located in the globular head region and the catalytic function of the NA enzyme is carried out through four major steps: The first step is the transformation of the carboxylate group from the axial position to the pseudo-equatorial position due to the charge-charge interactions to the arginine triad (118, 292 and 371) and steric constraints with tyrosine 406. The second step is the proton donation from solvent water molecules and the formation of endocyclic sialosyl cation transition-state intermediate. The third step of the catalytic pathway involved in the nucleophilic attack of tyrosine on the sialosyl cation. Final step involved the formation and release of sialic acid (Gong et al., 2007;Janakiraman et al., 1994). It was reported that in the early nineties, the catalytic role of sialic acid bound neuraminidase enzyme of influenza A and B viruses were experimentally studied, which provide the new opportunity to design more potential NA inhibitors for the flu virus (Burmeister et al., 1992;Varghese et al., 1992;von Itzstein et al., 1993). Numerous experimental studies reveal that the inhibiting potency of sialic acid based inhibitor against the NA enzyme by enzyme inhibition assay (Colombo et al., 2016;Ikematsu et al., 2012;Kim et al., 1998;Leang & Hurt, 2017).
M2 ion channel protein is another important enveloped protein, which maintains the pH value of the virion for uncoating of the virus during the early stages of the replication process (Balgi et al., 2013;Shu et al., 2011). Currently, H1N1 swine influenza and H5N1 avian Influenza are the emergent threats circulating in both animals and humans. Recently in India (January 2021), more than 50,000 poultry birds are killed due to the H5N1 influenza virus infection. Worldwide, M2 ion channel protein inhibitor (amantadine and rimantadine) and NA inhibitors (peramivir, oseltamivir and zanamivir) are used to treat Influenza viral infection (Bozdaganyan et al., 2014;Han et al., 2018). However, M2 ion channel protein inhibitor has limited applications due to their side effect. It was reported that all the NA inhibitors are effective only if they are given within the first two days of the manifestation of the symptoms (Glanz et al., 2018). There is no report on low pathogenic influenza virus resistance against NA inhibitor; on the other hand, mutation of influenza virus resists the developed antiviral drugs (Kannan & Kolandaivel, 2018). Hence the design of inhibitor molecules to HA and NA is quite crucial to avoid viral attachment, viral replication and viral infection. Molecular modeling methods such as Molecular docking and Molecular Dynamics (MD) provide the essential understanding of the processes and details about conformations and energies associated with the viral infection.
Molecular Docking is more commonly used in modeling the ligand molecule (small molecule) at the binding site of a macromolecular structure like protein. In molecular docking, docking score binding affinity values delineate the conformation of the ligand along with its position and orientation within the affinity of the protein-ligand complexes. MD simulation is a viable tool for rational drug design and this method is useful to study the protein folding, protein stability, conformational changes and conformation of carbohydrates, protein-carbohydrate interactions and ion transport in biological systems. It provides a wealth of information at the atomistic level of the dynamics of protein-carbohydrate interaction (Chang et al., 2011;Middleton et al., 2012;Parasuraman et al., 2015). In our earlier studies, we have carried out the sialic acid analog inhibitors for H5 hemagglutinin of H5N1 influenza virus (Jeyaram et al., 2019;. The present study is focused on investigating the role of N-acetylneuranamic acid (Neu5Ac, sialic acid, SIA) and zanamivir as effective inhibitors to the N1 neuraminidase of H5N1 influenza A virus by carrying out MD simulations which might help in designing the glycan based drug to inhibit influenza A virus.

Molecular docking and MD simulation methods
The protein data bank contains the 3-D structures of only group2 neuraminidase (N2, N9) complexed with N-acetylneuranamic acid. For group1 neuraminidase (N1) various sialic acid analog (zanamivir, oseltamivir, etc.) complex structures are available. But, no crystal structure is available for the N1 neuraminidase -N-acetylneuranamic acid complex system. Hence, a high resolution, crystal structure of N1 neuraminidase of H5N1 Influenza A virus in complex with ZMR is taken from the protein data (PDB id is 3CKZ). pdb4amber module of Amber 16 package is used to filter out the amino acid residues from the complex structure. Three-dimensional structure of SIA is taken from the 3dsdscar database (Veluraja et al., 2010). The pdbsum database and the online binding site prediction servers like Protein-Ligand Interaction Profiler, RaptorX are used to predict the binding site residues. Consurf web server is used for multiple sequence alignment and binding site prediction. Different parameters like the Ramachandran plot, QMEAN/ Z-score are used to check the quality of the initial structure. Autodock 4.2 (Morris et al., 2009) software is used for docking of SIA and ZMR into the active site residues of the N1 neuraminidase structure. With the aid of Autodock tools (ADT), the required polar hydrogen atoms and Kollman charges are added to the N1 protein structure and then non-polar hydrogen atoms are merged. The Gasteiger charges are added to the ligand atoms. The grid size is set to the range of 74 Â 74 Â 74 in xyz direction and the grid spacing set as the default range of 0.375 Å. The Lamarckian Genetic Algorithm (GA 4.2) is used as the docking parameter which gives the top ten binding score for the protein-ligand complex structure.
Based on analysis, the best binding pose of docked N1 neuraminidase of H5N1 Influenza A virus complexed with SIA and ZMR structures are taken as the initial structures to carry out the MD simulation by Amber 16 package (Case et al., 2016). The input parameter and topology files are generated using tleap module and the details are given in Table 1. Generalized amber force field (gaff) and ff14SB force fields are assigned for ligands (SIA, ZMR) and N1 protein respectively. Solvent water molecules from the TIP3P library and necessary counter ions are added to neutralize the complex system. The energy minimization of the system performed for 2000 minimization cycles in which the first 1000 cycles followed the steepest descent algorithm and then the conjugate gradient is used for the next 1000 cycles. The non-bonded cutoff interaction is set up with 8.0 Å. The system was heated and then the equilibration was carried out for 3,00,000 steps with the time step of 2 fs. SHAKE algorithm was used to constrain all bonds involving hydrogen and the system was equilibrated by NVT ensemble up to 300 K. The simulation temperature was fixed with 300 K and constant pressure (NPT ensemble) was used for 120 ns. Particle Mesh Ewald Molecular Dynamics (PMEMD) module (Salomon-Ferrer et al., 2013) is used to perform the MD simulation and the output trajectories are collected for every picosecond. To get insight into the interactions at the atomistic level, the collected trajectories are analyzed through Amber Tools17, Visual Molecular Dynamics (VMD), Nanoscale Molecular Dynamic Simulation (NAMD) and with the in-house developed FORTRAN programs. Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculation of CPPTRAJ module is used to calculate the binding free energy of the complex system by the formula DG bind ¼ DG complex -(DG N1 þDG lig ) . The pair interaction energy calculation between the active site residues of N1 NA and the ligands (SIA, ZMR) is computed by the NAMD software. Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) are calculated by the CPPTRAJ (Roe & Cheatham, 2013) module of Amber 16. Chimera (Pettersen et al., 2004), Molscript (Kraulis, 1991) and Raster (Hijmans & Etten, 2012) are used to rending the modeled structure.

Binding site prediction and molecular docking analyses
The N1-ZMR complex (3CKZ.pdb) is taken as the reference complex structure and its binding site residues are identified from the pdbsum database (Laskowski et al., 1997) Figure 1. RaptorX is a protein structure prediction server, the template based modeling used to predict the protein secondary and tertiary structures, contacts, solvent accessibility, disordered regions and binding sites. The RaptorX prediction is R118, E119, D151, R152, W178, I222, R224, S246, E276, E277, R292, N294, R371 and Y406. The Consurf server uses the phylogenic relation between the homologous sequences to estimate the conservation of amino acid/nucleic acid positions in a protein/DNA/RNA molecule. Multiple sequence alignment was done by using MAFFT and the homologous sequences were collected from the UNIREF90. The list of 150 homologues sequences is sampled to the query sequence. From this detailed sequence analysis, it is observed that most of the amino acid residues are structurally conserved as shown in Figure 1 and the predicted binding site residues are R118, E119, D151, R152, R156, W178, S179, I222, R224, E227, S246, E276, E277, R292, N294, G348, R371 and Y406. It is reported in literatures, that the active site of the NA enzyme contains 8 highly conserved amino acid residues (R118, D151, R152, R224, E276, R292, R371 and Y406) and 11 framework residues (E119, R156, W178, S179, D198, I222, E227, H274, E277, N294, and E425) in all influenza A and B viruses (Orozovic et al., 2014;Pizzorno et al., 2012;Tang et al., 2019;Yang et al., 2013). Molecular Docking of ligands (SIA, ZMR) into the active site of N1 neuraminidase of H5N1 influenza virus is performed using Autodock 4.2 software. The docking study provides the top ten structurally similar clusters with the increasing binding free energy, and the lowest binding free energy cluster is taken as the best binding pose of the ligand into the binding site of N1 protein. In the best binding pose of ZMR into N1, the binding site residues R118, E119, D151, R152, R156, W178, E277, R292, Y347, R371 and Y406 are interact with the ZMR which is consistent with the active site residues predicted by the pdbsum database and the online web servers (PLIP, RaptorX and Consurf). The substitution of guanidino group in O4 hydroxyl of Neu5Ac2en is known as ZMR and this guanidino group interacts with E119, D151 and W178 of viral neuraminidase, which is the same trend of interactions that are earlier reported in the previous experimental studies by Varghese et al. (1995) and Taylor et al. (1998). Moreover, the docked ZMR is bound to binding site residues in half-chair conformation, which is similar to the 4-guanidino-Neu5Ac2en (ZMR) inhibitor model evaluated by von Itzstein et al. (1993) and Varghese et al. (1995). They have reported that both 4-guanidino-Neu5Ac2en and SIA bind in half-chair conformation with the active site residues of NA of influenza virus. The metadynamics studies on the conformational free energy of SIA by Spiwok and Tvaroska (2009) reported that, in a planar system, the SIA based inhibitors zanamivir and oseltamivir bind in half-chair conformation and SIA binds in 4 S 2 /B 2,5 conformation with NA enzyme. The interacting residues of the binding site for the docked N1-SIA complex are similar to the docked N1-ZMR complex. Interestingly the docked SIA is binds in 2 C 5 chair conformation which is the minimum energy conformation of SIA (Priyadarzini et al., 2012a;Veluraja & Rao, 1980). The best binding pose of N1-SIA complex and N1-ZMR complex structure derived from molecular docking is schematically shown in Figure 2 and the inhibitory constant values and the estimated free energy of binding scores are given in Table 2. The experimental inhibitory constant value for N1-ZMR is 1.9 nM (Collins et al., 2008).

Quality check for the initial structure
Ramachandran plot, QMEAN are used to verify the quality of the experimental determination of the 3-dimensional coordinates of the protein structure. Ramachandran plot is a plot between the two torsional angle u and w of the peptide chain. PROCHECK webserver (Laskowski et al., 1993) is used to depict the torsional angles and given in Figure 3a. Ramachandran plot depicts that the conformation of amino acid residues fall in the most favoured regions (84.8%), allowed regions (14.9%) and generously allowed regions (0.3%). A good quality model is expected to have over 90% of the amino acid residues is in the most favoured regions. QMEAN is a Qualitative Model Energy Analysis used to describe the geometrical properties of the protein structure by generating a scoring function (Benkert et al., 2008). QMEANDisCo scoring method of Swiss-Model web server (Studer et al., 2020) is used for the QMEAN calculation and its global score is 0.95 ± 0.05. The local quality of the amino acid sequence is shown in Figure 3b.

ADME prediction
SWISS ADME web server (http://www.swissadme.ch/) (Daina et al., 2017) is used to evaluate the absorption, distribution, metabolism, excretion and toxicity of the ligand (SIA, ZMR) molecules. Both molecules have molecular weight of less than 500 g/mol and their lipophilicity, pharmacokinetics and synthetic accessibility values are close to each other. SIA obeys the Lipinski rules with the violation NH/OH >5 and ZMR passes three Lipinski rules with the violations of the sum of N&O > 10, NH/OH >5. The ADMET properties are given in Table 3, from this analysis it has been observed that each molecule can be used as a drug candidate for the inhibition of NA of influenza viruses.

Post simulation analysis
The docked structures of the N1-SIA complex and N1-ZMR complex are taken as the starting structure for the MD simulation. The MD simulation of duration 120 ns is carried out through PMEMD module of Amber 16 and the output trajectories are collected for every picosecond that leads to a total of 120000 frames. To get insight into the simulated trajectories, the analyses are done for each complex system. Measurements of RMSD and RMSF are used to identify the conformational stability of the complex structure. Hydrogen bonding interaction and hydrophobic contacts stabilize the complex structures. Binding affinity between N1 NA of H5N1 influenza virus and the ligand (SIA, ZMR) complex structures are computed by the use of NAMD pair interaction energy calculation (Phillips et al., 2005) and MM-PBSA binding free energy calculation (Miller et al., 2012).

Binding affinity measurement
The NAMD pair interaction energy is computed between SIA and the active site residues R118, E119, D151, R152, S246,    51-120 ns) trajectories of N1-SIA and BM1 of N1-ZMR. NAMD pair interaction energies and MM-PBSA binding free energies corresponding to BM1 and BM2 are listed in Table  2. The arginine triad (R118, R292 and R371) plays a vital role in binding specificity of the N1-SIA complex structure (Figure 4b). The relative binding energy difference between BM1 of N1-SIA and N1-ZMR is 121 Kcal/mol in NAMD and 24 Kcal/mol in MM-PBSA. The NAMD pair interaction energy calculations show a larger magnitude of energy difference in comparison with MM-PBSA. In our earlier study, we mentioned that the difference in the order of magnitude of energy is due to the non-consideration of water molecular interaction in NAMD pair interaction energy calculation (Jeyaram et al., 2019). However, both NAMD and MM-PBSA calculations confirm the higher binding affinity of BM1 of N1-SIA in comparison with N1-ZMR complex.

Measurement of RMSD and RMSF
RMSD and RMSF are the measures of quality of reproduction of the given structure and are used to check the conformational changes, protein stability and protein folding of the protein-ligand complex structures (Sargsyan et al., 2017). RMSD and RMSF of the N1 protein and ligand (SIA, ZMR) complex system are calculated using CPPTRAJ module of AMBER 16. RMSD measurement depends on the binding interaction and binding energy of the protein-ligand complex system and is calculated for the simulation time in picosecond. The RMSD of N1-SIA complex structure deviates from 0.75 Å to 1.9 Å and the N1-ZMR complex has the deviation from 0.75 Å to 3 Å. values which are measured for the Ca atom of each amino acid residue. From the RMSF plot, it is clear that, during the simulation period, the amino acid residues are remains in a same kind of orientation and it evident that the readjustment of the orientation of the ligand molecule to keeps the complex structure. The RMSD and RMSF plots are shown in Figure 5.

Hydrogen bonding analysis
Hydrogen bonding and hydrophobic interactions provide the structural stability of the protein-ligand complex system. The hydrogen bonding interactions and hydrophobic contacts of trajectories of N1-SIA complexed structures are analyzed using VMD software (Humphrey et al., 1996) and the representative projection of interaction pattern in BM1 and BM2 are shown in Figure 6. During simulation, slight changes in the orientation of binding site residues occur and their interactions with different hydroxyl position of SIA are observed. Interestingly N1 can hold the SIA in 2 C 5 chair confirmation at its binding site, the preferred conformation of SIA in its solution state. Raab et al. (2011)  The residential times of water molecules involved in water mediated hydrogen bonding interactions are calculated using in-house programs in randomly selected frames existing in different trajectories belonging to particular binding modes. Finally, each residential time of water interaction is averaged to calculate the average residential time of water interactions for a particular interaction. In the present work, the rigorous analysis shows that all the water mediating hydrogen bonds have an average residence time of more than 5 ps. Besides that, the clustering of hydrogen contributing to hydrophobic contacts is observed in the MD trajectories. The H3 and/or H4 of SIA forms hydrophobic contact with W403, I427, R430, P431 and K432. The methyl group of SIA interacts with S246, R292, N294, Y347, R371 and W403 contributing to hydrophobic interaction. The glycerol side chain hydrogens (H7, H8 and H9) are involved in forming a hydrophobic binding pocket with R152, I222, R224, S246, N247, N294, P326, A346, Y347, P431 and K432. It is interesting to note that in BM1, the binding pocket has enough space to accommodate the galactose (GAL) residue of SIAa(2!3)GAL and SIAa(2!6)GAL (as verified by in-house program) and in doing so, the actual binding site was able to form two or more hydrogen bonds. The projection of the interaction pattern is shown in Figure 7.
The second binding mode (BM2) is observed between 36-41 ns and 51-120 ns. Due to the flexibility of 150 loop, the constituent S370, S372, T401, W403 and K432 residues interact with SIA and form a maximum of 5 direct and 4 water mediated hydrogen bonds. The side chain hydroxyl group of S370 interacts with O1 carboxyl or O2 hydroxyl and O6 oxygen of SIA or glycerol side chain hydroxyl of O8/O9 of SIA. The arginine R371 backbone nitrogen and oxygen make the hydrogen bonds with O1 carboxyl & O2 hydroxyl group SIA. Side chain nitrogen of W403 or side chain hydroxyl of S372 binds with O2 & O6 of SIA. The backbone carbonyl group of T401 forms a direct hydrogen bond with O7 or O9 or O10 of SIA. Either the backbone carbonyl group or side chain nitrogen group of K432 forms direct or water mediated hydrogen bonds with different hydroxyl groups of SIA. The BM2 hydrogen bonding interactions are given in Figure 6b. The reduction in the number of hydrogen bonding interaction between SIA and active site residues in BM2, when   compared with BM1, is greatly reflected in the binding energy differences based on NAMD pair interaction energy and MM-PBSA binding free energy calculations. The hydrogen bonding interaction of N1-ZMR complex shows that the arginine triad (R118, R292, and R371) interacts with the carboxylate group of ZMR through direct or water mediated hydrogen bonding up to 5 ns, after that it loses its interaction with the carboxylate group and R118 interacts with guanidine group of ZMR. Side chains of S246, Y347 and Y406 partially interact with the carboxyl group. NH1, NH2 of ZMR hold the amino acid residues R156 and W178 through direct hydrogen bond up to 10 ns after that, the bond is mediated by a water molecule. In the entire simulation, the acetamide nitrogen N5 interacts with the side chains of E277 & E227, O7 and O10 of ZMR interact with the side chain of R152. The glycerol side chain (O8) interacts with N294, Y347/R224 & S246, and O9 of ZMR interacts with the R292 and N294 or N221, R224, G244 and S246. The details of hydrogen bonding interactions are given in Figure  8 and Table 5.

Protein allosteric prediction analysis
In a molecule, perturbation in one site by an effector causes a functional change in another site which is called as an allosteric site. These are the most common form of powerful mechanism for the regulation of protein function (Greener & Sternberg, 2015). In this study, the PARS web server (Panjkovich & Daura, 2014) is used to find out the allosteric ligand binding site of the protein structure. Eight ligand  binding sites are predicted by the PARS analysis and are given in Table 6. The first ranked site has the highest possibility of being an allosteric site and is displayed as a red ball in Figure 9. Cavities having the P-Value less than 0.05 Table 4. The hydrogen bonding interaction and water mediated hydrogen bonding interaction between N1 and SIA in the N1-SIA complex structure in binding mode   are considered as relevant to have the allostery and cavities whose structural conservation is more than 50% are considered as conserved. For the N1-ZMR complex, the amino acid residues D151, S246, E277, Y347 and Y420 and for in BM1 of N1-SIA complex, the amino acid residues E119, R152, E276, Y347, G348 and Y406 are found in the allosteric site as predicted by PARS and the allosteric site is given in red color in Figure 9. The cavities colored in cyan are structurally conserved and not having the allosteric site.

Discussion
The small molecule drug candidates SIA and ZMR are docked into the binding site of N1 neuraminidase of H5N1 Influenza A virus using Autodock 4.2 and each complex system is simulated over 120 ns by Amber16. Pdbsum database and few binding site prediction servers are used to find out the binding site residues, the active site residues of docked structures are compared and are similar to the result of prediction servers. The Ramachandran plot, QMEAN analyses confirming the good quality of the protein structure. The Swiss ADME analysis displays the physiochemical, pharmacological and drug like properties of the ligand molecules. The docking studies show that SIA is in 2 C 5 chair conformation and ZMR is in half-chair conformation in the binding site of N1 NA. Experimental studies show that ZMR is bound in half chair and SIA is bound to NA enzyme through boat or skewed boat or half chair conformation (Burmeister et al., 1992;Figure 8. Hydrogen bonding interactions between N1-ZMR complex structure.  Taylor et al., 1998;Varghese et al., 1995;von Itzstein et al., 1993), whereas our modeled SIA was bound in chair conformation. From the MD output trajectories, different parameters are used to find the binding affinity of each complex system and are given in Table 2. NAMD pair interaction energy calculation predicts that the BM1 of N1-SIA complex has a much better binding affinity than N1-ZMR complex. This is further validated by MM-PBSA binding free energy calculation. The RMSD and RMSF measurement ( Figure 5) shows that the binding stability of each complex structures and N1-SIA complex has lesser deviation and fluctuation than the N1-ZMR complex. Hydrogen bonding analysis shows that the interaction between the carboxylate group of SIA and the catalytic triad residues (Arginine) R118, R292, R371 and Y406 plays a dominant role in structural stabilization of N1-SIA complex structure; the same kind of structural stability was reported by Raab and Tvaroska (2011). It was previously reported that the negative charge on the carboxylate group of SIA and positively charged side chain of arginine triad (R118, R292 and R371) forms the strong charge-charge interaction (Gamblin & Skehel, 2010;Gong et al., 2007;Kim et al., 1999;Russell et al., 2006). Vavricka et al. (2013) has reported that the Y406 is a crucial factor for the catalytic activity and function as a nucleophile in the anomeric position of SIA. On the other hand, due to the very less negative charge on the carboxyl group of N1-ZMR complex, the arginine triad loosely binds to the carboxyl position, that too only for first 5 ns duration. Also, after 10 ns, R156 and W178 establish water mediated interaction with the guanidine group of ZMR. These may be the crucial factors for very less binding energy predicted by NAMD and MM-PBSA. The arginine triad (R118, R292 and R371) interaction with the carboxyl group of SIA is relatively stronger than the interactions of R152, N221, R224, G244 and S246 with the glycerol side chain of ZMR. The catalytic center of NA contains 19 highly conserved amino acid residues (R118, D151, R152, R224, E276, R292, R371, Y406, E119, R156, W178, S179, D198, I222, E227, H274, E277, N294, and E425) in all influenza A and B viruses (Orozovic et al., 2014;Tang et al., 2019;Yang et al., 2013). Amaro et al. (2007) and Han and Mu (2013) reported that N1 was able to adopt a wide variety of conformations due to the flexibility of 150 loop around the active site opens the formation of a large cavity adjacent to the active site. In our study, the amino acid residues N325, A346, G348, S370, S372, W403, and K432 interact along with the active site residues in BM1 of N1-SIA complex and N221 and G244 interact with the binding site of N1-ZMR complex. It is observed that the slight changes in the orientation of active site residues form at least 12 hydrogen bonds in each trajectory of BM1 of SIA. To understand the relative binding stability of the N1-SIA complex, the snapshots of the substrate and the binding site residues at an interval of 10 ns are given as the supplementary material. The reason for the changes in the orientation of the active site of NA is that the viral infection process is initiated by HA protein which is tightly fused with the SIA receptor and then it releases the NA. Hence the NA tends to readjust to bind and cleave the terminal SIA. The binding site of NA protein contains several charged amino acids to bind with the terminal SIA residue through hydrogen bonding or charge-charge interactions and also the active site residues are highly conserved in all NA subtypes (Kamali & Holodniy, 2013). In BM2, only six residues (S370, R371, S372, T401, W403 and K432) interact with SIA to form a maximum of 9 hydrogen bonds. All the analyses confirm that the BM1 of N1-SIA complex structure is the preferred binding mode. The PARS web server is used to find out the allosteric ligand binding site of the protein structure. Some residues (E119, R152, E276, Y347, G348 and Y406) in BM1 of N1-SIA and D151, S246, E277, Y347, Y406 in N1-ZMR found to have the allosteric regulation. Raab and Tvaroska (2011) studied the binding properties of the H5N1 influenza virus neuraminidase in complex with sialic acid and two trisaccharides (3SL & 6SL). They have started all the simulation with ligand molecules in boat conformation. During the simulation, ligand in N1-sialic acid complex changes its initial boat conformation into the 2 C 5 chair conformation; ligands in N1-trisaccharides complexes remain in boat conformation. In our studies, both the binding modes, N1 can accommodate the SIA in its minimum energy conformation ( 2 C 5 ) indicating that the NA need not spend any energy to accommodate SIA in its binding site. Besides that, the N1 binding pocket has the space to accommodate the galactose residues in a(2!3) and a(2!6) linkages along with the 2 C 5 chair conformation of sialic acid as checked by in-house Fortran program. As evident from the analysis of MD trajectories through RMSD, RMSF, NAMD pair interaction energy, MM-PBSA binding free energy, direct and water mediated hydrogen bonding interactions; SIA is preferred over ZMR by N1 NA of H5N1 influenza A virus.

Conclusion
In the present investigation, we have carried out 120 ns molecular dynamics simulations of N1-ligand (SIA, ZMR) complex structures. The binding affinities of SIA and ZMR with N1 NA are validated through binding free energy and pair interaction energy calculations. Furthermore, the analytical tools RMSD & RMSF measurements and hydrogen bonding analysis provide a much better understanding of the structural stability of the N1-SIA complex structure. All the analyses confirm that the BM1 of SIA to be the preferential binding mode. In both the binding modes of the N1-SIA complex, the N1 enzyme can accommodate the sialic acid in the 2 C 5 chair conformation and ZMR is bound in distorted chair conformation. Also, the NA binding site is able to accommodate the GAL residue of SIAa(2!3)GAL and SIAa(2!6)GAL. The inhibitory constant, binding affinity measurement by NAMD and MM-PBSA are the crucial factors, which will be compared to the kinetic properties of the further experimentation of N1-SIA complex structure. It can be concluded that SIA can be used as a drug candidate in inhibiting N1 NA of H5N1 Influenza A virus.