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 spatiotemporal 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 lowdimensional 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 lowviscosity 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 lowdimensional. 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 lowdimensional 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 coswimming 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}\equiv \int\nolimits_{0}^{1}a(s){\mathrm{d}}s\) and the wavelength \(\lambda \equiv \frac{\uppi L}{6}\int\nolimits_{0}^{1}(s1/2)\varphi (s){\mathrm{d}}s\), here we are interested in their arclength 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)/a_{0}, φ(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, 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 betweenconditions variance is about seven times larger than the withincondition variance. By contrast, the mean amplitude varies little between these conditions, although it shows large variations within conditions: the betweenconditions variance in amplitude is about onefifth the withinconditions 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 outerarm 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 axonemetoaxoneme 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 axonemetoaxoneme variance (Supplementary Fig. 1 and Supplementary Table 1). It is nonetheless possible that variability within a given condition is due to longterm 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 Qfactor of the oscillations can reach values as high as 150 (Fig. 3b, inset), confirming earlier observations^{44}. This Qfactor, 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 shortterm 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: \({\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 Ca^{2+}, 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 anticorrelated 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 lowdimensionality 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 nondimensional 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 nondimensionalizing 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 curvefitting 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 \({\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 lowfriction 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 \(({\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. Lowdimensional phenotypic spaces have been observed in other behaviours (Introduction), but the range of conditions probed in our study makes the lowdimensionality 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 (nondimensionalized 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}:
where ‘~’ denotes correlation. Thus, the different shapes correspond to different curvature sensitivities.
Wildtype 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, Ca^{2+} 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 outerarm 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 10fold 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 lowviscosity relation (dashed line). This extends our earlier finding that a low Ma accords with wildtype and mbo2 axonemes^{32}. Third, despite the wide variability in frequency and mean amplitude, the total range is ~10fold for both, there is an absence of highamplitude and highfrequency 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 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 lowviscosity 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, preadaptations, 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 motordriven 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 cilialike motility using heterologous motors and filaments^{33,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 10^{6} cells per millilitre. Cilia were isolated using dibucaine, then purified on a 25% sucrose cushion and demembranated in HMDEK buffer (30 mM HEPESKOH, 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 membranefree 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 easyclean glass and doublesided 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 ATPregeneration 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 phasecontrast microscopy, set up on an inverted Zeiss Axiovert S100 TV or Zeiss Observer Z1 microscope using a Zeiss 63× PlanApochromat NA 1.4 or a 40× PlanNeofluar 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 highspeed 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
Highprecision tracking of isolated axonemes
To track the shape of the axoneme in each movie frame with nanometre precision, the MATLABbased software tool FIESTA Version 1.03 was used^{35}. Prior to tracking, movies were backgroundsubtracted 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 signaltonoise ratio by a factor of three^{60}. Phasecontrast 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 twodimensional 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 arclength position in time. To study only the dynamic modes of the waveform we subtracted the static mode (the 0th Fourier mode)^{37} and Fourierdecomposed 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 arclengthdependent 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≥0} g_{n}ℓ_{n}(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 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 \(\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 f_{sl}(s) generated by motors to depend on the local curvature ∂_{s}ψ(s) and sliding \(\varDelta(s)\) = \(\varDelta_{\mathrm{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 (basetotip) 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 wildtype 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 crosssection 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: ∂_{s}ψ(0) = F_{b}, ∂_{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 nonlinear problem near the point of instability^{63}. Note that we have not considered the role of threedimensional 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 \(\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 conditions^{36,63}. We keep the solution with the lowest wavenumber, because for simple dynamical nonlinear 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, 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 \(\overline{\mathrm{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 lowMachinnumber regime, our analysis is quite insensitive to changes in these two parameters, which was confirmed by increasing and decreasing the Machin number 10fold. 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 \(\overline{\mathrm{Ma}}\to 0\). It is particularly convenient to work in terms of curvature ∂_{s}ψ, which can be shown to take the form
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
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 \(\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, nonzero basal compliance and nonzero 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.
Data availability
The data shown in the main figures are appended to this article in csv format. In addition, all tracked waveforms are available on Dryad (https://doi.org/10.5061/dryad.0gb5mkm2j).
Change history
23 February 2022
A Correction to this paper has been published: https://doi.org/10.1038/s41567022015529
References
Stephens, G. J., JohnsonKerner, B., Bialek, W. & Ryu, W. S. Dimensionality and dynamics in the behavior of C. elegans. PLoS Comput Biol 4, e1000028 (2008).
Helms, S. J. et al. Modelling the ballistictodiffusive transition in nematode motility reveals variation in exploratory behaviour across species. J. R. Soc. Interface 16, 20190174 (2019).
Jordan, D., Kuehn, S., Katifori, E. & Leibler, S. Behavioral diversity in microbes and lowdimensional phenotypic spaces. Proc. Natl Acad. Sci. USA 110, 14018–14023 (2013).
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).
Berman, G. J., Bialek, W. & Shaevitz, J. W. Predictability and hierarchy in Drosophila behavior. Proc. Natl Acad. Sci. USA 113, 11943–11948 (2016).
Mathis, A. et al. DeepLabCut: markerless pose estimation of userdefined body parts with deep learning. Nat. Neurosci. 21, 1281–1289 (2018).
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).
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).
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).
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).
O'Callaghan, C., Achaval, M., Forsythe, I. & Barry, P. W. Brain and respiratory cilia: the effect of temperature. Biol. Neonate 68, 394–397 (1995).
Rikmenspoel, R. Movements and active moments of bull sperm flagella as a function of temperature and viscosity. J. Exp. Biol. 108, 205–230 (1984).
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).
Brokaw, C. J. Effects of viscosity and ATP concentration on the movement of reactivated seaurchin sperm flagella. J. Exp. Biol. 62, 701–719 (1975).
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).
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).
Ishijima, S. & Witman, G. B. Flagellar movement of intact and demembranated, reactivated ram spermatozoa. Cell Motil. Cytoskeleton 8, 375–391 (1987).
Wan, K. Y. & Jékely, G. On the unity and diversity of cilia. Philos. Trans. R. Soc. B 375, 20190148 (2020).
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).
Summers, K. E. & Gibbons, I. R. Adenosine triphosphateinduced sliding of tubules in trypsintreated flagella of seaurchin sperm. Proc. Natl Acad. Sci. USA 68, 3092–3096 (1971).
Nicastro, D. et al. The molecular architecture of axonemes revealed by cryoelectron tomography. Science 313, 944–948 (2006).
Warner, F. D. Ciliary intermicrotubule bridges. J. Cell Sci. 20, 101–114 (1976).
Brokaw, C. J. Thinking about flagellar oscillation. Cell Motil. Cytoskeleton 66, 425–436 (2009).
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).
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).
Machin, K. E. Wave propagation along flagella. J. Exp. Biol. 35, 796–806 (1958).
Brokaw, C. J. Computer simulation of flagellar movement. VI. Simple curvaturecontrolled models are incompletely specified. Biophys. J. 48, 633–642 (1985).
Lindemann, C. B. A ‘geometric clutch’ hypothesis to explain oscillations of the axoneme of cilia and flagella. J. Theor. Biol. 168, 175–189 (1994).
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).
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).
RiedelKruse, I. H., Hilfinger, A., Howard, J. & Jülicher, F. How molecular motors shape the flagellar beat. HFSP J. 1, 192–208 (2007).
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).
Sanchez, T., Chen, D. T., DeCamp, S. J., Heymann, M. & Dogic, Z. Spontaneous motion in hierarchically assembled active matter. Nature 491, 431–434 (2012).
Pochitaloff, M. et al. Self organized wave like beating of actin bundles in a minimal actomyosin system of controlled architecture. Biophys. J. 114, 649a (2018).
Ruhnow, F., Zwicker, D. & Diez, S. Tracking single particles and elongated filaments with nanometer precision. Biophys. J. 100, 2820–2828 (2011).
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).
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).
Mitchell, D. R. & Rosenbaum, J. L. A motile Chlamydomonas flagellar mutant that lacks outer dynein arms. J. Cell Biol. 100, 1228–1234 (1985).
Kamiya, R., Kurimoto, E. & Muto, E. Two types of Chlamydomonas flagellar mutants missing different components of innerarm dynein. J. Cell Biol. 112, 441–447 (1991).
KatoMinoura, T., Hirono, M. & Kamiya, R. Chlamydomonas innerarm dynein mutant, ida5, has a mutation in an actinencoding gene. J. Cell Biol. 137, 649–656 (1997).
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).
Kubo, T., Yanagisawa, H., Yagi, T., Hirono, M. & Kamiya, R. Tubulin polyglutamylation regulates axonemal motility by modulating activities of innerarm dyneins. Curr. Biol. 20, 441–445 (2010).
Alford, L. M. et al. The nexin link and Btubule glutamylation maintain the alignment of outer doublets in the ciliary axoneme. Cytoskeleton 73, 331–340 (2016).
Foster, K. W., Vidyadharan, J. & Sangani, A. S. Evidence for a selforganized compliant mechanism for the spontaneous steady beating of cilia. Cytoskeleton 74, 260–280 (2017).
Kössl, M. & Russell, I. Basilar membrane resonance in the cochlea of the mustached bat. Proc. Natl Acad. Sci. USA 92, 276–279 (1995).
Xu, G. et al. Flexural rigidity and shear stiffness of flagella estimated from induced bends and counterbends. Biophys. J. 110, 2759–2768 (2016).
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).
Ling, F., Guo, H. & Kanso, E. Instabilitydriven oscillations of elastic microfilaments. J. R. Soc. Interface 15, 20180594 (2018).
Poole, R. J. The Deborah and Weissenberg numbers. Rheol. Bull. 53, 32–39 (2012).
Schaedel, L. et al. Microtubules selfrepair in response to mechanical stress. Nat. Mater. 14, 1156–1163 (2015).
Mondal, D., Adhikari, R. & Sharma, P. Internal friction controls active ciliary oscillations near the instability threshold. Sci. Adv. 6, eabb0503 (2020).
Lin, J. & Nicastro, D. Asymmetric distribution and spatial switching of dynein activity generates ciliary motility. Science 360, eaar1968 (2018).
Waddington, C. H. Canalization of development and the inheritance of acquired characters. Nature 150, 563–565 (1942).
Rutherford, S. L. & Lindquist, S. Hsp90 as a capacitor for morphological evolution. Nature 396, 336–342 (1998).
Eshel, I. & Matessi, C. Canalization, genetic assimilation and preadaptation: a quantitative genetic model. Genetics 149, 2119–2133 (1998).
Flatt, T. The evolutionary genetics of canalization. Q. Rev. Biol. 80, 287–316 (2005).
Prensier, G., Vivier, E., Goldstein, S. & Schrével, J. Motile flagellum with a ‘3 + 0’ ultrastructure. Science 207, 1493–1494 (1980).
Mooseker, M. S. & Tilney, L. G. Isolation and reactivation of the axostyle: evidence for a dyneinlike ATPase in the axostyle. J. Cell Biol. 56, 13–26 (1973).
Dutcher, S. K. Asymmetries in the cilia of Chlamydomonas. Philos. Trans. R. Soc. B 375, 20190153 (2020).
Alper, J., Geyer, V., Mukundan, V. & Howard, J. in Methods in Enzymology Vol. 524 (ed. Marshall, W. F) 343–369 (Elsevier, 2013).
Wakabayashi, K., Yagi, T. & Kamiya, R. Ca^{2+}dependent waveform conversion in the flagellar axoneme of Chlamydomonas mutants lacking the centralpair/radial spoke system. Cell Motil. Cytoskeleton 38, 22–28 (1997).
Schiff, P. B., Fant, J. & Horwitz, S. B. Promotion of microtubule assembly in vitro by taxol. Nature 277, 665–667 (1979).
Camalet, S. & Jülicher, F. Generic aspects of axonemal beating. New J. Phys. 2, 24 (2000).
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).
Sartori, P. Effect of curvature and normal forces on motor regulation of cilia. Preprint at https://arxiv.org/abs/1905.04138 (2019).
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 startup funds used to support this project.
Author information
Affiliations
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
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.
Additional information
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
Postprocessed 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
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Geyer, V.F., Howard, J. & Sartori, P. Ciliary beating patterns map onto a lowdimensional behavioural space. Nat. Phys. 18, 332–337 (2022). https://doi.org/10.1038/s41567021014462
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41567021014462
Further reading

The beat of isolated cilia
Nature Physics (2022)

Light moves artificial cilia to a complex beat
Nature (2022)