although 
 it is textbook knowledge that the functions of biomacromolecules 
 are strongly coupled to their conformational motions and fluctuations , computer simulation of such motions has been 
 a challenge for decades .
typically , distinct 
 algorithms are employed to estimate equilibrium quantities ( e.g. , 
 refs ( 3 ) and ( 4 ) ) and dynamical properties 
 ( e.g. , refs ( 510 ) ) . in principle , a single long dynamics trajectory would be sufficient 
 to determine both equilibrium and dynamical properties , but such simulations remain impractical for 
 most systems of interest . aside from straightforward simulations , 
 more technical approaches 
 that can yield both equilibrium and dynamical simulation , sometimes 
 under minor assumptions ,
a number 
 of approaches employ markov state models ( msms ) as part their overall 
 computational strategy . on the basis of replica exchange molecular 
 dynamics ( remd ) ,
it is possible to extract kinetic 
 information from continuous trajectory segments between exchanges 
 and thereby construct an msm .
the adaptive 
 seeding method ( asm ) similarly builds an msm based on trajectories 
 seeded from states discovered via remd or another of the so - called 
 generalized ensemble ( ge ) algorithms .
msms have also been used in combination with short , off - equilibrium 
 simulations to construct the equilibrium ensemble of folding pathways 
 of a protein .
another general strategy 
 is to employ a series of nonintersecting 
 interfaces that interpolate between states of interest selected in 
 advance .
milestoning generates and analyzes transitions between interfaces 
 assuming prior history does not affect the distribution of trajectories .
transition interface sampling ( tis ) and its variants 
 also analyze such transitions and can yield free energy barriers in 
 addition to rates while accounting for some history information .
forward flux sampling ( ffs ) again samples interface 
 transitions : it accounts for history information and can yield rates 
 and equilibrium information .
the  weighted 
 ensemble  ( we ) simulation strategy ( see figure 1 ) , which 
 has a rigorous basis as a path - sampling method , has also been suggested as an approach for computation 
 of both equilibrium and nonequilibrium properties .
although we was originally developed as a tool for characterizing 
 nonequilibrium dynamical pathways and rates ( e.g. , refs ( 5 ) , ( 2528 ) ) , the strategy 
 was extended to steady - state conditions including equilibrium .
the simultaneous computation of equilibrium 
 and kinetic properties using we was demonstrated with configuration 
 space separated into two states by a dividing surface and later for arbitrary states defined in advance of a 
 simulation .
in contrast to many other 
 advanced sampling strategies , we generates an ensemble of continuous 
 trajectories , all at the physical condition ( e.g. , temperature ) of 
 interest .
( a ) ensemble of trajectories 
 with arrow tips indicating the instantaneous configuration and tails 
 showing recent history in the space of two schematic coordinates q1 and q2 . states 
 a and b , shown in gray , are two arbitrary regions of phase space . 

( b ) dissection into two subsets based on whether a trajectory was 
 most recently in state a ( black solid arrows , the    
 steady state ) or state b ( red dashed , the    
 steady state ) .
( c ) statistically equivalent ensemble of weighted trajectories , 
 with arrow thickness suggesting weight .
configuration space has been 
 divided into cells (  bins  ) which each containing an 
 equal number of trajectories . here
, we further develop the capability of we simulation 
 to calculate 
 equilibrium and nonequilibrium quantities simultaneously in several 
 ways that may be important for future studies of increasingly complex 
 systems .
( i ) the approach described below permits the calculation 
 of rates between arbitrary states , which can be defined after a simulation has been completed . in a complex system ,
the most important 
 physical states , including intermediates , generally will not be obvious 
 prior to simulation .
further , the present approach opens up the possibility 
 to use rate calculations to aid in the state - definition process .
( ii ) 
 the non - markovian analysis described here enables unbiased rate calculations 
 in the typical case where  bins  used by we simulation 
 do not exhibit markovian behavior .
the analysis is general and can 
 be applied outside the we context , including the analysis of ordinary 
 long trajectories .
