Variable ligand- and receptor-binding hot spots in key strains of influenza neuraminidase.

Influenza A continues to be a major public health concern due to its ability to cause epidemic and pandemic disease outbreaks in humans. Computational investigations of structural dynamics of the major influenza glycoproteins, especially the neuraminidase (NA) enzyme, are able to provide key insights beyond what is currently accessible with standard experimental techniques. In particular, all-atom molecular dynamics simulations reveal the varying degrees of flexibility for such enzymes. Here we present an analysis of the relative flexibility of the ligand- and receptor-binding area of three key strains of influenza A: highly pathogenic H5N1, the 2009 pandemic H1N1, and a human N2 strain. Through computational solvent mapping, we investigate the various ligand- and receptor-binding "hot spots" that exist on the surface of NA which interacts with both sialic acid receptors on the host cells and antiviral drugs. This analysis suggests that the variable cavities found in the different strains and their corresponding capacities to bind ligand functional groups may play an important role in the ability of NA to form competent reaction encounter complexes with other species of interest, including antiviral drugs, sialic acid receptors on the host cell surface, and the hemagglutinin protein. Such considerations may be especially useful for the prediction of how such complexes form and with what binding capacity.

INTRODUCTION

The influenza A virus is a persistent public health threat that has the potential to cause extensive human and animal disease through both epidemic and pandemic events (Capua and Alexander, 2010; Medina and Garcia-Sastre, 2011). The two major glycoproteins on the surface of the influenza virus particle, neuraminidase (NA) and hemagglutinin (HA), have been well studied due to their roles as viral invasion machinery and major antiviral and vaccine targets, respectively. Yet, despite decades of investigation, many intriguing questions related to their structural dynamics and biophysical interactions during infection and treatment remain unanswered.

X-ray crystallographic structures provide critical information regarding the three dimensional structure(s) of the neuraminidase enzyme. However, they typically only provide one average snapshot of the protein among the ensemble of possible substrates that may be sampled. NA in particular has been shown to be an extraordinarily flexible enzyme, especially in the 150- and 430-loop regions (Russell et al, 2006; Amaro et al, 2007; Amaro et al, 2009). Although these two loops are not directly within the active site (i.e., where sialic acid (SA), the natural substrate of NA, binds and is cleaved), they line two other potentially important locations: the 150-cavity and the secondary sialic acid binding site (Laver et al, 1984). The composite “binding face” of NA is therefore comprised of the SA, secondary SA, and the 150- and 430-cavities (Figure 1). This region of NA is believed to complex with HA, host cell receptors that contain its natural substrate, and small molecule drugs. Four copies of the binding face would be exposed to incoming binding partners due to the tetramer oligomerization state of NA.

Figure 1.

NA binding face. A. The neuraminidase binding face is shown for a single monomer, with oseltamivir (yellow) bound in the active site. Blue indicates sialic acid binding pocket, green indicates 150-loop region, red indicates secondary sialic acid binding and 430-cavities. B. The tetramer configuration of the NA binding face, shown looking down at the top of the NA protein. C. Side view of the tetramer, with stalk region (green) and membrane (orange) indicated in cartoon.

The 150-cavity adjacent to the SA binding site has been recently shown to be accessible to small molecule compounds (Rudrawar et al, 2011). As some strains of NA, in particular the highly pathogenic avian flu H5N1 and the 2009 pandemic H1N1, have a high predominance of this cavity in their native structural ensemble, it has been proposed that this area may be targeted in the development of drugs that preferentially target those strains. This could be especially important for drug-resistant strains of N1, which periodically emerge in the population and are a constant public health threat. A more complete understanding of the binding site pockets available to different strains of NA may provide key strategic insights into the development of such compounds. Part of this understanding is evaluating the binding capacity of such newly revealed sites, which may be predictively assessed through computational solvent mapping experiments. In this work we employ FTMAP (Brenke et al, 2009) to carry out the solvent mapping experiments. Favorable binding regions of small organic probe molecules are determined via the following steps: (1) rigid body fragment docking using a fast Fourier transform approach, (2) minimization and rescoring of fragment-protein complexes, (3) clustering and ranking of low-energy fragment-protein complexes, (4) consensus site determination. Populated consensus sites found by FTMap have been shown to agree with ligand binding sites identified using experimental methods (Landon et al, 2007; Brenke et al, 2009; Landon et al, 2009).

