with 
 relentless advances in supercomputing , and rapid expansion 
 of the number of structurally characterized macromolecules , molecular 
 dynamics ( md ) has evolved into a standard tool for studying biomolecular 
 phenomena . while experimental techniques often provide 
 either a high - resolution static or a low - resolution dynamic picture 
 of the molecular phenomena , md simulations can violate this 
uncertainty 
 relation  between spatial and temporal resolutions and provide 
 a dynamic , yet detailed view at an atomic level .
conventional 
 md simulations , however , suffer from poor conformational sampling , 
 preventing one from achieving an accurate description of conformational 
 ensembles and free energy landscapes . despite the ever - increasing 
 capability of supercomputers ,
the time scale limitation remains a 
 great challenge ; the typical time scale of an atomistic md simulation 
 of biomolecular systems is much smaller than those required to describe 
 most biologically relevant molecular phenomena .
for instance , many 
 biological processes , e.g. , active membrane transport , which rely 
 on large - scale conformational changes of a protein , occur on time 
 scales of milliseconds or longer . over the past few decades , various enhanced sampling techniques 

these methods are often tested initially on relatively small representative 
 molecular systems , e.g. , dialanine peptide .
a handful of these methods are used routinely for the study of more 
 realistic biomolecular systems . from a practical perspective , however , 
 applying many of these advanced methods to complex biological systems 
 is challenging .
ideally , md simulations can be used to characterize 
 thermodynamic 
 and kinetic properties of a biomolecular system / process . from a numerical 
 perspective , these quantities require certain integrations over a 
 high - dimensional phase space  or configuration space , with making 
 certain assumptions  which in turn require a large set of independent 
 and identically distributed samples .
one can also reduce the many - dimensional 
 atomic coordinate space into a much smaller space simply by considering 
 the constraints associated with the steric factors ( e.g. , covalent 
 bonds ) . nonetheless , the curse of dimensionality 
 remains unbeatable unless additional assumptions or simplifications 
 are made .
one particular premise that many free energy calculation 
 and path 
 optimization methods rely on is the existence of a low - dimensional 
 manifold on which lie most of the relevant conformations visited by 
 a system during a structural transition , i.e. , an intrinsic manifold 
 embedded in the atomic coordinate space .
the intrinsic manifold premise relies mostly on the high cooperativity 
 of the biomolecules ; the motions of different parts of a protein , 
 for instance , are correlated along the transition pathway(s ) .
reaction coordinate 
 ( or the committor function ) which , in 
 principle , parametrizes the intrinsic manifold becomes particularly 
 important for inferring kinetic information from configurational integrations . 

note that one may need to make several assumptions with regard to 
 the committor function ( e.g. , velocity independence ) and the transition 
 tube ( e.g. , localization ) in order to connect a computable potential 
 of mean force ( pmf ) to transition rates . in the context of pmf calculations ,
the task of identifying 
 the 
 intrinsic manifold is often done in an ad - hoc manner . 
 in many free energy methods , e.g. , umbrella sampling ( us ) or
its nonequilibrium variations ( such as many 
 flavors of driven or adaptive - bias techniques ) , a low - dimensional collective variable is defined based on the a priori knowledge of the system .
the assumption is that 
 biasing the system along a  good  reaction coordinate 
 will result in samples that are distributed ( although not statistically 
 correct , due to the bias ) along the correct transition tube . in order to identify the intrinsic manifold more systematically
, 
 one may use a statistical learning method such as principal component 
 analysis ( pca ) , isomap , or diffusion map .
dimensionality reduction 
 techniques such as pca and its nonlinear counterparts are often used 
 to analyze md trajectories .
these have also been used in 
 conjunction with enhanced sampling techniques such as metadynamics , adaptive biasing force , and us .
many techniques take advantage of the 
 intrinsic manifold premise 
 without an explicit use of a conventional dimensionality reduction 
 technique .
for instance , path - optimization techniques ( e.g. , different 
 flavors of the string method ) or path - optimizing free energy techniques rely on the 
 existence of a localized transition tube .
these methods simplify the 
 choice of collective variables by allowing the use of many atomic 
 coordinates or collective variables since the sampling is assumed 
 to converge to the relevant regions of the configuration space ( i.e. , 
 the transition tube ) .
the approach discussed here can be categorized 
 as an ad - hoc dimensionality reduction method in that 
 the reaction coordinates 
 and biasing protocols are designed through an empirical search . moreover , 

our approach heavily relies on our knowledge of the system under study 
 in order to limit the conformational sampling to the relevant regions 
 of the phase space while keeping the calculations reliable and accurate . 

the approach and the underlying methodology are particularly tuned 
 toward systems that undergo large - scale and complex conformational 
 changes .
we take advantage of available structural information for 
 the system ( e.g. , the crystal structures of the end states even if 
 they only provide partial information on the atomic coordinates ) to 
 identify a set of system - specific reaction coordinates that can be 
 used to induce conformational transitions . by combining several techniques 
 within an empirical framework ,
the proposed approach could narrow 
 the gap between advanced sampling techniques and realistic applications . 

the general framework of this approach is based on two equally important 
 and novel elements : ( i ) empirical design of nonequilibrium driving 
 protocols and ( ii ) converting the optimized nonequilibrium protocols 
 into equilibrium biasing protocols for free energy calculations .
, we discuss our proposed approach with an example illustrating 
 the workings of the approach .
section 3 will 
 detail some of the technical aspects of the approach for the interested 
 practitioners , while the last section is reserved for a brief summary 
 and conclusion .
it is believed that functionally 
 important , large - scale conformational 
 motions of proteins are associated with a hierarchy of substates and 
 time scales . consider a complex conformational 
 transition in a rugged free energy landscape involving multiple metastable 
 states and several / many characteristic time scales . under these conditions , 

any md - based quantification of thermodynamic properties such as free 
 energy goes hand in hand with some simplification ( e.g. , selecting 
 a few reaction coordinates ) . even in methods that allow for simultaneous 
 application of many collective variables ( e.g. , string method ) ,
the choice of initial conformations / pathways is crucial to what 
 the final results converge to ( e.g. , free energy or transition pathway ) . recognizing these issues
, we propose an approach that takes into 
 account , in a balanced manner , all the factors determining the final 
 results including ( i ) the reaction coordinates used for biasing , ( ii ) 
 the initial conformations / pathways used to initiate the sampling , 
 ( iii ) the enhanced sampling techniques , and ( iv ) the posthoc analysis 
 methods used to quantify the results and estimate the uncertainties .
briefly , there are multiple mostly qualitative 
 ( i ) and mostly quantitative ( ii ) stages involved .
although 
 an accurate quantitative result such as a free energy profile is much 
 more desirable than a qualitative description , a reliable free energy 
 map is often too costly to be obtained . on the other hand
, the more 
 we know about a system and its approximate free energy landscape , 
 the better we can design an efficient free energy simulation . therefore , 
 one may improve the efficiency of free energy simulations by taking 
 advantage of qualitative ( but reliable ) information on the system . 

in particular , experimental clues combined with fast nonequilibrium 
 techniques provide an appropriate framework for this approach .
proposed approach 
 for the study of large - scale conformational changes 
 in proteins .
although the conclusions drawn 
 from the analysis of the results ( stages i.3 and ii.3 ) may be used 
 at any point ( depending on the assessed accuracy ) , iterating the process 
 improves the accuracy of such conclusions .
given the complexity of conformational free energy landscapes , 
 an iterative sampling approach seems necessary in which the results 
 of each iteration ( from imperfect sampling ) are used to design better 
 reaction coordinates , to generate better initial conformations , and 
 to better sample the relevant regions of the configuration space .
one of the fastest 
 ways to generate a physically meaningful pathway 
 between two known conformations for a molecular system is to use a 
 time - dependent external force .
this approach is employed most prominently 
 in simulation techniques such as steered and targeted md which use a harmonic 
 restraint with a center ( defined in a collective variable space ) moving 
 according to a given schedule .
an important disadvantage of 
 these methods is that the resulting trajectories are  far from 
 equilibrium ;  these methods often violate laws of equilibrium 
 thermodynamics as well as those laws formulated for local - equilibrium , 
 near - equilibrium , and steady - state nonequilibrium thermodynamics ( e.g. , 
 onsager reciprocity relations and fluctuation  dissipation theorem ) . in other words , inferring information with 
 regard to the behavior of the system at equilibrium from trajectories 
 generated by steered or targeted md is not straightforward . fortunately , 

externally driven , far - from - equilibrium systems too follow certain 
 laws such as  nonequilibrium work relations . 
