The Comparison of Substrate Stability in Neuraminidase Type 2 (N2) Active Site between A/Tokyo/3/67 and A/Pennsylvania/10218/84 with Heating Dynamics Simulation

A molecular dynamics simulation of two neuraminidase-sialic acid (NA-SA) complexes show a difference of the level of stability between sialic acid and neuraminidases that originated from viruses A/Tokyo/3/67 (Structure A) dan A/Pennsylvania/10218/84 (Structure B). Analyses of sialic acid RMSD and the change of torsional angles suggest that the sialic acid in Structure A is much more twisted and able to be influenced more by the binding of the neuraminidase functional residues than Structure B. Moreover, analyses upon hydrogen bond occupancy and binding free energy of both complexes showed that Structure A had more stable hydrogen bonds and each complex’s binding free energy were calculated to be –1.37 kcal/mol and 17.97 kcal/mol for Structure A and Structure B, respectively, further suggesting stability more dominant in Structure A than Structure B. Overall, Structure A has a more stable enzyme-substrate than Structure B.


Introduction
There are two pathogenicity levels of the avian influenza virus: the High Pathogenic Avian Influenza (HPAI) virus and the Low Pathogenic Avian Influenza (LPAI) virus.These levels of virus pathogenicity are classified based on the Intravenous Pathogenicity Index (IVPI) of a six-week-old chicken.An influenza virus is said to be highly pathogenic when it is able to kill >75% of chickens ages 4 -6 weeks in a ten-day window post-inoculation and features an IVPI larger than 1.2.The viruses that do not meet the HPAI criteria are called LPAI viruses [1].
Based on a few in vitro and in vivo studies [2][3][4][5][6][7][8], the enzymatic activity of neuraminidase (NA) participates in defining the pathogenicity level of avian influenza virus.Neuraminidase plays a major role in viral replication.During virus budding, neuraminidase cleaves the newly mature virion from its host cell.After the virion is detached from its host cell, it is then able to infect other host cells.In the previous studies, it was found that the neuraminidase of the HPAI viruses have the ability to cleave the sialic acid (SA) more effectively than the LPAI viruses.
Research efforts in neuraminidase inhibition throughout the years have resulted in three commercially available drugs (Oseltamivir, Laninamivir and Peramivir).However, at the rate in which the neuraminidase enzyme is mutating, these drugs could lose their potency over time.In other words, those three inhibitors are not stable in the confines of the neuraminidase binding pockets and may lead to the virus becoming resistant to the inhibitors.
Intricate studies in this area (NA-SA interaction) are very important, keeping in mind that the main principle of inhibition is to have a certain residue bind with the neuraminidase more than with the natural substrate, sialic acid.This study aims to shed light on neuraminidase so that structure-based drug design could solve the observed problem of neuraminidase resistance to drug molecules such as Oseltamivir, Zanamivir, Laninamivir and Peramivir in the future.
The computational approach of defining the correlation between pathogenicity and the neuraminidase's ability to bind to a substrate is still on the rise, and because of neuraminidase's rate of mutation, the knowledge in this field will complement experimental approaches to gain a more specific design to neutralize the neuraminidase activity more completely.While the standard intravenous pathogenicity index test can produce results in 10 days [9], a faster identification method through NASBA only focuses on hemagglutinin [10].There needs to be an alternative to these experimental tests that integrates all of the different methods and approaches in order to determine the viral pathogenicity in an accurate and timely manner.Computational approaches are valuable resources in this process, and this study serves as a starting point on how pathogenicity could be viewed on a molecular scale.
In this study, the structure of the A/Tokyo/3/67 avian influenza virus that was isolated in 1967 [11], during a time of prevalent infection, will be compared to the structure of the A/Pennsylvania/10218/84 avian influenza virus that is a non-pathogenic avian influenza virus [12].These avian influenza viruses were chosen to delegated two different level of pathogenicity.The relationship between the pathogenicity level of an avian influenza virus and the stability of the neuraminidase functional residues upon binding to sialic acid will be examined.Additionally, the natural substrate sialic acid will be examined in detail to observe its response to neuraminidase during molecular dynamics simulations.

Sequence Alignment and Template Searching
The amino acid sequences of both neuraminidases were obtained from the influenza database in NCBI [13] with accession code AAB05621 for the A/Tokyo/3/67 virus and BAF48360 for A/Pennsylvania/10218/84, which will be denoted as Structures A and B, respectively.A sequence alignment was executed using the fast pairwise alignment method that uses the BLOSUM 30 scoring matrix.
The x-ray crystallography structure of the neuraminidase-sialic acid (NA-SA) complex for Structure A [14] was obtained from the RCSB protein data bank [15] with accession code 2BAT.Furthermore, the crystal structure for Structure B was obtained by homology modeling as explained in the next section.

