Ciliary beating patterns map onto a low-dimensional behavioural space

Biological systems are robust to perturbations at both the genetic and environmental levels, although these same perturbations can elicit variation in behaviour. The interplay between functional robustness and behavioural variability is exemplified at the organellar level by the beating of cilia and flagella. Cilia are motile despite wide genetic diversity between and within species, differences in intracellular concentrations of ATP and calcium, and considerable environment fluctuations in temperature and viscosity. At the same time, these perturbations result in a variety of spatio-temporal patterns that span a rich behavioural space. To investigate this behavioural space we analysed the dynamics of isolated cilia from the unicellular algae Chlamydomonas reinhardtii under many different environmental and genetic conditions. We found that, despite large changes in beat frequency and amplitude, the space of waveform shapes is low-dimensional in the sense that two features account for 80% of the observed variation. The geometry of this behavioural space accords with the predictions of a simple mechanochemical model in the low-viscosity regime. This allowed us to associate waveform shape variability with changes in only the curvature response coefficients of the dynein motors. The beating of motile cilia arises from the collective action of hundreds of proteins. A study of the dynamics of cilia under different environmental and genetic conditions shows that the space of beating variations is low-dimensional.

B iological systems can function despite genetic and environmental perturbations in their molecular components. These perturbations, in turn, affect behaviour at higher levels, from organelles through cells to whole organisms. A key finding is that behavioural spaces, the spaces that embed the whole repertoire of behaviours, tend to be low-dimensional. Examples are the shapes of moving nematodes 1,2 , the swimming trajectories of ciliates and bacteria 3,4 , the postures of walking flies 5 and mice 6 and the activities of Hydra 7 , which can all be reduced to low-dimensional behavioural spaces. An important conclusion of these analyses is that only a few features define each behaviour, although there can be considerable variation of these features among individuals within that behaviour. A challenge is to map these behavioural spaces onto underlying molecular and biophysical mechanisms 8 .
Cilia are motile organelles that undergo complex shape changes. These changes drive motion relative to the surrounding fluid and allow cilia to respond to external cues. Cilia, therefore, display a simple form of behaviour that makes them a potentially good model system for relating behaviour to underlying molecular mechanisms.
A key feature of cilia is that they can beat in the presence of a wide range of genetic mutations of key proteins [9][10][11] and can operate over a wide range of temperatures 12 , viscosities 12,13 , ATP concentrations 14 and buffer conditions such as pH and calcium concentration [15][16][17] . Thus, cilia exemplify conserved function 18 , namely beating, in the face of large genetic and environmental perturbations. In this work we use variations in beat shapes elicited by these perturbations to construct the behavioural space of cilia.
The beating motion of cilia is powered by thousands of dynein motors belonging to several different classes 19 . These motors drive the sliding of adjacent doublet microtubules 20 within the axoneme, the structural core of cilia 21 . Sliding is converted to bending by proteins that constrain the shearing of doublets (for example, nexin links or radial spokes) 20,[22][23][24] and sliding at the base 25 . This biomechanical and structural knowledge forms the basis of mechanochemical models that can successfully recapitulate the oscillatory motion of cilia [26][27][28][29][30][31][32] and synthetic filament bundles 33,34 . The combination of mechanistic models, cited in this paragraph, with phenotypic variability in response to perturbations, cited in the previous paragraph, makes the cilium an ideal model system to study the connection between behavioural spaces and molecular mechanisms.
In this work, we have asked: what is the geometry of the behavioural space of beating cilia subject to a wide range of perturbations, and how do the individual components of this space relate to the underlying mechanochemistry.

Results
Quantification of the ciliary beat. We measured the waveforms of isolated and reactivated Chlamydomonas reinhardtii axonemes ( Fig. 1a) with high temporal and spatial precision (see ref. 32 and Methods). The waveforms, discretized in 25 points along the arc length of the axoneme, were tracked over time for up to 200 beat cycles using filament tracking software 35 (one beat cycle is shown in Fig. 1b).
The periodic beat of the axoneme was parameterized by the tangent angle ψ(s, t) (with s the arc length and t the time) relative to the co-swimming frame 36 (Fig. 1b). The power spectrum of the tangent angle ( Fig. 1d) showed that typically <10% of the power was in the modes n ≥ 2, which we therefore neglected. After subtracting the static mode n = 0, discussed in refs. 32,37 , we parameterized the beat dynamics as a travelling wave: where a(s) is the amplitude profile, φ(s) is the phase profile and f is the beat frequency. The arc length was normalized by the total length L of the axoneme, so that s ∈ [0, 1]. Although a(s) and φ(s) are often approximated by scalars, such as the mean amplitude a 0 ≡ ∫ 1 0 a(s)ds and the wavelength λ ≡ − πL