nonequilibrium work relations are powerful tools in describing 
 the behavior of systems driven far from equilibrium by the external 
 variation of a parameter . in an md simulation , a time - dependent biasing 
 potential defined in terms of a collective variable
can be used to 
 drive the system from an initial state toward a final one . here
, the 
 collective variable acts as a control parameter , varied according 
 to a protocol . methods such as steered and targeted md , therefore , fall into 
 the category of nonequilibrium driven simulations .
nonequilibrium 
 work relations have been used in many applications 
 to numerically estimate the free energies based on the nonequilibrium 
 work measurements ; however , the use of these relations 
 is not limited to free energy calculations . in principle , nonequilibrium 
 work relations can be used to estimate any equilibrium macroscopic 
 quantity from nonequilibrium driven trajectories . due to the sampling limitations , such a generalization is not necessarily 
 of practical use ; however , one may find certain quantities that can 
 be estimated at a modest computational cost .
as an example , 
 one may estimate the relative transition rates of 
 competing pathways using nonequilibrium work measurements .
the relative transition rates of different paths can be used to 
 estimate their relative importance without requiring an accurate estimate 
 of the whole free energy landscape .
suppose that 
 different  hypothetical  mechanisms associated with 
 a conformational transition can be attributed to distinct transition 
 tubes in a particular collective variable space .
the relative importance 
 of each tube can be estimated using nonequilibrium work measurements 
 based on biasing protocols that guide the system via these transition 
 tubes . along the same lines ,
one may measure the nonequilibrium 
 work along 
 different transition pathways , not for a quantitative description 
 of the transition rates but rather to simply assess , in a qualitative 
 and relative manner , the practicality of different nonequilibrium 
 protocols .
this semiquantitative ad hoc work analysis is one 
 of the novel features of our nonequilibrium approach that will be 
 discussed here .
however , we can only provide some general rules on 
 how to employ the method and how to improve the biasing protocols . 

the approach is intrinsically empirical and knowledge - based and will 
 not result in any improvements without some a priori knowledge of the molecular system under study . in principle , the probability of each transition 
 tube ( which is proportional to its associated transition rate ) can 
 be measured using work measurements given adequate sampling .
however , the presence of distinct work trends between the most relevant 
 transition tube ( which is dominant when the system is not biased ) 
 and other hypothetical transition tubes ( which are disfavored when 
 the system is not biased ) can simplify the calculations . note that due to the nonequilibrium nature of the 
 simulations , one can not make reliable statements based on single trajectories . 

however , the trend of the work ( determined from repeated simulations ) 
 can be used to compare different transition paths / mechanisms . in order to explore the transition tubes ,
one may define a set of 
 relevant reaction coordinates to reduce the configuration space to 
 a coarse coordinate space with a clear distinction between different 
 states of the system including initial , final , and different hypothetical 
 intermediate states .
let us assume that , using these reaction coordinates , 
 one designs two nonequilibrium driven protocols 1 and 2 , each sampling 
 a particular transition tube by driving the system from state a to b through transition tubes 1 and 2 , 
 respectively .
consider a scenario in which there is a clear trend 
 in the work profiles generated by the two protocols such that protocol 
 1 always results in a nonequilibrium work profile whose largest peak 
 is lower than the largest peaks of all work profiles generated by 
 protocol 2 .
one 
 may quickly identify which hypothetical transition mechanism is worth 
 investigating ( mechanism 1 in this case ) .
we note that any parameter 
 involved in the biasing protocol ( e.g. , 
 simulation time ) can influence the trend of the work .
one can simplify 
 the comparison by keeping some of these parameters constant between 
 different protocols associated with different mechanisms / pathways . 

in other words , the protocols should 
 be designed such that they all differ only in one parameter . here 

we consider the trend of the work as a  response  variable ; any parameter that varies between different 
 biasing protocols could be generally considered a candidate 
for instance , if the two protocols use two essentially different collective 
 variables ( e.g. , a distance versus an angle ) , the comparison will 
 be nontrivial ; the different work trends could be due to the way the 
 collective variables are defined ( not due to the difference of the 
 paths taken ) . using different collective variables exploring 
 similar pathways 

will be justified if it provides a way of empirically finding the 
 optimum protocol to sample a given pathway .
therefore , one may take 
 advantage of nonequilibrium driven md and nonequilibrium work relations 
 in an empirical manner , in order to set the stage for a more systematic 
 investigation . in section 3
, we will discuss 
 a few typical techniques for fast comparison of candidate protocols 
 that are not necessarily designed to represent different transition 
 tubes .
we show how nonequilibrium work measurements can be used to 
 assess the quality of a protocol and more importantly to compare it 
 to other protocols in order to define an optimal set of collective 
 variables and/or force constants . with a multidimensional free energy 
 landscape calculated , one can 
 identify and quantify the free energy minima and saddle points .
unfortunately , 
 this is currently too ambitious of a goal for an md - based study of 
 large - scale conformational changes .
most free energy calculations 
 thus effectively sample a localized region of the configuration space . 

using a knowledge - based approach in choosing the biasing protocol 
 involved in the sampling and employing an empirical work minimizer 
 in order to optimize the biasing protocol may increase the chance 
 of sampling along the  relevant  transition tube .
the 
 idea is not exclusive and may be combined with other approaches . for 
 instance
, our best nonequilibrium trajectory can be relaxed into an 
 even more optimized pathway using the string method .
note that the choice of the initial pathway is crucial to the convergence 
 of path - finding algorithms such as the string method .
with appropriate changes , a nonequilibrium driven scheme can 
 be converted into a free energy calculation method .
for instance , 
 repeating the nonequilibrium driven simulations ( preferably in a bidirectional 
 scheme ) may be used to reconstruct the pmf .
we recently proposed an adaptive - bias variation of nonequilibrium 
 driven md which takes advantage of an iterative adaptive biasing potential 
 in order to speed up the convergence .
nonetheless , here we will base our approach on the most popular 
 free energy calculation method , i.e. , us .
due to its similarity to steered md , an us protocol can be conveniently 
 designed based on a fine - tuned steered md protocol by replacing the 
 time - dependent driving potential with a time - independent biasing potential , 
 as will be discussed in more detail in section 3 .
us combined with the weighted histogram 
 analysis method ( wham ) is a standard 
 free energy calculation scheme for reconstructing the pmf along a 
 given reaction coordinate . employing the method to large - scale transitions , 
 however , is often challenging , and simple biasing protocols ( e.g. , 
 using the - root - mean - square deviation ( rmsd ) from a target structure 
 as the reaction coordinate ) usually produces unreliable estimates 
 for free energies . by using system - specific reaction
coordinates and 
 sampling around reliable transition pathways ( obtained using the approach 
 described above ) , one may significantly improve the sampling . sampling a continuous portion of the configuration space along a 
 reaction
coordinate becomes particularly more likely when us is used 
 in conjunction with a replica - exchange scheme , termed 
 here bias - exchange umbrella sampling ( beus ) ( also known as window - exchange or replica - exchange umbrella sampling ) . in replica - exchange md , each replica
is associated 
 with a different value of a given property whose periodic exchange 
 between the replicas based on an  exchange rule  accelerates 
 the exploration of the phase space .
temperature is the most commonly 
 used property to exchange between the replicas which accelerates the 
 sampling of all degrees of freedom somewhat blindly . alternatively
, 
 one can exchange , in a time - dependent or time - independent manner , biasing 
 potentials in a 
bias - exchange  scheme to specifically 
 accelerate the sampling of the degrees of freedom most relevant to 
 the transition of interest .
temperature- and bias - exchange simulations 
 could also be combined to yield a better convergence .
the mixing of the replicas in the bias - exchange method guarantees 
 the continuity of the conformational space sampled ( at least for each 
 individual replica ) , thereby yielding a more reliable free energy 
 profile for the process .
note that due to the large number of degrees 
 of freedom in most biomolecular systems , it is virtually impossible 
 to sample a continuous conformational space if the simulations were 
 to run independently as in a conventional us scheme .
the efficiency 
 of the beus simulations depends on ( i ) the definition 
 of collective variables , ( ii ) the choice of initial conformations , 
 and ( iii ) the choice of the window / umbrella parameters ( i.e. , centers 
 and force constants ) .
the choice of the collective variables ( and 
 initial conformations ) can be improved , e.g. , by using nonequilibrium 
 simulations as described above .
the umbrella parameters can be optimized 
 iteratively using short runs with the goal of achieving roughly similar 
 rates of exchange between neighboring replicas . with certain assumptions , 