Refinement, Homology Modeling and
Explicit Solvation Process The generation of Structure B was done by homology modeling.Homology modeling was initiated by alignment of the amino acid sequence of the target structure sequence with the template.For the template itself, the 2BAT structure was chosen.In the 2BAT model, there is a sialic acid ligand in the structure.During modeling, the sialic acid was considered to be a rigid structure.
The template was refined by removing unwanted water and calcium molecules, and the remaining NA-SA structure was then used to generate Structure B. The modeling of Structure B was completed by characterizing the structure with the charmm forcefield, which in turn added the missing hydrogen atoms to the structure.Following structure generation, both Structure A and Structure B were solvated in a TIP3P water box.

Minimization and Molecular Modeling
Both complexes were subjected to two steps of energy minimization.The first step was executed by the Steepest Descent algorithm with a targeted energy gradient of 0.5 kcal/mol and a 1,000,000-step maximum.The second step was executed by the conjugate gradient with a targeted energy gradient of 0.1 kcal/mol and a 1,000,000step maximum.For both complexes, a nonbond list radius of 14 Å was used, and a switching function was applied between 10 -12 Å for computational efficiency.To gain a long-range electrostatic energy contribution, it was visualized in a spherical cutoff mode.
The molecular dynamics simulation was executed with heating during the first 20 ps of the simulation; the temperature rose from 0 -300 K.The parameters used were 20,000 steps, 0.001 time step, 0 K initial temperature, 300 K target temperature, nonbond list radius identical to that in energy minimization, and the trajectory data were stored every 0.1 ps.

End-Point Energy Calculation
Calculation of the ligand-receptor complex was based on the equation below: G is the average Gibbs energy, solvation G is the electrostatic and nonpolar free energy from implicit solvation.In this particular study, Generalized Born with Molecular Volume (GBMV) was applied for implicit solvation.The last term (TS) is the temperature and entropy contribution, while the first term of the right hand side ( MM E ) is the energy term produced by the applied forcefield, which is the potential energy of the system [16]: The relationship between the ligand, receptor and complex energy is given in the next equation: A/Pennsylvania/10218/84 with Heating Dynamics Simulation Equation ( 1) is average Gibbs energy which constructed each component energy in Equation (3).Furthermore, all phases of the study described in this section from structure preparation to molecular dynamics simulation were conducted with Discovery Studio 2.1 (Accelrys).

Results
A sequence alignment indicated that the two neuraminidases had 91% amino acid similarity.This similarity allowed us to generate a homology model for Structure B based on the 2BAT template structure, since it was higher than the required minimum of 50% similarity [17].In addition, energy minimization of the solvated structures resulted in a decrease in energy to -711, 107.21 and -599,227.56kcal/mol for Structure A and Structure B, respectively.
Heating simulation was executed for the first 20 ps (20,000 steps) of the molecular dynamics simulation (from 0 -300 K for each complex) to raise the system temperature to room temperature.The heating process is illustrated in Figure 1(a) for both the high pathogenic and low pathogenic complexes.It can be seen that the increase in temperature from 0 -300 K occurs between 0 and 2.5 ps, while the rest of the molecular dynamics simulation continued in equilibrium at a 300 K average temperature.The increment of system temperature in 2.5 ps designed to computational efficiency.For both systems, the increase in temperature and the subsequent thermal stability of the system appeared to be similar.This can be seen by the overlapping curvatures in Figure 1(a).This could suggest that the interactions of both complexes are along similar energetic pathways.
To evaluate the stability of the systems during simulation, the root mean square deviation (RMSD) of the backbone and all atoms of the complexes were calculated for all conformations throughout the simulation.As a whole, the structural stability of the systems was wellmaintained and is illustrated in Figure 1(b).In that diagram, the RMSD of all atoms was below 1 Å and had very little fluctuation.The RMSDs for both systems were 0.89 Å and 0.93 Å for Structure A and Structure B, respectively.Furthermore, the movement and change in stability of the neuraminidase molecules were not that significant.This is indicated by the backbone RMSDs that were well below 0.6 Å (0.52 Å for Structure A and 0.51 Å for structure B).The RMSDs of all neuraminidase atoms suggest that both systems behaved in a similar manner.
A comparison of the structures of both backbones of neuraminidase and its substrate (sialic acid) at the end of the simulation is shown in Figure 2. The superimposed   structures of the complexes, including the neuraminidase (shown as a ribbon) and the sialic acid (shown as a tubular shape), suggest that they do not differ from each other significantly at the end of the molecular dynamics simulation.This does not mean that the non-bonding interactions in the complexes are the same.Therefore, the observation of the interactions between the substrate and the functional residues are necessary to determine the cause of the pathogenicity of avian influenza virus neuraminidase.
Calculation of the binding free energy at the end of the simulation with added implicit solvation using the Generalized Born with Molecular Volume (GBMV) method, resulting in the values of -1.37 kcal/mol for Structure A and 17.97 kcal/mol for Structure B. The resulting values indicate that ligand binding is more favorable in Structure A than in Structure B [18,19].

