Subharmonic oscillations of collective molecular motors

Molecular motors are force generators which can produce mechanical oscillations at the cellular scale and play a central role in muscle contraction and cell motility [1,2]. Insect fibrilar muscles of wasps and bees are known to exhibit a wing thrust which is asynchronous to the activating nervous impulses [3,4]. Skinned skeletal and cardiac muscle fibers have also been shown to exhibit spontaneous oscillations in vitro under various conditions [5,6]. Other nonmuscular examples are the periodic motion of cilia and flagella used for cell motility [7,8] and the amplification of weak stimuli via oscillatory instabilities of mechanosensory hair bundles [9,10]. Early models proposed to understand such systems were grounded on experimental evidences of muscle contraction [11–13] and were subsequently generalized to the study of eukaryotic flagella [14,15]. Two main approaches have been used to study instabilities in motor assemblies: two-state "crossbridge" and "stiff motor" models, which can be considered as two different limits of a more general description that considers the motor-filament interaction and the stiffness of the motors [16,17]. In both cases, oscillations are obtained for nonmonotonic velocity-force relationships in the presence of an elastic element. In this work, we study a generic model for the collective dynamics of molecular motors which generates spontaneous oscillations via a Hopf bifurcation [18,19], and successfully explained many experimental results on biological oscillating systems driven by molecular motors [10,20–22]. We show that the pair of integro-partial differential equations studied in ref. [18] can be reduced to a set of three ordinary differential equations which capture the main dynamics of the system and exhibit rich nonlinear behaviour. The robustness of this reduced description is validated through direct comparison with the full system of integro-partial differential equations. Finally, we compare our results with experimental oscillations of a minimal actomyosin system in vitro [22] and show that the unexplained secondary peaks observed in the power spectrum can be understood in the limit of high mobility and large elastic modulus as subharmonics resulting from nonlinearities in the system.

We consider a collection of N molecular motors which are rigidly attached to a common backbone filament as sketched in fig. 1. Motor heads can be found in two possible conformational states $\sigma=1,2$ subject to a l-periodic potential $W_{\sigma}(x)$. Depending on whether they are bound to or unbound from the filament, motors feel the potentials W₁,W₂, respectively. The backbone is coupled to a spring of elastic modulus KN, where K is the elastic modulus per motor. Motors can detach or attach with l-periodic rates $\omega_{1}(x)$, $\omega_{2}(x)$, respectively. In order to obtain spontaneous oscillations, detailed balance must be broken through the addition of localized "active sites" in the minima of the potential W₁. We can account for this effect by introducing an l-periodic function $\theta(x)$ which is proportional to the adenosine triphosphate (ATP) hydrolysis rate:

$$\begin{equation} \frac{\omega_1}{\omega_2}=\exp \left(\frac{\Delta W}{k_BT}\right) + \theta (x), \end{equation} \tag{ 1 }$$

where $\Delta W (x) \equiv W_1(x)-W_2(x)$. The first term on the right side of eq. (1) accounts for thermally driven excitations while the second term describes ATP-hydrolysis–driven excitations. Motors are spaced uniformly with distance q incommensurate to l. We define the probability density of finding a motor in state σ as $P_{\sigma}(x,t)$ and the backbone position as X(t). For simplicity, we will work in the limit of a large number of motors, where we can assume the total motor density distribution to be uniform and therefore $\sum_{\sigma} P_\sigma(x,t)=1/l$. Switching to the cyclic coordinate $\xi = x \, \text{mod} \, l$, the dynamics of the system reads [18]

$$\begin{eqnarray} \partial_t P_1 +v \partial_{\xi}P_1&=-\alpha P_1+\omega_2/l, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \lambda v &= - KX +F, \end{eqnarray} \tag{ 3 }$$

where $v \equiv \dot{X}$ is the backbone velocity, λ is the friction coefficient of the medium, $\alpha(\xi) \equiv \omega_1(\xi)+\omega_2(\xi)$ and $F= -\int_{0}^{l}\text{d}\xi P_1 \partial_{\xi} \Delta W$ is the active force generated per motor. Since P₁ and P₂ are constrained, we focus on the dynamics of one probability distribution (e.g. P₁). Equation (2) describes the dynamics of the probability density P₁ while eq. (3) accounts for force balance in the system. Following the work in ref. [23], we decompose $P_1(\xi,t)$ as an infinite sum of modes $a_n(t)$ and $b_n(t)$:

$$\begin{equation} P_1(\xi,t)=\sum_{n=0}^{\infty} [a_{n} \sin(2 \pi n \xi/l)+ b_n\cos(2 \pi n\xi/l)]. \end{equation} \tag{ 4 }$$