along with the factors mentioned above , the choice of the exchange 
 rules plays a role in the efficiency of the beus scheme as discussed 
 elsewhere .
the trajectories generated by ( be)us 
 md simulations are biased 
 by a known biasing scheme that can be used to reweight the samples , 
 recover the correct statistics , and extract ( by proper integration ) 
 information with regard to the transition mechanism
. however , a suboptimal 
 reaction coordinate for biasing can easily result in a poor sampling 
 and unreliable information .
it is thus important to assess the quality 
 of sampling before interpreting the results .
the 
 conventional us / wham approach often assumes a priori knowledge of a good reaction coordinate which is 
 used for both biasing and pmf reconstruction .
however , in the absence 
 of a perfect reaction coordinate , the sampling in the degrees of freedom 
 other than the one used for sampling is important for both determining 
 the reliability of the results and designing better reaction coordinates . 

the data collected from ( be)us simulations along an imperfect reaction 
 coordinate may not give us an accurate free energy profile , but it 
 can be used , particularly in an iterative manner , to arrive at better 
 reaction coordinates and thus reliable free energies . given 
 the computational cost of the ( be)us simulations , it is important 
 to recognize that a flexible analysis framework may be necessary in 
 order to ( i ) maximize the amount of acquired information on the transition 
 mechanism and ( ii ) detect sampling flaws and possibly identify better 
 reaction coordinates .
in both cases , the sampling along the degrees 
 of freedom orthogonal to the one used for biasing plays a key role . 

some of these degrees of freedom are too slow or too fast ; the former 
 can not be identified , and the latter will be sampled properly .
special 
 care must be taken , however , with regard to those events which occur 
 neither too slowly nor too quickly .
the degrees of freedom associated 
 with such events will be sampled poorly , thus needing to be identified 
 and taken into account in interpreting the results .
one of the 
 techniques that simplifies much of the analysis required 
 to derive information on the degrees of freedom other than the one 
 used for sampling is the nonparametric reweighting , which assigns a weight to each 
 individual conformation sampled .
a weighted histogram or other weighted 
 statistical analysis methods can thus be performed once the weights 
 are estimated .
the pmf ( or weighted histogram ) can be reconstructed 
 in terms of any arbitrary measurable quantity ; however , the pmf reconstructed 
 along a poorly sampled coordinate will not be accurate .
assessing 
 the quality of sampling along different coordinates is thus important 
 in order to ( i ) determine the reliability of the results and ( ii ) 
 diagnose the sampling . for a discussion on several variations of nonparametric 
 reweighting schemes ,
one 
 important set of reaction coordinates to investigate is the 
 knowledge - based set as defined in stage i.1 ( figure 1 ) .
the pmf along such reaction 
 coordinates is naturally only accurate within the well - sampled region . 

an unweighted histogram / pmf can be used to readily distinguish between 
 well - sampled and poorly sampled regions .
note that the importance 
 of the continuous sampling is evident here as a gap along any reaction 
 coordinate ( not just the one used for biasing ) indicates unreliable 
 statistics .
a systematic approach in detecting discontinuities ( or 
 outliers ) is clustering , e.g. , using pairwise rmsds .
a relatively 
 less costly clustering approach involves pca ( based on the biased beus samples ) .
a pairwise clustering can then 
 be performed using a metric defined only based on the significant 
 pcs ( i.e. , those with significant contributions to the variance ) . in an n - dimensional beus simulation ,
the first n pcs are likely to correlate well with the n reaction coordinates used for biasing . detecting more than n significant
pcs indicates the existence of slow degrees 
 of freedom orthogonal to the reaction coordinates used for sampling . 

note that for reliable sampling , all slow degrees of freedom must 
 be included in the biasing scheme .
the importance of the degrees of freedom 
 orthogonal to the order 
 parameter for efficient and robust sampling has been known for quite 
 a while , which is related to the so - called  hamiltonian lagging  
 problem .
our approach to addressing this 
 issue is to iteratively ( i ) detect hidden free energy barriers which 
 are in the conformational space orthogonal to the reaction coordinate 
 and ( ii ) incorporate the slow orthogonal degrees of freedom in the 
 biasing protocol .
orthogonal 
 space random walk  ( or more recently  orthogonal 
 space tempering  ) which simultaneously 
 enhances the sampling in the reaction coordinate space and its generalized 
 force space . in order to illustrate the effectiveness 
 and some of the technical aspects of the methodology introduced here 
 in the context of a concrete biomolecular system
, we will take advantage 
 of a set of simulations we have recently performed on a membrane transporter 
 in an explicit lipid / water environment .
we set out for an extensive search on conformational ensemble and 
 structural transitions of msba , a bacterial exporter of the atp - binding 
 cassette ( abc ) superfamily .
we note that 
 our recent report on this membrane transporter did not describe the methodological details of the approach 
 ( the subject of the present paper ) and only emphasized our findings 
 from a mechanistic perspective and in the context of the function 
 of abc exporters . here
, we provide in - depth discussion and analyses 
 that were mostly omitted in our recent report due to the difference 
 in the scope and emphasis .
an abc 
 exporter consists of two nucleotide - binding domains ( nbds ) and two 
 transmembrane domains ( tmds ) .
msba is a unique abc exporter in that its structure has been crystallized 
 in three distinct conformations , including an outward - facing ( of ) 
 and two inward - facing ( if ) conformations with varying degrees of cytoplasmic 
 opening , termed if - closed ( if - c ) and if - open ( if - o ) , respectively . 

the of structure was solved to 3.7  ( pdb : 3b60 ) , while both if 
 structures are solved at a lower resolution ( 4.5  ) , therefore 
 allowing for only c positions to be determined ( pdb : 3b5x and 3b5w ) .
experimental evidence supports an ifof transition 
 during 
 the transport cycle such that the transporter alternates 
 between if and of states (  alternating access  mechanism ) ; however , there is no consensus on what exactly 
 these two states are and what intermediate conformations are involved . 
 from a computational perspective
, conventional unbiased md simulations 
 can not be used directly to study the ifof transition , due 
 to their limited time scales .
however , the crystal structures of msba 
 provide a great framework for the study of the transition by designing 
 relevant collective variables and biasing protocols inducing the transition . as pointed 
 out earlier
, we can use any source of information on the end states 
 to design global reaction coordinates that best capture the molecular 
 motions involved in a transition . for the case of msba ,
although the crystal structures of msba ( if - c and 
 if - o , in particular ) are of low resolution , they can be used to identify 
 the bundling of the helices and their orientations in different states , 
 based on which we defined the following reaction coordinates :  , 
 the relative orientation of the tmd helices describing the cytoplasmic 
 opening ;  , the relative orientation of the tmd helices describing 
 the periplasmic opening ;  , the relative orientation of the 
 two nbds ( see figure 2 ) . note that  
 and  are associated with the tmd conformational changes and 
 can be used to induce opening of the cytoplasmic side and closing 
 of the periplasmic side , respectively . on the other hand ,  

is associated with the nbd conformational changes and can be used 
 to induce the twisting of the nbds .
cartoon representation of msba crystal 
 structures in of , if - c , 
 and if - o conformations along with the definitions of reaction coordinates 
  ,  , and . the panels on right show the projections 
 of the of ( square ) , if - c ( circle ) , and if - o ( triangle ) crystal structures 
 onto ( , ) , ( , ) , and ( , ) spaces . 

color code : nbd is colored yellow / green , and tmd bundles b1 ( tm1,2cis , tm4,5trans helices ) , b2 ( tm1,2trans , tm4,5cis ) , b3 ( tm3,6cis ) , and b4 ( tm3,6trans ) are colored blue , red , yellow , and green , respectively . in 
 the of conformation , the roll axes of bundles b1/b4 ( gcis ) and b2/b3 ( gtrans ) are colored blue and 
 red , respectively , to illustrate the definition of . in the 
 if - o conformation , the roll axes of bundles b1/b3 ( gcis ) and b2/b4 ( gtrans ) are colored blue and 
 red , respectively , to illustrate the definition of . the roll 
 axes of nbd in 
 the of and if - c conformations are colored yellow / green to illustrate 
 the definition of . triplet ( ,, ) describes the global conformational 
 features of the msba transporter in a 3d coarse coordinate space . 

one may think of many other collective variables to be used for biasing , 
 e.g. , the distance between the mass centers of nbds ( dnbd ) .
dnbd was actually used 
 in some of our simulations as well , but to simplify our discussion , 
 we will only focus on  ,  , and  reaction coordinates 
 here .
we used 
 the equilibrated of structure of apo msba to initiate 
 several nonequilibrium driven md simulations using different protocols 
 driving the protein toward the if - o state , including conventional 
 targeted md simulations and orientation - based driven simulations . 

