molecular 
 dynamics ( md ) simulations are commonly used for ( a ) studying 
 dynamics of protein structures , ( b ) protein structure prediction , 
 and ( c ) calculating thermodynamic properties such as free energies , 
 enthalpy , and entropy of conformational states of proteins .
all - atom cartesian md simulations is a classical mechanics based 
 toolkit using cartesian coordinates as degrees of freedom .
the cartesian 
 md method is used widely for calculating statistical and thermodynamic 
 properties of materials and biomaterials .
one of the attractive features 
 of the cartesian all - atom dynamics model which uses absolute coordinates 
 is its mathematical simplicity .
there has been extensive development 
 of all - atom cartesian dynamics algorithms for thermostats and simulation 
 ensembles such as constant temperature ( nvt ) , constant pressure ( npt ) , 
 and constant stress ( grand canonical ensemble ) .
the thermodynamic properties calculated from these ensemble 
 simulations can be directly compared with experimental measurements . 
 additionally , constraints and/or bias potentials are often used in 
 cartesian md to increase the time step size , enhance conformational 
 sampling , and simulate large - scale conformational changes .
the addition 
 of these constraints into the cartesian dynamics makes the equations 
 of motion differential algebraic , requiring differential - algebraic 
 equation solvers that can adversely impact simulation robustness and 
 complexity .
bond length , bond angle , and torsion angle ( bat ) 
 relative coordinates are more natural than cartesian absolute coordinates 
 for describing the bonded structure of proteins .
md simulations in 
 bat coordinates are referred to as internal coordinate md ( icmd ) methods . 
 in icmd models of proteins ,
degrees of freedom that are nonessential 
 for effecting large scale conformational changes in proteins , such 
 as high frequency bond length degrees of freedom , can be constrained 
 by simply excluding them from the model . not only do the resulting icmd models have 
 fewer number of degrees of freedom , but they also retain the simpler 
 structure of ordinary differential equations instead of the more complex 
 differential - algebraic structure required for constrained cartesian 
 models . other advantages of icmd are the following:1.the low - frequency 
 torsional coordinates 
 allow larger time steps for integration of the equations of motion.2.conformational search is 
 more effective 
 in the low frequency torsional degrees of freedom and leads to significant conformational changes.3.enhanced sampling methods are more 
 effective when performed in torsional space.4.the six translation and orientation 
 degrees of freedom for a molecule are explicit coordinates rather 
 than implicit as in cartesian models .
this is useful when calculating 
 the conformational entropy using quasi - harmonic analysis.5.they provide a large range 
 of options 
 for selecting and controlling the granularity of the dynamics model.6.fixman potential corrections 
 for constraint - induced biases in the partition 
 function are easier to apply .
the low - frequency 
 torsional coordinates 
 allow larger time steps for integration of the equations of motion .
conformational search is 
 more effective 
 in the low frequency torsional degrees of freedom and leads to significant conformational changes .
the six translation and orientation 
 degrees of freedom for a molecule are explicit coordinates rather 
 than implicit as in cartesian models .
they provide a large range 
 of options 
 for selecting and controlling the granularity of the dynamics model .
fixman potential corrections 
 for constraint - induced biases in the partition 
 function are easier to apply .
the ability of the icmd simulation methods to control 
 the granularity of the dynamics model makes them highly suitable for 
 the development of multiscale methods and strategies for the simulation 
 of protein macromolecular complexes and polymeric materials .
our vision 
 is to develop a robust and advanced icmd method that capitalizes on 
 the inherent advantages of icmd , and to use these methods to develop 
 a comprehensive icmd simulation toolkit for tackling important problems 
 in structural biology .
the coupled nature of the icmd coordinates 
 and the higher complexity of icmd models however demands far more 
 sophisticated mathematical techniques and algorithms .
torsional md is a well - known example of 
 an icmd method that constrains 
 all bond lengths and bond angles .
many research groups have used the 
 torsional md method for simulating the dynamics of peptides and proteins in addition to clarifying the stubborn challenges 
 associated with the torsional md methods .
the wider usage of torsional md methods has been 
 hampered by serious bottlenecks that slowed their progress and viability 
 of icmd methods as a structural biology tool .
two long - standing concerns 
 have been ( a ) the mathematical and computational complexity of the 
 equations of motion associated with the icmd models and ( b ) the increased 
 rigidity in the torsional md dynamics resulting from keeping the bond 
 lengths and angles constrained , affecting the transition barriers 
 and probability density functions of states in a protein .
