ion channels are the canonical example 
 of a protein nanomachine ; 
 they perform vital physiological functions , are exquisitely sensitive , 
 and can make solid - state electronic devices look crude by comparison .
they have been studied extensively for over 50 
 years , and a wealth of electrophysiological , structural , and biophysical 
 data have been accumulated .
we can not 
 yet calculate the conduction properties of an ion channel from an 
 experimental ( static ) structure .
instead , once a structure of a particular 
 ion channel becomes available , its conduction properties are explained 
 intuitively , perhaps by comparison to other channels .
calculating functional properties of ion channels from their 
 structures , directly linking structure and function , not only is desirable 
 in itself but also is an excellent test of our ability to accurately 
 model the behavior of membrane proteins in general .
we shall 
 focus on potassium ion channels , as the first structure 
 of an ion channel belongs to this family and it remains the family 
 with the most experimental structures .
all potassium ion channels 
 have a tetrameric central pore domain through which potassium ions 
 and waters pass . for example compare kcsa , a bacterial ph - gated k channel , and kvchim , a chimeric mammalian voltage - gated k channel ( figure 1 ) .
their central pore domains are remarkably similar despite 
 the two channels being gated by different stimuli and both having 
 additional large domains , namely , voltage sensors and an intracellular 
 tetramerization domain for kvchim versus a c - terminal intracellular 
 domain for kcsa .
the primary difference between the high - resolution 
 structures of the pore domains of kcsa and kvchim is that the intracellular 
 gate is closed in the former but open in the latter .
potassium ions can flow through in either direction ; 
 for convenience , 
 we shall consider a potassium ion approaching either channel from 
 the cytoplasm .
it first moves past the intracellular gate , if open , 
 and enters the central cavity .
the ion then passes through the narrowest 
 part of the pore , the selectivity filter , before exiting the cell . 

the selectivity filter is formed by four loops , one from each monomer , 
 each having the conserved sequence tvgyg .
the loops are arranged such 
 that their backbone carbonyls point toward the pore axis creating 
 a series of octahedral cages of oxygens in which potassium ions and 
 waters are stabilized ; these are conventionally labeled s0s4 
 ( figure 1 ) . before passing through the selectivity 
 filter
the dominant 
 view is that waters and ions are thought to enter the filter alternately 
 and move through in single file , with at least two ions bound at all 
 times .
we shall call this the ion  water  ion 
 ( kwk ) mechanism ( figure 1 ) . as a potassium 
 ion approaches the filter from the central cavity , a 
event occurs leading to a potassium ion being ejected from the top 
 of filter . as the name suggests , the selectivity filter is thought 
 to discriminate between potassium and sodium ions .
potassium ion channels ( a ) kcsa ( in 
 red ) , a ph - gated bacterial channel , and ( b ) kvchim ( in blue ) , a voltage - gated chimeric mammalian channel , 
 share a similar central pore domain .
ions pass through the channel along the 4-fold symmetric 
 pore axis , z , shown as a dotted gray line .
the c - terminal 
 domain of kcsa and tetramerization domain of kvchim are not shown . 

the narrowest part of the pore , the selectivity filter , is found above 
 the central cavity toward the extracellular side of the bilayer ; this 
 region is enclosed in a gray box .
the selectivity filters of ( c ) kcsa 
 and ( d ) kvchim have the same tvgyg amino acid sequence and are structurally 
 very similar .
the backbone oxygens , along with the side chain hydroxyl 
 of the threonine residue , form octahedral cages of oxygen within which 
 an ion or a water are stabilized .
we follow the nomenclature of bernche 
 and roux and label the three potassium 
 ions as shown .
the 
 position of the center of mass of the top two ions along the pore 
 axis , z , is therefore denoted by z12 , while the position of the lower ion is denoted by z3 ( see methods ) . we 
 shall describe the configuration of the selectivity filter as a string 
 of letters denoting the species present in each binding site , from 
 s0 to s4 , reading from left to right with the cavity species , where 
 appropriate , given in parentheses . the configuration shown is kwkwk . 

( e ) kcsa has a glutamate at position a. this is protonated , leading 
 to a hydrogen bond between gluh71 and asp80 ( drawn in green ) . in kvchim ,
the equivalent residues are val366 
 and asp375 , and therefore , no such hydrogen bond is possible . for 
 clarity ,
only two of the four monomers are drawn , and the coordinating 
 oxygens are red .
although the selectivity 
 filters of kcsa and kvchim have the same 
 tvgyg sequence ( figure 1c and d ) , kcsa has 
 a very low open probability due to c - type inactivation . replacing 
 the glutamate at position 71 of kcsa with alanine
kv channels typically have a valine 
 or isoleucine at the equivalent position to glu71 and so lack this 
 hydrogen bond . the e71a kcsa mutant therefore 
 behaves more like a kv channel , such as kvchim .
over half a 
 billion potassium ions pass through a single kcsa channel 
 each second at physiological voltages .
this is near the diffusive limit for potassium ions , which is remarkable 
 considering the high selectivity this channel has for potassium over 
 sodium ions . according to nernst 
planck theory estimates , this implies that the potassium ions must not 
 encounter any barriers greater than 5 kt ( 3 kcal / mol at 300 
 k ) , otherwise the conductance would be markedly reduced .
since 
 the first structure of kcsa was published , theoretical methods have been applied to investigate or 
 calculate the conductance of a range of potassium ion channels .
we shall consider here hierarchical approaches 
 that use a potential of mean force ( pmf ) .
the first step in the calculation of a pmf is to identify which 
 species to follow , and their corresponding reaction coordinates , zi , are understood as a set 
 of relevant collective variables that emerge naturally from the underlying 
 transport dynamics . in ion channels ,
this is straightforward ; one 
 wants to follow the motion of the n ions along the 
 pore axis , z. each point on the n - dimensional pmf landscape , ( z1 , z2 , .. zn ) , defines a unique configuration of 
 ions within the selectivity 
 filter , and the gradient yields the average force experienced by each 
 of the conducting species .
it has been shown that the total pmf , total(z1 , z2 , 
 .. zn ) , can be written 
 as the following:1where equilibrium(z1 , z2 , .. zn ) 
 is the potential of mean force measured 
 at equilibrium , i.e. , with a transmembrane potential of zero , q is the elementary charge and membrane(zi ) describes how the 
 transmembrane potential changes along the pore axis .
the effect of 
 the channel on the permeating species is collapsed down to a single , 
 multi - dimensional pmf , equilibrium(z1 , z2 , 
 .. zn ) , which , along with 
 the variation in 
 transmembrane potential along the pore axis , membrane(zi ) , may be used as 
 inputs into a higher - level model that predicts 
 the conductance of the channel .