targeted 
 md simulations based on the if - o structure . here , we used  two - stage 
targeted md simulations , 
 targeting the of state toward if - c in the first stage , and if - c state 
 toward if - o in the second stage .
figure 3 compares 
 the results of a typical one - stage targeted md protocol ( red ) to its 
 two - stage counterpart ( blue ) .
comparison of the results of one - stage ( red ) 
 and two - stage ( blue ) 
 targeted md protocols with our ( ,,)-based protocol 
 ( gray ) .
the three 160 ns trajectories are projected onto the ( rmsdif - o , rmsdif - c ) , ( , ) , 
 and ( , ) spaces in panels a , b , and c , respectively , 
 while panel d shows the nonequilibrium , transferred work measured 
 along the trajectories .
note that square , circle , and triangle represent 
 of , if - c , and if - o crystal structures , respectively . in the trajectories resulting from the one - stage 
 targeted md protocols 

the cytoplasmic opening occurs consistently prior to the closing of 
 the extracellular side resulting in a channel - like intermediate that 
 is simultaneously open to both cyto- and periplasm ( see figure 3b ) .
this clearly contradicts the alternating access 
 mechanism and questions the relevance of these trajectories ( knowledge - based 
 assessment ) .
in addition , the final conformation of one - stage targeted 
 md simulations is found to be unstable in subsequent equilibrium simulations , 
 during which the structure opens to both cyto- and periplasmic sides .
the two - stage targeted md simulations typically result in better 
 trajectories in that ( i ) the intermediate conformation does not violate 
 the alternating access mechanism ( see figure 3b ) and ( ii ) the equilibration of both final and intermediate conformations 
 ( i.e. , the ones at the end of first and second stages ) results in 
 locally stable conformations
. more interestingly , one can identify 
 a clear difference in the trend of nonequilibrium work of the two 
 protocols : two - stage simulations always require less work than one - stage 
 ones when the same simulation time is used ( see figure 3d for an example ) . introducing the if - c structure as 
 an intermediate in the targeted 
 md simulations proved to be helpful in improving the protocol both 
 qualitatively ( more consistent with the alternating access mechanism ) 
 and semiquantitatively ( based on nonequilibrium work ) .
nevertheless , 
 the amount of nonequilibrium work still appears to be too high , and 
 it is unlikely that the resulting trajectory provides a reliable pathway 
 that can be used for further quantitative studies , especially for 
 free energy calculations .
we thus use a more systematic approach to 
 sample the reaction - path ensemble of this complex transition .
we use different combinations of  ,  , and  
 reaction coordinates to induce the of  if transition .
assuming the major changes in  ,  , and  occur 
 in discrete stages , a systematic study will be feasible with only 
 six distinct classes of protocols ( all possible orders of these reaction 
 coordinates ) .
the results show that one particular biasing order (  
     ) consistently requires less amount 
 of work than others .
figure 3 compares examples of trajectories obtained from our optimum orientation - based 
 protocol and targeted md simulations .
although the nonequilibrium 
 work provides a semiquantitative tool 
 for comparing different protocols , a knowledge - based assessment of 
 the resulting trajectories is necessary , particularly if one decides 
 to perform free energy calculations based on these trajectories . for 
 instance
, it is quite likely that a protocol involving a short nonequilibrium 
 simulation results in internal conformational distortions such as 
 helix unwinding .
this artifact was observed in some of our shorter 
 simulations ( not shown in figure 3 ) . in addition 
 to the secondary structure
, one may monitor other conformational features 
 of the internal domains for a knowledge - based assessment ( see figure s3 for examples ) .
we note that certain 
 local conformational changes ( e.g. , a side - chain flipping ) may be 
 achieved only through further relaxation of conformations during an 
 unbiased or biased equilibrium simulation ( e.g. , during follow - up 
 beus simulations ) .
one may 
 use our optimum biasing protocol and its resulting nonequilibrium 
 trajectory to initiate a set of beus simulations to estimate the pmf . 

one may think of this process as transforming a nonequilibrium work 
 profile which is dominated by a dissipative term ( e.g. , the gray line 
 in figure 3 ) to an equilibrium free energy 
 profile . within the methodological framework introduced in this paper , 
 equilibrium free energy calculations
can be thought of as the final 
 step in the process of removing the dissipative term from work . in 
 other words ,
nonequilibrium protocol optimization reduces the amount 
 of dissipative work and makes the subsequent free energy calculations 
 more efficient and less costly ; however , the optimized work profiles 
 are likely to be dominated by a dissipative term . given the 
 multistage nature of our msba protocol , we set out to perform free 
 energy calculations only on the final stage ( i.e. , conformational 
 change along  between if - c and if - o conformations ) .
these calculations 
 are particularly useful in identifying the resting state of if msba 
 by characterizing the degree of cytoplasmic opening for the nucleotide - free apo msba in its if state .
figure 4 compares the pmf along  between the values associated with 
 if - c and if - o crystal structures ( i.e. , 
 if - c  16 and if - o  48 , respectively ) to the work profiles associated 
 with the nonequilibrium simulations discussed above ( i.e. , one - stage 
 and two - stage targeted md as well as the optimized protocol ) as plotted 
 against  and offset by the work at if - c  16. here , the pmf was estimated from beus md free 
 energy calculations along  as the reaction coordinate using 
 22 replicas each running for 24 ns , with degrees of opening ranging 
 from 13 to 49 in the  space .
nonequilibrium , transferred 
 work required by 160 ns md simulations 
 performed using one - stage ( red ) and two - stage ( blue ) targeted md as 
 well as the optimized protocol ( gray ) , as plotted against  
 and offset by the work at if - c  16. 
 the work profiles are obtained from the same simulations described 
 in figure 3 , and shown only within the  
 = 1648 (  if - cif - o ) range .
inset : nonequilibrium work obtained from 
 the optimized protocol ( gray ) compared to the pmf along  as 
 obtained from the beus simulations ( see figure 5 ) .
note that the circle 
 and triangle represent crystal structures if - c and if - o , respectively .
figure 5a shows 
 the reconstructed pmf along  within the entire sampled range . 

there exists a local minimum around   16 ( the 
 value associated with the if - c crystal structure ) , while no free energy 
 basin is discernible around   48 ( the value associated 
 with the if - o crystal structure ) , suggesting 
 that the large opening of the cytoplasmic end might be an artifact 
 of crystal contacts .
nonetheless , the deepest minima in the  
 space are in the 2632 range forming a basin that is 
 associated with an if conformation resembling the crystal structures 
 of a homologous transporter protein ( p - glycoprotein ) , which are obtained at higher resolutions . the overall picture emerging 
 from the pmf calculations along  represents msba as a fairly 
 flexible structure in its resting state in the absence of nucleotides 
 and substrates .
we note that while these observations support the 
 relevance of our pmf results , we can not rule out the possibility of 
 a more stable open conformation corresponding to the if - o crystal 
 structure which is not captured in our sampling .
pmf of apo msba in the if state along different 
 reaction coordinates ( a ,  ; b ,  ; c , pc1 ; and d , pc2 ) . 

the pmf is obtained from beus md simulations performed using an -based 
 biasing protocol and estimated using a gibbs sampler ( see supporting information ) .
principal components 
 were constructed from all c atoms of the protein 
 from all conformations used for the pmf estimations .
the mean values 
 and the error bars are estimated from 100 pmfs generated by bootstrapping .
although 
  is the single 
 most relevant reaction coordinate ( among the ones designed ) to capture 
 the conformational changes of msba associated with the opening of 
 the cytoplasmic side in the if state , more information can be extracted 
 from the beus simulations performed along  by reconstructing 
 the pmf along other reaction coordinates . among them ,  is of 
 particular interest , not only from a mechanistic perspective but also 
 to assess sampling efficiency .
this is due to a flexibility associated 
 with  ; diffusion along  is fast enough to be sampled 
 during 24 ns of beus simulations .
however , sampling along  
 is not homogeneous since it was not included in the biasing protocol . 

figure 5b shows the pmf along . noting 
 that in all initial conformations  is between 60 and 
 80 , it is not surprising to have the smallest estimated errors 
 around this region ( i.e. , better sampling ) while the largest estimated 
 errors are for the regions with  < 35 ( i.e. , poor 
 sampling ) . from a sampling efficiency assessment perspective
, 
 instead of choosing arbitrary reaction coordinates to analyze , one 
 may use pca .
we thus performed pca based 
 on all c atoms of the protein using all conformations 
 sampled by beus simulations .