Ciliary beating patterns map onto a low-dimensional behavioural space
Veikko F. Geyer 1 ✉ , Jonathon Howard 2 ✉ and Pablo Sartori 3 ✉ Biological systems are robust to perturbations at both the genetic and environmental levels, although these same perturbations can elicit variation in behaviour. The interplay between functional robustness and behavioural variability is exemplified at the organellar level by the beating of cilia and flagella. Cilia are motile despite wide genetic diversity between and within species, differences in intracellular concentrations of ATP and calcium, and considerable environment fluctuations in temperature and viscosity. At the same time, these perturbations result in a variety of spatio-temporal patterns that span a rich behavioural space. To investigate this behavioural space we analysed the dynamics of isolated cilia from the unicellular algae Chlamydomonas reinhardtii under many different environmental and genetic conditions. We found that, despite large changes in beat frequency and amplitude, the space of waveform shapes is low-dimensional in the sense that two features account for 80% of the observed variation. The geometry of this behavioural space accords with the predictions of a simple mechanochemical model in the low-viscosity regime. This allowed us to associate waveform shape variability with changes in only the curvature response coefficients of the dynein motors.
here we are interested in their arc-length dependence, which defines the shape of the beat. Our goal is to characterize variations in the set of parameters P = {a(s), φ(s), f, L}, which we define as the behavioural space of beating. A subset of these parameters, {a(s)/a 0 , φ(s)}, defines the shape. We explored the variability of P by observing groups of axonemes under several environmental perturbations (temperature, viscosity, different ATP concentrations, Ca 2+ and taxol) and genetic perturbations (oda1, ida3, ida5, mbo2 and tpg1) (Fig. 2). Our dataset comprises 498 axonemes undergoing a total of ~40,000 beat cycles. Under all conditions we found beating axonemes, with a percentage of reactivating axonemes greater than 70%. The range of conditions resulted in a wide range of beat shapes and lengths (Fig. 2a), beat frequencies ( Fig. 2b), amplitude profiles ( Fig. 2c) and phase profiles (Fig. 2d). Therefore, although the beating of cilia is robust against perturbations, the features of its waveform display substantial variations. We now set out to quantify these variations.
Variations in frequency, mean amplitude and wavelength. The beat frequency increases with ATP concentration and temperature and decreases with viscosity (Fig. 3a). This variation encompasses a frequency range of about one order of magnitude, and the between-conditions variance is about seven times larger than the within-condition variance. By contrast, the mean amplitude varies little between these conditions, although it shows large variations within conditions: the between-conditions variance in amplitude is about one-fifth the within-conditions variance. Thus, ATP concentration, temperature and viscosity lead to large changes in frequency with little effect on mean amplitude, in agreement with studies cited in the Introduction.
The mutations oda1 (which lacks outer-arm dyneins 38 ), ida3 and ida5 (which lack inner dynein arms f and a, c, d, e, respectively 39,40 ), mbo2 (which lacks microtubule inner proteins and has a symmetric beat 41 ) and tpg1 (which has reduced polyglutamamylation required for axonemal integrity 42,43 ), as well as addition of the ion Ca 2+ and the microtubule poison taxol, did not appreciably broaden the range of beat frequencies (Fig. 3b, 15 Hz to 160 Hz). We note, however, that the oda1 mutant has a twofold lower beat frequency over all temperatures, as reported in earlier studies (Introduction). The mbo2 mutation and taxol increased the range of mean amplitudes (Fig. 3b). Interestingly, in mbo2 the amplitude is high but the beat frequency is low, whereas in taxol (which has not been studied before) the amplitude is low but the beat frequency is high. This reciprocal variation of amplitude and frequency may reflect the energetics of the beat (Discussion).
The variations in the beat frequency and mean amplitude within individual axonemes are small compared to the axoneme-toaxoneme variation. Over the times of analysis, which were typically ~50 cycles, the variances of the frequency and the mean amplitude were considerably smaller than the axoneme-to-axoneme variance (Supplementary Fig. 1 and Supplementary Table 1). It is nonetheless possible that variability within a given condition is due to long-term dynamic effects, for example due to the sampling of axonemes in different states of rundown. Fluctuations in beat frequency over time are so small that the Q-factor of the oscillations can reach values as high as 150 (Fig. 3b, inset), confirming earlier observations 44 . This Q-factor, one of the highest in any biological system, is exceeded by the tuning of sonar responses in the inner ear of the moustached bat (Q up to 1,000 (ref. 45 )). Thus, the variation in beat frequency and amplitude is not due to short-term variation within individual axonemes, but rather due to differences between axonemes and their response to perturbations.
The lengths of axonemes vary twofold over the entire dataset, from 7.2 μm to 15.4 μm (Fig. 3c). This variation is primarily due to the shorter mbo2 axonemes and the longer oda1 axonemes. Remarkably, the wavelengths under all these conditions are almost equal to the lengths, with λ/L = 0.97 ± 0.08 (mean ± standard deviation for all 498 cilia) (see analysis relating to Fig. 5).