the first equilibrium pmf was 
 calculated over 10 years ago by bernche 
 and roux using the first low - resolution 
 structure of kcsa .
this seminal study 
 firmly established the technique , and the authors suggested a plausible 
 conduction mechanism . because the barriers were 23 kcal / mol 
 high ( 5 kt ) they were able to use this pmf to calculate , with 
 some success , the current  voltage relationship of kcsa using 
 brownian dynamics .
subsequent studies 
 by other groups calculated pmfs for different potassium channels suggested there might be conduction mechanisms involving vacancies , examined the selectivity , and even applied the methodology to bacterial outer membrane 
 pores .
following an early pioneering 
 study and the relentless increase in 
 computational power , it is now possible to directly simulate multiple conduction events at different magnitudes of the 
 electric field thereby permitting the direct in silico measurement 
 of the conductance .
a direct approach is intuitive but inevitably requires 
 far more computational resource than the equivalent pmf - based method , 
 and there remain difficulties applying an electric field in a simulation . 

calculating an appropriate pmf also allows other important questions 
 to be answered ; for example , one can study the stability of multi - ion 
 configurations and devise kinetic models to compute the escape rate 
 of an ion over an energetic barrier .
in this paper , we shall 
 calculate seven independent 2d pmfs that 
 encapsulate the conduction of potassium through two different channels , 
 kcsa and kvchim .
we shall ( 1 ) introduce the e71a mutation into kcsa 
 to facilitate comparison with kvchim , ( 2 ) determine the effect on 
 both channels of removing the charmm cmap forcefield correction , and 
 ( 3 ) examine the reproducibility of these types of calculation .
this 
 will permit us to discuss the current feasibility of calculating the 
 conductance of potassium ion channels from experimental structures .
we shall first focus on examining the sensitivity and reproducibility 
 of equilibrium pmfs , equilibrium(z1 , z2 , 
 .. zn ) , before later examining 
 how the transmembrane 
 potential varies along the pore axis , membrane(zi ) . following the method of 
 bernche and roux , we used a two - coordinate 
 system ( z12 , z3 ) to follow the motion of potassium ions 1 , 2 , and 3 in the selectivity 
 filter ( figure 1 and methods ) . like the majority of studies
, we assumed that the sites of the 
 selectivity filter are alternately filled by water molecules and potassium 
 ions ( kwk mechanism ) .
we produced a series of carefully equilibrated 
 structures of the pore domain of wild - type kcsa embedded in lipid bilayers as described in methods .
each structure had a different configuration of ions 
 and waters in the selectivity filter , all of which conformed to the 
 kwk mechanism .
these were then used to seed 577 simulations , each 
 of which was 0.2 ns long and sampled a different region of ( z12 , z3 ) space by 
 applying biasing forces to the potassium ions .
close inspection 
 of individual umbrella simulations revealed two complications concerning 
 the assumption of a simple kwk mechanism .
first , one or more val76 
 carbonyls on occasion flipped by 180 , reducing the oxygen 
 coordination number within sites s2 and s3 ( figure 2a ) and breaking the 4-fold symmetry of the selectivity filter
. 
 the flipping of carbonyls , especially those of the valine , has been 
 observed several times before and has been suggested to lead 
 to c - type inactivation .
second , 
 the kwk ordering was not always maintained . in common with previous 
 studies
, we have implicitly assumed that applying biasing forces to the potassium 
 ions would trap the water between the cations , as long as the cations 
 remain within the confines of the filter .
this has the advantage of 
 reducing the number of collective variables and , therefore , the dimensionality 
 of the pmf from five to two .
however , our simulations show that this 
 assumption does not always hold ( figure 2b ) . 

first , a water trapped between two potassium ions in the filter can 
 occasionally escape by squeezing past one of the ions leaving a vacancy 
 ( e.g. , kwwkk configuration in figure 2b ) .
second , 
 as ion 3 approaches the bottom of the selectivity filter , more than 
 one water may become trapped between it and ion 2 inside the filter . 
 reversing
this process leads to a similar behavior , i.e. , moving ion 
 3 from inside the filter into the central cavity because the trapped 
 water above ion 3 can escape as before .
we can not therefore assume 
 that applying biasing forces to only the potassium ions will strictly 
 maintain the alternating sequence of ions and waters necessary for 
 the kwk conduction mechanism .
this is an important result , and we 
 shall discuss it later . checking the assumptions that ( 1 ) the selectivity 
 filter remains 
 conductive and ( 2 ) that the ions and water proceed through the selectivity 
 filter in an alternating single file .
( a ) log histogram of the backbone 
  dihedral angle of val76 shows that the first assumption does 
 not always hold ; 2.7% of valines have a  > 0 that 
 corresponds 
 to a flipping of the backbone carbonyl oxygen and , as suggested by 
 bernche , leads to a selectivity 
 filter with altered conduction properties . because kcsa is a tetramer 
 and only a single flip
we remove these frames and those where the kwk mechanism 
 is not followed from all subsequent analyses .
( b ) to check the second 
 assumption , we have plotted the number of waters between adjacent 
 ions as a function of the z coordinate of the ion 
 furthest from the center of the selectivity filter .
the area of each 
 blue square is proportional to its percentage occurrence at that value 
 of z. we further define an ion to be within the selectivity 
 filter if it is within sites s1s4 ( gray area ) ; for a detailed 
 definition , see methods .
there is only a single 
 water between any two ions in the selectivity filter for most ( 98% ) , 
 but not all , of the time ( circled in red ) .
waters occasionally squeeze 
 past the potassium ions in the filter leading to unexpected configurations , 
 two examples of which are shown . to resolve these two problems
, we only kept data where there 
 was 
 exactly one water between any two adjacent ions ( as long as both ions 
 are  inside  the filter , see methods ) , and there were no carbonyl flips .
as a result , 12% of the data 
 were discarded ( table s1 , supporting information ) .
the resulting pmf ( and all others in this paper ) hence assumes 
 ( 1 ) the kwk conduction mechanism and ( 2 ) the filter is in a conducting 
 state .
we are therefore implicitly testing the classical kwk conduction 
 mechanism , because if it is true , we would expect the heights of the 
 barriers to be < 5 kt . following this filtering , each of the 577 umbrella 
 simulations was divided into 100 bins , each 2 ps wide .
a 2d pmf was 
 calculated for each using the weighted histogram analysis method and the minimum free energy path ( mfep ) found 
 as described in methods .
the height of the 
 first barrier was measured and plotted as a function of time .
calculating 
 the statistical inefficiency suggested 
 that we should discard the first 100 ps of the data to ensure equilibration 
 and that a conservative estimate of the correlation time for the remaining 
 100 ps is 20 ps ( figure s1 , supporting information ) .
because each umbrella simulation is 200 ps long , five independent 
 2d pmfs were calculated and averaged , with errors estimated in the 
 usual way ( figure s2 , supporting information ) .
