The Mechanism by which 146-N-Glycan Affects the Active Site of Neuraminidase

One of the most conserved glycosylation sites of neuraminidase (NA) is 146-N-glycan. This site is adjacent to the 150-cavity of NA, which is found within the active site and thought to be a target for rational drug development against the antiviral resistance of influenza. Here, through a total of 2.4 μs molecular dynamics (MD) simulations, we demonstrated that 146-N-glycan can stabilize the conformation of the 150-loop that controls the volume of the 150-cavity. Moreover, with 146-N-glycan, our simulation result was more consistent with crystal structures of NAs than simulations conducted without glycans. Cluster analysis of the MD trajectories showed that 146-N-glycan adopted three distinct conformations: monomer-bridged, dimer-bridged and standing. Of these conformations, the dimer-bridged 146-N-glycan was the most stable one and contributed to stabilization of the 150-loop conformation. Furthermore, our simulation revealed that various standing conformations of 146-N-glycan could block the entrance of the binding pocket. This result was consistent with experimental data and explained the relatively low activity of inhibitors with flexible substituents toward the 150-cavity. Together, our results lead us to hypothesize that rigid and hydrophobic substituents could serve as better inhibitors targeting the 150-cavity.


Introduction
Influenza viruses have extraordinary virulence, causing severe respiratory infection and even death. Every decade or so, a dangerous new strain appears and poses a threat to the public health [1]. Of the three types of influenza virus, type A infects a wide range of avian and mammalian species. Additionally, type A can be further classified into subtypes according to the serological reactivity of its surface glycoprotein antigens, hemagglutinin (HA) and neuraminidase (NA) [2]. Of all the serotypes, there are nine for NA (N1 to N9) and sixteen for HA (H1 to H16) found in avian and mammalian hosts [3]. The nine NA alleles have been divided into two groups based on phylogenetic analyses of the species (group-1: N1, N4, N5, N8; group-2: N2, N3, N6, N7, N9) [4]. Of these alleles, N1 and N2 are responsible for viral pandemics and recurrent annual epidemics in humans [5]. The influenza virus NA is a tetramer protein that catalyzes the cleavage of terminal α-ketosidically linked sialic acids from a large variety of glycoproteins, glycolipids, and oligosaccharides [6,7]. NA plays an important role in the final stages of influenza virus infection, removing sialic acid from both the infected cell surfaces and the newly formed virus. Simultaneously, NA facilitates the release of progeny viruses and thus furthers the spread of infection [8]. Given its importance in altering the progression of influenza infection, NA is an attractive target for structure-based antiviral drug design.
To date, only four anti-flu drugs have been approved by the Food and Drug Administration (FDA). Two (amantadine and rimantadine) target the M2 ion channel of the influenza virus [9]. The remaining two, oseltamivir (Tamiflu) and zanamivir (Relenza), both inhibit NA [10,11]. However, because these drugs are widely used, various NA subtypes have demonstrated resistance to the M2 channel inhibitors amantadine and rimantadine [12]. Fortunately, oseltamivir and zanamivir are still effective against new viruses. However, NA also faces selection pressure as mutations occur. Once drug-resistant strains arise, they will likely lead to the breakout of a novel flu, causing greater global public health concerns. Thus, to develop new drugs for NA, it is vital to proactively identify potential drug-resistance sites and identify the optimal binding modes of current drugs with NA.
Adjacent to the sialic acid binding site, a difference occurs between group-1 and group-2 NAs. In the crystal structure of NAs, the 150-loop has two distinct conformations, one open and one closed [13]. The majority of the crystal structures of group-1 NAs contain a 150-cavity that is rarely found in the crystal structures of group-2 NAs [13]. It has been hypothesized that targeting the 150-cavity may allow for the development of new antiviral agents with increased specificity and potency against group-1 enzymes [14]. Earlier computational studies revealed that a key salt bridge (D147-H150) predominantly controlled the closed conformation of the 150-cavity in group-2. The lack of this salt bridge provides flexibility to the 150-loop, which could be responsible for the open conformation in group-1 [13]. In these simulations, the open form of the 150-loop was noted as the main conformation in group-1 and was proposed as a new potential target for drug design [14]. The closed conformation previously observed in group-2 structures suggested that a slow conformational change may occur upon inhibitor binding to the NA of group-1 [15,16]. However, in the crystal structure 3NSS (group-1 NA), the closed conformation can also be observed without an inhibitor binding in the pocket [17]. Therefore, we hypothesize that there are additional factors that contribute to stabilizing the 150-loop of the NA in group-1. After investigating the 150-loop of NA, we found that none of these simulations considered the possible effects of N-glycan, which is located on residue N146 adjacent to the 150-loop.
In the past decades, numerous studies have demonstrated that glycosylation on HA and NA of influenza strains can affect host specificity, infectivity and virulence by altering the biological properties of HA and NA [18]. During the 1918 pandemic influenza virus, different N-glycan profiles may have caused viral resistance to proteinase digestion and lead to increased infectivity [18]. The glycosite located near the enzymatic active site on ASN146 of the NA N-terminus close to the 150-loop (residues 147-152) is highly conserved across approximately all strains [19]. To reveal the effects of the N-glycan on the conformational stability of the 150-loop and volume of the 150-cavity, we examined the flexibility of the 150-loop with and without 146-Nglycan using molecular dynamics (MD) simulations of the 09N1 crystal structure as well as other available structures of N1 and N2 alleles derived from human clinical samples. Additionally, we studied how 146-N-glycan affected the 150-cavity when the 150-loop was in an open conformation. Our results indicate that the 146-N-glycan could stabilize the 150-loop conformations and preferentially prevent transformation of the 150-loop from closed to open in both N1 and N2 enzymes. Conversely, the volume of the 150-cavity and the entrance to the 150-cavity could also be controlled and affected by the 146-N-glycan, which helps to explain the low activity of inhibitors targeting the 150-cavity of NA enzymes.