( iii ) the non - markovian analysis can improve the 
 efficiency of we simulations by yielding accurate estimates of observables 
 from shorter simulations .
the analysis is based on a previously suggested 
 decomposition of the equilibrium ensemble into two nonequilibrium 
 steady states .
we is easily parallelizable because it employs multiple 
 trajectories and was recently used with 3500 cores . because there is no need to  catch 
trajectories 
 at precise transition interfaces , we algorithms lend themselves to 
 a scripting - like implementation which has been employed to study a 
 wide range of stochastic systems via regular molecular dynamics , monte carlo , the 
 string strategy , and gillespie - algorithm 
 dynamics of chemical kinetic networks .
we simulation uses 
 multiple simultaneous trajectories , with weights 
 that sum to one , that are occasionally coupled by replication or combination 
 events every  units of time .
the 
 coupling events typically are governed by a static partition of configuration 
 space into  bins  ( figure 1c ) , 
 although dynamical / adaptive bins may be used . in the case of static bins ,
when one or more trajectories enters 
 an unoccupied bin , those trajectories are replicated so that their 
 count conforms to a ( typically ) preset value , m. 
 replicated  daughter 
if more than m trajectories 
 are found to occupy a bin , trajectories are combined statistically 
 in a pairwise fashion until m remain , with weight 
 from pruned trajectories assigned to others in the same bin .
these 
 procedures are carried out in such a way that dynamics remain statistically 
 unbiased .
this study does not adjust 
 weights according to previously developed reweighting procedures during the simulation .
rather , the we simulations 
 described here are long enough to permit relaxation to the equilibrium 
 state .
once 
 the equilibrium state is reached in a we simulation , meaning 
 that there is a detailed balance of probability flow between any two 
 states , equilibrium observables such as state populations or a potential 
 of mean force can be calculated simply by summing trajectory weights 
 in the corresponding regions of phase space .
estimation of observables . to calculate rates , the equilibrium 
 set of trajectories ( figure 1a ) is decomposed 
 into two steady states as shown in figure 1b : the  steady state consisting of trajectories more recently 
 in a than b , and the  steady state with those most recently 
 in b ; these were denoted
 ab  and 
  ba  steady states , respectively , in ref ( 31 ) .
trajectories are  labeled  
 according to the last state visited , i.e. , classified as  or 
  , during a we simulation or in a postsimulation analysis (  post - analysis  ) . 

the direct rate kab estimate 
 is computed from the probability arriving at the final state via1where mfpt is the mean - first - passage time , 
 flux(a  b| ) is the probability per unit time arriving 
 at state b in the  steady state , and p( ) 
 is the total probability in the  steady state . by construction p( ) + p( ) = 1 .
normalizing 
 by p( ) effectively excludes the reverse steady 
 state , and the rate calculation only  sees  the unidirectional 
  steady state as in ref ( 23 ) .
an expression analogous to eq 1 applies 
 for kba . also note that the effective 
 first order rate constant , defined by flux(a  b|)/paeq , can be determined from equilibrium we simulation
because paeq can be directly computed by summing weights in a. we note 
 that analogous direct calculation of observables can be 
 performed from an equilibrium ensemble of unweighted ( i.e. ,  brute 
 force  ) trajectories by assigning equal weights to each . beyond the direct estimates of observables based on trajectory 
 weights ,
we also generalize previous matrix formulations for nonequilibrium 
 steady states into an equilibrium formulation that explicitly accounts for the 
 embedded steady states ( as in figure 1b , c ) . 

these non - markovian matrix estimates are tested below and may prove 
 important for future we studies using shorter simulations , as described 
 in the discussion .
our matrix approach 
 explicitly uses the decomposition of the equilibrium population into 
  and  components for each bin i:2which implies p( ) 
 = ipi and p( ) = ipi. we called this a  labeled  
 analysis .
thus , with n bins , a set of 2n probabilities is required rather than n. similarly , 
 a 2n  2n rate matrix is required : kijv , where  and 
  can be either the  or  subsets of trajectories