Overall Substrate Stability
Figure 3 depicts the overall movement of the sialic acid in response to the hydrogen bonds that form between the sialic acid and neuraminidase and the long range interactions in the 14 Å spherical cutoff range.As shown in Figure 3, the substrate RMSD drastically increased from the starting point to the 6000 th step until it reached the 0.65 Å mark.This could be caused by the increase in kinetic energy of the atoms with the rise in temperature.After the heating phase, both substrates were observed to be relatively stable until the 9000 th step.Between the 9000 th and 16,000 th step, the RMSD of Structure A's sialic acid decreased to 0.45 Å.This suggests that there is higher electrostatic interaction towards the initial position of sialic acid in Structure A than Structure B. The following steps of the curvature showed that the RMSD fluctuation of the sialic acid was not that different in each structure.The final RMSD of the sialic acid in Structure A was measured to be 0.52 Å, while 0.64 Å was measured for Structure B.

Substrate movements
The motivation to understand the movement of the substrate is not only to see its response in active site of two neuraminidases with different pathogenicity, but also to give information about nature way of interaction of NA-SA.It is needed to improve and produce better inhibitor drugs in the future.In order to understand the substrate movement specifically, the substrate was divided into four torsional sections, with each section representing the movement of the binding pocket area that interacts with the functional residues of the neuraminidase active site.A few of the sialic acid atoms that form hydrogen bonds with the functional residues were characterized by Stoll et al. [20] into a number of sections.The order in which the torsional angles were arranged to the sialic acid carbon atoms was also defined similarly by Aruksakunwong et al. [21] and Mao Su et al. [22].
The first torsional angle was formed by atoms O6-C2-C1-O1A (see Figure 4).These atoms correspond to the first carboxylic binding pocket (C1) that is negatively charged and interacts with the positively charged 118-292-371 arginine triad subsite (S1).The bound sialic acid in Structure A changes its torsional angle 1 drastically.This is in contrast to the more subtle changes observed in Structure B. The interaction of S1 with C1 is stronger for Structure A than for Structure B (Figure 5).
The 1 of Structure A may be showed like unstable oscillation (Figure 6a).It is because we look it as two dimensional graphic.If we place the torsional change graphic in the cylinder, and bring to meet -180˚ with 180˚, it could clearly show that the torsional changes experiences only a 1-20˚ torsional angle change and not much different than Structure B torsional changes graphic.The differences in the torsional angle 1 from the beginning of the simulation indicated the rotation of the carboxylic head O1A and O1B upon energy minimization.mutated residues are D147G, V149I, I194V, K199R, V275I, T346N, Q347P, L370S, R403W, and K431P.In observing these residues in depth, the mutated residues do not correspond with the decrease in the torsional angle directly, but suggesting a repulsive force that may shift the neighboring backbone in the neuraminidase.
The process of the decrease of 2 was measured from the initial state to the 5000 th step of the simulation.This is thought to give a contribution to the sudden rise in the RMSD of all atoms of the substrate from step 0 -6000 th .
The third torsional angle (3) is defined as the angle formed by the atoms C5-C6-C7-C8 (Figure 4), while the fourth torsional angle (4) is formed by atoms O9-C9-C8-C7 (Figure 4).This region is chosen to observe the response of the C6 carbon pocket to the atoms that interact with subsite S5 (formed by A246 and E276).
The second torsional angle 2 was defined by atoms C4-C5-N5-C10, as shown in Figure 4.These atoms represent the area in which two carbon pockets interact with the functional residues (i.e., the C4 pocket that interacts with D151 and subsite S2, and the C5 carbon pocket that interacts with subsite S3).From Figure 6b, the 2 decreases from -114˚ to -145˚ for Structure A and from -130˚ to -153˚ for Structure B. This decrease is a response to the interaction of a negatively charged area (formed by D151 and subsite S2 composed of E119 and E227) with the C4 carbon pocket, which had been identified by Taylor et al. [23] as positively charged and indicative of a difference in strength in the S2 subsite.The difference in distance between D151 and subsite S2 from C4 in both complexes affects the attracting forces that the atoms in 2 are experiencing.This in turn could explain the difference in torsional angles from the initial to the final step of the simulation.By selecting the atoms at the 12 Å non-bonded spherical cutoff marks from the atoms in 2 there are a few residues (that were already mutated during homology modeling) that might experience non-bonded interactions, directly or indirectly, that could be observed.Those The differences in the torsional angles in the various figures suggest that the sialic acid in Structure A is much more twisted than that in Structure B. This is based on the large torsional angle changes but regular RMSD fluctuations experienced by the sialic acid in both complexes.

