2009 
 saw the first influenza pandemic of the 21st century , when 
 an h1n1 influenza a virus spread to over 200 countries , resulting 
 in more than 18  000 deaths globally .
the wide spread 
 of the virus resulted from its efficient human - to - human transmission .
transmission requires binding of the virus surface envelope protein , 
 hemagglutinin ( ha ) , to the host surface glycan receptors containing 
 terminal sialic acid ( sa ) .
the h1n1 human flu virus preferentially 
 binds to receptors with a 2,6 linkage between sa and galactose 
 ( gal ) . these 2,6-linked sialic acids ( 2,6-sa )
are predominantly 
 found in the epithelial cell surface of the human upper respiratory 
 tract .
in contrast , the has of avian - adapted virus preferentially bind 
 to 2,3-sa , which are found in the digestive and respiratory 
 tracts of birds , and on human lung alveolar cells . although the 
 overall case fatality rate of 2009 h1n1 was not much different from 
 seasonal flu , a specific mutation in ha at position 222 ( h1 numbering ) 
 from aspartic acid to glycine ( d222 g ) was more likely to be associated 
 with severe or fatal cases .
the presence of this d222 g mutation in both the ha from 1918 h1n1 
 and 2009 h1n1 viruses was reported to enhance ha binding to 2,3-sas , 
 potentially causing the severe cases observed early in the 2009 pandemic .
the results also showed that the d222 g mutation significantly weakened 
 binding to 2,6-sa for ha from 1918 h1n1 , but only modestly 
 affected the 2009 h1n1 virus . to explain this discrepancy , 
 we previously performed computer simulations 
 for the 2,6-sa analog 6sln binding with has from the 
 1918 h1n1 ( a / south carolina/1/18 , sc18 ) and 2009 h1n1 ( a / netherlands/602/2009 , 
 nl602 ) viruses with and without the d222 g mutation . for wild type sc18 ( sc18-wt ) , two binding modes for
these 
 modes differ in the interactions formed between 2,6-sa and 
 ha , which alter the position of the galactose sugar within the receptor 
 binding pocket .
structures with shorter distances between the ligand 
 and the residue at the 224 position are classified as mode 1 , and 
 those with longer distances are identified as mode 2 . by definition
, 
 the ligand in mode 1 is buried deep inside the pocket , while in mode 
 2 it is situated higher up . nl602-wt was also observed to adopt two 
 similar binding modes alongside many intermediate structures . for 

sc18 carrying the d222 g mutation ( sc18d222 g ) , almost no binding 
 mode 2 was observed , while for nl602d222 g , although a decrease 
 in binding mode 2 was seen , it was not entirely absent .
this contrasting 
 stability is a result of the additional interaction between 2,6-sa 
 and the side chains of 130k , 142k , and 224e , which can stabilize 2,6-sa 
 binding to ha .
we suggested that these interactions allow nl602d222 g 
 to maintain binding to 2,6-sa , explaining the fitness of this 
 mutant virus strain .
the binding pocket region distal 
 to the 
 viral membrane is highlighted for one monomer .
( b ) schematic illustration 
 of the two binding modes ; the 6sln glycan in mode 1 and mode 
 2 is colored blue and green , respectively . the distance between the 
 -c of 224a and the o3 atom of gal ( denoted as d224 )
other efforts have also 
 been made to understand how the d222 g mutation 
 changes the binding preference , in both experiments and computer simulations . 
 in 2013 ,
zhang et al . reported the crystal 
 structures of ha from both h1n1 a / california/04/2009 and sc18 flu 
 virus bound to 2,6-sa and 2,3-sa analogs . 

they suggested that the d222 g mutation results in a missing salt bridge 
 between 222d and 219k , loosening the 220-loop and enabling the key 
 residue 223q to interact with the avian receptor , switching the binding 
 preference . using nmr and molecular dynamics ( md ) simulation , elli 
 et al . also reported stable structures of 2,6-sa 
 binding to ha from 1918 flu before and after d222 g mutation .
the results 
 suggested that the interaction between 2,6-sa and the 220-loop 
 in wild type ha no longer exists after the mutation , which allows 
 the galactose sugar to move away from the protein surface .
also performed md simulations 
 investigating the binding of sialyldisaccharides to various has . using 
 neu5ac(26)gal as the analog of 2,6-sa and ha from the a / swine / iowa/30 
 virus , they also observed two binding 
 modes . due to the difference of the analog and ha model , these two 
 binding modes are not exactly the same as those reported in 2010 , most likely due to differences in the ha model 
 and glycan analogue employed .