In this work, we evaluate and compare the relative flexibility of three different N1 and N2 strains (A/Tokyo/3/67, A/Vietnam/1204/04, and A/California/04/2009) through all-atom molecular dynamics (MD) simulations, and investigate how such flexibility affects the binding site capabilities in different regions of the NA binding face, with particular emphasis on the 150-cavity and secondary sialic acid binding site area. Structural clustering performed on the residues that line the NA binding face provides information regarding the relative flexibility of this region among the strains. Subsequent computational solvent mapping experiments assess the capacity of these regions to bind to various interaction partners for NA, including host cell receptors, sialic acid molecules, drugs, and potentially HA.

MATERIALS AND METHODS

System setup

The system setup was performed as follows for all simulated systems (Table 1). Atomic coordinates were taken from 2HTY for A/Vietnam/1203/04 (VN04N1) (Russell et al, 2006), 3NSS for A/California/04/2009 (09N1) (Li et al, 2010), and 1NN2 for A/Aichi/3/67 (N2) (Varghese and Colman, 1991). Protonation states for histidines and other titratable groups were determined at pH 6.5 by the PDB2PQR (Dolinsky et al, 2004; Dolinsky et al, 2007) web server using PROPKA (Li et al, 2005) and manually verified. All crystallographically resolved water molecules and calcium ions were retained where possible. The system was setup using the program xLEaP from AMBERTools 1.5 (Case et al, 2010) using the AMBER99SB force field (Hornak et al, 2006). Disulfide bonds were enforced using the CYX residue notation in AMBER with the S-S bonds manually added in xLEaP. Each system was solvated in a orthorhombic box containing sufficient TIP3P (Jorgensen et al, 1983) water molecules to provide a minimum distance of 10Å between any solute atom and the edge of the box. Each system was neutralized by addition of Na+ or Cl- counter ions as appropriate and then additional Na+ and Cl- ions were added to reproduce experimental assay conditions of 20mM NaCl.

Table 1.

Description of simulated systems. System name, crystal structure PDB identifier, human strain isolate description, total simulation time for tetramer simulation, and total number of atoms, including solvent, are shown for each system.

System Name	Crystal Structure	Strain	Simulation Length (nsec)	No of Atoms
N2	1nn2	A/Tokyo/3/67	100	133,049
VN04N1	2hty	A/Vietnam/1203/04	100	112,311
09N1	3nss	A/California/04/2009	100	165,171

Molecular dynamics simulations

MD simulations were performed using a version of the PMEMD module from AMBER 11 that was custom performance tuned for each specific simulation and the NICS Cray XT4 and SDSC Trestles supercomputers. Each complex was minimized and equilibrated as follows: steric clashes caused by the addition of hydrogen atoms, water and ions were alleviated prior to performing molecular dynamics through a staged minimization protocol. Harmonic restraints, with an initial 5kcal mol-1 Å-2 force constant, on all non-hydrogen protein atoms, were slowly reduced over ∼40,000 combined steepest descent and conjugate gradient minimization steps.

Following minimization, the systems were linearly heated to 310K in the canonical NVT ensemble (constant number of particles, N; constant volume, V; constant temperature, T) using a Langevin thermostat, with a collision frequency of 5.0psec-1, and harmonic restraints of 4kcal mol-1 Å-2 on the backbone atoms. This was followed by three sequential 250ps long runs at 310K in the NPT ensemble, in which the restraint force constant was reduced by 1kcal mol-1 Å-2 each run. The pressure was controlled using a Berendsen barostat (Berendsen et al, 1984) with a coupling constant of 1ps and a target pressure of 1atm. A final equilibration was carried out with 250ps of NPT dynamics at 310K without restraints and a Langevin collision frequency of 2ps-1.

Production runs of 100ns were then conducted in the NVT ensemble at 310K. As with the heating, the temperature was controlled with a Langevin thermostat (but with a 1.0psec-1 collision frequency). The time step used for all stages was 2fsec and all bonds to hydrogen atoms were constrained using the SHAKE algorithm (Andersen, 1983). Long-range electrostatics were included on every step using the Particle Mesh Ewald algorithm (Darden et al, 1993) with a 4th order B-spline interpolation, a grid spacing of < 1.0Å, and a direct space cutoff of 8Å. For all trajectories, the random number stream was seeded using the wall clock time in microseconds. The production trajectories for each monomer of the tetramers were extracted and concatenated to approximate 400nsec of monomer sampling.

Clustering

