## Abstract

Cross-linked filament bundles, such as in cilia and flagella, are ubiquitous in biology. They are considered in textbooks as simple filaments with larger stiffness. Recent observations of flagellar counterbend, however, show that induction of curvature in one section of a passive flagellum instigates a compensatory counter-curvature elsewhere, exposing the intricate role of the diminutive cross-linking proteins at large scales. We show that this effect, a material property of the cross-linking mechanics, modifies the bundle dynamics and induces a bimodal *L*^{2} − *L*^{3} length-dependent material response that departs from the Euler–Bernoulli theory. Hence, the use of simpler theories to analyse experiments can result in paradoxical interpretations. Remarkably, the counterbend dynamics instigates counter-waves in opposition to driven oscillations in distant parts of the bundle, with potential impact on the regulation of flagellar bending waves. These results have a range of physical and biological applications, including the empirical disentanglement of material quantities via counterbend dynamics.

## 1. Introduction

The spontaneous generation of harmonic bending waves along a sperm flagellum has been a source of fascination since it was reported on in the late seventeenth century [1]. It was not until 1968, however, that the fundamental mechanism behind the flagellar wave propagation was unveiled [2–4]. ATP-induced inter-microtubule tangential motion is converted into transversal forces that are capable of bending the flagellar assembly altogether, laying the empirical basis for the sliding filament theory for eukaryotic flagellum.

Notably, almost one decade before the discovery of the interfilament sliding [2,3,5], the existence of such active elements along the sperm flagellum was theorized via a simple fluid–structure interaction model [6]. Machin demonstrated that the combined action of viscous and elastic dissipation experienced by a slender filament rapidly damps any driven oscillation along its length, thus requiring the action of contractile elements in order to sustain large waving amplitude [7]. Later, Wiggins & Goldstein [8] elegantly demonstrated that the elastohydrodynamics of any simple Euler–Bernoulli filament moving in a viscous fluid leads to a hyperdiffusive dissipation of bending, characterized by a bending penetration length ℓ_{ω}, which can be further exploited to extract material parameters from biological filaments in a wide range of length scales [9]. Hitherto, the elastohydrodynamics of active and passive filaments have generated a vast literature of analytical, computational and empirical studies across disciplines [10–19].

Despite the inherent complexity of filament bundles [20–31], as exemplified by the axonemal flagellum [23,24], with its 9 + 2 cross-linked microtubule doublets arranged in a cylindrical fashion [23], the textbook elastic bending stiffness has been estimated using a simplistic linear relation between bending moment and curvature [27–31], as derived from Euler–Bernoulli rod theory [32]. Incidentally, the inadequacy of classical rod theories, from Euler–Bernoulli to Timoshenko–Cosserat [32], emerged via paradoxical counterbend empirical responses, figure 1*a*, first observed by Lindemann and colleagues [33,37], and later captured via a geometrically exact mechanical model by Gadêlha *et al.* [34] (figure 1*a*). These studies revealed how the induction of curvature in one section of a passive sperm flagellum instigates compensatory counter-curvature elsewhere, namely the counterbend phenomenon [33,37]. They established the critical role of the diminutive elastic linking-proteins while instigating large-amplitude deformations, inherently coupling distant parts along the bundle assembly, despite their small slenderness ratio [33,34,37]. More recently, the counterbend phenomenon was also exploited in order to extract material quantities from *Chlamydomonas* flagella [38], despite the relatively short-length flagella. The dynamical response of the counterbend phenomenon in passive cross-linked bundles still remains unexplored in the literature.