structures of all initial end points . the names 
 of key residues 
 in sc18-wt binding mode 1 are labeled . the previous work described above mainly focused on the differences 
 between static structures
having identified alternative binding modes , one can 
 investigate not only their kinetic stability , both in terms of on / off 
 rates and the free energy of binding , but also how readily they interconvert . 

this information , especially the transition mechanism , the barriers , 
 and the key intermediates , could be helpful in identifying additional 
 therapeutic targets and provide a more complete understanding of the 
 effect of a particular mutation on viral fitness . in the present work , 
 we investigate interconversion pathways between two binding modes 
 to elucidate the contrasting affinity of the two different h1n1 virus 
 mutants for 2,6-sa .
a new intermediate network scheme is introduced 
 to accelerate the search for kinetically relevant paths in the potential 
 energy landscape . by comparing the distribution of energy barriers 
 and energy changes of key rearrangements ,
we find that the d222 g mutation 
 exchanges the energetic preference of the two binding modes in both 
 sc18 and nl602 has and that the effect is more dramatic in sc18 . by 
 analyzing the shape of the binding pocket
, we have identified differences 
 between the two binding modes for the 1918 and 2009 strains .
our results 
 indicate that changes at positions 183 and 224 generate a smaller 
 pocket in nl602 than in sc18 , permitting greater flexibility in the 
 2,6-sa transformation between the two modes in nl602 .
the amber ff99sb molecular mechanics force field in conjunction with 
 the gb generalized born implicit solvent method ( igb = 
 2 ) was used to represent the protein . 

the potential was modified to remove a discontinuity at the nonbonded 
 interaction cutoff , which was set to 16  .
the trisaccharide , 
 neuac2,6-gal14glcnac1 ( 6sln ) , 
 with a methyl cap attached to the o1 atom of the terminal sugar , was 
 used to represent 2,6-sa in our simulations , since it has been 
 seen experimentally to bind to a range of h1n1 flu viruses .
parameters for the three sugars that comprise this analog ( neuac , 
 gal , and glcnac ) were obtained from the glycam06 g force field , as in our previous study , where the binding 
 modes were first identified .
the initial 
 structure of the ha in nl602-wt was based upon the crystal structure 
 in the sc18 h1n1 .
the order parameter 
 used as the criterion to choose initial structures was the separation 
 of the center of mass of the backbone atoms of residue 224 and the 
 center of mass of the o2 and o3 atoms of the 2,6-sa galactose 
 sugar .
this order parameter was chosen because it was found to demarcate 
 two distinct binding modes , which had been observed in md snapshots 
 from earlier work .
once these snapshots 
 had been sorted into one of these binding modes ( or intermediate structures ) , 
 we selected representatives based on the following criteria : ( i ) structures 
 temporally separated from an observed transition between modes and 
 ( ii ) structures clearly assigned as mode 1 or mode 2 .
these choices 
 ensure that the structures are indeed representative , prior to minimization . for sc18-wt and nl602-wt ,
we then manually interpolated between the two end points , 
 removing those differences in binding mode 2 far away from the pocket , 
 which do not affect interconversion behavior . for sc18d222
g , 
 binding mode 2 was entirely absent from the md trajectory reported 
 previously .
hence , in the current work , 
 a structure for mode 2 was generated using the ligand from sc18-wt , 
 superimposing it into the pocket of sc18d222 g , incorporating 
 residue moves around the pocket corresponding to binding mode 2 , and 
 finally performing a local minimization . for nl602d222 g , although 
 some of the structures seen in previous work exhibit features of binding 
 mode 2 , the total number is very small . 
 in the present work , binding mode 2 was also generated using the ligand 
 from nl602-wt mode 2 , superimposing it into the binding pocket , then 
 minimizing .
a detailed description can 
 be found on our website and in the supporting information .
the 
 doubly nudged elastic band ( dneb ) method , as implemented 
 in the optim program , was used to generate 
 likely candidate structures for transition states along the pathway 
 between a pair of end point structures .
the two minima that each transition state connects were identified 
 by calculating approximate steepest - descent paths .
all the minima 
 and transition states characterized during the connection of each 
 pair of end points were collected in the initial database generated 
 by the intermediate network scheme , as described in the following 
 section .
the databases were then expanded using the discrete path 
 sampling ( dps ) framework , with several shortcut , shortcut 
 barrier , and untrap runs , as implemented in the 
 pathsample program .
