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.
Cilia and flagella play a crucial role in the survival, development, cell feeding and reproduction of micro-organisms . These lash-like appendages follow regular beating patterns which enable cell swimming in inertialess fluids . 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 . 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 . 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 . 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 . This mechanism leads to the propagation of bending undulations along the flagellum, as commonly observed during the movement of spermatozoa .
Many questions still remain unanswered on how dynein-driven sliding causes the oscillatory bending of cilia and flagella . 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  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 . In contrast, dynamic curvature regulation has been recently proposed to account for Chlamydomonas flagellar beating . 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 . 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 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 , with the orientation of the cross section at distance s along its length defined by the normal vector to the centreline , being ϕ≡ϕ(s,t) the angle between the tangent vector and the direction (taken along the x-axis). The subscripts (+) and (−) refer to the upper and lower filaments, respectively (figure 1c). The shape of the bundle is given at any time by the expression: 2.1The geometrical constraint of the filament bundle induces an arc length mismatch, Δ(s,t), denoted as sliding displacement: 2.2where ϕ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 as the total internal force density generated at 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 −f at the same point. Next, consider that N dyneins are anchored at each polar filament in a region lc around s, where lc 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 lc the ‘tug-of-war’ length. We define n±(s,t) as the number of bound dyneins in a region of size lc 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 2.3where ρ≡l−1c 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 , where , such that M(s,t) reads 2.4where , provided that ∂sϕ±≈∂sϕ for bundles characterized by b≪L. The combined bending stiffness of the filament bundle is given by Eb=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.5The binding/unbinding rates are given by π±=π0(N−n±), where ε0 and π0 are constant rates and fc 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 f0 and velocity at zero load v0, the loads are given by F±(s,t)=±f0(1∓Δt/v0). The stall force f0 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 2.6where 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=−Msss, where ζ⊥ is the normal drag coefficient [8,29]. Combining the last expression with equation (2.4), we have 2.7Hereinafter, we switch to dimensionless quantities while keeping the same notation. We non-dimensionalize 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 f0ρ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(ζ⊥/Ebτ0)1/4, μ≡Kb2L2/Eb, μa=bf0ρNL2/Eb and ζ≡b/v0τ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 . 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] 2.8where, in our case, the dimensionless internal force density f(s,t) takes the form: 2.9and Δ=ϕ−ϕ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: 2.10In the absence of dynein activity, expression (2.10) reduces to the dynamics of an elastic filament bundle with sliding resistance forces . 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 2.11where η≡π0/(π0+ε0) is the duty ratio of the motors and 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.
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 . The bending stiffness of the filament bundle has been reported to be Eb≃0.9×10−21 N m2 for sea-urchin sperm [4,25] and Eb≃1.7×10−21 N m2 for bull sperm . 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 , thus K≃2×103 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 . 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 ≃105 [31,32]. The stall force has been found in the range f0≃1−5 pN [33,34]. Following reference , we choose the characteristic unbinding force for dynein such that . Axonemal dynein is characterized by a low duty ratio estimated to be η≃0.14 and speeds at zero load in the range v0≃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 . Finally, we need to estimate ρ and N. Considering the length of the human sperm flagellum (L≃50 μm), we obtain ≃2×103 motors μm−1. In our description, we divide the axoneme into two regions with corresponding dynein teams. Therefore, we have ρN≃103 motors μm−1. In order to find ρ, we need to choose a criterion to decide the typical length scale or ‘tug-of-war’ length lc=ρ−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 lc∼1 μm and therefore ρ∼1 μm−1 and N∼103. Hence, we obtain Sp≃5−20, μ≃50−100, μa∼103 and ζ∼1. The motor activity is studied in a broad range μa≃(2−6)×103 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 . 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±=n0+δn±. Introducing the modulation δn≡δn+=−δn− around n0 and considering we obtain 4.1and 4.2where . We use the ansatzes and , where σ is a complex eigenvalue and c.c accounts for complex conjugate. From equation (4.1), we get , where χ′(σ) is a complex response function: 4.3Using equation (2.9) and considering , we obtain , where χ(σ) is a second complex response function: 4.4The latter response functions generalize the work in reference  for a complex eigenvalue σ and are equivalent to results presented in reference . With the ansatz in equation (4.2), we obtain the characteristic equation , where , and . Solving the characteristic equation, we obtain four possible roots qi, i=1,…,4 and the eigenfunctions read 4.5where . Once is known, where δN=|χ′| and . 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 . Indeed, the regulation of active forces by the time delay of the curvature was proposed as a mechanism to generate travelling bending waves . In order to find , we need to impose the four boundary conditions, obtaining a linear system of equations for , 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 , , and . 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 . 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 4.6where are free amplitude parameters. For λ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 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 . 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 . 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 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).
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 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 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 , 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 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 . 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 . For small ε, 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  to study sperm flagella. We discretize our flagella data with time-points ti, i=1,…,p and M=4×103 intervals corresponding to points sj=( 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=ϕ(sj,ti). This matrix represents a kymograph of the flagellar beat. We define the r×r covariance matrix as , where with all rows equal to , where is the mean tangent angle at si. 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 . 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 v1,…,vr and their corresponding eigenvalues d1,…,dr, such that Cvj=djvj. Without loss of generality, we can sort the eigenvalues in descending order d1≥⋯≥dr. 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=[ϕ(s1,ti),…,ϕ(sr,ti)] can be expressed now as a linear combination of the eigenvectors vk : 5.1where Bk are the shape scores computed by a linear least-square fit. In figure 4b (left), the two first eigenvectors v1,v2 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 v1,v2 (solid and dashed lines, respectively) and fitting the scores B1,B2. 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 . 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 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 (). Both populations decay exponentially with characteristic time to n0 and begin oscillating in anti-phase around this value, in a tug-of-war competition.
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 . 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 . 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)×103 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 . 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 . 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 . Note that a two-dimensional description would not hold for multifrequency oscillations, where an additional dimension is required . 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 . 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 . 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 . 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  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.
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 .
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.
We have no competing interests.
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.
↵† 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.