Waveform of free, hinged and clamped axonemes isolated from C. reinhardtii: influence of calcium

Figures 1(A)–(D) illustrates planar swimming motion of an exemplary free axoneme reactivated at 80 μM ATP, and in the absence of external calcium (see video 1). Bending waves initiate at the basal end (shown in blue in figures 1(A) and (B)) and travel toward the distal tip of the axoneme (shown in white in figures 1(A) and (B)) periodically in time, as shown in figure 1(C), with a frequency of f₀ = 18.17 ± 0.01 Hz. These traveling periodic curvature waves provide the necessary thrust to propel the axoneme in the surrounding water-like fluid. This thrust force is balanced with the viscous drag exerted by the fluid on the swimmer [44]. The swimming dynamics of axonemes differ from sperm flagellar propulsion primarily in that axonemes are shorter in length (∼10 μm compared to 50 μm of human spermatozoa) and beat with a non-zero static curvature (figure 2(A)), causing a circular swimming trajectory (figures 1(B) and (D)). To quantify the static curvature, we translate and rotate the axonemal configurations such that the basal end (s = 0) is at the origin of the coordinate system and the tangent vector at s = 0 is oriented along the $\hat{x}$ direction (see figure 1(A)). The filament in cyan color in figure 2(A) shows the time-averaged axonemal shape with mean curvature of κ₀ ∼ −0.21 μm⁻¹. The negative sign of κ₀ indicates a clockwise bend when moving from the basal end at s = 0 toward the distal tip at s = L (figure 2(A)). The base-to-tip bending deformations superimposed on this negative intrinsic curvature cause a counter-clockwise rotation of the axoneme (figure 1(B)). The wavelength λ of these bending waves is larger than the contour length. To calculate λ, we performed a Fourier transform of θ(s, t) in time at each position s along the contour length of the axoneme and obtained the phase ϕ(s) (figure 2(B)). The wavelength is calculated as λ ∼ 2πL/(ϕ(L) − ϕ(0)) = 15.95 μm [15], which is $\sim 30\%$ larger than the axonemal contour length L = 12.35 μm. We notice that the phase ϕ(0) is undefined since θ(s = 0, t) = 0. This is immaterial for the present purpose, as we are only interested in the phase gradient, which we determined by fitting the phase dependence along the axoneme, as indicated by the red line in figure 2(B).

Zoom In Zoom Out Reset image size

The beating frequency of axonemes depends on the ATP concentration following a Michaelis–Menten-type kinetics [15, 43]: it starts with a linear trend at small concentrations of ATP and saturates at higher ATP concentrations around 1 mM (see figure 2(C)). In our experiments, we measured a critical minimum of ATP concentration around 60 μM necessary to reactivate axonemes [15, 21]. Reactivated axonemes swim on circular trajectories effectively in 2D (see reference [45] for a study of the small out-of-plane beating components of isolated axonemes). An active axoneme undergoes planar shape deformations over time, but at any instant of time it may be considered as a solid body with translational and rotational velocities U_x , U_y and Ω_z , which we measure in the swimmer-fixed frame (figures 2(D)–(F)) [46]. These velocities oscillate in time, reflecting the fundamental beat frequency of the axoneme (∼18 Hz) and its higher harmonics.

We perform PCA to describe the time-dependent shape deformation of isolated axonemes [16, 34]. This analysis is based on the method introduced by Stephens et al [34], which was initially used to characterize waveforms in C. elegans. Using the x and y coordinates of the tracked flagellum (see section 4.2), we calculate θ(s, t), which is defined as the angle between the local tangent vector of the centerline of the tracked flagellum and the $\hat{x}$-axis (see figure 1(A)). We recall that we chose the basal end (s = 0) at the origin (0, 0), with the tangent vector at s = 0 along the x-axis (figure 2(A)). We then average the axonemal shapes over one beat cycle to obtain the mean shape of the axoneme ⟨θ(s,t)⟩_t , which corresponds to the arc-shaped filament colored in cyan in figure 2(A).

Next, we compute the covariance matrix θ_Cov(s, s') of the fluctuations of the mean-subtracted angle θ(s, t) − ⟨θ(s,t)⟩_t . Figure 3(A) shows the color map of the flagellum covariance matrix, which has an effective reduced dimensionality with a small number of non-zero eigenvalues. We note that four eigenmodes M_n (s) (n = 1, ..., 4) corresponding to the first four largest eigenvalues of θ_Cov(s, s') capture the shape of the flagellum with very high accuracy (see video 2). These four eigenvectors are plotted in figure 3(B). We reconstruct the axonemal shapes using the four dominant modes and the corresponding time-dependent motion amplitudes a_n (t) as:

$$\begin{equation}\theta (s,t)-{\langle \theta (s,t)\rangle }_{t}\approx \sum\limits _{n=1}^{4}\,{a}_{n}(t){M}_{n}(s).\end{equation} \tag{ 1 }$$

The reconstructed axonemal shapes are presented in figure 3(C), where the green dashed lines show the experimental data. The time-dependent coefficients a₁(t) and a₂(t) are shown in figure 3(D) (a₃(t) and a₄(t) are shown in figures S3 and S4). In order to compare the experimental data with the reconstructed shapes, figures 3(E) and (F) shows two exemplary configurations with the lowest and the highest root mean square errors of 0.006 and 0.085, respectively.