RMSD-based clustering was performed identically for each strain using the clustering algorithms implemented in the rmsdmat2 and cluster2 programs of the GROMOS++ analysis software (Christen et al, 2005). Tetramer conformations were sampled at 200psec intervals yielding a total of 500 conformations. Monomer conformations were then concatenated together, giving 2,000 conformations for each strain. To remove external translation and rotation, an alpha-carbon atom RMSD alignment to the first sampled conformation of chain A was performed for each sampled monomer conformation. Following alignment, clustering was carried out using the GROMOS++ clustering algorithm (Daura et al, 1998; Baron and McCammon, 2007; Amaro et al, 2008), implemented in GROMACS (Lindahl et al, 2001; Hess et al, 2008), using a cutoff of 2.2Å on the alpha carbon atoms of the following 70 binding-site residues: 117 to 119, 133 to 138, 146 to 152, 156, 178 to 180, 196 to 200, 223 to 227, 243 to 247, 276 to 278, 293 to 295, 325, 346 to 350, 368 to 371, 403 to 406, and 426 to 441.

Computational fragment mapping

Computational fragment mapping was performed on the cluster centroids of each strain, using the free FTmap web service (FT-Map http://ftmap.bu.edu) (Brenke et al, 2009). The resulting output included the structures with ranked consensus sites (CSs). An automated algorithm to assess the differences among the ensemble of structures was implemented in a script that takes these CSs as input. Overlapping CSs were grouped with increasing distance until reaching a user-specified cutoff. Initially, the CSs in each frame were placed along the axes of a two-dimensional, symmetric matrix. The cells of the matrix were populated with the measured distance between the centroids of each CS. The lowest value in the matrix not located in the main diagonal was iteratively extracted, thereby identifying the two closest CSs. Rows that contain each of these CSs were merged, as are the columns. Afterwards, new distances were calculated between the centroid of the new CS and the rest of the CSs. This process iterates until the lowest distance in the matrix is higher than the user-specified cutoff. This is a form of average link agglomerative hierarchical clustering. The primary purpose for its design is to allow one to feasibly see the differences among the identities and types of probes that have been docked to each CS by highlighting hot spots that exist across several enzyme mutants. It also seeks to outline the differences in the number of probes docked into each consensus site as well as to discern preferences for fragments to a specific CS.

All protein surface images were rendered with MSMS (Sanner et al, 1995) through VMD (Humphrey et al, 1996).

RESULTS

Different strains of N1 and N2 exhibit varying degrees of flexibility

The three strains studied here: N2, VN04N1, and 09N1 (Table 1), exhibit varying degrees of flexibility in the region constituting the NA binding face (Figure 1). RMSD-based clustering on all of the atoms lining the binding face region was utilized as an indicator of flexibility. The results of this analysis indicate that the N2 strain was the most flexible, with a total of 17 clusters required to represent the structural ensemble sampled. 09N1 was the second most flexible, with 12 clusters, and the apo VN04N1 strain was the least flexible overall, requiring only 8 clusters (Table 2).

Table 2.

Cluster results from molecular dynamics simulations. RMSD-based clustering results are presented for each system. System name, number of clusters total, and the individual percentages of each cluster are listed. Highlighted cluster representative structures are depicted in the accompanying figures.

System Name	N2	VN04N1	09N1
No of Clusters	17	8	12
Cl1 %	44.4	87.9	41.3
Cl2 %	17.7	4.9	28.8
Cl3 %	16.4	4.2	17.5
Cl4 %	8.7	1.5	6.2
Cl5 %	5.2	1.3	2.2
Cl6 %	3.9	0.2	1.7
Cl7 %	0.8	0.1	0.9
Cl8 %	0.7	0.05	0.6
Cl9 %	0.5	//	0.3
Cl10 %	0.5	//	0.25
Cl11 %	0.5	//	0.25
Cl12 %	0.3	//	0.02
Cl13 %	0.2	//	//
Cl14 %	0.15	//	//
Cl15 %	0.15	//	//
Cl16 %	0.05	//	//
Cl17 %	0.05	//	//

Hot spots in N2

The N2 strain is among the most flexible based on the clustering algorithm, and exhibits three predominant clusters that represent 44%, 18%, and 16% of the sampled ensemble, respectively. These clusters have variable hot spots that highlight how such flexibility can impact ligand and receptor binding (Figure 2). In the most predominant cluster, although the 150-loop is in the closed conformation (indicated in green in Figure 2), there are shallow ligand-binding hot spots that persist.

Figure 2.

Neuraminidase hot spots. Ligand- and receptor-binding hot spots are shown for the sialic acid cavity (blue), 150-cavity (green), and secondary sialic acid and 430-cavities (red) for the most predominant N2 (A), second most predominant N2 (B), third most predominant N2 (C), most predominant VN04N1 (D), second most predominant VN04N1 (E), third most predominant VN04N1 (F), most predominant 09N1 (G), second most predominant 09N1 (H), and third most predominant 09N1 (I) ensemble structures. Clusters of organic probes indicating actual hot spot/probe binding locations are shown in various colors in stick representation.

Hot spots in VN04N1

The VN04N1 strain shows remarkable rigidity in the overall NA binding face region, as compared to both N2 and 09N1. In the dominant first cluster of this strain, the 150-loop is in an open conformation; 88% of the ensemble falls into this first cluster (Table 2). This open cluster exhibits ligand and receptor binding hot spots in the secondary SA site, SA site, and the 150-loop region implying many favorable ligand binding areas (Figure 2). Cluster 2, which represents only ∼5% of the structural ensemble, has a closed 150-loop and therefore no cavity or ligand-binding hot spot is exhibited, in marked contrast to the predominant conformation. Cluster 3, representing ∼4% of the trajectory ensemble, lacks a hot spot area in the lower 430-cavity, but exhibits hot spots in both the 150-cavity and the secondary SA binding site.

Hot spots in 09N1

The 2009 pandemic H1N1 strain exhibits an intermediate degree of flexibility, as indicated by the total number of clusters (Table 2). Hot spots from the first three clusters, which represent 41%, 29%, and 17% of the structural ensemble, respectively, are presented in Figure 2. While this strain exhibited a closed 150-loop configuration in the crystal structure, all-atom MD simulations indicated that it would open to adopt a loop conformation similar to the VN04N1 strain (Amaro et al, 2011). In the most predominant cluster conformation, there is no hot spot in the 150-loop region.

DISCUSSION

The current study, in which we analyze the structural dynamics and ligand-binding capacities of the NA binding face, presents several novel insights regarding recognition events in NA. All of the dominant conformations for N2, VN04N1, and 09N1 exhibit ligand-binding hot spots in the sialic acid binding site, the 150-cavity, and the secondary sialic acid binding site, and indeed, many of the binding features are conserved among the strains. For example, the sialic acid binding cavity exhibits strong conservation of probe moiety and locations (Figure 3). However, variability of protein structure among the three strains, as sampled through the MD simulations, results in different exposed surfaces and probe binding preferences for all critical binding regions. Even in the highly conserved sialic acid binding site, subtle structural changes at residue position 165 alter the binding pose of hydrogen bonding probes (Figure 4). In VN04N1, residue 165 is a serine, which provides a hydrogen-bonding moiety for probes. In N2, residue 165 is an alanine, yet probes are still able to form hydrogen bonds with the backbone oxygen of this residue. In 09N1, residue 165 is a proline, however, neighboring residue 166, a serine, provides compensatory hydrogen bonding interactions for binding partners. The ability of NA to vary its structure and residue content while maintaining similar functional binding features appears to be a hallmark of this enzyme. Indeed, the degeneracy in NA structural features and binding capacity may have evolved specifically to facilitate the modulation of receptor binding events.

Figure 3.

Sialic acid site preference for aromatic group conserved across neuraminidase strains. Surface representation for the hot spots within the sialic acid cavity are shown for the most predominant RMSD cluster of VN04N1 (A), 09N1 (B), and N2 (C). For clarity, only aromatic probes (benzaldehyde, benzene, phenol) are shown in licorice representation.

Figure 4.

Variations between sialic acid site preferences for hydrogen-bond donating groups. Surface representation for hot spots for the sialic acid cavity are shown for the most predominant RMSD cluster of neuraminidase strain VN04N1 (A), 09N1 (B), and N2 (C). Only hydrogen-bond donating probes are shown in licorice representation (acetamide, methylamine, isobutanol, ethanol, phenol, isopropanol, urea). Residue position 165, indicated with a yellow circle, varies among the strains. In panel (C), residue position 166 (a serine) is indicated with an orange circle.

Despite the fact that the 150-loop of the N2 strain does open in a small fraction of the trajectory, the 150-cavity is much more shallow, compared to VN04N1, which has a wide open 150-loop and a very deep cavity for ligand binding. The solvent accessible surface of this deep cavity exposes many highly favorable locations for ligand-binding, as determined through computational solvent mapping experiments (Figure 5). The additional 150-cavity depth presented by VN04N1 may allow for additional or modified interactions with host cell receptors or inhibitors. This finding may help rationalize how VN04N1 has been able to infect humans in rare cases. The deep cavity presented by the VN04N1 strain seems particularly amenable to binding polar functional groups, such as those present on the host cell glycan receptors. The 09N1 strain was recently reported to adopt an open conformation for the 150-loop over the course of 100 ns MD simulations (Amaro et al, 2011). Yet, shallow or deep 150-cavity hot spots are completely lacking in clusters 1 and 3 for the 2009 pandemic strain (Figure 2G, I). This may suggest reduced capacity of this strain to bind to more varied host cell receptors, as compared to the VN04N1 strain; although the second cluster structure for this strain does indeed indicate a deep pocket available for binding in the 150-cavity, similar to the VN04N1 strain. The relative population of this deep-150-cavity structure is much less (∼17%), compared to the VN04N1 deep cavities (present in over 90% of the ensemble), which may also influence such recognition events. A recent article highlighting the critical balance between NA and HA activity as related to human infection underscores the need to further investigate the structural dynamics at the NA binding face for the 09N1 strain (Yen et al, 2011).

Figure 5.

Polar probes bound within the deep VN04N1 150 cavity. Surface representation of the 150 pocket in VN04N1 with polar probes (acetamide, acetonitrile, acetone, acetaldehyde, methylamine, benzaldehyde, isobutanol, NN-dimethylforamide, dimethyl ether, ethanol, phenol, isopropanol, urea) bound at the bottom of the site. Surface colored by electrostatic potential as calculated by DelPhi using the VMD plugin DelEnsembleElec.

The secondary SA binding site was consistently predicted to favorably bind solvent probes in all of the cluster representative frames investigated, across all of the NA subtypes studied. This site exhibits a high degree of structural variability, which is reflected and highlighted by the subtle differences in the identity and arrangement of the bound probes (Figures 2 and 6). In VN04N1, residue number 403 is an arginine, whereas in 09N1 and N2, residue 403 is a tryptophan. The loss of the positively charged group impacts the electrostatic binding properties of the secondary site, and correspondingly, the features of the solvent probes that are predicted to favorably bind to the area. The consistent prediction of this region to favorably bind solvent probes indicates its high propensity to form complexes with ligands. It therefore seems likely that this region participates in the host cell sialic acid receptor recognition events. Experimental evidence for the existance of this secondary sialic acid binding site has recently been presented using saturation-transfer-difference (STD) NMR spectroscopy (Lai et al, 2012). These studies validate the hot spot predictions at the secondary site, as presented in this work for the 2009 pandemic H1N1 strain. Unfortunately, the exact in vivo glycosaccharide composition of these receptors and their configurations is not well understood; such receptors are likely to be highly complex in terms of glycosaccharide content and include branched topologies, which makes predictive models of so-called “encounter complexes” more challenging. We propose that the information provided here could be utilized in complex modeling studies, in order to push the frontiers of understanding molecular recognition events and function of these critically important enzymes.

Figure 6.

Variations between secondary sialic acid binding sites. Surface representation for hot spots in the sialic acid cavity are shown for the most predominant cluster of VN04N1 (A), 09N1 (B), and N2 (C). Surface is colored by residue type: blue indicates positive residues; red, negative; green, polar; white, hydrophobic.

Considering the structural ensemble for NA enzymes reveals important insights relevant to its function. Analysis of X-ray crystal structures alone may be insufficient for a full understanding of these dynamic enzymes.

The 150-cavity exhibits varying topology and frequency among the three strains studied here: N2 from Tokyo in 1967, the highly pathogenic H5N1 from Vietnam in 2004, and the 2009 pandemic H1N1 from California; ligand-binding hot spots vary correspondingly.

The 09N1 strain exhibits an open 150-loop, but the 150-cavity it presents is much more shallow than VN04N1, which has a deep cavity. Both shallow and deep cavities exhibit persistent ligand-binding hot spots for these two strains.

The favorable binding of solvent probes to the secondary sialic acid binding site persists throughout the structural ensembles for all of the strains studied, despite the large degree of structural variability sampled in the trajectories. This suggests that this secondary site may play an important role in the complexation of NA with sialic acid receptors on the host cell.