Studying Rotary Motions in Proteins
with GROMACS: The »Flexible Axis« Approach

We describe a novel method to enforce rotation of a protein subunit in molecular dynamics (MD) simulations. Our »flexible axis« approach allows flexible adaptions of both the rotating subunit as well as the rotation axis during the simulation. For the example of F1-ATP synthase (Fig. 1) we show that the flexible method (apart from the rotation itself) imposes minimal constraints on the rotation group, and allows for conformational adaptions to the surrounding.


Rotary motion of a subunit is in many motor proteins an essential part of their function. Examples are

    • the Fo and F1 motors in F-ATPase,
    • the bacterial flagellar motor,
    • DNA helicases, and
    • DNA translocation into viral capsids.

    Such processes are typically too slow and/or too infrequent to be observed on molecular dynamics (MD) time scales. However, rotation can be induced and/or increased in rate with the help of external forces that are exerted on certain subunits. Previous studies adopted a fixed and stiff axis to enforce rotation of a structural part (Fig. 2AB). This approach, however, does not reflect situations such as F1-ATP synthase (Fig. 1), where the rotating part is flexible and adapts to the steric restraints of its bearing (green), see also Fig. 2C. To more realistically describe biomolecular rotations, we propose a »flexible axis« technique with that a protein part can be rotated also within a tight surrounding—like a pipe-cleaner inside a pipe.


    To enforce rotation, a group of N atoms (the »rotation group«) is subjected to a rotation potential V. Each atom with position xi gets assigned a reference (equilibrium) position yi(t), which is rotating at a constant angular rate ω about an axis v. Typically the initial positions of the rotation group provide the initial (t=0) set of reference positions yi0. The »isotropic« potential (Fig. 2 B, 3 A)

     \[  V^\text{iso} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t) ({\mathbf y}_i^0 - {\mathbf u}) - ({\mathbf x}_i - {\mathbf u})  \right]^2  \]

    with the rotation matrix Ω and the pivot u of the axis yields forces towards the (rotating) reference positions yi(t), thus effectively rotating xi. With the prefactors

    \[ w_i = N \, m_i/M \] with total mass \[ M = \sum_{i=1}^N m_i \]

    mass-weighting can be achieved. The fixed axis rotation movie (see download link on the right) shows the F1-ATPase γ subunit rotating inside the α3β3 bearing, applying Viso on the 272 Cα carbons of the γ subunit. Starting from the fixed axis potential Viso we will in 4 steps develop the flexible axis potential. Details can be found in the publication.

    Step 1: Detach the pivot

    A translation-invariant pivot is achieved by selecting the group’s center of mass xc as the pivot

    \[ \mathbf{x}_c   = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{x}_i  \] and \[ \mathbf{y}_c^0 = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{y}_i^0 \]

    resulting in the »pivot-free« potential

    \[ V^\text{iso-pf} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
    ({\mathbf y}_i^0 - {\mathbf y}_c^0) - ({\mathbf
    x}_i - {\mathbf x}_c) \right]^2 \]

    Step 2: Apply only angular restraints

    For Viso and Viso-pf the potential minimum is a single point at the position of the reference, thereby penalizing radial deviations from the reference conformation. The »radial motion« potential (Fig. 3B) does not constrain an atom if it already has the correct angular position:

    \[ V^\text{rm-pf} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[
    \frac{\hat{\mathbf{v}}\times \mathbf{\Omega}(t) ({\mathbf y}_i^0 - {\mathbf
    y}_c^0)} {\| \hat{\mathbf{v}}\times \mathbf{\Omega}(t) ({\mathbf y}_i^0 -
    {\mathbf y}_c^0)\|}
    \cdot({\mathbf x}_i - {\mathbf x}_c)
    \right]^2 \]

    or, alternatively (Fig. 3C),

    \[ V^\text{rm2-pf} =
    \frac{k}{2} \sum_{i=1}^{N} w_i\,
    \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c ))
    \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c) \right]^2}
    {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c) \|^2 +
    \epsilon'} \]

    since Vrm-pf still contains a small radial component, which is not the case for Vrm2-pf. Note that choosing a small positive ε’ yields a well-defined potential also in the vicinity of the pivot, compare Fig. 3C with 3D.

    Step 3: Segment into slabs

    Now the rotation group is subdivided into a set of equidistant slabs perpendicular to v (Fig. 2C). To avoid discontinuities in the potential, »soft slabs« are used by weighing the contributions to V from each slab by a Gaussian function (see Fig. 4)

    \[ g_n(\mathbf{x}_i) = \Gamma \ \mbox{exp} \left(
    -\frac{\beta_n^2(\mathbf{x}_i)}{2\sigma^2}  \right) \]

    centered at the midplane of slab n. σ is the width of the Gaussian, Δx the slab distance and \[ \beta_n(\mathbf{x}_i) := \mathbf{x}_i \cdot \hat{\mathbf{v}} - n \, \Delta x \]

    We define the instantaneous centers of the slabs (= local slab axis pivots, red in Fig. 1) as

    \[ \mathbf{x}_c^n =
    \frac{\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i \, \mathbf{x}_i}{\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i} \]

    and the reference slab centers as

    \[ \mathbf{y}_c^n =
    \frac{\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i \, \mathbf{y}_i^0}{\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i} \]

    Step 4: Sum up the contributions of all slabs

    For a smooth potential with continuous forces, Gaussian weighting is applied to Vrm-pf, yielding the »flexible« potential

    \[ V^\text{flex} =
    \frac{k}{2} \sum_{i=1}^{N} \sum_n w_i \, g_n(\mathbf{x}_i)
    \frac{\hat{\mathbf{v}} \times
    \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) }{ \| \hat{\mathbf{v}}
    \times \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \| }
     (\mathbf{x}_i - \mathbf{x}_c^n)
    \right]^2 \]

    or, alternatively, to Vrm2-pf, yielding the »flexible 2« potential

    \[ V^\text{flex2} =
    \frac{k}{2} \sum_{i=1}^{N} \sum_n w_i\,g_n(\mathbf{x}_i)
    \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c^n ))
    \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \right]^2}
    {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c^n) \|^2 +

    The flexible axis rotation movie at the bottom of this page shows the F1-ATPase γ subunit rotating inside the α3β3 bearing, applying Vflex2 on the 272 Cα carbons of the γ subunit.

    Results and Discussion

    We have implemented a versatile method to enforce rotation upon protein subunits for GROMACS 4.6. There are ten different rotation potentials available, of which the flexible ones impose the least constraints on the rotated subunit (see Fig. 5). With Viso, the rotated structure always stays close to its reference conformation, while for the flexible types the RMSD to the reference is about three times as large. Flexible rotation improves in several apects on the usually adopted fixed axis approach:

      • A curved axis fits the shape of a curved protein.
      • The impact of an arbitrary choice of the pivot is avoided.
      • The conformation of the rotation group is not artificially stabilized at the reference conformation.
      • The axis curvature continuously adapts to structural changes.

      The flexible axis approach makes it possible to simulate various rotation-related phenomena on MD time scales with minimal restrictions on the rotated parts.

      Exemplary GROMACS .mdp file entries for enforced rotation

      ; Enforced rotation: No or Yes
      rotation                 = Yes
      ; Output frequency for angle, torque and rotation potential energy for the whole group
      rot-nstrout              = 1
      ; Output frequency for per-slab data (angles, torques and slab centers)
      rot-nstsout              = 10
      ; Number of rotation groups
      rot-ngroups              = 1
      ; Rotation group name  
      rot-group0               = System
      ; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t
      rot-type0                = flex2-t
      ; Use mass-weighting of the rotation group positions
      rot-massw0               = yes
      ; Rotation vector, will get normalized
      rot-vec0                 = 1 0 0
      ; Pivot point for the potentials iso, pm, rm, and rm2 [nm]
      rot-pivot0               = 4.31852e+00  1.73201e+00  1.89800e+00
      ; Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)]
      rot-rate0                = 10.0
      rot-k0                   = 500.0
      ; Slab distance for flexible axis rotation [nm]
      rot-slab-dist0           = 1.5
      ; Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)
      rot-min-gauss0           = 0.001
      ; Value of additive constant epsilon' [nm^2] for rm2* and flex2* potentials
      rot-eps0                 = 0.0001
      ; Fitting method to determine angle of rotation group (rmsd, norm, or potential)
      rot-fit-method0          = norm
      ; For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated
      rot-potfit-nsteps0       = 21
      ; For fit type 'potential', distance in degrees between two consecutive angles
      rot-potfit-step0         = 0.25

      Go to Editor View