the resulting averaged 2d pmf for wild - type kcsa is shown in figure 3a .
the 2d pmf has three minima corresponding 
 to the following configurations : wwkwk(k ) , wkwkw(k ) , and kwkwk ( see 
 the legend of figure 1 for the definition of 
 this nomenclature ) .
the mfep between the wwkwk(k ) and kwkwk configurations 
 was calculated using the nudged elastic band method ( see methods ) and , as expected , it visits the wkwkw(k ) 
 configuration .
it appeared possible , by inspection , 
 that there are alternative 
 pathways between the wkwkw(k ) ( point 3 ) and kwkwk ( point 5 ) configurations . 
 calculating the mfep using a slightly different initial condition 
 bore out our suspicion and shows that there are two main pathways , 
 labeled path 1 and path 2 .
the latter is similar , but not identical , 
 to the secondary pathway previously seen by bernche and roux .
the height of the barrier encountered is similar 
 along both paths ( figure 3b , table 2 ) , and hence , they are degenerate . consider for 
 a moment the shape of the mfep if the conducting species were hard 
 spheres and the selectivity filter a cylinder  because the ions 
 and waters would smoothly move through the filter , the mfep would 
 be a straight line between the two bound configurations .
both paths 
 are instead curved in an approximately symmetrical fashion about the 
 idealized mfep . from path 1
, we infer that the ions and waters are 
  compressible  because ion 3 has to nearly enter s4 
 before the knock - on occurs , and from path 2 , we infer that , as expected , 
 the cations can interact at a distance .
there is therefore an equal 
 probability of the two ions in the filter  knocking - on  
 either spontaneously once ion 3 is immediately below s4 ( path 2 ) , 
 or after ion 3 has begun to enter s4 ( path 1 ) . following the 
 mfep allows us to suggest how potassium moves through the selectivity 
 filter of kcsa ( figure 3d ) . as ion 3 approaches 
 the filter from the central cavity , ion 1 moves from site s2 to s1 , 
 ion 2 moves from site s4 to s3 , and a water simultaneously enters 
 the filter to occupy site s4 .
as ion 3 approaches thr70 , a second  knock - on  
 event occurs .
this either occurs spontaneously with ion 3 still some 
 distance from entering s4 ( path 1 ) , thereby creating a 
partial 
 vacancy  , or it happens when ion 3 has begun to enter s4 .
after 
 the second  knock - on  , during which all the species 
 within the filter move one site closer to the extracellular side of 
 the channel , ion 3 enters s4 .
the net effect of these motions is the 
 displacement of a potassium ion from the central cavity , where it 
 is entirely coordinated by water , to site s0 , where it is only partly 
 coordinated by water .
this pmf therefore captures almost an entire 
 single conduction event and suggests a plausible conduction mechanism 
 that is consistent with the original proposals of hodgkin and keynes .
the quasi-1d pmf suggests that the kwkwk configuration , 
 which is the only one with three ions inside the filter , is more stable 
 than the other two by 610 kcal / mol , although this could be 
 because the intracellular gate is closed , destabilizing the wwkwk(k ) 
 configuration .
the heights of the kinetic barriers 
 between the configurations are 511 kcal / mol ( table 2 ) , which is significantly larger than the estimate 
 of 3 kcal / mol ( 5kt ) required by the electro - diffusion 
 model of bernche to reproduce conduction at the rates observed 
 for this channel .
there are three stable configurations of the 
 selectivity filter 
 of kcsa and three potassium ions .
( a ) the 2d pmf of wild - type kcsa 
 with the cmap correction has three minima .
the axes correspond to 
 the center of mass of the two potassium ions closest to the periplasm 
 ( ion 1 and ion 2 , z12 ) and the position 
 of the third potassium ion ( z3 ) .
contours 
 are drawn every 1 kcal / mol , and the minimum free energy path ( mfep ) 
 as calculated by the nudged elastic band method is drawn as a black 
 line .
a second , degenerate path between points 3 and 5 is drawn as 
 a gray line .
a dashed straight line between the wkwkw(k ) and kwkwk 
 configurations is drawn ; this is the mfep in the idealized case where 
 the permeating species are represented as hard spheres .
the smallest 
 value of the free energy is defined as zero by the wham procedure . 

( b ) the free energy along the mfep more clearly shows the 
 heights of the barriers between the identified stable states .
( c ) the relative 
 free energies of the three stable configurations and the heights of 
 the forward and reverse barriers between them is defined as shown 
 ( table 2 ) .
( d ) six 
 illustrative images of the positions of the ions and waters and the 
 structure of the selectivity filter .
these correspond to the points 
 labeled in a and b. ion 1 and ion 2 are blue , whereas ion 3 is cyan . 
 for clarity , only two of the four monomers are drawn , and the coordinating 
 oxygens are red
. we shall now examine how 
 sensitive this 2d pmf is to a number of factors .
an additional four 
 sets of umbrella simulations were run , kcsa-10 , kcsa-00 , kvchim-10 , 
 and kvchim-00 ( table 1 ) .
the data from these 
 1194  0.2 ns simulations were filtered as before ( table s1 , 
 figures s3 and s4 , supporting information ) to remove instances where either one or more carbonyl oxygens had 
 flipped or the ions and waters did not stay in an alternate single 
 file .
we assumed that , as shown for kcsa-11 , discarding the first 
 100 ps of each umbrella simulation would ensure equilibration , and 
 20 ps is a reasonable estimate of the correlation time for the remaining 
 100 ps .
this permitted us to calculate five independent 2d pmfs from 
 each set of umbrella simulations ( figure s5 , supporting 
 information ) and estimate errors in the usual way ( figure s6 , supporting information ) .
in most cases , the error 
 is 1.0 kcal / mol along the mfep , although it reaches 1.5 kcal / mol 
 in places .
the additional four filtered converged 2d pmfs , along with 
 that of wild - type kcsa ( kcsa-11 ) , are plotted in figure 4 , and illustrative snapshots of the configurations of the 
 selectivity filters can be found in figure s7 of the supporting information .
the kcsa-11 and kvchim-10 umbrella 
 simulations can be further decomposed into two independent sets of 
 simulations as shown , and hence , two 2d pmfs can be calculated from 
 each ( figure 5 ) .
the h - bond column details 
 whether a hydrogen bond can be formed at posa ( figure 1e ) in the channel .
the two digits in the name of the simulation 
 indicate whether the cmap correction was used ( 1 ) or not ( 0 ) and whether 
 the h - bond is present ( 1 ) or not ( 0 ) , respectively .
