Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Ciliary beating patterns map onto a low-dimensional behavioural space

### Subjects

A Publisher Correction to this article was published on 23 February 2022

## Abstract

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.

## Main

Biological 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 nematodes1,2, the swimming trajectories of ciliates and bacteria3,4, the postures of walking flies5 and mice6 and the activities of Hydra7, 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 mechanisms8.

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 proteins9,10,11 and can operate over a wide range of temperatures12, viscosities12,13, ATP concentrations14 and buffer conditions such as pH and calcium concentration15,16,17. Thus, cilia exemplify conserved function18, 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 classes19. These motors drive the sliding of adjacent doublet microtubules20 within the axoneme, the structural core of cilia21. 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 base25. This biomechanical and structural knowledge forms the basis of mechanochemical models that can successfully recapitulate the oscillatory motion of cilia26,27,28,29,30,31,32 and synthetic filament bundles33,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 software35 (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 frame36 (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:

$$\psi (s,t)=a(s)\cos (2\uppi ft+\varphi (s))\,,$$
(1)

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}\equiv \int\nolimits_{0}^{1}a(s){\mathrm{d}}s$$ and the wavelength $$\lambda \equiv -\frac{\uppi L}{6}\int\nolimits_{0}^{1}(s-1/2)\varphi (s){\mathrm{d}}s$$, 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 $${{{\mathcal{P}}}}=\{a(s),\varphi (s),f,L\}$$, which we define as the behavioural space of beating. A subset of these parameters, {a(s)/a0, φ(s)}, defines the shape.

We explored the variability of $${{{\mathcal{P}}}}$$ by observing groups of axonemes under several environmental perturbations (temperature, viscosity, different ATP concentrations, Ca2+ 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 dyneins38), ida3 and ida5 (which lack inner dynein arms f and a, c, d, e, respectively39,40), mbo2 (which lacks microtubule inner proteins and has a symmetric beat41) and tpg1 (which has reduced polyglutamamylation required for axonemal integrity42,43), as well as addition of the ion Ca2+ 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-to-axoneme 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 observations44. 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 (a0) equals the mean amplitude, the first order (a1) corresponds to the linear deviation from constant amplitude and is a measure of the proximal–distal asymmetry, and the second order (a2) 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: $${\bar{a}}_{1}={a}_{1}/{a}_{0}$$ and $${\bar{a}}_{2}={a}_{2}/{a}_{0}$$.

Figure 4d displays the space of variations of parabolicity and asymmetry of the amplitude. Most points cluster around $${\bar{a}}_{1}\approx 0.0$$ and $${\bar{a}}_{2}\approx 0.3$$, which corresponds to a symmetric and convex amplitude profile. The exceptions include taxol, ida5 and tpg1, which have $${\bar{a}}_{1} < 0$$ corresponding to decreasing amplitudes towards the distal tip, and Ca2+, which has $${\bar{a}}_{1} > 0$$ corresponding to increasing amplitude towards the distal tip. The oda1 waveforms have a small parabolicity $${\bar{a}}_{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}_{1}^{{{{\rm{pc}}}}}$$ and $${a}_{2}^{{{{\rm{pc}}}}}$$ (Fig. 4f, left boxes). Furthermore, the asymmetry of the amplitude profile, $${\bar{a}}_{1}$$, is anti-correlated with the parabolicity of the phase profile, φ2, whereas the parabolicity of the amplitude profile, $${\bar{a}}_{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.

### Motor response controls waveform variability

The low-dimensionality of shapes suggests that the perturbations, irrespective of their molecular mechanisms, affect only a few collective 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 $$\overline{\mathrm{Ma}}$$, which is a constant that comes from non-dimensionalizing Machin’s equation (see equation (4) in ref. 26). $$\overline{\mathrm{Ma}}=2\uppi f{\xi }_{\mathrm{n}}{L}^{4}/\kappa$$, where κ is the flexural rigidity of the axoneme and ξn is the normal friction coefficient. For Chlamydomonas, $$\overline{\mathrm{Ma}} \approx 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(si) and φ(si) 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 R2 > 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 $${\mathrm{Ma}}\equiv \overline{\mathrm{Ma}}/{(2\uppi )}^{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 Ca2+ data lie far apart from each other in the (β′, β″) space, just as in the $$({\bar{a}}_{1},{\bar{a}}_{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 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 ($${\bar{a}}_{1}$$ and φ1) and parabolicity ($${\bar{a}}_{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 model32. 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, $${\bar{a}}_{1}$$, $${\bar{a}}_{2}$$, φ1 and φ2:

$${\bar{a}}_{1} \sim \beta ^{\prime} \sim -{\varphi }_{2}\quad {\mathrm{and}}\quad {\bar{a}}_{2} \sim \beta ^{\prime\prime} \sim +{\varphi }_{1}\propto \lambda /L\,,$$
(2)

where ‘~’ denotes correlation. Thus, the different shapes correspond to different curvature sensitivities.

Wild-type cilia have symmetric and convex amplitude profiles, $${\bar{a}}_{1}\approx 0.0$$ and $${\bar{a}}_{2}\approx 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 $${\bar{a}}_{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, Ca2+ leads to $${\bar{a}}_{1} > 0$$, and so β′ > 0. The oda1 mutation has a flat profile, with $${\bar{a}}_{2}\approx 0$$ and $${\bar{a}}_{1}\approx 0$$, and so β′ ≈ 0 and β″ ≈ −12. This counters the general belief that outer-arm dyneins affect only beat frequency47. 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 mutant37,41 but is consistent with the observation that the mbo2 mutations do not affect dyneins47. 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 instability29,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 number49, is much smaller than unity (Ma  1). For a plane wave with wavenumber q = 2π/λ, we have

$${\mathrm{Ma}}=\frac{{\xi }_{\mathrm{n}}\omega }{\kappa {q}^{4}}=\frac{1}{{(2\uppi )}^{4}}\overline{\mathrm{Ma}}\,,$$
(3)

which ranges from 0.02 to 0.14 under all conditions. Note, however, that there is considerable uncertainty in κ, especially in the presence of taxol50 (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 axonemes32. 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}_{0}^{2}f$$ whereas the dissipation of viscous energy is proportional to $${a}_{0}^{2}{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 a0 ~ f−1/2 or a0 ~ 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 model32, as the larger Ma of sperm may allow for a different motor control mechanism31. 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 process53. Such processes can tolerate genetic variations, pre-adaptations, which can be selected for if large environmental changes occur54,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 three57 to as many as hundreds58. Furthermore, radial spokes and the central pair are missing in some motile cilia, and additional asymmetrically localized molecules lead to planar or asymmetric beat patterns59. 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 filaments33,34.

## Methods

### 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 106 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 MgSO4, 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 Ca2+ buffered reactivation solution with a concentration of free Ca2+ 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 one61. For the taxol, the standard buffer was augmented with 10 μM taxol, a concentration that stabilizes polymerized microtubules62, 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 used35. 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 three60. 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 $$\varphi (s)=\arg \psi (s)$$ profiles of the fundamental Fourier mode.

#### Polynomial decomposition

We use shifted Legendre polynomials, n(s), with normalization $$\int\nolimits_{0}^{1}{\ell }_{n}(s){\ell }_{m}(s){\mathrm{d}}s=\frac{1}{2n+1}{\delta }_{mn}$$, to characterize amplitude and phase profiles. Any function g(s) with support [0, 1] can be written as g(s) = ∑n≥0gnn(s), with $${g}_{n}=(2n+1)\int\nolimits_{0}^{1}g(s){\ell }_{n}(s){\mathrm{d}}s$$. With this definition, $${g}_{0}=\int\nolimits_{0}^{1}g(s){\mathrm{d}}s$$ 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 forces26,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

$$i\overline{\mathrm{Ma}}\psi ={\partial }_{s}^{4}\psi +{\partial }_{s}^{2}{f}_{\mathrm{sl}}\,,$$
(4)

where fsl is the dimensionless sliding force density and $$\overline{\mathrm{Ma}}=2\uppi f{\xi }_{\mathrm{n}}{L}^{4}/\kappa$$ is a dimensionless constant related to the Machin number by $$\overline{\mathrm{Ma}} = (2 \pi)^4{\rm{Ma}}$$. A global balance of sliding forces also holds, $$\int\nolimits_{0}^{1}{f}_{\mathrm{sl}}(s){\mathrm{d}}s={F}_{\mathrm{b}}$$, where the basal force is given by $$F_{\mathrm{b}}=\chi_{\mathrm{b}}{\varDelta_{\mathrm{b}}}$$ with $$\varDelta_{\mathrm{b}}$$ the basal sliding and $${\chi }_{\mathrm{b}}={\chi }_{\mathrm{b}}^{\prime} +i{\chi }_{\mathrm{b}}^{\prime\prime}$$ the complex basal response coefficient with $${\chi }_{\mathrm{b}}^{\prime},{\chi }_{\mathrm{b}}^{\prime\prime} > 0$$.

We take the sliding force density fsl(s) generated by motors to depend on the local curvature ∂sψ(s) and sliding $$\varDelta(s)$$ = $$\varDelta_{\mathrm{b}}$$ + ψ(s) − ψ(0) of doublets32. In particular,

$${f}_{\mathrm{sl}}=k{{\varDelta }}+\beta {\partial }_{s}\psi \,,$$
(5)

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 propagation64. 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’ model28 or to the coupling between twist and bending64. As boundary conditions we use torque and force balances: ∂sψ(0) = Fb, ∂sψ(1) = 0, $${\partial }_{s}^{2}\psi (0)={f}_{\mathrm{sl}}(0)$$, and $${\partial }_{s}^{2}\psi (1)={f}_{\mathrm{sl}}(1)$$. Equations (4) and (5) correspond to the solution of a non-linear problem near the point of instability63. Note that we have not considered the role of three-dimensional components to the beat, for example via twist64.

#### Model fitting

In previous works, we have used the Machin equation and its variants to reproduce the observed beating patterns of cilia31,64. For a given value of $$\overline{\mathrm{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 conditions36,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,\beta ^{\prime} ,\beta ^{\prime\prime} \}\to {\psi }_{\mathrm{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, R2(ψthe, ψexp), so that R2 = 1 corresponds to a perfect matching of theory and experiments36,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 $$\overline{\mathrm{Ma}}$$, we used κ = 580 pN μm2 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 μm2 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 parameters36. We found that such complications are greatly reduced in the limit $$\overline{\mathrm{Ma}}\to 0$$. It is particularly convenient to work in terms of curvature ∂sψ, which can be shown to take the form

$${\partial }_{s}\psi ={\mathrm{e}}^{cs}({\mathrm{e}}^{cd(s-1)}-{\mathrm{e}}^{-cd(s-1)})\,,$$
(6)

with c = β/2 and $$d=\sqrt{1+4k/{\beta }^{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

$$\tanh (cd)=-\frac{{\chi }_{\mathrm{b}}cd}{{\chi }_{\mathrm{b}}c+k }\,,$$
(7)

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

$${\partial }_{s}\psi =\exp \left(\beta ^{\prime} s/2\right)\sin \left(n\uppi (s-1)\right)\exp \left(i\beta ^{\prime\prime} s/2\right)\,,$$
(8)

with the sliding stiffness and curvature response coefficients related by $$\beta \sqrt{1+4k/{\beta }^{2}}=i2\uppi 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.

## References

1. Stephens, G. J., Johnson-Kerner, B., Bialek, W. & Ryu, W. S. Dimensionality and dynamics in the behavior of C. elegans. PLoS Comput Biol 4, e1000028 (2008).

2. Helms, S. J. et al. Modelling the ballistic-to-diffusive transition in nematode motility reveals variation in exploratory behaviour across species. J. R. Soc. Interface 16, 20190174 (2019).

3. Jordan, D., Kuehn, S., Katifori, E. & Leibler, S. Behavioral diversity in microbes and low-dimensional phenotypic spaces. Proc. Natl Acad. Sci. USA 110, 14018–14023 (2013).

4. Pleška, M., Jordan, D., Frentz, Z., Xue, B. & Leibler, S. Nongenetic individuality, changeability, and inheritance in bacterial behavior. Proc. Natl Acad. Sci. USA 118, e2023322118 (2021).

5. Berman, G. J., Bialek, W. & Shaevitz, J. W. Predictability and hierarchy in Drosophila behavior. Proc. Natl Acad. Sci. USA 113, 11943–11948 (2016).

6. Mathis, A. et al. DeepLabCut: markerless pose estimation of user-defined body parts with deep learning. Nat. Neurosci. 21, 1281–1289 (2018).

7. Han, S., Taralova, E., Dupre, C. & Yuste, R. Comprehensive machine learning analysis of Hydra behavior reveals a stable basal behavioral repertoire. eLife 7, e32605 (2018).

8. Wang, H. et al. From neuron to muscle to movement: a complete biomechanical model of Hydra contractile behaviors. Preprint at bioRxiv https://doi.org/10.1101/2020.12.14.422784 (2020).

9. Brokaw, C. J. & Kamiya, R. Bending patterns of Chlamydomonas flagella: IV. Mutants with defects in inner and outer dynein arms indicate differences in dynein arm function. Cell Motil. Cytoskeleton 8, 68–75 (1987).

10. Brokaw, C. J. & Luck, D. J. Bending patterns of Chlamydomonas flagella: III. A radial spoke head deficient mutant and a central pair deficient mutant. Cell Motil. Cytoskeleton 5, 195–208 (1985).

11. O'Callaghan, C., Achaval, M., Forsythe, I. & Barry, P. W. Brain and respiratory cilia: the effect of temperature. Biol. Neonate 68, 394–397 (1995).

12. Rikmenspoel, R. Movements and active moments of bull sperm flagella as a function of temperature and viscosity. J. Exp. Biol. 108, 205–230 (1984).

13. Yagi, T. et al. An axonemal dynein particularly important for flagellar movement at high viscosity. Implications from a new Chlamydomonas mutant deficient in the dynein heavy chain gene DHC9. J. Biol. Chem. 280, 41412–41420 (2005).

14. Brokaw, C. J. Effects of viscosity and ATP concentration on the movement of reactivated sea-urchin sperm flagella. J. Exp. Biol. 62, 701–719 (1975).

15. Bessen, M., Fay, R. B. & Witman, G. B. Calcium control of waveform in isolated flagellar axonemes of Chlamydomonas. J. Cell Biol. 86, 446–455 (1980).

16. Omoto, C. K. & Brokaw, C. J. Bending patterns of Chlamydomonas flagella: II. Calcium effects on reactivated Chlamydomonas flagella. Cell Motil. Cytoskeleton 5, 53–60 (1985).

17. Ishijima, S. & Witman, G. B. Flagellar movement of intact and demembranated, reactivated ram spermatozoa. Cell Motil. Cytoskeleton 8, 375–391 (1987).

18. Wan, K. Y. & Jékely, G. On the unity and diversity of cilia. Philos. Trans. R. Soc. B 375, 20190148 (2020).

19. Bui, K. H., Yagi, T., Yamamoto, R., Kamiya, R. & Ishikawa, T. Polarity and asymmetry in the arrangement of dynein and related structures in the Chlamydomonas axoneme. J. Cell Biol. 198, 913–925 (2012).

20. Summers, K. E. & Gibbons, I. R. Adenosine triphosphate-induced sliding of tubules in trypsin-treated flagella of sea-urchin sperm. Proc. Natl Acad. Sci. USA 68, 3092–3096 (1971).

21. Nicastro, D. et al. The molecular architecture of axonemes revealed by cryoelectron tomography. Science 313, 944–948 (2006).

22. Warner, F. D. Ciliary inter-microtubule bridges. J. Cell Sci. 20, 101–114 (1976).

23. Brokaw, C. J. Thinking about flagellar oscillation. Cell Motil. Cytoskeleton 66, 425–436 (2009).

24. Heuser, T., Raytchev, M., Krell, J., Porter, M. E. & Nicastro, D. The dynein regulatory complex is the nexin link and a major regulatory node in cilia and flagella. J. Cell Biol. 187, 921–933 (2009).

25. Satir, P. Studies on cilia: II. Examination of the distal region of the ciliary shaft and the role of the filaments in motility. J. Cell Biol. 26, 805–834 (1965).

26. Machin, K. E. Wave propagation along flagella. J. Exp. Biol. 35, 796–806 (1958).

27. Brokaw, C. J. Computer simulation of flagellar movement. VI. Simple curvature-controlled models are incompletely specified. Biophys. J. 48, 633–642 (1985).

28. Lindemann, C. B. A ‘geometric clutch’ hypothesis to explain oscillations of the axoneme of cilia and flagella. J. Theor. Biol. 168, 175–189 (1994).

29. Bayly, P. V. & Dutcher, S. K. Steady dynein forces induce flutter instability and propagating waves in mathematical models of flagella. J. R. Soc. Interface 13, 20160523 (2016).

30. Bayly, P. V. & Wilson, K. S. Equations of interdoublet separation during flagella motion reveal mechanisms of wave propagation and instability. Biophys. J. 107, 1756–1772 (2014).

31. Riedel-Kruse, I. H., Hilfinger, A., Howard, J. & Jülicher, F. How molecular motors shape the flagellar beat. HFSP J. 1, 192–208 (2007).

32. Sartori, P., Geyer, V. F., Scholich, A., Jülicher, F. & Howard, J. Dynamic curvature regulation accounts for the symmetric and asymmetric beats of Chlamydomonas flagella. eLife 5, e13258 (2016).

33. Sanchez, T., Chen, D. T., DeCamp, S. J., Heymann, M. & Dogic, Z. Spontaneous motion in hierarchically assembled active matter. Nature 491, 431–434 (2012).

34. Pochitaloff, M. et al. Self organized wave like beating of actin bundles in a minimal acto-myosin system of controlled architecture. Biophys. J. 114, 649a (2018).

35. Ruhnow, F., Zwicker, D. & Diez, S. Tracking single particles and elongated filaments with nanometer precision. Biophys. J. 100, 2820–2828 (2011).

36. Geyer, V. F., Sartori, P., Jülicher, F. & Howard, J. in Dyneins: Structure, Biology and Disease 2nd edn, Vol. 2 (ed. King, S. M.) 192–212 (Elsevier, 2018).

37. Geyer, V. F., Sartori, P., Friedrich, B. M., Jülicher, F. & Howard, J. Independent control of the static and dynamic components of the Chlamydomonas flagellar beat. Curr. Biol. 26, 1098–1103 (2016).

38. Mitchell, D. R. & Rosenbaum, J. L. A motile Chlamydomonas flagellar mutant that lacks outer dynein arms. J. Cell Biol. 100, 1228–1234 (1985).

39. Kamiya, R., Kurimoto, E. & Muto, E. Two types of Chlamydomonas flagellar mutants missing different components of inner-arm dynein. J. Cell Biol. 112, 441–447 (1991).

40. Kato-Minoura, T., Hirono, M. & Kamiya, R. Chlamydomonas inner-arm dynein mutant, ida5, has a mutation in an actin-encoding gene. J. Cell Biol. 137, 649–656 (1997).

41. Segal, R. A., Huang, B., Ramanis, Z. & Luck, D. J. Mutant strains of Chlamydomonas reinhardtii that move backwards only. J. Cell Biol. 98, 2026–2034 (1984).

42. Kubo, T., Yanagisawa, H., Yagi, T., Hirono, M. & Kamiya, R. Tubulin polyglutamylation regulates axonemal motility by modulating activities of inner-arm dyneins. Curr. Biol. 20, 441–445 (2010).

43. Alford, L. M. et al. The nexin link and B-tubule glutamylation maintain the alignment of outer doublets in the ciliary axoneme. Cytoskeleton 73, 331–340 (2016).

44. Foster, K. W., Vidyadharan, J. & Sangani, A. S. Evidence for a self-organized compliant mechanism for the spontaneous steady beating of cilia. Cytoskeleton 74, 260–280 (2017).

45. Kössl, M. & Russell, I. Basilar membrane resonance in the cochlea of the mustached bat. Proc. Natl Acad. Sci. USA 92, 276–279 (1995).

46. Xu, G. et al. Flexural rigidity and shear stiffness of flagella estimated from induced bends and counterbends. Biophys. J. 110, 2759–2768 (2016).

47. Porter, M. E. & Sale, W. S. The 9 + 2 axoneme anchors multiple inner arm dyneins and a network of kinases and phosphatases that control motility. J. Cell Biol. 151, 37–42 (2000).

48. Ling, F., Guo, H. & Kanso, E. Instability-driven oscillations of elastic microfilaments. J. R. Soc. Interface 15, 20180594 (2018).

49. Poole, R. J. The Deborah and Weissenberg numbers. Rheol. Bull. 53, 32–39 (2012).

50. Schaedel, L. et al. Microtubules self-repair in response to mechanical stress. Nat. Mater. 14, 1156–1163 (2015).

51. Mondal, D., Adhikari, R. & Sharma, P. Internal friction controls active ciliary oscillations near the instability threshold. Sci. Adv. 6, eabb0503 (2020).

52. Lin, J. & Nicastro, D. Asymmetric distribution and spatial switching of dynein activity generates ciliary motility. Science 360, eaar1968 (2018).

53. Waddington, C. H. Canalization of development and the inheritance of acquired characters. Nature 150, 563–565 (1942).

54. Rutherford, S. L. & Lindquist, S. Hsp90 as a capacitor for morphological evolution. Nature 396, 336–342 (1998).

55. Eshel, I. & Matessi, C. Canalization, genetic assimilation and preadaptation: a quantitative genetic model. Genetics 149, 2119–2133 (1998).

56. Flatt, T. The evolutionary genetics of canalization. Q. Rev. Biol. 80, 287–316 (2005).

57. Prensier, G., Vivier, E., Goldstein, S. & Schrével, J. Motile flagellum with a ‘3 + 0’ ultrastructure. Science 207, 1493–1494 (1980).

58. Mooseker, M. S. & Tilney, L. G. Isolation and reactivation of the axostyle: evidence for a dynein-like ATPase in the axostyle. J. Cell Biol. 56, 13–26 (1973).

59. Dutcher, S. K. Asymmetries in the cilia of Chlamydomonas. Philos. Trans. R. Soc. B 375, 20190153 (2020).

60. Alper, J., Geyer, V., Mukundan, V. & Howard, J. in Methods in Enzymology Vol. 524 (ed. Marshall, W. F) 343–369 (Elsevier, 2013).

61. Wakabayashi, K., Yagi, T. & Kamiya, R. Ca2+-dependent waveform conversion in the flagellar axoneme of Chlamydomonas mutants lacking the central-pair/radial spoke system. Cell Motil. Cytoskeleton 38, 22–28 (1997).

62. Schiff, P. B., Fant, J. & Horwitz, S. B. Promotion of microtubule assembly in vitro by taxol. Nature 277, 665–667 (1979).

63. Camalet, S. & Jülicher, F. Generic aspects of axonemal beating. New J. Phys. 2, 24 (2000).

64. Sartori, P., Geyer, V. F., Howard, J. & Jülicher, F. Curvature regulation of the ciliary beat through axonemal twist. Phys. Rev. E 94, 042426 (2016).

65. Sartori, P. Effect of curvature and normal forces on motor regulation of cilia. Preprint at https://arxiv.org/abs/1905.04138 (2019).

66. Bayly, P. et al. Propulsive forces on the flagellum during locomotion of Chlamydomonas reinhardtii. Biophys. J. 100, 2716–2725 (2011).

## Acknowledgements

V.F.G., P.S. and J.H. acknowledge support from the Max Planck Society during the early phases of this work. J.H. would also like to thank Yale University for start-up funds used to support this project.

## Author information

Authors

### Contributions

V.F.G. performed the experiments, P.S. and V.F.G. carried out data analysis, P.S. was responsible for modelling, P.S., V.F.G and J.H. conceptualized the work, and wrote the paper together.

### Corresponding authors

Correspondence to Veikko F. Geyer, Jonathon Howard or Pablo Sartori.

## Ethics declarations

### Competing interests

The authors have no competing interests.

## Peer review

### Peer review information

Nature Physics thanks Philip Bayly and Kirsty Wan for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

### Supplementary Information

Supplementary Fig. 1, Table 1 and description of data file.

### Supplementary Data

Post-processed waveform data. The first row contains descriptors of each column, whereas the subsequent 498 rows correspond to the axonemes analysed in this work. More details on the structure of this file are found in the Supplementary Information.

## Rights and permissions

Reprints and Permissions

Geyer, V.F., Howard, J. & Sartori, P. Ciliary beating patterns map onto a low-dimensional behavioural space. Nat. Phys. 18, 332–337 (2022). https://doi.org/10.1038/s41567-021-01446-2

• Accepted:

• Published:

• Issue Date:

• DOI: https://doi.org/10.1038/s41567-021-01446-2

• ### The beat of isolated cilia

• Kirsty Y. Wan

Nature Physics (2022)

• ### Light moves artificial cilia to a complex beat

• Dhanya Babu
• Nathalie Katsonis

Nature (2022)