Sialic Acid Interaction with Neuraminidase
Functional Residues The functional residues of the avian influenza neuraminidase those mediate binding with the sialic acid  through hydrogen binding were examined in more detail.This was done to compare the hydrogen bond contribution of both complexes towards the structural stability throughout the simulation.
There are a few differences in the hydrogen bond contributions between Structure A and Structure B. The carboxylic group that acts as a main attraction particularly for the S1subsite that initiates binding [24], is observed to form four simultaneous hydrogen bonds in the final step of the simulation of Structure A (Figure 7(a)).In contrast, Structure B only had one hydrogen bond for that region (Figure 7(b)).The formation of the hydrogen bond at R292 emphasizes the importance of this residue in sialic acid binding, which is also supported by studies that describe inhibitor resistance caused by an R292K mutation [25][26][27].The strong interaction with subsite S1 raises a few new questions concerning the mutated framework residues in the 12 Å radius of the sialic acid and S1 residues status as functional residues.One such area of inquiry is the observed hydrogen bond occupancy, where Y406 acts as a functional residue in the S1 of Structure B. We are currently studying this subject.
Another difference that is also important is the hydrogen bond formed by R152 in Structure A but not in Structure B. The interaction with R152 seemed to constrain D151 to move towards the sialic acid in Structure A, therefore sustaining the bond with E276.Additionally, the unstable interaction with R224 in Structure A raises a question concerning its role.By direct observation, D151 of Structure B is a focal point of the increase in the torsional angle 2, created when the number of hydrogen bonds decreases from two to one.
The difference in hydrogen bond influence could be observed by the occupancy, where those with at least 80% occupancy are considered to be strong hydrogen bonds [22].From Figure 5 it can be seen that there are four residues with at least 80% occupancy in Structure A, while there were three in Structure B. However, there were no strong hydrogen bonds in subsite S1 of Structure B, indicating that this could influence the rate in which a mature virion is released from the host cell.Moreover, the more variable 1 in Structure A could relate to the electrostatic attraction of the sialic acid in S1 that negates interaction with other subsites.
From the hydrogen bond occupancy, the different overall binding strength was also visible around E276.In Structure B, the hydrogen formed by E276 was disrupted and broke a few times during the simulation.In Structure B, it could also be seen that the hydrogen bond formed by R224 is much stronger than that in Structure A. This is may be related with the stability of the bond at E276.