table 1 shows 
 the fraction of variance explained by first four pcs along with the 
 correlation between the pcs with a select number of other reaction 
 coordinates . as expected , pc1 correlates well with  , which 
 was used for sampling ( correlation coefficient of 0.97 ) .
the pmf along 
 pc1 ( figure 5c ) resembles the pmf along  
 as well .
pc2 , 
 on the other hand , correlates strongly with  ( correlation coefficient 
 of 0.82 ) , and the pmf along pc2 ( figure 5d ) 
 resembles that along . pc2 explains about 22% of the variance .
pearson correlation coefficient 
 between a principal component ( column ) and a collective variable ( row ) .
the pca results suggest that 
 biasing the system along  alone 
 might not be enough to achieve good sampling in the pc2 space ( which 
 correlates strongly with  ) . meanwhile , the pmf along the ( , ) 
 space , although not accurate in some regions , may provide clues on 
 the mechanism of transition ( figure 6a ) .
 
 tends to  un - twist  in the intermediate  range , 
 a behavior also observed in our earlier equilibrium simulations .
( a ) contour plot of the pmf in the ( , ) space 
 ( in 
 kcal / mol ) as obtained from the -based beus simulations .
( b ) 
 same as a with a pathway generated by the lfep algorithm .
( c ) same as a but in the ( pc1,pc2 ) space ( as 
 defined in figure 5 ) .
inset : the fraction of 
 variance ( of unweighted beus - generated conformations )  explained  
 by the first five principal components ( bars ) and their cumulative 
 value ( gray line ) .
( d ) same as c with two pathways generated by the 
 lfep algorithm .
inset : free energy along the two lfep - generated pathways 
 ( black and gray ) compared to the pmf along  ( red ) . an if - c 
  if - o reaction coordinate is defined by transforming either 
  ( red ) or arc length of the path ( black and gray ) to make the 
 comparison between the three easier . in the next iteration , one may use a 2d beus md simulation 
 involving 
 both  and  in order to ensure good sampling along both 
 dimensions .
if the 2d beus is not feasible  it could be too 
 computationally demanding  the resulting approximate 2d pmf 
 provides a way of identifying important pathways , e.g. , the pathway 
 shown in figure 6b , found using the lowest 
 free energy path ( lfep ) algorithm , and 
 restraining the replicas along a 1d path .
if multiple competing pathways 
 coexist , one may perform several 1d beus simulations , instead of a 
 2d one .
obviously the 2d beus simulation will be more informative 
 but also computationally more costly .
better sampling along  
 not only improves the pmf estimate along  but also makes our 
 free energy estimate along  more accurate . however , note that 
 the definition of free energy along  will depend on the way 
  was incorporated into the protocol , i.e. , a second dimension 
 for sampling or a degree of freedom to be restrained .
the fraction of variance explained by the first five principal 
 components and their cumulative value are also shown .
the free energy 
 along the two lfep - generated pathways as obtained from our reconstructed 
 pmfs is compared to the pmf along pc1 .
the 1d pmf is significantly 
 lower than the free energy along either pathway , where the sampled 
 region is elongated along pc2 due to the presence of multiple metastable 
 states .
besides 
 pca and other analyses targeting the global features of the conformations , 
 one may be interested in more local conformational changes , particularly 
 those correlating with large - scale conformational changes .
a 2d pmf 
 can be constructed in terms of  ( or pc1 ) , representing the 
 global conformational change , and any quantity of interest representing 
 a more localized change . as an example , here we discuss a potentially 
 significant salt bridge between residues d252 and k299 .
figure 7 plots the 2d pmf in the ( ,d252k299 ) space in which d252k299 counts the number 
 of d252k299 salt bridges formed in the cis and trans monomers .
conformations associated with 
 several minima are presented in figure 7 . if - c 
 and if - c * are both similar to the if - c crystal structure ; however , 
 they differ in their number of salt bridges . if - o * and if - o are also similar to the if - o crystal structure ( although both less 
 open ) , exhibiting different numbers of salt bridges
. contour plot of the pmf 
 in the ( ,d252k299 ) space ( in kcal / mol ) 
 as obtained from the -based beus simulations . 

d252k299 smoothly counts the number of d252k299 
 salt bridges ( in the cis and trans monomers ) .
b is the distance between ( i ) the mass centers of 
 the side - chain oxygen atoms of residue a ( d252 ) and 
 ( ii ) the side - chain nitrogen atom of residue b ( k299 ) . 

inset : a smooth step function (r ) used to 
 quantify the formation / breakage of a salt bridge based on rd252k299 . select conformations associated with several minima 
 are presented in a cartoon representation .
tm5 ( green ) 
 and tm6 ( orange ) helices are highlighted along with d252 
 ( red ) and k299 ( blue ) residues in van der waals surface representation .
the 2d pmf reveals a strong ( nonlinear ) 
 correlation between the 
 two global and local quantities  and d252k299 .
the larger the  , the less likely it is for the salt bridges 
 to form ; however , it turns out that the largest basin ( associated 
 with if - s ) which spans over a large  range is associated with 
 an asymmetric salt - bridge pattern .
this implies that msba in the if 
 state ( at least in the apo form ) is more likely to 
 have only one of the salt bridges than both or none .
this observation 
 is in line with recent experimental work suggesting an asymmetrical 
 tmd arrangement in msba . beside 
 the mechanistic insight obtained from the 2d pmfs , e.g. , 
 the one shown in figure 7
, such multidimensional 
 pmfs can be used more quantitatively to arrive at a better estimate 
 for the transition barrier .
note that the pmf maximum along any collective 
 variable , if estimated accurately , only provides a lower bound to 
 the actual barrier .
this is due to the degeneracy associated with 
 collective variables which becomes evident in our example by comparing 
 the pmf along  with that in the ( ,d252k299 ) space .
there appears to be an approximately 2-kcal / mol barrier 
 around   20 in the 1d pmf ( figure 5a ) ; however , the   20 point is associated 
 with both a saddle point ( d252k299  
 1 ) and a basin ( d252k299  2 ) , implying 
 that the 2-kcal / mol estimate is dominated by the free energy of a 
 local minimum and not the transition state .
the 2d pmf indicates that 
 the free energy associated with the saddle point at   
 20 is around 4 kcal / mol rather than 2 kcal / mol .
although both 
 4-kcal / mol and 2-kcal / mol values can be considered as lower - bound 
 estimates for the free energy barrier , the 4-kcal / mol estimate is 
 a larger lower - bound , which makes it more informative .
note that , 
 however , there is a trade - off between the accuracy and precision since 
 including more dimensions in pmf reconstruction increases the uncertainty 
 of the results due to limited sampling .
in section 2 , we gave an overview of our 
 proposed scheme to combine several sampling techniques in order to 
 optimize the reliability and relevance of free energy calculations . 

here , we will discuss some of the more technical aspects of the approach 
 which could be of particular interest to md practitioners .
although 
 some of the techniques described here are already known to practitioners , 
 as referenced throughout the discussion , certain practical difficulties 
 arise  or become more serious  when one tries to combine 
 different techniques .
the emphasis here is on having a balanced approach 
 toward addressing such practical issues by taking into account the 
 entire process of designing , performing , and analyzing the simulations .
nonequilibrium 
 work measurements can be used for a qualitative comparison of different 
 protocols , as described in section 2 .
however , 
 such comparisons are limited to the sampled transition mechanisms 
 and can not rule out the possibility of alternative transition mechanisms 
 that are either not sampled or sampled only poorly due to the use 
 of a particular biasing protocol .
another important limitation arises 
 when multiple transition protocols are found with similar trends of 
 nonequilibrium work .
if there is a clear trend in work profiles of 
 different mechanisms favoring a number of paths over the others , one 
 can make general qualitative statements about their relative importance ; 
 however , the transition paths with similar trends of work can not necessarily 
 be considered similarly important . for
these cases , longer simulations 
 or more iterations might establish a difference in their trends of 
 nonequilibrium work . despite the considerations mentioned above
, one may take advantage 
 of nonequilibrium driven md simulations and nonequilibrium work measurements 
 in an empirical manner as an efficient way to optimize the biasing 
 protocols .
we note that before assessing the quality of a protocol 
 based on the amount of work , one needs to determine whether or not 
 the employed protocol is effective in inducing the transition of interest . 

this may be a trivial task for simple molecular systems , but in the 
 presence of many metastable states , it is quite likely for the system 
 to get trapped in an irrelevant metastable state without ever getting 
 close to the product .