Zoom In Zoom Out Reset image size

Let us focus on the first two modes M₁(s) and M₂(s), which in combination with the first two time-dependent motion amplitudes a₁(t) and a₂(t) capture 98.52% of the total variance (figure 3(G)). The black dashed lines in figure 3(B) represent the Fourier fits of the modes, imposing the boundary condition θ(s = 0, t) = 0:

$$\begin{equation}{M}_{1}(s)=-{b}_{1}+{b}_{1}\,\mathrm{cos}(ks)+{b}_{2}\,\mathrm{sin}(ks),\end{equation} \tag{ 2 }$$

$$\begin{equation}{M}_{2}(s)=-({c}_{1}+{c}_{3}+{c}_{5})+\sum\limits _{m=1}^{3}\left({c}_{2m-1}\,\mathrm{cos}(mks)+{c}_{2m}\,\mathrm{sin}(mks)\right),\end{equation} \tag{ 3 }$$

where k = 2π/λ and the parameters are: λ = 16.87 μm, b₁ = 0.136, b₂ = −0.055, c₁ = 2.3 × 10⁻⁵, c₂ = −0.239, c₃ = 0.003, c₄ = 0.061, c₅ = 0.011, and c₆ = 0.012. In addition, the cyan dashed lines in figure 3(D) fitting the time-dependent coefficients a₁(t) and a₂(t) in equation (1) are given by the Fourier modes as:

$$\begin{equation}{a}_{1}(t)={a}_{11}\,\mathrm{cos}({\omega }_{0}t)+{a}_{12}\,\mathrm{sin}({\omega }_{0}t)\quad \text{and}\;\quad {a}_{2}(t)={a}_{21}\,\mathrm{cos}({\omega }_{0}t)+{a}_{22}\,\mathrm{sin}({\omega }_{0}t),\end{equation} \tag{ 4 }$$

with the coefficients a₁₁ = 5.29, a₁₂ = −1.67, a₂₁ = 0.78, a₂₂ = 2.81 and ω₀ = 2πf₀ = 114.17 ± 0.06 rad s⁻¹. Since the mean-shape ⟨θ(s,t)⟩_t is subtracted, the time-average of both a₁(t) and a₂(t) is zero. Moreover, we note that, as shown in reference [16], the second harmonic can play a crucial role in the rotational speed of the axonemes, but at the level of our approximations, equation (4) provides a reasonable fit (see figure 3(D)). The spatial Fourier decomposition of the modes M₁(s) and M₂(s), as presented in equations (2) and (3), in combination with the time-dependent coefficients a₁(t) and a₂(t) in equation (4) describes the axonemal shapes with good accuracy, therefore the representation proposed here also accurately capture the boundary conditions at s = L. At the distal region of the axoneme, both θ(s = L, t) and its derivative oscillate at the beat frequency f₀, as shown in figure S5.

Having described the angle θ(s, t) in terms of the eigenmodes M₁(s) and M₂(s), we calculate next the curvature waves as:

$$\begin{equation}\kappa (s,t)=\frac{\mathrm{d}\theta (s,t)}{\mathrm{d}s}=\frac{\mathrm{d}{\langle \theta (s,t)\rangle }_{t}}{\mathrm{d}s}+{a}_{1}(t)\frac{\mathrm{d}{M}_{1}(s)}{\mathrm{d}s}+{a}_{2}(t)\frac{\mathrm{d}{M}_{2}(s)}{\mathrm{d}s}.\end{equation} \tag{ 5 }$$

The mean angle ⟨θ(s,t)⟩_t shows a linear trend with contour length s, indicating a constant static curvature κ₀ (cyan filament in figure 2(A)). Next, we define the dimensionless curvature as C = κλ/(2π) and use equations (2)–(4) to express the deviation of C(s, t) from its averaged value C₀ as a sum of propagating waves of wavenumber mk propagating forward and backward, respectively, with an amplitude C_m and ${C}_{m}^{\prime }$. The expressions for C_m and ${C}_{m}^{\prime }$ can be readily derived from equations (2)–(4), and the resulting explicit expressions are given in the supplementary material (table S1). Interestingly, combining the Fourier decomposition in space of M₁ (equation (2)) or M₂ (equation (3)) with the temporal dependence of a₁ or a₂ (equation (4)) results in standing waves. Specifically, only M₂ contributes to standing waves of wavenumbers 2k and 3k, resulting in ${C}_{2}={C}_{2}^{\prime }$ and ${C}_{3}={C}_{3}^{\prime }$. On the other hand, both M₁ and M₂ generate a standing wave at wavenumber k, which implies that the coefficients ${C}_{1}\ne {C}_{1}^{\prime }$. The expression for the spatial and temporal dependence of the curvature can then be arranged as:

$$\begin{equation}C(s,t)-{C}_{0}(s)\approx ({C}_{1}-{C}_{1}^{\prime })\,\mathrm{cos}(ks-{\omega }_{0}t+{\alpha }_{1})+\sum\limits _{m=1}^{3}\,{C}_{m}^{\prime }\left[\mathrm{cos}(mks-{\omega }_{0}t+{\alpha }_{m})+\mathrm{cos}(mks+{\omega }_{0}t+{\alpha }_{m}^{\prime })\right],\end{equation} \tag{ 6 }$$

