Skip to main content Accessibility help


  • Access
  • Open access



MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Swimming mediated by ciliary beating: comparison with a squirmer model
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Swimming mediated by ciliary beating: comparison with a squirmer model
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Swimming mediated by ciliary beating: comparison with a squirmer model
        Available formats
Export citation


The squirmer model of Lighthill and Blake has been widely used to analyse swimming ciliates. However, real ciliates are covered by hair-like organelles, called cilia; the differences between the squirmer model and real ciliates remain unclear. Here, we developed a ciliate model incorporating the distinct ciliary apparatus, and analysed motion using a boundary element–slender-body coupling method. This methodology allows us to accurately calculate hydrodynamic interactions between cilia and the cell body under free-swimming conditions. Results showed that an antiplectic metachronal wave was optimal in the swimming speed with various cell-body aspect ratios, which is consistent with former theoretical studies. Exploiting oblique wave propagation, we reproduced a helical trajectory, like Paramecium, although the cell body was spherical. We confirmed that the swimming velocity of model ciliates was well represented by the squirmer model. However, squirmer modelling outside the envelope failed to estimate the energy costs of swimming; over 90 % of energy was dissipated inside the ciliary envelope. The optimal swimming efficiency was given by the antiplectic wave; the value was 6.7 times larger than in-phase beating. Our findings provide a fundamental basis for modelling swimming micro-organisms.

1 Introduction

Swimming micro-organisms are ubiquitous, including oceanic micro-algae and human gut bacteria. To clarify swimming dynamics, several fluid mechanical models have been created to deal with the various geometries and swimming modes of natural micro-organisms. Lighthill (1952) introduced a simplified ciliate model, the so-called ‘squirmer’. The model was generalised by Blake (1971) and has been used to analyse various ciliate species. Brennen (1974) outlined a fluid mechanical model of self-propelled micro-organisms by introducing an oscillating boundary layer theory for ciliary propulsion. Keller & Wu (1977) built on the model by extending the squirmer to be prolate spheroidal in shape, and estimated the effect of shape on energy expenditure. Michelin & Lauga (2010) showed that the minimum energy dissipation at a given swimming speed was afforded by a neutral swimmer, thus neither a puller nor a pusher. Magar, Goto & Pedley (2003) investigated microbial nutrient uptake from the environment. Squirming enhanced uptake compared to that of a rigid sphere moving at the same speed. Ishikawa et al. (2016) extended the analysis to squirmer suspensions; nutrient uptake increased in a quadratic manner to volume fractions of up to 30 %. The squirmer model has also been used to study bioconvection in suspensions of upswimming cells, which generate gravitational instability (Kessler 1986; Ishikawa, Locsei & Pedley 2008; Evans et al. 2011), and to explore the coherent structures of non-gravitational active swimmers (Simha & Ramaswamy 2002; Saintillan & Shelley 2008). Ishikawa & Pedley (2007) computed the translational diffusion tensors of squirmers in semi-dilute suspensions; these can also be derived analytically at the dilution limits. Pedley, Brumley & Goldstein (2016) extended the Lighthill–Blake squirmer to incorporate an azimuthal swirl mode. The details can be found in a recent review by Pedley (2016).

The Lighthill–Blake squirmer is a sphere with a deformable surface reflecting periodic ciliary beating. The squirming surface is defined by the mean height of ciliary tips that exhibit low-amplitude oscillations. The deformable/stretchable spherical surface is considered to be an envelope marking the boundary of beating ciliary tips. When squirming velocity on the deformable surface is considered in both the radial and tangential directions, the body swims with a non-zero mean velocity. The squirmer model describes fluid motions outside the ciliary envelope; flow inside the envelope is ignored. In reality, however, fluid flow exists even within the envelope. Only a few studies (e.g. Keller & Wu 1977) have compared the energy expenditure within the ciliary layer to that outside the ciliary layer. Thus, flow generated by individual ciliary motions inside the envelope requires further attention.

Cilia-driven flow has been investigated by several groups (Blake 1972, 1974; Niedermayer, Eckhardt & Lenz 2008; Brumley et al. 2012; Eloy & Lauga 2012; Elgeti & Gompper 2013; Nasouri & Elfring 2016). In cilia-driven flows, metachronal wave propagation plays a role, and four metachronal wave patterns have been recognised (Sleigh 1962): a symplectic metachronal wave (the beat direction and that of wave transmission coincide); an antiplectic metachronal wave (the two directions are opposed); a dexioplectic metachronal wave (the effective beat moves to the right with respect to the direction of wave transmission); and a laeoplectic metachronal wave (the effective beat moves to the left with respect to the direction of wave transmission). These different wave patterns have been observed in various species. For example, symplectic waves are characteristic of Opalina (Sleigh 1962) and dexioplectic wave propagation has been observed in Paramecium swimming under normal conditions (Machemer 1972). Blake (1972) analysed the force, bending moment and rate of work exerted on a cilium with various types of metachronism. He reported that in antiplectic metachronism neighbouring cilia spread out during the effective stroke. Thus, the force exerted on a cilium and the rate of work were large during the effective stroke. Eloy & Lauga (2012) derived optimal ciliary beat waveforms with varying sperm number. Elgeti & Gompper (2013) numerically studied how ciliary geometrical parameters, such as ciliary spacing and beat direction, affected the properties of metachronal waves. The waves exhibited flow-induced emergence; development of an antiplectic metachronal wave increased the propulsion velocity by more than 3–10-fold compared to in-phase beating. Several studies have also found an antiplectic metachronal wave to be optimal for efficiency (Blake 1972; Michelin & Lauga 2010; Osterman & Vilfan 2011). To clarify how metachronal waves emerge, several self-sustained oscillator models have been developed. Brumley et al. (2012) presented such a model sustained by a linear elastic spring. It was concluded that orbit compliance mediated by spring elasticity facilitated fast robust synchronisation of two oscillators.

Although earlier studies on squirmers and cilia-driven flow yielded valuable insights into the dynamics of swimming ciliates and the propulsion velocity of metachronal waves, their relationship remains unclear, especially the difference between swimming mediated by surface squirming and ciliary motions. Here, we develop a three-dimensional ciliate model incorporating individual ciliary motions on cell surfaces. We analyse the flow field using a boundary element–slender-body coupling method. This methodology allows us to accurately calculate hydrodynamic interactions between cilia and the cell body under free-swimming conditions. The methodology is detailed in § 2. Ciliate swimming behaviours in terms of beat phase, cilia number and body aspect ratio are discussed in § 3. We use oblique metachronal waves to develop three-dimensional helical trajectories. In § 4, we compare our model ciliates with Lighthill–Blake squirmers, focusing on energy dissipation inside/outside the ciliary envelope. Section 5 provides conclusions.

2 Governing equations and numerical methods

2.1 Boundary integral equation with slender-body theory

Consider a ciliate immersed in an infinite Newtonian liquid of density $\unicode[STIX]{x1D70C}$ and viscosity $\unicode[STIX]{x1D707}$ ; the ciliate is propelled via individual ciliary motion. The inertial effect of fluid flow is negligible; the Reynolds numbers scaled by ciliary motion and swimming are markedly lower than 1 ( $Re\ll 1$ ). Thus, fluid flow around the ciliate is governed by the Stokes equation. The cell body is modelled as a rigid spheroid from which cilia emerge. As cilia are slender, slender-body theory (Tornberg & Shelly 2004) is applicable when analysing ciliary motion. We parametrise the ciliary centreline using arclength $s\in [0,L]$ , where $L$ is the ciliary length. The flow field at point $\boldsymbol{x}$ located on the $i$ th cilium, $\boldsymbol{x}\in s_{i}$ , is given by Pozrikidis (1992) and Tornberg & Shelly (2004) as