when 
 all the bond lengths and bond angles are constrained in torsional 
 md , the equations of motion in the dihedral angle space become computationally 
 expensive with the solution scaling as the cubic power of the number 
 of torsion degrees of freedom .
this has been a major bottleneck that 
 stunted the use and growth of icmd methods .
euler inverse mass operator ( gneimo ) method that is 
 based on the wealth of mathematical theory and analysis using spatial operator algebra techniques originally developed 
 for spacecraft and robot dynamics domains .
there has been a considerable amount of research over the years 
 on the development of the spatial operator algebra methods for the analysis of robot multibody 
 system dynamics .
the key insights that have been developed included 
 analytical techniques for the factorization and inversion of the mass 
 matrix for tree - topology systems .
these techniques have opened up 
 the opportunities to revisit and advance torsional md simulations 
 and make it computationally accessible to the wider community .
our 
 spatial operator algebra based icmd algorithm for solving the equations 
 of motion was the first to overcome the icmd computational cost bottleneck 
 by reducing costs to being just linearly , instead of cubically , proportional 
 to the number of degrees of freedom .
this low cost 
 algorithm has been adopted by other groups to implement torsional 
 md capability .
the gneimo method has been applied 
 to study a wide variety of structural biology problems such as protein 
 folding , protein structure refinement , and domain motion in proteins .
this algorithm is now being used routinely in refining 
 nmr structures ( in the software called cyana ) and x - ray crystal structures 
 ( in the software called nih - xplor ) . however , these two applications 
 require only short time scale dynamics and do not require calculation 
 of accurate thermodynamic properties from simulations .
although 
 the first challenge was eliminated making torsional md 
 simulations computationally feasible with the spatial operator algebra 
 based solution , the second bottleneck arising from the rigidity of 
 the model still persisted .
the increased rigidity of the dynamic model 
 stemming from freezing degrees of freedom led to fewer dihedral transitions 
 and to systematic errors in the probability density functions for 
 proteins and peptides .
researchers have observed that the use of rigid 
 constraints in both bond lengths and bond angles in the torsional 
 md simulations alters the potential energy surface and the free energy 
 surface of the system compared with unconstrained md simulations . addressing the rigidity of the dynamics model imposed 
 in the torsional md simulations , fixman in the 1970s proposed a compensating 
 potential that rigorously corrects for treating stiff and uncoupled 
 bond angles as rigid in the icmd model , and generates a partition 
 function that gives probability density functions closer to that of 
 all - atom cartesian simulations .
while the fixman potential removes such biases for those bond angles 
 that are stiff , it remained intractable and hence not tested for larger 
 branched molecules . by developing a spatial operator algebra based 
 general purpose and low - cost computational method to calculate the 
 fixman potential for all linear and branched molecules , we have solved the longstanding problem of using the fixman potential 
 in icmd simulations to remove constraint - induced biases in thermodynamic 
 properties and the probability density function .
our studies verify 
 that , as predicted , the inclusion of the fixman potential recovers 
 the equilibrium probability density function of the conformational 
 states , the transition barrier crossing rates , and the free energy 
 surface for serial and branched molecules . to the best of our knowledge , this is the first time that the fixman 
 potential has been used for icmd of large and realistic branched polymers . overcoming the two long - standing bottlenecks along with several 
 theoretical and algorithmic advances in the gneimo icmd method made 
 in the past five years
has led to a robust long time scale icmd simulation 
 method and an associated software package for use to study protein 
 structural dynamics .
the highlights of these methods 
 are described in section 2 , and the highlights 
 of the results for several protein simulation applications and their 
 comparison to experimental findings are given in section 3 .
the current state of the icmd methods and the 
 areas that need further development are discussed in section 4 .
the icmd 
 method is a molecular dynamics simulation method performed 
 in bat coordinates also known as internal coordinates .
the gneimo 
 method is a generalized icmd method based on spatial operator algebra , 
 for performing multibody dynamics of macromolecules .
constraints on 
 the high frequency bond lengths can be placed in the gneimo method 
 to perform icmd simulations with the bond angles and torsion angles 
 as degrees of freedom . if both the bond lengths and bond angles are 
 kept rigid , the resulting torsional md method is also supported by 
 gneimo . in the gneimo torsional md method ,
the bond length and bond 
 angle degrees of freedom are treated as rigid and the macromolecule 
 of tree topology is modeled as a collection of rigid bodies ( of varied 
 sizes ) connected by flexible hinges .
the rigid bodies also known as clusters can vary in size and shape ranging from a single 
 atom to a methyl group , or a helix or an entire domain of a protein . 

however , these hinges have one to six degrees of freedom , 
 with the one degree of freedom being just the torsion angle .