Asymmetry and parabolicity are dominant waveform traits.
To characterize the variation in the amplitude and phase profiles of the tangent angle, we decomposed both into shifted Legendre polynomials of increasing order (Methods and Fig. 4a-c). For the amplitude, the lowest order (a 0 ) equals the mean amplitude, the first order (a 1 ) corresponds to the linear deviation from constant amplitude and is a measure of the proximal-distal asymmetry, and the second order (a 2 ) measures the quadratic deviation. For the phase, we set φ 0 to zero (as it is arbitrary), the first order is φ 1 = −πL/λ by definition of λ (the linear component of the phase profile), and φ 2 corresponds to the quadratic deviation from the line. Because we are interested in shape, we normalized the amplitude coefficients: ā 1 = a 1 /a 0 and ā 2 = a 2 /a 0 . Figure 4d displays the space of variations of parabolicity and asymmetry of the amplitude. Most points cluster around ā 1 ≈ 0.0 and ā 2 ≈ 0.3, which corresponds to a symmetric and convex amplitude profile. The exceptions include taxol, ida5 and tpg1, which have ā 1 < 0 corresponding to decreasing amplitudes towards the distal tip, and Ca 2+ , which has ā 1 > 0 corresponding to increasing amplitude towards the distal tip. The oda1 waveforms have a small parabolicity ā 2 corresponding to their linear amplitude profiles. The ida5 waveforms have smaller parabolicity and negative asymmetry. Thus, despite the robustness of the beat, both genetic and environmental perturbations led to variations in the waveforms.
Parabolicity and asymmetry, the polynomial terms that together contribute 80% of the variance (Fig. 4e), are also the main features of the data: they strongly correlate with the first two principal components a pc 1 and a pc 2 (Fig. 4f, left boxes). Furthermore, the asymmetry of the amplitude profile, ā 1 , is anti-correlated with the parabolicity of the phase profile, φ 2 , whereas the parabolicity of the amplitude profile, ā 2 , is correlated with the asymmetry of the phase profile, φ 1 (Fig. 4f, right boxes). Thus, the dominant shape features are the parabolicity and asymmetry of the amplitude and phase profiles. properties of the axoneme. To explore this possibility, we asked whether a mechanical model of the ciliary beat, based on that in ref. 32 , could account for the diversity of waveforms. The model, described in the Methods, contains four non-dimensional coefficients: β′ and β″, which characterize the motors' dependence on the instantaneous curvature (instantaneous response) and the rate of change of curvature (dynamic response), respectively; k, which is the shear stiffness of the axoneme; and Ma, which is a constant that comes from non-dimensionalizing Machin's equation (see equation (4) in ref. 26 ). Ma = 2πfξ n L 4 /κ, where κ is the flexural rigidity of the axoneme and ξ n is the normal friction coefficient. For Chlamydomonas, Ma ≈ 50, using the parameters in the Methods. We fitted the mechanical model to the data (see example in Fig. 1e,f). The curve-fitting reduced the dimensionality of the waveform description from 48, a(s i ) and φ(s i ) with i = 1, …, 24 points along the length, to just 3, the model parameters (Fig. 5a). The agreement of the model with the data was excellent, with 92% of the axonemes having R 2 > 0.9 (Fig. 5b). Figure 5c shows that one of the parameters, the sliding stiffness k, is strongly correlated with another parameter, the dynamic motor response β″. The observed relationship between these parameters can be recovered analytically by noting that the ratio of friction to bending forces, which is proportional to the Machin number Ma ≡ Ma/(2π) 4 , is ≪1. Taking the limit Ma → 0 gives the curved dashed line in Fig. 5c (Methods). The good match between this 'dry friction' limit and the full theory supports the hypothesis that Chlamydomonas operate in a low-friction regime (Discussion). Therefore, the waveform is controlled by just the two motor response coefficients, which correlate with the asymmetry and parabolicity of the amplitude (Fig. 5c, inset).
In agreement with the above argument, Fig. 5d recapitulates the variability in waveform features observed in Fig. 4d. For example, the taxol and Ca 2+ data lie far apart from each other in the (β′, β″) space, just as in the (ā 1 ,ā 2 ) space. Furthermore, because amplitude and phase are correlated, the motor response coefficients also correlate with the phase features (Fig. 5d, inset). Parameters β′ and β″ not only control the asymmetry and parabolicity of the amplitude, respectively, but β″ also controls the wavelength. Specifically, in the limit Ma → 0, we can show that the curvature wavelength is λ = −4πL/β″, and therefore λ/L = 1, which is typically observed, requires β″ = −4π, in agreement with the average values in Fig. 5c (dotted line) and Fig. 5d. Thus, the space spanned by the motor response coefficients accurately recapitulates the space of amplitude and phase features, which suggests that all genetic and environmental  45 13 28 19 11 15 10 14 8 14 15 14 13 12  perturbations are buffered into changes of the curvature regulation of the motors.