The discovery of the counterbend phenomenon highlighted the current need to reassess both the established material measurements [20–22,25–31,35] and the resulting mechanical response, from statics to dynamics, of cross-linked filament bundles [27–31]. The former is crucial in a broad range of biological structures, from the cytoskeleton of eukaryotic cells to cellular division, cross-bridge muscle contraction and locomotion, via structures like the axoneme. A fundamental challenge, both experimentally and theoretically, is, therefore, to understand how this complex structure yields bulk material properties and overall cellular mechanical responses and, ultimately, function. In active bundles, the consequences of using inadequate material parameters, which have been used for the past 30 years to investigate flagellar waves [7,11,38–49], are still unknown. This is further confronted with an increasing number of, repeatedly contradicting, active control models for the flagellar wave coordination [7,11,36,40,41,44,45,47–52]. Paradoxically, in order to induce bending waves, flagellar control models rely on the implementation of filament-bundle deformations, in distinct material directions, that are yet to be scrutinized in isolation, from curvature [39,41,46,49] to interfilament sliding [7,11,36,45,52,53], and axial distortions [44,47,50]. This is aggravated by the strong coupling between the unknown activity, and the passive and dissipative components, leading to the non-identifiability of parameters when contrasted against experiments [54,55]. Without the disentanglement between the passive and active elements, and without the rationalization of the resultant mechanical response of cross-linked filament bundles, it is unclear, for example, which competing flagellar control hypothesis [7,11,36,40,41,44,45,47–51], if any, is able to provide a quantitative understanding of the flagellar regulation and, crucially, function of the internal mechanics and structure. Indeed, any comprehensive model of flagellar bending self-organization depends on reliable measurements of mechanical and material properties of the system in the absence of activity [34,38].

Here, we complement the seminal work by Machin [6] and Wiggins & Goldstein [8] on the dynamics of passive filaments, and demonstrate how the nanometric cross-linking proteins that are present in passive cross-linked filament bundles instigate novel dynamical counterbend phenomena. This is in contrast with previous models on flagellar wave coordination [7,11,36,40,41,44,45,47–52], which incorporate the cross-linking interaction in conjunction with molecular motor dynamics. We consider the dynamical situation in which only the structural passive elements are present. For axonemal filament bundles, this corresponds to the empirical situation in which molecular motors are rendered passive [33,37,38]; figure 1*a*. The filament-bundle elastohydrodynamical model unveils the occurrence of counter-travelling waves in distant parts of bundle, reducing the propulsive potential of driven oscillations, and even reversing the propulsive direction, from pushing to pulling hydrodynamics. We show that the interplay between the interfilament sliding at the base and cross-linking dissipation elsewhere gives rise to a bimodal *L*^{2} − *L*^{3} length-dependent material response that departs from canonical Euler–Bernoulli theory. Hence, the use of simpler rod theories to analyse experiments can result in paradoxical interpretations. Furthermore, the counterbend dynamics offers a robust way to measure material quantities empirically, bypassing cumbersome force-displacement experiments at the microscale [30,31,33,37,38]; figure 1*a*. These results further suggest that the dynamical counter-wave phenomena is likely to play a critical role on the waveform organization, and the subsequent wave direction, of long flagella [48,49,56].

## 2. Cross-linked filament-bundle elastohydrodynamics

We consider a planar representation of cross-linked filament bundles and flagellar axonemes, as depicted in figure 1*b*, used interchangeably hereafter, composed of two elastic, inextensible filaments that resist deformation with an elastic bending modulus *E* [7,34–36]. Each constituent filament is separated by a distance *a*, much smaller than the filament length *a* ≪ *L*, normal to the centre line **r**(*s*, *t*) at every point in arclength *s* and time *t*. Geometry constrains the normal vector to the plane, where *α* is the angle between the fixed frame *x*-axis and the tangent to the centre line . Like a rail-track [35], the constituent filaments travel distinct contour lengths forcing a geometrical arclength mismatch Δ(*s*, *t*) = Δ_{0}(*t*) + *a*(*α*(*s*, *t*) − *α*_{0}(*t*)), where Δ_{0} and *α*_{0} are the length mismatch and tangent angle at *s* = 0 (figure 1*b*). Points of equal contour length along the filament bundle are connected by elastic links which generate a shearing force, and thus an internal moment, proportional to the sliding displacement *f*(*s*, *t*) = *k*Δ(*s*, *t*) with an elastic sliding resistance *k*. At the basal end, the additional connecting compliance across the filaments, commonly found in spermatozoa and inhomogeneous bundles, is Hookean with a spring constant *κ*_{e} [34,36] (figure 1*b*).