each of the previously considered kij rate elements is thus decomposed 
 into four history - dependent elements which account for whether the 
 particular trajectory was last in state a or b and whether the trajectory 
 transitions between the  and  subsets .
the analysis 
 assumes states consist strictly of one or more bins , but this is always 
 possible in a post - analysis without a loss of generality . in other 
 words ,
given the flexibility we have when we define the bins , it is 
 not a real limitation that the states have to be strictly constituted 
 by bins . constructing a labeled rate matrix for unbiased calculations . for 
 purposes of illustration , here state
a consists solely of bin 1 and 
 state b solely of bin 3 . left : a traditional rate matrix with history - blind 
 elements .
the rate kij gives the conditional probability for transitioning from bin i to bin j in a fixed time increment , regardless 
 of previous history .
the element kijv is the conditional 
 probability for the i to j transition 
 for trajectories initially in the  subensemble which transition 
 to the  subensemble , where  and  are either  
 or . the labeled rate matrix correctly assigns the  
 and  subpopulations of each bin , whereas the traditional matrix 
 may not .
we wish to emphasize that this 
 analysis is  non - markovian  
 because we are explicitly including history information ( i.e. ,  
 and  labels ) in the new 2n  2n rate matrix .
once the matrix is built , the steady state 
 observables are obtained using the same mathematical formalism that 
 would be used in a regular markov model .
however , the matrix should 
 be seen as a tool of linear algebra and not as embodying any physical 
 assumptions .
for example , consider a bin in the  intermediate  
 region ( neither a nor b ) , such as bin 2 in figure 2 . in this region , an  trajectory can not change into 
 a  trajectory , nor vice versa ; hence rates for these processes 
 are zero .
similarly , an  trajectory in the intermediate region 
 which enters a bin in b must turn into a  trajectory , so the 
 rate will always be zero to the  components of bins in b. the non - markovian results below stem from the division into  
 and  steady states , but several steps are required . first , 
 rates among bins are estimated in a post - analysis as3where ij is the 
 probability flux , for a given iteration , from bin i to j of trajectories only with initial and final 
  labels   and  , respectively , while i is the population labeled as  which is initially in i. the subscript  2  in the numerator indicates 
 that the rate kijv is estimated to be nonzero 
 only when more than one transition is observed ; after the second event , 
 all events are included , from the first one , to avoid bias .
the requirement 
 for two transitions was found to greatly enhance numerical stability 
 in estimating fluxes and rates between macroscopic states : rates estimated 
 from single events exhibit large fluctuations .
notice that eq 3 is a ratio of averages and differs 
 from the average ratio ij/i , which might seem equally or 
 more  natural .
 however , our data show that eq 3 yields unbiased estimates , while the average ratio 
 may not ( data not shown ) .
the difference between the two estimators 
 indicates that transitions are correlated with trajectory weights . 

perhaps more importantly , the average ratio places less importance 
 on high weight transitions due to the instantaneous normalization  and 

that is , low - weight 
 transitions count as heavily as high - weight events , which evidently 
 biases the rate estimate . in the ratio of averages , high - weight events 
 appropriately count more . to obtain  macroscopic  
 rates between states consisting 
 of arbitrary sets of bins ( noting that arbitrary bins can be employed 
 in a post - analysis ) , we calculate  labeled 
fluxes for 
 use in eq 1 via4the labeled bin populations pi and pi are obtained from 
 the steady - state solution of the labeled rate matrix k = { kij}. a summary of the 
  labeled  or non - markovian matrix 
 procedure for estimating rates between arbitrary states is as follows . 

first , we obtain the labeled rate matrix k = { kijv } using eq 3 to 
 average interbin transitions .
second , we solve the matrix problem kpss = pss , yielding the steady state solution pss .
then , 
 the steady state solution pss along with 
 the labeled rate matrix elements are used to calculate the  
 flux entering state b and the  flux entering a ( eq 4 ) .
finally , the mfpt values are obtained from eq 1 . in the graphs below ,
each non - markovian estimate 
 shown is from the matrix solution using the kij rates calculated based on all data obtained until the given iteration 
 of the simulation .