Molecular Dynamics Simulations
To probe the effect of the N-glycan on the conformation of the 150-loop, we performed six separate 100 ns MD simulations with AMBER 12 [20] for three tetrameric NA enzymes with and without a 146-N-glycan: (1) A/Tokyo/3/67, a seasonal human H2N2 virus (N2, PDB ID: 1NN2) [21]; (2) A/Vietnam/1203/04, an avian-derived H5N1 virus isolated from a human (VN04N1, PDB ID: 2HTY) [14]; (3) A/California/04/2009, an H1N1 virus isolated early during the 2009 pandemic (09N1, PDB ID: 3NSS) [17]. The 1NN2 structure exclusively contained a salt bridge stabilizing the closed conformation of the N2 sample. The 2HTY structure was selected to represent an open conformation of N1 samples. The 3NSS structure was selected to represent a closed conformation of N1 samples.
Simulating the homo-tetramer configuration of NA led to the equivalent of approximately 0.5 μs (400 ns) of sampling for each NA monomer while simultaneously incorporating neighboring subunit (monomer) structural effects. The tetramer systems and individual monomer chains exhibited stability over 100 ns in backbone root mean square deviation (RMSD) plots (S1 Fig). Additionally, the experimental and simulation-derived B-factors were consistent (S2-S4 Figs).

146-N-glycan stabilizes predominant l50-loop conformations and prevents conformational transformation in both N1 and N2
Our simulations revealed that 146-N-glycan stabilizes predominant 150-loop conformations and prevents the conformational transformation from closed to open in both N1 and N2 (Fig 1). Surprisingly, the pandemic 09N1 with 146-N-glycan maintained its preferential conformation with a closed 150-cavity in normal solution dynamics. This finding is in contrast with both 09N1 without 146-N-glycan simulation (Fig 1) and simulations of other groups [13,14]. All six systems (three with the 146-N-glycan and three without) were simulated using 100 ns MD simulations. The snapshots of the four chains obtained every 40 ps were clustered with cpptraj analysis tool of AMBER 12   146-N-Glycan Affects the Active Site of Neuraminidase of the 150-loop conformations remained closed. For 1NN2 without 146-N-glycan, however, the percentage of closed conformations was 55.1%. Based on the six trajectories, for both N1 and N2 with 146-N-glycan, greater than 20% more of the initial conformational form (whether open or closed) was noted for both N1 and N2 with 146-N-glycan compared with N1 and N2 without 146-N-glycan. This result indicates that glycan can stabilize the predominant 150-loop conformations and prevent the conformational transformation from closed to open in both N1 and N2. As previous simulations have revealed, the 150-loop may not exhibit high levels of flexibility in the presence of glycan [13][14][15].

