## Abstract

The physical basis of flagellar and ciliary beating is a major problem in biology which is still far from completely understood. The fundamental cytoskeleton structure of cilia and flagella is the axoneme, a cylindrical array of microtubule doublets connected by passive cross-linkers and dynein motor proteins. The complex interplay of these elements leads to the generation of self-organized bending waves. Although many mathematical models have been proposed to understand this process, few attempts have been made to assess the role of dyneins on the nonlinear nature of the axoneme. Here, we investigate the nonlinear dynamics of flagella by considering an axonemal sliding control mechanism for dynein activity. This approach unveils the nonlinear selection of the oscillation amplitudes, which are typically either missed or prescribed in mathematical models. The explicit set of nonlinear equations are derived and solved numerically. Our analysis reveals the spatio-temporal dynamics of dynein populations and flagellum shape for different regimes of motor activity, medium viscosity and flagellum elasticity. Unstable modes saturate via the coupling of dynein kinetics and flagellum shape without the need of invoking a nonlinear axonemal response. Hence, our work reveals a novel mechanism for the saturation of unstable modes in axonemal beating.

## 1. Introduction

Cilia and flagella play a crucial role in the survival, development, cell feeding and reproduction of micro-organisms [1]. These lash-like appendages follow regular beating patterns which enable cell swimming in inertialess fluids [2]. Bending deformations of the flagellum are driven by the collective action of ATP-powered dynein motor proteins, which generate sliding forces within the flagellar cytoskeleton, named axoneme [3]. This structure has a characteristic ‘9+2’ composition across several eukaryotic organisms, corresponding to nine peripheral microtubule doublets in a cylindrical arrangement surrounding a central pair of microtubules [1]. Additional proteins, such as the radial spokes and nexin cross-linkers, connect the central to the peripheral microtubules and resist free sliding between the microtubule doublets, respectively. Each doublet consists of an A-microtubule in which dyneins are anchored at regular intervals along the length of the doublets, and a B-microtubule, where dynein heads bind in the neighbouring doublet [1,4]. In the presence of ATP, dyneins drive the sliding of neighbouring microtubule doublets, generating forces that can slide doublets apart if cross-linkers are removed [5]. In the presence of cross-linkers, sliding is transformed into bending. Remarkably, this process seems to be carried out in a highly coordinated manner, in such a way that when one team of dyneins in the axoneme is active, the other team remains inactive [6]. This mechanism leads to the propagation of bending undulations along the flagellum, as commonly observed during the movement of spermatozoa [4].

Many questions still remain unanswered on how dynein-driven sliding causes the oscillatory bending of cilia and flagella [7]. Over the last half a century, intensive experimental and theoretical work has been done to understand the underlying mechanisms of dynein coordination in axonemal beating. Different mathematical models have been proposed to explain how sliding forces shape the flagellar beat [8–16]. Coordinated beating has been hypothesized considering different mechanisms such as dynein’s activity regulation through local axonemal curvature [9–11,16], owing to the presence of a transverse force (t-force) acting on the axoneme [12] and by shear displacements [13–15]. Other studies also examined the dynamics of flagellar beating by prescribing its internal activity [17,18] or by considering a self-organized mechanism independent of the specific molecular details underlying the collective action of dyneins [13,14,19]. In particular, the latter approach, although general from a physics perspective, does not explicitly incorporate dynein kinetics along the flagellum, which has been shown to be crucial in order to understand experimental observations on sperm flagella [20–22]. Load-accelerated dissociation of dynein motors was proposed as a mechanism for axonemal sliding control, and was successfully used to infer the mechanical properties of motors from bull sperm flagella [15]. In contrast, dynamic curvature regulation has been recently proposed to account for *Chlamydomonas* flagellar beating [16]. In the previous studies, linearized solutions of the models were fit to experimental data; however, it is not clear that such results still hold at the nonlinear level. Recent studies also investigated the emergence and saturation of unstable modes for different dynein control models [23,24]; however, saturation of such unstable modes was not self-regulated, but achieved via the addition of a nonlinear elastic contribution in the flagellum constitutive relation. Nevertheless, predictions on how dynein activity influences the selection of the beating frequency, amplitude and shape of the flagellum remain elusive. Here, we provide a microscopic bottom-up approach and consider the intrinsic nonlinearities arising from the coupling between dynein activity and flagellar shape, regarding the eukaryotic flagellum as a generalized Euler-elastica filament bundle [25]. This allows a close inspection on the onset of the flagellar bending wave instability, its transient dynamics and later saturation of unstable modes, which is solely driven by the nonlinear interplay between the flagellar shape and dynein kinetics.