(2.1) $$\begin{eqnarray}\displaystyle \boldsymbol{v}(\boldsymbol{x}) & = & \displaystyle -\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\int _{cell}\boldsymbol{J}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{q}(\boldsymbol{y})\,\text{d}A(\boldsymbol{y})-\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\unicode[STIX]{x1D726}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{x})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\int _{cilia}[\boldsymbol{J}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{y})+\boldsymbol{K}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{x})]\,\text{d}s_{i}(\boldsymbol{y})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\mathop{\sum }_{j\neq i}^{N}\int _{cilia}[\boldsymbol{J}(\boldsymbol{x},\boldsymbol{y})+\boldsymbol{W}(\boldsymbol{x},\boldsymbol{y})]\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{y})\,\text{d}s_{j}(\boldsymbol{y}),\end{eqnarray}$$

where $\boldsymbol{q}$ is the viscous traction on the cell body, $\boldsymbol{f}$ is the force density per unit length and $N$ is the total number of cilia. The first integral on the right operates over the entire spheroidal cell surface, and the second and third integrals operate along the central ciliary lines. The Green’s function $\boldsymbol{J}$ is given by

(2.2) $$\begin{eqnarray}J_{ij}(\boldsymbol{x},\boldsymbol{y})=\frac{\unicode[STIX]{x1D6FF}_{ij}}{r}+\frac{r_{i}r_{j}}{r^{3}},\end{eqnarray}$$

where $r=|\boldsymbol{r}|$ and $\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{y}$ . In (2.1) $\unicode[STIX]{x1D726}$ and $\boldsymbol{K}$ are the local operators of the slender-body theory, which are given by (Tornberg & Shelly 2004)

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{ij}(\boldsymbol{x})=c[\unicode[STIX]{x1D6FF}_{ij}+t_{i}(\boldsymbol{x})t_{j}(\boldsymbol{x})]+2[\unicode[STIX]{x1D6FF}_{ij}-t_{i}(\boldsymbol{x})t_{j}(\boldsymbol{x})]\end{eqnarray}$$


(2.4) $$\begin{eqnarray}K_{ij}(\boldsymbol{x},\boldsymbol{y})=-\frac{\unicode[STIX]{x1D6FF}_{ij}+t_{i}(\boldsymbol{x})t_{j}(\boldsymbol{x})}{|s(\boldsymbol{x})-s(\boldsymbol{y})|},\end{eqnarray}$$

where $c=-\ln (\unicode[STIX]{x1D700}^{2}e)$ , $\boldsymbol{t}$ is the unit tangential vector to the centreline of each cilium and $\unicode[STIX]{x1D700}=a_{cilia}/L$ , with $a_{cilia}$ the ciliary radius. In terms of the practical ciliary radius and length ratio (Sleigh 1962), the slenderness value $\unicode[STIX]{x1D700}$ is set to $\unicode[STIX]{x1D700}=0.01$ throughout the present study. The slender-body kernel $\boldsymbol{W}$ is defined by

(2.5) $$\begin{eqnarray}W_{ij}(\boldsymbol{x},\boldsymbol{y})=\frac{(\unicode[STIX]{x1D700}L)^{2}}{2}\left(\frac{\unicode[STIX]{x1D6FF}_{ij}}{r^{3}}-3\frac{r_{i}r_{j}}{r^{5}}\right).\end{eqnarray}$$

When the observation point $\boldsymbol{x}$ is not on a cilium, $\boldsymbol{x}\notin s$ , the velocity is given by

(2.6) $$\begin{eqnarray}\displaystyle \boldsymbol{v}(\boldsymbol{x}) & = & \displaystyle -\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\int _{cell}\boldsymbol{J}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{q}(\boldsymbol{y})\,\text{d}A(\boldsymbol{y})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\mathop{\sum }_{j=1}^{N}\int _{cilia}[\boldsymbol{J}(\boldsymbol{x},\boldsymbol{y})+\boldsymbol{W}(\boldsymbol{x},\boldsymbol{y})]\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{y})\,\text{d}s_{j}(\boldsymbol{y}).\end{eqnarray}$$

Though the asymptotic accuracy of kernels $\unicode[STIX]{x1D726}$ and $\boldsymbol{K}$ is $O(\unicode[STIX]{x1D700}^{2}\ln \unicode[STIX]{x1D700})$ , kernel $\boldsymbol{W}$ is accurate only to the limit $O(\unicode[STIX]{x1D700})$ (Tornberg & Shelly 2004). Equations (2.1) and (2.6) are therefore accurate to the limit $O(\unicode[STIX]{x1D700})$ .

Figure 1. Problem setting. (a) Schematics of the body frames and the angles $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ . (b) Beat pattern of each cilium described by the parameters in table 1.

2.2 Ciliary motions

We first define an orthonormal frame for the cell body $\boldsymbol{e}_{i}$ with origin $\boldsymbol{x}_{c}$ , where $\boldsymbol{x}_{c}$ is the body’s centre of mass; $\boldsymbol{e}_{1}$ then reflects swimming orientation. To efficiently model ciliary motion on the cell surface, local vectors with an orthonormal basis $\boldsymbol{g}_{i}$ are defined as follows. A material point $\boldsymbol{x}_{b}$ , located at the base of a cilium on the cell surface, serves as the origin of the orthonormal body frame $\boldsymbol{g}_{i}$ (cf. figure 1 a). The basis vectors $\boldsymbol{g}_{1}$ and $\boldsymbol{g}_{2}$ are defined as $\boldsymbol{g}_{1}(\boldsymbol{x}_{b})=\boldsymbol{b}(\boldsymbol{x}_{b})\wedge \boldsymbol{n}(\boldsymbol{x}_{b})/|\boldsymbol{b}(\boldsymbol{x}_{b})\wedge \boldsymbol{n}(\boldsymbol{x}_{b})|$ and $\boldsymbol{g}_{2}(\boldsymbol{x}_{b})=\boldsymbol{n}(\boldsymbol{x}_{b})$ , where $\boldsymbol{b}=\boldsymbol{e}_{1}\wedge \boldsymbol{n}$ and $\boldsymbol{n}$ is the outward unit normal vector, respectively. The time-dependent profile of each ciliary motion is derived using the following mathematical formula of Fulford & Blake (1986):

(2.7) $$\begin{eqnarray}\boldsymbol{x}^{cilia}(\boldsymbol{x}_{b},s,t)=\unicode[STIX]{x1D709}^{1}(s,t)\boldsymbol{g}_{1}(\boldsymbol{x}_{b})+\unicode[STIX]{x1D709}^{2}(s,t)\boldsymbol{g}_{2}(\boldsymbol{x}_{b}),\end{eqnarray}$$


(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D709}^{i}(s,t)=\frac{1}{2}\unicode[STIX]{x1D6FC}_{0}^{i}(s)+\mathop{\sum }_{n=1}^{N_{0}}\unicode[STIX]{x1D6FC}_{n}^{i}(s)\cos n\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D6FD}_{n}^{i}(s)\sin n\unicode[STIX]{x1D714}t,\end{eqnarray}$$

with $\unicode[STIX]{x1D714}$ the angular beat frequency. The Fourier coefficients $\unicode[STIX]{x1D6FC}_{n}^{i}$ and $\unicode[STIX]{x1D6FD}_{n}^{i}$ are given by