how does the 2d pmf change 
 when the kvchim channel is used in place 
 of kcsa ? because the selectivity filter of kvchim is similar to the 
 e71a mutant of kcsa , let us first consider the effect of the e71a 
 mutation on kcsa ( kcsa-11 vs kcsa-10 ) .
almost a third of the data 
 from the kcsa-10 kcsa e71a umbrella simulations were removed , primarily 
 due to carbonyl flips of val76 ( table s1 , supporting 
 information ) .
the 2d pmf does not , as a result , sample the 
 wwkwk(k ) configuration ; the start of the mfep was instead placed close 
 to the wkwkw(k ) configuration .
both this and the kwkwk configurations 
 are identified , and the heights of the barriers are also similar to 
 the wild - type ( table 2 ) , although the difference in free energy between the two configurations 
 is slightly larger for the mutant .
again , two paths between the wkwkw(k ) 
 and kwkwk configurations can be identified , although they are more 
 similar to one another , and so do not represent significantly different 
 conduction mechanisms .
we can now compare the 2d pmfs of e71a 
 kcsa and kvchim ( figure 4b and d ) .
although 
 the selectivity filters of kcsa 
 and kvchim have the same sequence , one should not expect identical 
 2d pmfs because , as has been suggested before , the energetics of conduction 
 may be altered by the remainder of the channel . even if this effect is small
, one would still expect the energetics 
 of the wwkwk(k ) configurations to be different because the intracellular 
 gate that is directly below the central cavity is closed in kcsa but 
 open in kvchim .
the kvchim 2d pmf ( kvchim-10 ) 
 clearly identifies the kwkwk and wkwkw(k ) configurations , and again , 
 there are degenerate mfeps that connect the two
unfortunately , we can not determine if 
 this state exists for kcsa because , as mentioned above , most of the 
 data from this region of the 2d pmf were removed by the filtering . 
 in the kvchim pmf , the wkwkw(k )
configuration has a slight shoulder 
 at lower values of z12 , and the shape 
 of the well defining the kwkwk configuration is more elongated in 
 the z3 direction , indicating that the 
 free energy is less sensitive to the position of ion 3 than in e71a 
 kcsa , perhaps because the intracellular gate is closed in the latter 
 structure .
water selectivity filter configurations 
 are insensitive to the use of cmap or the presence of the gluh71asp80 
 hydrogen bond .
the relative stability of these states and the heights of 
 the barriers between the states , however , changes markedly .
( a ) the 
 result from figure 3 is shown in the top left ; 
 this is for wild - type kcsa ( i.e. , possessing the hydrogen bond behind 
 the selectivity filter ) with the cmap correction to the charmm forcefield . 

( b ) removing this hydrogen bond using the e71a mutant appears to increase 
 the relative stability of the kwkwk configuration . ( c )
when the cmap 
 correction is not used , the kwkwk and wkwkw(k ) configurations have 
 the same energy .
( d ) the 2d pmf for the pore domain of the paddle 
 chimera , which naturally does not have the hydrogen bond behind the 
 filter , is similar to the equivalent pmf of kcsa ( b ) .
( e ) not using 
 cmap , however , leads to quite a different 2d pmf . in all cases , the 
 same color scheme and labeling
is used as in figure 3 , and contours are drawn every 1 kcal / mol .
regions colored 
 white are where data were removed due to the kwk mechanism not being 
 followed in the simulations .
these gaps occasionally prevented the 
 full mfep being calculated , as in b and e. where possible , an alternate 
 mfep between points 3 and 5 ( path 2 ) is drawn in gray .
for comparison 
 5 kt = 3 kcal / mol 
 at 300 k. for a definition of the barrier heights , e.g. , g13 , see figure 3c . like the majority of previous studies , 
 we have used the charmm 
 forcefield because its potassium  carbonyl oxygen nonbonded 
 parameters have been altered to better reproduce the interaction energies 
 between potassium and both water and n - methylacetamide .
the charmm forcefield has undergone further 
 improvements in recent years , the most notable being the cmap correction .
this enhances the protein ( , ) 
 angular dependence of the potential energy and has been shown to improve 
 the simulated dynamics of soluble proteins .
the study by bernche and roux predated the development of the cmap correction , but it has subsequently 
 been shown to be required to simulate the cation selectivity of the 
 kv1.2 channel .
let us therefore examine 
 the effect of not using this forcefield modification for both the 
 e71a kcsa mutant ( figure 4b , c ) and wild - type 
 kvchim ( figure 4d , e ) . in both cases , we identify 
 all three expected stable configurations of the selectivity filter , 
 but not using cmap alters the heights of the kinetic barriers between 
 the configurations .
the wkwkw(k ) and kwkwk configurations of 
 the cmap e71a kcsa mutant 
 ( kcsa-00 ) have similar energies and a similar level of degeneracy 
 in the mfep compared to when the cmap correction is used ( kcsa-10 ) , 
 with the barrier between these configurations being reduced to 6 
 kcal / mol ( table 2 ) .
the wwkwk(k ) configuration 
 has a free energy of 9 kcal / mol more than the kwkwk configuration , 
 similar to what was seen in the 2d pmf of wild - type kcsa with the 
 cmap correction ( kcsa-11 ) .
all three stable configurations of 
 the filter are picked out by 
 the 2d pmf of kvchim when the cmap correction is not used ( kvchim-00 ) , 
 although the wwkwk(k ) configuration is only marginally stable .
part 
 of the 2d pmf between the wkwkw(k ) and kwkwk configurations is not 
 defined , making it difficult to determine the mfep and accurately 
 determine the height of the barrier .
this is despite that the filtering 
 only removing 7.5% of the data ( table s1 , supporting 
 information ) .
compared to kvchim-10 , the kwkwk configuration 
 is more stable than the other configurations .
this is the opposite 
 to what we observe when not using the cmap correction factor for e71a 
 kcsa .
( a ) the kcsa-11 2d pmf in figure 3 can be 
 decomposed into two independent 2d pmfs ; these are shown here .
note 
 that kcsa-11a and kcsa-11b were seeded from different structures , 
 and different amounts of data were discarded in ensure a kwk conduction 
 mechanism ( table s1 , figures s3 and s4 , supporting 
 information . )
( b ) the kvchim-10 2d pmf in figure 4 can also be decomposed into two independent 2d pmfs .
both 
 pmfs identify the same bound configurations of the selectivity filter 
 but differ significantly in their estimate of the free energy change 
 between the states .
kvchim-10a and kvchim-10b were seeded using the 
 same equilibrated structure of the pore of the paddle chimera , and 
 similar proportions of data were removed to ensure that the kwk mechanism 
 is maintained .
two of the five 2d 
 pmfs ( kcsa-11 and kvchim-10 ) were calculated by combining two independent 
 data sets of umbrella simulations .