where C₀(s) = κ₀(s)λ/(2π) is the dimensionless mean-curvature, κ₀ is the mean curvature of the cyan filament in figure 2(A). For our exemplary axoneme in figure 1, the values of these coefficients are listed in the first row of table S2. For comparison, the extracted coefficients of five other free axonemes, all reactivated at [ATP] = 200 μM, are also presented.

We summarize some important properties of equation (6). (i) To reconstruct the axonemal shapes with high precision, in addition to the main spatial mode k, higher modes 2k and 3k, must be included. (ii) The Fourier analysis of the PCA modes shows that the main base-to-tip traveling wave component (with amplitude ${C}_{1}-{C}_{1}^{\prime }$) co-exists with standing waves at wavenumbers mk, with amplitudes ${C}_{m}^{\prime }$. We note that these standing waves may result from modulations of the wave amplitude, as discussed in supplementary section S1.2. (iii) For the standing waves at wavenumber mk (m = 1, 2, 3) the waveform reduces to $2{C}_{m}^{\prime }\,\mathrm{cos}(\omega t+({\alpha }_{m}^{\prime }-{\alpha }_{m})/2)\mathrm{cos}(mks+({\alpha }_{m}+{\alpha }_{m}^{\prime })/2)$, with nodes at locations ${s}_{\text{node}}^{m}={(mk)}^{-1}[(2n+1)\pi /2-({\alpha }_{m}+{\alpha }_{m}^{\prime })/2]$, n being an integer.

We stress that the representation of the wavy motion provided by equation (6) is only an approximation, and that deviations due to higher order PCA modes M_n with n > 2, among others, are expected to lead to ${C}_{m}\ne {C}_{m}^{\prime }$ for m = 2 and 3. Still, given that the two modes considered here account for $\sim 98.52\%$ of the total variance of the signal, the corrections to equation (6) are expected to be small. We also insist that the biological implications of describing the motion in terms of a dominant base-to-tip propagating wave, modulated by a set of standing waves remain to be discussed.

To distinguish the role of different modes in reconstructing the axonemal shapes, we compared the original experimental data (without the mean-shape subtraction) with the shapes reconstructed using the combination of different modes as $\theta (s,t)={{\Sigma}}_{i}{a}_{i}^{\prime }(t){M}_{i}(s)$. Please note that the time-average of ${a}_{i}^{\prime }(t)$ is not zero and ${{\Sigma}}_{i}{\langle {a}_{i}^{\prime }(t)\rangle }_{t}{M}_{i}(s)$ gives the axonemal mean-shape (cyan filament in figure 2(A)). The coefficients ${a}_{i}^{\prime }$ and a_i are related by the simple relation ${a}_{i}^{\prime }(t)={a}_{i}(t)+{\langle {a}_{i}^{\prime }(t)\rangle }_{t}$. This comparison for different number of modes is presented in figure 4.

Zoom In Zoom Out Reset image size

Our PCA analysis was applied to 20 different axonemes. The corresponding eigenmodes exhibit some similar trends, as can be seen by superposing them after re-scaling the contour length s by the total length L or the wavelength λ (see figures S6–S8).

To understand the effect of various modes on the rotational and translational velocities of a freely beating axoneme, we analyzed analytically a simplified waveform. For this purpose, we set the Fourier coefficients c₃, c₄, c₅ and c₆ in equation (3) to zero. These coefficients capture the effect of the higher spatial modes 2k and 3k. The green line in figure S9 shows the deviation between the experimentally extracted mode M₂(s) and the corresponding fitted mode. We note that setting these coefficients to zero results to C'₂=C'₃=0 in equation (6), which leads to the simplified expression:

$$\begin{equation}C(s,t)\approx {C}_{0}+{C}_{1}\,\mathrm{cos}({\omega }_{0}t-ks+{\alpha }_{1})+{C}_{1}^{\prime }\,\mathrm{cos}({\omega }_{0}t+ks+{\alpha }_{1}^{\prime }).\end{equation} \tag{ 7 }$$

The total force and torque exerted on a freely swimming axoneme at low Reynolds number regime is zero. Given the prescribed form of curvature waves in equation (7), we used resistive force theory (RFT) [47, 48] to calculate propulsive forces and torques from the tangential and normal velocity components of each axonemal segment (see section 4.4). By imposing the force-free and torque-free constrains in 2D, we can uniquely determine the translational and rotational velocities of the axoneme. The resulting expressions for the translational and angular velocities, equations (S.3)–(S.5), reflect an important symmetry property of the system. Namely, for the axoneme to rotate, it is necessary to break the mirror symmetry, i.e. C₀ ≠ 0. If the wave pattern of the axoneme is transformed to its mirror-image after half a period (i.e. C(s, t) = −C(s, t + T/2)), then the flagella will swim on average in a straight trajectory. This is the case if C₀ = 0 in equation (7). The reason is that the mean rotations achieved in the first and second half-periods, have the same magnitude but opposite signs and the net rotation sums up to zero. Equations (S.3)–(S.5) also show that in the absence of even harmonics, ⟨U_y ⟩ = 0 and ⟨U_x ⟩ depends on the square of C₁ and ${C}_{1}^{\prime }$ [46, 49, 50]. One notices that in the case of standing waves, ${C}_{1}={C}_{1}^{\prime }$, the net propulsion vanishes. Finally, in the case of λ = L and with the value of η⁻¹ = ζ_⊥/ζ_∥ = 1.8 [48], equations (S.3)–(S.5) simplify to [20]:

$$\begin{equation}\frac{\langle {{\Omega}}_{z}\rangle }{{f}_{0}}\approx -2.23{C}_{0}({C}_{1}^{2}-{C}_{1}^{\prime 2}),\end{equation} \tag{ 8 }$$

$$\begin{equation}\frac{\langle {U}_{x}\rangle }{L{f}_{0}}\approx -0.93({C}_{1}^{2}-{C}_{1}^{\prime 2}),\end{equation} \tag{ 9 }$$

$$\begin{equation}\frac{\langle {U}_{y}\rangle }{L{f}_{0}}\approx 0.26{C}_{0}({C}_{1}^{2}-{C}_{1}^{\prime 2}).\end{equation} \tag{ 10 }$$

To investigate the effect of calcium ions on the flagellar waveform, we performed experiments at different calcium concentrations and systematically measured the static and dynamic components of the axonemal waveform. As shown in figure 5, we changed the calcium concentration from 10⁻⁴ mM to 1 mM and quantified the curvature waves of the reactivated axonemes. Our analysis confirms that the static component of the curvature waves C₀ decreases by $\sim 85\%$, as we increase [Ca²⁺] from 10⁻⁴ mM to 1 mM. The drop is relatively sharp at [Ca²⁺] around 0.02 mM, as shown in figure 5(A). This decrease in C₀ causes a transition from a circular swimming trajectory to a straight trajectory (see figures 5(A) and (B) and video 3), as expected based on equation (S.3). We also observed that the beat frequency of axonemes drops slightly at high calcium concentrations around 1 mM (figure 5(C)). In addition, calcium ions affect the amplitude of the main traveling wave component. To quantify this effect, we examined the wave amplitudes C₁ inferred from our PCA analysis. These results, which are plotted in figure 5(D), show a small decrease in C₁ at high [Ca²⁺].

Zoom In Zoom Out Reset image size

We notice that calcium concentrations lower than ∼0.02 mM do not qualitatively affect the picture drawn in the absence of calcium, as summarized in section 2.1. However, at high calcium concentrations, the rotational velocities of axonemes Ω_z (t) are small and highly noisy, and a direct measurement of ⟨Ω_z ⟩ is difficult. Therefore, in figure 5(B), ⟨Ω_z ⟩ are calculated using the full experimental shapes based on RFT simulations with η⁻¹ = ζ_⊥/ζ_∥ ∼ 1.8. To check the validity of RFT in our system, we used our experimental data at zero or very low concentrations of calcium, where axonemes swim on circular paths and a direct measurement of ⟨Ω_z ⟩ is possible. The comparison between the experimental measurements and RFT simulations is shown in figure 6(A). The black dashed line has the slope of 1 and passes through the origin, which provides a strong evidence for the validity of RFT in our experimental system. Furthermore, with respect to our PCA analysis discussed earlier, we mention that the mean axonemal rotational velocity determined with RFT depends on the number of modes considered for the shape reconstruction. Figure 6(B) shows that ⟨Ω_z ⟩/ω₀ converges for n ⩾ 4. The black dashed lines in figure 6(B) are the values of ${\langle {{\Omega}}_{z}\rangle }^{\text{RFT}}/{\omega }_{0}$ obtained using full experimental shapes, i.e. infinite number of modes.

Zoom In Zoom Out Reset image size

Finally, similar to figure 4, we performed PCA analysis to reconstruct axonemal shapes for an exemplary axoneme at calcium concentration of 0.1 mM. As shown in figure 7, the first four modes recapitulate the beat pattern of the axoneme with high accuracy.

Zoom In Zoom Out Reset image size

In our experiments, attachment of axonemes to the substrate occurred either from the basal end or the distal tip due to the non-specific axoneme-substrate interactions. These axonemes were either hinged from the basal or the distal side and rotated freely around the hinging point, as shown in figure 8 and videos 4–6, or they were clamped from the either side to the substrate and were unable to rotate, as shown in figure 9 and videos 7–10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us first focus on the hinged boundary condition and consider an exemplary case where the axoneme is attached from the basal end and rotates counter-clockwise at the mean rotational velocity of ⟨Ω_z ⟩ ∼ 15 rad s⁻¹ around the pinning point (figures 8(A) and (B)). Curvature waves initiate at the surface-attached proximal region (s = 0) and travel at the frequency of 44.34 ± 0.02 Hz toward the distal tip (s = L); see figure 8(C). This axoneme has an intrinsic curvature of ∼−0.17 μm⁻¹, which is comparable to the corresponding value for the free axoneme presented in figure 1. The amplitude of the traveling wave components, which are inferred from the PCA analysis, are listed in the top panel of table S3 (shapes reconstructed from the PCA analysis are illustrated in figure S10). Comparison of the wave amplitudes shows that the main wave amplitude C₁ is almost 5 times larger than the wave component ${C}_{1}^{\prime }$; see also table S4 for wave amplitudes of two other basal-hinged exemplary axonemes. For comparison, we also present a distal-hinged axoneme which rotates clockwise around the hinging point, as shown in figures 8(D)–(F). The hinging of axonemes from the distal tip to the substrate rarely occurred in our experiments. For this exemplary case, C₁ and ${C}_{1}^{\prime }$ extracted from the PCA analysis show that ${C}_{1}\sim 6{C}_{1}^{\prime }$ (see the second panel of table S3), which is comparable with the values obtained for the basal-hinged axoneme in panels (A)–(C).

