# 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.

**Introduction**

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

- the F
_{o}and F_{1}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 F_{1}-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.

**Methods**

To enforce rotation, a group of* N* atoms (the »rotation group«) is subjected to a rotation potential* V*. Each atom with position* x_{i}* gets assigned a reference (equilibrium) position

*, which is rotating at a constant angular rate ω about an axis*

**y**_{i}(t)**. Typically the initial positions of the rotation group provide the initial (**

*v**t=0*) set of reference positions

*. The »isotropic« potential (Fig. 2 B, 3 A)*

**y**_{i}^{0}\[ 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

*, thus effectively rotating*

**y**_{i}(t)*. With the prefactors*

**x**_{i}\[ 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 F_{1}-ATPase γ subunit rotating inside the α_{3}β_{3} bearing, applying *V ^{iso}* on the 272 C

_{α}carbons of the γ subunit. Starting from the fixed axis potential

*V*we will in 4 steps develop the flexible axis potential. Details can be found in the publication.

^{iso}**Step 1: Detach the pivot**

A translation-invariant pivot is achieved by selecting the group’s center of mass * x_{c}* 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 *V ^{iso}* and

*V*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:

^{iso-pf}\[ 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 *V ^{rm-pf}* still contains a small radial component, which is not the case for

*V*. Note that choosing a small positive ε’ yields a well-defined potential also in the vicinity of the pivot, compare Fig. 3C with 3D.

^{rm2-pf}**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 *V ^{rm-pf}*, yielding the »flexible« potential

\[ V^\text{flex} =

\frac{k}{2} \sum_{i=1}^{N} \sum_n w_i \, g_n(\mathbf{x}_i)

\left[

\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) \| }

\cdot

(\mathbf{x}_i - \mathbf{x}_c^n)

\right]^2 \]

or, alternatively, to *V ^{rm2-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 +

\epsilon'}\]

The flexible axis rotation movie at the bottom of this page shows the F1-ATPase γ subunit rotating inside the α_{3}β_{3} bearing, applying *V ^{flex2}* 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 *V*^{iso}, 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 `

; 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**