Discussion
The main purpose of this study was to compare the NA-SA interaction of two distinct viruses that possess different pathogenicity levels.In depth analyses of structure, energy, and electrostatic interaction of the substrate with the NA, functional residues were conducted to observe the correlation between pathogenicity and structural change throughout the simulation.And also, the similarity of both structures may enhance the curiousness about how could be the similar structures have different level of pathogenicity.
Comparison of the substrate movement indicated that there is a difference in sialic acid response to the neuraminidase, albeit not significant.For the highly pathogenic Structure A, the substrate that is bound by the neuraminidase possesses a lower average RMSD than  Structure B. The RMSD comparison here is compatible with multiple studies that infer that a higher substrate RMSD suggests the superior neuraminidase ability to reject an inhibitor [21,22,28].The torsional angle of the sialic acid in Structure A was far greater than its low pathogenic counterpart in Structure B, which suggests that the sialic acid was more twisted in the Structure A binding pocket than in the binding pocket of Structure B.
By comparing the functional residues, positions of those that are considered to be the functional residues in the neuraminidases are conserved and still satisfy the characterization of Stoll et al.There are many examples A/Pennsylvania/10218/84 with Heating Dynamics Simulation in which mutation of functional residues increases the potency of resistance against inhibitors [25][26][27][29][30][31][32][33][34][35][36].Evidence of this are the data shown in Figures 5 and 7 that depicts the residues principal to substrate binding and how mutation of certain residues alter the stability of the system.
A mutation that often occurs and is the main focus in type-2 neuraminidase (N2) studies is R292K.The mutation at this residue causes resistance to inhibitors.The level to which sensitivity to inhibition is lost also differs for each drug molecule [37], which indicates distinct viral pathogenicity.In this study, the R292 residue has relatively high hydrogen bond occupancy for the more pathogenic structure.This residue could very well play a very significant role in determining the pathogenicity level of a virus apart from other residues such as D151, E276 and R371 in N2.
Comparison of the binding free energy at the last step of the simulation based upon the procedure of Amaro et al. was used for computational efficiency [38].The ac- quired results were not far from the results of other studies that calculated the binding free energy of NA-inhibitor complexes with different approaches [39,40].
In another study, Masukawa et al. calculated the entire trajectory using MM/PBSA and successfully showed the correlation between the binding free energy with inhibition in an experimental study involving a type-9 neuraminidase (N9).The experimental results were in the form of an inhibition constant that was converted to G using the equation G = -RTlnIC 50 , and they were able to show agreement between the computational and experimental results.The binding free energy of the N9-SA complex in their study was calculated to be -1.69kcal/mol, and the conversion to an inhibition constant resulted in a value of -3.06 kcal/mol.The other drug molecules in the study had a binding free energy range between -15.2 kcal/mol and -6.91 kcal/mol that is similar to the experimental results.
Additionally, a similar study conducted by Smith et al. showed a table containing the inhibition constants from the study of von Itzstein et al. compared with the converted binding free energy of an interacting N2-inhibitor complex [41].It was calculated that the range of binding free energy was from -10.1 kcal/mol to -4.6 kcal/mol.In the same study, the values of binding free energy of an inhibitor were tabulated with inhibitors that have been modified.The range of the binding free energy was from -21.7 kcal/mol to 7.3 kcal/mol.The calculations from these different studies similarly suggested that the enzyme that bonded more tightly with the substrate possesses a more negative binding free energy value.Additionally, the order of the binding free energy values seems to agree with our study.

Conclusions
By heating dynamics simulation produced, we are able to determined different level of influenza virus pathogenicity.The more pathogenic influenza virus has more hydrogen bonds, higher percentage occupancy of hydrogen bond and more negative end-point binding free energy than the non-pathogenic.And also the different electrostatic interaction responded by different torsional angle changes of the same sialic acid.The experiment showed the different strength and stable interaction between NA and SA from the different level of influenza virus pathogenicity.

Figure 1 .
Figure 1.(a) Herlambang et al; temperature vs. simulation step graphic in both systems.

Figure 1 .
Figure 1.(b) Herlambang et al; RMSD all atoms and backbone NA vs. simulation step in both systems.

Figure 2 .
Figure 2. Herlambang et al; The superposition of neuraminidase A/Tokyo/3/67 complex (green) and A/Pennsylvania/10218/84 complex (gray).The flat ribbons show the backbone of both neuraminidases and the sialic acid bound in the active site "holy grail" shown as a stick molecule.
Focusing on 3, as shown in Figure6c, in the sialic acid of Structure A, there is widening of torsional angle from 149˚ until 168˚.This change is bigger than torsional angle change in Structure B which increases from 143˚ to 149˚.The different change in torsional angle is influenced by different interaction with E276 and also with the electrostatic attraction from the positively-charged side chain of R292.This could lead to the decrease in the angles formed by C8 and C7 and consequently the increasing 3.For 4 (Figure6(d)), the sialic acid of both complexes fluctuated from 50˚ -80˚.The sialic acid of Structure A was measured to be 64˚ and 63.5˚ by the end of the simulation.The sialic acid of Structure B was initially measured to be 65˚ then 77˚ by the end of the simulation.The difference of the torsional angles could be explained by the electrostatic interaction difference of E276 and E277.The electrostatic interaction of this region in Structure A is much stronger than in Structure B; this may be verified by the hydrogen bond occupancy (Figure 5).The O8 atom of the sialic acid that acts as a proton donor in a hydrogen bond with E276 is covalently linked with C8; thus, the change of 4 in Structure A was not of major significance.

Figure 5 .
Figure 5. Herlambang et al; the percentage of hydrogen bond occupancy of the functional residues during the simulation.

Figure 7 .
Figure 7.The functional residues which have a hydrogen bond (dashed line) in the last conformation of the NA-SA interaction in: a. Structure A (top); b.Structure B (bottom).