when 
 the hinge has six degrees of freedom , it allows the bond stretch , 
 bond angle bending , and torsion about the bond connecting the clusters . 

the default cluster model of proteins in the gneimo method is shown 
 in figure 1 . when the bond lengths and 
 bond angles are treated as rigid
, the 
 equations of motion in icmd become coupled and are shown in eq 11where  is the vector of the generalized 
 coordinates ( e.g. , torsional angles ) ,   denotes 
 the vector of generalized forces 
 ( e.g. , torques ) ,   denotes the mass matrix ( moment of inertia 
 tensor ) , and   includes 
 the velocity dependent coriolis 
 forces .
the dynamics of motion is obtained by solving eq 1 for the  acceleration and integrating them 
 to obtain new velocities and coordinates . with conventional algorithms , 

the equations of motion involve calculating the inverse of the dense 
 mass matrix that scales as the cubic power of the number of degrees 
 of freedom .
the gneimo method uses a spatial 
 operator algebra based method to derive an analytical expression for 
 the inverse of the mass matrix followed by the following expression 
 for 2the ,  , , etc .
, terms 
 in the above expression are 
 associated with mass matrix related factorizations and are described 
 in detail in refs ( 5 and 19 ) .
standard clustering 
 scheme used for torsional md in gneimo shown 
 for the tripeptide ala - tyr - ala .
the rigid clusters are shown in a 
 single color , and these clustered atoms move together as a single 
 unit . the gray rods and arrows represent hinges connecting two clusters . 

the expression on the 
 right can be evaluated using recursive algorithms 
 whose cost scales just linearly with the number of degrees of freedom 
 and thus provides a computationally tractable method for solving the 
 equations .
this recursive algorithm described in detail in ref ( 5 ) avoids the computationally 
 intensive inversion of the dense mass matrix shown in the matrix equation 
 in eq 1 .
we further extended this method to 
 simulate the canonical n , v , t ensemble with constant temperature dynamics using the 
 nos  hoover thermostat method . 

hoover 
 icmd method are given below.34 here   is the additional 
 frictional force term 
 due to the canonical ensemble which is 
 dependent on  , the dynamic variable representing the thermostat , 
  is the mass parameter of the thermostat , t is 
 the instantaneous temperature , and tb is the 
 thermostat temperature . since the gneimo nos 
hoover 
 equations of motion shown above are similar to the equations of motion 
 in eq 1 , all the spatial operator equations 
 and factorization discussed in ref ( 5 ) hold good for these equations as well . the nos 
hoover 
 thermostat mass parameter  was optimized to be 10 times the 
 time step size for torsional md simulations .
the accuracy and stability of the gneimo torsional md simulations 
 were measured by the conservation of the total hamiltonian and the 
 extent of fluctuations in the temperature .
other 
 research groups such as brooks and co - workers and schweiters and co - workers 
 adapted the gneimo algorithm into charmm and nih - xplor programs for 
 icmd simulations .
however , the systematic 
 biases in the probability density functions 
 calculated from icmd simulations arising from the rigidity of the 
 model remain to be solved .
to rigorously correct for the systematic 
 biases in the probability density function caused by treating stiff 
 bond angles as rigid , fixman proposed a compensating potential of 
 the form5where   denotes the 
 mass matrix in the full bat 
 coordinates , q0 the coordinates for the 
 frozen degrees of freedom , k the boltzmann constant , 
 and t the temperature . while the fixman potential removes such biases , it was computationally 
 intractable for generalized branched molecules and hence remained 
 untested and unused for several decades for icmd simulations , except 
 for small model systems such as c4 or c5 .
we derived a spatial operator algebra based 
 algorithm that calculates the fixman potential with just 24% additional 
 cost in the computation time .
this method which we 
 call gneimo - fixman is rigorous , general purpose , low - cost , and applicable 
 to general linear and branched systems with little additional cost . 

more significantly , the gneimo - fixman method is also able to compute 
 the partial derivatives of the fixman potential .
these partial derivatives 
 define the additional forces , i.e. , the fixman torque , to be applied 
 within constrained md simulations .
this makes it feasible to use the 
 fixman potential during the torsional md simulations and test the 
 effect of the fixman potential on recovering the accuracy of the probability 
 density functions and the dihedral transition rates . to study the effectiveness of 
 the fixman potential in recovering the thermodynamic probability density 
 function of torsion angles , we performed torsional md simulations 
 with langevin force for various linear and branched molecular systems 
 without force fields .