We first derive the governing nonlinear equations using a load-accelerated feedback mechanism for dynein along the flagellum. The linear stability analysis is presented, and eigenmode solutions are obtained similarly to references [15,24], to allow analytical progress and pedagogical understanding. The nonlinear dynamics far from the Hopf bifurcation is studied numerically and the resulting flagellar shapes are further analysed using principal component analysis [26,27]. Finally, bending initiation and transient dynamics are studied subject to different initial conditions.

## 2. Continuum flagella equations

We consider a filament bundle composed of two polar filaments subjected to planar deformations. Each filament is modelled as an inextensible, unshearable, homogeneous elastic rod, for which the bending moment is proportional to the curvature and the Young modulus is *E*. The filaments are of length *L* and separated by a constant gap of size *b*, where *b*≪*L* (figure 1*c*). We define a material curve describing the shape of the filament bundle centreline as ** r**(

*s*,

*t*). The positions of each polar filament forming the bundle read

*s*along its length defined by the normal vector to the centreline

*ϕ*≡

*ϕ*(

*s*,

*t*) the angle between the tangent vector

*x*-axis). The subscripts (+) and (−) refer to the upper and lower filaments, respectively (figure 1

*c*). The shape of the bundle is given at any time by the expression:

*s*,

*t*), denoted as sliding displacement:

*ϕ*

_{0}≡

*ϕ*(0,

*t*). For simplicity, we have set any arc length incongruity between the two filaments at the base to zero, and we consider the filaments clamped at the base. A similar approach can be used to include basal compliance and other types of boundary conditions at the base (e.g. pivoting or free swimming head). Here we centre our study on the nonlinear action of motors along the flagellum. We aim to study the active and passive forces generated at each point along the arc length of the filament bundle. We define

*s*at time

*t*on the plus-filament owing to the action of active and passive forces (figure 1). By virtue of the action–reaction law, the minus-filament will experience a force density −

**at the same point. Next, consider that**

*f**N*dyneins are anchored at each polar filament in a region

*l*

_{c}around

*s*, where

*l*

_{c}is much smaller than the length of the flagellum

*L*and much larger than the length of the regular intervals at which dyneins are attached along the microtubule doublets. We shall call

*l*

_{c}the ‘tug-of-war’ length. We define

*n*

_{±}(

*s*,

*t*) as the number of bound dyneins in a region of size

*l*

_{c}around

*s*at time

*t*which are anchored in the plus- or minus-filament, respectively. We consider a tug-of-war at each point

*s*along the flagellum with two antagonistic groups of

*N*dyneins. The elastic sliding resistance between the two polar filaments exerted by nexin cross-linkers is assumed to be Hookean with an elastic modulus

*K*. Thus, the internal force density

*f*(

*s*,

*t*) reads

*ρ*≡

*l*

^{−1}

_{c}is the density of tug-of-war units along the flagellum and

*F*

_{±}(

*s*,

*t*) is the load per motor each group of dyneins experiences owing to the action of the antagonistic group. The stresses on the filament bundle are given by a resultant contact force

**(**

*N**s*,

*t*) and resultant contact moment

**(**

*M**s*,

*t*) acting at the point

**(**

*r**s*,

*t*). The internal force density,