we will now look at these data 
 sets separately ( kcsa-11a and kcsa-11b and kvchim-10a and kvchim-10b , 
 table 1 ) to test the reproducibility of the 
 2d pmf maps .
the kcsa-11a data set was produced using structures of 
 kcsa embedded in a 7:3 pope / popg lipid bilayer , whereas a simple popc 
 lipid bilayer was used in the simulation used to seed the kcsa-11b 
 data set and all other simulations .
although anionic lipids have been 
 shown to modulate the opening probability of kcsa , they do not change the conductance , and therefore , we do 
 not expect to observe any differences between kcsa-11a and kcsa-11b . 

the umbrella spring constant k12 in the 
 kcsa-11b set was twice as large as the kcsa-11a set , which in turn 
 necessitated tighter spacing between umbrella simulations along the z12 coordinate , and hence , 381 umbrella simulations 
 contribute to the kcsa-11b pmf , whereas only 196 were required for 
 the kcsa-11a pmf .
only 6% of frames were removed from the kcsa-11a 
 simulations to ensure an kwk conduction mechanism with no carbonyl 
 flips , whereas 15% were removed from the kcsa-11b data set ( table 
 s1 , supporting information ) . despite these 
 differences ,
the kcsa-11a and kcsa-11b pmfs and their minimum free 
 energy paths , including the degenerate path 2 , are remarkably similar 
 ( figure 5a ) .
the only difference was the spacing of the 
 umbrella simulations , and a few additional umbrella simulations with 
 a higher value of k12 were required for 
 kvchim-10b .
both pmfs identify the wkwkw(k ) and kwkwk configurations 
 and suggest that the kwkwk configuration is the most stable ( figure 5b ) .
unfortunately , removing the data necessary to 
 ensure the kwk mechanism is maintained created an unsampled region 
 in the kvchim-10b 2d pmf that one would expect to lie on the mfep 
 between the wkwkw(k ) and kwkwk configurations .
the error is also estimated 
 to be very high , especially in the region around where the wwkwk(k ) 
 configuration would be ( figure s8 , supporting 
 information ) .
these difficulties prevent us from drawing any 
 conclusions about the reproducibility of the kvchim-10 2d pmf .
so far we have focused on the feasibility of calculating equilibrium(z1 , z2 , 
 .. zn ) from eq 1 .
now , 
 let us briefly examine how the transmembrane potential changes along 
 the pore axis ; this is the membrane(zi ) term in eq 1 .
the transmembrane potential is usually assumed to drop linearly 
 across the selectivity filter . solving 
 the poisson 
boltzmann equation numerically for both the ( closed ) 
 kcsa and ( open ) kvchim structures shows that only 45% and 
 60% , respectively , of the voltage drop occurs across the selectivity 
 filter ( figure 6 ) .
the value for kcsa is likely to be an underestimate because the 
 intracellular gate is closed , and we note that it agrees with a previous 
 calculation of how much the potential drops across the low resolution 
 closed structure of kcsa . of more concern 

is that not all the voltage drop is predicted to take place across 
 the selectivity filer of kvchim . the poisson 
boltzmann method 
 assumes the lipid bilayer can be represented by a dielectric slab 
 and is therefore less accurate than , for example , simulating the electric 
 field more directly using a double - bilayer arrangement .
this is complex and beyond the scope of this 
 work but deserves further investigation .
boltzmann 
 calculations on the crystal structure 
 suggest that around 60% of the voltage drop occurs across the selectivity 
 filter for the open structure of kvchim , and this is approximately 
 linear . as described in methods
, these calculations 
 were adjusted to take account of the presence of the membrane .
we 
 have investigated the conduction of potassium ions through two 
 potassium ion channels , focusing mainly on a bacterial channel ( kcsa ) but also examining a mammalian chimeric ( kvchim ) channel .
this was done by calculating 2d potential 
 of mean force ( pmf ) maps that capture the energetics of almost an 
 entire conduction event .
the minima in all seven 2d pmf maps corresponded 
 to the arrangements of ions and waters bound in the selectivity filter 
 required for a single - file  knock - on  mechanism ( table 
 s2 , supporting information ) .
these are 
 the wwkwk(k ) , wkwkw(k ) , and kwkwk configurations ( for a definition 
 of the nomenclature , see the legend of figure 1 ) .
because we have enforced the kwk conduction mechanism , this might 
 at first glance appear unsurprising ; however , this does not guarantee 
 that these configurations would be identified as the most stable . 

our simulations predict a qualitative similarity in the conduction 
 of potassium ions by kcsa and kvchim , which is not unexpected given 
 the sequence and structural similarity of both selectivity filters 
 and the surrounding pore domains ( figure 1 ) . 

the functional similarity has been demonstrated by the sensitivity 
 of kcsa to toxins that have evolved to block mammalian voltage - gated 
 k channels .
the most stable configuration 
 of the filter ( kwkwk ) was obtained 
 when three potassiums were bound ( at sites s0 , s2 , and s4 ) with waters 
 at the intervening positions ( s1 and s3 ) in six of the seven pmf maps . 
 with one exception ,
this configuration was 59 kcal / mol more 
 stable than the wkwkw(k ) configuration ( g53 , table 2 ) and 714 
 kcal / mol more stable than the wwkwk(k ) configuration ( g51 ) .
this is consistent with other published 
 computational studies ; this three - ion configuration is estimated to 
 be more stable than the other two - ion states from 0.5 kcal / mol to 5.5 kcal / mol for kv1.2 and from 1 kcal / mol to 2.5 
 kcal / mol or 4 kcal / mol for kcsa .
our simulations indicate that the heights 
 of all the forward and 
 reverse energetic barriers between the stable configurations are 513 
 kcal / mol ( table 2 ) , with three exceptions where 
 the barrier height was only 23 kcal / mol .
barriers of this 
 magnitude are too large to reproduce using a standard electro - diffusion 
 model ; the high conductances ( 97 ps for kcsa at 0 v ) were measured by single - channel electrophysiology . of 
 the published pmf studies , only two determined the heights of all 
 the barriers to be less than 3 kcal / mol .
for the remainder , 
 the heights of the barriers were in the range of 38 kcal / mol , 
 with only a few barriers <3 kcal / mol .
relatively few of these studies estimated error , which is an important 
 consideration given the sensitivity of the conductance to the barrier 
 heights .
one was reproducible 
 despite large differences in the setup of the simulations , and the 
 results of the second set were inconclusive ( figure 5 ) .
finally , not using the cmap correction to the charmm forcefield 
 resulted in substantial changes to the 2d pmfs of both channels , with 
 the relative stabilities of the filter configurations and the heights 
 of the energetic barriers both changing .
this is important , as it 
 suggests that such 2d pmfs are more sensitive to initial conditions 
 and methods than might be anticipated .