To calculate the mean rotational velocity of a model hinged axoneme, we used RFT and the simplified waveform given by equation (7). As noted above, a hinged axoneme was free to rotate but not to translate, thus the total torque exerted on the axoneme is zero while the total force is non-zero. Following the procedure outlined in section 4.5, we calculate the mean rotational velocity of the axoneme up to the first order in C₀ and second order in C₁ and ${C}_{1}^{\prime }$ we arrive at equation (S.6), where in the limit of λ = L and with η⁻¹ = ζ_⊥/ζ_∥ = 1.8, it further simplifies to

$$\begin{equation}\frac{\langle {{\Omega}}_{z}\rangle }{{f}_{0}}\approx -0.63{C}_{0}({C}_{1}^{2}-{C}_{1}^{\prime 2}).\end{equation} \tag{ 11 }$$

Comparison with the pre-factor in equation (8) shows that with the same intrinsic curvature C₀ and wave amplitudes C₁ and ${C}_{1}^{\prime }$, the hinged axoneme rotates on average slower than a free one.

Finally, we present the experimental data of two exemplary non-rotating clamped axonemes where part of the axonemal length, either at the basal or the distal tip, was firmly attached to the substrate. Figures 9(A)–(C) shows a basal-clamped axoneme that beats at the frequency of f₀ = 20.31 ± 0.29 Hz with a reduced intrinsic curvature κ₀ ∼ −0.05  μm⁻¹, compared to the free and hinged axonemes. The amplitude of the main traveling wave component C₁, extracted form PCA analysis is also reduced compared to the free and hinged axonemes. Moreover, the amplitude of the standing wave ${C}_{1}^{\prime }$ is three times smaller than the main base-to-tip propagating wave component C₁ (see the third panel of table S3). Notably, for the case where the axoneme is clamped from the distal tip, the magnitude of the standing wave component becomes larger $({C}_{1}/{C}_{1}^{\prime }\sim 2)$, indicating possible wave reflection from the clamped side of the axoneme (see figures 9(D)–(F) and the forth panel of table S3). In intact C. reinhardtii cells, the attachment of the flagellum to the cell body is close to a basal-clamped boundary condition, so we repeated our PCA analysis for a flagellum in an intact cell (see figures S1(A) and (B) and S11). This analysis gives ${C}_{1}/{C}_{1}^{\prime }\sim 3$, which is comparable to the values mentioned above for the basal-clamped isolated axonemes.

We used wild-type C. reinhardtii cells, strain SAG 11-32b, to isolate flagella using dibucaine following the protocols in references [36, 65]. First, we grew the cells axenically in TAP (tris-acetate-phosphate) medium on a 12 h/12 h day–night cycle. Cells release their flagella upon adding dibucaine, which we purify on a 25% sucrose cushion, and demembranate using detergent NP-40 in HMDEK solution (30 mM HEPES-KOH, 5 mM MgSO₄, 1 mM DTT, 1 mM EGTA, 50 mM potassium acetate, pH 7.4) supplemented with 0.2 mM Pefabloc. The membrane-free axonemes were resuspended in HMDEK buffer plus 1% (w/v) polyethylene glycol (M_w = 20 kg mol⁻¹), 0.2 mM Pefabloc and used freshly after isolation. To perform reactivation experiments, we diluted axonemes in HMDEKP reactivation buffer (HMDEK plus 1% PEG) supplemented with 0.2 mM Pefabloc and different concentrations of ATP. Reactivation solution plus axonemes was infused into 100 μm deep flow chambers, built from cleaned glass and double-sided tape. To avoid axoneme attachment to the substrate, glass slides are coated with 2 mg mL⁻¹ casein solved in HMDEKP. For calcium experiments, HMDEKP reactivation buffer is supplemented with calcium at different concentrations.

We used high-frame phase-contrast microscopy to analyze fast beating dynamics of axonemes (100× objective, 1000 fps). First, we inverted phase-contrast images and subtracted the mean-intensity of the time series to increase the signal to noise ratio [15]. Next, we applied a Gaussian filter to smooth the images. Tracking of axonemes is done using GVF technique [41, 42]. For the first frame, we select a region of interest that should contain only one actively beating axoneme (see figure S2). Then, we initialize the snake by drawing a line polygon along the contour of the axoneme in the first frame. This polygon is interpolated at N equally spaced points and used as a starting parameter for the snake. The GVF is calculated using the GVF regularization coefficient μ = 0.1, as defined in equation (10) of reference [41], with 20 iterations. The snake is then deformed according to the GVF where we have adapted the original algorithm by Xu and Prince for open boundary conditions [41, 42]. We obtain positions of N points along the contour length s of the filament so that s = 0 is the basal end and s = L is the distal side, where L is the total contour length of filament. The position of filament at s_i is denoted by r (s_i ) = (x(s_i ), y(s_i )).