*f*(

*s*,

*t*), contributes only to the internal moment of the bundle

*M*(

*s*,

*t*) reads

_{s}

*ϕ*

_{±}≈∂

_{s}

*ϕ*for bundles characterized by

*b*≪

*L*. The combined bending stiffness of the filament bundle is given by

*E*

_{b}=2

*EI*, where

*I*is the second moment of the area of the external rods.

Dynein kinetics is modelled by using a minimal two-state mechanochemical model with states *k*=1,2, corresponding to microtubule bound or unbound dyneins, respectively. Because the sum of bound and unbound motors at *s* remains constant at all times, we study only the plus and minus bound motor distributions *n*_{±}(*s*,*t*). Dyneins bind with rates *π*_{±} and unbind with rates *ε*_{±} (figure 1*b*). The corresponding bound motor population dynamics reads
*π*_{±}=*π*_{0}(*N*−*n*_{±}), *ε*_{0} and *π*_{0} are constant rates and *f*_{c} is the characteristic unbinding force. Motivated by previous studies on the collective action of molecular motors [15,28], we assume an exponential dependence of the unbinding force on the resulting load. By considering that dyneins fulfil a linear velocity–force relationship with stall force *f*_{0} and velocity at zero load *v*_{0}, the loads are given by *F*_{±}(*s*,*t*)=±*f*_{0}(1∓Δ_{t}/*v*_{0}). The stall force *f*_{0} is defined as the absolute value of the load a motor experiences at stall (Δ_{t}=0). Substituting the different definitions, the internal force density *f*(*s*,*t*) reads
*n*′≡*n*_{+}+*n*_{−}. For simplicity, we derive the equations governing the tangent angle *ϕ* in the limit of small curvature (but possibly large amplitudes) such that tangential forces can be neglected. The derivation for arbitrary large curvature is also presented in the electronic supplementary material. Using resistive force theory in the limit of small curvature, we consider only normal forces along the flagellum obtaining *ζ*_{⊥}*ϕ*_{t}=−*M*_{sss}, where *ζ*_{⊥} is the normal drag coefficient [8,29]. Combining the last expression with equation (2.4), we have
*L*, time with respect to the correlation time of the system *τ*_{0}=1/(*ε*_{0}+*π*_{0}), motor number with respect to *N*, internal force density with respect to *f*_{0}*ρN* and sliding displacement with respect to *b*. The correlation time defines how fast the motors will respond to a change in load. We also define Sp=*L*(*ζ*_{⊥}/*E*_{b}*τ*_{0})^{1/4}, *μ*≡*Kb*^{2}*L*^{2}/*E*_{b}, *μ*_{a}=*bf*_{0}*ρNL*^{2}/*E*_{b} and *ζ*≡*b*/*v*_{0}*τ*_{0}. The sperm number, Sp, characterizes the relative importance of bending forces to viscous drag. The parameter *μ* measures the relative importance of the sliding resistance compared with the bending stiffness [25]. On the other hand, the parameter *μ*_{a} denotes the activity of dyneins, measuring the relative importance of motor force generation compared with the bending stiffness of the bundle. Finally, *ζ* denotes the ratio of the bundle diameter and the characteristic shear induced by the motors. The dimensionless sperm equation in the limit of small curvature reads [2,13,18]
*f*(*s*,*t*) takes the form:
*ϕ*−*ϕ*_{0}. Because the flagellar base is clamped, without loss of generality, we set *ϕ*_{0}=0. Combining equations (2.8) and (2.9), we obtain the nonlinear dynamics for the tangent angle:
*n*_{±} reads
*η*≡*π*_{0}/(*π*_{0}+*ε*_{0}) is the duty ratio of the motors and

## 3. Parameter values

We present the choice of parameters based on experimental studies on sperm flagella and on the green algae *Chlamydomonas*. We first discuss the passive properties of a flagellum. The typical length of a human flagellum is *L*≃50 μm, and the axonemal diameter is found to be *b*≃200 nm [4]. The bending stiffness of the filament bundle has been reported to be *E*_{b}≃0.9×10^{−21} N m^{2} for sea-urchin sperm [4,25] and *E*_{b}≃1.7×10^{−21} N m^{2} for bull sperm [15]. On the other hand, the interdoublet elastic resistance from demembranated flagellar axonemes of *Chlamydomonas* yields an estimated spring constant 2×10^{−3} N m^{−1} for 1 μm of axoneme [30], thus *K*≃2×10^{3} N m^{−2}. Finally, typical medium viscosities for sperm flagella range from *ζ*_{⊥}≃10^{−3} Pa s in low viscous media to *ζ*_{⊥}≃1 Pa s in high viscous media [4,18].

Next, we discuss the mechanochemical parameters associated with axonemal dynein. Axonemal dynein are subdivided into inner and outer arms depending on its position in the axoneme, and can be found in heterodimeric and monomeric forms [22]. For the sake of simplicity, we consider identical force generating dynein motor domains acting along the flagellum. The total number of motor domains in a beating flagellum has been estimated to be ≃10^{5} [31,32]. The stall force has been found in the range *f*_{0}≃1−5 pN [33,34]. Following reference [15], we choose the characteristic unbinding force for dynein such that *η*≃0.14 and speeds at zero load in the range *v*_{0}≃5–7 μm s^{−1} [33,35]. The characteristic time scale of dynein kinetics sets the order of magnitude of the beating frequency in our model. We choose an estimated value of *τ*_{0}≃50 ms [24,35], which corresponds to a frequency of ≃10 Hz, comparable to the case of human sperm [4]. Finally, we need to estimate *ρ* and *N*. Considering the length of the human sperm flagellum (*L*≃50 μm), we obtain ≃2×10^{3} motors μm^{−1}. In our description, we divide the axoneme into two regions with corresponding dynein teams. Therefore, we have *ρN*≃10^{3} motors μm^{−1}. In order to find *ρ*, we need to choose a criterion to decide the typical length scale or ‘tug-of-war’ length *l*_{c}=*ρ*^{−1} in our coarse-grained description. The typical length scale in the system is given by *l*_{c}∼1 μm and therefore *ρ*∼1 μm^{−1} and *N*∼10^{3}. Hence, we obtain Sp≃5−20, *μ*≃50−100, *μ*_{a}∼10^{3} and *ζ*∼1. The motor activity is studied in a broad range *μ*_{a}≃(2−6)×10^{3} because it plays the role of the main control parameter in our study.

## 4. Linear stability analysis

In this section, we perform a linear stability analysis of equations (2.10) and (2.11). The non-moving state is characterized by *ϕ*=0 and *x*-axis, and the number of plus and minus bound motors is constant in space and time. For the linear stability analysis, we consider the perturbed variables around the base state as *ϕ*=*δϕ* and *n*_{±}=*n*_{0}+*δn*_{±}. Introducing the modulation *δn*≡*δn*_{+}=−*δn*_{−} around *n*_{0} and considering *σ* is a complex eigenvalue and c.c accounts for complex conjugate. From equation (4.1), we get *χ*′(*σ*) is a complex response function:
*χ*(*σ*) is a second complex response function:
*σ* and are equivalent to results presented in reference [24]. With the ansatz *q*_{i}, *i*=1,…,4 and the eigenfunctions read
*δN*=|*χ*′| and *δn* is the same as for *ϕ* except for a phase shift *Δθ* and an overall change in the amplitude *δN*, which depends on *χ*′(*σ*). This result indicates the presence of a time delay between the action of motors and the response of the flagellum. Time delays commonly arise in systems where molecular motors work collectively [36]. Indeed, the regulation of active forces by the time delay of the curvature was proposed as a mechanism to generate travelling bending waves [9]. In order to find *j*=1,…,4. The procedure to obtain the set of boundary conditions for a clamped head is detailed in the electronic supplementary material. The boundary conditions read *σ*_{n}, with the corresponding growth rates λ_{n}=Re[*σ*_{n}] and frequencies *ω*_{n}=Im[*σ*_{n}], which satisfy the boundary conditions, where _{n+1}<λ_{n}, such that the first one has the largest growth rate λ_{1}. Defining ** u**=(

*ϕ*,

*δn*)

^{T}, the general solution of the system reads

_{n}<0, ∀

*n*, solutions decay exponentially to the non-moving state. On the other hand, when λ

_{1}becomes positive, in the range of parameters studied, the system undergoes a Hopf bifurcation and solutions follow an exponential growth, oscillating with frequency

*ω*

_{1}. Next, we study the marginal stable solutions, i.e. when the maximum growth rate equals zero (λ

_{1}=0). For this, we define the critical frequency of oscillation as

*ω*

_{c}≡|

*ω*

_{1}|. Travelling waves propagate from tip to base, a feature already reported for the clamped-type boundary condition [13,14,24]. In figure 2

*c*, the marginal stability curve in phase space is shown. Intuitively, as Sp is increased the travelling instability occurs for higher motor activity

*μ*

_{a}and the critical frequency of oscillation

*ω*

_{c}follows a non-monotonic decrease (see electronic supplementary material, figure S1). For low viscosity (Sp=5), the wave propagation velocity is slightly oscillatory, whereas for high viscosity (Sp=10), it becomes more uniform (figure 2

*a*,

*b*, bottom panels). These results are in agreement with studies on migrating human sperm, where in the limit of high viscosity waves propagated approximately at constant speed [37]. For high viscosity, curvature tends to increase from base to tip, finally dropping to zero owing to the zero curvature boundary condition at the tail (figure 2

*d*and electronic supplementary material). This modulation is consistent with experimental studies on human sperm, which show viscosity modulation of the bending amplitude [37]. In the latter study, however, the effect is more pronounced possibly owing to external elastic reinforcing structures found along the flagellum of mammalian species, as well as other nonlinear viscoelastic effects. Defining

*d*, inset). Similar results can be obtained using other definitions for