the non - markovian matrix formulation exhibits 
 a number of desirable 
 properties : ( i ) unlike with unlabeled ( i.e. , implicitly markovian ) 
 analysis , kinetic properties will be unbiased as shown below .
( ii ) 
 solution of both the  and  steady states is performed 
 simultaneously via a standard markov - state - like analysis of the kij rate matrix .
by contrast , if 
 the  and  steady states are independently solved within 
 a markov formalism , there can be substantial ambiguity in how to assign 
 feedback from the target to initial state when the initial state consists 
 of more than one bin .
( iii ) the labeled formulation guarantees , by 
 construction , the flux balance intrinsic to equilibrium , namely , flux(a 

( iv ) the analysis 
 can be performed using arbitrary bins ( and states defined as sets 
 of these bins ) .
it is not necessary to employ the bins originally 
 used to run the we simulation because a post - analysis can calculate 
 rates among any regions of configuration space .
( v ) the analysis is 
 equally applicable to ordinary brute - force simulations . for reference
, 
 we also perform a traditional markov analysis of 
 the trajectories , which will prove to yield biased rate estimates 
 because most divisions of configuration space ( e.g. , we bins ) are 
 not true markovian states . the markov analysis proceeds without labeling the trajectories .
elements of the rate 
 matrix are estimated as5where the subscript  2  again 
 means that we only estimate a rate as nonzero once at least two transitions 
 from i to j have occurred .
bin populations 
 are then computed by solving for the steady - state solution of the 
 markov matrix with elements kij .
the computation of an mfpt requires the use of source 
 ( a ) and sink 
 ( b ) states .
hence , we determine markovian macroscopic 
 rates by substituting the markovian kij for all nonzero elements of the kij. we emphasize that this is merely an accounting 
 trick to establish sources and sinks and simultaneously measure both 
 a - to - b and b - to - a fluxes / rates . we perform a smoothing operation 
 on the macroscopic markovian rates 
 because otherwise the data are fairly noisy .
the mfpt results shown 
 for the markovian matrix analysis are running averages based on the 
 last 50% of the estimates ( where each estimate is from the matrix 
 solution using kij estimates 
 from all data obtained until the particular iteration ) .
we confirmed 
 numerically that such smoothing did not contribute bias to any of 
 the mfpt estimates .
weighted ensemble simulations were performed on two 
 systems : the 
 alanine tetrapeptide ( ala4 ) solvated implicitly and a pair of explicitly 
 solvated methane molecules .
all simulations were performed at 300 
 k with a stochastic thermostat ( langevin thermostat ) .
friction constants 
 of 5.0 and 1.0 ps were used for ala4 and methane 
 systems , respectively . the molecular dynamics time step used for all 
 systems was t = 2 fs .
an iteration is defined 
 to be the simultaneous propagation of all trajectories in the ensemble 
 for some amount of time , . in these studies , a value of  
 = 2500t is used for ala4 and  = 250t for the methane  methane system . for ala4 , 

the all - atom amber ff99sb force field with implicit gb / sa solvent and no cutoff for the evaluation 
 of nonbonded interactions was simulated using the amber 11 software 
 package .
the hawkins , cramer , and truhlar pairwise generalized born model is used , with parameters described 
 by tsui and case ( option igb=1 in amber 
 11 input file ) . the progress coordinates were selected and  binned  
 using a 10  10 partition of a 2d space .
a dihedral distance d = ( ( 1/n)idi )    with respect to a reference set 
 of torsions is used in the first dimension , where n is the number of torsional angles considered and di is the circular distance between the 
 current value of the ith angle and our reference , 
 i.e. , the smaller of the two arclengths along the circumference .
this 
 dimension was divided every 14 from 0 to 126 and then 
 a final partition covering the space ( 126,180 ] ) . in the second dimension , 
 a regular rmsd , using only heavy atoms , is measured with respect to 
 an -helical structure . in this case , the space was divided 
 every 0.4  from 0 to 3.6  and then a final partition covering 
 the space [ 3.6, ) .