Monomer-bridged, dimer-bridged and standing conformations of the 146-N-glycan
To study the mechanism by which the 146-N-glycan stabilizes the conformations of the 150-loop, conformations of the 146-N-glycan of three systems with glycan were clustered. A relative, stable reference plane (defined by CA atom of S145 and two β-sheet end CA atoms which have the following B-factors: S145:CA = 11.77, L223:CA = 8.07 and P301:CA = 8.10) was chosen to define the dihedral angle of the 146-glycan (S6 Fig). The glycan conformations on the protein were divided into three clusters: monomer-bridged, standing and dimer-bridged (Fig 2). Monomer-bridged is the conformation in which the glycan attaches to the same chain to which 146-N-glycan is bound. Dimer-bridged is the conformation in which the glycan bridges to the adjacent chain. The standing conformation is that in which the glycan stands perpendicular to the dimer interface, as presented in Fig 2. Both monomer-bridged and dimer- bridged conformations form structures that are more stable than those formed by the standing conformation by interacting with a monomer or adjacent chain or glycans. In three MD systems, the percentages of conformations among the clusters were different, as presented in Table 1. In 1NN2, 94.30% of the glycan conformations were standing or monomer-bridged. In 2HTY, 94.25% of the glycan conformations were dimer-bridged or standing. In 3NSS, 73. To further study the stability of the system, interaction energies between 146-N-glycan and the other parts of 3NSS NA were also calculated using AMBER 12 software (Fig 3; calculation details were provided in S7 Fig and S1 Text). As shown in Fig 3, the binding energies of monomer-bridged conformations were approximately 20 to 40 kcal/mol reduced compared with standing conformations. Additionally, the dimer-bridged conformation energies were approximately 40 to 80 kcal/mol lower than those of the standing conformations, which indicates that monomer-bridged and dimer-bridged conformations are more stable. Specifically, the dimer-bridged conformation, which displayed the lowest binding energy and    Fig  4B. The filled conformation (9%) represents the conformation in which 146-N-glycan occupies the 150-cavity, whereas the covered conformation (7%) represents the conformation in which the 146-N-glycan covers the entrance of the 150-cavity. Together, these two conformations are defined as the blocked conformation (Fig 4C), whereas the other conformation is defined as the open conformation. Based on these results, we hypothesized that the activities of inhibitors designed to target the 150-cavtiy might be affected by blocking the 146-N-glycan. To test the hypothesis, five inhibitors ( Fig 5) were selected for the following study (inhibitors 1 and 2 only target orthosteric sialic acid binding sites, whereas compounds 3, 4, 5 also target the 150-cavity). These inhibitors form co-crystal structures when used to soak crystals of N8 [A/duck/ Ukraine/1/63(H3N8)] (group-1 NA) [26], and they display experimental Ki values in vitro with H5N1 [23]. The binding energies of crystal conformations (S8 Fig) of inhibitors 1, 2, 3, 4,  Table 2, the calculated binding energies of inhibitors 3, 4 and 5 were similar to those of inhibitors 1 and 2, which indicates that all five inhibitors should exhibit similar inhibition activity. However, the experimental Ki values of inhibitors 1 and 2 (only binding to the sialic acid cavity) were considerably reduced compared with inhibitors 3, 4, and 5. One primary reason for the difference in activity between these two sets of inhibitors is that inhibitors 3, 4, and 5 might be blocked from binding to the entrance to the 150-cavity by the 146-N-glycan in 16% of conformations (as in Fig 4B). In contrast, for inhibitors 1 and 2, only a small portion of the entrance to the sialic acid cavity might be blocked by the 146-N-glycan. Furthermore, for inhibitors 3, 4, and 5, the most flexible inhibitor (inhibitor 3) displayed the lowest activity. This finding may be caused by the increased probability that a flexible substitute will interact with 146-N-glycan. Therefore, to prevent an inhibitor from interference with 146-N-glycan, the substitute should be longer and more rigid to reach the bottom of the 150-cavity where the hydrophobic residues V116, I117, A138 and V149 (Fig 4D) are located. Recently, in Yuanchao Xie's work [25], two series of inhibitors were synthesized, and the ability of these molecules to inhibit H5N1 was verified by experiments. Seven inhibitors were selected from two series (Fig 6) from Yuanchao Xie's work [25] and docked into the N1 protein (PDB ID: 2HU0). After a 2 ns MD simulation, the binding energies of the seven inhibitors were calculated using MM/GBSA in AMBER 12 ( Table 2). The binding conformations are provided in S9 and S10 Figs. The results obtained are consistent with our conclusions. Table 3 shows that the calculated binding energies of inhibitors 17h, 17l, and 20m (with flexible substituents) do not correlate with their low experimental activities. Alternatively, inhibitors 17e and 20l (with rigid substituents), which may not interact with the 146-glycan site, exhibit better experimental activities. Inhibitors 20l, which can interact with the hydrophobic zone of the 150-cavity (shown in Fig 7), exhibits a particularly high activity (IC 50 = 1.9 nm). Based on this model, we suggest that to avoid being blocked by the 146-N-glycan, inhibitors targeting the 150-cavity should have long, rigid substituents that can bind to the hydrophobic zone, which is surrounded by residues V116, I117, A138 and V149 of the 150-cavity. Furthermore, in 1993, Shengqiang Li's work [27] revealed that NA activity of Influenza A/WSN/33 can be reduced by 10-fold by a glycosylate at position 130 (corresponding to the 146-N-glycan of 3NSS). This result indicates that if the 146-N-glycan is shifted to block the entrance of the NA active site, it can reduce the activity of NA. Therefore, it may be useful to design compounds to prevent the 146-N-glycan from forming a dimer-bridge to provide more standing conformations. Such a compound could represent a new type of NA inhibitor.

