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₁ 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₁-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 y_i(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 y_i⁰. 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 y_i(t), thus effectively rotating x_i. 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 F₁-ATPase γ subunit rotating inside the α₃β₃ bearing, applying V^iso on the 272 C_α carbons of the γ subunit. Starting from the fixed axis potential V^iso 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 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^iso-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 V^rm-pf still contains a small radial component, which is not the case for V^rm2-pf. Note that choosing a small positive ε’ yields a well-defined potential also in the vicinity of the pivot, compare Fig. 3C with 3D.

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 α₃β₃ 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.