Nonlinear amplitude dynamics in flagellar beating

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.


Introduction
Cilia and flagella play a crucial role in the survival, development, cell feeding and reproduction of micro-organisms [1]. These lashlike appendages follow regular beating patterns which enable cell swimming in inertialess fluids [2].  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][9][10][11][12][13][14][15][16]. Coordinated beating has been hypothesized considering different mechanisms such as dynein's activity regulation through local axonemal curvature [9][10][11]16], owing to the presence of a transverse force (t-force) acting on the axoneme [12] and by shear displacements [13][14][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][21][22]. Loadaccelerated 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.

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 1c). 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 r ± = r ± (b/2)n, with the orientation of the cross section at distance s along its length defined by the normal vector to the centrelinen = − sin φî + cos φĵ , being φ ≡ φ(s, t) the angle between the tangent vectorŝ ≡ ∂ s r ≡ r s and theî direction (taken along the x-axis  The geometrical constraint of the filament bundle induces an arc length mismatch, (s, t), denoted as sliding displacement: where φ 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 f (s, t) = f (s, t)ŝ as the total internal force density generated at s at time t on the plusfilament 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 −f at the same point. Next, consider that 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

M(s, t) = Mk,
wherek =ĵ ×î, such that M(s, t) reads where F(s, t) = L s f (s , t) ds , provided that ∂ s φ ± ≈ ∂ s φ for bundles characterized by b L. The combined bending stiffness of the filament bundle is given by E b = 2EI, 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 1b). The corresponding bound motor population dynamics reads (2.5) The binding/unbinding rates are given by π ± = π 0 (N − n ± ), ε ± = ε 0 n ± exp(±F ± /f c ) where ε 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 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 wheren ≡ n + − n − and 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 Hereinafter, we switch to dimensionless quantities while keeping the same notation. We nondimensionalize the arc length with respect to the length scale 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 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] where, in our case, the dimensionless internal force density f (s, t) takes the form: 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: (2.10) In the absence of dynein activity, expression (2.10) reduces to the dynamics of an elastic filament bundle with sliding resistance forces [25]. Note that this expression is obtained considering the sliding mechanism and a linear velocity-force relationship for dyneins, but it is independent of dynein kinetics. On the other hand, the dimensionless form of the bound motor population dynamics n ± reads where η ≡ π 0 /(π 0 + ε 0 ) is the duty ratio of the motors andf ≡ f 0 /f c dictates the sensitivity of the unbinding rate on the load. In the electronic supplementary material, we include a table with a summary of all the variables and parameters in the model with their corresponding symbols.

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 thatf = 2. Axonemal dynein is characterized by a low duty ratio estimated to be η 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 (see Linear stability analysis). Using the previous parameters, we get 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.

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 n ± = n 0 ≡ π 0 /(π 0 + ε 0 ef ). This means that the flagellum is aligned with respect to the 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 consideringf ζ φ t 1, we obtain and whereτ ≡ n 0 /η. We use the ansatzes φ =φ(s) e σ t + c.c and δn = δñ(s) e σ t + c.c, where σ is a complex eigenvalue and c.c accounts for complex conjugate. From equation (4.1), we get δñ = χ (σ )φ, where χ (σ ) is a complex response function: Using equation (2.9) and considering f =f (s) e σ t + c.c, we obtainf = χ (σ )φ, where χ (σ ) is a second complex response function: The latter response functions generalize the work in reference [15] for a complex eigenvalue σ and are equivalent to results presented in reference [24]. With the ansatzφ ∼ δñ ∼ e iqs in equation (4.2), we obtain the characteristic equation q 4 −χq 2 +σ = 0, whereχ ≡ μ a χ ,σ ≡ σ Sp 4 andχ ,σ ∈ C. Solving the characteristic equation, we obtain four possible roots q i , i = 1, . . . , 4 and the eigenfunctions read j=1Φ j e iq j s , (4.5) where q j ,Φ j ∈ C. Onceφ is known, δñ(s) = δNφ(s) exp(i θ ), where δN = |χ | and θ = arg(χ ). Therefore, the evolution of δ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 φ(s), we need to impose the four boundary conditions, obtaining a linear system of equations forΦ j , j = 1, . . .  (1). By setting the determinant of the system to zero, we find the set of complex eigenvalues σ n , with the corresponding growth rates λ n = Re[σ n ] and frequencies ω n = Im[σ n ], which satisfy the boundary conditions, where λ n , ω n ∈ R. We order the set of different eigenvalues according to its growth rate λ 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 u(s, t) = n A n φ n δñ n e λ n t e iω n t + c.c, (4.6) where A n ∈ C are free amplitude parameters. For λ n < 0, ∀n, solutions decay exponentially to the nonmoving 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 2c, 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 2a,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 2d 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 k ≡ max |q i |/2π as the characteristic wavenumber, we obtain that it increases almost linearly with Sp (figure 2d, inset). Similar results can be obtained using other definitions for k, for example using the covariance matrix (see Principal component analysis section).

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 3a (left) corresponds to a case where the system is found close to the Hopf bifurcation, whereas figure 3a (right) corresponds to a regime far from the bifurcation. We note that the marginal solution obtained in the linear stability analysis ( figure 2b, 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 3a). The colour code in figure 3a indicates the value of the semi-difference of plus-and minusbound 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 3a are provided in the supplementary material, movies S1 and S2. In figure 3b, the time evolution of φ and δn is shown at s = 3 4 for the cases in figure 3a, 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 3b, right). This behaviour is typical of molecular motor assemblies working far from the instability threshold [38]. Despite the signals S in figure 3b  (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 3c, we study how the amplitude and frequency of the oscillations vary with the relative distance from the bifurcation point ε ≡ (μ a − μ c a )/μ c a . For small ε, the maximum absolute value of the tangent angle seems to follow a squareroot 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.

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 C = (Φ −Φ) T (Φ −Φ), whereΦ = [φ; . . . ;φ] with all rows equal toφ = [φ(s 1 ), . . . ,φ(s r )], wherē φ(s i ) is the mean tangent angle at s i . The covariance matrix C is shown for μ a = 4300 (figure 4a, left) and μ a = 5600 (figure 4a, right). In figure 4a (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 Cv 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 singlefrequency 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]: where B k are the shape scores computed by a linear least-square fit. In figure 4b (left), the two first eigenvectors v 1 , v 2 are shown for μ a = 4300. In figure 4b (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 4c, we show the shape space trajectories beginning with small amplitude eigenmode solutions. While close to the bifurcation the limit cycle is elliptic ( figure 4c, left), far from the bifurcation the limit cycle becomes distorted (figure 4c, 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.

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 5a,b the spatio-temporal transient dynamics are shown for the case of an initial eigenmode solution corresponding to the maximum eigenvalue (figure 5a) and an initial sine perturbation in φ, with equal constant bound motor densities (figure 5b). In case (b), travelling waves initially propagate in both directions and interfere at t = t (figure 5b,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 5c, the transient dynamics for case (b) are shown for plus-and minus-bound dynein populations close to the tail (s = 9 10 ). Both populations decay exponentially with characteristic timeτ to n 0 and begin oscillating in anti-phase around this value, in a tug-of-war competition.

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][14][15][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  (c) Flagellar dynamics in a reduced two-dimensional shape space for μ a = 4300 (left) and μ a = 5600 (right). Elliptic limit cycles are rescaled to better appreciate the distortion owing to the nonlinear terms. Sp = 10, μ = 50, η = 0.14, ζ = 0.  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 2d, 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].