Conclusion
Our MD simulations highlight the important effect of the 146-N-glycan on the 150-cavity of N1 and N2 systems across a dynamic ensemble of NA proteins. To our knowledge, this is the first study to employ long time-scale MD simulations of NA tetramer proteins complexed with glycan. Our simulation revealed that the 146-N-glycan could stabilize predominant l50-loop conformations and prevent conformational transformation in both N1 and N2. In particular, when simulating the 09N1 protein 3NSS, simulations with glycan showed the molecule to have the opposite preferential conformation to that without glycan. Cluster analysis of the MD trajectories captured three distinct conformational clusters of glycan on proteins: monomerbridged, dimer-bridged, and standing conformations, which correlate to the closed, open, and open+closed forms of 150-loop, respectively. These conformations therefore control the volume of the 150-cavity. Our simulations also revealed that in a small portion of conformations, glycan can block the 150-cavity. This finding could be a primary reason that existing inhibitors Seven inhibitors design to target 150-cavity. These inhibitors were come from Yuanchao Xie's work [25]. oseltamivir carboxylate (OS-C) and oseltamivir derivatives (9, 17e, 17h, 17l, 20m, 20l). Inhibitors 9 and OS-C have no substituents bound to 150-cavity. Inhibitors 17h, 17l, 20m have flexible substituents while inhibitors 17e and 20l have rigid substituents.
b MM/GBSA in AMBER 12 was used to calculated the binding energy.
c Binding conformation of OS-C is the co-crystal conformation (PDB ID: 2HU0).
targeting the 150-cavity exhibit lower bio-activities relative to those that bind to the original sialic binding cavity. Based on this model, we suggest that to avoid being blocked by the 146-N-glycan, inhibitors targeting the 150-cavity should have long, rigid substituents that can bind to the hydrophobic zone, which is surrounded by residues V116, I117, A138 and V149 of the 150-cavity. Furthermore, we propose that compounds designed to prevent 146-N-glycan from forming a dimer bridge might result in a greater proportion of the 146-N-glycans in the standing conformation, which could reduce the probability of binding to sialic acid. Molecules designed in this manner could represent a new form of NA inhibitor.