the effect of the fixman potential can be most 
 clearly seen by comparing the probability density functions calculated 
 from torsional md simulations with the fixman potential to the probability 
 density functions calculated from unconstrained cartesian md simulations . 

figure 2 shows the joint probability distribution function 
 of the two main chain torsion angles in alanine dipeptide .
it is seen 
 from the plots that the torsional md introduces biases in the joint 
 probability density function that are removed when the torsional md 
 simulations are performed with the fixman potential and the torques . 

the magnitude of the fixman potential is much smaller in comparison 
 to the potential energy from the all - atom force field .
there are some 
 bond angles in the protein that are weakly coupled to the dihedrals 
 and the nonbond forces while other bond angles such as backbone bond 
 angles that show strong coupling to the dihedrals and the nonbond 
 forces .
although the fixman potential corrects for the errors in the 
 probability density function stemming from treating stiff and uncoupled 
 bond angles as rigid , it does not do so for soft bond angles that 
 couple to the torsion angles or the nonbond interactions .
there 
 are two approaches for eliminating the bias in the potential energy 
 surface imposed by treating the bond angles that are coupled to torsions 
 and nonbond interactions , as rigid .
the first approach is to open 
 up some of the bond angles and treat them as movable degrees of freedom 
 in the icmd simulations .
this approach reduces rigidity 
 and the error in probability density function and transition rates , 
 with little impact on time step size .
the second approach is to use 
 correction torsional potentials that compensate for the rigidity in 
 the coupled bond angles to recover the probability density function . 

examples of this approach are the ecepp family of force fields or icmff that correct the torsional angle potential 
 for all of the bond angles .
the icmff corrects the force constants 
 used for the torsional angle potentials by refitting the torsional 
 energy curve obtained from the rigid model to the torsional energy 
 curve calculated using the flexible model .
the force constants thus 
 obtained will reproduce the torsion energy barriers of the flexible 
 cartesian model .
chen et al . showed that this method enhances the 
 number of dihedral transitions during the torsional md simulations . 

the possible caveat with this method is that it is not rigorous and 
 is system specific .
the icmff has been tested for alanine dipeptide 
 and remains to be tested for larger systems . in summary , it is necessary 
 to treat the bond angles that show strong coupling to torsion angles 
 as flexible degrees of freedom in order to achieve thermodynamic accuracy 
 in the probability density function . however , this is not required 
 if the torsional md simulations are used for structure prediction 
 applications .
this is because the goal in structure prediction applications , 
 such as refining the loop structures in proteins , is to enrich the 
 native ensemble by widening the conformational sampling and generally 
 not in recovering the accuracy of the partition function of the cartesian 
 md simulations .
besides the two 
 major challenges ( described above ) that have hampered icmd method 
 development , there were other developments that were required to make 
 the icmd method robust and suitable for long time scale simulations . 

several research groups have shown that constrained icmd models are 
 not the limiting case of stiff cartesian models . to address the differences between constrained icmd models and cartesian 
 models of the molecule
, we showed that the application of the conventional 
 equipartition theorem that is based on the cartesian model to the 
 icmd model in internal coordinates does not yield an equipartition 
 principle analogous to that for cartesian models .
instead , the ensemble 
 averages involve configuration dependent coupled coordinates , that 
 are not easy to interpret or use .
therefore , we introduced a coordinate 
 transformation to modal coordinates that transforms the system s 
 kinetic energy into a decoupled form similar to that for cartesian 
 models .
additionally , we showed that using the  modal velocity  
 coordinates one can derive an equipartition principle for the icmd 
 model that is analogous to the one for cartesian models .
this principle 
 holds even though the modal coordinates are not canonical coordinates 
 in the hamiltonian sense . the equipartition principle based on modal 
 coordinates provides a method for the thermodynamically correct initialization 
 of velocities in icmd simulations .
joint probability 
 distribution function of the two backbone torsion 
 angles in alanine dipeptide : ( a ) cartesian simulations ; ( b ) icmd simulations ; 
 ( c ) icmd + fixman simulations .
the bin size for each axis is d 
 = 18. material : reprinted with permission from ref ( 18 ) .
we have 
 developed a modular software package called gneimosim for efficient implementation of icmd simulations .
gneimosim is quite possibly the first software package 
 to include advanced features such as the icmd equipartition principle 
 and methods for including the fixman potential to eliminate systematic 
 statistical biases introduced by the use of hard constraints .
moreover , 
 gneimosim s architecture allows it to be extended and easily 
 interfaced with third party force field packages for icmd simulations . 

the availability of a comprehensive 
 python interface to the underlying c++ classes and their methods provides 
 a powerful and versatile mechanism for users to develop simulation 
 scripts to configure the simulation as well as to control the simulation 
 flow .