the rule of thumb is that once the system is 
 steered along the designed reaction coordinate , followed by a careful 
 equilibration , the product must be locally stable while verifiably 
 resembling the target structure .
the careful equilibration here includes 
 ( i ) restraining the system around the final target , ( ii ) slow removal 
 of the restraint , and ( iii ) unrestrained equilibration . in cases in 
 which
the target is not a well - known state , our low - resolution knowledge 
 of the target may be used ; however , the stability of the product should 
 be established before performing costly computations such as free 
 energy calculations .
the amount 
 of nonequilibrium 
 work might be sensitive to the choice of force constant when a harmonic 
 bias is used . the force constant needs to be large enough to induce 
 the transition .
if the force constant is too large , however , the simulation 
 could become unstable , and the molecular system could undergo deformation 
 or distortion . within these limits ,
any force constant may be used 
 to induce the transition of interest ; however , choice of the force 
 constant will influence the amount of dissipation .
nonequilibrium 
 work measurements may be used in order to find an optimum value for 
 force constant in an empirical manner .
we take advantage of the fact 
 that different force constants result in different work distributions ; 
 thus , one may find an optimum force constant which typically ( or on 
 average ) results in a lower work .
see section 
 s2 for a more detailed discussion and an example for this use 
 of work measurements .
suppose that one 
 is interested in inducing a particular conformational transition such 
 as an interdomain orientational change .
for instance , it should induce the desired conformational 
 change without distorting the system , that is , without taking the 
 system through irrelevant molecular configurations ( e.g. , unlikely , 
 high - energy states ) .
if several collective variables satisfying this 
 criterion are available , one may use nonequilibrium work measurements 
 in a semiquantitative manner to choose the best one .
thus , one can 
 try different definitions and compare the results , first , by a qualitative 
 knowledge - based assessment of the conformations sampled and , if there 
 is no clear preference , by comparing the work profiles .
similarly , 
 one may use qualitative comparison of work profiles ( along with knowledge - based 
 comparison of the generated trajectories ) in order to optimize the 
 choice of atoms directly affected by the biasing potential .
for instance , 
 either all heavy atoms , all backbone atoms , or only c atoms can be subjected to a biasing potential to induce an orientational 
 change in a helical domain .
nonequilibrium work measurements again 
 provide a basis for optimizing the details of a protocol in a systematic 
 fashion . in order to induce a conformational transition , 
 it is often relevant 
 to define a set of collective coordinates such that by applying appropriate 
 forces on the system , one can vary these collective variables and 
 change the conformation of the system .
any measure / metric distinguishing 
 between different conformations can be used for the analysis of a 
 trajectory ; however , it does not necessarily make a practical collective 
 variable suitable for biasing / driving a system from one state toward 
 another .
metrics used for the analysis are often chosen to quantify 
 certain properties in an intuitive manner , but they do not need to 
 be , for instance , differentiable . on the other hand ,
a collective 
 variable has to be well - behaved to be used in a biasing potential 
 since its gradient and time derivative , for instance , are needed to 
 calculate the biasing force and work , respectively .
one particular 
 collective variable that is widely used in the context of structural 
 transition of proteins is rmsd .
although using the rmsd from a target 
 structure as a collective variable ( e.g. , in a targeted md simulation ) has proven useful , the method has 
 its own pitfalls and limitations .
rmsd is associated with both extreme 
 degeneracy and large entropy loss ( for 
 large and small values of rmsd , respectively ) .
the trajectory generated 
 by targeted md represents a pathway along which the rmsd decreases 
 almost monotonically and nearly linearly .
another limitation that 
 reduces the flexibility of the method to a great extent is its high 
 sensitivity to the quality / resolution of the target structure . in 
 addition
, a targeted md simulation typically requires a large amount 
 of work to induce a transition , thus making the interpretation of 
 its results difficult in the context of nonequilibrium work relations .
other conventional collective variables such as distance and radius 
 of gyration have their own limitations in inducing complex structural 
 transitions .
collective variables such as rmsd and radius of gyration 
 can be considered special cases of a  generalized distance  
 which has been studied with regard to some of its features ( e.g. , 
 its entropy loss ) when used in a biasing potential .
one particular feature that seems to best describe 
 a variety of 
 large - scale conformational changes in proteins is semirigid - body domain 
 orientational changes .
there are several ways of defining a collective 
 variable that quantifies an orientation - based conformational change . 
 among them ,
the orientation quaternion technique has proven successful as a well - behaved , 
 flexible method for defining system - specific collective variables , 
 specifically aimed at inducing interdomain orientational changes .
note that here we use the term system - specific reaction coordinates 
 or collective variables for definitions that rely on our knowledge 
 of the conformational changes involved in a transition . if the conformational 
 changes are spatially localized ( e.g. , a side - chain flipping or a 
 salt - bridge formation / breakage ) , one may use relatively simple techniques 
 to define a reaction coordinate on a select number of atoms .
if the 
 conformational change is global , more advanced collective variables 
 such as the orientation quaternion may be needed to induce the transition 
 of interest .
the selection of atoms ( grouping of domains , bundling 
 of helices , etc . ) in the definition of the collective variables is 
 again a knowledge - based component .
we note that when a local conformational 
 change triggers a global conformational change , using collective variables 
 defined on a limited number of atoms may prove more useful in inducing 
 the transition .
using combinations of ( generalized ) distance - based 
 and/or orientation - based collective variables may be necessary in 
 order to induce a complex global transition .
suppose 
 that the relative orientation of two molecular domains change during 
 a transition .
the angle between the two domains can be defined as 
 a simple geometric angle based on the mass centers of three groups 
 of atoms ( three - group definition ) .
similarly the arccosine of the 
 dot product of principal axes ( usually the roll axes ) of two groups 
 of atoms can be used to define the orientation angle ( vector - based 
 definition ) .
unfortunately , using a simple three - group definition 
 of an angle to induce a global conformational change often results 
 in undesired deformations of the protein . using the principal axes 

is associated with certain technical difficulties as well ; e.g. , the 
 three principal axis components may interchange during the simulation 
 due to conformational changes .
we will discuss an alternative technique 
 to define orientation - based collective variables , namely , the orientation 
 quaternion which has been recently implemented in several md packages by fiorin 
 et al .
this technique allows for a flexible 
 and reliable definition of system - specific collective variables , specifically 
 aimed at inducing interdomain orientational changes .
the orientation 
 quaternion , often 
 used for  optimal superposition  in computational biology , is a tool to deal with the so - called  absolute 
 orientation  problem .
n ) , we have two different sets of measurements : { xk } and { yk}. to simplify the problem , we assume both sets have been already 
 shifted to bring their barycenters to the origin ( optimum translation ) . 
 to find the optimum rotation to superimpose { yk } on { xk }
a quaternion is a four - component vector that can be 
 considered as a composite of a scalar and an ordinary vector , or as 
 a complex number with three different imaginary parts .
a quaternion 
 whose scalar part is zero is called pure ( reminiscent of pure imaginary 
 numbers ) .
the optimal rotation can be parametrized by a unit quaternion , q that minimizes qxkq *  yk 
 in which  denotes an average over k , q * is the conjugate of q , and 
 q  qq * ( see ref ( 101 ) for 
 more details ) .
the optimal rotation unit quaternion ( or orientation 
 quaternion ) q can be written as ( cos(/2 ) , 
 sin(/2) ) in which  and  ( a unit vector ) are the optimum angle and axis 
 of rotation , respectively . as a collective variable
, an orientation 
 quaternion can be used 
 not only to monitor rotational changes but also to apply force ( which 
 is proportional to the gradient of the quaternion ) on the system in 
 a practical way to induce the desired rotational changes .
suppose 
 that we are interested in inducing a particular rotation  given 
 by its axis of rotation ( unit vector  ) and 
 its target angle of rotation target  on a 
 particular segment of a biomolecule , e.g. , part of a helix , an entire 
 helix , or a bundle of helices .
this can be done by using a time - dependent 
 harmonic potential similar to steered md in spirit:1here , qref({xk } ) is the optimum orientation 
 quaternion to superimpose { xk } on a reference set { xkref}. the reference could 
 be the initial , target , or any other structure . here , to simplify 
 the notations
q(t )  ( cos((t)/2 ) , sin((t)/2) ) is a unit quaternion that is varied externally , providing the center 
 of the harmonic potential at time t during a simulation 
 ( 0  t  t ) .
if the 
 reference is the same as the initial structure , (0 ) and (t ) can be set to 0 and target , respectively . 

once we have q(0 ) and q(t ) , we can use different interpolation methods to determine q(t ) .
a simple method is varying (t ) linearly , which is a special case of the spherical linear 
 interpolation ( slerp ) method .