we are forced to conclude 
 that it is currently too challenging 
 to accurately calculate the conductance of a potassium ion channel 
 from an experimental structure using classical forcefields , as the 
 heights of the barriers derived from the pmfs are consistently too 
 high and the pmfs themselves are too sensitive to the method employed . 

the last point is important because , as pointed out by bernche 
 and roux , the predicted ion flux is very 
 sensitive to small changes in the predicted heights of any energetic 
 barriers .
one way to avoid the conflation of errors involved 
 in a hierarchical 
 modeling approach is to simply measure 
 the current through a pore domain by 
worryingly , the 
 extensive simulations of the pore domain kv1.2 carried out by jensen 
 et al . dramatically underestimated the 
 conductance for positive membrane potentials , and no conduction events 
 at all were observed for negative values of the membrane potential . 

the membrane potentials were underestimated , a fact corrected in their more recent papers , further exacerbating the lack of agreement with experiment .
this 
 further hints that classical forcefields may not yet be able to accurately 
 model the conduction of ions through a potassium channel selectivity 
 filter .
there are five mutually compatible main reasons that 
 could account 
 for the failure to accurately capture the kinetics of the conductance 
 of a potassium ion channel : ( 1 ) the description of the atomic interactions 
 is inaccurate .
( 2 ) there is insufficient sampling . ( 3 ) the standard 
 electro - diffusion theory is unable to predict current 
( 4 ) the simulations fail to take 
 into account some important factor present in the electrophysiological 
 experiments .
( 5 ) there are other mechanisms in addition to kwk that 
 contribute to conduction . of course
, any combination of these five 
 reasons could also account for the failure to accurately capture the 
 kinetics of the conductance of a potassium ion channel .  
first , as mentioned 
 above , it is likely that the current classical forcefields do not 
 describe well enough the high - field environment within the selectivity 
 filter .
simulating conduction through potassium channels will always 
 be a very hard problem because the ions completely dehydrate upon 
 entering the narrow selectivity filter , and therefore , it is likely 
 that polarization , at the very least , and maybe quantum mechanical 
 effects also , will need to be incorporated into any forcefield before 
 we can accurately model the key interactions between the conducting 
 species and the carbonyl oxygens of the selectivity filter .
second , despite the reproducibility of the kcsa-11 
 2d pmf , it remains possible that our simulations have insufficiently 
 sampled conformational space .
longer simulations and running repeats 
 will also reduce the statistical errors inherent in these calculations . 

third , an applied external electrical field is a nonequilibrium perturbation 
 and could result in structural changes within the pore domain , thereby 
 making eq 1 invalid .
such changes could either 
 make the channel more conductive or alter how the electrical potential 
 changes across the selectivity filter , affecting the heights of the 
 barriers between configurations in complex ways .
alternatively , it 
 is usually assumed that all of the entire transmembrane voltage drops 
 across the selectivity filter ; our poisson  boltzmann calculations 
 suggest this is not the case .
standard electro - diffusion theory also 
 ignores defects arising from the discrete nature of ions .
fourth , it has been shown that the conduction properties 
 of kvchim are sensitive to changes in the pressure across a patch 
 of cell membrane ; therefore , the conditions 
 we have simulated may not match those of the published electrophysiological 
 data .
finally , let us now reflect on how likely it is that there are 
 other mechanisms in addition to kwk that contribute to the conductance .
consider an oversimplified model of the selectivity filter where 
 it is represented as a hollow cylinder , and the ions and waters are 
 modeled as hard spheres with no long - range interactions .
ions and 
 waters would smoothly progress through the filter , resulting in a 
 straight mfep ( figures 3 and 7a ) .
because we observe stable configurations , there must be 
 regions in the filter that can stabilize ions and water , which is 
 obvious from the structures .
the mfeps are curved , from which we infer 
 that the distances between the species in the filter , i.e. , ions and 
 waters do not remain constant during conduction , and hence , the conducting 
 species can be described as  compressible  ( figure 7a ) .
a single file of ions and 
 waters would smoothly move through the selectivity filter if they 
 were hard spheres , leading to a straight line on the mfep .
if the 
 ions and waters are  softer  , then one would expect 
 the mfep to be curved as the approach of ion 3 causes the species 
 in the selectivity filter to be squeezed together .
( b ) if the waters 
 between potassium ions remain trapped when in the selectivity filter , 
 i.e. , their persistence time is infinite , then each conduction event 
 is governed by a single conduction mechanism . if , as we have observed , 
 the persistence time is long but finite , then the situation is more 
 complex because the water in a kwk configuration may leave the filter , 
 leading to other configurations ( kk and k_k ) . in both panels , ions 
 and waters
one usually assumes that 
 a water between two potassium ions in 
 the selectivity filter would remain trapped and therefore have an 
 infinite persistence time ( of course , in practice the water would 
 eventually escape as the ions move through the filter ) .
if true , then 
 each conduction event can be completely described by the order of 
 the conducting species ( e.g. , kwk ) .
waters can occasionally 
 squeeze past a potassium ion and thereby escape the selectivity filter , 
 suggesting that the water persistence time is long ( nanoseconds perhaps ) 
 but not infinite .
second , a variable number of waters can enter the 
 filter ahead of ion 3 , as has been in seen in direct simulations of 
 potassium conductance .
we infer that 
 not only is it possible for the selectivity filter to start in configurations 
 with variable numbers of waters between ions but also the order of 
 the conducting species itself can change during a conduction event 
 ( figure 7b ) .
clearly this complicates matters 
 because a filter that starts in a kwk configuration could end up with 
 a vacancy ( kk , k_k ) or two intervening waters ( kwwk ) ( figure 7b ) .
the implications of this are as follow : ( 1 ) 
 there is no unique simple conduction mechanism .
( 2 ) while it is likely 
  kwk - like  conduction mechanisms dominate , other mechanisms 
 will contribute toward conductance .
( 3 ) it is incorrect to assume , 
 as done in previous studies , that the intervening waters remain trapped 
 and so only the positions of the potassium ions along the pore axis 
 have to be treated as collective variables .
one must either filter 
 the data , as we have done , or also make the positions of waters along 
 the pore axis collective variables and thereby explicitly control 
 their dynamics .
the second option would substantially increase the 
 complexity and cost of any pmf calculation and will almost certainly 
 require the use of more efficient adaptive free energy methods .
finally , the apparent stability of the three - ion 
 configuration kwkwk suggests that four - ion knock - on mechanisms may 
 exist , as has been suggested for kirbac1.1 .
what is the evidence that potassium ions and waters alternate 
 in 
 the selectivity filter and move through it in a single file ( the kwk 
 mechanism ) ?
hodgkin and keynes deduced 
 that ions were conducted through channels in a single file ; however , 
 the evidence that they must alternate with water molecules is less 
 convincing .