(2.9a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{n}^{i}(s)=\mathop{\sum }_{m=1}^{M_{0}}A_{mn}^{i}s^{m},\quad \unicode[STIX]{x1D6FD}_{n}^{i}(s)=\mathop{\sum }_{m=1}^{M_{0}}B_{mn}^{i}s^{m}.\end{eqnarray}$$

The coefficients $A_{mn}^{i}$ and $B_{mn}^{i}$ are summarised in table 1. In this study, the wavenumbers $N_{0}$ and $M_{0}$ are set to 3. The beat pattern described by the parameters in table 1 is shown in figure 1(b).

Table 1. Fourier coefficients for ciliary beating, which are taken from Fulford & Blake (1986).

2.3 Boundary element method

In order to simulate free swimming of the ciliate model, force-free and torque-free conditions are taken into account:

(2.10) $$\begin{eqnarray}\int _{cell}\boldsymbol{q}\,\text{d}A+\mathop{\sum }_{i=1}^{N}\int _{cilia}\boldsymbol{f}\,\text{d}s_{i}=0\end{eqnarray}$$


(2.11) $$\begin{eqnarray}\int _{cell}\boldsymbol{q}\wedge \hat{\boldsymbol{r}}\,\text{d}A+\mathop{\sum }_{i=1}^{N}\int _{cilia}\boldsymbol{f}\wedge \hat{\boldsymbol{r}}\,\text{d}s_{i}=0,\end{eqnarray}$$

where $\hat{\boldsymbol{r}}=\boldsymbol{x}-\boldsymbol{x}_{c}$ , and $\boldsymbol{x}_{c}$ is the centre of mass of the cell body.

Assuming that the cell body shows a rigid motion, the velocity on the spheroidal cell body surface and on the cilia can be expressed by

(2.12) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\boldsymbol{v}(\boldsymbol{x})=\boldsymbol{V}+\unicode[STIX]{x1D734}\wedge \hat{\boldsymbol{r}}(\boldsymbol{x}), & \quad \boldsymbol{x}\in \text{cell body},\\ \boldsymbol{v}(\boldsymbol{x})=\boldsymbol{V}+\unicode[STIX]{x1D734}\wedge \hat{\boldsymbol{r}}(\boldsymbol{x})+\boldsymbol{v}^{cilia}(\boldsymbol{x}), & \quad \boldsymbol{x}\in \text{cilia},\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{V}$ is the translational and $\unicode[STIX]{x1D734}$ the angular velocity, and $\boldsymbol{v}^{cilia}$ is the ciliary velocity with respect to the body frame $\boldsymbol{e}_{i}$ (i.e.  $\boldsymbol{v}^{cilia}=\unicode[STIX]{x2202}\boldsymbol{x}^{cilia}/\unicode[STIX]{x2202}t$ ). We then solve the resistance problem with respect to the unknowns $\boldsymbol{V}$ and $\unicode[STIX]{x1D734}$ , and derive the viscous tractions $\boldsymbol{q}$ and $\boldsymbol{f}$ .

The cell body is modelled as a rigid spheroid, and the body surface is discretised into 5120 triangular mesh elements with 2562 nodal points. Each cilium is discretised into 16 nodes that are interpolated using the centripetal Catmull–Rom spline method (Catmull & Rom 1974). All physical quantities are computed at each discretised point. The boundary integrals of (2.1) and (2.6) are computed with the aid of numerical Gaussian integration. When an observation point $\boldsymbol{x}$ is located on the cell body, the following linear algebraic equation can be derived from (2.6):

(2.13) $$\begin{eqnarray}\displaystyle [\boldsymbol{v}^{b}]=[{\mathcal{J}}^{bb}][\boldsymbol{q}]+[{\mathcal{J}}^{bc}][\,\boldsymbol{f}]. & & \displaystyle\end{eqnarray}$$

When $\boldsymbol{x}$ is on the cilia, on the other hand, we have the following equation from (2.1):

(2.14) $$\begin{eqnarray}\displaystyle [\boldsymbol{v}^{c}]=[{\mathcal{J}}^{cb}][\boldsymbol{q}]+[{\mathcal{J}}^{cc}][\,\boldsymbol{f}]. & & \displaystyle\end{eqnarray}$$

The vector sizes of $[\boldsymbol{v}^{b}]$ and $[\boldsymbol{q}]$ are $3N_{b}$ , while $[\boldsymbol{v}^{c}]$ and $[\,\boldsymbol{f}]$ have the size $3N_{c}$ , where $N_{b}$ is the number of nodes on the cell body and $N_{c}$ is the total number of nodes on the cilia. The matrix sizes of $[{\mathcal{J}}^{bb}]$ and $[{\mathcal{J}}^{bc}]$ are $3N_{b}\times 3N_{b}$ and $3N_{b}\times 3N_{c}$ , whereas $[{\mathcal{J}}^{cb}]$ and $[{\mathcal{J}}^{cc}]$ are $3N_{c}\times 3N_{b}$ and $3N_{c}\times 3N_{c}$ , respectively. In a similar manner, discretised forms of the force–torque conditions, (2.10) and (2.11), can be written as

(2.15) $$\begin{eqnarray}\displaystyle [{\mathcal{F}}^{b}][\boldsymbol{q}]+[{\mathcal{F}}^{c}][\,\boldsymbol{f}]=[\mathbf{0}] & & \displaystyle\end{eqnarray}$$


(2.16) $$\begin{eqnarray}\displaystyle [{\mathcal{T}}^{b}][\boldsymbol{q}]+[{\mathcal{T}}^{c}][\,\boldsymbol{f}]=[\mathbf{0}]. & & \displaystyle\end{eqnarray}$$

Considering the boundary condition of (2.12), the system can be expanded to

(2.17) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}cccc@{}}{\mathcal{J}}^{bb} & {\mathcal{J}}^{bc} & {\mathcal{V}}^{b} & {\mathcal{A}}^{b}\\ {\mathcal{J}}^{cb} & {\mathcal{J}}^{cc} & {\mathcal{V}}^{c} & {\mathcal{A}}^{c}\\ {\mathcal{F}}^{b} & {\mathcal{F}}^{c} & \mathbf{0} & \mathbf{0}\\ {\mathcal{T}}^{b} & {\mathcal{T}}^{c} & \mathbf{0} & \mathbf{0}\end{array}\right]\left[\begin{array}{@{}c@{}}\boldsymbol{q}\\ \boldsymbol{f}\\ \boldsymbol{V}\\ \unicode[STIX]{x1D734}\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathbf{0}\\ \boldsymbol{v}^{cilia}\\ \mathbf{0}\\ \boldsymbol{ 0}\end{array}\right]. & & \displaystyle\end{eqnarray}$$

The matrix components ${\mathcal{V}}^{b}$ and ${\mathcal{A}}^{b}$ are of sizes $3N_{b}\times 3$ , whereas ${\mathcal{V}}^{c}$ and ${\mathcal{A}}^{c}$ are both $3N_{c}\times 3$ . The dense matrix (2.17) is solved using the lower–upper (LU) factorisation technique. Given the translational and angular velocities, all material points are updated using the second-order Runge–Kutta method. Validations of the numerical method are shown in appendix A.

2.4 Parameter setting