Discussion
The dimensionality of the ciliary beat is low. By using a wide variety of environmental and genetic perturbations, we have constructed the behavioural space of ciliary swimming waveforms. Asymmetry (ā 1 and φ 1 ) and parabolicity (ā 2 and φ 2 ) of the amplitude and phase profiles suffice to describe about 80% of the variation between and within conditions. The dimensionality of swimming behaviour for isolated axonemes is therefore very low: beat frequency, mean amplitude, shape asymmetry, shape parabolicity and axoneme length. Low-dimensional phenotypic spaces have been observed in other behaviours (Introduction), but the range of conditions probed in our study makes the low-dimensionality even more remarkable. It remains open how the dimensionality of isolated swimming axonemes relates to that of intact cells.
The shape space is spanned by the curvature sensitivity of the dyneins. The geometry of the shape space (the amplitude and phase profiles) accords with a simple biomechanical model of ciliary beat, which is an extension of a previous model 32 . In this model, dynein motors respond to curvature with an instantaneous coefficient β′ and a dynamic coefficient β″. The model also includes a passive stiffness to the sliding of doublets, k. Note that the typical value k = 20 (non-dimensionalized from Fig. 5c), which corresponds to 120 pN, is in good agreement with the measured value of 80 pN from ref. 46 . The connection between the shape and biomechanical spaces enabled a remarkable finding: the curvature coefficients β′ and β″ strongly correlate with the shape parameters, ā 1 , ā 2 , φ 1 and φ 2 :ā where '~' denotes correlation. Thus, the different shapes correspond to different curvature sensitivities. Wild-type cilia have symmetric and convex amplitude profiles, ā 1 ≈ 0.0 and ā 2 ≈ 0.3, corresponding to β′ = 0 and β″ ≈ −10. Using equation (2), we can ascribe deviations from this typical behaviour to changes in motor properties. For example, taxol, tpg1 and ida5 have ā 1 < 0, corresponding to a beat whose amplitude decays distally. Despite the different nature of these perturbations, we find that all three perturbations decrease the instantaneous response coefficient, β′, and increase the dynamic response coefficient, β″. This means that the sensitivity to instantaneous curvature is increased, whereas the sensitivity to dynamic curvature is decreased. In other words, the phase of the curvature response is altered.
Conversely, Ca 2+ leads to ā 1 > 0, and so β′ > 0. The oda1 mutation has a flat profile, with ā 2 ≈ 0 and ā 1 ≈ 0, and so β′ ≈ 0 and β″ ≈ −12. This counters the general belief that outer-arm dyneins affect only beat frequency 47 . The ida3 mutation has a high parabolicity, that is, a larger dip in the middle, which reflects a less negative value of the dynamic curvature coefficient. The mbo2 mutation has a waveform very similar to that of wild type, which is remarkable given the absence of a static mode in the waveform of this mutant 37,41 but is consistent with the observation that the mbo2 mutations do not affect dyneins 47 . Thus, changes in the instantaneous and dynamic curvature sensitivity of the dyneins can account for the effects of environmental and genetic perturbations, although a causal connection between the perturbations and dynein activity has not been proven. Studying the relationship between shape and model parameters for alternative motor mechanisms (for example, 'geometric clutch' 28 or compression instability 29,48 ) may give further insight in this direction.
Elasticity dominates Chlamydomonas swimming. Our results contain three arguments that viscous forces are smaller than elastic forces for Chlamydomonas axonemes. First, the Machin number, which is the ratio of viscous to elastic forces and is related to the Weissenberg number 49 , is much smaller than unity (Ma ≪ 1). For a plane wave with wavenumber q = 2π/λ, we have which ranges from 0.02 to 0.14 under all conditions. Note, however, that there is considerable uncertainty in κ, especially in the presence of taxol 50 (even if the flexural rigidity of the axoneme were reduced 10-fold by taxol, our models shows that there would be only a small effect on the shape). Second, the data collapse of the fitted k and β′′ in Fig. 5c almost coincides with that predicted by the low-viscosity relation (dashed line). This extends our earlier finding that a low Ma accords with wild-type and mbo2 axonemes 32 . Third, despite the wide variability in frequency and mean amplitude, the total range is ~10-fold for both, there is an absence of high-amplitude and high-frequency beats. The dissipation of elastic energy is proportional to a 2 0 f whereas the dissipation of viscous energy is proportional to a 2 0 f 2 (ref. 26 ). Thus, if elastic or viscous energy were limited by the energy available from the hydrolysis of ATP by the motors, then a 0 ~ f −1/2 or a 0 ~ f −1 , respectively. The former provides a better bound on the data than the latter (Fig. 3b dashed and dotted curves, respectively), as expected when elastic dissipation dominates.
A puzzling finding is that the Machin number measured experimentally is relatively insensitive to changing viscosity because increasing viscosity decreases the beat frequency, compensating for much of the expected change in Ma. This highlights our lack of understanding of how the beat frequency is selected. The low-viscosity regime is also supported by the literature. In ref. 33 , the ATPase rate of axonemes was shown to increase in proportion to beat frequency, as predicted for elastic dissipation, and in ref. 51 measurements of the flow field around Chlamydomonas cilia showed that friction forces are smaller than the estimated bending forces. Furthermore, a low Ma may also explain the different relationship between curvature and sliding force reported in ref. 52 relative to that of our model 32 , as the larger Ma of sperm may allow for a different motor control mechanism 31 . We conclude that Chlamydomonas swims in the regime of low Ma, that is, low friction. Because the Reynolds number Re, which is the ratio of inertial to viscous forces, is also small, we have Re ≪ Ma ≪ 1.
Perspective. The functional robustness of the ciliary beat, evidenced by its persistence in the face of genetic and environmental perturbations, indicates that the ciliary beat is a highly canalized process 53 . Such processes can tolerate genetic variations, pre-adaptations, which can be selected for if large environmental changes occur [54][55][56] . This may underlie the diversity of axonemal forms. For example, the number of doublet microtubules in motile axonemes ranges from as few as three 57 to as many as hundreds 58 . Furthermore, radial spokes and the central pair are missing in some motile cilia, and additional asymmetrically localized molecules lead to planar or asymmetric beat patterns 59 . Axonemes can be long or short, can be solitary or arrayed and can drive cell motility or fluid flow. This diversity of axonemal structure and function likely originates, at least in part, in the strong tendency of motor-driven bending of cytoskeletal filaments to produce oscillations when the bending feeds back on the activity of the motors, as we postulate for dynein. The propensity for oscillation likely underlies the successful reconstitution of cilia-like motility using heterologous motors and filaments 33,34 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41567-021-01446-2. for the phase profile. Fitting to a mechanical model, we map these data to a three-dimensional mechanical space, with axes the two response coefficients of the motors and the shear stiffness. b, Histogram of the R 2 of the fits. More than 90% of the fits have R 2 > 0.9. c, Variability in shear stiffness k and motor response β′′. The values of k and β′′ are tightly correlated, which corresponds to the low-viscosity limit (dashed line). The dotted line corresponds to −4π and the crossed colored lines correspond to the means of each dataset. The inset shows the correlation between curvature response coefficients and phase asymmetry and parabolicity. d, Variability in curvature response coefficients captures waveform variability, with low negative β′ corresponding to low negative ā 1 and thus decaying waveforms. The dotted line corresponds to −4π. The inset shows the correlation between the curvature response coefficients and the asymmetry and parabolicity of the phase.  Table of correlation coefficients of amplitude polynomial coefficients with amplitude principal components and phase polynomial coefficients. Note that ā 1 and ā 2 are strongly correlated with the first two principal components, a pc 1 and a pc 2 , respectively, and are anti-correlated and correlated with φ 2 and φ 1 , respectively. The green backgrounds are scatter plots of data.