the alpha version of the gneimosim software is available 
 for use for the research community at http://dartslab.jpl.nasa.gov/gneimo/index.php . in figure 3
, the blue line shows the 
 average run time for a simulation for the start - up and solution of 
 the equations of motion in the gneimo algorithm without using a force 
 field .
it shows the effective linear scaling of the dynamics equations 
 solver cost for the gneimo method in gneimosim .
the average run time 
 for calculating the forces using openmm with gbsa ( on a gpu ) and using 
 lammps and rosetta force fields on a cpu are shown in yellow , red , 
 and green lines , respectively .
the cost of running the equations of 
 motion solver in gneimo is equivalent to the force field cost of using 
 the amber force field with the generalized born solvation module on 
 gpu using the openmm force field engine .
however , for larger systems , 
 the cost of the gneimo equations of motion solver will have a linear 
 dependence while the cost of the calculation of forces will show a 
 higher order increase in computational time .
additionally , the increase 
 in computational time for icmd is compensated to an extent by the 
 ability to take larger time steps than are possible in a cartesian 
 simulation .
we have implemented 
 in gneimosim several methods that are routinely used in md simulations 
 such as1.advanced 
 integrators including runge 
kutta , 
 lobatto , adaptive cvode , and verlet integrators for stable long time 
 scale ( microseconds ) torsional md simulations.2.an adaptation 
 of the generalized born 
 solvation method ( gbsa ) for implicit solvation.3.support for multiple molecules of any 
 type , including explicit solvent.4.a temperature - based replica - exchange 
 ( remd ) method in which temperatures may 
 be switched randomly or probabilistically using the metropolis algorithm.5.periodic boundary conditions 
 for simulations 
 with explicit water.6.support for standard cartesian simulations 
 for comparison.7.user 
 defined harmonic distance restraints 
 between pairs of atoms.8.more recently , we have added langevin 
 dynamics and accelerated md ( amd ) methods to gneimosim .
kutta , 
 lobatto , adaptive cvode , and verlet integrators for stable long time 
 scale ( microseconds ) torsional md simulations .
an adaptation 
 of the generalized born 
 solvation method ( gbsa ) for implicit solvation .
a temperature - based replica - exchange 
 ( remd ) method in which temperatures may 
 be switched randomly or probabilistically using the metropolis algorithm .
more recently , we have added langevin 
 dynamics and accelerated md ( amd ) methods to gneimosim .
average 
 run times per step for the gneimosim simulation using different 
 force fields .
the 
 bat coordinates are natural coordinates to describe a bonded system 
 like proteins and polymers .
they also lend themselves readily to freezing 
 and thawing any internal coordinate degree of freedom .
for example , 
 many torsional degrees of freedom can be kept constrained during dynamics 
 such as freezing a whole helixin a protein structure , if needed .
this 
 freeze and thaw can be done at the start of the simulation or during 
 the course of the md simulations in the gneimo icmd method .
we refer 
 to the freeze and thaw technique on the fly during the dynamics simulations 
 as  dynamic clustering  .
the dynamic clustering scheme can be applied automatically on the 
 basis of user specified criteria to freeze and thaw internal coordinate 
 degrees of freedom .
the default clustering scheme for the gneimo torsional 
 md treats all the torsions of the backbone and main chain as degrees 
 of freedom .
all - torsion md simulations may not be adequate for simulating 
 long time scale events such as large conformational changes in proteins . 

adaptable clustering strategies that allow change in the cluster model 
 of the protein during the simulations are more suitable for simulating 
 large conformational changes .
wagner et al . used criteria to freeze secondary 
 structure elements when performing folding simulations at high temperature 
 replicas and allow them to thaw at lower temperature replicas .
however , it should be noted that , since the 
 freeze and thaw mechanism can alter the pathway of the process being 
 simulated , it is advisable to use this method for conformational search 
 purposes only . other criterion to freeze or thaw a degree of freedom 

is by using an upper threshold for the velocity , a hinge can be treated 
 as flexible if its velocity stays closer to the threshold velocity . 

this ensures that torsions that undergo significant motion are not 
 treated as rigid . additionally , monitoring the accumulation of stress 
 forces at a hinge that is treated rigid
we have 
 used the freeze and thaw strategy to fold four proteins starting from 
 their respective extended structure .
the secondary structure elements 
 were predicted using sequence information and built into the extended 
 structure .
gneimo torsional md with the replica exchange method was 
 used to run the folding simulations . the adaptive time step cvode integrator using the adams 