shortest - path algorithm with appropriate 
 edge weights was used to determine the discrete 
 path that makes the largest contribution to the rate constant between 
 the end points of interest ( the fastest path ) . efficiently 
 finding
a connected pathway between two significantly different end 
 points can be very difficult if attempted in a single step due to 
 limitations in how intelligently we can interpolate between them . 

we found that the main differences between the two end points are 
 mainly rotations of 6sln and the side chains of residues on 
 the edge of the binding pocket .
the main body of the complex , including 
 the position of the center of mass of 6sln in the pocket , 
 is almost the same .
this observation inspired a new intermediate network 
 interpolation technique to optimize refinement of the kinetic transition 
 network and ensure that we have included as many paths as possible 
 in the initial database .
this scheme is illustrated in figure 3 and consists of the following steps:(1)for two end points a and b , identify all the differences of local 
 conformations and 
 label these differences as ai and bi , ( i = 1 , n ) .
the end point a can then be denoted as the set a = { a1 , a2 , ... , an}. for example , in figure 3 , the molecular fragment in question is represented by the 
 blue square , which should be similar for a and b. four differences are identified and labeled with small 
 circles , which are white in a and red in b.(2)new configurations 
 are generated for 
 end point a by transplanting the different local b conformations into the original a minimum . 

we adopt the notation i to denote the 
 minimum obtained by relaxing the structure after the local conformation 
 replacement , where  is a binary string , identifying the origin 
 of each local configuration . those originally from a and b are denoted as  0  and  1  , 
 respectively .
the leftmost digit represents a1 or b1 . for example , in figure 3 , i0110 = { a1 , b2 , b3 , a4 } corresponds to the minimum 
 obtained by relaxing the structure with the b2 and b3 conformations replacing a2 and a3 . hence i00
= a. we note 
 that i11 can differ in detail 
 from b because the local relaxations after transplanting 
 conformations can produce minor changes in nearby side chains . in 
 principle
the total number of ia is 2 , including a. however , some minima obtained after relaxation may be the same , 
 because certain substitutions may have to follow a particular sequence 
 to be favorable.(3)double - ended 
 connection attempts are 
 then run for all the possible i and i corresponding to a single local conformation 
 change between  and . here we require the pathway search 
 algorithm to locate multiple transition states , ensuring that any 
 given intermediate pairs will eventually be connected , as implemented 
 in the  missing connection  approach , which is coded in optim . in 
 figure 3 , we have shown several possible connections . 
 in principle
, the total number of connections covered by this network 
 could be n  2 + 1 , where the extra + 1 corresponds to the case when i11  b. the number will be less than this upper 
 bound if 
 some of the minima generated in step 2 are the same . for two end points a and b , identify all the differences of local 
 conformations and 
 label these differences as ai and bi , ( i = 1 , n ) .
the end point a can then be denoted as the set a = { a1 , a2 , ... , an}. for example , in figure 3 , the molecular fragment in question is represented by the 
 blue square , which should be similar for a and b. four differences are identified and labeled with small 
 circles , which are white in a and red in b. new configurations 
 are generated for 
 end point a by transplanting the different local b conformations into the original a minimum . 

we adopt the notation i to denote the 
 minimum obtained by relaxing the structure after the local conformation 
 replacement , where  is a binary string , identifying the origin 
 of each local configuration . those originally from a and b are denoted as  0  and  1  , 
 respectively .
for example , in figure 3 , i0110 = { a1 , b2 , b3 , a4 } corresponds to the minimum 
 obtained by relaxing the structure with the b2 and b3 conformations replacing a2 and a3 .
hence i00 = a. we note 
 that i11 can differ in detail 
 from b because the local relaxations after transplanting 
 conformations can produce minor changes in nearby side chains . in 
 principle
the total number of ia is 2 , including a. however , some minima obtained after relaxation may be the same , 
 because certain substitutions may have to follow a particular sequence 
 to be favorable .
double - ended 
 connection attempts are 
 then run for all the possible i and i corresponding to a single local conformation 
 change between  and . here we require the pathway search 
 algorithm to locate multiple transition states , ensuring that any 
 given intermediate pairs will eventually be connected , as implemented 
 in the  missing connection  approach , which is coded in optim . in 
 figure 3 , we have shown several possible connections . 
 in principle
, the total number of connections covered by this network 
 could be n  2 + 1 , where the extra + 1 corresponds to the case when i11  b. the number will be less than this upper 
 bound if 
 some of the minima generated in step 2 are the same .