*k*, for example using the covariance matrix (see Principal component analysis section).

## 5. Nonlinear flagellar dynamics

In this section, we study the nonlinear dynamics of the flagellum in the limit of small curvature by numerically solving equations (2.10) and (2.11) using a second-order accurate implicit-explicit numerical scheme (see electronic supplementary material). The unstable modes presented in §4 follow an initial exponential growth and eventually saturate at the steady state owing to the nonlinearities in the system. In figure 3, two different saturated amplitude solutions are shown. Figure 3*a* (left) corresponds to a case where the system is found close to the Hopf bifurcation, whereas figure 3*a* (right) corresponds to a regime far from the bifurcation. We note that the marginal solution obtained in the linear stability analysis (figure 2*b*, top panel) gives a very good estimate of the nonlinear profile close to the bifurcation point, although it does not provide the magnitude of *ϕ* or *δn*. Frequencies are approximately 10 Hz and maximum amplitudes are found to be small, around ≃4% of the total flagellum length. However, the oscillation amplitude for high motor activity is more than double in respect to the case of low activity (figure 3*a*). The colour code in figure 3*a* indicates the value of the semi-difference of plus- and minus-bound motors *δn*. Plus-bound motors are predominant in regions of positive curvature (*ϕ*_{s}>0) along the flagellum and vice versa. Despite the low duty ratio of dynein motors [35], approximately 2% bound dyneins along the flagellum are sufficient to produce micrometre-sized amplitude oscillations. The full flagella dynamics corresponding to figure 3*a* are provided in the supplementary material, movies S1 and S2.