In a Stokes flow regime, the fluid viscosity is simply a multiplier of both force and traction. The viscosity $\unicode[STIX]{x1D707}$ can be assumed to be unity, without loss of generality. We further assume that all cilia exhibited identical beat frequencies, and beat periodically during computation. To express phase differences among ciliary beats on the cell surface (i.e. metachronal waves), the initial ciliary beat phase $\unicode[STIX]{x1D713}^{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ was defined as

(2.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}^{0}/2\unicode[STIX]{x03C0}=\sin (k\unicode[STIX]{x1D703}/2)+\sin (\unicode[STIX]{x1D708}\unicode[STIX]{x1D719}/4), & & \displaystyle\end{eqnarray}$$

where $k$ and $\unicode[STIX]{x1D708}$ are the wavenumbers in the $\unicode[STIX]{x1D703}$ - and $\unicode[STIX]{x1D719}$ -directions, respectively (cf. figure 1 a). The angle $\unicode[STIX]{x1D703}=[0,\unicode[STIX]{x03C0}]$ is that between the orientation vector $\boldsymbol{e}_{1}$ and $\hat{\boldsymbol{r}}_{b}$ ; and $\unicode[STIX]{x1D703}=\cos ^{-1}(\boldsymbol{e}_{1}\boldsymbol{\cdot }\hat{\boldsymbol{r}}_{b})$ , where $\hat{\boldsymbol{r}}_{b}=(\boldsymbol{x}_{b}-\boldsymbol{x}_{c})/|\boldsymbol{x}_{b}-\boldsymbol{x}_{c}|$ . The angle $\unicode[STIX]{x1D719}=[-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$ is defined as $\unicode[STIX]{x1D719}=\pm \cos ^{-1}(\boldsymbol{e}_{2}\boldsymbol{\cdot }\hat{\boldsymbol{r}}_{b})$ ; the sign is determined by the sign of $\boldsymbol{e}_{3}\boldsymbol{\cdot }\hat{\boldsymbol{r}}_{b}$ .

In the following sections, the wavenumber $k$ varies from $-2.0$ to 2.0, and $\unicode[STIX]{x1D708}$ is fixed at 0 except in § 3.5. If $k$ is positive, an antiplectic metachronal wave is in play; a symplectic metachronal wave is triggered by a negative $k$ . When $k$ is set to zero, all cilia beat in phase. In § 3.5, we consider oblique wave propagation, i.e. dexioplectic/laeoplectic metachronal waves; we set $\unicode[STIX]{x1D708}\neq 0$ .

The cell body is modelled by a rigid spheroid of major radius $a_{0}$ and minor axis $b_{0}$ . To explore the effects of aspect ratio, we considered two cell shapes. First, the cell shape was controlled to hold the volume constant regardless of the aspect ratio $\unicode[STIX]{x1D6FC}_{v}$ . Next, the major radius was fixed in various aspect ratio, represented as $\unicode[STIX]{x1D6FC}_{l}$ . For all computations, the time interval $\unicode[STIX]{x0394}t$ was set to $\unicode[STIX]{x0394}t/T=0.01$ , where $T$ is the ciliary beat period.

3 Swimming motions induced by ciliary beating

3.1 Spherical ciliates featuring in-phase ciliary beating

We first investigated the swimming of spherical ciliates featuring in-phase ciliary beating. The cell radius $a_{0}$ was set to $a_{0}/L=3.0$ , similar to the value exhibited by freshwater ciliates Tetrahymena. The ciliary number was set to $N=160$ . To reproduce in-phase ciliary beating, the wavenumbers $k$ and $\unicode[STIX]{x1D708}$ in (2.18) were set to zero. The orientation vector was initially set to $\boldsymbol{e}_{1}=(1.0,0,0)$ ; the ciliate swam in the $x$ -direction.

Ciliate configurations during a single ciliary beating phase are shown in figure 2 (see also the supplementary movie). All cilia exhibited periodic effective and recovery strokes; the latter strokes occurred at $0.25\leqslant t/T\leqslant 0.9$ . When the strokes were effective, the ciliate swam forwards, but the direction of motion was reversed during recovery strokes. Thus, the cell was gradually propelled in the $x$ -direction via periodic back-and-forth movement.

To explore swimming velocity in detail, time changes of swimming velocity over a single period are shown in figure 3(a). In the first quarter period $0\leqslant t/T\leqslant 0.25$ , ciliary beating was effective, and swimming velocity was maximal at $t/T\sim 0.1$ . The maximum swimming velocity was approximately $V_{x}T/L\simeq 3.0$ , equivalent to a half-body-length per period. During the recovery phase, the velocity was negative, and the ciliate reversed its progress (thus in the minus $x$ -direction). Although back-and-forth movement was in play, the total displacement $u_{x}=\int _{0}^{T}V_{x}\,\text{d}t$ was positive; the ciliate moved forwards.

Figure 2. Swimming of the spherical ciliate during one period ( $N=160$ , $a_{0}/L=3.0$ and in-phase beating). The colour on the cilia indicates the propulsion force of individual cilia calculated by $\int f_{x}\,\text{d}s$ . A movie can be seen in the supplemental material available at

Figure 3. (a) Swimming velocity of the ciliate as a function of time. (b) Propulsion force of individual cilia in three regions calculated by $\sum _{i}\int f_{x}\,\text{d}s_{i}$ . Regions A, B and C are posterior, middle and anterior regions, as shown in the inset.

In figure 2, all cilia beat in-phase, but the contributions of individual cilia to swimming velocity depended on ciliary position. The contour colours of the ciliary surface indicate propulsion afforded by ciliary motion, i.e. the $\int f_{x}\,\text{d}s$ values. Red suggests that the cell is pushed forwards by ciliary motion, and blue that the cell is forced back by such motion. For example, at $t/T=0.0$ , the anterior cilia are blue, suggesting that locomotion was compromised. Both posterior and central cilia can trigger locomotion. In figure 3(b), the propulsion forces generated in three regions are shown. Region A represents the posterior area (encompassing one-third of the body from the posterior end). Regions B and C are the middle and anterior regions (each encompassing one-third of the body length, as shown in the inset). In figure 3(b), it is apparent that the locomotion curves of different ciliary positions vary, and that the locomotion peak shifts from region A to C over time. The total force curve is similar to the swimming velocity curve; the total ciliate drag is not affected by ciliary configuration.

3.2 Effect of metachronal wave

Elgeti & Gompper (2013) investigated the emergence of metachronal waves in ciliary arrays active on a flat plane, and concluded that such waves increased propulsion velocity. Thus, even when the cell body is spherical, metachronal waves may enhance swimming velocity. We thus investigated the effects of such waves on cellular locomotion. Here, the aspect ratio was set to 1.0, the ciliary radius to $a_{0}/L=3.0$ , and the ciliary number to $N=160$ .

Figure 4. Swimming of the ciliate with different initial phase difference $k$ . By setting positive $k$ , the ciliate shows antiplectic metachronal waves, while negative $k$ represents symplectic metachronal wave. The other parameters are the same as in figure 2.

We first examined the effects of an antiplectic metachronal wave. The wavenumbers $k$ were set to $k=0,~0.5,~1.0$ and 2.0; and $\unicode[STIX]{x1D708}$ to zero. The ciliate configurations at $k=1.0$ are shown in figure 4. The effective beat is propagated from posterior to anterior over time. The temporal swimming velocities during periods with various wavenumbers $k$ are shown in figure 5(a). When $k=0$ (i.e. during in-phase beating), the swimming velocities (both forwards and backwards) are maximised. When an antiplectic metachronal wave forms ( $k>0$ ), the peak is lower than that associated with in-phase beating, but the effective phase time tends to be longer, as shown in figure 5(a). Time-averaged velocities $\bar{V}_{x}$ at various $k$ -values are shown in figure 5(b); $\bar{V}_{x}$ is maximum when $k\simeq 1.0$ , suggesting that an antiplectic metachronal wave with $k=1.0$ is optimal in this cellular configuration.

Figure 5. Swimming velocity with various $k$ . (a) Time change of swimming velocity, and (b) time-averaged swimming speed. When $k>0$ , the ciliate shows antiplectic metachronal waves, whereas it shows symplectic metachronal waves with $k<0$ . Movies of $k=1$ and $k=-1$ are available in the supplementary material.

Figure 6. Swimming of spheroidal ciliates: (a) oblate-type ciliate with aspect ratio $\unicode[STIX]{x1D6FC}_{v}=0.5$ and (b) prolate-type ciliate with $\unicode[STIX]{x1D6FC}_{v}=2.0$ . In both cases, ciliary beats are in-phase.

We also investigated the effect of a symplectic metachronal wave (i.e.  $k<0$ ); the results are also shown in figure 5. Unlike what was seen when an antiplectic wave was in play, swimming velocity was not enhanced by a symplectic metachronal wave. With such a wave, ciliary stroke effectiveness is likely to be compromised by cilia of the posterior side, as shown in figure 4. On the other hand, when an antiplectic wave forms, the distance between the fore and aft cilia engaging in effective strokes becomes larger, and fluid flows freely in the ciliary layer. This agrees well with the findings of earlier studies (Blake 1972; Guo et al. 2014). Thus, an antiplectic metachronal wave enhances the swimming velocity of spherical ciliates.

3.3 Effect of aspect ratio

In previous sections, the cell considered was spherical, like the micro-alga Volvox. In nature, however, micro-organisms assume various forms. For example, Tetrahymena is ellipsoidal with an aspect ratio of approximately 2. We thus investigated the effects of cellular shape by varying the aspect ratio. As explained in § 2.4, we defined two aspect ratios $\unicode[STIX]{x1D6FC}_{v}$ and $\unicode[STIX]{x1D6FC}_{l}$ . The swimming behaviours of an oblate ciliate ( $\unicode[STIX]{x1D6FC}_{v}=0.5$ ) and a prolate ciliate ( $\unicode[STIX]{x1D6FC}_{v}=2.0$ ) are shown in figure 6. For an oblate ciliate, the effective strokes of the anterior and posterior cilia are probably perpendicular to the direction of swimming. Thus, the average swimming velocity is smaller than that of prolate ciliates for both the $\unicode[STIX]{x1D6FC}_{v}$ and $\unicode[STIX]{x1D6FC}_{l}$ cases, as shown in figures 7(a) and 8. When the aspect ratio was varied while keeping the volume constant, the drag coefficient of a rigid spheroid became the minimum when $\unicode[STIX]{x1D6FC}_{v}\sim 2$ , as shown in figure 7(b). Thus, large velocities were observed around $\unicode[STIX]{x1D6FC}_{v}=2$ in the antiplectic mode. The highest velocities appear at the aspect ratio of approximately 1.5 with $k=0.5$ –1.5. Hence, the velocity is not simply a function of aspect ratio, but hydrodynamic interactions between cilia play an important role. For $\unicode[STIX]{x1D6FC}_{l}$ (i.e. the major axis is constant), the cellular volume becomes smaller as the aspect ratio increases and the drag on the cell body thus decreases. Accordingly, the swimming speed increases with $\unicode[STIX]{x1D6FC}_{l}$ in the moderate-aspect-ratio regime, as shown in figure 8.

Figure 7. Effect of aspect ratio $\unicode[STIX]{x1D6FC}_{v}$ so as to keep the volume constant. (a) Average swimming speed of ciliate. For all cases, cilia beat as in-phase and number of cilia $N=160$ . (b) Drag coefficient of a rigid spheroid without cilia as a function of $\unicode[STIX]{x1D6FC}_{v}$ . Here $U$ is the fluid velocity relative to the rigid spheroid.

Figure 8. Effect of aspect ratio $\unicode[STIX]{x1D6FC}_{l}$ so as to keep the major axis constant.

We compared our results with the swimming behaviour of a real ciliate Tetrahymena. In Ishikawa & Kikuchi (2018), the aspect ratio of Tetrahymena was found to be approximately 2.1. The body length was approximately $65~\unicode[STIX]{x03BC}\text{m}$ and the swimming velocity was $440~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ (i.e. 6.8 body lengths per second). In figure 7(b), the average velocity when $\unicode[STIX]{x1D6FC}_{l}=2.0$ and $k=1.0$ was approximately 0.35, corresponding to 0.058 body lengths per beat. The beat frequency of Tetrahymena was reported to be approximately 30 Hz (Ishikawa & Kikuchi 2018). Our ciliate model yields a swimming velocity of 1.74 body lengths per second, thus approximately 3.9-fold less than the actual swimming velocity of Tetrahymena. The real Tetrahymena has more than 1000 cilia, whereas our model features only 160 cilia. We thus explored the effect of cilia numbers on swimming velocity.

Figure 9. Swimming velocity with different number of cilia and metachronal waves: (a) in-phase, (b) antiplectic and (c) symplectic waves. In panels (df), the velocities are normalised by each maximum swimming velocity. For all the cases, we set the aspect ratio equal to 1.0.

3.4 Effect of cilia numbers

Here, we investigate the effect of cilia numbers. The aspect ratio was again set to 1.0 and the radius to $a_{0}/L=3.0$ . The temporal swimming velocities for $N=10,~40,~160$ and 320 are shown in figure 9. The wavenumber $k$ was set to $k=-1,~0$ and 1, and $\unicode[STIX]{x1D708}$ was set to zero. In the case of $k=0$ , i.e. in-phase beating, the cell showed larger back-and-forth movements as $N$ increases; while in the antiplectic waves, we see a longer period of propelled motions. The maximum instantaneous swimming velocity attains $V_{x}T/L\simeq 3.7$ when $N=320$ at in-phase beating, and was $V_{x}T/L\simeq 3.0$ when $N=160$ . Thus, the maximum swimming velocity was not simply proportional to the number of cilia. However, when the velocity was normalised by the maximum velocity, all curves almost converged onto one curve regardless of the number of cilia. Thus, the shape of velocity over time is independent of $N$ , but it depends on the metachronal waves. We next explored average ciliate velocity as a function of the number of cilia $N$ (cf. figure 10). The swimming velocity increased monotonically with $N$ for all $k$ , although the slope decreased as $N$ increased. Again, the wavenumber $k=1.0$ was associated with the maximum velocity regardless of the $N$ -value. Thus, the optimal wavenumber is not greatly affected by cilia numbers.

Figure 10. Average swimming velocity with different number of cilia $N$ .

If hydrodynamic interactions among cilia are omitted, the cilia-driven velocity and the thrust force may be proportional to the number of cilia $N$ . Using resistive force theory, the force generated by a single cilium $F$ can be scaled as (Osterman & Vilfan 2011)

(3.1) $$\begin{eqnarray}F\sim c_{N}\unicode[STIX]{x1D714}L^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}$ is the angular velocity of ciliary beat, $L$ is the length of the cilium, $c_{N}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}/(\ln (2L/r)-0.5)$ is the drag coefficient of a slender body (Blake 1974), $\unicode[STIX]{x1D707}$ is the viscosity and $r$ is the radius of the cilium. The swimming velocity of a spherical ciliate without hydrodynamic interaction can then be scaled as