schematic illustration of the intermediate network scheme . the 
 blue square represents end points a and b. four small circles represent the variable fragments in a and b , which are white ( a1a4 ) in a and red ( b1b4 ) in b ( b  i1111 ) .
we consider this procedure to be converged when pathways 
 have been 
 found between all of the states i. one advantage of the 
 scheme is that all the minimizations in step 2 and the connection 
 attempts in step 3 can be run in parallel independently .
the corresponding 
 pathways can have multiple transition states with additional local 
 minima stabilized by nonlocal interactions .
however , they are precisely 
 the paths that are likely to have the lowest barriers and make the 
 largest contribution to the overall rate constant . for the ideal case 

that i11 = b , and any pair of intervening minima have only one possible connection , 
 corresponding to a single transition state , then in principle , this 
 network will cover all the possible paths that connect the two end 
 points . in practice , we may have to group more than one residue into 
 a single fragment so that one single connection in step 3 still contains 
 several elementary steps and alternative possibilities .
hence , it 
 is still worthwhile to further refine the initial database that includes 
 all the minima and transition states from the intermediate network 
 connection searches .
the glycan binding 
 pocket is located on the top of each monomer of 
 the ha trimer . as shown in figure 1b , it consists 
 mainly of three secondary structure features : the 190-helix ( residues 
 187 to 195 ) , 220-loop ( residues 218 to 225 ) , and 130-loop ( residues 
 132 to 135 ) .
the aspartic acid residues located at the 222 and 187 positions 
 contribute most to the binding of 2,6-sa . by definition ,
the 
 fastest path makes the largest contribution to the rate constant between 
 the end points of interest .
analyzing this path helps us to understand 
 how displacements of individual residues contribute to the overall 
 interconversion mechanism between the two binding modes .
the fastest 
 paths in each of the four strains all involve loss of the interaction 
 between neuac and 187d and formation of a new hydrogen - bond between 
 gal and the 220-loop .
several motions are shared by the fastest paths 
 of all four strains , namely ( 1 ) the rotation of the 4-hydroxyl of 
 gal toward the backbone of 222d / g ( g4 movement ) , ( 2 ) the 
 3-hydroxyl of gal binding to the side chain of 222d of wild type has , 
 or moving to where it would have bound to the side chain of 222d in 
 mutant has with 222 g ( g3 movement ) , and ( 3 ) the rotation 
 of the side chain of 187d toward the outside of the binding pocket 
 ( r187 movement ) . for the fastest path in sc18-wt , 
 we can distinguish four sets of movements .
neuac first loses interactions 
 with 187d by rotating the 7-hydroxyl and forms new intramolecular 
 hydrogen - bonds , after which the g4 and g3 movements 
 occur in succession .
the 3-hydroxyl of gal points toward the 2-hydroxyl 
 of gal in the binding mode 1 , which is unique to the sc18-wt strain . 

these two steps together involve a highest barrier of 7.8 kcal / mol 
 and make the largest contribution to the distance change between the 
 -c of 224a and the o3 atom of gal ( denoted as d224 ) , which distinguishes the two binding modes .
the final 
 change is the r187 motion , where 187d loses most of its 
 interactions with the ligand , leaving only one hydrogen - bond with 
 glcnac .
further details of the fastest path can be found in the supporting information . in sc18d222 g , 
 the position of the ligand is first changed 
 by the rotation of the 8-hydroxyl of neuac , which breaks its hydrogen - bond 
 with 187d and forms new hydrogen - bonds within neuac itself .
this step 
 has a high barrier of 14.1 kcal / mol and is endothermic by 11.7 
 kcal / mol .
next are the g4 and 
 g3 movements , which are similar to those in sc18-wt , except 
 that the 3-hydroxyl of gal interacts with the backbone of 222 g in 
 binding mode 1 .
the g3 movement still results in the largest 
 change in d224 but the corresponding energy 
 changes are much smaller than those in sc18-wt .
the first 
 steps involve breaking the hydrogen - bond 
 between 224e and the 2-hydroxyl of gal , followed by conformational 
 changes of the 224e side chain .
this step does not exist in sc18 since 
 the 224a side chain in sc18 is too short to interact with the ligand . 

the 9-hydroxyl of gal then rotates , forming a hydrogen - bond with 224e , 
 which shifts the whole 220-loop toward the ligand .
the g4 , g3 , and r187 movements that follow have two 
 interesting differences from previous pathways : ( 1 ) for both sc18-wt 
 and sc18d222 g strains , we observed translational motion of 
 the ligand toward the 220-loop when the g4 and g3 movements happen , which is not seen in nl602-wt .