in the subsequent structures of a range of potassium ion 
 channels , electron density 
 is observed at well - defined positions along the pore axis , but it 
 is not possible to uniquely infer what mixture of filter configurations 
 leads to the observed electron density because this is an under - determined 
 problem .
for example zhou and mackinnon successfully fitted a simple model to electron density observed 
 along the pore axis , but because the model assumed that either an 
 ion or a water must be present in each site , it can not prove that 
 other mechanisms do not exist .
the kwk mechanism is usually justified 
 through an appeal to electrostatics ( i.e. , 
 it is intuitive that two adjacent potassium ions should be higher 
 in energy than having two potassiums separated by a water molecule ) , 
 but this does not constitute strong evidence given the presence of 
 surrounding oxygen atoms from the filter .
hodgkin and keynes also suggested that vacancies ( empty sites ) 
 may arise during conduction , and this was subsequently investigated 
 by furini and domene .
with hindsight , the calculation of a conduction 
 pmf over 10 years 
 ago was a tour de force , given the scaling 
 of molecular dynamics codes and speed and availability of high performance 
 computers at the time .
codes are more optimized and scale better , and high performance 
 computers are more powerful and numerous .
it has become feasible to 
 repeat the calculation of a pmf map several times ( or try to reproduce 
 a published result ) , and our results 
 show that this is instructive . despite this ,
most published studies 
 do not repeat the simulations that underlie the pmf , although there 
 are some notable exceptions .
neither 
 do most published studies explicitly test the convergence of the resulting 
 pmf nor calculate the statistical error .
this state of affairs is 
 somewhat unsatisfactory , especially where only additional analysis 
 is required . as we have shown , examining the sensitivity of a pmf 
 to the method employed , such as the precise forcefield used , is a 
 useful test of the universality of a result .
adopting these , and no 
 doubt other , best practices will help improve the 
 accuracy , utility , and hence adoption of computational free energy 
 methods in molecular biology .
the historical focus on calculating 
 pmfs of ions conducting through 
 potassium channels and thereby predicting the conductance was primarily 
 due to the determination of high resolution crystal structures of 
 potassium ion channels .
now that structures of other ion channels 
 are being solved , it is likely that sodium channels in particular will be more amenable to this approach . unlike potassium channels , 
 not only
do the ions pass through hydrated but also the conduction 
 mechanism appears to be simpler with fewer ions involved .
indeed , reasonable 
 barriers to conduction have already been calculated for na channels , and even good estimates of the 
 single - channel conductance have been measured by  brute - force  

approaches . in conclusion , assuming 
 the kwk mechanism dominates , we have successfully 
 predicted the expected stable configurations of the selectivity filter 
 but have been unable to recover energy barriers sufficient low to 
 predict the conductance of either kcsa or the pore domain of kvchim . 
 in future work
, we will examine the detailed structural changes that 
 occur during a conduction event .
crystallographic structures 
 of kcsa ( pdb:1k4c ) 
 and kvchim ( pdb : 2r9r ) were obtained . 

the voltage sensors of kvchim were removed , leaving the pore domain 
 ( residues 321417 ) , and missing atoms were added as described 
 in fowler and sansom .
the psfgen package 
 of vmd was used to add the missing atoms 
 to arg117 of kcsa and to protonate glu71 or mutate it to alanine , 
 where necessary .
potassium ions were then placed at the s0 , s2 , and 
 s4 positions with waters at s1 and s3 .
the resulting pore domains were embedded in simple lipid bilayers 
 by conversion into a coarse - grained representation and subsequent 
 self - assembly of a lipid bilayer around the protein as described elsewhere .
all simulations used a simple popc bilayer , with the exception of 
 kcsa-11a , which used a mixed pope / popg bilayer in the ratio 7:3 .
simulations 
 of kcsa had 234 lipids ( 66 popg and 168 pope for kcsa-11a ) and simulations 
 of the pore domain of kvchim had 227 popc lipids .
the small patch 
 of lipid bilayer was solvated with explicit water using the solvate 
 plugin to vmd and neutralizing sodium 
 and chloride ions added making a total of 70 26870 889 atoms 
 for kcsa and 72 354 atoms for kvchim .
the charmm forcefield with cmap 
 correction was used , and therefore , all 
 hydrogens are explicitly modeled .
the energy of the resulting simulation 
 unit cells was minimized for 1000 steps using namd2.7 , and then the temperature was gradually increased 
 from 100 to 300 k in steps of 20 k , running for 40 ps at a time . during 
 this procedure
, the positions of the heavy atoms of the lipid headgroups 
 ( atom names n c11 c15 o2 p1 o3 c1 ) were constrained in the direction 
 of the membrane normal with a spring constant of 10 kcal / mol /  .
the constraints on the lipids were then reduced 10 fold , 
 and a further 60 ps of md was run at 310 k. this and all subsequent 
 simulations were run in the npt ensemble .
the temperature was held 
 at 310 k using a langevin thermostat with a damping coefficient of 
 1.0 ps . the pressure was maintained at 1 atm using 
 either a berendsen barostat ( applied 
 anisotropically with a compressibility of 4.46  10 bar and a relaxation time of 0.2 ps ) or a langevin 
 piston barostat ( with an oscillation 
 period of 200 fs and a damping time constant of 100 fs ) .
electrostatic 
 forces were calculated between charged atoms using the particle mesh 
 ewald method with a grid spacing of 1 
  and a real space cutoff of 13.5  .
van der waals forces 
 were truncated at 12  with a switching function applied from 
 10  .
settle , and shake used to constrain the length of all bonds involving 
 a hydrogen atom , permitting an integration time step of 2 fs .
the 
 dynamics of the system were then simulated until the pore domain had 
 equilibrated , as measured by rmsd .
we follow the method of bernche and roux and define the pore axis as follows .
the origin 
 of the axis is the center of mass of the heavy backbone atoms ( atom 
 names c ca n ) of the residues tvgy in the selectivity filter .
the 
 direction of the axis is defined by a unit vector pointing from the 
 center of mass of the residues tv to the center of mass of the residues 
 gy .
the positions of ion 1 , ion 2 , and ion 3 along this axis ( z1 , z2 , z3 ) were then calculated using linear algebra . as per bernche 
 and roux ,
we simplified the calculation 
 by only following the center of mass of ion 1 and ion 2 , ion 12 ( z12 ) , in addition to ion 3 ( z3 ) .
we assumed then that the configuration of the selectivity 
 filter can be described by a point in ( z12 , z3 ) space . during the equilibration 
 process described above
, harmonic potentials with spring constants 
 10 kcal / mol /  were applied to ion 12 and ion 3 to 
 prevent the selectivity filter changing configuration .
the coordinates 
 of the ions were calculated , and the harmonic forces updated every 
 time step using the tclforces module in namd .