(3.2) $$\begin{eqnarray}U\sim \frac{c_{N}\unicode[STIX]{x1D714}L^{2}}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a_{0}}N.\end{eqnarray}$$

The above scaling can be used only for small numbers of cilia and is qualitatively different from the numerical results in figure 10. These results illustrate that hydrodynamic interactions between cilia and the cell body considerably affect swimming velocity.

3.5 Helical swimming with an oblique metachronal wave

Oblique metachronal waves (dexioplectic or laeoplectic waves) are evident in the natural world. To determine the effects of such waves, we considered wave propagations in both the $\unicode[STIX]{x1D703}$ - and $\unicode[STIX]{x1D719}$ -directions on a spherical cellular surface. By setting $\unicode[STIX]{x1D708}=1.0$ and $k=1.0$ , we established dexioplectic wave propagation; the trajectory of the centre of mass over 10 beats is shown in figure 11. The trajectory is helical, because the ciliary beat is no longer instantaneously axisymmetric, although the time-averaged ciliary propulsion force remains axisymmetric. Thus, helical swimming was reproduced using an oblique metachronal wave even when both the cellular body and time-averaged ciliary propulsion force were axisymmetric.

Figure 11. Trajectory of the centre of the spherical ciliate with dexioplectic metachronal wave ( $k=1.0$ and $\unicode[STIX]{x1D708}=1.0$ ). Black dot indicates starting point, while red dot is the end point. A movie can be seen in the supplementary material.

4 Comparison with the Lighthill–Blake squirmer

4.1 Swimming velocity