values and coordinates for the references 
 used to compute the order parameters are given in the supporting information ( si ) .
the methane 
 molecules were simulated using the gromacs 4.5 software 
 package with the united - atom gromos 45a3 
 force field and dodecahedral periodic 
 box of tip3p water molecules ( about 900 
 water molecules in a 34  34  24  box ) .
van der waals 
 interactions were switched off smoothly between 8 and 9  ; real - space 
 electrostatic interactions were truncated at 10  .
long range 
 electrostatic interactions were calculated using particle mesh ewald 
 ( pme ) summation .
the single progress coordinate was the distance r between the two methane molecules , following ref ( 28 ) .
the coordinate r  [ 0, )  was partitioned with a bin 
 spacing of 1  from 0 to 16  and a last bin covering the 
 space r  [ 16, )  . for 
 the post analysis of methane ,
the coordinate r  [ 0, )  was partitioned so that the first bin 
 is the space r  [ 0,5 )  , then a bin 
 spacing of 2  was used from 5 to 17 , while the last bin covers 
 the space r  [ 17, )  .
the results shown below include all data generated 
 in all trajectories : no transient or relaxation period has been omitted .
for ala4 , 
 populations and mfpts 
 are estimated using we and compared to independent measurements based 
 on ordinary  brute force  ( bf ) simulation .
rates are 
 estimated in both directions between the two sets of states a1,b1 
 and a2,b2 shown in figure 3 ( see si to visualize representative structures ) .
the 
 second set is less populated and consequently expected to be more 
 difficult to sample .
figure 3 also shows the 
 bin definitions used in the post - analysis , which were the same as 
 those used during the we simulation .
however , as we shall see in our 
 second system , we can use any partition of the space for the post 
 analysis . the ala4 free energy surface .
the surface is projected onto two 
 coordinates : d = ( ( 1/n)idi )    from one reference 
 structure ( see si ) and the rmsd with respect 
 to an ideal -helix .
the surface was computed using 3.0 s 
 of ordinary  brute force  simulation .
the set of states 
 a1,b1 is highlighted in green , while the second set a2,b2 is highlighted 
 in red .
the grid shows bins that were used both for we simulation 
 and for the post - analysis calculation of observables via the non - markovian 
 matrix formulation .
the data shown below 
 are based on the same total simulation times 
 in bf and we .
the bf estimates and confidence intervals are based 
 on a single long trajectory of 3.0 s where thousands of transitions 
 between states were observed .
five independent we simulations were 
 run , each employing a total of 3.0 s accounting for all the 
 trajectories . the use of independent we runs permits straightforward 
 error analysis for comparison with bf . as described above ,  direct  

we measurements sum 
 trajectory weights for population and flux calculations . figures 4 and 5 show direct estimates 
 for both equilibrium and kinetic quantities for both sets of states . 

we estimates as a function of simulation time are compared to 95% 
 confidence intervals for bf simulation .
direct we estimates for 
 populations and mean first passage times 
 ( mfpts ) for ala4 states a1,b1 from figure 3 .
five independent we runs are shown , each based on 3.0 s 
 of total simulation time .
dashed lines indicate roughly a 95% confidence 
 interval based on 3.0 s of brute force simulation .
each nanosecond 
 of molecular ( single - trajectory ) time corresponds to approximately 
 200 ns of we simulation including all trajectories in a single run .
direct we estimates for populations and mean 
 first passage times 
 for ala4 states a2,b2 from figure 3 .
five independent 
 we runs are shown , each based on 3.0 s of total simulation 
 time .
dashed lines indicate roughly a 95% confidence interval based 
 on 3.0 s of brute force simulation .
each nanosecond of molecular 
 time corresponds to approximately 200 ns of we simulation accounting 
 for all trajectories in a single run . as with all observables , data from five
is the estimate 
 using all data from the run and thus is based on a total simulation 
 time equal to that of bf ( 3 s ) .
the spread of the rightmost 
 we data points therefore can be compared with the bf confidence interval 
 to gauge statistical quality .