In figure 3*b*, the time evolution of *ϕ* and *δn* is shown at *a*, respectively. As mentioned in §4, the tangent angle *ϕ* is delayed with respect to *δn*, and the time delay is not considerably affected by motor activity. Close to the instability threshold, both signals are very similar, because the system is found near the linear regime; however, far from threshold, both signals greatly differ. For high motor activity, both the tangent angle and the fraction of bound dyneins at certain points along the flagellum exhibit cusp-like oscillations (figure 3*b*, right). This behaviour is typical of molecular motor assemblies working far from the instability threshold [38]. Despite the signals S in figure 3*b* (right) being nonlinear, they conserve the symmetry S(*t*+*T*/2)=−S(*t*), where *T* is the period of the signal. This is a consequence of the plus and minus motor populations being identical, a property also found in spontaneous oscillations of motor assemblies [36,38]. Finally, in figure 3*c*, we study how the amplitude and frequency of the oscillations vary with the relative distance from the bifurcation point *ε*, the maximum absolute value of the tangent angle seems to follow a square-root dependence, characteristic of supercritical Hopf bifurcation; however, in the strongly nonlinear regime the curve deviates from this trend. On the other hand, the beating frequency decreases for increasing activity. The only possibility to increase *μ*_{a} keeping the other dimensionless parameters fixed is by increasing *ρN*; hence, the larger each dynein team becomes, the larger the activity. Consequently, the necessary time to unbind a sufficient number of dynein motors to drive the instability increases, leading to a lower beating frequency.

