Internal friction controls active ciliary oscillations near the instability threshold

TL;DR

This study demonstrates that internal friction, rather than external fluid resistance, is the primary dissipation mechanism driving ciliary oscillations near the instability threshold, revealing a new understanding of active filament dynamics.

Contribution

It shows that internal elastic stress dominates dissipation in ciliary beating, contrasting with previous assumptions of external fluid friction being primary.

Findings

Internal friction surpasses external fluid friction in ciliary oscillations.

Active stresses leading to oscillations are strain softening and shear thinning.

A combined experimental and theoretical approach elucidates ciliary instability mechanisms.

Abstract

Ciliary oscillations driven by molecular motors cause fluid motion at micron scale. Stable oscillations require a substantial source of dissipation to balance the energy input of motors. Conventionally, it stems from external fluid. We show, in contrast, that external fluid friction is negligible compared to internal elastic stress through a simultaneous measurement of motion and flow field of an isolated and active Chlamydomonas cilium beating near the instability threshold. Consequently, internal friction emerges as the sole source of dissipation for ciliary oscillations. We combine these experimental insights with theoretical modeling of active filaments to show that an instability to oscillations takes place when active stresses are strain softening and shear thinning. Together, our results reveal a counterintuitive mechanism of ciliary beating and provide a general experimental and theoretical methodology to analyze other active filaments, both biological and synthetic ones.

Tables

Table 1. Table S1 : Possible values of shear friction coefficient from the literature. Shear friction coefficient from existing literature and reasons for neglecting/accepting them.

Reference	${\Gamma }_u$ [$\times {10}^{-6}$]	Comparison of shear friction with elasticity
	Ns/m
Brokaw Brokaw (1972)	0.06-0.51	negligible compared to elastic terms
Murase-Shimizu Murase and Shimizu (1986)	30	very high compared to elastic terms,
		overdamping the system
Bayly-Wilson Bayly and Wilson (2014)	0.5	negligible compared to elastic terms
Bayly-Dutcher Bayly and Dutcher (2016)	0.16	negligible compared to elastic terms
Ingmar Riedel Riedel (2005)	10	comparable and competing with elastic terms, ✓
Minoura-Yagi-	365	very high compared to elastic terms,
Kamiya Minoura et al. (2001)		overdamping the system*

Equations

$${\partial }_s\mathbf{F}+\mathbf{f}=0,{\partial }_s\mathbf{M}+\mathbf{t}\times \mathbf{F}+\mathbf{m}=0$$

$$\mathbf{v}(\mathbf{r})=-{\int }_0^L\mathbf{G}(\mathbf{r},\mathbf{R}(s))\cdot {\mathbf{f}}^v(s)ds$$

$${\mathbf{f}}^v(s)=-\mathbf{\gamma }\cdot \dot{\mathbf{R}}(s),\mathbf{\gamma }={\gamma }_n\mathbf{n}\otimes \mathbf{n}+{\gamma }_t\mathbf{t}\otimes \mathbf{t}$$

$$EI{\partial }_s\kappa -aku-{\gamma }_n{g}_n^{\prime }+m^{\mathcal{A}}=0$$

$${\partial }_s^2\Delta \theta -\Delta \theta -g_n+m^{\mathcal{A}}=0$$

$${\partial }_s^2\Delta \theta +{\partial }_s^2{\partial }_t\Delta \theta -\Delta \theta -\frac{{\nu }_{\kappa }}{{\nu }_u}{\partial }_t\Delta \theta +m^{\mathcal{A}}=0$$

$${\partial }_t\left[\begin{array}{c}({\partial }_s^2-{\nu }_{\kappa }/{\nu }_u)\Delta \theta \\ m^{\mathcal{A}}\end{array}\right]=\left[\begin{array}{cc}1-{\partial }_s^2& -1\\ b_3& -b_1\end{array}\right]\left[\begin{array}{c}\Delta \theta \\ m^{\mathcal{A}}\end{array}\right]$$

$$\Delta \text{uptheta}(s,t)=\int \frac{dz}{2\pi }\sum _n\widetilde{\Delta \text{uptheta}}(q_n,z)e^{i(q_{n}s-zt)},\text{[similarly for mA(s,t)→mA(qn,z)]}$$