Here we compare the swimming velocity calculated in the present study to that of the squirmer model. The Lighthill–Blake squirmer (Blake 1971) has been widely used in earlier analytical and numerical studies, as mentioned above in § 1. Using spherical polar coordinates, fluid motion caused by squirming in the radial and azimuthal directions is given by

(4.1) $$\begin{eqnarray}\displaystyle v_{r}(r,\unicode[STIX]{x1D703}) & = & \displaystyle -U\cos \unicode[STIX]{x1D703}+A_{0}\frac{a^{2}}{r^{2}}P_{0}+\frac{2}{3}(A_{1}+B_{1})\frac{a^{3}}{r^{3}}P_{1}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{n=2}^{\infty }\left[\left(\frac{1}{2}n\frac{a^{n}}{r^{n}}-\left(\frac{1}{2}n-1\right)\frac{a^{n+2}}{r^{n+2}}\right)A_{n}P_{n}+\left(\frac{a^{n+2}}{r^{n+2}}-\frac{a^{n}}{r^{n}}\right)B_{n}P_{n}\right]\qquad\end{eqnarray}$$


(4.2) $$\begin{eqnarray}\displaystyle v_{\unicode[STIX]{x1D703}}(r,\unicode[STIX]{x1D703}) & = & \displaystyle U\sin \unicode[STIX]{x1D703}+\frac{1}{3}(A_{1}+B_{1})\frac{a^{3}}{r^{3}}V_{1}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{n=2}^{\infty }\left[\left(\frac{1}{2}n\frac{a^{n+2}}{r^{n+2}}-\left(\frac{1}{2}n-1\right)\frac{a^{n}}{r^{n}}\right)B_{n}V_{n}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{1}{2}n\left(\frac{1}{2}n-1\right)\left(\frac{a^{n}}{r^{n}}-\frac{a^{n+2}}{r^{n+2}}\right)A_{n}V_{n}\right],\end{eqnarray}$$

where $a$ is the squirming surface, i.e. the mean radius of the ciliary envelope, $U$ is the swimming speed of the squirmer, $P_{n}$ is the Legendre polynomial and

(4.3) $$\begin{eqnarray}V_{n}=\frac{2\sin \unicode[STIX]{x1D703}}{n(n+1)}P_{n}^{\prime }.\end{eqnarray}$$

The velocity on the surface of sphere $r=a$ is given by

(4.4a,b ) $$\begin{eqnarray}\displaystyle v_{r}(a,\unicode[STIX]{x1D703})=\mathop{\sum }_{n=0}^{\infty }A_{n}(t)P_{n}(\cos \unicode[STIX]{x1D703}),\quad v_{\unicode[STIX]{x1D703}}(a,\unicode[STIX]{x1D703})=\mathop{\sum }_{n=1}^{\infty }B_{n}(t)V_{n}(\cos \unicode[STIX]{x1D703}). & & \displaystyle\end{eqnarray}$$

The swimming velocity $U$ and the coefficients $A_{1}$ and $B_{1}$ satisfy the following relation:

(4.5) $$\begin{eqnarray}\displaystyle U={\textstyle \frac{1}{3}}(2B_{1}-A_{1}). & & \displaystyle\end{eqnarray}$$

If the ciliary envelope is set outside the individual cilia, equation (4.5) should describe the swimming velocity of our present model. However, if the ciliary envelope overlaps with individual cilia, a propulsion force is generated even outside the envelope and (4.5) is no longer valid. We thus calculated the velocities $v_{r}$ and $v_{\unicode[STIX]{x1D703}}$ at $r=a_{0}+L+\unicode[STIX]{x1D716}$ , where $a_{0}$ is the cell radius and $L$ is the length of a cilium. Here $\unicode[STIX]{x1D716}$ is an offset used to avoid the singularity and to increase the accuracy of numerical integration; $\unicode[STIX]{x1D716}/L$ was set to 0.1. We averaged $v_{r}$ and $v_{\unicode[STIX]{x1D703}}$ in the circumferential direction, and estimated $A_{n}$ and $B_{n}$ of our present ciliate model using a least-squares method to fit the data to (4.4). We used 50 points in the $\unicode[STIX]{x1D703}$ -direction to estimate $A_{1}$ and $B_{1}$ , and the swimming velocity was then calculated employing (4.5).

Figure 12. Comparison with the squirmer. Solid lines express present numerical results, whereas symbols are the theory (Blake 1971). (a) Effect of the cilia number $N$ on the temporal swimming speed with in-phase wave. (b) Temporal swimming speed with antiplectic and symplectic waves $(N=160)$ . (c) Difference of maximum swimming speed between the present numerical results (in-phase beating) and the squirmer.

The time changes in swimming velocities outlined in this study and squirmers with three different types of metachronism are shown in figure 12. In all cases, there is quantitative agreement between the two models. We also compared maximum swimming speeds between our numerical model and the squirmers. As shown in figure 12(c), we found quantitative agreement, and the difference between the two models was less than 1 % when $N>10$ . We thus confirmed that swimming velocity can be predicted by the Lighthill–Blake theory when the ciliary envelope is set outside the individual cilia.

4.2 Energy dissipation

In the squirmer model, the instantaneous rate of working of the stress at the squirming surface $a$ is given by

(4.6) $$\begin{eqnarray}\displaystyle P=2\unicode[STIX]{x03C0}a^{2}\int _{0}^{\unicode[STIX]{x03C0}}(v_{r}\unicode[STIX]{x1D70E}_{rr}+v_{\unicode[STIX]{x1D703}}\unicode[STIX]{x1D70E}_{r\unicode[STIX]{x1D703}})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D748}$ is the viscous stress. In an isolated system, the above $P$ -value reflects the rate of viscous energy dissipation via fluid motion. In the squirmer model, fluid motion inside the ciliary envelope is not considered; thus, viscous energy dissipation is only considered outside the envelope. We calculated energy dissipation inside the ciliary envelope using our model.

Figure 13. Power generated by ciliary beat with three different metachronal waves.

The powers associated with movement of the cell body and cilia are given by

(4.7) $$\begin{eqnarray}\displaystyle P_{cell}=\int _{cell}\boldsymbol{q}\boldsymbol{\cdot }(\boldsymbol{V}+\unicode[STIX]{x1D734}\wedge \hat{\boldsymbol{r}})\,\text{d}A & & \displaystyle\end{eqnarray}$$


(4.8) $$\begin{eqnarray}\displaystyle P_{cilia}=\mathop{\sum }_{i=1}^{N}\int _{cilia}\boldsymbol{f}\boldsymbol{\cdot }(\boldsymbol{V}+\unicode[STIX]{x1D734}\wedge \hat{\boldsymbol{r}}+\boldsymbol{v}^{cilia})\,\text{d}s_{i}. & & \displaystyle\end{eqnarray}$$

The rate of energy dissipation occurring outside the envelope can be given by

(4.9) $$\begin{eqnarray}\displaystyle P_{out}=\int _{a}^{\infty }\unicode[STIX]{x1D748}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{v}\,\text{d}V, & & \displaystyle\end{eqnarray}$$

where $a=a_{0}+L+\unicode[STIX]{x1D716}$ is the squirming surface. Considering the energy conservation, the rate of energy dissipation inside the envelope $P_{in}$ is given by

(4.10) $$\begin{eqnarray}P_{in}=P_{cell}+P_{cilia}-P_{out}.\end{eqnarray}$$

To calculate the infinite integral (4.9), we split $P_{out}$ into two parts,

(4.11) $$\begin{eqnarray}\displaystyle P_{out}=P_{out}^{near}+P_{out}^{far}=\int _{a}^{r_{c}}\unicode[STIX]{x1D748}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{v}\,\text{d}V+\int _{r_{c}}^{\infty }\unicode[STIX]{x1D748}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{v}\,\text{d}V, & & \displaystyle\end{eqnarray}$$