this was also used to 
 write out every 200 fs the coordinates of the heavy atoms of the selectivity 
 filter and the three bound potassium ions and intervening waters . 

the positions of the ions was written to disc every 20 fs and coordinates 
 of the entire system every 1020 ps . each umbrella simulation 
 sampled a region of ( z12 , z3 ) space .
they were either arranged on a rectangular or 
 trigonal grid ( we found the latter to be more efficient ) . during the 
 umbrella simulations , harmonic potentials with spring constants k12 and
k3 were applied 
 to ion 12 and ion 3 , and their coordinates recorded as described above . 

the spring constant k3 was 
 set to 20 kcal / mol /  , whereas a range of magnitudes 
 for k12 was used in the range 2040 
 kcal / mol /  .
occasionally , additional simulations 
 were run with k12 set to 80 kcal / mol /  to improve the sampling of the umbrella simulations in some 
 regions of ( z12 , z3 ) space ( table 1 ) .
the exception was 
 kcsa-11b , where k12 was set to 80 kcal / mol /  for all simulations ; this necessitated reducing the spacing 
 of the simulations in the z12 direction 
 to 0.25  .
a set of initial structures with the selectivity filter 
 in a range of configurations was generated from the equilibrated structure 
 of the channel by permuting the positions of the ions and water molecules . 
 to prevent distortion of the filter ,
the magnitude of the umbrella 
 force was capped at 5.0 kcal / mol /  for the first 
 2050 ps of each simulation . in total ,
1 771 umbrella simulations 
 were run , each for 200 ps , making a total of 0.354 s of simulations .
all the umbrella simulations were then analyzed , and the number 
 of waters between each potassium ion was counted ( figure s3 , supporting information ) using vmd .
a potassium ion was defined as being inside 
 the filter if its pore axis coordinate was within the bounds set by 
 the centers of mass of the hydroxyl oxygen of thr75 ( kcsa or thr371 
 - kvchim ) and the backbone carbonyl oxygens of tyr78 ( kcsa or tyr374 , 
 kvchim ) .
any frame where there was not a single water between both 
 potassium ions and were themselves inside the filter was discarded . 
 to be conservative ,
all subsequent frames from a single umbrella simulation 
 were also discarded because the escape of an intervening water is 
 a slow process .
if either potassium ion was deemed to be outside the 
 filter by the above definition , any number of intervening waters was 
 permitted .
the  backbone 
 dihedral angle of val-76 ( kcsa ) or val-371 ( kvchim ) was also calculated 
 ( figure s4 , supporting information ) using 
 mdanalysis for all 1 771 umbrella simulations . 

the resulting seven data sets ( table 1 ) were first divided into 100 bins , each 2 ps wide . 

the bias was removed , and a 2d potential of mean force for each bin 
 was calculated using an implementation of the weighted histogram analysis method ( wham ) .
the minimum free energy path was found using an implementation 
 of the nudged elastic band method .
the 
 height of the barrier between the first two configurations was measured , 
 as shown for kcsa-11 in figure s1b of the supporting 
 information .
plotting this against time gives a qualitative 
 check of when the calculation has converged ( figure s1c , supporting information ) .
the statistical 
 inefficiency , s(n ) , is defined aswhere x is the 100 ( presumably ) 
 correlated values of the barrier height , which are then binned with 
 width n to produce a second data set of xn values .
if the data set is converged 
 then the ratio of the variances multiplied by the bin width , n , should tend toward a limit as n tends 
 to . we followed the method of yang et al . , and
by penetrating the reversed data set x to different degrees and calculating s(n ) , one can estimate both the convergence and subsequent 
 correlation times .
this is shown in figure s1d and e of the supporting information . from this , we estimate 
 that the first 100 ps should be discarded , and 20 ps is a conservative 
 estimate of the correlation time for the remaining data .
five independent 2d pmfs can therefore be produced for each data set 
 ( figure s6 , supporting information ) and 
 combined in the usual way to estimate errors ( figure s7 , supporting information ) .
all graphs were produced 
 using gnuplot , and all images were rendered using vmd we implemented 
 the nudged elastic band method in python , including several of enhancements suggested by the authors , such 
 as gradually turning on the perpendicular component of the spring 
 force when the path is kinky .
the initial 
 path was formed by connecting 64 beads in a straight line between 
 two points : one just below the wkwkw(k ) configuration ( 1.9,9.4 ) , 
 and the other just above the kwkwk configuration ( 4.3,4.1 ) . 

we investigated whether a second degenerate pathway existing by connecting 
 48 beads between points 3 and 5 taken from the first mfep .
the initial 
 path taken was not a simple straight line but ran via a waypoint to 
 the right of point 3 .
this was usually sufficient for the algorithm 
 to find the second path , path2 .
the poisson  boltzmann 
 ( pb ) equation was solved without and with applied voltage at an ionic 
 strength of 150 mm .
(0 v ) of 
 the potential along the pore axis represents the contribution from 
 the external potential to the total potential .
calculations were performed 
 on a grid of dimensions 100  100  110 grid points and 
 a spacing of 1  .
the position of the transmembrane domain inside 
 the lipid bilayer was extracted from simulations and a dielectric 
 slab of thickness 45  with dielectric coefficient  = 
 2 was set up to mimic the low - dielectric environment of the protein . 

the protein itself was set to  = 4 and the aqueous environment 
 to 80 .
the channel pore was also modeled as aqueous solution . the 
 linear poisson  boltzmann was solved , and the solution used 
 as a starting point for solving the nonlinear pb equation with under - relaxation 
 and  = 0.6 .
 potassium ion channels form pores 
 in cell membranes , allowing potassium 
 ions through while preventing the passage of sodium ions . despite 
 numerous high - resolution structures , 
 it is not yet possible to relate 
 their structure to their single molecule function other than at a 
 qualitative level . over the past decade 
 , there has been a concerted 
 effort using molecular dynamics to capture the thermodynamics and 
 kinetics of conduction by calculating potentials of mean force ( pmf ) . 
 
 these can be used , in conjunction with the electro - diffusion theory , 
 to predict the conductance of a specific ion channel . here , we calculate 
 seven independent pmfs , thereby studying the differences between two 
 potassium ion channels , the effect of the charmm cmap forcefield correction , 
 and the sensitivity and reproducibility of the method . thermodynamically 
 stable ion 
  water configurations of the selectivity filter can 
 be identified from all the free energy landscapes , but the heights 
 of the kinetic barriers for potassium ions to move through the selectivity 
 filter are , in nearly all cases , too high to predict conductances 
 in line with experiment . 
 this implies it is not currently feasible 
 to predict the conductance of potassium ion channels , but other simpler 
 channels may be more tractable . 
Introduction
Results
Discussion
Methods