$$\left[\begin{array}{cc}(1+q_n^2)-iz(q_n^2+{\text{upnu}}_{\text{upkappa}}/{\text{upnu}}_u)& -1\\ b_3& iz-b_1\end{array}\right]\left[\begin{array}{c}\widetilde{\Delta \text{uptheta}}\\ \widetilde{m^{\mathcal{A}}}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

$$\big{[}(1+q_{n}^{2})-iz(q_{n}^{2}+\upnu_{\upkappa}/\upnu_{u})\big{]}(iz-b_{1})+b_{3}=0$$

$$(-q_n^2-i\text{upomega}{q}_n^2-1-i\text{upomega}{\text{upnu}}_{\text{upkappa}}/{\text{upnu}}_u)\widetilde{\Delta \text{uptheta}}+\widetilde{m^{\mathcal{A}}}=0$$

$$\underset{\text{passive elastic}}{\underbrace{-(1+q_n^2)\widetilde{\Delta \text{uptheta}}}}\underset{\text{passive viscous}}{\underbrace{-i\text{upomega}(q_n^2+{\text{upnu}}_{\text{upkappa}}/{\text{upnu}}_u)\widetilde{\Delta \text{uptheta}}}}+\underset{\text{active elastic + viscous}}{\underbrace{(G^{\prime }+i\text{upomega}{G}^{\prime \prime })\widetilde{\Delta \text{uptheta}}}}=0$$

$$k_{off}(f_{+})=k_{0}\exp\bigg{[}\frac{f_{+}}{f_{c}}\bigg{]}=k_{0}\exp\bigg{[}\frac{\bar{f}-f^{\prime}v_{d}}{f_{c}}\bigg{]}$$

$$k_{off}(f_{+})=\bar{k}_{off}\bigg{[}1-\frac{f^{\prime}v_{d}}{f_{c}}\bigg{]}$$

Full text

Introduction

Cilia and flagella are prototypical engines of micron-scale motility, used by the biological world in myriad contexts Goldstein (2015); Gaffney et al. (2011). They are classic examples of nonequilibrium active filaments that undergo spontaneous oscillations by converting stored or ambient energy into mechanical motion Marchetti et al. (2013). Their elemental structure, called axoneme, consists of cross-linked microtubule (MT) doublets and dynein motors which apply forces on MT, through adenosine triphosphate (ATP) hydrolysis, to cause periodic bending of the whole structure Nicastro et al. (2006). In addition to its importance in the cellular context, ciliary beating has been mimicked in synthetic filaments for applications in cargo transport, microfluidics and drug delivery Dreyfus et al. (2005); Paxton et al. (2004).

Naturally, there have been a number of studies devoted to understand the microscopic mechanisms of ciliary oscillations Brokaw (1971); Rikmenspoel (1978); Blum and Hines (1979); Riedel-Kruse et al. (2007); Gaffney et al. (2011); Bayly and Wilson (2015); Sartori et al. (2016). Most of them analyze beating based on passive elastic stresses calculated from measured filament waveform, models of active drive and dissipation stemming from external fluid. The dissipation is usually computed using slender body integral equations since filament waveform alone cannot be used to validate its form. Intuitively, hydrodynamics should play an essential role as cilia operate in the regime of low Reynolds number ($\sim 10^{-3}$). However, its contribution to synchronization and collective behavior of cilia is highly debated because most of these in vivo experiments are conducted on live cells where it is difficult to delineate hydrodynamics from other coupling mechanisms such as cell body rocking or intracellular means Goldstein et al. (2009); Geyer et al. (2013); Quaranta et al. (2015); Wan and Goldstein (2016); Klindt et al. (2016); Elgeti and Gompper (2013). In addition, biopolymers often have substantial internal friction Poirier and Marko (2002). Hence, to determine the dominance or lack of hydrodynamics in ciliary oscillations, an accurate measurement of the external viscous drag and its competition with the other elastic stresses for an in vitro system of cell-free isolated cilium is required, where cell body rocking and intracellular basal body coupling are eliminated. This can be accomplished by measuring its flow field along with the waveform.

Here, we present the first simultaneous measurement of the bending waveform and flow field of an isolated and active Chlamydomonas axoneme, beating near the critical ATP concentration at which oscillations set in, at high spatio-temporal resolution to address the role of external fluid in its beating. Our measurements demonstrate that hydrodynamic dissipation, accurately described by resistive force theory (RFT), is negligible compared to internal elastic stresses. Consequently, a dissipation mechanism internal to the filament is essential for stable driven oscillations, in contrast to the widely-held view that fluid friction is the only important source of dissipation Blum and Hines (1979); Sartori et al. (2016). We combine these insights with a theoretical model of filament motion that includes a minimal spring-dashpot form of active stress. We show there exist critical values of active stress beyond which the model exhibits oscillations, namely, active stresses should be strain softening and shear thinning.

Results

Experimental system

Isolated and reactivated axonemes of $\sim 11~{}\upmu$m length and $\sim 0.2~{}\upmu$m diameter Nicastro et al. (2006), purified from the unicellular algae Chlamydomonas reinhardtii, are clamped at one end on a glass cover slip. Their oscillations, with beat frequency $\nu_{b}\sim 16.22$ Hz, are approximately planar with an average centreline height $h\sim 0.9~{}\upmu$m from the surface. Passive microspheres are introduced into the suspension as tracers for measuring flow of the ambient fluid using particle tracking velocimetry (PTV) (see Materials and Methods). We capture motion of both the axoneme and tracers at high spatio-temporal resolution using a 60X phase objective coupled with a high speed camera (Fig. 1A and movie S1). The detailed experimental procedure is described in Material and Methods. The position of the axoneme is sampled at many points over its length to construct a global Chebyshev polynomial interpolant of the parametric form $\boldsymbol{R}(s)=\sum\boldsymbol{a}_{n}(t)T_{n}(s)$ where $0\leq s\leq L$ is the arc-length and $\boldsymbol{a}_{n}$ is a vector of Chebyshev coefficients for the $x$ and $y$ components of the position (Fig. 1B).

Filament geometry and mechanics

We perform experiments near the critical ATP concentration at which the axoneme transitions from quiescent to oscillatory state through an instability threshold (Fig. 1C) Camalet and Jülicher (2000). Near this threshold, details of the internal structure of the axoneme become irrelevant. Therefore, we model it as an inextensible, but shearable, active rod of uniform circular cross section of diameter $a$ and length $L$, with a centreline described by the curve $\boldsymbol{R}(s)$ to calculate stresses averaged over its cross-section. We attach an orthogonal Frenet-Serret frame to each point of the curve. Planar motion of the filament results in the curve tangent to be parametrised by the tangent angle $\theta$ as $\boldsymbol{t}(s)=[\cos\theta(s),\sin\theta(s)]$, which automatically satisfies the inextensibility constraint, $\boldsymbol{t}\cdot\boldsymbol{t}=1$. The position of the curve is obtained in terms of the tangent angle by integrating $\partial_{s}\boldsymbol{R}=\boldsymbol{t}$. The shear strain of the rod is $u(s)$ which is assumed to vanish at the clamped base, $u(0)=0$. Inextensibility then requires the shear strain to accommodate the filament bending as $u(s)=a[\theta(s)-\theta(0)]\equiv a\Delta\theta(s)$ and the kinematics of the filament is completely specified by $\theta$ or equivalently by shear angle, $\Delta\theta$ Brokaw (1971). Notably, the unobservable shearing of the filament can be inferred from its observable curvature, $\kappa=\partial_{s}\theta\equiv\partial_{s}\Delta\theta$. There is no assumption of small curvature in this kinematic description. Internal active stresses cause the filament to shear and, by the kinematic constraint, to curve.

We assume that the rod supports an internal stress whose stress and moment resultants on the cross section at $s$ are, respectively, $\boldsymbol{F}(s)$ and $\boldsymbol{M}(s)$ and that it is acted upon by forces and moments whose sum per unit length are, respectively, $\boldsymbol{f}$ and $\boldsymbol{m}$. Then, in the absence of inertia, the balance equations for force and torque are Ericksen and Truesdell (1957)

[TABLE]

Internal moments included in $\boldsymbol{m}$, which can exist only if the internal stress is antisymmetric, are generally omitted in the elasticity of rods, but are essential to our model Ericksen and Truesdell (1957). The above equations are closed by identifying the relevant forces, moments and their constitutive equations in terms of the kinematic variables. Integrating the force equation, the stress resultant can be expressed as $\boldsymbol{F}(s)=\int_{s}^{L}\boldsymbol{f}(s)\thinspace ds+\boldsymbol{F}(L)$, where $\boldsymbol{F}(L)=0$ at the free terminus of the rod. The only force per unit length relevant here is the external hydrodynamic drag, $\boldsymbol{f}^{v}$.

External viscous drag

A rod moving through a viscous fluid experiences a drag $\boldsymbol{f}^{v}(s)$ and creates a flow $\boldsymbol{v}(\boldsymbol{r})$. In the experimentally relevant limit of slow viscous flow and a slender rod, $L\gg a,$ the integral representation of Stokes equation gives

[TABLE]

This represents a distribution of Stokeslets of strength $\boldsymbol{f}^{v}(s)$ where $\boldsymbol{G}$ is a Green’s function of the Stokes equation. The matching of the fluid flow to the velocity of the rod at its surface yields the slender body integral equation whose formal solution is $\boldsymbol{f}^{v}(s)=-\int_{0}^{L}\boldsymbol{\gamma}(s,s^{\prime})\cdot\dot{\boldsymbol{R}}(s^{\prime})\thinspace ds^{\prime}$. The friction kernel is often approximated by a local form with a constant friction coefficient, $\boldsymbol{\gamma}(s,s^{\prime})=\boldsymbol{\gamma}\delta(s-s^{\prime})$. In this RFT limit, the drag is

[TABLE]

where $\gamma_{n}$ and $\gamma_{t}$ are the normal and tangential components of the friction coefficient and $\dot{\boldsymbol{R}}=\dot{R}_{n}\boldsymbol{n}+\dot{R}_{t}\boldsymbol{t}$ is the centreline velocity of the rod in terms of its normal and tangential components.

With the viscous drag thus determined in terms of the centreline velocity, we now validate RFT using experimentally measured instantaneous flow fields (Fig. 2, A and B). We compare these experimental flows with theoretically computed ones using Eqs. 2 and 3 (Fig. 2, C and D), wherein $\dot{\boldsymbol{R}}$ is determined from the measured waveform, $\gamma_{n}=4\pi\eta/ln(4h/a)=4.35$ mPa.s, $\gamma_{t}=\gamma_{n}/2$ Howard (2001) (fluid viscosity, $\eta=1$ mPa.s) and $\boldsymbol{G}$ is the Lorentz-Blake tensor for flow near a no-slip wall Blake (1971). Representative cuts along the experimental and theoretical flows show that there is a good agreement between the two (Fig. 2, E and F). A more comprehensive comparison is given by root mean square deviation of the flows, $RMSD=\sqrt{\sum_{i=1}^{NS}(v^{expt}_{i}-v^{th}_{i})^{2}/NS}$, where $NS$ is the number of grid points, $v^{expt}_{i}$ and $v^{th}_{i}$ are the experimental and theoretical flow magnitudes at the $i$-th grid point, respectively. RMSD of $v_{x}$, $v_{y}$ and $|\boldsymbol{v}|$ in Fig. 2 are (A and C) 4.32, 8.28 and 7.47; (B and D) 7.26, 7.89 and 9.81 $\upmu$m/s, respectively, all of which are within the Brownian noise regime implying a good match. Therefore, we have now verified by direct measurement that the hydrodynamic drag force is unambiguously given by RFT. This form of drag is used to evaluate the stress resultant, $\boldsymbol{F}$ and its normal component, $F_{n}$, is used in the torque balance equation below.

Scalar equation of motion of the filament

In slow viscous flow as is shown above, the drag acts in the plane of motion and the stress resultant remains in that plane. Therefore, the couple resultant $\boldsymbol{M}$ is normal to the plane of motion, along the frame binormal $\boldsymbol{b}$, and vanishes at the free terminus, giving $M(L)=0$. As the only non-zero components of $\boldsymbol{t}\times\boldsymbol{F}=F_{n}\boldsymbol{b}$ and $\boldsymbol{m}$ are normal to the plane, torque balance reduces to a scalar equation, $\partial_{s}M+F_{n}+m=0$. The dissipation in this equation is contained in the normal component of the stress resultant due to external viscous drag as $F_{n}(s)=\boldsymbol{n}(s)\cdot\int_{s}^{L}\boldsymbol{f}^{v}(s^{\prime})ds^{\prime}=-\gamma_{n}g^{\prime}_{n}(s)$ where g^{\prime}_{n}(s)=\boldsymbol{n}(s)\cdot\int_{s}^{L}\big{[}\dot{R}_{n}(s^{\prime})\boldsymbol{n}(s^{\prime})+\dot{R}_{t}(s^{\prime})\boldsymbol{t}(s^{\prime})/2\big{]}ds^{\prime} using Eq. 3. Herein, the filament velocity components can be expressed in terms of the tangent angle, $\theta$, using $\dot{\boldsymbol{R}}(s)=\int_{0}^{s}\dot{\boldsymbol{t}}(s^{\prime})ds^{\prime}$. The simplest elastic contributions for an inextensible and shearable active rod are that due to bending, $M=EI\kappa$, and shear, $m=-aku$ Brokaw (1971). The source of filament motion is included as the active moment per unit length, $m^{\mathcal{A}}$. Material parameters in these constitutive relations are bending rigidity, $EI=600$ pN$\upmu$m2 Xu et al. (2016); Sartori et al. (2016) and shear stiffness, $k=2000$ pN/$\upmu$m2 Minoura et al. (2001); Xu et al. (2016). Torque balance, closed by the constitutive relations and the expression for $F_{n}$, yields the following dynamical equation,

[TABLE]

All the passive terms in this equation can be expressed completely in terms of the tangent angle, $\theta$. Experimental data reveals that the tangent angle can be parametrised by a traveling waveform, $\theta(s,t)=\theta_{0}(s)\sin[\omega t-\phi(s)]+\bar{\theta}(s)$, where $\theta_{0}$, $\omega$, $\phi$ and $\bar{\theta}$ are the amplitude, angular frequency, phase and offset respectively (Fig. 3) Rikmenspoel (1978); Sartori et al. (2016). As we are interested in oscillations about the time averaged shape of the beat, $\bar{\theta}$, we define $\theta^{\prime}(s,t)=\theta(s,t)-\bar{\theta}(s)$. In the following, we focus on $\theta^{\prime}$ and drop the prime such that $\theta(s,t)\equiv\theta_{0}(s)\sin[\omega t-\phi(s)]$ which represents the dynamic oscillatory beat of the filament about the mean shape. Therefore, we use this parametric form, instead of the Chebyshev interpolant, to estimate all angle-dependent quantities in the dynamical equation.

The dimensionless equation of motion with $l_{\kappa}=\sqrt{EI/a^{2}k}=2.74~{}\upmu$m as the curvature penetration length scale and $\nu_{h}=EI/\gamma_{n}l_{\kappa}^{4}=2447$ Hz as the hydrodynamic relaxation frequency scale of the system is,

[TABLE]

Here, $g_{n}=g^{\prime}_{n}/l_{\kappa}^{2}\nu_{h}$ and $m^{\mathcal{A}}$ is rescaled by $l_{\kappa}^{2}/EI$. Using $\theta$ from the traveling wave parametrisation of experimental data and constitutive parameters from the literature, the space-time variation of elastic and viscous terms of Eq. 5 are plotted in Fig. 4 (A to C). It shows that the nonlinear viscous dissipation has a standing wave component in contrast to the linear elastic terms Rikmenspoel (1978). On comparing their colorbars, we conclude that the hydrodynamic dissipation, $g_{n}\sim\mathcal{O}$(0.01), is negligible compared to the elastic forces which are of $\mathcal{O}$(1).

The axoneme consumes energy in the form of ATP and exhibits stable oscillations. Such a continuously driven system can be in a dynamical steady state only when the elastic stresses due to the drive are balanced by some dissipation. Since external fluid friction cannot account for this dissipation, consistency demands that the internal stress has, in addition to an elastic component, a dissipative frictional component. Each kinematic degree of freedom, i.e. bending and sliding, can contribute to dissipation. Notably, bending friction in MTs has been experimentally demonstrated to be dominant over hydrodyamic friction for length scales smaller than $20~{}\upmu$m and attributed to slow structural rearrangements Poirier and Marko (2002); Brangwynne et al. (2007); Janson and Dogterom (2004). The bending friction coefficient of a 11 $\upmu$m long axoneme, $\Gamma_{\kappa}=1.6~{}\rm{pN}\upmu\rm{m^{2}s}$, is same as that of a single MT since it is an intensive quantity Brangwynne et al. (2007); Janson and Dogterom (2004). Though there is no experimental measure of the shear friction coefficient of an axoneme, $\Gamma_{u}$, several experimental studies suggest the presence of inter-MT sliding friction Minoura et al. (2001). We consider $\Gamma_{u}=10~{}\rm{pNs/}\upmu\rm{m}$ from estimates of nexin protein friction Riedel (2005) (see table S1). Earlier work introduced similar terms for internal viscous stresses for either stabilization of the numerical simulations of bend propagation in active filaments Brokaw (1972); Murase and Shimizu (1986); Hines and Blum (1978); Blum and Hines (1979) or on the basis of theory alone Bayly and Wilson (2015); Bayly and Dutcher (2016). However, here we have shown experimentally that such terms are necessary to completely account for dissipation in ciliary motion.

Dynamical equation without external friction

We, therefore, neglect the external viscous drag and include the internal viscous stresses to re-write the scalar torque balance equation, $\partial_{s}M+m=0$. Hence, the constitutive equations are modified as $M=EI\kappa+\Gamma_{\kappa}\dot{\kappa}$ and $m=-aku-\Gamma_{u}\dot{u}+m^{\mathcal{A}}$. The modified dynamical equation in terms of the shear angle, $\Delta\theta$, is $EI\partial_{s}^{2}\Delta\theta+\Gamma_{\kappa}\partial_{s}^{2}\partial_{t}\Delta\theta-a^{2}k\Delta\theta-a\Gamma_{u}\partial_{t}\Delta\theta+m^{\mathcal{A}}=0$. The negligible role of fluid friction leads to a second-order reaction-diffusion equation for the shear angle, rather than fourth-order partial differential equations that are commonly obtained when fluid friction is retained Bayly and Wilson (2015); Riedel-Kruse et al. (2007); Murase and Shimizu (1986). Identifying the two frequency scales of the system, a curvature relaxation frequency scale, $\nu_{\kappa}=EI/\Gamma_{\kappa}$ and a sliding relaxation frequency scale, $\nu_{u}=ak/\Gamma_{u}$, we note that $\nu_{\kappa}/\nu_{u}=9.38$ i.e both kinematical degrees of freedom contribute to dissipation. The dimensionless equation of motion with $\nu_{\kappa}$ as the frequency scale is,

[TABLE]

where $m^{\mathcal{A}}$ is rescaled by $l_{\kappa}^{2}/EI$. Figure 4 (D and E) show the variation of the internal viscous stresses over three beat cycles along the filament length. The colorbars in Fig. 4 (A to E) indicate that the internal viscous stresses due to bending friction of $\mathcal{O}(0.1)$ and shear friction of $\mathcal{O}(1)$ compete with the elastic stresses of $\mathcal{O}(1)$, unlike the negligibly small external drag of $\mathcal{O}(0.01)$.

The internal passive stresses being completely defined, a constitutive relation for the active moment in terms of kinematic variables is required. The axoneme having both elastic and viscous material parameters, the motor activity will induce a viscoelastic response. Hence, we assume a minimal spring-dashpot form for the dynamics of $m^{\mathcal{A}}$, $\partial_{t}m^{\mathcal{A}}+b_{1}m^{\mathcal{A}}=b_{3}\Delta\theta$ Camalet et al. (1999); Camalet and Jülicher (2000), whose Fourier representation describes the spring-like elastic and dashpot-like viscous dampening response of the system under oscillating shear as shown later. Here, the active stress is parametrised by two constants: $b_{1}$ controls the autonomous dynamics of $m^{\mathcal{A}}$ and $b_{3}$ controls the amount of feedback it receives from the sliding kinematics. We have chosen this simplest form of constitutive relation as it is first order in time, lowest order in wavenumber and linear in $\Delta\theta$. This minimal constitutive relation is most relevant for a coarse-grained description of axoneme as used here. This relation for $m^{\mathcal{A}}$ together with Eq. 6 form a pair of coupled equations,

[TABLE]

which emphasizes that the passive and active parts in this model are independent, but coupled, degrees of freedom.

These dissipative linear partial differential equations can sustain stable oscillations only in the presence of non-linearities whose choice will become crucial far from the instability threshold Camalet and Jülicher (2000); Oriola et al. (2017); Bayly and Wilson (2015). We now focus on the linear regime and seek to identify the threshold for the onset of oscillations and its frequency near threshold. This is relevant to our experimental results because the axoneme in this study beats at 60 $\upmu$M ATP, near the critical ATP concentration of 45 $\upmu$M at which oscillations set in (Fig. 1C). Hence, the axoneme is beating near the instability threshold, where the nonlinearity is weak and the oscillation frequency of the limit cycle is that of the linear analysis evaluated at the threshold Camalet and Jülicher (2000); Oriola et al. (2017).

We perform linear stability analysis on the coupled Eq. 7 in the Fourier domain (see section S2 for detailed analysis). The dispersion relation is quadratic in complex frequency $z$, \big{[}(1+q_{n}^{2})-iz(q_{n}^{2}+\nu_{\kappa}/\nu_{u})\big{]}(iz-b_{1})+b_{3}=0. The clamped head-free end boundary conditions of the filament discretize the wavenumbers to $q_{n}=\dfrac{n\pi}{2(L/l_{\kappa})}$ for odd $n$. Unstable oscillations are possible only when both $b_{1}<0$ and $b_{3}<0$ (Fig. 4F). These parameters can alternatively be related to the elastic, $G^{\prime}$, and viscous, $G^{\prime\prime}$, response coefficients of the active stress through its fundamental Fourier mode of oscillation frequency $\omega$ as $\widetilde{m^{\mathcal{A}}}=(G^{\prime}+i\omega G^{\prime\prime})\widetilde{\Delta\theta}$, where $G^{\prime}=b_{1}b_{3}/(b_{1}^{2}+\omega^{2})$ and $G^{\prime\prime}=-b_{3}/(b_{1}^{2}+\omega^{2})$ Riedel-Kruse et al. (2007); Camalet and Jülicher (2000). The sign of response coefficients determines the nature of active stress. If $G^{\prime},G^{\prime\prime}<0$, the system’s passive spring constant and friction coefficient get renormalized by the ATP dependent dynein activity. In our case, $b_{1},b_{3}<0\implies G^{\prime},G^{\prime\prime}>0$, i.e. activity works against the material response and leads to strain softening ($G^{\prime}>0$) and shear thinning ($G^{\prime\prime}>0$) in the system (see section S3 for detailed discussion). The experimental beat frequency of the axoneme ($\omega_{expt}=2\pi\nu_{b}/\nu_{\kappa}=0.272$), operating close to the threshold, constrains the parameters of the constitutive relation at hand, namely, $b_{1}=-0.121,~{}b_{3}=-0.85$ for the fundamental oscillatory mode (Fig. 4G). Elastoviscous response coefficients of the active stress computed with these values, $G^{\prime}=1.16,~{}G^{\prime\prime}=9.6$, imply that the viscous response of the active stress dominates the elastic one in our experiment. On comparing the nature of active stress so obtained with a microscopic model of load dependent detachment of motors in the experimentally relevant dominantly linear regime of the post-threshold dynamics Riedel-Kruse et al. (2007); Camalet and Jülicher (2000); Oriola et al. (2017), we infer that axonemal dyneins are low duty ratio motors (see section S4 for details), which agrees with previous studies Howard (2001).

We emphasize again that the existence and dominance of internal friction over hydrodynamic drag in isolated ciliary dynamics is borne out of experimental measurements alone and not through detailed modeling. Simultaneous measurement of flow field and waveform of an active cilium gives us the unambiguous measure of the hydrodynamic drag force i.e. the external fluid friction. The passive internal elastic stresses are calculated from the experimental waveform using minimal constitutive relations for bending and shear elasticity which are widely accepted in the literature for elasticity of rods. Comparing the experimentally measured fluid friction with that of the passive elasticity led to the conclusion that fluid friction is not enough to counteract the elastic stresses due to the active drive and consequently, stable ciliary oscillations need internal friction to reach their dynamical steady state. The above theoretical analysis of the filament equation of motion using a minimal constitutive relation for the active drive is simply to illustrate that oscillations exist in presence of internal friction too.

Discussion

In conclusion, to the best of our knowledge, this is the first direct experimental evidence of the negligible role of external fluid friction in ciliary oscillations near the instability threshold, thereby suggesting that internal dissipation mechanisms are essential to have a self-consistent understanding of ciliary beating. Our results are complementary to studies conducted on live cells/sperms Qin et al. (2015) and reactivated cilia at high ATP conditions Brokaw (1975) wherein changing the fluid viscosity alters both waveform and beat frequency of cilia by about 20%. At low ATP concentration investigated in this study, ciliary beat frequency and waveform is almost independent of the viscosity of ambient fluid Brokaw (1975). It is crucial to study this regime near the instability threshold because, here, the system is governed by linear terms that are generic and do not depend on the minor structural details of the system, leading to an understanding that is universal in nature Camalet and Jülicher (2000). As in this study, our results imply that one should include the contribution of internal friction with the obvious source of external fluid friction in constructing balance equations for active filaments in viscous fluids to correctly understand the dynamical regimes in which they operate.

This result will inspire further studies on the role of fluid friction in ciliary beating far from the threshold and in collective behavior of cilia. In addition, the measured flow field can distinguish internally active filaments from those driven by surface forces, for example, in phoretic chains (see fig. S1) Laskar et al. (2013). Therefore, a measurement of the flow field of an oscillating cilium provides vital information on its mechanisms of operation that cannot be obtained from the measurement of motion alone. Our forward approach of using the experimental insights from simultaneous waveform and flow measurements to build a theoretical model of active filament is in contrast to the existing reverse approaches of starting with microscopic models of active stress to match the observed macroscopic waveform of the filament Sartori et al. (2016); Bayly and Wilson (2015); Riedel-Kruse et al. (2007). Lastly, our approach can serve as a paradigm for analysis and regulation of any active slender body, both biological and synthetic one, in a viscous fluid.

Materials and Methods

Axoneme purification and reactivation: Wild-type Chlamydomonas reinhardtii cells (strain: CC1690) are synchronously grown in 12:12 hours light:dark cycle in TAP+P medium (tris acetate phosphate) Alper et al. (2013). They are collected 2-3 hrs after the beginning of light cycle, at OD750 (optical density at 750 nm) $\approx 0.15-0.25$, followed by washing in HES buffer [10 mM HEPES (pH 7.4), 1 mM EGTA, 4% Sucrose] thrice. The plasma membrane covering the cells is disrupted by adding 0.15% of a non-ionic detergent IGEPAL CA-630 (I8896, Sigma-Aldrich) in HMDEK buffer [30 mM HEPES (pH 7.4), 5 mM MgSO4, 1 mM dithiothreitol (DTT), 1 mM EGTA, 50 mM K-acetate] to the cell pellet. The cells in IGEPAL are then kept in ice for $\sim 5-7$ hours. These cells are devoid of cell membranes and are called cell models. Some of them shed their demembranated flagella, called axonemes, due to weakening of the attachment to the cell body. Isolated axonemes are then separated from the cell models by centrifugation. Axonemes are then mixed with 30% saturated sucrose, flash freezed and stored at -80*∘*C for long term usage.

This method of purification, by long-term exposure to detergent, yields nonsticky axonemes, in contrast to the commonly used dibucaine procedure Alper et al. (2013); Wakabayashi and Kamiya (2015), hence essential for flow field measurements. These non-motile axonemes regain their motility in presence of 45 $\upmu$M to 1 mM ATP in HMDEKP buffer (HMDEK + 1% PEG-20k) and are said to be reactivated. There was no distinction in the nature of beating of the frozen axonemes, when thawed, from those that were not frozen. We use an ATP regeneration system, to hold the ATP concentration constant within the sample for approximately 40 mins, enough for imaging approximately 6-7 isolated axonemes in one sample. The ATP regeneration system composes of 6 mM sodium creatine phosphate (27920, Sigma-Aldrich) and 40 U/ml creatine phosphokinase (CPK) from bovine heart (C7886, Sigma-Aldrich). The commercially available CPK is typically oxidized, hence we reduce it by DTT at room temperature for efficient reactivation Favre and Muellhaupt (2005).

Tracer particles: The tracer particles chosen, for flow field measurements, are neutrally buoyant (polystyrene microspheres with a density of $\rm 1.055~{}g/cc$) and small (diameter of 200 nm, the lowest size that can be used with diffraction-limited optical imaging) so that their motion is nearly identical to the fluid in which they reside. Commercially available charge-stabilized microspheres stick to each other and to the axoneme due to screening of electrostatic interactions by the divalent ions and salts present in the reactivation buffer. We graft long chains of block copolymer PLL-PEG20k [poly-L-lysine (PLL), P7890, 15 to 30 kDa, Sigma-Aldrich; mPEG-SVA-20k, NANOCS], onto 200-nm negatively charged sulfate latex beads (S37491,Thermo Fisher Scientific) to impart additional steric stabilization.

Surface modification of glass surfaces and sample preparation: The coverslips and slides are cleaned with a hot soap solution (1% Hellmanex III), followed by rinsing with ethanol and 100 mM potassium hydroxide. We graft polyacrylamide brush on these clean glass surfaces to suppress depletion interaction of beads and filament with the surfaces. The coverslips are further modified with randomly located sticky patches of 0.05% PLL over the polyacrylamide brush, to clamp the axonemes at one end. We introduce the reactivation buffer containing axonemes and tracer particles inside a sample chamber made up of a glass slide and coverslip sandwich with double sided tape of thickness 65 $\upmu$m as the spacer.

Imaging and depth of focus: The sample is mounted on an inverted microscope (Olympus IX83), equipped with 60x oil immersion phase objective [0.65 to 1.25 numerical aperture (NA), UPlanFL N, PH3] and connected to a high speed complementary metal oxide semiconductor (CMOS) camera (frame rate, 1200 frames per second; Phantom Miro C110, Vision Research) for imaging. We use an intermediate NA between 0.65 and 1.25 for the 60x variable NA phase objective, to capture most of the filament beat in focus. We have measured the depth of focus ($\Delta Z$) of the objective at this intermediate NA to be $\rm 1.4~{}\upmu m$. We only image time lapse sequences of those axonemes which are clean, clamped at one end with proper beating, having $80$ to $90\%$ of the filament in focus at a $\Delta Z$ of $\rm 1.4~{}\upmu m$ and do not have another reactivated filament in the surrounding area of $\rm 30~{}\upmu m~{}\times~{}30~{}\upmu m$. The frame rate in the high speed imaging is suitably chosen to capture approximately 1000 beat cycles per axoneme, with $\rm\sim 55$ to $75$ conformations per beat cycle, for example 1200 frames per second for the axoneme in movie S1, whose beat frequency is approximately 16.63 Hz. Furthermore, the theoretically computed flow field is also depth averaged over this $\Delta Z$ for appropriate comparison with the experimental flow field.

Extracting filament conformation from images: We manually identify position coordinates of the axoneme for three beat cycles from the recorded time lapse sequences, which are then smoothened along $s$ and $t$ by using Savitzky-Golay filter of orders 3 and 5, respectively. The filament positions $\boldsymbol{R}(s)$, arc length $s$ are nondimensionalized with the filament length $L$ and time $t$ with beat period $T$, followed by setting the clamped end of the filament to (0,0). These nondimensionalized experimental positions are used to construct their Chebyshev interpolant as mentioned in the Results section. Higher order derivatives of the interpolation lose accuracy at the ends of the axoneme and data from those parts are, accordingly, discarded (fig. S2). Hence, we consider $s/L\in[0.1,0.9]$ when computing derivatives higher than first order and neglect the end values (as shown in Fig. 4, A to E).

Particle Tracking Velocimetry (PTV): The recorded movies are background subtracted, before tracking the tracer displacements, to consider only those tracers which are not stuck to the coverslip and are following the flow. We track the tracer displacements in between two filament conformations with an infinitesimal time gap, $\Delta t$ (for example, $\Delta t\approx 3.33$ ms for power stroke waveforms and $\Delta t\approx 4.98$ ms in recovery stroke waveforms). As the filament has a nonuniform beat frequency $\nu_{b}=16.63\pm 0.62$ Hz, over $\sim 940$ beat cycles, each conformation will not exactly repeat itself after one beat period. Hence, we cross-correlate each filament conformation for a given beat period with the whole recorded sequence of $\sim$ 940 beat cycles. The correlation peaks indicate the frames that have similar conformation. The stack of matched conformations is checked manually to delete any conformation with more than $10\%$ dissimilarity. The displacement of tracers in between these two conformations is calculated using standard MATLAB tracking routines and velocity vectors are obtained from $\sim 500$ to $700$ beat cycles. The resulting velocity vectors are placed on $29\times 29$ mesh grid, each of size $\rm 0.74~{}\upmu m~{}\times~{}0.74~{}\upmu m$ over the image, and the mean at each grid point is computed. The gridded velocity vectors are further smoothened using a $2\times 2$ point averaging filter.

Supplementary Materials

It includes the following:

Section S1. Tension forces in the filament

Section S2. Linear stability analysis

Section S3. Why strain softening and shear thinning facilitates the instability to oscillations?

Section S4. Comparison with microscopic load dependent detachment model of dynein motors

Section S5. Estimates of $\Gamma_{u}$ from literature

Movie S1. Movie of an isolated and reactivated axoneme in presence of tracers.

Table S1. Possible values of shear friction coefficient from the literature.

Fig. S1. Active filaments driven by slip at boundary vs driven internally by motors.

Fig. S2. Chebyshev differentiation of traveling wave parameters.

Fig. S3. Tension forces in the filament.

Fig. S4. Existence of oscillations with varying sliding friction coefficient.

References Prasath et al. (2016); Chakrabarti and Saintillan (2019); Yokota and Mabuchi (1994); Bayly and Wilson (2014).

Acknowledgements: D.M. and P.S. thank J. Sutradhar, K. Ghosh, M. Jain, G. Prasath, R. Govindarajan, T. G. Murthy, and D. Ghosh for useful discussions; B. J. Rao and S. Upadhyaya for help in growing Chlamydomonas; K.-i. Wakabayashi for help in purifying cilia; and S. Ramaswamy for important discussions and critical reading of our manuscript. D.M. and P.S. thank A. Baskaran for crucial inputs during the late stages of this work and continuous support during peer review process. R.A. thanks R. E. Goldstein for helpful discussions and suggesting appropriate references. We thank A. Laskar for assistance in the initial stages of this work and for sharing codes and data on phoretic particle chains.

Funding: This work was supported by the DBT/Wellcome Trust India Alliance Fellowship (grant number IA/I/16/1/502356) awarded to P.S. R.A. acknowledges support from the Isaac Newton Trust.

Author Contributions: R.A. and P.S. designed the research, D.M. performed the research, and all authors wrote the paper.

Competing Interests: The authors declare that they have no competing interests.

Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Supplementary Materials

Section S1. Tension forces in the filament

The tangential component of the stress resultant is the tension within the filament Prasath et al. (2016), which we can compute from filament velocity as $F_{t}(s)=\mathbf{t}(s)\cdot\int_{s}^{L}\mathbf{f}^{v}(s^{\prime})ds^{\prime}=-\mathbf{t}(s)\cdot\int_{s}^{L}\mathbf{\gamma}\cdot\dot{\mathbf{R}}(s^{\prime})ds^{\prime}=-\gamma_{n}g_{t}(s)$, where, $g_{t}(s)=\mathbf{t}(s)\cdot\int_{s}^{L}[\dot{R}_{n}(s^{\prime})\mathbf{n}(s^{\prime})+(\dot{R}_{t}(s^{\prime})/2)\mathbf{t}(s^{\prime})]ds^{\prime}$. The tension force in the filament is also too small (fig. S3) compared to the internal elastic forces.

Section S2. Linear stability analysis

We Fourier transform the coupled equations of motion (Eq. 7 in main text) in space and time with the following convention

[TABLE]

where $q_{n}$ is the discretized wavenumber of the $n$th mode due to finite length of the filament and $z$ is the complex frequency. In matrix form,

[TABLE]

For non-trivial solution, the determinant of the above matrix must be zero. Hence the dispersion relation is

[TABLE]

This is a quadratic equation in the complex frequency $z$ of the form $A(q_{n})z^{2}+B(q_{n})z+C(q_{n})=0$ whose coefficients are $A(q_{n})=q_{n}^{2}+\upnu_{\upkappa}/\upnu_{u}$, B(q_{n})=i\big{[}(1+q_{n}^{2})+b_{1}(q_{n}^{2}+\upnu_{\upkappa}/\upnu_{u})\big{]} and $C(q_{n})=b_{3}-b_{1}(1+q_{n}^{2})$. Roots are z_{1,2}=\big{[}-B(q_{n})\pm\sqrt{B^{2}(q_{n})-4A(q_{n})C(q_{n})}\big{]}/2A(q_{n}). As the time dependent component in the solution is $e^{-izt}$, existence of the real part of the root will imply that the solution is oscillatory and sign of the imaginary part of $z$ will decide if the solution is growing or decaying. Hence, the conditions for unstable oscillations are $Im[z]>0$ and $Re[z]\neq 0$. The frequency of the oscillations is therefore given by $\upomega=-Re[z]$.

Section S3. Why strain softening and shear thinning facilitates the instability to oscillations?

We Fourier transform the dynamical equation of motion (Eq. 6 in the main text) in space and time under oscillating shear of fundamental frequency $\upomega$ i.e. $\Delta\uptheta\sim\widetilde{\Delta\uptheta}e^{i(q_{n}s+\upomega t)}$ where $\upomega=-Re(z)$ is real and $q_{n}$ is the discretized wavenumber of the $n$th mode.

[TABLE]

Replacing the fundamental Fourier mode of active stress as $\widetilde{m^{\mathcal{A}}}=(G^{\prime}+i\upomega G^{\prime\prime})\widetilde{\Delta\uptheta}$, where where $G^{\prime}$ and $G^{\prime\prime}$ corresponds to elastic and viscous response of the system, respectively and are related to $b_{1},b_{3}$ such that $b_{1},b_{3}<0\implies G^{\prime},G^{\prime\prime}>0$ (see main text), we obtain

[TABLE]

Hence, the passive elastic and viscous terms in the equation of motion are negative i.e. they resist the sliding caused by the active dynein motors. Now, if $G^{\prime},G^{\prime\prime}<0$, the system’s passive spring constant and friction coefficient get renormalized by the ATP dependent dynein activity and thus, both the active and passive components of the system resist the sliding ultimately leading to a quiescent stable state of the filament. On the other hand, if $G^{\prime},G^{\prime\prime}>0$, dynein motors work against the material response so that the system becomes unstable and undergoes oscillations. Here, $G^{\prime}>0$ indicates that elastic stresses reduce within the axoneme as motor activity increases which is called ‘strain softening’ and $G^{\prime\prime}>0$ indicates that the axoneme becomes less viscous with increasing motor activity which is called ‘shear thinning’.

We note that nonlinear viscoelastic effects are not needed for an active material, such as the axoneme, to shear thin/strain soft, because the nonequilibrium active stresses generated by the dynein motors can produce structural rearrangements within the axoneme by binding and unbinding the dynein cross-bridges across microtubule doublets. This is in contrast to passive equilibrium systems which have to be intrinsically nonlinear to strain soft/shear thin under large external shear as thermal energy alone is insufficient to drive structural rearrangements in such systems.

Section S4. Comparison with microscopic load dependent detachment model of dynein motors

In the microscopic load dependent detachment model, the motor detachment rate ($k_{off}$) is assumed to increase exponentially with increasing load (i.e. single motor force, $f_{+}$), which in turn decreases linearly with the dynein sliding speed ($v_{d}$) by the force-velocity relationship of the motors Riedel-Kruse et al. (2007); Oriola et al. (2017); Chakrabarti and Saintillan (2019).

[TABLE]

where $\bar{f}$ is the dynein stall force, $f^{\prime}$ is the slope of the dynein force-velocity curve, $f_{c}$ is the characteristic unbinding force generally given by $f_{c}\approx\bar{f}/2$ Riedel-Kruse et al. (2007); Oriola et al. (2017). The force-velocity slope is given by $f^{\prime}=\bar{f}/v_{0}$, where $v_{0}$ is the dynein velocity at zero load. In the limit of low sliding speed i.e. $\dfrac{f^{\prime}v_{d}}{f_{c}}<<1$, the above exponential relation linearizes to

[TABLE]

where $\bar{k}_{off}=k_{off}(\bar{f})=k_{0}\exp(\bar{f}/f_{c})$ is the motor detachment rate at stall (refer to Eq. B4 in Appendix B of Riedel-Kruse et al. (2007)). Now let us see if our experiments near the critical ATP concentration satisfy the linearizing condition, $\dfrac{f^{\prime}v_{d}}{f_{c}}<<1\implies\dfrac{v_{d}}{v_{0}/2}<<1$.

Axoneme being a cross-linked filament, the angular speed of the axoneme ($\partial_{t}\Delta\uptheta$) is related to the sliding speed per dynein motor as $v_{d}=a_{MT}\partial_{t}\Delta\uptheta/L\uprho\bar{p}$ Riedel-Kruse et al. (2007)). Here $a_{MT}=24$ nm is the MT interdoublet spacing in which the dyneins work Nicastro et al. (2006), $L$ is the filament length at which angular speed is calculated from experiments, $\uprho=198~{}\upmu$m*-1* is motor density Nicastro et al. (2006) and $\bar{p}=0.02$ is the fraction of motor domains that are attached to MT, equivalent to the duty ratio Howard (2001); Oriola et al. (2017). Therefore, $L\uprho\bar{p}$ is the total number of motors bound to a single MT along the length of the filament. At 60 $\upmu$M ATP, $\partial_{t}\Delta\uptheta\approx 80$ rad/s (from Fig. 4E of main text) at $L\approx 9~{}\upmu$m. Using these values $v_{d}\approx 54$ nm/s. Earlier experiments have measured the zero load dynein velocity at 60 $\upmu$M ATP to be $v_{0}\approx 2~{}\upmu$m/s Howard (2001); Yokota and Mabuchi (1994). Hence the ratio of dynein sliding speed in our experiments to the half of its zero load velocity $\dfrac{v_{d}}{v_{0}/2}\approx 0.05<<1$. To summarize, the axoneme beating near the instability threshold at 60 $\upmu$M ATP is in the dominantly linear regime (equivalently weakly nonlinear regime) of the post-bifurcation dynamics. Therefore, our choice of linear constitutive relationship for the active moment and the associated linear stability analysis is valid near the instability threshold (also refer to Fig. 3a and associated text of Oriola et al. (2017)).