We describe the shape of each flagellum by its unit tangent vector $\hat{\mathbf{t}}(s)$ and the unit normal vector $\hat{\mathbf{n}}(s)$ at distance s along the contour [66, 67]. Instantaneous deformation of flagellum is described by curvature κ(s, t), which using Frenet–Serret formulas is given by [34]

$$\begin{equation}\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s}=\hat{\mathbf{t}}(s)\quad \text{and}\quad \frac{\mathrm{d}\hat{\mathbf{t}}(s)}{\mathrm{d}s}=\kappa (s)\hat{\mathbf{n}}(s).\end{equation} \tag{ 12 }$$

Let us define θ(s) to be the angle between the tangent vector at distance s and the x-axis, then κ(s) = dθ(s)/ds. For shape analysis, we translate and rotate each flagellum such that basal end is at (0, 0) and the orientation of the tangent vector at s = 0 is along the $\hat{x}$-axis i.e. θ(s = 0, t) = 0. Following Stephens et al [34], we performed PCA by calculating the covariance matrix of mean-subtracted angles θ_m (s, t) = θ(s, t) − ⟨θ(s,t)⟩_t , defined as θ_Cov(s, s') = ⟨(θ_m (s, t) − ⟨θ_m ⟩)(θ_m (s', t) − ⟨θ_m ⟩)⟩, where ⟨θ_m ⟩ is the spatial average of θ_m (s, t) at a given time t. The eigenvalues λ_n and the corresponding eigenvectors M_n (s) of the covariance matrix are given by Σ_s' θ_Cov(s, s')M_n (s') = λ_n (t)M_n (s). Figure 3 shows that the superposition of four eigenvectors corresponding to the four largest eigenvalues can describe the flagellum shape with high accuracy:

$$\begin{equation}{\theta }_{m}(s,t)=\theta (s,t)-{\langle \theta (s,t)\rangle }_{t}\approx \sum\limits _{n=1}^{4}\,{a}_{n}(t){M}_{n}(s).\end{equation} \tag{ 13 }$$

Here the four variables a₁(t), ..., a₄(t) are the amplitudes of motion along different principal components and are given by a_n (t) = ∑_s M_n (s)θ_m (s, t). The fractional variance of flagellum's shape captured by n eigenvectors is calculated as ${\sigma }_{n}^{2}={{\Sigma}}_{i=1}^{n}\,{\lambda }_{i}/{\sigma }^{2}$ where ${\sigma }^{2}={{\Sigma}}_{i}^{N}{\lambda }_{i}$. Here N is the total number of eigenvectors. Figure 3(G) shows that already two modes capture 98.52% and four modes capture 99.85% of the total variance.

Biological microorganisms swim with flagella and cilia in the world of 'low Reynolds number' where they experience viscous forces many orders of magnitude larger than inertial forces [44, 68]. In this world where inertia is negligible, Newton's law becomes an instantaneous balance between external and fluid forces and torques exerted on the swimmer, i.e. ${\boldsymbol{\mathcal{F}}}_{\text{ext}}+{\boldsymbol{\mathcal{F}}}_{\text{fluid}}=0$ and ${\boldsymbol{\mathcal{T}}}_{\text{ext}}+{\boldsymbol{\mathcal{T}}}_{\text{fluid}}=0$. The force ${\boldsymbol{\mathcal{F}}}_{\text{fluid}}$ and torque ${\boldsymbol{\mathcal{T}}}_{\text{fluid}}$ exerted by the fluid on the axoneme can be written as

$$\begin{equation}{\boldsymbol{\mathcal{F}}}_{\text{fluid}}={\int }_{0}^{L}\mathrm{d}s\enspace {\boldsymbol{\mathcal{F}}}_{\text{Axoneme}}(s,t),\end{equation} \tag{ 14 }$$

$$\begin{equation}{\boldsymbol{\mathcal{T}}}_{\text{fluid}}={\int }_{0}^{L}\mathrm{d}s\enspace \boldsymbol{r}(s,t)\times {\boldsymbol{\mathcal{F}}}_{\text{Axoneme}}(s,t),\end{equation} \tag{ 15 }$$

where the integrals over the contour length L of the axoneme calculate the total hydrodynamic force and torque exerted by the fluid on the axoneme. ATP-reactivated axonemes show oscillatory shape deformations. At any given time, we consider an axoneme as a solid body with unknown translational and rotational velocities U(t) and Ω(t), yet to be determined. ${\boldsymbol{\mathcal{F}}}_{\text{fluid}}$ and ${\boldsymbol{\mathcal{T}}}_{\text{fluid}}$ can be separated into propulsive part due to the relative shape deformation of the axoneme in the body-fixed frame and the drag part [69]

$$\begin{equation}\left(\begin{matrix}\hfill {\boldsymbol{\mathcal{F}}}_{\text{fluid}}\hfill \\ \hfill {\boldsymbol{\mathcal{T}}}_{\text{fluid}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\boldsymbol{\mathcal{F}}}^{\,\text{prop}}\hfill \\ \hfill {\boldsymbol{\mathcal{T}}}^{\,\text{prop}}\hfill \end{matrix}\right)-{\mathcal{D}}_{\mathrm{A}}\left(\begin{matrix}\hfill \boldsymbol{U}\hfill \\ \hfill \mathbf{\Omega }\hfill \end{matrix}\right),\end{equation} \tag{ 16 }$$

where 6 × 6 geometry-dependent drag matrix of the axoneme ${\mathcal{D}}_{\mathrm{A}}$ is symmetric and nonsingular (invertible). We also note that a freely swimming axoneme experiences no external forces and torques, thus ${\boldsymbol{\mathcal{F}}}_{\text{fluid}}$ and ${\boldsymbol{\mathcal{T}}}_{\text{fluid}}$ must vanish. Further, since swimming effectively occurs in 2D, ${\mathcal{D}}_{\mathrm{A}}$ is reduced to a 3 × 3 matrix and equation (16) can be reformulated as

$$\begin{equation}\left(\begin{array}{c}{U}_{x}\\ {U}_{y}\\ {{\Omega}}_{z}\end{array}\right)={\mathcal{D}}_{\mathrm{A}}^{-1}\left(\begin{array}{c}{\mathcal{F}}_{x}^{\;\text{prop}}\\ {\mathcal{F}}_{y}^{\;\text{prop}}\\ {\mathcal{T}}_{z}^{\;\text{prop}}\end{array}\right)={\left(\begin{array}{ccccccccc}{d}_{11}& {d}_{12}& {d}_{13}\\ {d}_{21}& {d}_{22}& {d}_{23}\\ {d}_{31}& {d}_{32}& {d}_{33}\end{array}\right)}^{-1}\left(\begin{array}{c}{\mathcal{F}}_{x}^{\;\text{prop}}\\ {\mathcal{F}}_{y}^{\;\text{prop}}\\ {\mathcal{T}}_{z}^{\;\text{prop}}\end{array}\right),\end{equation} \tag{ 17 }$$

which we use to calculate translational and rotational velocities of the swimmer after determining the drag matrices ${\mathcal{D}}_{\mathrm{A}}$ and the propulsive forces and torque $({\mathcal{F}}_{x}^{\,\text{prop}},{\mathcal{F}}_{y}^{\,\text{prop}},{\mathcal{T}}_{z}^{\,\text{prop}})$.

We calculate ${\mathcal{F}}_{x}^{\,\text{prop}}$, ${\mathcal{F}}_{y}^{\,\text{prop}}$ and ${\mathcal{T}}_{z}^{\,\text{prop}}$ which are propulsive forces and torque due to the shape deformations of the axoneme in the swimmer-fixed frame by selecting the basal end of the axoneme as the origin of the swimmer-fixed frame. As shown in figure 1(A), we define the local tangent vector at contour length s = 0 as $\hat{X}$ direction, the local normal vector $\hat{n}$ as $\hat{Y}$ direction, and assume that $\hat{z}$ and $\hat{Z}$ are parallel. Let us define θ₀(t) = θ(s = 0, t) as the angle between $\hat{x}$ and $\hat{X}$ which gives the velocity of the basal end in the laboratory frame as ${U}_{x}^{\text{Basal}-\text{Lab}}=\mathrm{c}\mathrm{o}\mathrm{s}\,{\theta }_{0}(t){U}_{x}+\mathrm{sin}\,{\theta }_{0}(t){U}_{y}$ and ${U}_{y}^{\text{Basal}-\text{Lab}}=-\mathrm{s}\mathrm{i}\mathrm{n}\,{\theta }_{0}(t){U}_{x}+\mathrm{cos}\,{\theta }_{0}(t){U}_{y}$. Furthermore, we note that the instantaneous velocity of the axoneme in the lab frame is given by u = U + Ω × r(s, t) + u', where u' is the deformation velocity of flagella in the swimmer-fixed frame, U = (U_x , U_y , 0) and Ω = (0, 0, Ω_z ) with Ω_z = dθ₀(t)/dt.

To calculate ${\mathcal{F}}_{x}^{\,\text{prop}}$, ${\mathcal{F}}_{y}^{\,\text{prop}}$ and ${\mathcal{T}}_{z}^{\,\text{prop}}$ for a given beating pattern of axoneme in the body-fixed frame, we used classical framework of RFT which neglects long-range hydrodynamic interactions between different parts of the axoneme [47, 48]. In this theory, axoneme is divided into small cylindrical segments moving with velocity u'(s, t) in the body-frame and propulsive force ${\mathcal{F}}^{\,\text{prop}}$ is proportional to the local centerline velocity components of each segment in parallel and perpendicular directions

$$\begin{equation}\begin{aligned}\hfill {\boldsymbol{\mathcal{F}}}^{\text{prop}}(s,t)& ={\zeta }_{{\Vert}}{\mathbf{u}}_{{\Vert}}^{\prime }(s,t)+{\zeta }_{\perp }{\mathbf{u}}_{\perp }^{\prime }(s,t),\hfill \\ \hfill {\mathbf{u}}_{{\Vert}}^{\prime }(s,t)& =[\dot{\mathbf{r}}(s,t).\mathbf{t}(s,t)]\mathbf{t}(s,t),\hfill \\ \hfill {\mathbf{u}}_{\perp }^{\prime }(s,t)& =\dot{\mathbf{r}}(s,t)-{\mathbf{u}}_{{\Vert}}^{\prime }(s,t),\hfill \end{aligned}\end{equation} \tag{ 18 }$$

where ${\mathbf{u}}_{{\Vert}}^{\prime }$ and ${\mathbf{u}}_{\perp }^{\prime }$ are the projections of the local velocity on the directions parallel and perpendicular to the axoneme. The friction coefficients in parallel and perpendicular directions are defined as ζ_∥ = 4πμ/(ln(2L/a) + 0.5) and ζ_⊥ = 2ζ_∥, respectively. This anisotropy indicates that to obtain the same velocity, one would need to apply a force in the perpendicular direction twice as large as that in the parallel direction [47, 48]. Axoneme is a filament of length L ∼ 10 μm and radius a ∼ 100 nm. For a water-like environment with viscosity μ = 0.96 pN ms μm^{− 2}, we obtain ζ_∥ ∼ 2.1 pN ms μm^{− 2}.

Consider an axoneme which is attached from the basal end to the substrate and rotates around the hinging point (see figure 8(A)). We select the hinging point of the axoneme as the origin of the swimmer-fixed frame, and define tangent and normal vectors at s = 0 as the coordinate system in the swimmer-fixed frame. Since the axoneme is free to rotate but not to translate, the total torque ${\mathcal{T}}_{z}^{\text{A}}$ exerted on the axoneme is zero but the total force ${\mathcal{F}}_{x}^{\text{A}}$ and ${\mathcal{F}}_{y}^{\text{A}}$ is non-zero. Reformulating in terms of equation (16), we have

$$\begin{equation}\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\text{A}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\text{A}}\hfill \\ \hfill {\mathcal{T}}_{z}^{\text{A}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{T}}_{z}^{\,\text{prop}}\hfill \end{matrix}\right)+\left(\begin{matrix}\hfill {d}_{11}\hfill & \hfill {d}_{12}\hfill & \hfill {d}_{13}\hfill \\ \hfill {d}_{21}\hfill & \hfill {d}_{22}\hfill & \hfill {d}_{23}\hfill \\ \hfill {d}_{31}\hfill & \hfill {d}_{32}\hfill & \hfill {d}_{33}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {U}_{x}\hfill \\ \hfill {U}_{y}\hfill \\ \hfill {{\Omega}}_{z}\hfill \end{matrix}\right)\end{equation} \tag{ 19 }$$