Since the functions $W_1(\xi),\omega_1(\xi), \omega_2(\xi)$ are periodic, they can be generally described in terms of Fourier series. We will consider the simplest approximation by keeping only the first Fourier mode and defining $g(\xi)=1+\cos(2 \pi \xi/l)$. We choose a symmetric sinusoidal potential $W_1(\xi)=(U/2)g(\xi+l/2)$, where U is the amplitude of the potential. We also take α to be constant, such that the sum of the transition rates is independent of the spatial coordinate ξ. This approximation is convenient when the transition rates are almost uniform. The transition rates are given by

$$\begin{equation} \begin{array}{rcl} \omega_1(\xi)&=& (\beta/2) g(\xi),\\ \omega_2(\xi)&=&\alpha - \omega_1(\xi), \end{array} \end{equation} \tag{ 5 }$$

where $\alpha, \beta$ are unknown rates. Without loss of generality, we assume W₂ to be constant in the detached state. Since biological motors work far from equilibrium, $\omega_1/\omega_2 \simeq \theta(\xi)$, and thus

$$\begin{equation} \frac{\omega_1}{\omega_2} \simeq \frac{\Omega g(\xi)}{1-\Omega g(\xi)}, \end{equation} \tag{ 6 }$$

where $\Omega \equiv \beta/(2\alpha)$ is a dimensionless ATP hydrolysis amplitude and $\Omega \in [0,1/2)$. Using the orthogonality of the Fourier modes, the active force per motor is simply given by $F=-U \pi a_1/2$. Hence, we obtain an infinite set of ordinary differential equations for the modes $a_n(t),b_n(t)$, $n \geq 0$ and X(t). Remarkably, for constant α, the modes $a_1(t),b_1(t)$ and X(t) are formally decoupled, that is, they evolve independently from the rest (b₀ and a_n, b_n for n > 1). The last choice of $\alpha(\xi)$ is known as the constant rate approximation [23]. The evolution of b₀ is also decoupled from the rest and relaxes exponentially to its steady-state value. Nondimensionalizing with respect to the length scale l, time scale $1/\alpha$ and force density $\lambda \alpha$, the dimensionless governing equations read

$$\begin{eqnarray} \begin{array}{rcl} \dot{\bar{a}}_1&=& - \bar{a}_1 - 2 \pi \rho \bar{X} \bar{b}_1 - \mu \bar{a}_1\bar{b}_1,\\ \dot{\bar{b}}_1&=& -\Omega- \bar{b}_1 + 2 \pi \rho \bar{X} \bar{a}_1 + \mu \bar{a}_1^2,\\ \dot{\bar{X}}&=&-\displaystyle\rho \bar{X} - \frac{\mu}{2\pi} \bar{a}_1 \end{array} \end{eqnarray} \tag{ 7 }$$

and for n > 1

$$\begin{eqnarray} \begin{array}{rcl} \dot{\bar{a}}_n&=- \bar{a}_n - 2\pi n \rho \bar{X} \bar{b}_n - n \mu \bar{a}_1 \bar{b}_n,\\ \dot{\bar{b}}_n&=-\bar{b}_n + 2 \pi n \rho \bar{X} \bar{a}_1 + n \mu \bar{a}_1\bar{a}_n, \end{array} \end{eqnarray} \tag{ 8 }$$

where Ω plays the role of the control parameter, $\rho \equiv K/ \lambda \alpha$ sets the damping rate of oscillations and $\mu \equiv \pi^2 U/\lambda l^2 \alpha$ is an effective mobility. It is interesting to remark that a similar decoupling occurs for the Markus-Lorenz waterwheel equations [24,25] which share some similarities with eqs. (2) and (3). It is noteworthy that all nonlinear terms depend inversely on the friction λ, which plays a major role on the nonlinear nature of the system. Defining the dimensionless active force per motor as $f \equiv F/\lambda \alpha l$ we notice that f is proportional to the first asymmetric mode through $f=- \mu \bar{a}_1/2 \pi$. The system has the symmetry $(\bar{X},-f) \rightarrow (-\bar{X},f)$, given that $W_1(\xi)$ is symmetric. Furthermore, the backbone position $\bar{X}$ far from the initial value is uniquely given by the evolution of f, from (7):

$$\begin{equation} \bar{X}(\tau)= \int_0^{\infty} \text{d} \tau^{\prime} e^{-\rho \tau^{\prime}} f(\tau-\tau^{\prime}), \end{equation} \tag{ 9 }$$

where $\tau^{\prime}$ is the time-lag. Hence, $\bar{X}$ represents the memory of f, i.e. its exponentially weighted past evolution. Through a suitable linear change of variables $\{2\pi \rho \bar{X} + \mu \bar{a}_1, \bar{a}_1, \bar{b}_1+\Omega\} \rightarrow \{x,y,z\}$ the set of equations (7) can be rewritten as

