there is an urgent 
 need to better understand the molecular factors 
 that govern influenza a virus ( iav ) virulence .
seasonal ( human - adapted ) 
 iav is a major cause of morbidity and mortality each year , but the 
 pandemic potential of iav strains that originate in animal hosts poses 
 an even more serious threat to public health .
there have been four 
 influenza pandemics in the past century , in 1918 , 1957 , 1968 , and 
 2009 .
the precise origin of the 1918 pandemic strain , which killed 
 roughly 50 million people , is unclear .
however , the three subsequent 
 pandemic strains resulted from the process of genetic reassortment , 
 in which gene segments from different viruses mix during host coinfection 
 to produce novel viral strains . there is much concern that a new pandemic strain might result from 
 reassortment involving the highly pathogenic h5n1 and h7n9 avian influenza 
 viruses currently circulating in asia .
the majority of these viruses 
 undergo a deletion in the stalk of the neuraminidase surface protein , 
 a characteristic that is associated with virulence when causing outbreaks 
 in domestic poultry .
it is thus vitally 
 important to understand the effect of stalk deletion on viral function .
influenza a strains are named according to the alleles coding for 
 their two surface antigens , the hemagglutinin ( ha ) and neuraminidase 
 ( na ) glycoproteins .
hemagglutinin 
 ( ha , 18 antigenically distinct alleles identified to date ) facilitates 
 viral adhesion to human host cells by binding to cellular sialylated - oligosaccharide 
 receptors , and neuraminidase ( na , 11 antigenically distinct alleles ) 
 facilitates viral release by cleaving sialic acid linkages .
for example , hyperglycosylation 
 and active - site mutation tend to reduce ha / sialic acid binding . na 
 activity is generally proportional to the length of the highly variable 
 na stalk , and recent neuraminidase stalk deletions often produce highly pathogenic 
 avian influenza viruses such as h5n1 .
influenza virulence is determined in large part by the balance 
 between na and ha activities at the virion surface .
recent viral strains 
 with stalk - deletion na glycoproteins ( nadel ) are more virulent 
 because the reduced nadel activity alters this balance . 
 for example , the virulence of both a pandemic 2009 h1n1 virus and 
 a more recent 2013 h7n9 strain was enhanced when deletion mutations 
 containing 20 and 5 amino acids , respectively , were introduced into 
 the corresponding na stalks . in these strains ,
reduced na activity due to a stalk deletion may 
 enhance viral infectivity by reducing oligosaccharide - receptor cleavage 
 and thus increasing the number of receptors available for ha binding . 

this strategy is particularly advantageous when ha activity is itself 
 compromised ( e.g. , due to hyperglycosylation , active - site mutations , 
 etc . ) , leading many to express concerns that a current h5n1 strain 
 might obtain a long - stalk na through reassortment with a human - adapted 
 virus .
optimal viral replication is achieved only when the activities 
 of na and ha are ideally complementary .
the mechanism by which reduced 
 na stalk height impacts na activity 
 remains uncertain .
