1. Introduction
The unicellular green alga Chlamydomonas reinhardtii has been widely used as a model organism for studying various problems in fluid dynamics and cell biology (Goldstein Reference Goldstein2015), such as tactic behaviours (Leptos et al. Reference Leptos, Chioccioli, Furlan, Pesci and Goldstein2023), flagellar dynamics (Quaranta, Aubin-Tam & Tam Reference Quaranta, Aubin-Tam and Tam2015; Sartori et al. Reference Sartori, Geyer, Scholich, Jülicher and Howard2016) and single-cell flow fields (Drescher et al. Reference Drescher, Goldstein, Michel, Polin and Tuval2010; Wei et al. Reference Wei, Dehnavi, Aubin-Tam and Tam2019, Reference Wei, Dehnavi, Aubin-Tam and Tam2021). The wild-type cell of C. reinhardtii has two flagella anchored at the anterior side of the spheroidal cell body. During swimming, the two flagella beat approximately in a two-dimensional plane, with a breaststroke-like gait consisting of distinct power and recovery strokes. It is now well established that flagellar beating is driven by distributed molecular motors, a characteristic conserved across eukaryotic cells. Studies of the swimming mechanisms of green algae have inspired numerous designs of artificial flagella and microswimmers (see e.g. Diaz et al. Reference Diaz, Robinson, Aydin, Aydin, Goldman and Wan2021; Moreau et al. Reference Moreau, Walker, Poon, Soto, Goldman, Gaffney and Wan2024).
The analysis of eukaryotic flagellar dynamics often assumes a fixed flagellar base in space. Machin (Reference Machin1958) solved the bending waves of an elastic beam actuated at the boundary. Although this boundary-driven assumption is a significant simplification of realistic flagellar dynamics, later developments reveal the effects of various physical elements on flagellar propulsion, such as effective flexibility (Wiggins & Goldstein Reference Wiggins and Goldstein1998; Yu, Lauga & Hosoi Reference Yu, Lauga and Hosoi2006; Peng, Elfring & Pak Reference Peng, Elfring and Pak2017), swimmer geometry (Lauga Reference Lauga2007), flagellar curvature (Liu et al. Reference Liu, Qin, Zhu, Yang and Luo2020; Hu, Zhang & Shelley Reference Hu, Zhang and Shelley2022) and hydrodynamic interactions (HIs) (Elfasi, Elimelech & Gat Reference Elfasi, Elimelech and Gat2018; Hu & Meng Reference Hu and Meng2024). On the other hand, studies that combine mechanical models of flagellar deformation with molecular motor dynamics focus on how beating waveforms are regulated (see the review by Gilpin, Bull & Prakash Reference Gilpin, Bull and Prakash2020).
Under dynamic swimming conditions, the modelling of microorganisms involves different complexities, especially in the coupling between the cell body and flagella. In addition to mechanical coupling at the flagellar base, flagella are also coupled to the cell body through long-range HIs. Therefore, combining flagellar beating with a moving cell body is crucial for understanding the swimming mechanisms. Simulations and experimental measurements on various microorganisms have reported a decrease in swimming speed with increasing body size, including spermatozoa (Higdon Reference Higdon1979), monotrichous bacteria (Kamdar et al. Reference Kamdar2023; Liu, Chen & Xu Reference Liu, Chen and Xu2025), and the ciliate Volvox (Omori, Ito & Ishikawa Reference Omori, Ito and Ishikawa2020). In particular, body–flagella HIs have been shown to influence both the swimming speed and efficiency of bacteria (Liu et al. Reference Liu, Chen and Xu2025). For flagellated microalgae, we expect more complex and stronger body–flagella HIs due to significantly time-varying and closer spacing between them. Over a beating period, the cell body moves back and forth, entraining a considerable volume of fluid. At low Reynolds number (
$\textit{Re}$
), this body-induced flow decays slowly in space and can significantly influence flagellar propulsion. Meanwhile, the flow induced by the flagella also modifies the drag force acting on the cell body. As a minimal swimming model, the three-sphere model, in which the flagellar beating is represented by spheres moving along closed orbits, provides insights into flagellar synchronisation (Friedrich & Jülicher Reference Friedrich and Jülicher2012), run-and-tumble behaviour (Bennett & Golestanian Reference Bennett and Golestanian2013), and three-dimensional helical trajectories (Cortese & Wan Reference Cortese and Wan2021). Optimal gaits for swimming and feeding have been identified using numerical optimisation (Tam & Hosoi Reference Tam and Hosoi2011), where the swimming dynamics is modelled via flow singularities placed along the flagella and at the centre of the cell body. Although these studies reveal different dynamic behaviours of microalgal swimming, the effect of body–flagella HIs on swimming performance remains elusive. In recent work (Hu & Meng Reference Hu and Meng2024), we modelled the swimming dynamics of multiflagellated microalgae, and demonstrated that body–flagella HIs determine the hydrodynamic advantage of multiflagellarity. In particular, swimming efficiency increases monotonically with the number of flagella only when they are attached to the anterior side of the cell body. However, whether an optimal body size exists for microalgal swimming is still an open question.
In this work, we use numerical simulations to investigate the effect of body–flagella HIs on the swimming performance of C. reinhardtii. We prescribe the flagellar waveform based on experimental measurements, and incorporate both body–flagella and interflagellar HIs using a boundary element method. We find that body–flagella HIs can significantly enhance the swimming performance and lead to an optimal body size for both speed and efficiency. Further analysis based on the three-sphere model reveals that body–flagella HIs are effectively non-reciprocal.
2. Simulation model
Our swimming model of C. reinhardtii consists of two flagella of length
$L$
, and a spheroidal cell body with semi-minor axis
$a$
and semi-major axis
$b$
. Planar waveforms of cells held stationary are measured experimentally via high-speed videography (figure 1
a) and serve as inputs for numerical simulations. Physical parameters of a sample cell are given in table 1. We parametrise each flagellum by the arc length
$s\in [0,L]$
along its centreline. From the measured waveform, we evaluate the velocities of the two flagella relative to the cell body,
$\boldsymbol{U}_{\!i}(s,t)$
for
$i = 1,2$
, at each recorded time instant. In this work, we focus on the normal breaststroke gait. Although a small phase lag (approximately
$0.1\pi$
) exists between the cis and trans flagella (Wan, Leptos & Goldstein Reference Wan, Leptos and Goldstein2014), only a sufficiently large phase difference (
$\gtrsim \pi /2$
) can lead to a non-negligible rotation of the cell body (Geyer et al. Reference Geyer, Jülicher, Howard and Friedrich2013). Therefore, the two flagella are assumed to beat mirror-symmetrically, and the swimming model undergoes only translational motion. We further neglect the non-planar components of flagellar beating, as the resulting rotation frequency of the cell is much smaller than the flagellar beating frequency (Cortese & Wan Reference Cortese and Wan2021).
Parameters of a representative C. reinhardtii cell measured in experiments. The value of flagellum radius is taken from Sager & Palade (Reference Sager and Palade1957).
We simulate the dynamics of the model swimmer using a hybrid boundary element and regularised Stokeslet method (Smith Reference Smith2009). The regularised Stokeslet is the exact solution to the Stokes equation subjected to a smoothed point force. With a specific smooth function, it can be expressed as (Cortez, Fauci & Medovikov Reference Cortez, Fauci and Medovikov2005)
\begin{align} \unicode{x1D642}_{\epsilon } (\boldsymbol{x}_0, \boldsymbol{x}) = \frac {\unicode{x1D644} \big(R^2 + 2 \epsilon ^2 \big) + \boldsymbol{R}\boldsymbol{R}}{R_{\epsilon }^3}, \end{align}
where
$\epsilon$
is the regularisation parameter,
$\boldsymbol{R} = \boldsymbol{x}_0-\boldsymbol{x}$
, and
$R_{\epsilon } = \sqrt {R^2 + \epsilon ^2}$
.
Model of swimming C. reinhardtii. (a) Flagellar waveform measured from experiments. Blue indicates recovery stroke, and red indicates power stroke. (b) Time-averaged disturbance flow field generated by the swimming model.
We represent the motions of the cell body and the flagella by distributions of regularised Stokeslets on the body surface and along the flagellar centrelines, respectively. As verified by a test problem on a rigid scallop in Appendix A, this flagellar representation yields swimming dynamics that closely matches that obtained from a non-local slender body theory (Johnson Reference Johnson1980). Denote the body surface by
$S$
and the flagella centrelines by
$C_i$
with
$i=1,2$
. The disturbance velocity in the fluid domain generated by the motion of the cell body is given by
\begin{align} \boldsymbol{u}_{{b}}(\boldsymbol{x}_0) = \frac {1}{8\pi \mu } \int _{S} \unicode{x1D642}_{\epsilon }(\boldsymbol{x}_0, \boldsymbol{x}) \boldsymbol{\cdot } \boldsymbol{f}_{\!{b}}(\boldsymbol{x})\, {\rm d}\boldsymbol{x}, \end{align}
and that by the flagella is given by
\begin{align} \boldsymbol{u}_{i}(\boldsymbol{x}_0) = \frac {1}{8\pi \mu } \int _{C_i} \unicode{x1D642}_{\epsilon }[\boldsymbol{x}_0, \boldsymbol{x}_i(s)] \boldsymbol{\cdot } \boldsymbol{f}_{i}(s)\, {\rm d}s, \quad \text{for} \quad i = 1,2, \end{align}
where
$\mu$
is the fluid viscosity,
$\boldsymbol{x}_i(s)$
is the centreline position, and
$\boldsymbol{f}_{\!{b}}$
and
$\boldsymbol{f}_{i}$
are the hydrodynamic force densities exerted by the body and flagella on the fluid, respectively. When the evaluation points
$\boldsymbol{x}_0$
lie on the boundaries, the total disturbance velocity due to the flagella and cell body satisfies the following constraints imposed by the no-slip condition:
\begin{align}\boldsymbol{u}_{{b}}(\boldsymbol{x}_0) + \sum _{i} \boldsymbol{u}_{i}(\boldsymbol{x}_0) = \boldsymbol{U}_{\!{b}}(t), \quad \text{for } \quad \boldsymbol{x}_0 \in {S},\end{align}
\begin{align}\boldsymbol{u}_{{b}}(\boldsymbol{x}_0) + \sum _{j} \boldsymbol{u}_{\!j}(\boldsymbol{x}_0) = \boldsymbol{U}_{\!{b}}(t) + \boldsymbol{U}_{\!i}(s,t), \quad \text{for } \quad \boldsymbol{x}_0 \in {C_i} \quad \text{ and } \quad i,j = 1,2,\end{align}
where
$\boldsymbol{U}_{\!{b}}(t)$
is the instantaneous velocity of the cell body. Equations (2.4) and (2.5) incorporate both the interflagellar and the body–flagella HIs. Finally, the total hydrodynamic force and torque are zero:
\begin{align} \int _{S}\boldsymbol{f}_{\!{b}}(\boldsymbol{x})\, {\rm d}\boldsymbol{x} + \sum _i \int _{C_i} \boldsymbol{f}_{\kern-1pt i}(s) \, {\rm d}s = 0, \end{align}
\begin{align} \int _{S} (\boldsymbol{x}-\boldsymbol{x}_{{b}}) \times \boldsymbol{f}_{\!{b}}(\boldsymbol{x})\, {\rm d}\boldsymbol{x} + \sum _i \int _{C_i} [\boldsymbol{x}_i(s)-\boldsymbol{x}_{{b}}] \times \boldsymbol{f}_{\kern-1pt i}(s) \, {\rm d}s = 0, \end{align}
where the torques are evaluated relative to the body centre
$\boldsymbol{x}_{{b}}$
.
We scale lengths by
$L$
, time by the beating period
$T$
, velocities by
$U_0 = L/T$
, and forces by
$F_0 = \mu L^2/T$
. The body aspect ratio
$b/a$
is kept constant when varying the semi-major axis
$b$
. The integral formulations (2.2)–(2.7) are discretised into a system of linear equations for the unknowns
$\boldsymbol{U}_{\!{b}}$
,
$\boldsymbol{f}_{\!{b}}$
and
$\boldsymbol{f}_{\kern-1pt i}$
. Since the Stokeslet
$\unicode{x1D642}_{\epsilon }(\boldsymbol{x}_0, \boldsymbol{x})$
often varies more rapidly near the source position
$\boldsymbol{x}$
than the hydrodynamic forces, the evaluation of the Stokeslet is performed on a finer grid than that of the forces (Smith Reference Smith2009). Equations (2.2) and (2.3) are approximated as
\begin{align} \boldsymbol{u} (\boldsymbol{x}_0) = \frac {1}{8\pi \mu } \sum _{n=1}^N \boldsymbol{f}^n\boldsymbol{\cdot } \int _{\delta S^n (\delta C^n)} \unicode{x1D642}_{\epsilon }(\boldsymbol{x}_0, \boldsymbol{x}) \,{\rm d}\boldsymbol{x}, \end{align}
where
$N$
is the number of discrete surface elements
$\delta S^n$
(line elements
$\delta C^n$
) on the cell body (on each flagellum) with force density
$\boldsymbol{f}^n$
. The integral of
$\unicode{x1D642}_{\epsilon }$
over each element is evaluated using Gauss–Legendre quadrature with
$m\times m$
points on
$\delta S^n$
(
$m$
points on
$\delta C^n$
). For most simulations, we use
$N \approx 200$
on the sphere,
$N = 25$
on each flagellum, and
$m=6$
. The regularisation parameter
$\epsilon$
is taken to be the same as the flagellum radius
$r_0$
(table 1).
3. Results and discussion
3.1. Optimal cell body size for swimming
Figure 1(b) shows the disturbance velocity field of the model swimmer averaged over a beating period. There is one vortex on each side of the cell body, and a flow stagnation point exists in front of the swimmer along its swimming direction. These features closely resemble the measured flow field around free-swimming cells (Drescher et al. Reference Drescher, Goldstein, Michel, Polin and Tuval2010). Using the parameters of the sample cell (table 1), we compute the time-averaged swimming speed
$\langle U_{{b}} \rangle \approx 70.4$
$\unicode{x03BC}$
m s−1, in agreement with the experimentally measured speed (Buchner et al. Reference Buchner, Muller, Mehmood and Tam2021). Slight variations in the anchored angle
$\phi$
of the flagella (figure 1
a) have little effect on the swimming performance.
Effect of cell body size on swimming performance: (a)
$\langle U_{{b}}\rangle /U_{{0}}$
as a function of
$b/L$
; (b) swimming efficiency
$\eta$
as a function of
$b/L$
. The shaded areas indicate the distribution of
$b/L$
from experimental measurements. The vertical dashed lines mark the value of
$b/L$
of the sample cell. Other parameters are given in table 1.
We first evaluate the dependence of
$\langle U_{{b}} \rangle$
on cell body size by varying the semi-major axis length
$b$
while using the same flagellar waveform data. As shown in figure 2(a),
$\langle U_{{b}} \rangle$
varies non-monotonically with
$b/L$
, the ratio of body size to flagellar length. As the cell body becomes more slender (larger values of
$b/a$
), the maximum of
$\langle U_{{b}} \rangle$
shifts to larger values of
$b/L$
. From a collection of approximately 40 cells, we measure an average
$b/L = 0.38 \pm 0.08$
, as indicated by the shaded regions in figure 2. For the sample cell (
$b/a = 1.53$
), the maximum occurs at
$b/L \approx 0.21$
, which is lower than its actual geometric ratio (0.425; see table 1) and lies more than one standard deviation below the measured average.
We further compute the swimming efficiency
$\eta$
using a definition similar to that in Lighthill (Reference Lighthill1952):
\begin{align} \eta = \frac {\zeta _{{b}} \langle U_{{b}} \rangle ^2}{\left \langle \int _{S} \boldsymbol{f}_{\!{b}} \boldsymbol{\cdot } \boldsymbol{U}_{\!{b}} \, {\rm d}\boldsymbol{x} \right \rangle + \sum _i \left \langle \int _{C_i} \boldsymbol{f}_{\kern-1pt i} \boldsymbol{\cdot } (\boldsymbol{U}_{\!{b}} + \boldsymbol{U}_{\!i})\,{\rm d}s \right \rangle }, \end{align}
where
$\zeta _{{b}}$
is the friction coefficient of a spheroid translating parallel to its long axis (Kim & Karrila Reference Kim and Karrila2013). Equation (3.1) expresses the ratio of the work required to drag a spheroid moving steadily at speed
$\langle U_{{b}} \rangle$
to the work dissipated by the motions of the flagella and cell body. Figure 2(b) shows that
$\eta$
also exhibits a pronounced optimum with respect to
$b/L$
. Unlike
$\langle U_{{b}} \rangle$
, for the sample cell, the optimal value of
$b/L$
for
$\eta$
lies within one standard deviation of the measured average.
Removing the interflagellar HIs, which amounts to neglecting the disturbance velocity from the other flagellum in (2.5), has only a minor effect on the swimming performance (see the right-hand images in figure 2
a,b). This is due to hydrodynamic screening by the cell body. Further removing the body–flagella HIs, i.e. neglecting
$\boldsymbol{u}_i$
in (2.4) and
$\boldsymbol{u}_{{b}}$
in (2.5), results in a much lower
$\langle U_{{b}} \rangle$
and a monotonically decreasing trend with
$b$
. Therefore, the surprising increase of
$\langle U_{{b}} \rangle$
at small values of
$b$
is attributed to increasingly strong body–flagella HIs as
$b$
increases. Meanwhile, the intrinsic viscous drag of the cell body
$\zeta _{{b}} U_{{b}}$
, also increases with
$b$
and eventually overwhelms the flagellar propulsion, leading to a decrease in
$\langle U_{{b}} \rangle$
at sufficiently large
$b$
.
Effect of body-flagella HIs on (a) instantaneous swimming speed
$U_{{b}}(t)$
and (b) instantaneous displacement
$y_{{b}}(t)$
. Insets show the differences between the results obtained with and without body–flagella HIs for (a)
$U_{{b}}(t)$
and (b)
$y_{{b}}(t)$
over a beating period. Red and blue regions indicate power and recovery strokes, respectively.
In figure 3, we compare the swimming speed
$U_{{b}}(t)$
and the instantaneous displacement
$y_{{b}}(t)$
of the cell body, obtained with and without body–flagella HIs. Our key observation is that body–flagella HIs not only increase the forward speed during the power stroke, but also reduce the backward speed during the recovery stroke (figure 3
a); as a result, the swimmer achieves a larger forward displacement during the power stroke, and a smaller backward displacement during the recovery stroke (figure 3
b). Over one period, the swimmer accumulates a net displacement of approximately
$0.1L$
, which is approximately twice that achieved without body–flagella HIs.
The appearance of an optimal body size for swimming is robust and not sensitive to variations in flagellar waveforms. In Appendix B, we present additional simulation results using waveforms measured from different sample cells. Despite large variations in the static curvature, wavelength and amplitude, the swimming performance consistently varies non-monotonically with
$b/L$
. This suggests that body–flagella HIs – reflecting how they influence each other’s drag forces – are predominantly determined by their relative large-scale approaching and receding motions on the order of
$L$
during the power and recovery strokes. The robustness also implies that our results are not likely to be affected by slight flagellar adaptations in vivo arising from flexibility and active forcing.
3.2. Three-sphere model of Chlamydomonas swimming
To reveal the mechanisms underlying the enhancement of swimming performance by body–flagella HIs, we consider a three-sphere model consisting of a central sphere of radius
$r_{{b}}$
, representing the cell body, and two side spheres of radius
$r_{{\kern-1.5pt f}}$
, representing the flagella (figure 4
a). The flagellar spheres are assumed to be much smaller than the cell body,
$r_{{\kern-1.5pt f}} \ll r_{{b}}$
. Unlike previous works (Friedrich & Jülicher Reference Friedrich and Jülicher2012; Bennett & Golestanian Reference Bennett and Golestanian2013), we explicitly account for the no-slip boundary condition on the body surface by using the spherical image system
$\unicode{x1D642}_{{im}}$
(Kim & Karrila Reference Kim and Karrila2013). The two flagellar spheres, located at
$\boldsymbol{x}_i$
(
$i = 1,2$
), are prescribed to move with velocity
$v_0 \hat {\boldsymbol{t}}_i$
relative to the body, and are mirror-symmetric with respect to the model centreline. Due to the force-free condition, the cell body, located at
$\boldsymbol{x}_{{b}}$
, develops an instantaneous translational velocity
$\boldsymbol{U}_{\!{b}}$
. The motion of the cell body induces a disturbance flow
$\boldsymbol{v}_{{b}}$
in the surrounding fluid,
\begin{align} \boldsymbol{v}_{{b}}(\boldsymbol{x}_0, \boldsymbol{x}_{{b}}; \boldsymbol{U}_{\!{b}}, r_{{b}}) = \frac {3}{4} r_{{b}} \boldsymbol{U}_{\!{b}} \boldsymbol{\cdot } \bigg (1 + \frac {r_{{b}}^2}{6}{\nabla} ^2 \bigg) \unicode{x1D642}(\boldsymbol{x}_0, \boldsymbol{x}_{{b}}), \end{align}
where
$\unicode{x1D642}$
denotes the free-space Stokeslet.
The total drag force on the flagellar spheres is
$\boldsymbol{F}_{\kern-1pt i} = -6\pi \mu r_{{\kern-1.5pt f}} (v_0 \hat {\boldsymbol{t}}_i + \boldsymbol{U}_{\!{b}}) + \boldsymbol{F}_{{b} \to i}$
, where
$\boldsymbol{F}_{\!{b}\to i}$
is the contribution to
$\boldsymbol{F}_{\kern-1pt i}$
by the cell body, i.e. the body-to-flagella (B-to-F) force, given by
\begin{align} \boldsymbol{F}_{\!{b}\to i} = 6\pi \mu r_{{\kern-1.5pt f}}\! \left [\boldsymbol{v}_{{b}}(\boldsymbol{x}_i, \boldsymbol{x}_{{b}}; \boldsymbol{U}_{\!{b}}, r_{{b}}) + \frac {1}{8\pi \mu } \unicode{x1D642}_{{im}}(\boldsymbol{x}_i,\boldsymbol{x}_{{b}}; r_{{b}})\boldsymbol{\cdot } \boldsymbol{F}_{\kern-1pt i} \right ]\!. \end{align}
In (3.3),
$\boldsymbol{v}_{{b}}$
represents the dynamic component of the B-to-F force, and exists only when the cell body is in motion. The term involving
$\unicode{x1D642}_{{im}}$
represents the static component of the B-to-F force, and captures the effect of a nearby no-slip boundary on the flagellar force. The total drag force on the cell body is
$\boldsymbol{F}_{{b}} = -6\pi \mu r_{{b}} \boldsymbol{U}_{\!{b}} + \sum _i \boldsymbol{F}_{i \to {b}}$
, where
$\boldsymbol{F}_{i \to {b}}$
is the contribution to the drag force acting on the cell body by the flagella-induced flow, i.e. the flagella-to-body (F-to-B) force. Using the strength of the image Stokeslet in
$\unicode{x1D642}_{{im}}$
,
$\boldsymbol{F}_{i \to {b}}$
can be expressed as (Higdon Reference Higdon1979; Kim & Karrila Reference Kim and Karrila2013)
\begin{align} \boldsymbol{F}_{i \to {b}} = -\left [\left (\frac {3 r_{{b}}}{2R_i} - \frac {r_{{b}}^3}{2 R_i^3} \right ) \hat {\boldsymbol{R}}_i\hat {\boldsymbol{R}}_i + \left (\frac {3 r_{{b}}}{4 R_i} + \frac {r_{{b}}^3}{4 R_i^3}\right )(\unicode{x1D644} - \hat {\boldsymbol{R}}_i\hat {\boldsymbol{R}}_i)\right ] \boldsymbol{\cdot } \boldsymbol{F}_{\kern-1pt i}, \end{align}
where
$\boldsymbol{R}_i = \boldsymbol{x}_i {-} \boldsymbol{x}_{{b}}$
and
$\hat {\boldsymbol{R}}_i = \boldsymbol{R}_i/|\boldsymbol{R}_i|$
. The body velocity
$\boldsymbol{U}_{\!{b}}$
can be obtained from the force-free condition
$\boldsymbol{F}_{{b}} + \sum _i \boldsymbol{F}_{\kern-1pt i} = 0$
.
3.3. Effective non-reciprocal interactions between cell body and flagella
The three-sphere model (§ 3.2) enables us to separately evaluate the effects of the F-to-B force and the dynamic and static components of the B-to-F force. We consider the time instant when the three spheres form an isosceles right triangle. The left flagellar sphere moves in the direction
$\hat {\boldsymbol{t}}_1 = (\cos \alpha , \sin \alpha , 0)$
, where we set the angle relative to the
$x$
axis as
$\alpha = 5\pi /4$
. Reversing the direction of motion by setting
$\alpha = \pi /4$
corresponds to switching from the power stroke to the recovery stroke, and reverses the sign of
$U_{{b}}$
. Figure 4(b) shows that the F-to-B force alone can significantly reduce
$U_{{b}}$
compared with that without any HIs. This implies that the reduction in swimming speed during the recovery stroke (figure 3) results from the impedance of the cell body by the flagella-induced flow. The dynamic component of the B-to-F force increases
$U_{{b}}$
, while the static component has the opposite effect but decays more rapidly with the distance from the body surface
$h$
. When
$h \gtrsim r_{{b}}$
, which is approximately the average separation between the flagella and cell body (figure 1
a), the static component has a negligible effect on
$U_{{b}}$
. Therefore, the slight increase in
$U_{{b}}$
during the power stroke (figure 3
a) may arise from the dominance of the dynamic component of the B-to-F force.
(a) Schematic of the three-sphere model. (b) Effect of different components of body–flagella HIs on the body speed
$U_{{b}}$
in the three-sphere model. The radius of the flagellar sphere is
$r_{{\kern-1.5pt f}} = 0.1\ r_{{b}}$
. (c) Average of the normalised drag forces
$\langle \boldsymbol{F}_{i\to {b}} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{{b}}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
and
$\langle \boldsymbol{F}_{\!{b}\to i} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{\kern-1pt i}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
(see text) over the power and recovery strokes.
To verify the predictions from the three-sphere model, we evaluate how the cell body and flagella contribute to the drag forces on each other in the flagellated swimming model (§ 2). We approximate the spheroidal cell body as a sphere of radius
$r_{{b}} = (a+b)/2$
, and compute the F-to-B force
$\boldsymbol{F}_{i\to {b}}(t)$
by integrating (3.4) along the arc length
$s$
, with
$\boldsymbol{F}_{\kern-1pt i}$
replaced by the flagellar force density
$\boldsymbol{f}_{\kern-1pt i}(s)$
. The dynamic component of the B-to-F force is computed as
\begin{align} \boldsymbol{F}_{{b} \to i}(t) = \int _{C_i} [\zeta _\parallel \boldsymbol{p}_i\boldsymbol{p}_i + \zeta _\perp (\unicode{x1D644} - \boldsymbol{p}_i\boldsymbol{p}_i)] \boldsymbol{\cdot } \boldsymbol{v}_{{b}}[\boldsymbol{x}_i(s,t), \boldsymbol{x}_{{b}}; \boldsymbol{U}_{\!{b}}, r_{{b}}] \, {\rm d}s, \end{align}
where
$\zeta _\parallel$
and
$\zeta _\perp$
are the parallel and perpendicular drag coefficients of the slender flagella (Cox Reference Cox1970), and
$\boldsymbol{p}_i$
is the unit tangent vector. To compare the relative strength of the effects of
$\boldsymbol{F}_{i\to {b}}$
on the cell body and
$\boldsymbol{F}_{{b} \to i}$
on the flagella, we normalise
$\boldsymbol{F}_{i\to {b}}$
by
$\boldsymbol{F}_{{b}}^0(t) = -6 \pi \mu r_{{b}}\, \boldsymbol{U}_{\!{b}} (t)$
, and
$\boldsymbol{F}_{{b} \to i}$
by
\begin{align} \boldsymbol{F}_{\kern-1pt i}^0 (t) = -\int _{C_i} [\zeta _\parallel \boldsymbol{p}_i\boldsymbol{p}_i + \zeta _\perp (\unicode{x1D644} - \boldsymbol{p}_i\boldsymbol{p}_i)] \boldsymbol{\cdot } [\boldsymbol{U}_{\!{b}}(t) + \boldsymbol{U}_{\!i} (s, t)] \, {\rm d}s. \end{align}
Here,
$\boldsymbol{F}_{{b}}^0$
and
$\boldsymbol{F}_{\kern-1pt i}^0$
are the drag forces on the cell body and flagella, respectively, due to their own translational motions in a quiescent fluid. We normalise the components of
$\boldsymbol{F}_{i\to {b}}$
and
$\boldsymbol{F}_{{b} \to i}$
along the swimming direction
$\hat {\boldsymbol{y}}$
, and compute the average over the power and recovery strokes separately.
Figure 4(c) shows that the body–flagella HIs are effectively non-reciprocal during both the power and recovery strokes. During the power stroke,
$\langle \boldsymbol{F}_{\!{b}\to i} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{\kern-1pt i}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
is larger than
$\langle \boldsymbol{F}_{i\to {b}} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{{b}}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
, and the effect of the B-to-F force is stronger than that of the F-to-B force, leading to a slightly larger forward speed (figure 3
a). During the recovery stroke,
$\langle \boldsymbol{F}_{i\to {b}} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{{b}}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
is larger than
$\langle \boldsymbol{F}_{\!{b}\to i} \boldsymbol{\cdot } \hat {\boldsymbol{y}} / \boldsymbol{F}_{\kern-1pt i}^0 \boldsymbol{\cdot } \hat {\boldsymbol{y}}\rangle$
, and the effect of the F-to-B force dominates, leading to a smaller backward speed. The reversal of non-reciprocity when switching from the power stroke to the recovery stroke (figure 4
c) arises from the time-non-reciprocal deformation of the flagella. During the power stroke, the flagella are extended straight and oriented perpendicular to the swimming direction, resulting in a stronger B-to-F effect. During the recovery stroke, the flagella are coiled, more aligned with the swimming direction, and positioned closer to the cell body, resulting in a stronger F-to-B effect.
4. Conclusion
We modelled the swimming dynamics of C. reinhardtii using numerical simulations based on a boundary element method and experimentally measured flagellar waveform. We find that the swimming performance is significantly enhanced by the body–flagella HIs. As the cell body size increases, an optimal body size for swimming appears due to the competition between enhanced body–flagella HIs and increased viscous drag on the cell body. The measured average body size of C. reinhardtii is close to the computed optimal value that maximises the swimming efficiency. Using the three-sphere model, we demonstrate that the body–flagella HIs are effectively non-reciprocal: the body affects the flagella more strongly than vice versa during the power stroke, while the reverse holds during the recovery stroke. As a result, the forward speed is larger than that without HIs during the power stroke, and the backward speed is smaller during the recovery stroke.
Our results provide a hydrodynamic interpretation for the characteristic body size of microalgae. The existence of an optimal cell body size also suggests that the relation between the swimming performance and cargo size may be non-monotonic for certain swimmers composed of multiple interacting components, such as the three-sphere swimmer, which moves through sequential changes in the linear separations between the spheres (Nasouri, Vilfan & Golestanian Reference Nasouri, Vilfan and Golestanian2019). This optimality could inform laboratory designs of biohybrid microrobots (Zhang et al. Reference Zhang, Li, Chen, Luan, Fang, Zhang and Wang2024).
Non-reciprocity generally emerges in interactions mediated by velocity fields. Another example of non-reciprocal HIs is found in dense cilia arrays (Hickey, Golestanian & Vilfan Reference Hickey, Golestanian and Vilfan2023). Two other geometric features of eukaryotic microorganisms, namely the ratio of flagellar wave amplitude to wavelength and the body aspect ratio, also align with hydrodynamic predictions (Lisicki, Rodrigues & Lauga Reference Lisicki, Rodrigues and Lauga2024). Finally, beyond hydrodynamic effects, various factors can shape cell morphology. In particular, biological constraints naturally impose a lower limit on body size, as essential organelles such as the nucleus and chloroplasts require a minimum volume to function.
Acknowledgements
We thank J. Zhang and Y. Man for helpful conversations.
Funding
S.H. acknowledges support from the National Natural Science Foundation of China (NSFC grant no. 12504234) and the Fundamental Research Funds for the Central Universities. D.W. acknowledges support from NSFC grant no. 12574238. Z.L. acknowledges support from NSFC grant no. 12202438.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All the data supporting this work and the simulation codes that reproduce the dynamics of the rigid scallop model using both the hybrid BEM method and the non-local slender body theory are available at https://github.com/shiyuanhu/microalgal_swimming.
Appendix A
As a test problem for our implementation of the hybrid boundary element method (BEM) and regularised Stokeslet method, we consider a rigid scallop consisting of two slender rigid filaments hinged at a common end. This model does not swim on average, but oscillates. The orientation of the upper filament is prescribed as
$\theta (t) = \theta _0[\sin (\omega t) + 1]$
, where
$\theta _0$
is the amplitude, and
$\omega$
is the frequency. The lower filament is mirror-symmetric.
Comparison of the numerical results for a rigid scallop model, obtained from a non-local slender body theory and the hybrid BEM and regularised Stokeslet method used in this work. The two filaments are separated by a small distance at the hinge point to avoid overlapping. The number of line elements is
$N=100$
, and
$m=6$
.
From a non-local slender body theory (Johnson Reference Johnson1980), the velocity of the upper filament
$\partial \boldsymbol{x}_1/\partial t$
is related to the force densities
$\boldsymbol{f}_{\kern-1pt i}$
(for
$i=1,2$
) through
\begin{align} 8\pi \mu\! \left [\frac {\partial \boldsymbol{x}_1(s,t)}{\partial t} - \boldsymbol{v}_{2\to 1}(s,t)\right ] = \left [c (\unicode{x1D644} + \boldsymbol{p}_1\boldsymbol{p}_1) + 2(\unicode{x1D644} - \boldsymbol{p}_1\boldsymbol{p}_1)\right ]\boldsymbol{\cdot } \boldsymbol{f}_{\!1}(s) + \unicode{x1D646}[\boldsymbol{f}_{\!1}](s), \end{align}
where the slenderness parameter is
$c = |{\ln (\delta ^2 e)}|$
, with
$\delta$
the aspect ratio, the unit tangent vector is
$\boldsymbol{p} = (\cos \theta , \sin \theta )$
, and the integral operator
$\unicode{x1D646}$
captures the non-local effects within the filament. The HIs between the two filaments are captured by
\begin{align} \boldsymbol{v}_{2\to 1}(s) = \frac {1}{8\pi \mu }\int _0^L \frac {\unicode{x1D644} + \hat {\boldsymbol{R}}(s,s')\,\hat {\boldsymbol{R}}(s,s')}{R(s,s')} \boldsymbol{\cdot } \boldsymbol{f}_{\!2}(s')\,{\rm d}s', \end{align}
where
$\boldsymbol{R}(s,s') = \boldsymbol{x}_1(s){-} \boldsymbol{x}_2(s')$
. The filament velocity can be decoupled into a translation of the hinged point with velocity
$U(t)\,\hat {\boldsymbol{x}}$
and a rotation with angular velocity
$\dot {\theta }$
.
We discretise the filament centreline using Chebyshev points. The operator
$\unicode{x1D646}[\boldsymbol{f}](s)$
is first regularised (Tornberg & Shelley Reference Tornberg and Shelley2004) and then evaluated via Clenshaw–Curtis quadrature with spectral accuracy (Trefethen Reference Trefethen2000). Equations (A1) and (A2) form a dense linear system for the unknowns
$\boldsymbol{f}_{\kern-1pt i}$
and
$U(t)$
. Figure 5 shows that the result obtained from the hybrid BEM method using a line distribution of regularised Stokeslets agrees closely with the non-local slender body theory.
Appendix B
We quantify the geometric properties of waveforms used in our simulations by separating the tangent angle
$\theta (s,t)$
into static and dynamic components (Geyer et al. Reference Geyer, Sartori, Friedrich, Jülicher and Howard2016). The static component is defined as the time average of the tangent angle over a beating period,
$\bar {\theta }(s) = (1/T)\int _0^T \theta (s,t)\,{\rm d}t$
, and the dynamic component is defined as
$\tilde {\theta }(s,t) = \theta (s,t){-}\bar {\theta }(s)$
. The static component
$\bar {\theta }(s)$
is approximately linear in
$s$
, which corresponds to a circular shape. This static shape is characterised by its curvature
$C = |\bar {\theta }(L) - \bar {\theta }(0)|/L$
. The dynamic component
$\tilde {\theta }(s,t)$
is close to a travelling wave with amplitude
$A$
and wavelength
$\lambda$
. These geometric parameters are summarised in table 2. Although the four waveforms differ significantly in curvature
$C$
, length
$L$
, and wavelength
$\lambda$
, the computed swimming speed
$\langle U_{{b}} \rangle$
and efficiency
$\eta$
consistently display pronounced optima as a function of the body size
$b$
when HIs are incorporated (figure 6).
Variations in experimental waveforms explored in simulations. Dataset 1 is used in the main text.
Effect of cell body size
$b/L$
on the swimming speed
$\langle U_{{b}}\rangle /U_0$
and efficiency
$\eta$
, obtained from three additional waveforms measured in experiments (table 2): (a) dataset 2, (b) dataset 3, (c) dataset 4. The vertical dashed lines mark the values of
$b/L$
of each sample cell.