( 2 ) after the r187 motion , the ligand in sc18-wt and sc18d222 g moves 
 away from the base of the pocket , losing interactions with 187d , while 
 in nl602-wt , most of the hydrogen - bonds remain after the r187 displacement .
the r187 displacement 
 comes first , together with hydrogen - bond breaking between 
 224e and the 2-hydroxyl of gal . as discussed for nl602-wt ,
then , following the g4 and g3 movements , d224 increases 
 from 5.2 to 6.4  .
the final two steps correspond to motion of 
 residues 180 h and 191l toward the base of the pocket . by comparing the 
 fastest path that connects the two binding modes 
 in each case , we find that the shift of the 3-hydroxyl of gal , either 
 binding to the side chain of 222d of wild type has , or moving to where 
 it would have bound to the side chain of 222d in mutant has with 222 g 
 ( the g3 movement ) , is the key elementary step differentiating 
 the two binding modes
although the g3 displacement is common 
 to all four cases , the fastest transition paths exhibit systematic 
 differences . in sc18-wt
, the 3-hydroxyl does not bind to the backbone 
 first before binding to the side chain of 222d . in nl602-wt
this change , 
 and the binding of the 4-hydroxyl of gal to the backbone of 222d , 
 occur simultaneously . to compare
further we collect all the transition 
 states corresponding to steps where the 3-hydroxyl of gal either binds 
 to the side chain of 222d in wild type has , or moves to the place 
 where it would have bound to the side chain of 222d in mutant has 
 with 222 g .
the barrier and energy change distributions for the four 
 cases are illustrated in the contour plots shown in figure 4 . in sc18-wt ,
the peak value of the distribution 
 is around ( 0.7 , 1.7 ) , corresponding to an energy barrier of 0.7 kcal / mol and an energy difference 
 of 1.7 
 kcal / mol .
the pathways that contribute to this peak all involve breaking 
 the hydrogen - bond between the 3-hydroxyl of gal and the backbone of 
 222d .
the transition states that contribute to the shoulder peak around 
 ( 0.5 , 2.5 ) are broadly similar , with insignificant differences 
 in nearby side chains , for example , the position of the 7-hydroxyl 
 of neuac .
the step in which the 3-hydroxyl of gal breaks hydrogen - bonds 
 with non-222d residues mainly contributes to the small peak at ( 3.9 , 
 1.2 ) , as we observed in the fastest path .
contour plots illustrating 
 the distribution of the barriers and 
 the energy changes corresponding to the 3-hydroxyl of gal , either 
 binding to the side chain of 222d in ( a ) sc18-wt and ( c ) nl602-wt , 
 or moving to where it would have bound to the side chain of 222d in 
 ( b ) sc18d222 g and ( d ) nl602d222 g with 222 g .
the x axis corresponds to the height of the barrier and the y axis corresponds to the energy difference . in sc18d222 g
, 126 out of 3475 transition 
 states in the 
 database correspond to the key elementary step .
the main peak shifts 
 to ( 0.6 , 0.5 ) , which still mainly correlates with interactions of 
 the 3-hydroxyl of gal and the backbone of 222 g in the initial state . 

the energy barrier remains almost the same as in the wild type ha , 
 while the overall energy difference increases by at least 2.2 kcal / mol , 
 as a result of the d222 g mutation .
the small peak located at around 
 ( 2.8 , 2.3 ) corresponds to cases where the 3-hydroxyl of gal is not 
 bound to the backbone of 222 g in the initial state .
the structure 
 of the corresponding initial state is very similar to the one reported 
 by zhang et al .
, where the 3-hydroxyl 
 binds to 219k and the 4-hydroxyl of gal binds to the backbone of 222 g . in nl602-wt ,
332 out of 3719 transition states in the database 
 correspond to the 3-hydroxyl motion .
the main peak in figure 4c is around ( 0.1 , 3.4 ) and represents paths 
 where the 3-hydroxyl breaks the hydrogen - bond to the backbone of 222d , 
 while simultaneously binding to the side chain of 222d .
compared to 
 sc18-wt , this step is significantly faster and more exothermic in 
 nl602-wt .
the same mechanism also contributes to the broad peak ranging 
 from ( 2.7 , 5.5 ) to ( 3.6 , 4.9 ) , and in these paths , 
 the 4-hydroxyl of gal is already bound to the backbone of 222d in 
 the initial state .