the mean values of the direct 
 estimates are in agreement with bf 
 confidence intervals in all cases
. in some cases , the spread of we 
 estimates is significantly less than that for bf prior to the full 
 extent of we simulation .
each nanosecond of  molecular time  
 in figures 4 and 5 ( i.e. , 
 single - trajectory time ) corresponds to approximately 200 ns of total 
 simulation in a single we run accounting for all trajectories .
hence , 
 in some cases , considerably less we simulation is required for an 
 estimate of the same statistical quality as resulted from the full 
 bf simulation of 3.0 s .
we 
 also show results of the non - markovian matrix analysis for select 
 observables .
figure 6 shows that the non - markovian 
 analysis yields unbiased estimates of the same equilibrium and nonequilibrium 
 properties calculated with direct estimates .
( results for other observables , 
 like the population of a1 and the a1b1 mfpt , not shown , exhibit 
 qualitatively similar agreement . )
the agreement contrasts with a purely 
 markovian matrix formulation , which does not account for the  labeling  
 described above , which can yield statistically biased estimates for 
 kinetic quantities ( see methane results , below ) .
unbiased matrix - based 
 estimates are important when reweighting is used in we as noted in the discussion .
reweighting was not used in the present study , however . population 
 of a2 and mean first passage time for ala4 from a2 to 
 b2 , estimated by the non - markovian matrix analysis of we data .
dashed 
 lines indicate roughly a 95% confidence interval from brute force 
 simulation , as in figures 4 and 5 .
the states are defined in figure 3 . in the methane system , we 
 simulation is used to measure first - passage times based on a range 
 of state definitions . for a complex system , analyzing
the sensitivity 
 of the mfpt to state definitions could aid in the definition of states .
the mfpt was estimated directly , as well as by both non - markovian 
 and markovian matrix analysis . to assess statistical uncertainty , 

the bins used 
 for post - analysis differ from those used in the original we simulation , 
 as a matter of convenience  underscoring the flexibility of 
 the approach .
figure 7 shows passage 
 times measured as 
 a function of the boundary position for the unbound state .
the boundary 
 of the bound state a was held fixed at a separation of 5  while 
 the definition of the unbound state was varied from 5 to 17  . 

the passage times were measured in increments of 2  and compared 
 with bf results as shown in figure 7 .
the bf 
 confidence intervals are based on a single long trajectory of 0.4 
 s , the same total simulation time used in each we simulation .
the mean 
 first passage time for methane association ( b to a ) and 
 dissociation ( a to b ) measured  directly  and from the 
 non - markovian matrix analysis from we simulation as a function of 
 the boundary of state a. the inset displays the pmf along with the 
 definitions of the unbound and bound states , indicated by b and a , 
 respectively .
dashed lines indicate roughly a 95% confidence interval 
 based on 0.4 s of brute force simulation .
figure 7 shows that both direct and 
 non - markovian 
 matrix estimates are in agreement with bf confidence intervals .
for fixed state definitions , figure 8 shows 
 the evolution of state populations mfpts , as was done for ala4 .
we 
 fix the movable boundary position in figure 7 ( inset ) , defining state b as all configurations with r > 11  .
direct and 
 non - markovian 
 we estimates for populations and mean first passage times ( mfpts ) 
 are plotted vs molecular time .
five independent we runs are shown , 
 each based on 0.4 s of total simulation time .
dashed lines 
 indicate roughly a 95% confidence interval based on 0.4 s of 
 brute force simulation .
each nanosecond of molecular time corresponds 
 to approximately 80 ns of we simulation accounting for all trajectories 
 in a single run .
the bound state ( a ) is defined by distances less 
 than 5  , and b is defined by distances greater than 11  . the performance of the non - markovian 
 matrix estimates are particularly 
 noteworthy in figure 8 .
the matrix estimates 
 converge faster than direct estimates to the exact results for the 
 state populations . presumably , this is because the direct approach 
 requires relaxation of the full probability distribution to equilibrium , 
 whereas the matrix approach requires only relaxation of the distribution 
 with each bin ( in order to obtain accurate interbin rates kij ) .