$$\begin{eqnarray} \begin{array}{rcl} \dot{x}&=& -\rho x + \mu \dot{y},\\ \dot{y}&=& x(\Omega-z)-y,\\ \dot{z}&=&-z+xy, \end{array} \end{eqnarray} \tag{ 10 }$$

where x describes the position of an overdamped particle in an effective system under a forcing term $\mu \dot{y}$. Now the previous symmetry condition is manifested as $(x,y) \rightarrow (-x,-y)$. The dynamics of y and z is the same as in the Lorenz model (for the particular choice of $\beta=1$ in ref. [26]); however the first equation in (10) is different. The system has a single fixed point at the origin. A linear stability analysis around this point gives one negative eigenvalue −1 and two complex conjugated eigenvalues $\gamma \pm i \omega$, where $\gamma = \mu(\Omega-\Omega_c)/2$, $\omega=\sqrt{\rho-\gamma^2}$. $\Omega_c=(1+\rho)/\mu$ is defined as the critical value of Ω where a Hopf bifurcation takes place, i.e. when $\gamma=0$. In view of the similarities between (10) and the Lorenz equations, we investigated the possibility of chaos in the system. We considered the attractor shown in fig. 2, where nearby trajectories tend to locally diverge when found in the two spirals and to locally converge when switching between spirals. The computation of the Lyapunov characteristic exponents for the attractor in fig. 2 was studied numerically by using a Mathematica package [27] based on the algorithms presented in refs. [28,29]. The maximum Lyapunov exponent $\lambda_1$ was found to be nearly zero and the remainder exponents were $\lambda_2 \gtrsim \lambda_3 \simeq -0.02$, indicating that the limit set is a periodic orbit and consistently with the system being dissipative. In conclusion, no clear signs of chaotic behaviour were found in the system. Despite the similarities between (10) and the Lorenz equations, the different terms in the first equation of (10) crucially affect the dynamics of the system. Interestingly, as we will discuss next, the system can generally exhibit subharmonic oscillations. Hereinafter, we will focus on the dynamics of $\bar{X}$ rather than x to allow a physical interpretation compatible with empirical observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We choose the model parameters consistent with the minimal actomyosin system in ref. [22] which was shown to be successfully described by the two-state model in ref. [18]. Typical values of KN are found in the range of $10^{-2}\text{-}10^{-1}\ \text{pN/nm}$ and the number of motors are $N \sim 10\text{-}100$, which gives $K \simeq 10^{-4}\text{-}10^{-2}\ \text{pN/nm}$, while the friction coefficient is found to be $\lambda \simeq 10^{-2}\text{-}1\ \mu \text{N}\cdot \text{s}\cdot \text{m}^{-1}$ [18,22,30]. We take $\alpha=10\ \text{s}^{-1}$ and choose Ω much larger than the critical value. By considering $l = 6\ \text{nm}$, $U = 10 k_BT$ we find $\rho \simeq 10^{-2}\text{-}10^{2}$ and $\mu \simeq 10^3\text{-}10^5$. Therefore, while the dimensionless damping rate ρ is found to be around the characteristic rate of the system (i.e. $\sim 1$), the mobility μ is specially high, thus the system works in a regime where nonlinearities are particularly important. In fig. 3 we study two characteristic types of nonlinear oscillations by solving the reduced system (7) (top row) and the complete system (eqs. (2) and (3)) (bottom row) for two different set of parameters corresponding to small (fig. 3(a)) and large (fig. 3(b)) elastic modulus. Numerical solutions for the complete system are carried out taking $W_1(\xi)$ as a symmetric saw-tooth potential of amplitude U, a constant value for W₂, a sinusoidal form of the hydrolysis amplitude modulation $\theta(\xi)=(\Omega/2) g(\xi)$ and a constant attachment rate $\omega_2=\alpha$. The case (a) corresponds to $\sim 1\ \text{Hz}$ cusp-like oscillations with peak-to-peak amplitude of $\sim 100\ \text{nm}$. On the other hand, case (b) corresponds to $\sim 10\ \text{Hz}$ oscillations with peak-to-peak amplitude of $\sim 10\ \text{nm}$. Notice that despite of using a different set of functions to solve the reduced and the complete system, the shape of the oscillations is very similar. Both cases qualitatively agree with the observed amplitude and frequency measurements reported for spontaneous oscillations in ref. [22]. In fig. 4, subharmonic oscillations are shown by solving (7) (fig. 4(a), (b), (c) and (d) (left)) and the complete system (fig. 4(d) (right)). The periodic movement of the backbone is characterized by the formation of one (fig. 4(c) (left)) or two (fig. 4(c) (right)) subharmonics, where $\omega_0$ is the fundamental frequency of the signal. Remarkably, similar results where found experimentally in ref. [22] for large optical trap stiffness where the fundamental frequency was 2.2 Hz and a clear subharmonic peak at approximately $2\omega_0/3$ was observed. Although it was argued that molecular details of the actomyosin interaction could be the reason for this effect, we find that for high mobility values μ and in the limit of large elastic modulus (i.e. $\rho \sim 1\text{-}10$), the system is expected generally to exhibit periodic motion with two subharmonics and less frequently with only one, as shown in figs. 5 and 6. We identify eight types of staircase-shaped oscillations, which exemplify the complex bifurcation scenario. The final steady state of the system is sensitive to the initial conditions as shown in fig. 6, where different basins of attraction are identified. Two types of oscillations exhibit two subharmonics (dark-blue and grey diamonds) and one type exhibits a single subharmonic (light-blue diamonds). It is worth noting that subharmonic oscillations lose the symmetry property $\bar{X}(\tau+\mathcal{T}/2) =-\bar{X}(\tau)$, where $\mathcal{T}$ is the fundamental period of the signal. This property is fulfilled for oscillations in fig. 3 since W₁ is symmetric; however this can be lost for asymmetric W₁ [18]. In this case, a novel symmetry breaking property emerges dynamically due to subharmonic bifurcations of the system. Additionally, the time average of the backbone position in the steady state $\langle \bar{X} \rangle$ is different from zero for subharmonic oscillations and introduces an overall shift, as shown in figs. 4(a) and (d). The latter reflects asymmetries on the time the backbone spends for positive and negative displacements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work we have studied in detail a minimal three-variable description for the spontaneous oscillations of collective molecular motors, based on a generic two-state model coupled to an elastic spring. We find that this reduced description captures the dynamics of the full system, which, in turn, unveils different types of nonlinear oscillations which can generally exhibit one or two subharmonics in the limit of large elastic modulus and high mobility. Although the reduced description is only exact for a specific set of functions, direct numerical simulations indicate that the dynamics captured by the reduced system is generic and thus shared by other sets of functions, in particular concerning the occurrence of the referred subharmonic bifurcations. This suggests that a general three-dimensional reduction is inherent to the system, although in general, the corresponding set of variables may not coincide with the explicit ones used in our case. In addition, we found that the importance of nonlinearities is enhanced for smaller friction, a regime that is relevant for myosin motors in vitro. The emergence of two subharmonics on the spectrum of the backbone signal for large elastic modulus is in remarkable accordance with similar subharmonic peaks observed experimentally in analogous conditions [22]. It is worth stressing that the experiments are subject to strong noise sources, mainly due to the stochastic nature of myosin binding kinetics. This fact is likely to modify in a nontrivial way the amplitude and shape of the subharmonic peaks, by inducing transitions between the different oscillatory regimes, and presumably decreasing the time spent in the subharmonic oscillations in favor of the fundamental frequency oscillations. Furthermore, the studied model is not expected to yield an accurate description of the physical system at molecular level, since the precise forms of the potentials and transition rates are essentially unknown. Nevertheless, the mere fact that the subharmonic peaks are observable under the experimental conditions constitutes by itself a strong evidence of the robustness of this phenomenon. Other types of complex nonlinear behavior of molecular motor assemblies have also been reported in the literature for spontaneous sarcomere dynamics [31]. In that case, the assumptions of the model accounted for different physical ingredients and led to notably different nonlinear dynamics, including excitable behaviour and a Ruelle-Takens route to chaos.

The results presented here suggest that hydrodynamic friction could be tuned in biological systems in order to suppress undesired multifrequency oscillations. For example, mechanosensory hair bundles are self-tuned close to a Hopf bifurcation in order to obtain a nonlinear amplification gain [10]. Thus, it is essential for hair cells to have a characteristic frequency at which they are most sensitive. Likewise, flagellated organisms, such as spermetozoa, rely on the propagation of harmonic bending waves along the flagellum to achieve locomotion [32], and thus high motor mobility could slow down their motion and lead to uncoordinated bending waves. Finally, we hope that our results will stimulate future empirical tests on minimal in vitro systems to further scrutinize the role of nonlinearities in molecular-motor–driven oscillations.

We acknowledge the support under Projects No. FIS2010-21924-C02-02 and No. 2009-SGR-014. DO acknowledges a FPU grant from de Spanish Government. HG is an Oxford University Hooke Fellowship Holder and is supported by Award KUK-C1-013-04 from King Abdulah University of Science and Technology.