where $r_{c}$ is the domain truncation. As the stress $\unicode[STIX]{x1D748}$ is proportional to $r^{-3}$ in the far field, the energy $\unicode[STIX]{x1D748}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{v}$ is proportional to $r^{-6}$ (which can be numerically confirmed (data not shown)). We then assumed that $\unicode[STIX]{x1D748}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{v}\propto r^{-6}$ applied in the far field. In (4.11) $P_{out}^{far}$ was analytically calculated at $r=r_{c}$ using numerical data. The domain truncation $r_{c}$ was set to $r_{c}/L=10$ , which is sufficiently large to allow computation of the near-field contribution. In the near field, $P_{out}^{near}$ in (4.11) was numerically calculated via Gaussian numerical integration. The details are given in our earlier work (Omori et al. 2017).

We first computed the power of the ciliary beat over a cycle with three different metachronal waves (in-phase, antiplectic and symplectic), and the results are shown in figure 13. The power generated was highest when cilia were in-phase. However, due to the increased duration of the recovery stroke, the rate of working over the cycle was small. In the case of the antiplectic wave, on the other hand, the effective time was longer, but the peak decreased, as did the swimming velocity (cf. figure 9). In the antiplectic mode, the rate of working of each cilium was large during the effective stroke and the wave propagated to neighbouring cilia as in Blake (1972). Thus, the total period of effective power was longer than for in-phase beating.

Figure 14. Effect of cilia number on the energy dissipation of a ciliate model. (a) Energy dissipations occurred outside and inside the ciliary envelope. (b) Efficiency of the swimming.

We then calculated energy dissipation occurring inside and outside the ciliary envelope, and results with in-phase beating are shown in figure 14(a). The energy loss inside the envelope, $P_{in}$ , rapidly increases with $N$ , and approximately 90 % of the energy is dissipated within the envelope when $N=320$ . Thus, most energy dissipation occurs within the envelope; ciliary motion, rather than fluid motion outside the envelope, explains most of the dissipation. Keller & Wu (1977) estimated energy dissipation inside the ciliary envelope analytically and found that over 95 % of dissipation occurred inside the envelope. They assumed that each cilium worked at the same rate and did not consider hydrodynamic interactions among cilia, unlike the current study. If we neglect hydrodynamic interactions, the total power due to ciliary beating should be proportional to the number of cilia. To clarify the effects of hydrodynamic interactions, we additionally calculated the power of a single cilium located on a flat plane; the result of this assessment was $\bar{P}_{single}T^{2}/\unicode[STIX]{x1D707}L^{3}=9.74$ . Multiplying by $N=320$ , we have $N\bar{P}_{single}T^{2}/\unicode[STIX]{x1D707}L^{3}=3116.8$ . This value is approximately twice as large as the numerical result in figure 14(a). This implies that hydrodynamic interactions within ciliary layers diminish the rate of working. Although the magnitude of ciliary power differed, the ratio of energy dissipation inside the envelope was almost the same in the present study as in the study by Keller & Wu (1977). This may be because both swimming velocity and energy dissipation outside the envelope decrease as a reflection of the reduced power of individual cilia.

When the energy dissipation is known, it is possible to calculate swimming efficiency. Swimming efficiency is defined as the time-average swimming speed and the power, e.g.  $6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a_{0}\bar{V}^{2}/\bar{P}$ for a spherical ciliate. In the case of in-phase beating, the efficiency of the ciliate model is 10-fold lower than that of the squirmer model, as shown in figure 14(b).

We then investigated the effect of the aspect ratio on efficiency. As stated in former studies (Pironneau 1973; Vilfan 2012), the optimal profile for minimum dissipation is a prolate rather than an oblate ellipsoid, and the drag coefficient of a rigid spheroid is minimised when $\unicode[STIX]{x1D6FC}_{v}\sim 2$ (cf. figure 7 b). Swimming efficiency, on the other hand, has a different tendency due to the effect of the metachronal wave, as shown in figure 15. The efficiency of a symplectic wave increases with the aspect ratio, whereas an antiplectic wave is maximal at $\unicode[STIX]{x1D6FC}_{v}\sim 1$ .

Figure 15. Swimming efficiency with various aspect ratios $\unicode[STIX]{x1D6FC}_{v}$ ( $N=160$ ).

Figure 16. Effect of metachronal wave on the efficiency.

Finally, we investigated the effect of a metachronal wave on swimming efficiency. The ciliate shape was again set to a sphere with radius $a_{0}/L=3.0$ , and the result is shown in figure 16. The optimal efficiency was given by the antiplectic wave, and the value monotonically increased with the number of cilia. The efficiency of the antiplectic wave was approximately 6.7 times larger than that of the in-phase wave when $N=320$ . Several studies (Blake 1972; Michelin & Lauga 2010; Osterman & Vilfan 2011) also found that an antiplectic metachronal wave resulted in optimal efficiency. Thus, we confirmed qualitative agreement between the present and former studies. We also found that the ratio of energy dissipation inside and outside of the ciliary envelope was not significantly modified by the ciliary beat configuration. Hence, the antiplectic wave becomes the most efficient not only by considering the dissipation outside the envelope but also by considering both inside and outside dissipations.

Blake (1974) investigated the optimal spacing of cilia for efficiency and found efficiency to be maximum at $\unicode[STIX]{x0394}x/L=0.25$ for antiplectic waves, where $\unicode[STIX]{x0394}x$ is the spacing. In our simulations, swimming efficiency increased monotonically with $N$ when cilia were beating in antiplectic mode (cf. figure 16), and the minimum spacing was $\unicode[STIX]{x0394}x/L=0.59$ with $N=320$ . This value is approximately 2.3 times larger than Blake’s estimation. Osterman & Vilfan (2011) numerically calculated ciliary beating efficiency by considering the losses inside the ciliary layer, and reported that the maximum efficiency was given by an antiplectic wave with an inter-ciliary distance $\unicode[STIX]{x0394}x\sim 0.25L$ . The resulting optimal swimming efficiency was 0.016, which is approximately 10 times larger than our numerical result of the antiplectic wave with $N=320$ (cf. figure 16). The discrepancy may arise from the difference in the number of cilia between the two studies: at most 320 cilia were used in the present study, whereas an infinite number of cilia with periodic boundary was set in Osterman & Vilfan (2011). Since the efficiency of an antiplectic wave monotonically increases with the number of cilia in the present study, the difference may be reduced by further increasing $N$ . Another possibility for the discrepancy is the difference in the ciliary beat pattern between the two studies. The amplitude of the ciliary beat in Osterman & Vilfan (2011) is larger than ours, which may make their efficiency larger. We could not employ the same beat pattern in this study, because the cilia tend to collide frequently with each other, and non-hydrodynamic forces start to play a role.

Katsu-Kimura et al. (2009) also experimentally estimated the energetic efficiency of swimming in Paramecium as 0.078 %, which includes losses in metabolism. This value is still above the results of the present study. We think the ciliary configurations used in this study, such as number and density of cilia, and ciliary beat pattern, caused the quantitative difference from the experimental result of Katsu-Kimura et al. (2009). We would like to address these effects in our future studies by overcoming numerical difficulties caused by the collision of cilia.

5 Conclusions

We developed a ciliate model using the slender-body–boundary integral coupling method. We found that an antiplectic metachronal wave of wavenumber $k=1.0$ optimised the swimming speed when the cellular body assumed various aspect ratios. We used oblique wave propagation to reproduce the helical swimming trajectories of even spherical ciliates exerting axisymmetric, time-averaged ciliary forces. We also compared our numerical model to the squirmer model of Lighthill and Blake. The swimming speed afforded by our discrete ciliary model agreed with that derived using Lighthill–Blake theory. However, the squirmer model considers only events outside the ciliary envelope, and thus fails to accurately estimate the energy cost of swimming; over 90 % of the energy was dissipated inside the ciliary envelope. The results of this study show that swimming efficiency was optimal with an antiplectic wave, which is consistent with former studies. The efficiency was 6.7 times larger than that of in-phase beating when the number of cilia was $N=320$ . Our findings provide a fundamental basis for modelling swimming micro-organisms.


