molecular 
 dynamics ( md ) simulations of detailed atomic models is 
 a powerful tool to study the properties of complex biomolecular systems . however , while simulations based on realistic all - atom ( aa ) models 
 arguably offer the most detailed information , such models evolve on 
 a complex and rugged energy surface , and their dynamics is often burdened 
 by a host of slow processes . for this reason , achieving an adequate 
 sampling of all the relevant configurations of a system from straight 
 md simulations is often challenging . most proposed approaches to enhance 
 sampling by accelerating the exploration of configurational space 
 are either not guarantied to yield boltzmann equilibrium or only applicable 
 to a small subset of degrees of freedom ; see refs ( 1417 ) for reviews . the metropolis monte 
 carlo ( mc ) algorithm , which consists of generating a random walk in 
 configuration space from a set of proposed moves that are attempted 
 and then accepted or rejected , also offers a powerful 
 method to generate a boltzmann equilibrium distribution that is not , 
 in principle , limited by dynamical processes . in practice , its application 
 is only limited by the richness of the set of proposed moves that 
 are attempted for generating a random walk in configurational space . 
 nevertheless , attempts to sample large conformational changes with 
 mc remain completely ineffective for complex molecular systems in 
 the presence of explicit solvent due to the rejection of all proposed 
 new configurations . an appealing alternative approach is to 
 consider a simplified coarse - grained 
 ( cg ) model that relies on an effective many - body potential of mean 
 force ( pmf ) associated with a reduced set of degrees of freedom . because the many - body pmf governing a cg model is expected to be 
 generally smoother than the energy surface of its aa counterpart , 
 its dynamics is often intrinsically simpler and faster . however , except 
 for the simplest situations , the effective many - body pmf associated 
 with the cg degrees of freedom is not known exactly . commonly , an 
 approximate pmf is constructed , which affects many thermodynamic and 
 kinetic results . force - matching , whereby information extracted from 
 the fine - grained ( fg ) model via the mapping function is transferred 
 toward the cg model , provides a formal route to improve the many - body 
 pmf underlying a cg model . however , the accuracy 
 of the resulting cg model remains inherently limited by the functional 
 form chosen to represent the many - body pmf . even assuming that 
 a long trajectory of a cg model has been generated , 
 exploiting this information to enhance the configurational sampling 
 of the fg model is not straightforward . in principle , a set of configurations 
 for the fg model could be generated by reconstructing , via a operation , the missing degrees of freedom 
 onto snapshots extracted from a cg trajectory . unfortunately , rcg 
 generally involves mostly ad hoc empirical modeling 
 procedures that do not yield a boltzmann equilibrium distribution . 
 in this sense , more rigorous approaches allowing the information to flow back 
 and forth between a fg model and its associated cg model have sought 
 to couple the two levels via a resolution - exchange strategy . in particular , resolution - exchange with incremental coarsening relies 
 on a potential energy that interpolates between the fg and cg models 
 according to a parameter . the hope of resolution - exchange is to combine the efficiency of 
 cg simulation and the accuracy of a fg model by swapping configurations 
 between the two representations , although achieving considerable gains 
 in computational efficiency remains a challenge . one proposed route 
 to improve the efficiency of the algorithm has been to first relax 
 configurations before an exchange is attempted , but this procedure 
 yields the correct canonical sampling only if the rate of exchanges 
 is sufficiently low . an alternative avenue 
 is the multiscale enhanced sampling ( mses ) method of zhang and chen , 
 in which a small set of cg auxiliary particles evolving on a simplified 
 potential energy surface is coupled to an all - atom system via harmonic 
 springs of variable strength within a replica - exchange simulation 
 framework . despite these previous 
 efforts , a robust framework that exploits 
 the information from a cg model to generate the proper boltzmann equilibrium 
 distribution for a fg model is presented by the recent extensions to the 
 mc algorithm that combines the strength of nonequilibrium molecular 
 dynamics ( nemd ) and mc . in hybrid nemd - mc , the value of some chosen variable is altered 
 gradually in a time - dependent controlled manner , while the remaining 
 degrees of freedom are allowed to evolve freely according to the dynamical 
 equations of motion . the configuration generated by the nonequilibrium 
 switching process is then treated as a candidate that must be either 
 accepted or rejected via a metropolis criterion to generate the equilibrium 
 boltzmann distribution . this hybrid nemd - mc framework has most notably 
 been used to formulate a constant - ph simulation algorithm . here , this idea is pursued 
 further to design a novel multiscale 
 method that exploits the evolution of a cg model to help generate 
 target candidate configurations that are then used to guide the fg model during the nemd switching trajectories . monte carlo 
 ( cg - guided hybrid nemd - mc ) algorithm , comprises six steps , illustrated 
 schematically in figure 1 : ( 1 ) a fg configuration of an atomic system is dynamically propagated 
 for some period of time using equilibrium md ; ( 2 ) the resulting fg 
 configuration is mapped onto a simplified cg model ; ( 3 ) the cg model 
 is propagated for a brief time interval to yield a new cg configuration ; 
 ( 4 ) the resulting cg configuration is used as a target to guide the 
 evolution of fg system ; ( 5 ) the fg configuration ( from step 1 ) is 
 driven via a nonequilibrium md ( nemd ) simulation toward the cg target ; 
 ( 6 ) the resulting fg configuration at the end of the nemd trajectory 
 is then accepted or rejected according to a metropolis criterion before 
 returning to step 1 . a symmetric two - ends momentum reversal prescription 
 is used for the nemd trajectories of the fg system to guarantee that 
 the cg - guided hybrid nemd - mc algorithm obeys microscopic detailed 
 balance and rigorously yields the equilibrium boltzmann distribution . it is shown that the cg - guided hybrid nemd - mc algorithm can be carefully 
 engineered to achieve a reasonably high acceptance probability , even 
 when using fairly short nemd switching trajectories . more importantly , 
 because the transition probabilities are constructed to satisfy detailed 
 balance , the cg - guided hybrid nemd - mc is guaranteed to yield the equilibrium 
 boltzmann distribution . in the next section , the performance of the method is 
 illustrated with applications to simple model systems and solvated 
 polypeptides . let the total energy of a system be e(x ) = u(r ) + k(p ) , where u is the potential energy , k is the kinetic energy , and x represents all 
 of the degrees of freedom , including the coordinates r and the momenta p. in metropolis monte carlo , we seek 
 to generate a stochastic random walk process that moves the system 
 from a configuration x to a configuration x and ensure that the random walk will obey detailed balance 
 for the boltzmann distribution:1where (x ) = q exp[e(x ) ] is the equilibrium probability (  = 1/kbt ) , and   is the 
 transition probability . one common 
 approach to construct such a random walk is to separate the transition 
 probability into two stages : ( 1 ) the probability to generate a proposed 
 move and ( 2 ) the probability to accept ( or reject ) this move,2which leads 
 to the condition3the probability of a proposed move 
 describes 
 the probability of reaction coordinates moving toward the target value . the probability to accept or reject the proposed move , however , is 
 determined after the switch . the most common approach is when the 
 transition probability of the proposed move is inherently symmetric , 
 i.e. ,   ( x  x ) = ( x  x ) , leading to the 
 condition for the acceptance probability4that is formally satisfied by the 
 metropolis 
 criterion,5examples include constant - ph simulation or 
 biminima simulation , where the reaction coordinates of the nemd - mc 
 can only choose predetermined values . for example , in a hybrid md / nemd - mc 
 constant - ph simulation , the reaction coordinate for nemd - mc is the 
 protonation state . therefore , the target value is always the deprotonated 
 state when the initial value is the protonated state and vice versa . 
 as a result ,   is always 1 . 
 in the biminima simulation , the reaction coordinates are usually shifted for a fixed amount , and 
 therefore   is also 1 . by 
 substitution , it can be readily 
 shown that6and therefore the metropolis construct satisfies eq 1 . in this case , if the 
 constraint schedule is time reversible and the momentum is properly 
 treated , the nemd - mc will generate the correct boltzmann distribution . if the long time scale 
 relaxation of the system is known to be 
 dominated by the dynamics along a set of coarse - grained ( cg ) variables , r , it is possible to exploit this information to increase 
 the efficiency of the hybrid nemd - mc algorithm . first , it is assumed 
 that at any time , the cg variables r(t ) are uniquely defined from x(t ) through 
 a mapping function as , r = m[x ] . the general idea is to separate the transition probability ( x  x ) 
 into two distinct steps,7where   represents the transition probability for 
 a set of cg variables r , and latter   is the 
 transition probability for the fg 
 variables x , conditional on the transition r  r taking place . here , the 
 first step only involves changes in the cg variables . the transition 
 probability , , is constructed such that it obeys the 
 cg detailed balance relationship,8where (r ) 
 is the equilibrium probability of the r coordinates associated 
 with the cg model . assuming that the cg model is constructed on the 
 basis of a potential of mean force w(r ) , the probability ratio is,9 in practice , a 
 number of methods could be 
 used to fulfill these conditions while propagating the cg coordinates 
 ( e.g. , brownian dynamics , langevin dynamics , metropolis mc , etc . ) . in the second step , the transition r  r in cg space must be used to guide the changes in the remaining 
 atomic fg degrees of freedom . to formalize this idea , it is useful 
 to rewrite the transition probability   as the 
 product of the probability of a proposed move , and the probability 
 to accept 
 or reject the proposed move , 10 in principle , 
 generating the transition probability   for the proposed move in x , conditional on the transition r  r , could be a daunting task . even simply mapping the set of 
 fg variables x that is consistent with the cg variables r , a problem that is commonly referred to as however , 
 all issues with such a reverse mapping problem are resolved naturally 
 and rigorously if the system s degrees of freedom x(t ) are generated from a nemd propagation under 
 the influence of a time - dependent constraint dragging the cg variables 
 from their initial r to the final r 
 value . during the nemd dynamical r  r switching process , the evolution of the cg coordinates r(t ) is externally controlled and follows 
 a fixed time - dependent schedule over a time interval nemd , while the remaining degrees of freedom are propagated freely11the 
 cg coordinates start in the initial state r at the beginning 
 of the nemd switch and are gradually altered 
 in a time - dependent manner to reach the final state r at a time interval nemd later . this process 
 can be implemented by evolving the system under a time - dependent holonomic 
 constraint , or by applying stiff harmonic restraints centered on r(t ) . finally , the cg - guided scheme must 
 satisfy the global detailed balance relationship,12to guarantee detailed balance , the time - dependent r(t ) for the forward and backward switching , r  r and r 
 r , must be consistent . this is automatically 
 satisfied if r(t ) is constructed from 
 the linear form eq 11 . a nonlinear form could also be used for the switching , as long 
 as it is symmetric and time - reversible . assuming that the transition 
 probability of the proposed move is symmetric,13which is verified if the dynamical 
 propagation 
 used to generate the nemd trajectory is deterministic and reversible 
 ( e.g. , using a symplectic integrator ) , we obtain,14it follows that in the 
 cg - guided hybrid nemd - mc 
 scheme , the metropolis acceptance probability , , is15 in the present development , it is assumed 
 that the cg model is propagated via a thermalized dynamics satisfying 
 the detailed balance condition eq 8 . this is the reason why the energy difference between 
 the cg configurations ( w(r ) w(r ) ) appears in the acceptance criterion eq 15 . alternatively , one 
 could carry out the dynamics of the cg system with a propagator that 
 conserves energy , consistent with a microcanonical ensemble . in this 
 case , the energy difference of cg configurations would not appear 
 in the metropolis criterion . it is also worth pointing out that the 
 cg coordinates r correspond to a subspace of the fg coordinate x that is generated via the mapping function m[x ] through eq 11 . in this sense , the acceptance criterion operates only in 
 the fg space and this is the reason why the acceptance criterion does 
 not involve a joint distribution in terms of ( x , r ) . as illustrated in figure 1 , the cg - guided hybrid nemd - mc algorithm 
 comprises the following 
 steps : ( 1 ) propagate the fg system with equilibrium md using the potential 
 energy u(x ) for a period of time md ; ( 2 ) extract the cg coordinates r from x via the mapping function m ; ( 3 ) propagate the 
 simple cg on the free energy surface w(r ) ( for a predetermined time cg or until a chosen stopping criterion is met ) to yield the new cg 
 configuration r ; ( 4 ) save the final cg coordinates r to use as a target for the nemd trajectory ; ( 5 ) 
 propagate the fg model from x to x 
 under the time - dependent constraint that m[x ] changes linearly from r to r 
 over an interval nemd ; and ( 6 ) 
 accept or reject the final configuration of the fg system according 
 to the metropolis probability eq 15 before returning to step 1 . by construction , the cg - guided hybrid algorithm yields the correct 
 boltzmann equilibrium distribution for the fg model , regardless of 
 the underlying cg model that is chosen to generate the attempted moves . nonetheless , the choice of an optimal cg model to achieve the highest 
 efficiency is an important issue that needs further consideration . generally , important sampling techniques aim at achieving variance 
 reduction in computer simulations by sampling high and low probability 
 regions with equal frequency while recovering a correct distribution 
 by assigning a proper statistical weight to the different regions . the cg - guided hybrid algorithm can be designed according to the same 
 guiding principles by recognizing that the potential energy driving 
 the cg model actually governs the statistical weight attributed to 
 the coordinates r. for a given fg model , the  exact  
 pmf with respect to the cg degrees of freedom is defined as,16where c is some constant . if w(r ) were known , 
 then using it in the algorithm would result in an optimal cg model . this is because the cg - guided nemd simulation would sample all regions 
 with equal probability , as demonstrated by the following development:17this 
 relationship is valid for any r and r 
 within the cg space . in addition , the 
 average acceptance probability for nemd transitions between r and r , which also affects the simulation 
 efficiency , is18when increasing w(r ) relative to w(r ) , (x  x)r increases while (x  x)r decreases and vice versa ; eq 18 reaches a maximum if the average acceptance 
 probability is roughly equal , when using w(r ) for the cg model . in practice , the exact pmf with 
 respect to the cg degrees of freedom is not known , but the above argument 
 shows that , in an importance sampling sense , the optimal efficiency 
 will be achieved if one chooses a reasonable approximation to w(r ) that accurately captures 
 the dominant effects . a one - dimensional 
 linear chain system with 12 particles linked together 
 in a periodic potential was designed to test and illustrate the enhanced 
 sampling from the cg - guided hybrid nemd - mc . the fg system comprises 
 12 particles , each with a mass of 1 , and each pair of adjacent particles 
 is connected by a bond . the only two types of forces acting on each 
 particles are the bond force and the potential force . the potential 
 force , shown in figure 3 , is defined by u(i ) = cos(2 xi ) cos(10 xi)/2 , where x is the coordinate , and i is the index 
 of the particle ( going from 1 to 12 ) . the first part of u(i ) defines the periodicity . the bond force is defined 
 by ebond(i , j ) = 20(rij  
 1)(|i  j| 1 ) , where rij is the distance of the two particles , and the function  indicates 
 that only adjacent particles are connected by a bond . the periodic 
 potential energy divides the configuration space into infinite number 
 of wells along the x axis , with the width of each 
 well being 1 and the energy barrier between adjacent wells being around 
 2.6 . since the 
 optimal bond length is also 1 , the most stable configuration of the 
 molecule is to occupy 12 adjacent wells , each with one particle . to 
 diffuse along the x axis , a couple of particles , 
 if not all , have to cross the energy barrier at the same time . therefore , 
 the diffusion constant from md simulations is extremely low ( see figure 6a ) . the cg system 
 contains only three effective particles , each with a mass of 4 . the bond force 
 for the cg system is given by ebondcg(i , j ) = 20 ( rij  4 ) (|i  j|  1 ) . the cg particles 1 , 2 , and 3 correspond 
 to the centers of mass of fg particles 14 , 58 , and 
 912 , respectively . the simulation of this system is carried 
 out using an in - house code written in c++ . the temperature kbt is set 
 to 1 ; diffusion constant to 0.5 ; time step to 0.005 ; and collision 
 rate to 2.0 . the lengths of equilibrium md , cg , and nemd simulations 
 are all 1000 steps . flowchart of 
 cg - guided hybrid nemd - mc for a 12 particle model system . big ( small ) spheres represent the particles in the cg ( fg ) system . configuration of the cg ( fg ) system is represented 
 using r ( x ) . the solid lines between the fg and cg system represent 
 the center - of - mass constraints applied on the fg system using the 
 position of the cg particles . the arrows between the fg and cg system 
 represent the mapping function m. ( a ) from the initial 
 fg structure x , the cg structure r is built . ( b ) dynamical propagation is performed on the cg model , generating 
 a new configuration r. ( c ) dynamical propagation 
 is performed on the fg model , generating a new configuration x. during this propagation , the position of the cg 
 model is used to constrain the center - of - mass for the fg system . its analytic form is u = cos(2x ) cos(10x)/2 . the potential energy 
 is periodic and extends to infinity . the algorithm 
 was tested on four different aa peptide systems with explicit solvent : 
 trialanine(ala3 ) , penta - alanine(ala5 ) , deca - alanine(ala10 ) , and ( aaqaa)3 . ( aaqaa)3 is a 15 residue 
 peptide , with residues 3 , 8 , and 13 as glutamine and the rest as alanine . , two different initial structures were prepared , 
 one -sheet like and another -helix like , as shown in figure 4 . ala3 and ( aaqaa)3 have an acetylated n - terminus and an n - methylamide c - terminus . ala5 and ala10 have unblocked , charged termini . the potential energy of the solvated 
 peptide systems is represented by the charmm36 force field and tip3 water potential . ala3 is solvated in a 24  cubic box : ala5 , 30  ; ala10 , 40  ; and ( aaqaa)3 , 60  . the solvent for the ala5 and ala10 systems is a 1 m kcl aqueous solution . the cg systems for 
 polyalanines and ( aaqaa)3 systems include only a single 
 atom type for the backbone  carbon ( cgc ) . unless mentioned otherwise , the bond , 
 angle and dihedral energy terms for cgc are empirically set as ebond = 100 ( r  3.83 ) , eangle = 94.84 (   
 100.0 ) , and edihedral = 0.1 
 (   55.0 ) , where r ,  , 
 and  are the bond length , angle , and dihedral , respectively . there is no water and no periodic boundary conditions ( pbc ) for the 
 cg system . the red one has an -helix 
 for the first 10 residues and -sheet for the rest of the 5 
 residues . the charmm program version c36b1 is 
 used for md simulations for both the aa and cg system . in - house python 
 scripts are used for the general control flow of the algorithm ( generating 
 charmm input files , controlling mc acceptance , etc . ) . input files 
 are generated at each attempted move with the newest parameters and 
 constraints . unbiased brute - force md simulations were also performed 
 with the same systems to provide a comparison . for the aa simulations 
 with explicit solvent ( md and nemd ) , particle - mesh ewald ( pme ) summation is used to treat the electrostatic interactions , 
 with a real - space cutoff set to 14  and grid spacing smaller 
 than 0.5  . the lennard - jones ( lj ) interactions are smoothly 
 truncated with a switching function from 10 to 12  . the equations 
 of motion are integrated with a time - step of 2 fs , and shake is used to constrain covalent bonds involving 
 hydrogen atoms . for md , the peptide is kept near the center of the 
 pbc box with a weak global center - of - mass restraint . the leapfrog 
 varlet integrator is used with constant temperature and pressure control 
 ( cpt ) based on berendsen s method . temperature is controlled 
 at 300 k and pressure at 1 atm . for nemd , the initial constraint position is mapped from 
 the aa structure , and the final constraint position is generated by 
 cg simulations . theoretically , there could be a small  tracking  
 error between r and m[x ] when 
 harmonic restraints are used instead of holonomic constraints , but 
 the tracking error is expected to be small with strong restraints , 
 and does not affect the accuracy or the foundation of the theory . 
 to propagate the cg system , a langevin dynamics is used with a temperature 
 of 300 k. the overall rotation and transition are removed from the 
 cg configuration . in each cycle , the initial velocities of the cg 
 model were generated according to the maxwell distribution . in the present implementation , the md and nemd trajectories have 
 a fixed length while the propagation of the cg coordinates is controlled 
 by a  stopping criterion  . choosing carefully the stopping 
 criterion to control the magnitude of the displacement r  r can greatly increase the efficiency 
 of the method ( see the discussion section ) . 
 here , two different stopping criteria were tested . the first criterion 
 is based on the root - mean - square - deviation ( rmsd ) of the set of cg 
 coordinates ; the cg simulation is stopped when the difference r the second criterion is based 
 on the deviation of angles within the set of cg coordinates and works 
 also with a preset maximum allowed value , . taking the ( aaqaa)3 system as an example , the cg system 
 comprises 15 atoms and therefore 13 angles formed by adjacent atoms . the initial angles are recorded , the absolute change in degrees ( between 
 180 to 180 ) for each is monitored , and the cg simulation is 
 stopped when the average absolute change exceeds the preset maximum 
 allowed value . for ala3 , ala5 , ala10 , and ( aaqaa)3 , both were tested . lastly , 
 to prevent sampling of the cis - peptide bond conformation , 
 all target configurations generated by cg simulation with bond lengths 
 shorter than 3.5  were discarded . the general validity of the cg - guided hybrid 
 nemd - mc was first 
 ascertained using a simple model of a linear chain of linked particles 
 in a one - dimensional periodic potential ( figure 2 ) . in figure 5 , the distribution of each single particle , the center - of - mass 
 of each particle group , and of the entire linear chain are plotted . 
 for all degrees of freedom , the cg - guided hybrid nemd - mc and the equilibrium 
 md generate the exact same distributions . the global diffusion constant 
 of the linked chain is a good indicator of the overall sampling efficiency 
 of the simulation . for each scheme , the evolution of the center - of - mass 
 for the entire chain is shown in figure 6 . the diffusion constant 
 is 0.53 per 1000 rounds for equilibrium md and 43 for cg - guided hybrid 
 nemd - mc . according to this estimator , the cg - guided hybrid nemd - mc 
 gives a speedup of 80 times for this system . while the cg - guided hybrid 
 nemd - mc algorithm can rigorously yield the proper boltzmann equilibrium 
 distributions , it is absolutely critical to satisfy all conditions 
 of microscopic reversibility for a valid algorithm . for example , using 
 a time - dependent constraint that is not symmetric with respect to r and r for the switching schedule fails 
 to reproduce the correct equilibrium distributions of each particle 
 and particle group ( results not shown ) . black lines show the distribution simulated using equilibrium 
 md simulation ; red lines show the cg - guided hybrid nemd - mc . the subplot 
 ( 112 ) presents the distribution for each atom ; ( 1315 ) 
 for the center - of - mass of each atom group ; and ( 16 ) for the center - of - mass 
 of the entire molecule . this figure presents the relative position of the atom or center - of - mass 
 inside each well , no matter which well it is residing in . evolution of the center - of - mass 
 of the entire linked chain molecule . each round contains md , cg simulation , and nemd - mc , each for 1000 
 steps . the linear regression generates a coefficient of determination r value above 0.99 for either scheme . the overall performance of the 
 cg - guided hybrid nemd - mc algorithm 
 was then examined for realistic atomic models of biomolecules with 
 explicit solvent . those include tri- , penta- , and deca - alanine peptides 
 solvated in water . to ascertain the formal correctness of the method , 
 the population of the various conformations was calculated . the results 
 of the simulation are summarized in table 1 . according to this analysis , cg - guided hybrid 
 nemd - mc and equilibrium md generate the same distributions . for a 
 valid comparison , however , it is important that the equilibrium md 
 simulations be sufficiently long to reach proper convergence . the 
 accuracy of the method was further examined for ala3 using 
 a set of simulation parameters for the cg - guided hybrid nemd - mc simulations 
 ( md = 2ps , nemd = 4ps , and rmsd = 1.5 ) . the histogram 
 of distances and angles of adjacent carbon- is shown in figure 7 . the pmf along  
 or  is shown in figure 8 , and the 2d - pmf along  and  angles is shown 
 in figure 9 . according 
 to this analysis , the histograms and pmfs generated by cg - guided hybrid 
 nemd - mc are very similar to those obtained from long equilibrium md 
 simulations . the similarity of the pmfs along  or  angles 
 with previous results from mu et al . provides 
 additional confirmation of the validity of the cg - guided hybrid nemd - mc . using rmsd for 
 ala3 and  for others . this simulation of ( aaqaa)3 had an initial structure with an -helix ( figure 4 , right ) . the following row 
 presents result for cg - guided hybrid nemd - mc of ( aaqaa)3 with an initial structure with a -sheet ( figure 4 , left ) . is defined as 
 160 
 <  < 20 , 120 <  < 50 ( figure 9b , blue box ) ;  
 180 <  < 90 , 50 <  < 180 , 
 or 180 <  < 90 , 180 <  
 < 120 , or 160 <  < 180 , 110 <  < 
 180 ( figure 9b , cyan 
 box ) ; ppii 90 <  < 20 , 50 <  
 < 180 , or 90 < < 20 , 180 
 <  < 120 ( figure 9b , magenta box ) ; l- 0 <  < 120 , 
 30 <  < 90 ( figure 9b , yellow box ) . for ala3 and ( aaqaa)3 , termini are blocked , and the averages of all residues are 
 shown ; for ala5 and ala10 , averages of all but 
 two terminal residues are shown . trajectories are too short to generate 
 useful statistics . for ala3 and ala5 , only transitions of the middle residue 
 are counted . for 
 ala10 , the average of transitions of residues 38 
 is counted . for ( the black line presents the results from equilibrium 
 md ; the red line from cg - guided hybrid nemd - mc ; and the green line 
 the theoretical distribution calculated from cg parameters . subplot a presents 
 a histogram of the angle of ca1-ca2-ca3 ; b the distance between ca1 
 and ca2 ; and c the distance between ca2 and ca3 . the black 
 line presents the results from equilibrium md and the red line from 
 cg - guided hybrid nemd - mc . subplots a and d present pmf along the  
 and  angles of residue 1 ; b and e of residue 2 ; and c and f 
 of residue 3 . c present the pmf calculated from 
 the simulation results of equilibrium md ; d  f of cg - guided 
 hybrid nemd - mc . subplots a and d present the pmf of residue 1 ; b and 
 e of residue 2 ; and c and f of residue 3 . the color key is shown in 
 a. definitions of  ( blue box ) ,  ( cyan box ) , ppii ( magenta 
 box ) , and l- ( yellow box ) are shown in b. transitions between the  ,  , ppii , and l- 
 backbone 
 conformers are expected to be representative of the slowest motions 
 of polyalanine peptides in solvent . to quantify the kinetic acceleration 
 gained by the cg - guided hybrid nemd - mc algorithm , we define the speedup 
 factor ,  , as the ratio of the number of transition events from 
 the cg - guided hybrid nemd - mc simulation , relative to that from the 
 equilibrium md simulation . for the sake of the comparison , the total 
 simulation length ascribed to the cg - guided hybrid nemd - mc algorithm 
 comprises the length of both the equilibrium md and all of the attempted 
 nemd switches ( whether they are accepted or rejected ) . to facilitate 
 the discussion , it is useful to distinguish  global transitions  , 
 corresponding to interconversions between the  , l- , 
  , and ppii conformers ( not including interconversions between 
  and ppii ) , and transitions 
 from  to  are opposed by a large energy barrier . however , 
 the energy barrier for transitions from ppii to  is relatively 
 small . accordingly , we determined the global and local speedup factors , 
 global and local , from the simulation 
 data . on the basis of this analysis , the sampling efficiency appears 
 to vary with the length of md , md , length of nemd , 
 nemd , and the maximum allowed displacements  or rmsd . in almost all of the cases 
 examined here , cg - guided hybrid nemd - mc has higher efficiency than 
 straight equilibrium md , going up to 15-fold acceleration . according 
 to eq 17 , the accuracy 
 of the underlying cg model affects the overall efficiency 
 of the cg - guided hybrid nemd - mc algorithm . specifically , optimal efficiency 
 is achieved if the effective energy surface of the cg model is a reasonably 
 good approximation to the exact pmf with respect to those degrees 
 of freedom . to illustrate this important point , we examined the impact 
 of an improved cg model on the sampling efficiency of the ala5 peptide . the energy terms of the improved cg model were set 
 to ebond = 400 ( r  
 3.83 ) , eangle = 30 (  
 115.0 ) , and edihedral = 0.1 (   55.0 ) , in order to better match 
 the distribution from unbiased md simulation ( figure 10 ) . the other parameters were the same as 
 that in the first row for ala5 in table 1 . under these conditions , the average acceptance 
 ratio increases from 35% to 45% , demonstrating the additional gain 
 in sampling efficiency with an improved cg model . the black solid line presents the distribution defined by ebond = 100 ( r  3.83 ) , eangle = 94.84 (   
 100.0 ) , and edihedral = 0.1 
 (   55.0 ) . the black dotted line presents 
 the distribution defined by ebond = 400 
 ( r  3.83 ) , eangle = 30 (   115.0 ) , and edihedral = 0.1 (   55.0 ) . subplot a presents the angle between the carbon- of residues 
 2 , 3 , and 4 . finally , the accuracy and sampling efficiency of the cg - guided 
 hybrid nemd - mc algorithm was examined for the solvated ( aaqaa)3 peptide ( figure 4 ) . compared with equilibrium md , an estimated 
 speedup on the order of 21 times is achieved with the cg - guided hybrid 
 nemd - mc . the evolution of  or  angles for residues 1 , 
 4 , and 8 is shown in figure 11 . for residue 1 , the number of total transitions are comparable 
 for equilibrium md simulation and cg - guided hybrid nemd - mc . for residues 
 4 and 8 , no global or local transitions were observed using equilibrium 
 md , while many were observed using cg - guided hybrid nemd - mc . the reason 
 is that the -helix of the initial structure did not unfold 
 during equilibrium md simulations but unfolded and refolded with cg - guided 
 hybrid nemd - mc . the evolution of average population of  ,  , 
 ppii , l- , and helix is shown in figure 12 . the cg - guided hybrid nemd - mc simulations 
 starting from different conformations appear to quickly converge toward 
 unique results . the black line presents the results from 
 equilibrium md and the 
 red line from cg - guided hybrid nemd - mc . both simulations had an initial 
 structure with an -helix ( figure 4 , right ) . subplots a and d present the evolution 
 of  and  angles of residues 1 ; b and e of residue 4 ; 
 and c and f of residue 8 . the solid 
 lines present the results from a cg - guided hybrid nemd - mc ; the dashed 
 lines are from equilibrium md . the red and 
 blue lines present simulations where the initial structure has an 
 -helix ( figure 4 , right ) ; magenta  sheet ( figure 4 , left ) . subplots a d present the 
 percentage of  ,  , ppii , and l-. the average percentage 
 of residues 610 is shown . the structure contains an -helix secondary 
 structure when three successive residues are in  conformation . one process observed in the cg - guided 
 hybrid nemd - mc simulations 
 of the solvated peptides is the cis  trans isomerization 
 of the peptide linkages . this was somewhat unexpected because this 
 process , opposed by a high energy barrier , is commonly never observed 
 in equilibrium md simulations at room temperature . in retrospect , 
 however , it is clear that the occurrence of some cis configuration is entirely consistent with the boltzmann equilibrium 
 distribution of the polypeptide system . in this sense , long equilibrium 
 md simulations only explore a fraction of all possible configurational 
 space ; they generate the equilibrium distribution of a polypeptide 
 that is kinetically trapped in the trans state . in 
 contrast , cg - guided hybrid nemd - mc simulations aim at achieving a 
 full configurational sampling of the system and thus are not limited 
 by slow kinetic processes that are opposed by large energy barriers . this observation illustrates vividly the power unleashed by a simulation 
 method that truly seeks to enhance configurational sampling . however , 
 there are certainly circumstances where one would realistically want 
 to restrict the conformational sampling to the subspace that is explored 
 by equilibrium md . in particular , in the present simulations we explicitly 
 prevented cis more generally , one could restrict the conformational 
 sampling to the subspace accessible to equilibrium md in various ways , 
 for example , by introducing confinement potentials on different degrees 
 of freedom . the present set of simulations allows us to draw 
 some general principles 
 regarding the main factors affecting the performance of the cg - guided 
 hybrid nemd - mc algorithm . the magnitude of the displacement r  r that should be taken from 
 the cg model to generated proposed mc moves for the aa system is one 
 particularly important factor to consider . generally , one should avoid 
 using a target configuration r corresponding 
 to an exceedingly large conformational transition since the proposed 
 move is unlikely to be accepted during the metropolis step . however , 
 using a target configuration r corresponding 
 to a minuscule conformational transition is clearly unproductive and 
 should also be avoided . the magnitude of the displacement can be controlled 
 by using a fixed propagation time of the cg model ( cg ) or by using a  stopping criterion  based on a maximum 
 allowed displacement . one of the stopping criteria used here was based 
 on a maximum allowed value for rmsd of the cg coordinates , rmsd . for all of the solvated peptide systems examined here , 
 the sampling performance was better with a rmsd - based stopping criterion 
 than with a fixed propagation time ( results not shown ) . the reason 
 is that the change in conformation can vary substantially for cg simulation 
 propagated by a fixed amount of time . however , rmsd - based stopping 
 criteria can also become suboptimal for long peptides like ala10 or ( aaqaa)3 because global transitions for residues 
 near the middle of the chain are not well sampled . the latter do not 
 occur during the brief cg simulations as they correspond to rmsd values 
 that greatly exceed the stopping criterion rmsd . this observation led us to introduce an angle - based stopping criterion 
 with a maximum allowed change in the angle , . the cg - guided hybrid nemd - mc simulations of deca - alanine and ( aaqaa)3 using the angle - dependent stopping criterion outperformed 
 the ones with rmsd - dependent stopping criterion , by at least 2-fold 
 ( results not shown ) . the angle - based stopping criterion succeeds in 
 enhancing the sampling of these transitions because it treats the 
 conformational change within each residue more uniformly than a rmsd - based 
 stopping criterion . consequently , a broad range of motions is accelerated , 
 including large conformational change transition involving residues 
 that are in the middle of the chain . it might be possible to design 
 alternative stopping criterion to further accelerate global transition 
 near the middle of the chain . ultimately , the performance of 
 the cg - guided hybrid nemd - mc algorithm 
 is sensitive to both the magnitude of the cg displacement , r  r , and the length of the nemd switching 
 trajectory , nemd . a large nemd is definitely required to obtain a reasonably high acceptance probability 
 for large proposed mc moves . the simulations reported in table 1 were designed so 
 that a large nemd was always used when a large displacement 
 is allowed . in the final analysis , while the optimal combination of 
 nemd and maximum allowed cg displacement is probably 
 system - dependent , it should be possible to adaptively refine the value 
 of these parameters depending on the type of conformational change 
 that needs to be enhanced . nevertheless , in trying to optimize the 
 performance of the method , it is crucial to make sure that the microscopic 
 reversibility of the cg simulation imposed by eq 8 is rigorously maintained . a typical aa model of a complex bimolecular 
 system often contains 
 a very large number of degrees of freedom . straight unbiased md simulations 
 progress very slowly in this high dimensional space and are often 
 inefficient to adequately sample all of the meaningful configurations . the cg - guided hybrid nemd - mc simulation described here offers a promising 
 departure from equilibrium md . it aims to increase the sampling efficiency 
 by using the evolution of the simpler cg model as a guide to drive 
 any chosen motions with the aa system that are thought to be intrinsically 
 slow . the cg - guided hybrid nemd - mc method produces a boltzmann 
 equilibrium 
 sampling that is rigorously valid , for any reasonable choice of effective 
 energy surface for the cg model . technically , the cg model is used 
 only as a guide to generate proposed mc moves that are then accepted 
 or rejected via a metropolis probability eq 15 . thus , a key feature of the cg - guided hybrid 
 nemd - mc method is that the imperfections and inaccuracies of the cg 
 model do not formally affect the ultimate outcome of the simulation . for example , the parameters of the cg model of polyalanine were not 
 optimized , and the latter is clearly imperfect and inaccurate , as 
 shown from the differences between the distributions with respect 
 to the cg coordinates ( figure 7 ) . notwithstanding these imperfections of the cg model , the 
 cg - guided hybrid nemd - mc nonetheless yields the correct distributions , 
 identical to those obtained from equilibrium md of the aa model . of 
 particular importance , the cg - guided hybrid nemd - mc method does not 
 rely on a multicopy replica - exchange framework , which can become very 
 costly from a computational point of view . a wide range of avenues could be explored to build on the 
 present 
 cg - guided hybrid nemd - mc framework . of paramount importance is the 
 choice of cg degrees of freedom and how they relate to slow motions 
 within the aa system . to investigate a novel aa system with unknown 
 properties , a sound strategy would be to first try to detect the relevant 
 slow modes from the limited information from aa simulation data in 
 order to construct a meaningful cg model , which could then be progressively 
 refined . another practical aspect that also affects the efficiency 
 of the algorithm is the magnitude of the motion from the cg simulation 
 that is utilized to generate proposed mc moves . clearly , no substantial 
 gain in sampling efficiency can be expected if the cg model is used 
 only to generate very small displacements . however , being overly ambitious 
 has also some drawbacks ; the acceptance probability may become vanishingly 
 small if very large displacements of the cg coordinates are used to 
 drive the aa system with nemd simulations . to resolve this issue , 
 one may conceive of an adaptive procedure aimed at maximizing the 
 acceptance probability and the sampling efficacy of the cg - guided 
 hybrid nemd - mc . in practice , the motion of the cg variables may be 
 controlled either via the length of the cg trajectory or be determined 
 on the basis of some maximum displacement criterion . for the solvated 
 ala5 , ala10 , and ( aaqaa)3 systems , 
 the best results were obtained with an angle dependent criterion ; 
 the cg simulations were stopped when the average difference of the 
 angle was larger than a certain cutoff . this was chosen such that 
 all residues were treated equivalently in terms of transition size . ultimately , while the cg - guided hybrid nemd - mc is rigorously valid 
 for any cg model , the overall efficiency will be better if the effective 
 energy surface of the cg model is reasonably accurate . for example , 
 the average acceptance ratio in figure 10 and table 1 shows the additional gains in sampling efficiency 
 with an improved cg model . however , if the cg model always generates 
 proposed moves that are energetically forbidden for the aa model then 
 most will be rejected , and the method becomes very inefficient . our 
 formal analysis based on eq 17 shows that the optimal choice is formally achieved when the 
 effective energy surface of the cg model corresponds to the exact 
 many - body pmf computed with respect to the cg degrees of freedom within 
 the aa system . while the exact pmf with respect to the cg degrees 
 of freedom is generally unknown , this analysis offers a direct route 
 to improve efficiency by progressively constructing a reasonable approximation 
 that captures the dominant effects . to achieve maximum sampling efficiency 
 with the cg - guided hybrid nemd - mc algorithm , one could adopt an adaptive 
 procedure in which the information accumulated from the aa simulation 
 would serve to progressively improve the parameters of the cg model . for example , one could adjust the function w(r )  on - the - fly  during the 
 simulation using a force - matching algorithm .
molecular 
 dynamics ( md ) trajectories based on a classical equation 
 of motion provide a straightforward , albeit somewhat inefficient approach , 
 to explore and sample the configurational space of a complex molecular 
 system . while a broad range of techniques can be used to accelerate 
 and enhance the sampling efficiency of classical simulations , only 
 algorithms that are consistent with the boltzmann equilibrium distribution 
 yield a proper statistical mechanical computational framework . here 
 , 
 a multiscale hybrid algorithm relying simultaneously on all - atom fine - grained 
 ( fg ) and coarse - grained ( cg ) representations of a system is designed 
 to improve sampling efficiency by combining the strength of nonequilibrium 
 molecular dynamics ( nemd ) and metropolis monte carlo ( mc ) . 
 this cg - guided 
 hybrid nemd - mc algorithm comprises six steps : ( 1 ) a fg configuration 
 of an atomic system is dynamically propagated for some period of time 
 using equilibrium md ; ( 2 ) the resulting fg configuration is mapped 
 onto a simplified cg model ; ( 3 ) the cg model is propagated for a brief 
 time interval to yield a new cg configuration ; ( 4 ) the resulting cg 
 configuration is used as a target to guide the evolution of the fg 
 system ; ( 5 ) the fg configuration ( from step 1 ) is driven via a nonequilibrium 
 md ( nemd ) simulation toward the cg target ; ( 6 ) the resulting fg configuration 
 at the end of the nemd trajectory is then accepted or rejected according 
 to a metropolis criterion before returning to step 1 . 
 a symmetric 
 two - ends momentum reversal prescription is used for the nemd trajectories 
 of the fg system to guarantee that the cg - guided hybrid nemd - mc algorithm 
 obeys microscopic detailed balance and rigorously yields the equilibrium 
 boltzmann distribution . 
 the enhanced sampling achieved with the method 
 is illustrated with a model system with hindered diffusion and explicit - solvent 
 peptide simulations . 
 illustrative tests indicate that the method can 
 yield a speedup of about 80 times for the model system and up to 21 
 times for polyalanine and ( aaqaa)3 in water .