the cvode adaptive 
 integrator ensures stable simulations while allowing rapid conformational 
 sampling of clustered rigid bodies , by adjusting the time step size 
 to avoid clashes that destabilize the md simulations . using the md 
 trajectory of the b domain of staphylococcal protein a , we calculated 
 the population density histogram shown in figure 4 .
this figure also shows the corresponding representative 
 structures from each of the conformational clusters .
the x - axis is the root - mean - square deviation ( rmsd ) of the backbone atom 
 coordinates from the crystal structure .
we used 12 temperature replicas ranging from 
 300 to 1050 k. the folding process begins with the extended structure 
 with just the predicted helices using secondary structure prediction 
 methods .
initially , a small number of interhelical contacts are made 
 that collapses the structure to a rmsd range of 1216  . 
 at 811  ,
slightly incorrect protein topologies are explored 
 that finally fold to the correct topology below 7  .
this example 
 demonstrates the rapid conformational sampling realized using the 
 dynamic clustering strategy .
large scale domain 
 motion caused by the low frequency modes in a protein are often difficult 
 to capture with cartesian md simulations .
the conformational transitions 
 between the substates are often affected by long time scale and rare 
 events difficult to capture within the relatively short ( hundreds 
 of nanoseconds ) time scale md simulations .
enhanced sampling methods 
 are therefore often used to simulate the conformational transitions 
 in proteins .
some of the enhanced sampling 
 methods require explicit knowledge of the conformations for which 
 the transitions are being simulated , and therefore not readily applicable 
 to unknown systems .
the temperature based remd method when used with 
 cartesian md simulations often leads to unraveling and reforming of 
 secondary structure motifs .
thus , the transitions are wasted in high 
 frequency modes and the extra thermal energy is not focused on searching 
 in the low frequency modes .
we studied the use of remd with gneimo 
 torsional md simulations to see how effective they are in conformational 
 search in two highly flexible proteins , namely , fasciculin and calmodulin .
calmodulin is a calcium signaling protein that 
 belongs to the ef - hand family of proteins .
the three - dimensional structure 
 of calmodulin consists of two domains , namely , the amino terminus 
 ( n - terminus ) domain and the carboxy - terminus ( c - terminus ) domain , 
 that are connected by a long helical stretch .
calmodulin is a dynamic 
 protein that exhibits major conformational changes in response to 
 calcium binding as revealed by nmr studies . upon binding 
 to calcium ,
calmodulin samples a large ensemble of structures that 
 allows it to bind to a wide range of proteins .
there are several experimental 
 studies on the dynamics and conformational changes in calmodulin . 

these studies show that there are two major steps involved in the 
 conformational changes : ( 1 ) upon removal of the bound calcium , the 
 central helix that connects the carboxy terminal domain and the amino 
 terminal domain collapses followed by ( 2 ) the dynamics of the n - terminus 
 domain relative to the static c - terminus domain .
cartesian md simulations in explicit solvent showed the collapse 
 of the central helix connecting the two domains but could not map 
 the entire dynamics ensemble obtained from the nmr study .
the gneimo - remd simulations with no bias potential were performed 
 on the entire calmodulin protein .
the trajectory showed conformational 
 sampling of the collapse of the central helix , and the sampling of 
 various conformations resulting from the relative reorientation of 
 the two domains .
although the collapse was observed in the cartesian 
 md , the flexibility of the n - terminus domain was observed only in 
 the gneimo - remd simulations .
comparison of the ensemble of conformations 
 of the carboxy terminus domain to the ensemble of conformations from 
 nmr experiments is shown in figure 5 .
it is 
 seen that the gneimo conformations ( shown in blue in figure 5 ) cover most of the nmr conformations .
backbone crmsd 
 histogram of a  dynamic clustering  
 replica - exchange simulation of 1bdd , beginning from an extended conformation 
 containing predicted secondary structure elements with representative 
 structures .
overlay of the conformations 
 sampled in the gneimo - remd torsional 
 md simulations ( shown in blue ) of calmodulin without calcium , to the 
 nmr structures pdb i d : 1dmo ( shown in gold ) .
quantitative comparison of the average hydrogen bond distances 
 between residues showing hydrogen bonds in the nmr structures ( pdb i d : 1dmo ) to those in the gneimo - generated trajectories 
 shows that about half of the average distances calculated fall within 
 one standard deviation of the corresponding distances in the nmr structures . 

gneimo torsional md simulations on proteins such as crambin and bpti 
 at 310 k have shown that the correlations in the torsional motions 
 of residues that are farther apart in space can be captured in a shorter 
 simulation time scale than with cartesian md simulations .
