rotary motions , in particular , are essential for the function of many motor proteins .
these nanomotors use the free energy of chemical reactions or ion concentration gradients to generate mechanical torque .
rotary mechanisms were unequivocally demonstrated for three molecular engines , the fo and f1 motors in f - atp synthase ( f - atpase ) and the bacterial flagellar motor.(3 ) recently , rotary motion was also shown for the v1 portion of the prokaryotic homologue of the vacuolar atpase ( v - atpase).(4 ) other motor proteins that are assumed to be rotary include dna helicases(5 ) and proteins that translocate viral dna into preformed capsids .
the molecular mechanisms by which chemical reactions or transmembrane gradients drive protein rotary motions are in most cases not understood in full detail.(9 ) also , these often quite complex motions are typically too slow or infrequent to be accessible to equilibrium molecular dynamics ( md ) simulations . to overcome this limitation
, techniques have been developed to exert external forces or torques to certain subunits to induce rotation and/or to increase its rate without severely perturbing the nature of the involved structural changes .
this approach has also been used to simulate experiments in which biomolecules , such as proteins or dna , are mechanically driven to rotate by externally applied torques by single molecule manipulation techniques.(16 ) in one impressive example , the f1 portion of atp synthase ( f1-atpase ) has been shown to produce atp when the  subunit is enforced to rotate using magnetic tweezers.(17 ) with exceptions,(18 ) in most simulations involving external torque , a fixed ,  stiff  rotation axis has been used so far ( dashed line in figure 1a ) . as shown in the figure
, this approach does not properly describe situations such as f1-atpase , where the rotating part flexibly adapts ( dotted lines ) to the steric restraints set by the bearing ( gray ) . to more realistically describe biomolecular rotations ,
we have therefore developed a flexible axis rotation technique that ( i ) exerts torque with a curved axis that flexibly fits the shape of an arbitrarily shaped cavity ( figure 1a ) , ( ii ) avoids any impact or bias previously introduced by the necessary choice of the pivot for the axis , ( iii ) perturbs the internal dynamics and flexibility of the rotated structure as little as possible , and ( iv ) allows the curvature of the axis to adapt to structural changes of the bearing . in summary , a rotated fragment such as the  subunit inside the atpase 33 stator should deform like a rotating pipe - cleaner .
( a ) rotating the sketched shape inside the white tubular cavity creates severe artifacts when a conventional fixed rotation axis ( dashed ) is used .
more realistically , the shape would revolve like a flexible pipe - cleaner ( dotted ) inside the bearing ( gray ) .
( b ) fixed rotation around an axis v with a pivot point specified by the vector u. ( c ) subdividing the rotating fragment into slabs with separate rotation axes (  ) and pivot points (  ) for each slab allows for the required flexibility .
the distance between two slabs with indices n and n + 1 is x . to clarify notation and to explain the basic ingredients needed for the flexible technique , we start with a recapitulation of the established fixed axis rotation , as implemented , e.g. , in namd(21 ) or ego.(22 ) from these notions , several more complex potentials will be developed and characterized , and the resulting forces will be derived
. we will then motivate and describe in detail the flexible axis approach , for which we present two different variants . after outlining details of our gromacs implementation
, we will apply flexible axis rotation to the f1-atpase molecular motor and test if our approach is indeed capable of providing more accurate torque or free energy profiles .
in the established fixed axis approach ( figure 1b ) , torque on a group of n atoms with positions xi ( denoted  rotation group  ) is applied by rotating a reference set of atomic positions 
usually their initial positions yi0at a constant angular velocity  around an axis defined by a direction vector v and a pivot point u. to that aim , each atom with position xi is attracted by a  virtual spring 
potential to its moving reference position yi = (t ) ( yi0  u ) , where (t ) is a matrix that describes the rotation around the axis . in the simplest case ,
the  springs  are described by a harmonic potentialwith optional mass - weighted prefactors wi = nmi / m with total mass m = i=1nmi .
the rotation matrix (t ) iswhere vx , vy , and vz are the components of the normalized rotation vector v and  : = 1  cos ( t ) . as illustrated in figure 2a for a single atom j , the rotation matrix (t ) operates on the initial reference positions yj0 = xj(t0 ) of atom j at t = t0 . at a later time t
, the reference position has rotated away from its initial place ( along the blue dashed line ) , resulting in the forcewhich is directed toward the reference position .
all four potentials v ( color coded ) are shown for a single atom at position xj(t ) .
( a ) isotropic potential v , ( b ) radial motion potential v and flexible potential v , ( c , d ) radial motion 2 potential v and flexible 2 potential v for  = 0 nm ( c ) and  = 0.01 nm ( d ) .
the rotation axis is perpendicular to the plane and marked by . the light gray contours indicate boltzmann factors e in the xj plane for t = 300 k and k = 200 kj/(molnm ) . the green arrow shows the direction of the force fj acting on atom j ;
we first address the bias introduced by an arbitrary choice of the pivot vector u. this arbitrariness is avoided by defining as the pivot the center of mass xc of the rotation groupwhich yields the  pivot - free 
potentialwith forceswithout mass - weighting , the pivot xc is the geometrical center of the group . obviously ,
the forces generated by the isotropic potentials ( eqs 1 and 4 ) also contain components parallel to the rotation axis and thereby restrain motions along the axis of either the whole rotation group ( in case of v ) or within the rotation group ( in case of v ) . for cases where unrestrained motion along the axis is preferred ,
we have implemented a  parallel motion  variant by eliminating all components parallel to the rotation axis for the potential .
this is achieved by projecting the distance vectors between reference and actual positions : onto the plane perpendicular to the rotation vectoryieldingand similarly replacing in eq 8 the fixed pivot u with the center of mass xc yields the pivot - free variant of the parallel motion potential .
withthe respective potential and forces are in the above variants , the minimum of the rotation potential is either a single point at the reference position yi ( for the isotropic potentials ) or a single line through yi parallel to the rotation axis ( for the parallel motion potentials ) . as a result ,
the two subsequent types of rotation potentials , v and v , drastically reduce or even eliminate this effect .
the first variant , v ( figure 2b ) , eliminates all force components parallel to the vector connecting the reference atom and the rotation axiswiththis variant depends only on the distance pi(xiu ) of atom i from the plane spanned by v and (t)(yi0u ) .
the resulting force is proceeding similar to the pivot - free isotropic potential yields a pivot - free version of the above potential .
withthe potential and force for the pivot free variant of the radial motion potential read as seen in figure 2b , the force resulting from v still contains a small , second - order radial component . in most cases , this perturbation is tolerable ; if not , the following alternative , v , fully eliminates the radial contribution to the force , as depicted in figure 2c , where a small parameter  has been introduced to avoid singularities . for  = 0 nm ,
the equipotential planes are spanned by xi  u and v , yielding a force perpendicular to xi  u , thus not contracting or expanding structural parts that moved away from or toward the rotation axis .
we note that this variant is particularly suitable for free energy calculations via umbrella sampling techniques,(25 ) because the radial orientation of the equipotential planes shown in figure 2c guarantees statistically consistent sampling of adjacent umbrella windows , as required for a consistent definition of the free energy profile via subspace projection . to see why this is actually the case , note that consistent umbrella sampling requires that for adjacent umbrella windows the  stack  of ( 3n  1 dimensional ) configurational subspaces defined by the values of the chosen reaction coordinate agrees , subspace by subspace , with the one defined by the values of the umbrella potential .
this in turn requires that the equipotential planes shown in figure 2 coincide with those of a rotated potential , which is obviously the case for figure 2c , but not for figure 2a or b. choosing a small positive  ( e.g. ,  = 0.01 nm , figure 2d ) in the denominator of eq 19 yields a well - defined potential and continuous forces also close to the rotation axis , which is not the case for  = 0 nm ( figure 2c ) .
withthe force on atom j reads the pivot free variant of the above potential iswiththe force on atom j reads
as sketched in figure 1a , b , the rigid body behavior of the fixed axis rotation scheme is a drawback for many applications .
in particular , deformations of the rotation group are suppressed when the equilibrium atom positions directly depend on the reference positions . to avoid this limitation , eqs 18 and 24
will now be generalized toward a  flexible axis  , as sketched in figure 1c .
this will be achieved by subdividing the rotation group into a set of equidistant slabs perpendicular to the rotation vector , and by applying a separate rotation potential to each of these slabs .
figure 1c shows the midplanes of the slabs as dotted straight lines and the centers as thick black dots . to avoid discontinuities in the potential and in the forces
soft slabs  by weighing the contributions of each slab n to the total potential function v by a gaussian functioncentered at the midplane of the nth slab . here ,  is the width of the gaussian function , x the distance between adjacent slabs , and a most convenient choice is  = 0.7x andwhich yields a nearly constant sum , essentially independent of xi ( dashed line in figure 3 ) , i.e.,with |(xi)| < 1.3  10 .
this choice also implies that the individual contributions to the force from the slabs add up to unity such that no further normalization is required .
gaussian functions gn centered at nx for a slab distance x = 1.5 nm and n  2 .
gaussian function g0 is highlighted in bold ; the dashed line depicts the sum of the shown gaussian functions . to each slab center xcn ,
all atoms contribute by their gaussian - weighted ( optionally also mass - weighted ) position vectors gn(xi)xi . the instantaneous slab centers xcn are calculated from the current positions xiwhile the reference centers ycn are calculated from the reference positions yi0due to the rapid decay of gn , each slab will essentially involve contributions from atoms located within 3x from the slab center only .
for the first variant , the slab segmentation procedure with gaussian weighting is applied to the radial motion potential ( eq 18/figure 2b ) , yielding as the contribution of slab nand a total potential functionnote that the global center of mass xc used in eq 18 is now replaced by xcn , the center of mass of the slab .
withthe resulting force on atom j readsnote that for v , as defined , the slabs are fixed in space and so are the reference centers ycn . if during the simulation the rotation group moves too far in the v direction , it may enter a region where  due to the lack of nearby reference positions  no reference slab centers are defined , rendering the potential evaluation impossible .
we therefore have included a slightly modified version of this potential that avoids this problem by attaching the midplane of slab n = 0 to the center of mass of the rotation group , yielding slabs that move with the rotation group .
this is achieved by subtracting the center of mass xc of the group from the positionssuch that to simplify the force derivation , and for efficiency reasons , we here assume xc to be constant , and thus xc/x = xc/y = xc/z = 0 .
the resulting force error is small ( on the order of o(1/n ) or o(mj / m ) if mass - weighting is applied ) and can therefore be tolerated . with this assumption ,
the forces f have the same form as eq 37 . in our second variant ,
slab segmentation is applied to v ( eq 24 ) , resulting in a flexible axis potential without radial force contributions ( figure 2c)withthe force on atom j reads applying transformation 38 yields a translation - tolerant version of the flexible 2 potential , v. again , assuming that xc/x , xc/y , and xc/z are small , the resulting equations for v and f are similar to those of v and f.
gromacs 4 distributes the atoms among the parallel processors by domain - decomposing(24 ) the simulation box and assigning each domain to a processor . depending on van der waals and coulomb cutoff settings , positions of atoms near the domain boundaries are communicated such that each processor can compute the forces assigned to its domain .
however , the calculation of some of the proposed potentials and forces requires atom positions not present on the local processor .
for instance , the pivot free potentials require the center of mass of the rotation group , while the flexible potentials require all n positions of the rotation group .
the required coordinates are therefore distributed to all processors before the force calculations , which entails one extra communication step in the rotation module .
further , repeated expressions such as the last terms in eqs 18 and 28 are precalculated whenever possible .
for the efficient computation of the forces f , the inner sum of the last term of eq 37is rewritten assuch that the repeated termsare also precomputed for each relevant slab n and then used for the calculation of each fj term . likewise , for f , the terms s ( eq 45 ) of eq 46 are precalculated .
moreover , for the flexible potentials , only significant contributions to v and f are computed , defined by a cutoff value of gn(x )  gnmin with a default value gnmin = 0.001 , which is checked according to a simple distance criterion . also , the atoms of the rotation group are sorted according to their position along the rotation vector such that for each slab n , a first and a last index i between gn(xi )  gnmin i  [ ifirst ... ilast ] is stored , and all contributions outside that range can safely be ignored .
, the appropriate periodic image for each of the particles of the rotation group has to be chosen such that groups are not split . for fixed axis rotation , each atom
is put closest to its current reference position . for the flexible and pivot - free radial motion potentials ,
each atom is put next to its position at the previous time step , thereby ensuring the integrity of all rotation fragments .
as a sample application of our flexible axis approach , and to compare results obtained by fixed and flexible axis rotation , a series of all - atom md simulations was performed in which the  subunit of f1-atpase was enforced to rotate with respect to its stator part , 33 ( figure 4 ) .
corners , the full protein structure ( 33 ) is shown in a side ( top ) view .
subunit color - coding is  , red ;  , green ;  , cyan ;  , magenta ; and  , orange .
the central panel illustrates the initial orientation of the rotor domain (  ) with respect to the stator ( 33 ) ; for the sake of simplicity , only the  and two  subunits are shown .
the 3-fold symmetry axis of 33 that was used as a rotation axis in v is shown in magenta .
the red spheres and yellow arrows depict slab centers and local rotation axes as used by the flexible potentials .
the left and right side panels show the orientation of the rotor after 120 of enforced rotation using v and v , respectively .
the two orange spheres denote harmonic restraints applied to the n - terminal tags of the  subunits .
this is to prevent co - rotation of the 33 stator in close resemblance to single - molecule force probe experiments , in which the stator is immobilized by attaching the protein to the surface via his tags attached to one subunit type ( usually the  chains ) .
figure prepared with vmd.(38 ) f1-atpase is the soluble domain of the fof1-atp synthase , a rotary motor protein that synthesizes atp from adp using the electrochemical proton gradient across the membrane as its energy source.(26 ) the mitochondrial f1-atpase is an oligomeric protein consisting of nine polypeptide chains , 33.(27 ) in synthesis direction , f1-atpase is driven by the membrane - embedded proton - translocating fo motor while the f1 mobile subunit ,  , rotates clockwise ( seen from the membrane ) within the bearing formed by the hexagonally arranged  and  chains .
the energy transmitted mechanically via the rotating subunit is subsequently used at the catalytic sites of 33 for atp synthesis . to prevent co - rotation ,
the 33 hexamer is connected to the membrane - embedded fo motor by a peripheral linker stalk . despite numerous theoretical and simulation studies ,
the molecular mechanism of energy transmission between the rotor subunit and the ligand binding sites in the stator is still not fully understood.(37 ) the initial configuration of the f1 motor was based on the x - ray structure of bovine f1-atpase determined at 2.4  resolution(39 ) ( protein data bank entry 1e79 ) .
the covalently bound inhibitor as well as the glycerol and sulfate molecules were removed , leaving only mgatp and mgadp ligands in their respective binding sites .
two five - residue - long loops missing from the  subunit were modeled with tconcoord.(40 ) protonation states of ionizable groups were set according to the pka shifts calculated with the delphi(41 ) interface of whatif.(42 ) the protein structure was solvated with 87  321 water molecules in a 16.7  13.8  13.8 nm rectangular unit cell . to neutralize the system and to obtain physiological ionic strength , 261 na and
first , all heavy atoms of the protein and the protein s ligands were kept fixed ; subsequently , all atoms in the system were allowed to relax .
all simulations were performed with gromacs 4.0(24 ) in which the potentials v , v , and v were implemented . for convenience , we here also describe the newer 4.5 version , which produces the same results for the v , v , and v potentials but includes nine additional rotation potentials . for the protein as well as its ligands and ions ,
the opls / aa force field was used , and tip4p(45 ) was used for the water .
all production runs were carried out in the npt ensemble at 300 k and 1 bar .
temperature and pressure were controlled by noshoover ( coupling constant t = 0.5 ps ) and parrinellorahman ( p = 2.0 ps ) schemes , respectively . to avoid severe density oscillations , the first 5 ns of the npt equilibration run were performed with berendsen weak coupling(50 ) for temperature and pressure .
periodic boundary conditions were applied in 3d , and electrostatic forces were calculated with the particle mesh ewald ( pme ) method using a real - space cutoff of 1 nm and an fft grid density of 10 nm .
covalent bond lengths in the protein and ligand were constrained to their reference values with p - lincs.(53 ) settle was used to constrain the water geometry.(54 ) equations of motion were integrated using the leapfrog scheme with a time step of 2 fs . prior to enforcing the rotor movement ,
the system was equilibrated for 10 ns at the target temperature and pressure . during the first 1 ns of this run ,
all protein heavy atoms were harmonically restrained to their initial positions . to mimic the effect exerted on the f1 subunit by the rotation of the fo motor , a potential of the form v ( eq 1 ) , v ( 34 ) , or v ( 40 ) was applied during the production runs .
the longest principal axis of the 33 stator , i.e. , the eigenvector of the inertia tensor of 33 corresponding to the largest eigenvalue , was used as a rotation vector v. for the fixed variant , the pivot vector u of the axis was placed at the center of mass of the 33 units , thus defining the 3-fold pseudosymmetry axis of the stator subunit ( figure 4 ) .
for the flexible axis runs , a slab distance of x = 1 nm , a gaussian function cutoff of gnmin = 0.001 , and  = 0 nm were chosen .
the  reference positions were rotated counter - clockwise around v at an angular rate of  = 0.021/ps over 6 ns of the simulation time , yielding a 120 rotation of the  domain .
due to its symmetry , this covers a complete synthesis cycle , as also seen from the observed stepped motion of the  domain.(28 ) to examine the effect of the chosen spring constant k , for each of the three potentials , five runs were performed with k values ranging from 100 to 800 kj/(molnm ) . in each case , all heavy atoms of the n - terminal six - residue sequences of each  subunit were harmonically restrained to their initial positions using a force constant of 1500 kj/(molnm ) . additionally , for k = 600 kj/(molnm ) ,
a complete 360 rotation of the  domain was simulated . we demonstrate that , due to the optimizations described in the implementation section , simulation performance is nearly unaffected for typical setups such as the f1-atpse system , where only a small fraction of all atoms are subjected to a rotation potential .
the described atpase example with the implemented rotation types was benchmarked ( figure 5 ) on a cluster of intel xeon l5430 nodes connected by a ddr infiniband network .
an intel mpi 3.2.1 was used with the intel 11.1 compiler and the fftw 3.2 library . for the benchmarks , coulomb and van der waals cutoffs
were set to 0.9 nm and the fourier grid to 144  120  120 points , yielding a grid spacing of less than 0.12 nm in each dimension .
separation into long - range ( pme - only ) and short - range ( particleparticle ) processes was allowed .
the optimal number of pme - only processes was derived with the g_tune_pme(55 ) tool using 2000 equilibration steps for the dynamic load balancing , with run times taken from 2000 subsequent steps .
gromacs 4.5 performance for various rotation potentials ( colors ) compared to a simulation without rotation .
the system comprises 401  152 atoms in total , of which 2116 are subjected to the rotation potential . for the flexible potentials , slab distance x = 1 nm and gnmin = 0.001 have been chosen .
the thin black line denotes ideal scaling . with the n = 272 c atoms of the  subunit as the rotation group
, none of the potentials significantly reduced the md performance . to be able to analyze the scaling behavior ( figure 5 ) ,
the rotation group was therefore enlarged to contain all n = 2116 atoms of the  subunit .
as seen , the overall performance decreases only slightly compared to the case without rotation . for the most computationally demanding flexible potentials , on eight processors ,
to verify that the proposed methods properly control the motion of the rotary subunit , we first determined the time evolution of the rotor angle . the actual rotation angle (t ) of the  subunit was determined by a mass - weighted root - mean - square deviation ( rmsd ) fit to the initial (  = 0 ) configuration of the  backbone .
time evolution of the  rotor angle with respect to the 33 symmetry axis for the f1-atpase motor enforced to rotate in the synthesis direction using the potentials v ( a ) , v ( b ) , and v ( c ) with spring constants k of 100800 kj/(molnm ) .
the results show that in all 6-ns - long enforced rotation runs the rotor changes its orientation with respect to the stator by the expected 120. the angle increases nearly linearly with time , with the slope reflecting the constant angular velocity of 0.021/ps , at which the reference is rotated . for fixed axis rotation ,
the subunit closely follows the reference for all tested force constants k = 100800 kj/(molnm ) .
in contrast , for both flexible variants , a less regular evolution is observed , as indicated by the large fluctuations of  for k = 100 kj/(molnm ) . these result from conformational changes of the rotor that occur because the flexible method allows for structural relaxations and adaptations to the bearing .
additionally , at high rotation velocities , frictional forces occur , which cause further conformational changes .
movies illustrating the effect of the fixed and flexible axis methods have been included within the supporting information . in the movies ,
a v and a v rotation potential with k = 600 kj/(molnm ) is applied to all c atoms of the  subunit . for a quantitative comparison of the  subunit internal deformation , figure 7
shows the time evolution of the rmsd of the  backbone atoms from their initial configuration .
relatively small rmsd variations are observed for the fixed method , confirming nearly rigid - body like rotation .
in contrast , both flexible axis methods allow for structural rearrangements particularly for small k values .
a secondary structure analysis shows that for the f1-atpase flexibly rotating at 0.021/ps the force constant k should be 200 kj/(molnm ) or larger to preserve the rotor coiled - coil conformation of the crystal structure ( figure 4 ) . for any of the rotation potentials ,
generally , higher rotation rates will require larger force constants that stabilize the rotation group with the help of a stronger coupling to its reference .
yet , the decrease in conformational freedom with increasing k ( figure 7 ) shows that when using the flexible axis approach one can optimize the tradeoff between structural flexibility and mechanical resistance of the rotary subunit .
note that the 120 rotations , in principle , can not perfectly reproduce the starting configuration of the f1-atpase , as in our simulations the rotor motion is not accompanied by occupancy changes of the active sites .
rmsd of the  subunit backbone atoms with respect to the x - ray structure as a function of time for the f1 motor driven to rotate in the synthesis direction using the potentials v ( a ) , v ( b ) , and v ( c ) with spring constants k of 100800 kj/(molnm ) . because for the flexible potentials the local rotation axis adapts dynamically , it is interesting to monitor the evolution of the f1 rotor angle  also with respect to a variable axis .
figure 8 shows the time dependence of  computed in the same manner as previously but now with the instantaneous ( longest ) principal axis of the  subunit used as the reference axis .
significantly smoother variation of  with time is seen in figure 8 compared to using a fixed symmetry axis ( figure 6b , c ) .
this result illustrates the ability of the flexible methods to adapt the rotation geometry to the structure and conformational changes of the stator .
time evolution of the angular position of the  rotor computed as the best - fit angle with respect to the  longest principal axis for the f1 motor enforced to rotate in the synthesis direction using v ( a ) and v ( b ) .
we will now characterize the different rotation methods in terms of torque and energetics . because the efficiency of the chemomechanical energy transmission in the f1-atpase , when studied in single molecule measurements , is close to 100%,(28 ) and
due to the implied tight coupling between the mechanical reaction coordinate ( e. g. , the  angle ) and conformational changes in the catalytic subunit , the work necessary to enforce a 120 rotation of f1-atpase in the synthesis direction should approach the free energy of 5070 kj / mol required for atp synthesis .
figure 9 therefore compares the torque profiles along the mechanical reaction coordinate  ( eq 53 ) for the three methods considered above .
as can be seen , for both flexible axis potentials , the average torque along  is about 5 times smaller than that for the fixed axis potential . assuming that for infinitely slow rotation the ( equilibrium ) torque curve is smaller than the observed torques , and , further , that this difference is due to dissipation or other nonequilibrium effects , this result implies that the flexible axis approach reduces the dissipated energy by at least a factor of 5 with respect to the fixed axis potential .
the angular dependence of the driving torque for the  subunit enforced to rotate in the synthesis direction using v ( fixed axis ) , v ( flexible axis ) , and v ( flexible axis 2 ) , using five different spring constants k = 100800 kj/(molnm ) .
all torque profiles were smoothed using a running average window of 8. integrated over 120 , the corresponding work is 5900  300 , 1400  100 , and 1490  80 kj / mol , for the fixed , flexible , and flexible 2 potential , respectively . due to the large angular velocity applied as well as the resulting nonequilibrium nature of this process ,
this work is still much larger than the free energy of atp synthesis but clearly shows the dramatic reduction by the flexible axis method . for much lower velocities of 0.00042/ps
, the integrated work reduces further to 350  50 kj / mol ( data not shown ) . already for the short simulations , the dependence of the torque on the angular position of the rotor reveals details of the free energy landscape governing the f1 rotation . in all simulations with a flexible axis
potential ( figures 9 and 10 ) , only small variations of the average torque with respect to the rotor angle are observed .
this suggests that the underlying energy landscape is smooth and nearly linear , which is in agreement with recent experiments .
however , due to the much larger angular velocity employed here , the torque values calculated from our simulations are at least one order of magnitude larger than the f1-generated torque measured under a viscous load ( 4050 pnnm ) . in addition ,
the calculated torque profiles indicate that the free energy landscape is steepest at   40. the increase of the torque for  > 90 that shows up in the fixed type points to shortcomings of this particular method , which will be discussed in the next section .
evolution of the driving torque for the  subunit enforced to rotate in the synthesis direction using v ( red ) and v ( green ) with k = 600 kj/(molnm ) ( a ) .
rmsd of the  rotor backbone atoms ( solid ) and of the 33 stator backbone atoms ( dotted ) with respect to their respective x - ray structure ( b ) .
since in the simulations the motion along the mechanical reaction coordinate is not synchronized with changes in chemical occupancy of the stator binding sites , we do not expect the torque to drop to zero after 120 of rotation . to examine
how far our nonequilibrium simulations are from the reversible limit , figure 10a shows torque profiles of full 360 rotations applying v and v with k = 600 kj/(molnm ) . as can be seen , the torque determined for the fixed axis case increases strongly , whereas for the flexible case , after a small increase up to   40 , the torques decrease toward considerably smaller values .
this result underscores that the flexible potential perturbs the system to a much lesser extent , such that it remains much closer to equilibrium than for a fixed axis rotation .
the fixed axis potential induces structural changes almost exclusively in the bearing ( stator in figure 10b ) while in the flexible axis case , the structural changes are distributed rather equally among the rotating subunit and its bearing .
moreover , in the flexible axis case , both structures nearly approach the starting structure (  = 0 ) at the end of a whole 360 turn with an rmsd below 2.5  , which is not seen for the fixed axis simulation .
when using a simple fixed axis rotation potential , the rotating part behaves like a rigid body . in combination with the fixed axis
, this behavior can cause unphysical close contacts and strong torques between the rotor and the bearing , which may cause extensive artificial structural changes of the bearing .
the flexible axis approach , in contrast , keeps the system closer to the equilibrium for two reasons .
first , the self - adjusting local rotation axis ensures an overall optimal position of the pivot ; second , the built - in flexibility allows for structural relaxation of the rotating part and thus locally minimizes sterical hindrances . as the f1-atpase motor components are strongly coupled and leave only little room for the rotating subunit inside the bearing , both reasons allow for the necessary tight adaption of the  rotor to the 33 bearing . to quantify this effect , figure
11 displays the enforced conformational changes of the bearing , in terms of stator rmsd with respect to its x - ray structure as a function of time for the 2 subunit .
this subunit was chosen because it interacts most closely with the  rotor throughout the whole runs .
it is evident that the structural changes induced in 2 are considerably larger for fixed axis rotation than for the flexible potentials .
secondary structure analysis reveals that in the former case the structural motifs exposed to the center of the 33 hexamer are distorted by the rotating  subunit .
also in figure 4 one can notice partial disruption of the helices in the c - terminal part of the 3 subunit ( the bottom part of  on the left side of  , red dashed circle ) when it is pressed upon by the rotor driven to rotate around the fixed axis .
the torque increase for angles  > 90 ( figure 9 , left , and figure 10 ) reflects this effect , which is mainly due to wrapping of the 3 c - terminal domain around the  subunit .
rmsd of the 3 subunit backbone atoms with respect to the x - ray structure when driving  subunit rotation using v ( a ) and v ( b ) . for comparison ,
the corresponding rmsd evolution for five independent free md runs is also shown ( c ) .
we have developed , implemented , and tested a new method to enforce the rotation of protein subunits that allows for ( i ) a flexible rotation axis and ( ii ) structural adaptions of the rotated subunit to its environment . for  subunit rotation in f1-atpase , we have shown that our flexible axis method reduced the frictional dissipation of the  subunit within the 33 bearing by more than a factor of 5 . as a result ,
also the induced torque was 5-fold smaller compared to the one using a fixed axis . concerning the use of the flexible axis potentials developed here
the first caveat is due to the fact that , while the pivot vector is free to adapt flexibly , the orientation of the direction vector is fixed . for systems where the subsystem subjected to the rotation potential is embedded within a curved  bearing  , the flexible adaptation will work properly only as long as the angle between the orientations of the bearing axis and the direction vector is not too large , i.e. , for not too strong bending of the bearing .
in extreme cases such as a complete u - shaped bearing , artificial structural changes of the bearing similar to those induced by fixed axes may occur .
this problem can be addressed by subdividing the system into several parts and using a separate flexible rotation axis for each of these parts , with orientation vectors locally adapted to the respective part of the bearing .
if chosen too small , only a few atoms will be assigned to each slab , thus compromising the averaging that defines the pivot vector of each slab . in contrast , if chosen too large , the slabs might stretch over regions that would require changing pivot vectors , in which case the enforced rotation would induce , albeit to a lesser extent , the artifacts caused by fixed axis approaches . obviously , in the limit of just a single slab covering the complete rotating subsystem , the fixed axis potential is recovered . with these limitations and caveats in mind , our flexible axis potentials are applicable to a broad range of quite diverse biomolecular systems , processes , and functions . apart from mimicking molecular rotary motors , it can also serve to restrain the orientation of a protein or ligand , or , in combination with umbrella sampling , to calculate the preferred orientation of transmembrane proteins or membrane - active agents within a lipid bilayer . further
, the method is expected to yield more accurate free energy profiles along circular reaction coordinates via umbrella sampling . in the long run , our flexible axis approach might prove useful for the study and design of synthetic nanodevices with rotating elements , such as those considered for molecular nanotechnology.(57 )
 we describe a versatile method to enforce the rotation of subsets of atoms , e.g. , a protein subunit , in molecular dynamics ( md ) simulations . 
 in particular , we introduce a  flexible axis  technique that allows realistic flexible adaptions of both the rotary subunit as well as the local rotation axis during the simulation . a variety of useful rotation potentials were implemented for the gromacs 4.5 md package . 
 application to the molecular motor f1-atp synthase demonstrates the advantages of the flexible axis approach over the established fixed axis rotation technique . 
Introduction
Fixed Axis Rotation
Flexible Axis Rotation
GROMACS Implementation
Application to F
Results
Conclusions