Preparation of the initial structures
Six simulation systems were created and designated as 3NSS-No-Glycan, 3NSS-With-Glycan, 2HTY-No-Glycan, 2HTY-With-Glycan, 1NN2-No-Glycan and 1NN2-With-Glycan. The atomic coordinates of 3NSS-No-Glycan and 2HTY-No-Glycan were extracted from 3NSS and 2HTY, respectively, by removing the glycans. The atomic coordinates of 1NN2-With-Glycan were extracted from 1NN2, whereas the atomic coordinates of 1NN2-No-Glycan were extracted from 1NN2 by removing the glycans. The glycans of 3NSS-With-Glycan and 2HTY-With-Glycan were created using the corresponding atomic coordinates of the fractional 146-N-glycan in 3NSS as a template, using a complex glycan structure type command that was obtained from the Glyprot Server (http://www.glycosciences.de/modeling/glyprot/) [28].
The structure of the NA tetramer was created from monomer structures by applying a VMD 1.9 matrix transformation using Tcl script [29]. Histidine residues were protonated at pH 6.5 using the PDB2PQR server 6 [30]. Disulfide linkages were labeled with the proper AMBER notation and the resultant files were processed using XLeap from the AMBER 12 [20]  146-N-Glycan Affects the Active Site of Neuraminidase program, using the AMBER99SB and the GLYCAM06 force fields. A truncated octahedral box of pre-equilibrated TIP3P water was added to solvate each system with a pad of 10 Å. Counter ions were added to neutralize each system.

Molecular dynamics simulations
MD simulations of all systems were conducted using a version of the PMEMD module from AMBER 12. For each system containing glycans at position N146, the minimization and equilibration procedures were as follows. To maintain the glycan direction as it appears in the crystal structure, the structures were minimized in four stages. First, glycans and protein were all restrained. Second, the first three glycosyl units of the glycans and the protein were restrained. Third, only the protein was restrained. Finally, all systems were unrestrained. The restraints were reduced in the four stages over 20000 steps of steepest descent and conjugate gradient minimization.
After minimization, the system was gradually heated from 0 K to 300 K in a canonical constant volume ensemble using Langevin dynamics with the collision frequency gamma_ln = 2.0 and with position restraints on the solute molecules. The pressure was then equilibrated in five steps at 300 K in the NPT ensemble, during which a decreasing harmonic force constraint on the 146-glycans of 50, 25, 10, 5 and 1 kcal•mol -1 •Å -1 was applied at a target pressure of 1 atm. Then, 1 ns of NPT MD simulation without restraints was run before production runs at 300 K using Langevin dynamics with the collision frequency gamma_ln = 2.0. Finally, 100 ns MD simulations for each system were conducted at 300 K and at 1.0 atm. of pressure. Electrostatics were handled using the particle mesh Ewald (PME) algorithm with a 10.0 Å direct space non-bonded cutoff. All bonds involving hydrogen atoms were constrained using the SHAKE algorithm with a time step of 2.0 fs. The coordinates of the trajectories were saved every 10 ps during the entire MD run.

Docking
The protein structure (2HU0) was downloaded from the PDB Databank as a receptor. Autodock-vina was applied to dock the inhibitors into the receptor, with residue R156 being defined as flexible. Exhaustiveness was set at 16. The optimal docking conformation of inhibitor 20l was defined as a reference conformation. Finally, conformations that had the same binding mode as the reference conformation were selected for binding energy calculations for the remainder of the five inhibitors.

Binding energy calculations
The ligand binding energies and the interaction energies between the 146-N-glycan and the NA tetramer were calculated using MM/GBSA [30]. For calculations of the interaction energy, only the E vdw and E coul were considered. For calculations of the ligand binding energy, E MM and G solv were considered. The entropy contributions were neglected because the same receptor was used, and the normal mode calculations are computationally expensive and introduce errors that cause significant uncertainty in the results.
where E MM is a gas-phase MM energy. The free energy for each species (ligand, receptor, and complex) was decomposed into a gas-phase MM energy that contained polar and nonpolar solvation terms and an entropy term (which was neglected). E MM comprises E bat (the sum of bond, angle, and torsion terms in the force field), a Van der Waals term, E vdw and a coulombic term, E coul . G solv.p is the polar contribution to the free energy of solvation, often computed via the Generalized-Born (GB) approximation. G solv.np is the nonpolar free energy of solvation, usually computed as a linear function of the solvent-accessible surface area (SASA). A 2 ns MD simulation was performed on each complex. The final 50 snapshots with an interval of 10 ps were extracted to calculate the binding energy.

Volume calculation
The volume of the 150-cavity was calculated using the POVME1.0 program [31]. A 3D-grid with 0.5 Å spacing that covered the 150-cavity was created in the initial structure (see S11 Fig). Snapshots were extracted from the trajectory files every 40 ps and superimposed onto the reference structure (the initial crystal structure). Grid points overlapping with the protein atoms in the snapshots were systematically removed. The volume of the 150-cavity was then calculated by counting the remaining grid points.

Cluster analysis
To identify the structural representations of the 150-loop conformations that were generated by the MD simulations, cluster analysis was performed using cpptraj in AMBER 12 [20]. Snapshots were clustered every 10 ps using an average-linkage [32] clustering algorithm with an epsilon value of 2.0.
To reveal the relationship between the conformations of the 150-loop and 146-N-glycan, the conformers from the MD snapshots were projected onto the essential space (RMSD and dihedral) [33], and the molecular potential was given by Eq (2) [34].
where E represents the molecular potential, ρ is the conformer density and RT is the thermal energy. Three local minima of energies were selected and the relative conformations of 146-N-glycan were selected as representative cluster conformations. Finally, cluster analysis was performed.  Table. Activities of inhibitors of 6~12 to H5N1-1220 and H9N2-S2 (DOCX) S1 Text. Calculation and Statistics of interaction Energy between each glycan and the other part of complex (DOCX)