For asymptotically slender filament bundles, the hydrodynamic forces experienced by an infinitesimal element is anisotropic and linearly related to the local velocity , where *ζ*_{⊥}, *ζ*_{∥} are the lowest order resistive coefficients derived from inertialess hydrodynamics [27,28]. Contact forces are not defined constitutively due to the inextensibility constraint [32]. The filament-bundle elastohydrodynamics is governed by the balance of contact forces and contact moments
2.1simplified here for small curvatures [8,9,11,36]. The filament-bundle shape is given by the initial value problem for an arclength *s* and time *t*. Boundary conditions ensure the total balance of forces, *F*(*s*) = − *E**α*_{ss} + *af*(*s*, *t*), and torques, , acting on the bundle [32,34], as detailed in the electronic supplementary material, S1. The resulting cross-linking mechanics couples distant parts along the bundle via the total momentum balance, now modified non-locally by
(the electronic supplementary material, S1). Dissipation from different material directions are mediated by the hydrodynamic drag. The filament-bundle dynamics is dictated by the interplay between the elastohydrodynamic hyperdiffusion [8,9] and the cross-linking diffusion [7,35,36]. These boundary moments alter the hyperdiffusion balance in equation (2.1) and instigate novel long-range phenomena as we explore below.

### 2.1. The counterbend dynamics: angular actuation