in contrast to the unbiased 
 mfpt estimates obtained by both direct 
 and non - markovian analysis , the markov analysis can be significantly 
 biased for the mfpt .
figure 9 shows that applying 
 the markovian analysis ( section 2.3 ) leads 
 to mfpt estimates clearly outside the bf confidence interval .
data 
 in the si show that the use of a more sophisticated 
 model such as a maximum - likelihood estimator for reversible markov 
 models yields similar results and does 
 not correct the bias .
populations of a ( r < 5  ) and 
 b ( r > 11  ) and mfpts for the methane system , 
 estimated 
 by the non - markovian matrix analysis and the markovian analysis without 
 history information . dashed lines indicate roughly a 95% confidence 
 interval from brute force simulation based on 0.4 s of total 
 simulation time .
equilibrium properties , 
 however , can be estimated without bias 
 in a markovian analysis because history dependence is immaterial . 

figure 9 also illustrates correct ( equilibrium ) 
 population estimates based on the markovian analysis .
to our knowledge , this is the first weighted 
 ensemble ( we ) study 
 using the original huber and kim algorithm to simultaneously calculate both equilibrium and nonequilibrium 
 quantities . the present study estimates observables ( populations and 
 mfpts ) based on arbitrary states defined in a postsimulation analysis , 
 permitting the examination of different state definitions and their 
 effects on observables .
two qualitatively different estimation schemes 
 were examined , including a non - markovian rate - matrix formulation which 
 shows promise for reducing transient initial - state bias ( a bias which 
 is intrinsic to direct estimation of observables based on weights ) . 

both schemes showed substantial efficiency gains for some observables 
 even in the test systems which appear to lack significant energy barriers 
 in their configurational landscapes .
nevertheless , 
 as described below , the present data do point to further challenges 
 likely to be exhibited by larger , more complex systems
. one key feature of the 
 we implementation studied here is the ability to investigate a range 
 of state choices . as
computer simulations tackle systems of growing 
 complexity , it seems increasingly unlikely that states chosen prior 
 to a study will prove physically or biochemically relevant .
indeed , 
 it is already the case that specialized algorithms are invoked to 
 identify physical states , separated by the slowest time scales , from 
 existing trajectories . with we
simulation , as suggested 
 by our methane data , one can adjust state boundaries to minimize the 
 sensitivity of rates to those boundaries . a possible concern 
 with postsimulation state construction
is the need to store a potentially 
 large set of coordinates to ensure sufficient flexibility in post 
 analysis .
as an illustration , storage of { x , y , z } coordinates for 1000 heavy atoms in 
 a we run of 1000 iterations using 1000 trajectories would require 
 10 gb .
the estimation of both equilibrium 
 and kinetic properties from 
 relatively short simulations is an important goal of current methods 
 development , including for we . here
, we have demonstrated 
 as a  proof of principle  that we simulation can do 
 this efficiently ( compared to brute force simulation ) , without bias , 
 in parallel , and with flexibility in defining states .
given the relatively 
 fast time scales ( nanosecond scale ) characterizing the present systems , 
 it is somewhat surprising that we is better than brute - force simulation 
 for some of the observables and never worse .
previous studies suggest 
 that we has the potential for greater efficiency in more complex systems .
once a configuration space is discretized 
 ( e.g. , bins in we simulation ) , one expects in general that transitions 
 among such discrete regions will not be markovian . to take the simplest 
 example ,
in a 1d system , whether a trajectory enters a finite - width 
 bin from the left or right will affect the probability to make a transition 
 in a given direction .
so generally , discretized systems are non - markovian , 
 even when the underlying continuous dynamics are markovian .
this study compared 
 estimation of equilibrium and nonequilibrium observables using the 
 original we algorithm and via post - analysis . as mentioned in the introduction
, the occasional rescaling of weights 
 to match an equilibrium or nonequilibrium steady - state condition was not used to avoid any potential complications .