one prevailing theory suggests that nadel glycoproteins have reduced sialidase activity relative to wild type 
 ( nawt ) because their diminished height hinders access to 
 cellular sialylated - oligosaccharide receptors ( i.e. , the 
the 
 implication is that a towering canopy of ha effectively buries the 
 short - stalk nadel such that cell - bound oligosaccharides 
 simply can not reach the na enzymatic pockets .
while there is merit 
 to the limited - access theory , especially in cases where the stalk 
 length is profoundly reduced , experiment and computational modeling 
 suggest it may not fully explain the reduced nadel activity .
we hypothesize that , in addition to any limited - access effects , 
 na stalk length also reduces na enzymatic activity by altering the 
 binding affinity of sialic acid , the native substrate , to the enzymatic 
 pocket .
given that the stalk is roughly 40  from the active 
 site , we propose an indirect mechanism by which stalk length impacts 
 binding - pocket dynamics at a distance .
other studies have demonstrated 
 that the na enzymatic pocket is highly flexible and so samples many 
 conformations .
for example , x - ray crystallography and molecular dynamics 
 simulations of the isolated neuraminidase catalytic domain suggest 
 that key loops surrounding the active site ( i.e. , the 150- and 430-loops ) 
 are highly mobile , permitting large cavities to form immediately adjacent 
 to the main pocket .
given that the enzyme is known to adopt a variety of conformational 
 states , it is reasonable to hypothesize that stalk length may have 
 a quasi - allosteric impact on ligand binding affinity by shifting the 
 binding - pocket conformational population toward a subtly different 
 distribution . to properly study the impact of stalk length
, we built a full 
 model 
 of 2009 h1n1 neuraminidase as well as a model that included a 20 amino 
 acid stalk - deletion mutation similar to one known to enhance virulence 
 in a recent strains .
, who 
 replaced the long - stalk nawt of a 2009 h1n1 pandemic virus 
 with a short - stalk h5n1 nadel .
they attributed the resulting 
 decrease in transmissibility at least in part to a reduction in the 
 ability of the short - stalk nadel to cleave tethered substrates 
 and penetrate the host mucus barrier . rather than exactly replicating 
 the work of blumenkrantz et al .
, who inserted a foreign na ( from h5n1 ) 
 into a h1n1 virus , we made our deletion directly in the h1n1 nawt stalk itself .
our models were constructed using data from 
 both x - ray crystallography ( used to model the enzymatic head ) and 
 also protein  protein docking and modeling ( used to model the 
 stalk , transmembrane , and intravirion domains , which have never been 
 crystallized ) .
all na models were embedded in complex - mixture lipid 
 bilayers ( figure 1 ) .
equilibrated 
 neuraminidase models used in the current study : ( a ) 
 wild - type ( long - stalk ) nawt and ( b ) stalk - deletion nadel .
subsequent atomistic molecular 
 dynamics simulations of these systems 
 support the theory that the na stalk length impacts the dynamics of 
 the na binding pocket and therefore the affinity for the endogenous 
 substrate .
furthermore , these simulations reveal previously uncharacterized 
 druggable  hotspots  that may open up additional avenues 
 for structure - based drug design .
previous simulations have included 
 only the na head ; our models , which include additional na domains 
 as well as the associated lipid bilayer , allow us to take the next 
 steps toward more realistic atomic simulations of viral components
to generate a 
 complete model of long - stalk ( wild - type ) apo nawt , a monomer 
 of the tetrameric na head was built by homology using schrdinger s 

prime and the na sequence of the 2009 h1n1 pandemic reference strain 
 a / california/04/2009 ( genbank accession fj966084 ) , which codes for 
 a na protein identical to that of the a / england/195/09 virus used 
 by blumenkrantz et al .
default prime 
 parameters were used . the monomeric homology model , based on the 3nss 
 crystal structure , was missing a single 
 lysine residue on the c terminus , which was added manually . the appropriate 
 tetramer was formed by aligning copies of the homology - model monomer 
 to chains a  d of the 2hu4 crystal structure using multiseq , as implemented 
 in vmd .
the resulting tetrameric head 
 was then processed with schrdinger maestro s protein 
 preparation wizard and subjected to repeated minimization in the maestro 
 environment .
both sspro 
 4.5 , a neural - network - based algorithm 
 for predicting secondary structure from primary sequence , and abtmpro , a svm algorithm for determining whether a given 
 transmembrane segment is  helical or  barrel , suggested 
 that the stalk has an -helical structure .
we therefore used 
 the amber xleap module to generate a single  helix with the 
 appropriate long - stalk sequence .
rosettadock 3.4 was used to assemble four identical copies of this helix into a 
 four - helix bundle , one helix at a time . at each step ,
the top 200 best - scoring structures were 
 clustered by root - mean - square distance ( rmsd ) using cutoffs ranging 
 from 0.75 to 3.0  , and representative dockings were examined 
 visually .
the top representative docked pose that positioned the newly 
 introduced helix in the expected parallel conformation was retained . 
 when no such pose was identified by clustering , all docked poses were 
 examined sequentially , starting with the pose that had the highest 
 score , until an acceptable pose was identified .
the resulting four - helix 
 bundle was then processed with schrdinger maestro s 
 protein preparation wizard and subjected to repeated minimization 
 in the maestro environment .
we assumed that the cysteine residues located in the extra - virion 
 region of the stalk likely participate in interchain disulfide bonds , but that the cysteine residues located in the 
 transmembrane region likely do not .
the tetrameric - head and stalk / transmembrane - region models 
 were next positioned near each other using vmd .
the combined models 
 were subsequently loaded into maestro and again subjected to multiple 
 minimizations .
finally , maestro was used to build and position 
 the residues of 
 the short intravirion topological domain .
based on the assumption 
 that this domain lies at the interface of the inner membrane leaf 
 and bulk intravirion water , all hydrophobic side chains were manually 
 pointed upward toward the bilayer , and all polar side chains were 
 pointed downward . to generate a model of apo na with a virulence - enhancing 
 stalk 
 deletion (
hereafter called nadel ) , the appropriate mutant 
 amino acid sequence was first identified .
these residues were manually removed from the nawt model , and the flanking na segments were juxtaposed and fused using 
 vmd and schrdinger s maestro , followed by multiple 
 minimizations as needed .
our efforts 
 to construct a 
 complex - mixture lipid - membrane model were guided by a recent paper 
 describing the composition of the influenza viral envelope .
we identified and classified all listed membrane 
 components with mol % > 1.5 , taking care to keep the composition 
 of 
 the inner and outer leaflets separate .
as our ultimate goal was to 
 perform molecular dynamics simulations of this system , we were limited 
 to the lipids of the charmm 36 force field .
consequently , we had to substitute some experimentally identified 
 lipids with available lipid models .
3-palmitoyl-2-oleoyl - d - glycero-1-phosphatidylethanolamine ( pope ) was used for any phosphatidylethanolamine 
 ( pe ) variant , 2,3-distearoyl - d - glycero-1-phosphatidylserine 
 ( dsps ) for phosphatidylserine ( ps ) , 3-palmitoyl-2-oleoyl - d - glycero-1-phosphatidylglycerol ( popg ) for the forssman glycolipid 
 hapten , and 3-palmitoyl-2-oleoyl - d - glycero-1-phosphatidylcholine 
 ( popc ) for sphingomyelin ( sm ) . cholesterol and phosphatidylcholine , 
 also major components of the envelope membrane , are included in the 
 charmm 36 force field .
the na models were then positioned within this generated bilayer , 
 and clashing lipid residues were deleted with pymolecule . in preparation for molecular dynamics 
 simulations , we resolved multiple steric clashes in both the long - stalk 
 and stalk - deletion systems using serial iterations of minimization 
 and geometry optimization in schrdinger s maestro , 
 coupled with manual atomic and molecular manipulation in vmd .
partial 
 charges were 
 assigned to all atoms according to the specifications of the charmm22 all - hydrogen force field for proteins and the charmm36 all - hydrogen force field for lipids .
the vmd plug - in cionize(47 ) was used to position sodium and chloride ions as required to bring 
 the systems to electrical neutrality and to simulate a 20 nm solution . 

limited manual adjustments to the positions of some ions that had 
 been placed far from any protein or lipid molecule were made .
for both the nawt and nadel systems , amber s 
 xleap module was used to generate a water 
 box extending 15  beyond any atom .
water molecules were subsequently 
 removed if they were not positioned directly above the bilayer , so 
 that the 15  margin was ultimately retained only along the z 
 axis , perpendicular to the bilayer itself .
the nawt and 
 nadel systems had 366  293 and 351  699 atoms , 
 respectively .
psfgen 
 was used to fully 
 parametrize each system according to the charmm22 , charmm36 , and tip3p force fields for proteins , lipids , and water molecules , 
 respectively .
namd 2.9 running on the 
 stampede supercomputer at the texas advanced computing center was 
 used for all md simulations .
periodic boundary conditions were employed 
 with the particle mesh ewald method to account for electrostatic effects 
 ( smoothing cutoff : 14  ) . because of observed instabilities when 
 our standard minimization protocols were used , we subjected the initial 
 nawt system to only 1000 steps of unrestrained conjugate 
 gradient minimization , with the first and max initial steps for the 
 line minimizer reduced to 1.0  10 and 1.0 
  10 , respectively .
the structure was 
 minimized in four distinct steps : hydrogen atoms were first relaxed 
 for 5  000 steps ; hydrogen atoms , water molecules , and ions 
 were next relaxed for 5  000 steps ; hydrogen atoms , water molecules , 
 ions , protein side chains , and lipid tails were then relaxed for 10  000 
 steps ; and , finally , all atoms were relaxed for 25  000 steps . following minimization ,
each system was subjected to a constrained 
 equilibration using an npt ensemble at 310 k. the 
 equilibration was realized in four steps , using stepwise harmonic - constraint 
 force constants of 4 , 3 , 2 , and 1 kcal / mol /  on the 
 protein backbone . for each force
constant ( 1 fs time step ) , 250  000 
 steps of md simulation were performed .
langevin dynamics were applied 
 to maintain the temperature , and a modified langevin piston nos 
the nawt model underwent an additional 5 million 
 2 fs steps of unconstrained simulation , followed by another 21 
 million 1 fs steps .
the nadel model underwent an additional 
 5.5 million 2 fs steps of unconstrained simulation , followed 
 by another 36.5 million 1 fs steps .
satisfied that both systems 
 were sufficiently equilibrated , five 
 distinct 100 ns productive runs were performed for each ( 100 
 million 1 fs steps ) using distinct random seeds in order to sample 
 diverse protein configurations .
( 5 simulations )  ( 100 ns per simulation ) 
  ( 2 systems )  ( 4 monomers per system ) = 4 s of 
 monomeric simulation total .
5000 and 7000 frames were extracted from 
 the last 50% and 70% of the tetrameric nawt and nadel simulations , respectively , for subsequent analysis . for comparison purposes
, we also analyzed an archived simulation 
 of the isolated apo na head ( hereafter called nahead ) , 
 described previously . in brief , the na 
 head ( pdb i d : 3nss(26 ) ) was simulated for 100 ns using the 
 amber99sb force field , yielding 400 ns of monomeric trajectory .
for 
 the current study , we extracted 5000 frames from this simulation for 
 subsequent analysis .
we used root - mean - square - distance ( rmsd ) calculations 
 to judge 
 the geometric similarity between the primary na enzymatic pockets 
 of crystallographic holo and simulated apo models .
we first identified 
 active - site residues as those that came within 5.0  of a crystallographic 
 oseltamivir ligand ( pdb i d : 2hu4:a , h5n1 ) .
these residues are arg118 , glu119 , asp151 , 
 arg152 , arg156 , trp178 , ser179 , ile222 , arg224 , glu227 , ser246 , glu276 , 
 glu277 , arg292 , asn294 , tyr347 , gly348 , arg371 , and tyr406 . for the 
 purposes of our rmsd analysis , we discarded tyr347 because that tyrosine 
 is an asparagine in h1n1 na .
hydrogen atoms were also removed prior 
 to performing all rmsd calculations ; rmsd alignment was performed 
 only on the remaining heavy atoms .
we performed principal 
 component analysis ( pca ) on these same active - site heavy atoms .
after 
 rmsd aligning the nawt trajectories by these atoms using 
 vmd , we used the python modules scikit - learn and numpy to calculate their first two principal components 
 and to project each trajectory frame onto those components .
the frames 
 of the nadel trajectory were similarly projected onto the 
 same two nawt principal components to facilitate comparison .
the povme algorithm was used to measure the volume of the primary 
 sialic - acid - binding pocket over the course of the nawt , 
 nadel , and nahead simulations . to calculate 
 the volume of the binding pocket for a given trajectory frame , we 
 monitored the space within a pocket - centered sphere of radius 16.0 
  . in all cases ,
the pocket volume was calculated from the portion 
 of the sphere that was unoccupied by protein atoms . only volumetric 
 regions contiguous with the pocket
grid spacing and 
 padding parameters were set to 1.0 and 1.09  , respectively .
all - atom , residue root - mean - square 
 fluctuations ( rmsf ) were calculated for 
 each monomer of each simulation . for the nawt and nadel simulations ,
for the nahead simulation , 
 there were 4 associated monomeric simulations ( 1 run  4 monomers ) . 

in all cases , the monomeric trajectories were first aligned by the 
  carbons of residues ser90 to asp469 ( i.e. , the residues of 
 the na head )
. the residue center - of - mass rmsf values were calculated 
 for each monomeric simulation separately . to identify druggable 
 hotspots , we first identified representative protein conformations 
 from the md simulations . the heads of the aligned monomeric trajectories 
 were subjected to rmsd clustering , as 
 implemented in the gromacs computer package .
we adjusted the rmsd cutoff until the trajectory frames were separated 
 into five clusters .
the required cutoffs were 1.72 , 1.74 , and 1.74 
  for the nawt , nadel , and nahead simulations , respectively .
the centroid members of each of these 
 clusters were considered representative of the corresponding simulated 
 model .
the 15 representative na structures ( 3 models  
 5 cluster centroids ) were then submitted to the ftmap server to identify druggable hotspots . to simplify the analysis , all the 
 docked molecular probes associated with each simulation ( nawt , nadel , and nahead ) were considered collectively . 
 to identify druggable hotspots that were specific to a given model
, 
 the associated docked probes were superimposed , and any probe atom 
 associated with one model that came within 5  of any probe atom 
 associated with the other was deleted .
the remaining probe atoms congregated 
 at hotspots that were unique to their respective models .
as described 
 in the introduction , the limited - access theory 
 holds that stalk - deletion neuraminidase glyocproteins ( nadel ) have reduced sialidase activity relative to nawt because 
 their diminished height hinders access to cellular sialylated - oligosaccharide 
 receptors . the fact that na enzymatic activity varies depending on 
 ligand bulkiness is perhaps the most compelling evidence in favor 
 of this theory .
studies have shown that short - stalk nadel glycoproteins have significantly reduced activity when a bulky substrate 
 ( e.g. , fetuin ) is used , presumably because 
 the large substrate can not physically reach the na enzymatic sites . 

in contrast , na activity is largely stalk - independent when cleaving 
 the small molecule munana ( essentially a single sialic acid with a 
 4-methylcoumarin fluorophore attached to the 2-carbon ) , apparently because this smaller substrate is not subject to the 
 same large - scale steric hindrance .
however , several lines of 
 evidence suggest that the limited - access theory may not fully explain 
 the reduced nadel activity .
binding to 
 these receptors , which are extended , large , and cell bound , is likely 
 diffusion limited ; in contrast , munana is a relatively small molecule 
 with an affinity that is likely little dependent on kon .
while munana may be a fine substrate for measuring 
 relative na enzymatic activity in the context of a small - molecule 
 screen , it is arguably an inadequate surrogate for the actual receptor .
second , fetuin ( -2-hs - glycoprotein ) is also very different 
 than the endogenous receptor , and the theory that its bulk alone prevents 
 it from accessing the nadel enzymatic sites is hardly certain . 

we note that all available monoclonal antibodies known to target na 
 antigenic sites bind equally well to virions with nadel and nawt glycoproteins , despite 
 the fact that antibodies are much larger than fetuin ( 160 
 and 1000 kda for igg and igm , respectively , vs 40 
 kda for fetuin ) .
furthermore , 
 recent cryo - electron microscopy studies of whole - virion particles 
 demonstrate that the na glycoproteins are not evenly interspersed 
 among the more prevalent ha s as the limited - access theory 
 suggests ; rather , the na particles tend to cluster together on the 
 viral surface .
even when stalk mutations 
 substantially reduce the na height relative to ha , this clustering 
 effect alone should in theory preserve much access to the na enzymatic 
 head .
computational evidence likewise suggests that simple steric 
 considerations 
 are not solely responsible for reduced nadel activity . 

our models of 2009 h1n1 glycoproteins demonstrate that each stalk 
 amino acid contributes only 1.2  to the total na height 
 ( data not shown ) .
the 2013 h7n9 five - amino - acid stalk deletion with 
 enhanced virulence mentioned in the introduction therefore reduced the na height by only 6  ( roughly 
 4% ) .
more substantial stalk deletions 
 do not necessarily force the na enzymatic head to retreat entirely 
 behind a towering canopy of ha glycoproteins either . in our models 
 of h1n1
ha and na , we found the height of nawt was nearly 
 equal to that of ha ( 149 vs 154  , respectively ) .
when a virulence - enhancing 
 20 amino acid na stalk deletion was modeled 
 ( nadel ) , the height of the na fell to only 125  , 
 still 81% of the ha height .
given that the na head is itself a boxlike 
 structure 80  square , the resulting 24  height reduction 
 seems unlikely to substantially interfere with access .
we suggest 
 that , beyond limited - access effects , actual differences 
 in the affinities of nawt and nadel for the 
 endogenous sialylated - oligosaccharide receptors contribute to the 
 differences in enzymatic activity . as support for our reduced - affinity 
 hypothesis
, we analyze atomistic molecular dynamics simulations of 
 three na models : five 100 ns simulations of 2009 h1n1 nawt ; five 100 ns simulations of 2009 h1n1 na with a 20-amino - acid deletion 
 in the stalk ( nadel ) ; and 
 one 100 ns simulation of the isolated h1n1 na head , nahead ( figure 1 ) .
these 
 simulations indicate that stalk length does in fact impact the dynamics 
 of the enzymatic sialic acid binding pocket , suggesting a potential 
 mechanism for the proposed stalk - length - dependent differences in affinity .
to 
 verify that the simulations of the three na systems ( i.e. , nawt , nadel , and nahead ) had properly equilibrated , 
 we rmsd aligned each tetrameric trajectory by its c atoms using vmd .
the rmsd was then 
 calculated between each trajectory frame and a corresponding reference 
 structure ( figure s1 ) . despite extensive 
 pre - equilibration ,
the rmsd values of the whole - glycoprotein nawt and nadel simulations continued to drift slightly 
 throughout the productive runs .
to further investigate the source 
 of this drift , we calculated separate c rmsd values 
 for the residues of the na head ( figure 2 ) , the lower ( membrane - embedded ) stalk , and 
 the upper stalk .
the dynamics of the head and lower stalk were in 
 fact stable in both the nawt and nadel simulations . 

it was the upper stalk that was dynamically unstable . given that the 
 upper stalk serves as a linker connecting two fairly ridged bodies 
 ( the lower stalk and the head ) , it is reasonable to suppose that this 
 region is legitimately hypermobile .
rmsds of all trajectory monomers to their 
 respective reference 
 structures . in all cases ,
only the c of the head 
 domain ( i.e. , the c of all residues homologous to 
 the 388 residues present in the 3nss crystal structure ) were considered . 

( a ) nawt ( 5 tetrameric trajectories  4 monomers per 
 tetramer = 20 runs ) ; ( b ) nadel ( 5 tetrameric trajectories 
  4 monomers per tetramer = 20 runs ) ; and ( c ) nahead ( 1 tetrameric trajectory  4 monomers per tetramer = 4 runs ) .
subsequent analyses ( described 
 below ) focused primarily on the 
 residues of the na binding pocket , located on the enzymatic head . 

given the instability of the upper stalk , however , we chose to discard 
 the initial 50% and 30% of the nawt and nadel trajectories , respectively , to ensure that subsequent analysis focused 
 on the most equilibrated portion of the simulations .
in contrast , 
 nahead , based on the 3nss crystal structure , equilibrated far faster .
we therefore discarded 
 only the first 20% of the nahead simulation ( figures 2 and s1 ) .
we next 
 sought to characterize the effects of the stalk on binding - pocket 
 dynamics by comparing the geometries of our simulated apo binding 
 sites to that of holo ( ligand - bound ) na .
we calculated the rmsd values 
 between the simulated binding - pocket conformations and that of a reference 
 crystallographic structure , 2hu4 ( h5n1 ) , 
 which shares in common with h1n1 a total of 18 out of 19 binding - pocket 
 residues . in calculating these rmsd values , the single differing residue 
 was ignored .
a histogram of the rmsd values suggests that the 
 long - stalk nawt pocket exists in two conformational states 
 ( figure 3a ) , with rmsds 
 to the apo ( reference ) conformation of roughly 1.8 and 2.7  , 
 respectively .
in contrast , the pocket conformations sampled during 
 the nadel and nahead simulations were not binary , 
 suggesting the pocket is less flexible in these systems .
( a ) simulated apo binding - pocket 
 conformations are compared to the crystallographic ( 2hu4 ) holo conformation . 

geometric similarity was judged by calculating the rmsds between simulated 
 and crystallographic pocket residues .
. 
 geometric similarity was judged by calculating the rmsds between simulated 
 and crystallographic triad conformations .
( c ) representative 
 nawt arginine - triad conformations sampled from the 2.4 
  state ( in yellow ) and the 3.4  state ( in red ) .
the crystallographic 
 position of the triad residues is shown in thick , white licorice . 

an oseltamivir molecule has been positioned within the binding pocket 
 for reference ( in orange ) , though the simulations did not include 
 any ligand .
( d ) the simulated conformations of nawt and 
 nadel active - site heavy atoms , in black and red , respectively , 
 projected onto the first two principal components of the nawt simulation .
the nawt site assumes a greater variety of 
 conformations , in harmony with the rmsd analysis .
an arginine triad formed by arg118 , arg292 , and arg371 is 
 a key 
 feature of the enzymatic pocket .
this triad forms a positively charged 
 arginine nest that surrounds the negatively charged carboxylate groups 
 of the endogenous substrate as well as all fda - approved na inhibitors . 
 as with the whole pocket ,
the average rmsd between the atom positions 
 of the simulated apo triad configurations and the holo 2hu4 conformation were 2.5 and 2.4  for the 
 nawt and nadel simulations , respectively .
a t - test led us to reject the null hypothesis that these two 
 averages are statistically equal ( p = 0.0 ) .
a histogram of the arginine - triad rmsd values suggests that , unlike 
 the nadel triad , the nawt triad again assumes 
 two distinct conformations ( figure 3b ) .
to further characterize these two nawt states , we visually inspected representative conformations from 
 each ( figure 3c ) . in 
 the state with rmsd =  2.4  , arg371 is near its crystallographic 
 conformation , and arg292 is displaced only slightly in the direction 
 opposite of arg371 .
in contrast , in the state with rmsd = 3.4 
  , arg371 is generally displaced in the direction of arg118 , 
 and arg292 is even further from arg371 .
arg118 possesses roughly the 
 same conformation in both states . to further verify that the 
 nawt pocket assumes a greater 
 variety of conformational states than the nadel pocket
, 
 we next performed a principal component analysis ( pca ) of the nawt active - site heavy atoms ( figure 3d ) .
pocket conformations from both the nawt and nadel simulations were projected onto the 
 first two principal components . the nawt pocket ( in black ) 
 sampled the same region of pca space as the nawt simulation 
 ( in red ) , but the nawt pocket additionally sampled other 
 conformations not seen in the nadel simulation , confirming 
 that the nawt pockets adopted a greater number of conformational 
 states .
having studied the impact of stalk length on binding - pocket 
 dynamics generally , we next considered the dynamics of specific binding - pocket 
 residues .
we calculated the associated all - atom , residue root - mean - square 
 fluctuations for each of the 20 monomeric 
 trajectories associated with the nawt and nadel simulations . for each residue
, we compared the ensembles of nawt and nadel rmsf values using a t - test . a number of residues exhibited statistically significant average 
 rmsf differences when a liberal significance level was used ( t - test ,  = 0.10 ) .
most interestingly , the rmsf values 
 of two of the three arginine residues of the key active - site triad 
 ( arg371 and arg118 ) differed significantly ( table 1 ) .
this analysis again implicates the arginine 
 triad and 371 loop as contributing to differential , stalk - length - dependent 
 binding - pocket dynamics .
data from the 
 corresponding nahead ( stalk - absent ) simulations are also 
 included . the average 
 rmsf value is shown plus or minus the standard deviation .
the residue - specific 
 populations of all nawt and nadel rmsf values 
 were compared using t - tests .
to better visualize differences in active - site dynamics , 
 we colored 
 the residues of a nawt conformation by rmsf = rmsfwt  rmsfdel ( figure 4 ) .
the 371 , 406 , and 430 loops were generally 
 more flexible in the nawt simulation ( shown in red ) . however , 
 a portion of the 150 loop was more flexible in the nadel simulation ( shown in blue ) .
na active - site residues colored by rmsf 
 = rmsfwt  rmsfdel , from 0.2 
 to 0.2 .
residues in 
 red and blue were more flexible in the nawt and nadel simulations , respectively .
a crystallographic oseltamivir molecule 
 is shown for reference , though no ligand was included in the simulations .
finally , we also used our 
 molecular dynamics simulations to study 
 na druggability . to get an initial impression of the full extent of 
 pocket dynamics , we first performed a volumetric analysis of the primary 
 sialic - acid - binding pocket using the povme algorithm .
the trajectory frames with both the largest and smallest 
 volumes from the nawt , nadel , and nahead simulations were compared separately ( figure 5 ) . in all cases , the pocket fluctuated between 
 a closed conformation ( figure 5a ) and an open conformation ( figure 5b ) , suggesting a highly dynamic pharmacophore . 

interestingly , the nawt pocket tended to be larger than 
 the nadel pocket ( figure 5c , p = 0.0 per a t - test ) , implying possible differences in the pharmacophore distributions 
 of these two na variants that may have important implications for 
 drug discovery . to further study variant - dependent differences in 
 na - pocket druggability .
ftmap seeks to replicate in silico the experimental 
 multiple solvent crystal structures ( mscs ) method developed by mattos 
 et al . in brief , the ftmap algorithm uses computer 
 docking to flood the surface of a protein with small organic molecular 
 probes .
these probes are then clustered to identify regions of the 
 protein surface that are most likely to bind small molecules .
( a ) the smallest and ( b ) largest binding - pocket 
 conformations sampled by the nawt simulation .
the smallest 
 and largest conformations from the nadel and nahead simulations were similar ( data not shown ) .
( in licorice ) have been placed in each pocket for reference , though 
 the simulations included no ligand .
( ) histograms of the volumes sampled 
 over the course of the nawt and nadel simulations , 
 in black and gray , respectively .
we used ftmap to identify druggable hotspots on the surfaces 
 of 
 five representative na structures extracted from each of our molecular 
 dynamics simulations .
for the most part , the druggable regions of 
 the enzymatic site did not differ substantially between the nawt and nadel simulations .
however , ftmap did position 
 some molecular probes near a 430-loop - adjacent region of the nawt pocket that was not deemed  druggable  when 
 the nadel conformations were considered , suggesting subtle 
 differences that may be pharmacologically relevant ( figure 6a ) .
unique druggable hotspots 
 ( outlined ) . the 371 , 430 , and 150 loops 
 are shown in mauve , pink , and green , respectively .
( a ) a druggable 
 hotspot was identified near the 430 loop of the nawt pocket . 

( b ) a second druggable hotspot was identified beneath the 371 loop 
 of both the nawt and nadel pockets .
( c ) a large 
 druggable hotspot was identified near the base of the nawt and nadel heads .
the nawt and nadel hotspots were then 
 merged 
 and similarly compared to the nahead hotspots . surprisingly , 
 in both stalk - inclusive na simulations ( nawt and nadel ) , the 371 loop occasionally flipped upward , opening a novel 
 druggable hotspot that was not apparent when the representative structures 
 of the nahead simulation were considered ( figure 6b ) .
this pocket has been described 
 previously , but the nawt and nadel simulations confirm its importance . to date , no known na 
 inhibitor exploits this pocket , and it has not been the focus of previous 
 cadd efforts .
finally , another large druggable hotspot located at 
 the base of the na head was unique to the stalk - inclusive ( i.e. , nawt and nadel ) simulations ( figure 6c ) .
factors 
 that alter the sialic acid binding / association mechanism 
 may profoundly alter na affinity and activity in additional ways not 
 explained by the  theory of limited access . 
the length 
 of the stalk appears to influence the dynamics of the 150 , 430 , 371 , 
 and 406 binding - pocket loops .
the 371 loop and the critical arginine 
 catalytic triad are particularly affected . as a consequence
, 
 the apo nawt pocket may adopt a greater 
 variety of conformational states than does that of nadel .
specifically , the nadel pocket tends to sample a single 
 conformational state that is generally similar to that of the holo 
 conformation .
in contrast , nawt additionally samples more 
 distant pocket geometries ( figure 3 ) .
while the impact of stalk length on binding - pocket 
 dynamics may 
 well be partially responsible for the reduced nadel affinity , 
 we note also that stalk length likely has a profound influence on 
 the long - range electrostatics of the whole - virion surface , which could 
 also impact kon .
future efforts will focus 
 on elucidating stalk - length - dependent differences in this effect as 
 well . finally , the predicted druggable hotspots evident in the 
 stalk - inclusive 
 simulations also suggest novel and urgently needed opportunities for 
 drug discovery . prior to the 2009 h1n1 pandemic , seasonal h1n1 in 
 the united states was almost entirely oseltamivir resistant . fortunately , that pandemic strain itself remained 
 largely susceptible to oseltamivir , but 
 given the viral proclivity for rapid mutation and adaptation , multidrug 
 combinatorial therapies may be necessary in the future .
these novel hotspots also provide means 
 for experimentally testing 
 some of the hypotheses generated in the current study .
virtual screening 
 into the na pocket conformations sampled over the course of the nawt and nadel simulations could identify compounds 
 that are uniquely suited to the corresponding pocket geometries .
if 
 subsequent experimental validation confirms binding , x - ray crystallography 
 could be used to determine if the novel ligand(s ) stabilize one of 
 the predicted na conformations .
 deletions in the 
 stalk of the influenza neuraminidase ( na ) surface 
 protein are associated with increased virulence , but the mechanisms 
 responsible for this enhanced virulence are unclear . here 
 we use microsecond 
 molecular dynamics simulations to explore the effect of stalk deletion 
 on enzymatic activity , contrasting na proteins from the a / swine / shandong / n1/2009 
 strain both with and without a stalk deletion . by modeling and simulating 
 neuraminidase apo glycoproteins embedded in complex - mixture lipid 
 bilayers , 
 we show that the geometry and dynamics of the neuraminidase 
 enzymatic pocket may differ depending on stalk length , with possible 
 repercussions on the binding of the endogenous sialylated - oligosaccharide 
 receptors . 
 we also use these simulations to predict previously unrecognized 
 druggable  hotspots  on the neuraminidase surface that 
 may prove useful for future efforts aimed at structure - based drug 
 design . 
Introduction
Materials
and Methods
Results
and Discussion
Conclusion