### 5.1. Principal component analysis

In this section, we study the obtained nonlinear solutions using principal component analysis [26,27]. This technique treats the flagellar shapes as multi-feature datasets, which can be projected to a lower dimensional space characterized by principal shape modes. Here we analyse the numerically resolved data following reference [26] to study sperm flagella. We discretize our flagella data with time-points *t*_{i}, *i*=1,…,*p* and *M*=4×10^{3} intervals corresponding to points *s*_{j}=( *j*−1)/*M*, *j*=1,…,*r* along the flagellum, with *r*=*M*+1. We construct a measurement matrix *Φ* of size *p*×*r* for the tangent angle where *Φ*_{ij}=*ϕ*(*s*_{j},*t*_{i}). This matrix represents a kymograph of the flagellar beat. We define the *r*×*r* covariance matrix as *s*_{i}. The covariance matrix C is shown for *μ*_{a}=4300 (figure 4*a*, left) and *μ*_{a}=5600 (figure 4*a*, right). In figure 4*a* (left), we find negative correlation between tangent angles that are a distance λ/2 apart. Hence, a characteristic wavelength λ can be identified in the system, which manifests as a long-range correlation in the matrix C. On the other hand, strong positive correlations around the main diagonals correspond to short-ranged correlations mainly owing to the bending stiffness of the bundle [26]. The number of local maxima along the diagonals in C decreases from *μ*_{a}=4300 to *μ*_{a}=5600, and at the same time λ′>λ. Hence, an increase in motor activity slightly increases the characteristic wavelength while decreasing the number of local maxima, which is related to the characteristic wavenumber. Employing an eigenvalue decomposition of the covariance matrix, we can obtain the eigenvectors *v*_{1},…,*v*_{r} and their corresponding eigenvalues *d*_{1},…,*d*_{r}, such that C*v*_{j}=*d*_{j}*v*_{j}. Without loss of generality, we can sort the eigenvalues in descending order *d*_{1}≥⋯≥*d*_{r}. We find that the first two eigenvalues capture >99% variance of the data. This fact indicates that our flagellar waves can be suitably described in a two-dimensional shape space, because they can be regarded as single-frequency oscillators. Each flagellar shape *Φ*_{i}=[*ϕ*(*s*_{1},*t*_{i}),…,*ϕ*(*s*_{r},*t*_{i})] can be expressed now as a linear combination of the eigenvectors *v*_{k} [26]:
*B*_{k} are the shape scores computed by a linear least-square fit. In figure 4*b* (left), the two first eigenvectors *v*_{1},*v*_{2} are shown for *μ*_{a}=4300. In figure 4*b* (right), the flagellar shape at a certain time (thick solid line) is reconstructed (white line) by using a superposition of the two principal shape modes *v*_{1},*v*_{2} (solid and dashed lines, respectively) and fitting the scores *B*_{1},*B*_{2}. Finally, in figure 4*c,* we show the shape space trajectories beginning with small amplitude eigenmode solutions. While close to the bifurcation the limit cycle is elliptic (figure 4*c*, left), far from the bifurcation the limit cycle becomes distorted (figure 4*c*, right). Elliptic limit cycles were also found experimentally for bull sperm flagella [26]. Hence, as found in §5, motor activity in the nonlinear regime significantly affects the shape of the flagellum when compared with the linear solutions, which only provide good estimates sufficiently close to the Hopf bifurcation.

### 5.2. Bending initiation and transient dynamics

Finally, we study bending initiation and transient dynamics for two different initial conditions, in order to understand the selection of the unstable modes. In figure 5*a*,*b* the spatio-temporal transient dynamics are shown for the case of an initial eigenmode solution corresponding to the maximum eigenvalue (figure 5*a*) and an initial sine perturbation in *ϕ*, with equal constant bound motor densities (figure 5*b*). In case (b), travelling waves initially propagate in both directions and interfere at *t*=*t*′ (figure 5*b*,*d* and electronic supplementary material, movie S3). However, in the steady state, both the eigenmode and sine cases reach the same steady-state solution, despite the sinusoidal initial condition being a superposition of eigenmodes. This result provides strong evidence that the fastest-growing mode is the one that takes over and saturates in the steady state. In figure 5*c,* the transient dynamics for case (b) are shown for plus- and minus-bound dynein populations close to the tail (*n*_{0} and begin oscillating in anti-phase around this value, in a tug-of-war competition.

## 6. Discussion

In this work, we presented a theoretical framework for planar axonemal beating by formulating a full set of nonlinear equations to test how flagellar amplitude and shape vary with dynein activity. We have shown how the nonlinear coupling of flagellar shape and dynein kinetics in a ‘sliding-controlled’ model provides a novel mechanism for the saturation of unstable modes in flagellar beating. Our study advances understanding of the nonlinear nature of the axoneme, typically studied at the linear level [13–16].

The origin of the bending wave instability can be understood as a consequence of the antagonistic action of dyneins competing along the flagellum [15,39]. The instability is then further stabilized by the nonlinear coupling between dynein activity and flagellum shape, without the need to invoke a nonlinear axonemal response to account for the saturation of the unstable modes, in contrast to previous studies [23,24]. Moreover, the governing equations (equations (2.10) and (2.11)) contain all the nonlinearities in the limit of small curvature, and they are not the result of a power expansion to leading nonlinear order [19]. Far from the Hopf bifurcation, linearized solutions fail to describe the flagellar shape, and nonlinear effects arise in the system solely owing to motor activity. At the nonlinear level, both the tangent angle and dynein population dynamics exhibit relaxation or cusp-like oscillations at some regions along the flagellum. Similar cusp-like shapes for the curvature have also been reported in sea urchin sperm [40]. This phenomenology is characteristic of motor assemblies working in far-from-equilibrium conditions and has been found in other biological systems such as in the spontaneous oscillations of myofibrils [38,41]. Interestingly, despite the low duty ratio of axonemal dynein, a fraction of approximately 2% bound dyneins along the flagellum is sufficient to drive micrometre-sized amplitude oscillations.

Angular deflections are found to be approximately 0.1 *rad* in the experimentally relevant activity range *μ*_{a}≃(2−6)×10^{3} and the order of magnitude seems not to be crucially affected by viscosity nor other parameters in the system. Hence, despite of the fact that our description provides an intrinsic mechanism for amplitude saturation, it is only able to generate small deflections, typically an order of magnitude smaller than the ones reported for bull sperm flagella [15]. Other structural constraints, such as line tension, are likely to influence the amplitude saturation, owing to the elastohydrodynamic coupling with motor activity (see electronic supplementary material). The presence of tension on the self-organized beating of flagella was previously investigated at leading nonlinear order; however, deflections were also found to be small [19]. Hence, we conclude that a ‘sliding-controlled’ mechanism may not be sufficient to generate large deflections. Our work adds to other recent studies were the ‘sliding-controlled’ hypothesis seems to lose support as the main mechanism responsible for flagellar beating [16,24]. Basal dynamics and elasticity are also likely to influence the amplitude saturation, and substantial further research is still needed to infer whether sliding-controlled regulation is the responsible mechanism behind flagellar wave generation.

Principal component analysis allowed us to reduce the nonlinear dynamics of the flagellum in a two-dimensional shape space, regarding the flagellum as a single-frequency biological oscillator [31]. Note that a two-dimensional description would not hold for multifrequency oscillations, where an additional dimension is required [36]. Interestingly, we found that as activity increases, the characteristic wavenumber of the system slightly decreases. Thus, dynein activity has an opposite effect on wavenumber selection when compared with the medium viscosity (figure 2*d*, inset). We also showed that the steady-state amplitude is selected by the fastest-growing mode under the influence of competing unstable modes, provided that the initial mode amplitudes are sufficiently small.

An important aspect which is not studied explicitly in this work is the direction of wave propagation. For simplicity, we used clamped boundary conditions at the head which are known to induce travelling waves which propagate from tip to base in the sliding-controlled model [13,24]. It is beyond the scope of our study to assess the effects of different boundary conditions and the role of basal compliance at the head of the flagellum, which are known to crucially affect wave propagation [15]. This work is also restricted to the case of small curvatures; however, the full nonlinear equations including the presence of tension could, in principle, be numerically solved as in previous studies [18,19] (see electronic supplementary material). Finally, a real flagellum is subject to chemical noise owing to the stochastic binding and unbinding of dynein motors. Recent studies have provided insights into this problem by investigating a noisy oscillator driven by molecular motors. However, their approach was not spatially extended [31]. Our approach could be suitably extended to include chemical noise in the system through equation (2.11) by considering a chemical Langevin equation for the bound dynein populations including multiplicative noise [42]. In particular, it can be easily deduced from our study that by considering a force-independent unbinding rate, fluctuations of bound motors around the base state have mean *Nη* and variance *Nη*(1−*η*), in agreement with the results in reference [31] where a different model was considered.

The possibility to experimentally probe the activity of dyneins inside the axoneme is one of the most exciting future challenges in the study of cilia and flagella. These studies will be of vital importance to validate mathematical models of axonemal beating and the underlying mechanisms coordinating dynein activity and flagellar beating.

## Data accessibility

Electronic supplementary information for this article includes a supplementary text and three supplementary movies. The data from the supplementary movies are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.qs65q [43].

## Authors' contributions

D.O. carried out analytical work, numerical simulations, data analysis and drafted the manuscript. H.G. and J.C. conceived of the study, coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

J.C. and D.O. acknowledge financial support from the Ministerio de Economía y Competitividad under projects FIS2010-21924-C02-02 and FIS2013-41144-P, and the Generalitat de Catalunya under projects 2009 SGR 14 and 2014 SGR 878. D.O. also acknowledges a FPU grant from the Spanish Government with award number AP-2010-2503 and an EMBO Short Term Fellowship with ASTF number 314-2014. H.G. acknowledges support by the Hooke Fellowship, University of Oxford.

## Footnotes

↵† Present address: Max Planck Institute of Molecular Cell Biology and Genetics, Pfotenhauerstraße 108, 01307 and Max Planck Institute for the Physics of Complex Systems, Nöthnitzerstraße 38, 01187 Dresden, Germany.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3694462.

- Received September 14, 2016.
- Accepted February 2, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.