Experiments. Preparation and reactivation of axonemes. Axonemes from
Chlamydomonas reinhardtii cells (received from http://chlamy.org) were purified and reactivated. The procedures described here are detailed in ref. 60 .
Chemicals were purchased from Sigma Aldrich unless stated otherwise. In brief, cells were grown in tris-acetate-phosphate medium with additional phosphate (TAP+P medium) under conditions of continuous illumination (from two 75 W fluorescent bulbs) and of air bubbling at 24 °C over the course of 2 days, to a final density of 10 6 cells per millilitre. Cilia were isolated using dibucaine, then purified on a 25% sucrose cushion and demembranated in HMDEK buffer (30 mM HEPES-KOH, 5 mM MgSO 4 , 1 mM dithiothreitol, 1 mM EGTA and 50 mM potassium acetate, at pH 7.4) augmented with 1% (v/v) IGEPAL detergent and 0.2 mM Pefabloc SC protease inhibitor. The membrane-free axonemes were resuspended in HMDEK plus 1% (w/v) polyethylene glycol (molecular weight 20 kDa), 30% sucrose and 0.2 mM Pefabloc, and were stored at −80 °C. Prior to reactivation, axonemes were thawed at room temperature, then kept on ice. Thawed axonemes were used for up to 2 h.
Reactivation was performed in flow chambers with a depth of 100 μm built from easy-clean glass and double-sided sticky tape. Thawed axonemes were diluted in HMDEKP (HMDEK augmented with 1% (w/v) PEG 20.000) reactivation buffer. Unless stated otherwise, a standard reactivation buffer containing 1 mM ATP and an ATP-regeneration system (5 units per millilitre creatine kinase and 6 mM creatine phosphate) was used to maintain a constant ATP concentration. The axoneme dilution was infused into a glass chamber, which was blocked using casein solution (from bovine milk, 2 mg mL −1 ) for 10 min and then sealed with vacuum grease. Prior to imaging, the sample was equilibrated on the microscope for 5 min and data were collected for a maximum time of 20 min.
Special reactivation conditions. For the temperature series, the temperature was controlled using an objective heater from Bioptech. Unless stated otherwise, the sample temperature was kept constant at 24 °C. The temperature series was acquired by increasing the temperature in 2 °C steps and letting the system equilibrate for 10 min. After equilibration, the target temperature was checked using an inbuilt reference thermistor. For the ATP series, the standard buffer (without ATP) was augmented with different amounts of ATP (50, 66, 100, 240, 370, 500, 750 and 1,000 μM). For the viscosity series, the standard buffer was augmented with Ficoll 400 (1%, 5% 10% (w/v)), and then axonemes were added to this solution. For the calcium, we used a Ca 2+ buffered reactivation solution with a concentration of free Ca 2+ of 100 μM (calculated with the program Maxchelator https://somapp.ucdmc. ucdavis.edu/pharmacology/bers/maxchelator). This concentration was chosen as it converts the asymmetric beat into a symmetric one 61 . For the taxol, the standard buffer was augmented with 10 μM taxol, a concentration that stabilizes polymerized microtubules 62 , and then axonemes were added to this solution.
Imaging of axonemes. The reactivated axonemes were imaged by phase-contrast microscopy, set up on an inverted Zeiss Axiovert S100 TV or Zeiss Observer Z1 microscope using a Zeiss 63× Plan-Apochromat NA 1.4 or a 40× Plan-Neofluar NA 1.3 Phase 3 oil lens in combination with a 1.6× tube lens and a Zeiss oil condenser (NA 1.4). Movies were acquired using an EoSens 3CL CMOS high-speed camera. The effective pixel size was 139 nm or 219 nm per pixel. Movies of up to 3,000 frames were recorded at a frame rate of 1,000 fps.

Data analysis.
High-precision tracking of isolated axonemes. To track the shape of the axoneme in each movie frame with nanometre precision, the MATLAB-based software tool FIESTA Version 1.03 was used 35 . Prior to tracking, movies were background-subtracted to remove static inhomogeneities arising from uneven illumination and from dirt particles. The background image contained the mean intensity in each pixel calculated over the entire movie. This procedure increased the signal-to-noise ratio by a factor of three 60 . Phase-contrast images were inverted. For tracking, a segment size of 733 nm (approximately 5 × 5 pixels) was used, corresponding to program settings of a full width at half maximum of 750 nm and of a 'reduced box size for tracking especially curved filaments' of 30%. Along the arc length of each filament, 25 equally spaced segments were fitted using two-dimensional Gaussian functions.
Fourier analysis of tangent angle. Using the x, y positions of the filament shape we calculated the tangent angle ψ(s, t) at every arc-length position in time. To study only the dynamic modes of the waveform we subtracted the static mode (the 0th Fourier mode) 37 and Fourier-decomposed the tangent angle. The spectrum of dynamic modes of the tangent angle shows a peak at the fundamental frequency (Fig. 1d, first Fourier mode). Because the fundamental mode accounts for >90% of the total power, we neglected the higher harmonics (n = 2, 3, 4, …) in all further analysis (see ref. 32 for more detail) and considered only the first Fourier mode of ψ(s, t), which is denoted as ψ(s) and will be used hereafter for all dynamic quantities. We calculated the frequency f = ⟨f(s)⟩ and the arc-length-dependent amplitude a(s) = |ψ(s)| and phase φ(s) = arg ψ(s) profiles of the fundamental Fourier mode.
Polynomial decomposition. We use shifted Legendre polynomials, ℓ n (s), with normalization ∫ 1 0 ℓn(s)ℓm(s)ds = 1 2n+1 δmn, to characterize amplitude and phase profiles. Any function g(s) with support [0, 1] can be written as g(s) = ∑ n≥0 g n ℓ n (s), with gn = (2n + 1) ∫ 1 0 g(s)ℓn(s)ds. With this definition, g0 = ∫ 1 0 g(s)ds is the mean over the arc length. Theory. Mechanical model of the axoneme. The dynamics of the axoneme is characterized by a balance of hydrodynamic, elastic and internal sliding forces 26,63 . The sliding forces can be active, exerted by motors, or passive, associated with structures such as nexin links and radial spokes in addition to the motors. To linear order and in frequency space, the force balance equation, which we call the Machin equation, is given by where f sl is the dimensionless sliding force density and Ma = 2πfξ n L 4 /κ is a dimensionless constant related to the Machin number by Ma = (2π) 4 Ma. A global balance of sliding forces also holds, where the basal force is given by We take the sliding force density f sl (s) generated by motors to depend on the local curvature ∂ s ψ(s) and sliding ∆(s) = ∆ b + ψ(s) − ψ(0) of doublets 32 . In particular, where k > 0 is a passive sliding stiffness and β = β′ + iβ′′ is a complex curvature response coefficient. For β′′ < 0 there is forward (base-to-tip) wave propagation 64 . Note that, unlike in ref. 64 , we allow for instantaneous curvature response via β′, which is key for characterizing amplitude asymmetry (in ref. 32 , β′ = 0 was sufficient because only the symmetric wild-type and mbo2 axonemes were studied). The model does not speak to the mechanism of curvature sensing, which must be very sensitive given that axonemal curvatures are low. Such curvature sensitivity might be realized by strains that distort the cross-section of the axoneme, for example due to the 'geometric clutch' model 28 or to the coupling between twist and bending 64 .
As boundary conditions we use torque and force balances: (4) and (5) correspond to the solution of a non-linear problem near the point of instability 63 . Note that we have not considered the role of three-dimensional components to the beat, for example via twist 64 .
Model fitting. In previous works, we have used the Machin equation and its variants to reproduce the observed beating patterns of cilia 31,64 . For a given value of Ma, estimated from experiments, and values for the sliding response coefficients k, β′ and β′′, there exists a discrete set of solutions to the system of equations posed by equations (4) and (5), the global sliding balance and the boundary conditions 36,63 . We keep the solution with the lowest wavenumber, because for simple dynamical non-linear motor models it has been shown to be the first one to be excited (see Section 5.2.2 in ref. 65 ), so that for a given set of motor response parameters we can determine a unique theoretical waveform, {k, β ′ , β ′′ } → ψ the (s). The procedure to fit the model to the data has been extensively described in ref. 36 . In brief, we define a score between the theoretically predicted waveform and the experimentally observed one, R 2 (ψ the , ψ exp ), so that R 2 = 1 corresponds to a perfect matching of theory and experiments 36,64 . We then maximize this score with respect to the three parameters {k, β′, β′′} using a principal axis algorithm with three different initial points: {3π 2 , 0, −4π}, which corresponds to the plane wave approximation, and two perturbations around it, {1.7, −2.1, −5.8} and {3π 2 , 0.15, −9.6}. In addition, to calculate Ma, we used κ = 580 pN μm 2 and ξ n = 0.0049 pN s μm −2 (ref. 64 ). There is a considerable range of values of these parameters in the literature. For example, in ref. 66 the normal friction coefficient was estimated to be ξ n = 0.0015 pN s μm −2 based on measurements in intact cells, and in ref. 46 a bending modulus of κ = 800 pN μm 2 was obtained using optical tweezers. Because we are in the low-Machin-number regime, our analysis is quite insensitive to changes in these two parameters, which was confirmed by increasing and decreasing the Machin number 10-fold. In the Ficoll datasets we adjusted ξ n following our own measurements by factors of 1.1, 1.6 and 3.2 for dilutions of 1%, 5% and 10%, respectively.
Low friction limit. Although equations (4) and (5) are linear, the boundary conditions result in complicated relationships between the exponential scales in the solution and the motor parameters 36 . We found that such complications are greatly reduced in the limit Ma → 0. It is particularly convenient to work in terms of curvature ∂ s ψ, which can be shown to take the form ∂sψ = e cs (e cd(s−1) − e −cd(s−1) ) , with c = β/2 and d = √ 1 + 4k/β 2 where the sliding stiffness, k, is now allowed to be complex, with the imaginary component corresponding to sliding friction. The general condition for the existence of a solution in this limit is simply As the observed values for the basal stiffness arising from the fits are small, we explore the limit χ b → 0. In this limit equation (7) becomes tanh(cd) = 0, which implies cd = inπ, with n = 1, 2, … From this we obtain with the sliding stiffness and curvature response coefficients related by β √ 1 + 4k/β 2 = i2πn. Note that, in the case that k is real, we have β′ = 0. The dashed line in Fig. 5 corresponds to the previous equation when n = 1, which provides good agreement with the fits for the full model (finite viscosity, non-zero basal compliance and non-zero asymmetry). Note that in equation (8) the amplitude asymmetry is directly related to β′, with β′ < 0 corresponding to decreasing amplitude. Furthermore, the wavelength can be directly read off from this expression, and is given by λ = −4πL/β′′. This is also in agreement with observations from fitting the full model.