for example , 
 in the case of bpti , analysis of correlation in the backbone torsion 
 angles of residues in two loops connected 
 by a disulfide bond was shown to be strong in the millisecond level 
 simulations of shaw and co - workers .
we 
 have shown that we could capture the torsional angle correlations 
 between residues 9 to 18 and residues 35 to 40 in 100 ns of nvt gneimo 
 torsional md simulations at 310 k. the 
 gneimo algorithm has been implemented in the software package 
 cyana for nmr structure refinement and 
 in nih - xplor for x - ray crystal structure refinement . brooks and co - workers have tested torsional md simulations 
 based on the gneimo algorithm for few proteins . however
, this method has not been tested for refining protein 
 homology models of lower accuracy compared to their respective crystal 
 structures . due to the large number of protein structures of similar 
 proteins available in the protein data bank ,
homology or template 
 based modeling of related protein sequences has been shown to be the 
 most feasible and accurate method of predicting protein structures . 

however , the target sequence has to have a > 60% 
 sequence similarity to the template crystal structure .
the accuracy 
 of the model derived using homology modeling methods depends on the 
 sequence similarity between the template and the target .
there is 
 an increasing need for high throughput computational methods that 
 refine low accuracy homology models .
using 
 cartesian md simulation methods leads to refined structures only when 
 used with restraints  with the choice of restraints remaining 
 arbitrary and system dependent and hence difficult to automate .
starting from the homology based models , robust structure refinement 
 methods are important to derive high useful accuracy models .
an efficient conformational sampling method and an accurate energy 
 function to identify the native structure are the two important components 
 of a structure refinement algorithm .
we have evaluated the usefulness 
 of the conformational sampling afforded by gneimo - remd with the generalized 
 born solvation method in enriching the refined structures compared 
 to the starting decoy .
the gneimo - remd method 
 leads to a refinement of up to 1.3  for 28 out of 30 casp target 
 proteins .
these torsional md simulations using gneimo were done without 
 any experimental information as restraints .
figure 6 shows the contact map for one target ( t0453 ) as an example 
 of the refinement using the gneimo method .
figure 6 shows the contact map for refinement of the target protein 
 t0453 . the contact map shows 
 far is a given residue in the gneimo refined model from its crystal 
 structure .
the white color regions indicate the residue positions 
 that are closer to the native structure than the residues in the deep 
 red regions .
it is evident from figure 6 that 
 the gneimo refined structure shows significant refinement in the packing 
 of the loop structure against the core of the protein .
the long - range 
 contacts between the loop residues 30 to 40 to those residues between 
 40 and 60 improve remarkably from 14 to 16  to 2 to 4  , 
 as seen in figure 6 .
the 
 contact map shows the distance to the decoy structure , and the contact 
 map in the third column shows how far each residue is refined corresponding 
 to the native structure .
while the gneimo method leads to substantial refinement in 
 the 
 loop packing , the unconstrained cartesian remd method shows unraveling 
 and reformation of secondary structure elements and hence is less 
 effective in conformational search .
thus , starting from homology models 
 of varying accuracy , the gneimo method shows improvement in structure 
 refinement without unraveling the structures .
the overall extent of 
 refinement across many targets was modest and needs further improvement . 

the icmd simulation 
 method is not a replacement for the cartesian 
 md method .
we envision developing an icmd simulation method to perform 
 md simulations in bond , angle , and torsion ( bat ) coordinate systems . 

such a md simulation method is highly suitable for selecting and controlling 
 the granularity of the dynamics model of proteins and polymers , and 
 is a foundation for the development of multiscale simulation methods 
 for the simulation of protein macromolecular complexes .
the multiscale 
 icmd simulation method can also be used for real space refinement 
 of structural model fitting to low resolution crystallography or electron 
 microscopy data of large protein complexes .
currently , 
 studies of the large scale dynamics that govern protein function commonly 
 use all - atom cartesian md simulations , often with geometric constraints . 

such geometric constraints can adversely impact the integration time - step 
 size and the stability of the dynamics . on the other hand ,
our recent theoretical developments based 
 on spatial operator algebra have led to a robust long time scale icmd 
 simulation method gneimo and associated toolkit known as gneimosim . 
 we have demonstrated the power of the gneimo icmd technique with applications 
 to protein domain motion studies and homology model refinement studies . 

the results of these applications have shown great promise in the 
 wider use of this icmd technique .
however , there remain extensions 
 of the methods that need to be addressed to harness the advantages 
 of icmd methods to their fullest .