Now, that we have shown our experiments are consistent with the condition for linearizing the exponential dependence of motor detachment rate on sliding speed, we connect our constitutive equation for active moment to this linearized version of microscopic motor dynamics model (sliding control motor coordination model) proposed by Riedel-Kruse et al. Riedel-Kruse et al. (2007). In Riedel-Kruse et al. (2007), the active shear force per unit length is related to the shear displacement by a response function, $\upchi=K+i\upomega\uplambda$, as per their notation. The sign convention of $G^{\prime}$ and $G^{\prime\prime}$ is opposite to the elasto-viscous response coefficients ($K$ and $\uplambda$) of Riedel-Kruse et al. (2007) (also of Camalet and Jülicher (2000); Camalet et al. (1999)), as active drive in these references has opposite sign convention. The equivalence of our response coefficients to their microscopic model of load dependent detachment of motors (refer to Eq. B9 in Appendix B of Riedel-Kruse et al. (2007)) are as follows:

(a) $G^{\prime}\equiv-a^{2}K=2a^{2}\uprho\bar{f}\frac{f^{\prime}}{f_{c}}\bar{p}(1-\bar{p})\frac{\upomega^{2}\bar{\uptau}}{1+(\upomega\bar{\uptau})^{2}}$. This quantity is always positive.