The post-transient behaviour of the shape dynamics is captured by single frequency solutions of the form . After convenient rescaling, the eigenvalue problem reduces to *r*^{4} − *μ**r*^{2} − i Sp^{4} = 0, with eigenfunctions , where the coefficients *C*_{j} are determined by non-local boundary moments (the electronic supplementary material, S1). The dimensionless filament-bundle compliance parameter, also referred to as sperm number, *Sp* = *L*(*ζ*_{⊥}*ω*/*E*)^{1/4}, captures the battle between elastic and viscous forces [8], while the sliding resistance parameter, *μ* = *a*^{2}*L*^{2}*k*/*E*, compares bending stiffness with the cross-linking resistance [34,36]. The basal compliance is given by *γ* = *kL*/(*kL* + *κ*_{e}), and varies from *γ* = 0, corresponding to no interfilament sliding at the base, to *γ* = 1, for free basal sliding [34,36].

The emergence of the non-local, counterbend dynamics is depicted in figure 2 for a sinusoidal angular actuation of the proximal end (electronic supplementary material, Movie S1) with **r**_{t}(0, *t*) = 0, and zero force and torque condition at the distal end (the electronic supplementary material, S1). An angular amplitude of 0.4362 rad is used to limit the maximum radius of curvature to 0.1. Figure 2 contrasts Machin's original solutions, figure 2*a*,*b*, with the filament-bundle post-transient dynamics, figure 2*c*,*d*. Travelling waves originating from the distal end, indicative of the non-local counterbend effect [34], can be clearly seen in figure 2*c*,*d*. A relatively small sliding displacement Δ (overlaid colour in figure 2*c*,*d*), equivalent to only half of the bundle diameter, is capable of deforming the bundle non-locally with the same magnitude of the imposed actuation.

The amplitude modulation *A*(*s*), phase *Φ*(*s*) and velocity of propagation *v*_{p}(*s*) = 1/∂_{s}*Φ* are depicted in figure 3, following suitable transformation to solutions of the form *α*(*s*, *t*) = *A*(*s*)cos(*t* − *Φ*(*s*)). The abrupt change in wave direction is triggered by the loss in monotonicity of the phase. Non-local counter-waves propagate in the opposite direction from the distal end with a non-uniform decaying magnitude along the arclength due to the high-order dissipation. This is in contrast with the one-directional wave propagation of Euler–Bernoulli filaments, *μ* = 0 in figure 3. The sharp change in wave direction coincides with reduced wave amplitudes at this point in arclength (figure 3). Interestingly, the non-local actuation of cross-linking moments at distal parts of the bundle is delayed by the overdamped dynamics, as indicated by the proximal–distal phase difference. Higher *Sp* causes larger proximal–distal phase mismatch, faster decay of the counterbend wave speed and amplitude, indicative of a destructive interference between proximal and distal waves in figures 2*d* and 3*b*, demonstrated by the sharp jump in phase in figure 3.

The counterbend dynamics impacts significantly the resulting hydrodynamic propulsion. The time-averaged propulsive force exerted by fluid on the filament bundle is obtained by integrating the *x* component of the local drag force, , along the length of the bundle and averaging over one period of oscillation [8,9,14,15,57], as detailed in the electronic supplementary material, S1. Owing to the symmetry of the oscillations around the *x*-axis, the *y* component of the propulsive force averages to zero. Following [8,15], we conveniently rescale length by (*E*/*ω**ζ*_{⊥})^{1/4} and time by 2*π*/*ω*, so that the averaged propulsive force in the *x*-direction may be written as
with a common elastohydrodynamic factor multiplying the force scaling function *ϒ*_{x}(*Sp*, *μ*, *γ*), which is now modified by the cross-linking mechanics. The force scaling function *ϒ*_{x}(*Sp*, *μ*, *γ*) depicted in figure 4 captures the comparative propulsive potential between different points in the parameter space (*Sp*, *μ*, *γ*), where is the force exerted by the fluid on the filament bundle [8,15,57]. For a simple Euler–Bernoulli filament (*μ* = 0), dashed curve in figure 4, the propulsive force is always positive, except for effectively stiff filaments (*Sp* = 0), in which no propulsion can be generated *Υ*_{x} = 0 [8,15,58]. This is in accordance with Purcell's scallop theorem, where no net propulsion can be achieved via time-reversible motion in inertialess fluids [8,15,58]. As *Sp* increases, the filament becomes effectively more flexible, breaking Purcell's scallop theorem, thus inducing an increasingly large hydrodynamic propulsion, characterized by the emergence of a maximum value (dashed curve in figure 4). Increased filament floppiness, large *Sp*, causes the hydrodynamic force to decay and plateau. In this limit, the main contribution of the propulsive force arises from the tangential component of the hydrodynamic drag. Floppy filaments are unable to induce sizeable deformations, causing the filament to move effectively back and forth tangentially, along its length in the *x* direction. The sweet spot of the hydrodynamic propulsion depicted in figure 4 (dashed curve) is the signature of the competition between viscous drag and bending stiffness, captured by the so-called penetration length ℓ_{ω} = (*E*/*ω**ζ*_{⊥})^{1/4} [8,9,15]; first predicted by Wiggins & Goldstein [8] and empirically observed by Yu *et al.* [15].

A very distinct scenario is observed for filament bundles. As the basal compliance becomes stiffer, by reducing *γ*, the cross-linking mechanics switch from a mostly local contribution with small counterbend deformations (*γ* = 1 in figure 4), to non-local effects, with increasingly large amplitudes at the distal end (*γ* = 0 in figure 4, inset). Ultimately, this causes the propulsive force to vanish (circles in figure 4), and even switch the propulsive direction (negative values of *ϒ*_{x}), thus equivalent to a backward net motion. Imposed oscillations are now counteracted by waves travelling in opposition at distal parts of the bundle (compare top and bottom net plots of each inset indicated by an asterisk in figure 4). Small values of *γ* increases the non-local counterbending amplitude and causes the propulsive force to switch direction (negative values) as *Sp* increases (*γ* = 0, 0.4 and 0.6 in figure 4). However, similarly to the Euler–Bernoulli case, increased floppiness of the filament bundle, i.e. large *Sp*, causes the amplitude of deformation to decay (both amplitude of imposed oscillation and counter-waves), reducing in this way the magnitude of the propulsive force which is now characterized by the appearance of a minimum value (see *γ* = 0 in figure 4).

The zero values of *ϒ*_{x} for *γ* = 0, 0.4 and 0.6, circles in figure 4, are therefore a genuine manifestation of the perfect hydrodynamic balance between the imposed oscillation at the proximal end and the counter-waves at the distal end, with amplitudes modulated by the effective bundle stiffness *Sp*. Increased bundle floppiness causes the counter-wave amplitude to vanish, thus the change in sign to positive values of the propulsive force before it plateaus to a constant positive level, similarly to the Euler–Bernoulli filaments. The lower plateau value for filament bundles is due to the added effective stiffness arising from the cross-linking mechanics. The separatrix in figure 5 captures the region in parameter space where the local extrema of *ϒ*_{x}(*Sp*_{m}) changes sign, the non-trivial roots of *ϒ*_{x} (circles) in figure 4. This indicates the region where a significant influence of non-local counterbend effect is predicted, and illustrates how the triad (*Sp*, *μ*, *γ*) may be conveniently tuned to achieve zero, forward or backward propulsion (figure 5). Reversal in swimming direction may be achieved by simply increasing the frequency of oscillation for instance, recalling that *Sp* = *L*(*ζ*_{⊥}*ω*/*E*)^{1/4}. It is worth noting that the cross-linking dissipation does not affect the main elastohydrodynamic characteristics, governed by the penetration length ℓ_{ω} = *L*/*Sp* [8], as indicated by the positive asymptote for both Euler–Bernoulli filaments and filament bundles in figure 4.

### 2.2. The bimodal length-dependent relaxation dynamics

The cross-linking mechanics introduces a diffusion-like time scale, *L*^{2}*ζ*_{⊥}/*a*^{2}*k*, with a somewhat weaker *L*^{2} geometrical dependence, in addition to the high-order, hyperdiffusive scaling *L*^{4}*ζ*_{⊥}/*E*. Indeed, the cross-linking resistance *μ* = *a*^{2}*L*^{2}*k*/*E* contrasts the elastohydrodynamic and cross-linking time scales. The cross-linking resistance depends on geometrical aspects of the bundle, as *μ* measures the ratio between the natural cross-linking elastic length relative to total axial length via *μ* = (*L*/ℓ)^{2}. The cross-linking elastic length ℓ has an important biophysical interpretation: it is the dimensional length by which cross-linking effects become prevalent. The cross-linking mechanics become increasingly important when *L* > ℓ, while the opposite is found for or smaller. The latter entails the possibility of studying relaxation counterbend phenomena by only varying the length of the filament bundle, and motivates rescaling *μ* relative to the reference length *L*_{0}, so that *μ*(L) = L^{2}*μ*_{0}, *μ*_{0} = *a*^{2}*L*^{2}_{0}*k*/*E* and now L = *L*/*L*_{0}, with smallest dimensionless length L = 1.

The non-local, counterbend dynamics decaying from initial data, with an amplitude *a*_{n}, is depicted in figure 6 for filament bundles that are clamped at the proximal end (the electronic supplementary material, S1). This is characterized by an effective relaxation constant *λ*^{−4}_{n} via with *S*_{n} = *C*_{1} sin(*q*_{1n}*s*) + *C*_{2} cos(*q*_{1n}*s*) + *C*_{3} sinh(*q*_{2n}*s*) + *C*_{4}cosh(*q*_{2n}*s*) for a given participating mode *n*. The triad of dissipative contributions acts as an effective dispersion medium for both bending and cross-linking deformations, dictated by the same dispersion relation. However, the mode shape is captured by the wavenumber–eigenvalue coupling for *l* = 1, 2, which is not only influenced by the effective relaxation constant, reminiscent of pure high-order elastohydrodynamic dissipation, but also by the cross-linking diffusion. Boundary conditions define the transcendental solvability condition for *λ*_{n}, which depend implicitly on both *μ* and *γ*, with an infinite number of mode solutions for each parameter set. Curve-fitting expressions for *λ*_{1} obtained from numerical solutions are presented in the electronic supplementary material, S1.

The non-local cross-linking diffusion introduces a bimodal length-dependent material response, as illustrated in figure 6 for the relaxation time of the fundamental mode
The relaxation dynamics of Euler–Bernoulli filaments decay with a characteristic *L*^{4}-dependence associated with the filament elastohydrodynamic hyperdiffusion [8]. When *μ* = 0, *λ*_{1} is a constant solely defined by boundary conditions, thus *τ*_{1} ∝ *L*^{4} as discussed above. This canonical case is shown by the black straight curve in figure 6, also indicated by the ‘∝*L*^{4}’ at the top right corner of figure 6. The relaxation dynamics of filament bundles, on the other hand (depicted by the light blue and dark blue curves), deviate from the canonical *L*^{4}-dependence (straight black curve) as *L* increases, or equivalently *μ* (the top horizontal axis in figure 6 shows the associated *μ* values via the relationship *μ*(*L*) = *L*^{2}*μ*_{0}). This is due to the fact that, for filament bundles, *λ*_{1} = *λ*_{1}(*L*), thus inducing an additional length dependence on the denominator of *τ*_{1}. This causes the relaxation time to gradually deviate from the classical *L*^{4}-dependence, with greater magnitude as *L* increases, depending on the details of the basal compliance. For *γ* = 0, the deviation from the Euler–Bernoulli case starts approximately from *μ* > 1, while, for *γ* = 1, it occurs for approximately *μ* > 10. An *L*^{2} asymptote arises for long filament bundles with *γ* = 0, while for *γ* = 1 the transition is to an *L*^{3} behaviour instead. These long-length asymptotes are indicated by the dashed grey curves labelled by ‘∝*L*^{3}’ and ‘∝*L*^{2}’ at the bottom left corner of figure 6.

The length-dependent transition between *L*^{2} and *L*^{3} asymptotes for long filament bundles is governed by the basal compliance. This *γ*-dependence is depicted in the inset (*a*) of figure 6, which plots the exponent *ζ* of the relationship *τ*_{1}∝*L*^{ζ} in the limit of large *L*. For asymptotically long filament bundles, the exponent of *τ*_{1}∝*L*^{ζ} is quadratic (*ζ* = 2), and remains nearly quadratic until *γ* approaches 1, at which *ζ* = 3. Such bijection of the material response entails that simultaneous measurement of the bundle mechanical properties, in different material directions, can be extracted from simple relaxation experiments. In particular, increased inter-filament sliding, concentrated towards the clamped end, induces curvature reversal for long filament bundles (figure 6 inset (*b*)), reminiscent of the counterbend phenomenon [34]. Boundary conditions require zero contact forces and torques at the free end, while clamped constraint facilitates the accumulation of cross-linking sliding towards the proximal end. Both bending and cross-linking deformations relax towards the reference configuration with the same effective rate *λ*_{1}. Figure 6 inset (*b*) shows the simultaneous hyperdiffusive relaxation of bundle curvature (bundle shape) and sliding displacement (overlaid colour) towards the equilibrium straight configuration. Figure 6 thus summarizes the impact of the cross-linking mechanics on the relaxation time across the parameter space (*γ*, *μ*).

## 3. Discussion

We studied the transient and post-transient dynamics of overdamped filament bundles that are interconnected by linking elastic proteins. Deformations in distinct material directions, arising from the cross-linking interfilament sliding and pure bending deformation, are coupled with local slender-body hydrodynamics. This leads to an effective dispersion mechanism governed by the superposition of short- and long-range dissipation mechanisms. Cross-linking stresses are transmitted to distant parts of the bundle via boundary balance of moments. The cumulative moments are able to surpass the high-order elastohydrodynamic dissipation, and shape the bundle structure non-locally, with increased influence for long filament bundles, or equivalently, large *μ*.

The delicate interplay between the interfilament sliding at the base and the rest of the bundle results in a bimodal dynamic response, which departs from the classical Euler–Bernoulli theory [6,8,32]. When the basal sliding is permitted, cross-linking diffusion is mostly local, and acts to effectively reinforce the bundle structure. Long-range curvature-reversal events, however, are magnified when the basal sliding is constrained [34]. The counterbend dynamics generate spontaneous travelling waves in opposition to driven oscillations, which are capable of suppressing the propulsive potential, and even reverting the direction of propulsion (figure 4). Curvature perturbations diffuse more rapidly, a hundred times faster than Euler–Bernoulli hyperdiffusion with an equivalently higher bending rigidity (figure 6). Relatively small cross-linking deformations, up to only 30% of the bundle diameter, are capable of exciting large counterbend modes (figure 6 inset (*b*)), and induce a bimodal *L*^{2} − *L*^{3} length-dependent deviation from the *L*^{4}-dependence of canonical filaments. Paradoxical measurements may arise if standard Euler–Bernoulli theory is used to interpret experiments [29,31,33,37], as exemplified by the paradoxical length-dependent bending stiffness in microtubules [59]. Indeed, the length deviation predicted here may be mistakenly interpreted as an effective length-dependent bending rigidity via *L*^{4}*ζ*_{⊥}/*E*(*L*) [8], if rather the Euler–Bernoulli theory is used; de facto, the Euler–Bernoulli theory is traditionally used since the first measurements of flagellar bundles [27,29,31].

Static, force-displacement experiments that are often used to probe flagellar material quantities [29,31,33,38] are cumbersome (figure 1*a*). They require high-precision force-calibrated probes and micro-manipulators, and often rely on the rare attachment of the filament's tip to the coverslip to microprobe actuation [29,31,33]. This proximity of the filament bundle to the coverslip can interfere with the interfilament sliding due to surface adhesion, biasing in this way force and shape measurements [34]. The counterbend dynamics provides a simpler and robust empirical route for the disentanglement of material parameters. This includes measurements of the basal interfilament elasticity, despite being deeply embedded at the connecting piece of the bundle (figure 1). As a result, standard microfluidic designs may be explored to induce shape changes dynamically [18]. Likewise, the dynamical counter-wave phenomenon may also inspire the design of artificial swimmers [60] that are able to reverse the swimming direction by simply increasing, for instance, the frequency of oscillation (figure 4).

The counter-wave phenomena becomes increasingly important for bundles longer than , typically 5 μm for flagella [2,7,36,49,61,62]. Interestingly, the majority of eukaryotic flagella exceed ℓ by a few orders of magnitude, from approximately 30 μm for *Chlamydomonas* and sea urchin sperm to almost 200 μm for quail sperm [61,62]. Cross-linking effects may also become increasingly important during flagellar growth, and influencing in this way the wave coordination during flagellar reconstitution in *Chlamydomonas* [63]. Molecular motors organization thus may operate differently for *L* < ℓ and *L* > ℓ. Indeed, local flagellar control models [7,49] recently gained empirical support when tested against short flagella experiments [49], a regime where counterbend phenomenon may be negligible (*L* < ℓ). This is despite of the well-known negative support of curvature control models [36,41], tested instead against long flagella (*L* > ℓ). The recurrent contradictions between sliding and curvature control models [7,11,36,41,49,52] may suggest the occurrence of distinct length-dependent flagellar regimes.

Linear models coupling the molecular motor reaction kinetics with interfilament sliding spontaneously propagate waves along the axoneme via a Hopf bifurcation [11,36,42,45,48]. The resulting wave train is observed to move from tip to base, i.e. in the direction opposite to what one would expect from a local dissipation theory (when the interfilament sliding is prevented at the base) [11,49,56]. Incidentally, the basal compliance was observed to influence the direction of wave propagation in flagella self-organization models [36,56], demonstrating the sensitivity of the direction of the travelling wave to details of the connecting piece (basal part) and boundary conditions. This is in agreement with the bimodal response predicted here (figure 2), in which counterbend is maximized when *γ* = 0 [34]. The wave direction is influenced non-locally by cross-linking effects whose magnitude is regulated by the basal interfilament sliding (figures 2 and 6). This might explain the surprising significance of the basal compliance during the flagellar wave coordination observed recently in empirical studies [61], and even flagellar synchronization that may arise without recurring to hydrodynamical coupling [64]. Nevertheless, the mechanisms by which the flagellar wave direction is selected is poorly understood. Previous studies were reduced to the linear level, and at the nonlinear regime, the dynamical instability generates unstable travelling waves that can propagate in both directions, with potential for multi-frequency modes [54,65]. The boundary conditions and basal mechanics assist the mode selection nonlinearly, and thus the direction of propagation, emphasizing how the flagellar dynamics is critically dependent on the underlying structural mechanics of the axoneme. Non-local hydrodynamic interactions [62], transversal axonemal deformations [44,47] and geometrical non-linearities [10,39] are also likely to affect the emergence of self-organization in flagellar systems.

The high-order diffusive interaction, intrinsic to elastohydrodynamic systems in equation (2.1), is observed throughout nature. In non-dilute systems, particles are affected by density variations beyond the nearest neighbours via, for example, biharmonic interactions *u*_{t} = *D*_{1}∇^{2}*u* − *D*_{2}∇^{4}*u* [66,67], thus closely related to equation (2.1). Despite the relatively short-range influence, such higher-order diffusion instigates complex spatio-temporal dynamics and self-organization in all fields of science [68,69]. They drive instabilities and even mediate the coexistence of spatial patterns and temporal chaos, as observed in Kuramoto–Sivashinsky systems [69]. Other exemplars of local, higher-order diffusion are found in Ginzburg–Landau superconductors, spatial patterning in Cahn–Hilliard and biochemical systems, plus generalized Fisher–Kolmogorov models, water waves and continuum mechanics systems, among others [66,68,69].

In contrast with canonical high-order diffusion systems [66,68,69], the flagellar scaffold, or equivalently, any cross-linked filament-bundle immersed in a viscous fluid is governed by high-order dispersion medium that is inherently *non-local*. Here, the biharmonic diffusion arises instead via a *local* elastohydrodynamic dissipation [6,8], while the Fickian-like interaction arises through the *long-range* coupling reflecting the bundle mechanics [7,11,34,35], effectively connecting distant parts of the system via boundary bending moments. This unveils the potential for rich long-range phenomena via reaction–diffusion interactions [11,54,56], from non-local pattern formation to long-range synchronization of auto-oscillators, that are yet to be fully explored in the realm of mathematical biology.

We hope that these results will inspire theoreticians and experimentalists to study the dynamical effects of the counterbend phenomenon in the filament bundle as found throughout nature, including prospects for counterbend reaction–diffusion systems in flagellar dynamics, effectively bridging, non-locally, dynamical systems and PDEs.

## Data accessibility

This work is based on a theoretical development and generated no data.

## Author's contributions

R.C. carried out the mathematical modelling, analysed the result and drafted the manuscript; H.G. carried out the modelling, conceived the study, designed the study, coordinated the study and drafted the manuscript. All authors gave their final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

R.C. thanks Cambridge Bridgwater Summer Research Programme. H.G. acknowledges support by the Hooke Fellowship, University of Oxford, and WYNG Fellowship, Trinity Hall, Cambridge.

## Acknowledgments

The authors thank Dr E.A. Gaffney for enlightening discussions. We dedicate this work in memory of Prof. John R. Blake, whose work and devotion will continue to inspire future generations of scientists.

## Footnotes

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

- Received January 27, 2017.
- Accepted May 3, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.