the 4-hydroxyl subsequently binds to the backbone 
 of 222d , as in the fastest path .
another small peak located at around 
 ( 0.1 , 1.0 ) corresponds to paths where the 3-hydroxyl of gal 
 does not interact with the backbone of 222d in the initial state . in nl602d222 g ,
only one peak located at 
 ( 0.3 , 0.1 ) can be seen in figure 4d , corresponding 
 to paths with the 3-hydroxyl of gal bound to the backbone of 222 g 
 in the initial state .
for wild type has , we see that binding 
 the 3-hydroxyl of gal to 
 the side chain is generally an energetically favorable process .
when 
 the hydroxyl shifts from interacting with the backbone to the side 
 chain , the barrier is always lower than 1 kcal / mol . as a result , mode 
 2 corresponds to an energy minimum in this transformation , increasing 
 the equilibrium population and making it accessible during md simulations , 
 as observed in previous work . for 
 mutant
has , mode 2 is generally energetically unfavorable compared 
 to mode 1 . in particular , for sc18d222 g mode 2 is so unfavorable 
 that the reverse process ( mode 2 to mode 1 ) has a barrier of only 
 0.1 kcal / mol . according to the energy profile , other processes , such 
 as the rotation of the 187d residue ,
are also more endothermic in 
 sc18d222 g than in nl602d222 g . as a result ,
the equilibrium 
 occupation probability of mode 2 in sc18d222 g is rather low , 
 which is consistent with the previous observation that mode 2 is rarely 
 observed during md simulations . for nl602d222 g
, 
 the paths are similar , but the energy difference is lower . hence , 
 binding mode 2
can still be observed in the nl602d222 g mutant 
 ha , albeit with decreased probability .
a more detailed analysis of 
 this key elementary step will be addressed in future work one reason 
 why the binding of the 3-hydroxyl of gal to the side chain of 222d 
 is more exothermic and has a lower barrier in nl602 may be the differences 
 in the shape of the binding pocket .
when the ligand binds to 222d / g 
 or unbinds from 187d in sc18 , we always observe a significant shift 
 of the gal and glcnac sugars toward the 222 position .
we examined the shape of the 
 pocket for the four strains using the povme program and illustrate the differences observed in figure 5 .
structures were aligned based on the 190-helix , 130-loop , 
 and 220-loop using vmd . two major differences 
 between sc18-wt and nl602-wt can be observed in figure 5a : ( 1 ) the pocket in sc18-wt is larger than that in nl602-wt 
 on the side of 220-loop and 190-helix , where the majority of the interactions 
 between ligand and pocket are . ( 2 ) the pocket in nl602-wt is larger 
 on the side of the 130-loop and also between the 220-loop and 130-loop . 
 comparing d222 g mutants in figure 5b , the blue 
 region is larger than that in figure 5a .
however , 
 the extra volume is located mainly at the top of the pocket , which 
 may arise due to measurement error , since the pockets are open at 
 the top .
pocket shape differences between ( a ) sc18-wt ( cyan ) and nl602-wt 
 ( pink ) and ( b ) sc18d222 g ( cyan ) and nl602d222 g ( pink ) . 

the blue and red regions correspond to areas where the pocket in sc18 
 is larger and than that in nl602 , respectively .
to analyze the pocket shape differences further , in figure 6 we show a triangular region defined by 187d , 219k , 
 and 222d /
the shape of this region provides a simple representation 
 of the 
 220-loop side of the binding pocket , where the majority of the rearrangements 
 between the two binding modes are localized .
we find that the 219k 
 residue in sc18 is about 12  further away from 
 187d but slightly closer to 222d / g compared to 
 nl602 . as a result ,
224a 
 interacts with the 4-hydroxyl of gal in binding mode 1 and 224e in 
 nl602 interacts with the 3- and 4-hydroxyl of gal . hence , 224a
/ e mainly 
 determines how far to the left ( from the viewpoint of figure 6 ) and how deep 2,6-sa binds within the pocket . 
 in sc18 ,
the distance between 224a and the line connecting 187d and 222d / g 
 is 2.8  in the wild type ha and 2.5  in the 222 g mutant 
 ha . in nl602 ,
the position of 224e is close to the line connecting 
 187d and 222d / g .
the distance between 224e and the line connecting 
 187d and 222d / g is 1.1  in the wild type ha and 1.4  in 
 the mutant ha , and the distance between 224a and 187d in sc18 is also 
 longer than that in nl602 . since 222d