(b) G^{\prime\prime}\equiv-a^{2}\uplambda=-2a^{2}\uprho f^{\prime}\bar{p}\big{[}1-\frac{\bar{f}}{f_{c}}\frac{(1-\bar{p})}{\{1+(\upomega\bar{\uptau})^{2}\}}\big{]}. This quantity can be positive or negative.

In the above expressions, $\bar{\uptau}$ is the relaxation time of motor attachment/detachment and all other variables are already defined in the preceding paragraph. The sign in (b) depends on $\bar{\uptau}$ and $\bar{p}$ as: (i) $G^{\prime\prime}>0$ for $\upomega\bar{\uptau}<<1$ and $\bar{p}\sim 0$ and (ii) $G^{\prime\prime}<0$ for $\upomega\bar{\uptau}>>1$ and $\bar{p}\sim 1$. For our active filament $G^{\prime\prime}>0$ i.e. $b_{3}<0$ asserts that $\upomega\bar{\uptau}<<1$ and $\bar{p}\sim 0$. This means that the axonemal dynein motors are short lived with low duty ratio, which agrees with Howard (2001).

Section S5. Estimates of $\Gamma_{u}$ from literature

All parameters, except the shear friction coefficient, have been experimentally measured for an axoneme or at least for microtubules in the existing literature as mentioned in the main text. There is uncertainty in the value of $\Gamma_{u}$. Following (table S1) are the values used for this coefficient for constructing active filament models in the literature, except for the last entry which is an experimental study on shear elasticity.

We note that the exact value of this coefficient does not affect the existence of oscillation in our theoretical model (fig. S4) rather it modulates the magnitudes of the viscoelastic response coefficients of the active stress. The sliding friction coefficient $\Gamma_{u}$ appears in the dimensionless dynamical equation as the ratio $\upnu_{\upkappa}/\upnu_{u}$. Variation of $\upnu_{\upkappa}/\upnu_{u}$ implies variation in $\Gamma_{u}$ because all parameters except $\Gamma_{u}$ are experimentally known for an axoneme/MT. Figure S4 shows that unstable oscillations exists for $\upnu_{\upkappa}/\upnu_{u}\in[0.1,~{}50]$. Correspondingly, the varation of $\Gamma_{u}$ in this range is $[10^{-7},~{}0.5\times 10^{-4}]$ Ns/m.

Movie caption

Movie S1. Movie of an isolated and reactivated axoneme in presence of tracers. High speed phase contrast movie of a reactivated and clamped Chlamydomonas axoneme at 60 $\upmu$M ATP in presence of 200 nm tracer particles. Scale bar, 5 $\upmu$m.