the particular 
 method discussed here  implemented by fiorin et al.is based on the linear interpolation 
 of the quaternion q(t ) , using the 
 current and final centers , followed by its normalization at each time step .
finally (p,q ) is the length of the geodesic between two points 
 on the unit sphere , transformed by p and q from an arbitrary point on the unit sphere .
an 
 approximate estimate for (p,q ) is known to be arccos(|pq| ) in which pq is the inner product of p and q. the nonequilibrium work along 
 a trajectory generated by the quaternion - based 
 biasing protocol ( eq 1 ) can be measured via2which is the accumulated work at time t. the nonequilibrium , transferred work can also be measured by subtracting ub(qref({xk } ) , t ) 
ub0 from the 
 accumulated work , in which ub0 is the biasing potential measured 
 at t = 0 .
one can collect the biasing potential ( ub(qref , t ) ) and its partial time derivative ( ub(qref , t)/t ) based on the instantaneous qref at time t. for the particular qref schedule , q(t ) that comes from the linear interpolation of the quaternion ( q(t + t ) = q(t ) + ( q(t)q(t))t/(t  t ) ) followed 
 by its normalization ( q(t + t ) = ( q(t + t)/q(t + t) ) , one can show:34 we use our transporter 
 example to demonstrate how orientation quaternions 
 can be used to induce rotational changes .
let us recall the definition 
 of  and  ( defined on the tmd helices ) as well as  
 ( defined on the nbds ) .
suppose tmicis and tmitrans denote the ith transmembrane 
 helix of the two monomers ( labeled cis and trans , arbitrarily ) and consider four relatively rigid bundles , 
 b1 ( tm1,2cis , tm4,5trans helices ) , b2 ( tm1,2trans , tm4,5cis helices ) , b3 ( tm3,6cis helices ) , 
 and b4 ( tm3,6trans helices ) , colored in figure 2 blue , red , yellow , and green , respectively .
 
 describes the angle between two groups of bundles b1/b3 ( or gcis ) and b2/b4 ( or gtrans ) . on the other hand , 
  describes the angle between b1/b4 and 
 b2/b3 ( or gcis and gtrans , respectively ) .
the bundling of the helices 
 and the grouping of the bundles as explained here are clearly inspired 
 by the crystal structures of msba . in 
 order to induce conformational changes along  ,  , or 
 
, one may define several orientation quaternions based on ( i ) 
 individual helices , ( ii ) different bundles of helices , or ( iii ) segments 
 of helices and even combine them in different ways .
we have discussed 
 this empirical process in section s2 ( see supporting 
 information ) .
our empirical process led us to use six orientation 
 quaternions as collective variables to describe the three - dimensional 
 ( ,, ) space including ( qcis , qtrans , qcis , qtrans , qcis , qtrans ) . for each angle
qcis and qtrans are 
 the orientations of the two groups of bundles gcis and gtrans , respectively .
similarly , qcis and qtrans are 
 the orientations of gcis and gtrans , respectively .
the connection 
 between harmonic - based nonequilibrium driven md and harmonic - based 
 ( be)us is quite evident . one can easily design a number of time - independent 
 biasing potentials along the time - dependent biasing potential as the 
 umbrella window potentials and use the generated nonequilibrium conformations 
 close to the window centers as initial conformations .
the force constants 
 and the window centers may be adjusted in order to achieve sufficient 
 overlapping of the sampled windows and , in the beus scheme , reasonable 
 mixing of the replicas .
it is shown that placing the window 
 centers equidistantly defined with a proper thermodynamic metric ( i.e. , 
 along the geodesic ) minimizes the variance in the estimate of free 
 energies and optimizes the mixing of replicas ( within a given protocol ) .
similarly , within the linear response regime , 
 a nonequilibrium driven protocol can be optimized by minimizing the 
 dissipation or the thermodynamic divergence by moving along the geodesic 
 with a constant speed .
the two are indeed 
 related and one may use the information obtained from one ( e.g. , the 
 thermodynamic metric ) to optimize the other . in other words , the overlap 
 ( mixing ) required for the convergence of ( be)us simulations is related 
 to the average dissipative work measured from the driven simulations . 
 roughly speaking , one may conclude that the more dissipation that 
 is observed in a driven protocol , the more expensive the free energy 
 calculations from ( be)us simulations would be .
suppose 
 that we are interested in sampling the transition between a and b state along the protocol . 
 we propose the following steps to prepare and perform beus simulations .
1 . nonequilibrium driven md . in order to decrease 
 the risk of discontinuous sampling in beus simulations , we propose 
 generating the initial conformations from a continuous trajectory . 

we assume 
 we have already optimized the protocol . starting from a , we drive the system along  toward b and generate m conformations along the protocol .
efficiency of the sampling 
 in beus simulations relies on both ( i ) mixing of the replicas and 
 ( ii ) overlap between the windows .
suppose that the exchange rate between 
 replica i and i + 1 is ri , i+1 ( for i = 1 , ... , n  1 , in which n is the number of replicas ) .
the optimum exchange rate r = ri , i+1 depends on several factors , e.g. , the exchange scheme .
however , the flatness of the ri , i+1 function ( i.e. , ri , i+1 = r , where r is constant ) is particularly 
 important to achieve an efficient diffusion along the reaction coordinate 
 ( assuming no phase transition involved ) . with regard to the overlap 
 between the windows , a critical criterion ( particularly when weak 
 force constants are applied )
note that having an equal exchange 
 rate between the replicas does not guarantee a flat sampling along 
 the reaction coordinate .
therefore , when adjusting the parameters , 
 one needs to first make sure that there are no nonexchanging / nonoverlapping 
 windows .
one may use the following procedure prior to production 
 runs of 
 beus simulations to prepare the initial conformations and umbrella 
 potentials using a trajectory generated in a prior nonequilibrium 
 simulation : ( a ) take n initial conformations 
 from the m nonequilibrium conformations . at first
( b ) based on each conformation i ( selected above ) , 
 identify a center i along the reaction 
 coordinate . for each umbrella
i = 1 , ... , n , design a biasing potential ( or umbrella potential ) ui restraining  around 
 i with harmonic constant ki .
( c ) perform short 
 beus md simulations starting with the initial 
 conformations and using the umbrella potentials obtained from steps 
 a and b , respectively .
( d ) iterate steps a to c with different 
 initial conformations ( n can be varied as well ) and 
 harmonic constants until ( i ) 
 the exchange rate between any two neighboring replicas was estimated 
 to be in a given range ( e.g. , 2040% ) and ( ii ) the  
 space ( in a given continuous range ) is expected to be sampled without 
 any gap .
this process may be done empirically ; however , a more automatic 
 process can also be implemented .
the final set of initial conformations 
 and umbrella potentials generated above can be used to perform longer 
 beus simulations to achieve convergence . while the mixing 
 of the replicas in a beus scheme decreases the likelihood of discontinuity
the continuity of the regions sampled 
 by each replica is evident ; however , there could be discontinuity 
 introduced by the initial conformations such that some of the replicas 
 never cross the regions sampled by other replicas .
suppose one picks 
 the initial conformations of a beus simulation from several independent 
 simulations or a single simulation much longer than the individual 
 beus simulations .
this introduces a high risk of nonoverlapping replicas 
 which all sample perfectly well along the reaction coordinate ( used 
 for sampling ) with a reasonable mixing ( in terms of exchange rate ) , 
 but they never actually cross each other along a degree of freedom 
 orthogonal to the reaction coordinate . although clustering will detect 
 such problems , it is important to avoid risky choices for initial 
 conformations .
we suggest using a single nonequilibrium simulation 
 to generate the initial conformations . in the conventional us 
 scheme ,
using weak force constants may result in a higher risk of 
 discontinuous sampling . in the beus simulations , however , the discontinuity 
 in the sampling is less likely and the use of weaker force constants 
 is justified as long as they are stiff enough to ensure a continuous 
 sampling free of any gaps . the convergence , however , is typically 
 slower than that in a stiff - spring setting .
note that using 
 stiff harmonics which result in narrow gaussian 
 distributions is a straightforward approach to simplify the optimization 
 of us / beus parameters .
. the narrower the distribution , 
 the larger the number of windows needed for reasonable exchange rate 
 and overlap . on the other hand
, a broader distribution may converge 
 slower ; thus , the trade - off between the convergence time and number 
 of windows must be considered in choosing the force constants .
finally , we note that for detecting insufficient sampling and hysteresis 
 along the degrees of freedom orthogonal to the reaction coordinate , 
 one may use techniques such as a pairwise consistency test between 
 the probability distribution of adjacent windows as recently introduced 
 by zhu and hummer . for other techniques 
 used to characterize the sampling efficiency , we refer the reader 
 to refs ( 71 ) , ( 73 ) , ( 74 ) , and ( 108 ) .
conventional us is a 
 time - independent biasing scheme ; thus , the samples generated by each 
 umbrella are assumed to satisfy the detailed balance criterion and 
 have the correct biased boltzmann distribution .
however , the initial 
 conformations are not necessarily selected from the correct distribution . 
 the initial part of us simulations , thus , can not be used to construct 
 the biased ( and subsequently unbiased ) distributions and must be discarded . 

note that the equilibration process may be longer for weaker biasing 
 potentials since a broader configuration space is accessible .
in general , 
 the time scale associated with this conformational relaxation is not 
 known a priori and must be determined a posteriori .
for instance , one may discard the samples from an initial period 
 several times longer than the autocorrelation time .
alternatively , 
 one may take advantage of the fact that the sampling during the equilibration 
 period is often inconsistent with the rest of the data and discard 
 ( from the beginning ) as many samples as possible in order to minimize 
 the posterior error . in the beus scheme
, the conformations sampled 
 in a particular umbrella window may belong to different replicas . 
 unlike the extended ensemble ( i.e. , all replicas / windows combined ) , 
 neither individual replica trajectories nor individual ( reconstructed ) 
 window trajectories satisfy the detailed balance criterion .
however , 
 conformations sampled by any given umbrella have the correct biased 
 boltzmann distribution , assuming the extended ensemble is equilibrated . 

this implies that in addition to conformational relaxation , another 
 relaxation associated with the diffusion of replicas within the extended 
 ensemble scheme is required .
similar to the conformational relaxation , 
 the replica / window relaxation time scale may be determined a posteriori .
note that one may equilibrate the initial 
 conformations ( which may come from nonequilibrium simulations ) in 
 a conventional us setting prior to performing the beus simulations . 
 in this case , the relaxation of the replicas in the window space is 
 still required .
an estimate for the relaxation time associated 
 with mixing of the 
 replicas can be obtained from the exchange rates . a lower bound for the relaxation time can be estimated 
 from 2  t/(1  

2 ) in which t is the effective 
 time between exchange attempts and 2 is the second 
 largest eigenvalue of an empirical transition matrix constructed using 
 the exchange rates between the windows . by considering an ideal exchange rate between the windows ,
for instance , in a 1d 
 beus simulation including n windows / replicas , if 
 the exchange is only attempted between the neighboring windows , one 
 can easily show a flat exchange rate of r results 
 in 2  ( n + 
 1)/(10r)t .
we note that we 
 have assumed a stochastic choice of odd  even pairs while the 
 deterministic choice ( as was used in our example ) is expected to result 
 in faster relaxation .
let us examine 
 our beus md trajectories of msba transporter used 
 for sampling along  ( see section 2 ) 
 to estimate the relaxation and correlation times .
note that the two 
 are not the same but related ; the former determines the portion of 
 data to be discarded due to nonequilibrium relaxation effects , while 
 the latter is needed for error estimation and convergence assessment . by constructing the empirical transition matrix , 2
was estimated to be about 184 ps , which is close to the estimate 
 from 2  ( n + 
 1)/(10r)t  180 ps 
 ps .
we also explicitly estimated the autocorrelation times associated 
 with different quantities including the window index , select collective 
 variables , and several principal components ( table 2 ) .
we estimated the autocorrelation times from the correlation 
 functions using the scheme discussed in ref ( 73 ) .
for a given window , the 
 exchange of the conformations in the beus scheme results in apparent 
 decorrelation of the data as compared to conventional us .
the autocorrelation 
 time for different quantities was estimated to be around 50 to 500 
 ps .
however , one may notice that the actual trajectories associated 
 with the replicas have much longer autocorrelation times when compared 
 to reconstructed window - based trajectories ( see table 2 ) .
the correlations were estimated 
 along ( i ) an actual trajectory associated with a given replica and 
 ( ii ) a reconstructed trajectory associated with a given window .
the mean values and standard deviations 
 are given based on 22 trajectories / windows .
pc1 and pc2 
 may be considered the most significant biased and unbiased principal 
 components , respectively .
pc2 roughly represents the slowest degree 
 of freedom orthogonal to the reaction coordinate used for biasing . 

the largest autocorrelation time , which is incidentally associated 
 with pc2 , may be used to estimate the statistical inefficiency of 
 the entire data set as g  ( 2ac * )  3.4 ns .
alternatively , g can be estimated using the so - called structural decorrelation 
 time which is related to the block - averaging 
 approach . defining a metric in the ( pc1,pc2 ) 

space , we estimated dec  3.2 and 0.5 ns 
 for replica- and window - based trajectories , respectively
. with correlation 
 times on the order of nanoseconds , reaching a verifiable convergence 
 may require simulation times on the order of microseconds .
the presence 
 of multiple copies of the system and the mixing of the replicas , however , 
 may speed up the convergence considerably .
nonetheless , given these 
 highly correlated data sets , typical of conformational sampling simulations , 
 it is important to use an appropriate technique for error analysis , 
 e.g. , bayesian block bootstrapping as 
 discussed in section s1 in the supporting information .
in summary , 
 here we introduce a computational recipe for efficient 
 description of large - scale conformational transitions in complex molecular 
 systems .
the approach combines a set of ( i ) driven nonequilibrium 
 and ( ii ) biased equilibrium simulations with similar ( but not identical ) 
 ( i ) time - dependent and
the biasing protocols are designed using a knowledge - based approach 
 taking advantage of system - specific collective variables and empirically 
 optimized parameters .
the process is knowledge - based in that any structural 
 information available on the end states is used both in the design 
 of the biasing protocols and in their quality assessment .
the process 
 is also empirical in that both ( i ) nonequilibrium and ( ii ) equilibrium 
 biasing protocols are optimized iteratively by varying the protocol 
 parameters such as collective variables or centers of umbrella windows 
 and assessing the results , e.g. , based on nonequilibrium work or rate 
 of exchange . the approach discussed here can be viewed as an ad - hoc dimensionality reduction method in that the sampling 
 is adjusted empirically to restrain the system in a particular region 
 of the configuration space representing relevant transition tube(s ) 
 and intrinsic manifold .
several ideas / techniques described here 
 have been previously discussed / employed 
 by us and others elsewhere , some of which are referenced throughout 
 the paper .
the emphasis here is on presenting a coherent picture of 
 the entire process of performing , analyzing , and improving md simulations 
 aimed at describing functionally relevant , large - scale conformational 
 transitions in complex , realistic biomolecular systems .
the process 
 is indeed challenging with diversified , multifaceted issues involved , 
 some of which less appreciated than others .
we propose a balanced 
 strategy to combine several techniques within an empirical framework , 
 to address the practical obstacles involved in sampling large - scale 
 conformational changes .
moreover , the current work calls for a more 
 concerted attempt to develop strategies aimed at improving the applicability 
 of advanced sampling techniques to realistic biological problems .
 characterizing 
 large - scale structural transitions in biomolecular 
 systems poses major technical challenges to both experimental and 
 computational approaches . on the computational side , 
 efficient sampling 
 of the configuration space along the transition pathway remains the 
 most daunting challenge . 
 recognizing this issue , we introduce a knowledge - based 
 computational approach toward describing large - scale conformational 
 transitions using ( i ) nonequilibrium , driven simulations combined 
 with work measurements and ( ii ) free energy calculations using empirically 
 optimized biasing protocols . the first part is based on designing 
 mechanistically relevant 
 , system - specific reaction coordinates whose 
 usefulness and applicability in inducing the transition of interest 
 are examined using knowledge - based , qualitative assessments along 
 with nonequilirbrium work measurements which provide an empirical 
 framework for optimizing the biasing protocol . 
 the second part employs 
 the optimized biasing protocol resulting from the first part to initiate 
 free energy calculations and characterize the transition quantitatively . 
 using a biasing protocol fine - tuned to a particular transition not 
 only improves the accuracy of the resulting free energies but also 
 speeds up the convergence . 
 the efficiency of the sampling will be 
 assessed by employing dimensionality reduction techniques to help 
 detect possible flaws and provide potential improvements in the design 
 of the biasing protocol 
 . structural transition of a membrane transporter 
 will be used as an example to illustrate the workings of the proposed 
 approach . 
Introduction
Combining Nonequilibrium
and Equilibrium Approaches
Technical Considerations
Conclusion