protrudes slightly from the 
 edge of the pocket , the position of 224a in sc18 allows the ligand 
 to bind deeper and closer to the left side .
the variations in the 
 distance between the 220-loop and 190-helix confirm the shape differences 
 observed from the povme measurement : the pockets in sc18-wt and sc18d222 g 
 are larger than for nl602-wt and nl602d222 g on the side of 
 the 220-loop .
( a  d ) triangle formed by 187d , 219k , and 222d / g 
 in mode 
 1 of ( a ) sc18-wt , ( b ) sc18d222 g , ( c ) nl602-wt , and ( d ) nl602d222 g , 
 respectively . the distances between residues are calculated between 
 the backbone nitrogen atoms , and the values in brackets are the corresponding 
 distances in mode 2 .
binding pocket changes during the interconversion between the two 
 modes in four different strains .
the blue and red regions correspond 
 to volumes lost and gained after a particular displacement , respectively . 

( c ) rotation 
 of the 8-hydroxyl of gal and g4 and g3 shifts 
 in sc18d222 g ( m1m9 in figure s3 ) .
( d ) r187 movement in sc18d222 g ( m9m11 
 in figure s3 ) .
( e ) rotation of the 9-hydroxyl 
 of gal and g4 shifts in nl602-wt ( m4m6 in figure s5 ) .
( h ) g4 , g3 , 
 180h , and 191l movement in nl602d222 g ( m7m11 in figure s7 ) .
we also used povme to identify changes in pocket shape during 
 the 
 interconversion .
for sc18-wt , the first nontrivial shape change of 
 the pocket results from the rotation of the 4-hydroxyl of gal ( figure 7a ) , which shrinks the volume between the 220-loop 
 and the 130-loop .
this reduction in size is caused by both the conformational 
 change of the 223q side chain and the movement of the 220-loop toward 
 the 130-loop .
another significant shape change results from the rotation 
 of 187d ( figure 7b ) , which reduces the size 
 of the region in front of 187d .
the observations in the case of sc18d222 g 
 are similar to those for sc18-wt .
the first few movements , including 
 the rotation of the 8-hydroxyl of gal and the movement of the 4- and 
 3-hydroxyl of gal , contribute to the shrinking of the pocket between 
 the 220- and 130-loops , as shown in figure 7c .
the side chain rotation of 187d then removes the region in front 
 of 187d ( figure 7d ) .
the changes in pocket 
 shape for nl602-wt differ greatly from those 
 for sc18-wt .
the rotation of the 9-hydroxyl of gal shrinks the pocket 
 between the 220-loop and 190-helix and enlarges the volume between 
 the 190-helix and 130-loop , as shown in figure 7e .
the rotation of 187d , however , contributes less to the shape change 
 than it does for sc18-wt , which is consistent with the previous observation 
 that the rotation of 187d does not remove many interactions with the 
 ligand .
the rotation of 187d slightly shrinks 
 the region in front of 187d , and enlargement of the pocket next to 
 224e occurs via movement of the 224e side chain ( figure 7 g ) .
the next few steps , including the movement of the 4- and 
 3-hydroxyl and residues 180 h and 191l , slightly shrink the pocket 
 between the 220-loop and 190-helix , as for nl602-wt . in summary , 

the two strains from 1918 south carolina are more flexible 
 between the 220-loop and 130-loop but more rigid between the 220-loop 
 and 190-helix when compared to the 2009 netherlands strains .
this 
 observation is consistent with the distance changes among key residues 
 ( figure 6 ) .
the maximum distance variation 
 in sc18-wt is only 0.3  for 187d and 222d .
however , in nl602 , 
 the distance from 187d to residues 219k , 222d , and 224a is reduced 
 by roughly 1  in binding mode 2 .
since the ligand mainly interacts 
 with the 190-helix and the 220-loop , the flexibility of this region 
 plays an important role in the two binding modes .
when the transformation 
 between the two modes occurs , due to the shape and the size of the 
 pocket in nl602 , the 220-loop can adapt its position to reflect the 
 position of the ligand , which allows most interactions with the edge 
 of the pocket to be conserved .
however , for sc18 , because the pocket 
 is larger , the ligand must move from one side to the other during 
 the transformation to maintain these interactions , and the 220-loop 
 can not adapt to the ligand position as flexibly as for nl602 . as a 
 result , the binding of 2,6-sa to sc18 relies more heavily on 
 the side chain of 222d than in nl602 , since the ligand has limited 
 access to the edge of the pocket .