some of the possible extensions 
 are the following:1.to calculate accurate conformational 
 entropy using the gneimo torsional md simulation trajectories with 
 the quasi - harmonic analysis ( qha ) method to improve the accuracy while 
 correctly including traditionally difficult metric tensor correction 
 terms.2.the second advancement in icmd simulations 
 will be to extend the icmd simulation method to allow movement of 
 certain bond angles that show strong coupling to the dihedral angles . 

for example , it has been shown that the bond angle hinged on the c 
 atoms of a protein have strong coupling to the backbone dihedral angles .
we will extend the gneimo icmd method to free up desired bond angle 
 degrees of freedom .
this will lead to a comprehensive suite of hybrid 
 icmd models that bridge the wide gap between coarse grain torsional 
 md and fine grain all - atom md models .
these hybrid icmd models , with 
 all bond angles open , can also be used for accurately calculating 
 conformational entropy as discussed here.3.gneimosim software has enhanced sampling methods such as langevin 
 dynamics , replica exchange method implemented .
gneimo can be readily 
 combined with other enhanced sampling techniques such as steered md , 
 umbrella sampling , or accelerated md ( amd ) to name a few . to calculate accurate conformational 
 entropy using the gneimo torsional md simulation trajectories with 
 the quasi - harmonic analysis ( qha ) method to improve the accuracy while 
 correctly including traditionally difficult metric tensor correction 
 terms .
the second advancement in icmd simulations 
 will be to extend the icmd simulation method to allow movement of 
 certain bond angles that show strong coupling to the dihedral angles . 

for example , it has been shown that the bond angle hinged on the c 
 atoms of a protein have strong coupling to the backbone dihedral angles .
we will extend the gneimo icmd method to free up desired bond angle 
 degrees of freedom .
this will lead to a comprehensive suite of hybrid 
 icmd models that bridge the wide gap between coarse grain torsional 
 md and fine grain all - atom md models .
these hybrid icmd models , with 
 all bond angles open , can also be used for accurately calculating 
 conformational entropy as discussed here .
gneimosim software has enhanced sampling methods such as langevin 
 dynamics , replica exchange method implemented .
gneimo can be readily 
 combined with other enhanced sampling techniques such as steered md , 
 umbrella sampling , or accelerated md ( amd ) to name a few . with these improvements
, 
 the gneimo method will be useful for applications 
 to large protein dynamics simulations with explicit solvent .
the gneimo 
 method can be used with any type of force field such as an all - atom 
 force field or a coarse grain force field .
we also envision using 
 the gneimo torsional md method combined with the torsional monte carlo 
 method for protein structure prediction applications . in summary , 

the gneimo method offers a comprehensive icmd simulation method that 
 allows easy coarse graining of the dynamics model ranging from all - atom 
 to large domains of proteins being treated as clusters .
 internal coordinates such as bond 
 lengths , bond angles , and torsion 
 angles ( bat ) are natural coordinates for describing a bonded molecular 
 system . however , the molecular dynamics ( md ) simulation methods that 
 are widely used for proteins , dna , and polymers are based on cartesian 
 coordinates owing to the mathematical simplicity of the equations 
 of motion . 
 however , constraints are often needed with cartesian md 
 simulations to enhance the conformational sampling . 
 this makes the 
 equations of motion in the cartesian coordinates differential - algebraic , 
 which adversely impacts the complexity and the robustness of the simulations . 
 on the other hand , 
 constraints can be easily placed in bat coordinates 
 by removing the degrees of freedom that need to be constrained . 
 thus , 
 the internal coordinate md ( icmd ) offers an attractive alternative 
 to cartesian coordinate md for developing multiscale md method . 
 the 
 torsional md method is a special adaptation of the icmd method , where 
 all the bond lengths and bond angles are kept rigid . 
 the advantages 
 of icmd simulation methods are the longer time step size afforded 
 by freezing high frequency degrees of freedom and performing a conformational 
 search in the more important low frequency torsional degrees of freedom . 
 
 however , the advancements in the icmd simulations have been slow and 
 stifled by long - standing mathematical bottlenecks . in this review 
 , 
 we summarize the recent mathematical advancements we have made based 
 on spatial operator algebra , in developing a robust long time scale 
 icmd simulation toolkit useful for various applications 
 . we also present 
 the applications of icmd simulations to study conformational changes 
 in proteins and protein structure refinement . 
 we review the advantages 
 of the icmd simulations over the cartesian simulations when used with 
 enhanced sampling methods and project the future use of icmd simulations 
 in protein dynamics . 
Introduction
Methods
Results
and Discussion
Conclusions