This research was supported by JSPS KAKENHI grants, numbers 17H00853 and 18K18354.

Supplementary movies

Supplementary movies are available at

Appendix A. Validation of the numerical method

In this appendix, we explain the accuracy of our numerical method. We validated the numerical method by calculating the viscous drag of a rigid sphere or a slender ellipsoid. We first calculated the drag of a sphere in a uniform flow by using a boundary element method. We compared the numerical results with the analytical solution of Stokes’ law, $6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aU$ , where $U$ is the velocity, $\unicode[STIX]{x1D707}$ is the viscosity and $a$ is the radius of the sphere. We confirmed that the error became less than 0.1 % when we used 5120 triangular elements or larger meshes. We then calculated the drag of a slender ellipsoid in a uniform flow. We set two different flow directions, as shown in figure 17(a), and compared the results using slender-body theory and the boundary element method. To estimate numerical accuracy, we defined the force difference between the slender-body theory (SBT) and the boundary element method (BEM) as

(A 1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}F=\frac{|F^{SBT}-F^{BEM}|}{F^{BEM}}. & & \displaystyle\end{eqnarray}$$

The results are shown in figure 17(b). The force differences are less than 1 % when the slenderness $\unicode[STIX]{x1D700}=r/L$ is smaller than $10^{-2}$ in both the normal and tangential directions ( $\unicode[STIX]{x0394}F_{n}$ and $\unicode[STIX]{x0394}F_{t}$ ), where $r$ is the semi-minor axis and $L$ is the length. In the case of a cilium, a typical length is $10~\unicode[STIX]{x03BC}\text{m}$ and a typical radius is 100 nm; thus $\unicode[STIX]{x1D700}$ can be estimated as $10^{-2}$ . This confirms that our numerical method has sufficient accuracy to simulate the swimming of ciliates.

Figure 17. Comparison with a slender-body theory (SBT) and a boundary element method (BEM). (a) Uniform flow around a slender body with $\unicode[STIX]{x1D700}=10^{-2}$ . We set two different flow directions. (b) Force difference between the slender-body theory and the boundary element method.


Blake, J. 1971 A spherical envelope approach to ciliary propulsion. J. Fluid Mech. 46, 199208.
Blake, J. 1972 A model for the micro-structure in ciliated organisms. J. Fluid Mech. 55, 123.
Blake, J. 1974 Hydrodynamic calculations on the movements of cilia and flagella I. Paramecium . J. Theor. Biol. 45, 183203.
Brennen, C. 1974 An oscillating-boundary-layer theory for ciliary propulsion. J. Fluid Mech. 65, 799824.
Brumley, D. R., Polin, M., Pedley, T. J. & Goldstein, R. E. 2012 Hydrodynamic synchronization and metachronal waves on the surface of the colonial alga Volvox carteri . Phys. Rev. Lett. 109, 268102.
Catmull, E. & Rom, R. 1974 A class of local interpolating splines. Computer Aided Geometric Design. Academic Press.
Elgeti, J. & Gompper, G. 2013 Emergence of metachronal waves in cilia arrays. Proc. Natl Acad. Sci. USA 110, 44704475.
Eloy, C. & Lauga, E. 2012 Kinematics of the most efficient cilium. Phys. Rev. Lett. 109, 038101.
Evans, A. A., Ishikawa, T., Yamaguchi, T. & Lauga, E. 2011 Orientational order in concentrated suspensions of spherical microswimmers. Phys. Fluids 23, 111702.
Fulford, G. R. & Blake, J. R. 1986 Muco-ciliary transport in the lung. J. Theor. Biol. 121, 381402.
Guo, H., Nawroth, J., Ding, Y. & Kanso, E. 2014 Cilia beating patterns are not hydrodynamically optimal. Phys. Fluids 26, 091901.
Ishikawa, T., Kajiki, S., Imai, Y. & Omori, T. 2016 Nutrient uptake in a suspension of squirmers. J. Fluid Mech. 789, 481499.
Ishikawa, T. & Kikuchi, K. 2018 Biomechanics of Tetrahymena escaping from a dead end. Proc. R. Soc. Lond. B 285, 20172368.
Ishikawa, T., Locsei, J. T. & Pedley, T. J. 2008 Development of coherent structures in concentrated suspensions of swimming model micro-organisms. J. Fluid Mech. 615, 401431.
Ishikawa, T. & Pedley, T. J. 2007 Diffusion of swimming model micro-organisms in a semi-dilute suspension. J. Fluid Mech. 588, 437462.
Katsu-Kimura, Y., Nakaya, F., Baba, S. A. & Mogami, Y. 2009 Substantial energy expenditure for locomotion in ciliates verified by means of simultaneous measurement of oxygen consumption rate and swimming speed. J. Expl Biol 212, 18191824.
Keller, S. R. & Wu, T. Y. 1977 A porous prolate-spheroidal model for ciliated micro-organisms. J. Fluid Mech. 80, 259278.
Kessler, J. O. 1986 Individual and collective fluid dynamics of swimming cells. J. Fluid Mech. 173, 191205.
Lighthill, M. J. 1952 On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers. Commun. Pure Appl. Maths 5, 109118.
Machemer, H. 1972 Ciliary activity and the origin of metachrony in Paramecium: effects of increased viscosity. J. Expl Biol. 57, 239259.
Magar, V., Goto, T. & Pedley, T. J. 2003 Nutrient uptake by a self-propelled steady squirmer. Q. J. Mech. Appl. Maths 56, 6591.
Michelin, S. & Lauga, E. 2010 Efficiency optimization and symmetry-breaking in a model of ciliary locomotion. Phys. Fluids 22, 111901.
Nasouri, B. & Elfring, G. J. 2016 Hydrodynamic interactions of cilia on a spherical body. Phys. Rev. E 93, 033111.
Niedermayer, T., Eckhardt, B. & Lenz, P. 2008 Synchronization, phase locking, and metachronal wave formation in ciliary chains. Chaos 18, 037128.
Omori, T., Sugai, H., Imai, Y. & Ishikawa, T. 2017 Nodal cilia-driven flow: development of a computational model of the nodal cilia axoneme. J. Biomech. 61, 242249.
Osterman, N. & Vilfan, A. 2011 Finding the ciliary beating pattern with optimal efficiency. Proc. Natl Acad. Sci. USA 108, 1572715732.
Pedley, T. J. 2016 Spherical squirmers: models for swimming micro-organisms. IMA J. Appl. Maths 81, 488521.
Pedley, T. J., Brumley, D. R. & Goldstein, R. E. 2016 Squirmers with swirl: a model for Volvox swimming. J. Fluid Mech. 798, 165186.
Pironneau, O. 1973 On optimum profiles in Stokes flow. J. Fluid Mech. 59, 117128.
Pozrikidis, C. 1992 Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press.
Saintillan, D. & Shelley, M. 2008 Instabilities, pattern formation, and mixing in active suspensions. Phys. Fluids 20, 123304.
Simha, R. A. & Ramaswamy, S. 2002 Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles. Phys. Rev. Lett. 89, 058101.
Sleigh, M. A. 1962 The Biology of Cilia and Flagella. Pergamon Press.
Tornberg, A.-K. & Shelly, M. J. 2004 Simulating the dynamics and interactions of flexible fibers in Stokes flows. J. Comput. Phys. 196, 840.
Vilfan, A. 2012 Optimal shapes of surface slip driven self-propelled microswimmers. Phys. Rev. Lett. 109, 128105.