our data clearly show that a standard markovian analysis of we 
 simulation is inadequate ( figure 9 ) , since 
 we bins typically are not markovian
. additional information  history 
 dependence , as embodied in the / labeling scheme  is 
 needed to obtain unbiased results .
inclusion of history information 
 in the matrix analysis means it is intrinsically  non - markovian  
 regardless of the linear algebra employed .
future work will 
 incorporate the rate estimation and non - markovian 
 matrix schemes developed here , as well as possibly the simpler markovian 
 scheme shown in section 2.3 .
our data ( figure 8) suggest that these could be very successful in 
 bringing a we simulation closer to a specified steady state .
but it 
 is an open question whether reweighting simulations will prove superior 
 to the type of post - analysis suggested here .
importantly , data presented 
 here indicate that some rate estimators could lead to biased estimates 
 for populations , which , in turn , would bias a reweighted simulation .
one practical future approach , suggested by the work of darve and 
 co - workers , could be to define preliminary 
 states in advance to aid sampling transitions in both directions and 
 then to subject the data to the same post analysis performed here 
 to examine additional state definitions besides the initial choices .
the present study has not 
 addressed some of the intrinsic limitations of the we approach , which 
 are the related issues of correlations among trajectories ( due to 
 the replication and merging events ) and sampling  orthogonal  
 coordinates not divided up by we bins . in the systems examined here , 

there was sufficient sampling in orthogonal dimensions to obtain excellent 
 agreement with brute force results in all cases .
however , significant 
 future effort will be required to address correlations and orthogonal 
 sampling , the latter being a problem common to methods which preselect 
 coordinates such as multiple - window umbrella sampling and metadynamics .
in this proof - of - principle 
 study , the parallel weighted ensemble 
 ( we ) approach has been applied to measure equilibrium and kinetic 
 properties from a single simulation in small but nontrivial molecular 
 systems .
importantly , populations and rates could be measured for 
 arbitrary states chosen after the simulation . for all tested observables , 
 unbiased estimates
were obtained , as validated by independent brute - force 
 simulations . in a number of instances , we was significantly more efficient  yielding 
 estimates of a given statistical quality in less overall computing 
 time compared to simple simulation , including all trajectories . in 
 this sense , not only
is we a parallel method but it can exhibit  super - linear 
 scaling ; 
e.g. , 100 cores can yield desired information more 
 than 100 times faster than single - core simulation .
we also developed 
 a non - markovian matrix approach for analyzing 
 we or brute - force trajectories , capable of yielding unbiased results , 
 sometimes faster than direct estimates of observables from we .
the 
 non - markovian formulation also yields simultaneous estimates of equilibrium 
 and nonequilibrium observables based on an arbitrary division of phase 
 space , which is not possible in a standard markovian analysis .
the approaches tested here will need to be further developed and 
 tested in more complex systems .
 equilibrium formally can be represented 
 as an ensemble of uncoupled 
 systems undergoing unbiased dynamics in which detailed balance is 
 maintained . 
 many nonequilibrium processes can be described by suitable 
 subsets of the equilibrium ensemble . here , we employ the  weighted 
 ensemble  ( we ) simulation protocol [ huber and kim , biophys . 
 j.1996 , 70 , 97110 ] 
 to generate equilibrium trajectory ensembles and extract nonequilibrium 
 subsets for computing kinetic quantities . 
 states do not need to be 
 chosen in advance . 
 the procedure formally allows estimation of kinetic 
 rates between arbitrary states chosen after the simulation , along 
 with their equilibrium populations . 
 we also describe a related history - dependent 
 matrix procedure for estimating equilibrium and nonequilibrium observables 
 when phase space has been divided into arbitrary non - markovian regions , 
 whether in we or ordinary simulation . in this proof - of - principle study , 
 
 these methods are successfully applied and validated on two molecular 
 systems : explicitly solvated methane association and the implicitly 
 solvated ala4 peptide . 
 we comment on challenges remaining in we calculations . 
Introduction
Theoretical Formulation
Model Systems and Simulation
Details
Results
Discussion
Conclusions