which by imposing the constrains of U_x = 0, U_y = 0 and ${\mathcal{T}}_{z}^{\text{A}}=0$ simplifies to

$$\begin{equation}\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\text{A}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\text{A}}\hfill \\ \hfill 0\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{T}}_{z}^{\,\text{prop}}\hfill \end{matrix}\right)+\left(\begin{matrix}\hfill {d}_{11}\hfill & \hfill {d}_{12}\hfill & \hfill {d}_{13}\hfill \\ \hfill {d}_{21}\hfill & \hfill {d}_{22}\hfill & \hfill {d}_{23}\hfill \\ \hfill {d}_{31}\hfill & \hfill {d}_{32}\hfill & \hfill {d}_{33}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill {{\Omega}}_{z}\hfill \end{matrix}\right).\end{equation} \tag{ 20 }$$

Using this equation, we can calculate the instantaneous rotational velocity of a hinged axoneme as

$$\begin{equation}{\mathcal{T}}_{z}^{\,\text{prop}}=-{d}_{33}{{\Omega}}_{z},\end{equation} \tag{ 21 }$$

where

$$\begin{equation}{\mathcal{T}}_{z}^{\,\text{prop}}(t)={\int }_{0}^{L}\mathrm{d}s\enspace \mathbf{r}(s,t)\times \mathbf{f}(s,t).\end{equation} \tag{ 22 }$$