the additional interactions between 
 the ligand and the pocket in nl602 that were reported previously can 
 also been explained in the light of this observation .
key residues in ( a ) sc18 and ( b ) nl602 that determine the relative 
 positions of the 190-helix and 220-loop .
interestingly , as shown in figure 8a , the junction structures 
 between the 190-helix 
 and 220-loop , which determine the relative position of the two regimes , 
 are different in the two has . in sc18
183p is in front of the 190-helix 
 and does not interact with other residues nearby .
the hydrogen - bonds 
 between 215a and 182p , together with 217r and 224a , determine the 
 relative positions of the 190-helix and the 220-loop .
the residues 
 at the junction between the 190-helix and 220-loop change from 183p 
 and 224a in sc18 , to 183s and 224e in nl602 . due to these differences , 

the shape of the pocket in nl602 is affected by the interaction between 
 183s and 213e , as well as the interaction between 183s , 224e and 226r , 
 as shown in figure 8b .
there have been 
 several experimental results reported since 2010 
 that the s183p mutation and the e224a mutation with or without other mutations , 
 are correlated with the pathogenesis and transmission changes observed 
 for h1n1 . since the majority of the shape difference between sc18 
 and nl602 during the binding mode interconversion occurs between the 
 190-helix and 220-loop , the accumulated mutations at positions 183 
 and 224 may explain the differences in pocket shape flexibility during 
 binding mode interconversion observed here .
in the present work , we have investigated the 
 transformation mechanism 
 between two binding modes for 2,6-sa in both the wild type 
 and mutant h1n1 hemagglutinin proteins from 1918 and 2009 pandemic 
 flu using the 6sln ligand .
we introduced an intermediate network 
 scheme to accelerate the search process of pathways connecting the 
 two binding modes .
the advantage of this procedure is that searches 
 can be systematically conducted in parallel for all the single changes 
 in identifiable local conformational states . by comparing the paths 
 predicted to dominate the kinetics , we have identified the key elementary 
 step that distinguishes the two binding modes , namely the 3-hydroxyl 
 of gal either binding to the side chain of 222d in wild type has or 
 moving to the place where it would have bound to the side chain of 
 222d in mutant has with 222 g .
we find that the different residues 
 at the 183 and 224 positions in the two has change the connecting 
 pair between the 190-helix and 220-loop , and hence the shape of the 
 pocket .
the smaller size of the binding pocket in nl602 allows the 
 ligand to interact with more residues in nl602 than in sc18 and also 
 makes the 220-loop more flexible .
the d222 g mutation exchanges 
 the energy preference between the 
 two binding modes in both sc18 and nl602 hemagglutinins .
however , 
 the transition between mode 1 and mode 2 raises the energy more in 
 sc18 than in nl602 , probably because of the different shape of the 
 binding pocket .
the different shape , the flexibility of the 
 pocket during binding 
 mode interconversion , and the resulting energy profile changes provide 
 insight into the contrasting pathogenicities induced by the d222 g 
 mutation in the sc18 and nl602 viruses .
these results may help to 
 explain the difference between the impact of the 1918 and 2009 h1n1 
 strains .
a change in preference from 2,6- to 2,3-linked 
 sialic acid is associated with a shift in pathogenicity .
however , 
 in the 2009 d222 g mutant , the 2,6 affinity is not weakened as much 
 as for the corresponding 1918 mutant .
understanding such distinctions 
 at an atomic level of detail should be helpful in predicting the potential 
 emergence of strains that may pose serious threats to human health .
 in 
 2009 , the d222 g mutation in the hemagglutinin ( ha ) glycoprotein 
 of pandemic h1n1 influenza a virus was found to correlate with fatal 
 and severe human infections . 
 previous static structural analysis suggested 
 that , unlike the h1n1 viruses prevalent in 1918 , the mutation did 
 not compromise binding to human 2,6-linked glycan receptors , 
 allowing it to transmit efficiently . 
 here we investigate the interconversion 
 mechanism between two predicted binding modes in both 2009 and 1918 
 has , introducing a highly parallel intermediate network search scheme 
 to construct kinetically relevant pathways efficiently . 
 accumulated 
 mutations at positions 183 and 224 that alter the size of the binding 
 pocket are identified with the fitness of the 2009 pandemic virus 
 carrying the d222 g mutation . 
 this result suggests that the pandemic 
 h1n1 viruses could gain binding affinity to the 2,3-linked 
 glycan receptors in the lungs , usually associated with highly pathogenic 
 avian influenza , without compromising viability . 