Here f(s, t) is calculated using equation (18) based on the simplified waveform defined in equation (7).

To calculate mean rotational velocity for a hinged axoneme, we average over one beat cycle

$$\begin{equation}\left\langle {{\Omega}}_{z}\right\rangle /{\omega }_{0}=-\frac{1}{2\pi }{\int }_{0}^{1/{f}_{0}}\mathrm{d}t\enspace \frac{{\mathcal{T}}_{z}^{\,\text{prop}}(t)}{{d}_{33}(t)}.\end{equation} \tag{ 23 }$$

Please note that for a clamped boundary condition, since the mean translational and rotational velocities of the axoneme are zero, the total torque $({\mathcal{T}}_{z}^{\text{A}})$ and the total force (${\mathcal{F}}_{x}^{\text{A}}$ and ${\mathcal{F}}_{y}^{\text{A}}$) are non-zero, and are equal to the propulsive force and torque

$$\begin{equation}\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\text{A}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\text{A}}\hfill \\ \hfill {\mathcal{T}}_{z}^{\text{A}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\mathcal{F}}_{x}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{F}}_{y}^{\,\text{prop}}\hfill \\ \hfill {\mathcal{T}}_{z}^{\,\text{prop}}\hfill \end{matrix}\right).\end{equation} \tag{ 24 }$$