Swimming mediated by ciliary beating: comparison with a squirmer model

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.


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 Swimming mediated by ciliary beating 775 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(Blake , 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.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.

Boundary integral equation with slender-body theory
Consider a ciliate immersed in an infinite Newtonian liquid of density ρ and viscosity µ; 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 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 ∈ [0, L], where L is the ciliary length.The flow field at point x located on the ith cilium, x ∈ s i , is given by Pozrikidis (1992) and Tornberg & Shelly (2004) as where q is the viscous traction on the cell body, 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 J is given by where r = |r| and r = x − y.In (2.1) Λ and K are the local operators of the slenderbody theory, which are given by (Tornberg & Shelly 2004) and where c = − ln(ε 2 e), t is the unit tangential vector to the centreline of each cilium and ε = a cilia /L, with a cilia the ciliary radius.In terms of the practical ciliary radius and length ratio (Sleigh 1962), the slenderness value ε is set to ε = 0.01 throughout the present study.The slender-body kernel W is defined by (2.5) When the observation point x is not on a cilium, x / ∈ s, the velocity is given by v (2.6) Though the asymptotic accuracy of kernels Λ and K is O(ε 2 ln ε), kernel W is accurate only to the limit O(ε) (Tornberg & Shelly 2004).Equations (2.1) and (2.6) are therefore accurate to the limit O(ε).

Ciliary motions
We first define an orthonormal frame for the cell body e i with origin x c , where x c is the body's centre of mass; e 1 then reflects swimming orientation.To efficiently model ciliary motion on the cell surface, local vectors with an orthonormal basis g i are defined as follows.A material point x b , located at the base of a cilium on the cell surface, serves as the origin of the orthonormal body frame g i (cf.figure 1a).The basis vectors g 1 and g 2 are defined as , where b = e 1 ∧ n and n is the outward unit normal vector, respectively.778 H. Ito, T. Omori and T. Ishikawa  (1986).
The time-dependent profile of each ciliary motion is derived using the following mathematical formula of Fulford & Blake (1986): where with ω the angular beat frequency.The Fourier coefficients α i n and β i n are given by

9a,b)
The coefficients A i mn and B i mn 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).

Boundary element method
In order to simulate free swimming of the ciliate model, force-free and torque-free conditions are taken into account: (2.11) where r = x − x c , and 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 x ∈ cilia, (2.12) where V is the translational and Ω the angular velocity, and v cilia is the ciliary velocity with respect to the body frame e i (i.e.v cilia = ∂x cilia /∂t).We then solve the resistance problem with respect to the unknowns V and Ω, and derive the viscous tractions q and 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 x is located on the cell body, the following linear algebraic equation can be derived from (2.6): When x is on the cilia, on the other hand, we have the following equation from (2.1): (2.14) The vector sizes of ] and [ 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 [J bb  ] and ] are 3N c × 3N b and 3N c × 3N c , respectively.In a similar manner, discretised forms of the force-torque conditions, (2.10) and (2.11), can be written as and Considering the boundary condition of (2.12), the system can be expanded to (2.17) The matrix components V b and A b are of sizes 3N b × 3, whereas V c and A c are both 3N c × 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.

Parameter setting
In a Stokes flow regime, the fluid viscosity is simply a multiplier of both force and traction.The viscosity µ 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 ψ 0 (θ , φ) was defined as where k and ν are the wavenumbers in the θand φ-directions, respectively (cf.figure 1a).The angle θ = [0, π] is that between the orientation vector e 1 and rb ; and ; the sign is determined by the sign of e 3 • rb .
In the following sections, the wavenumber k varies from −2.0 to 2.0, and ν 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 ν = 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 α v .Next, the major radius was fixed in various aspect ratio, represented as α l .For all computations, the time interval t was set to t/T = 0.01, where T is the ciliary beat period.

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 ν in (2.18) were set to zero.The orientation vector was initially set to 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 t/T 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 backand-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 t/T 0.25, ciliary beating was effective, and swimming velocity was maximal at t/T ∼ 0.1.The maximum swimming velocity was approximately V x T/L 3.0, equivalent to a halfbody-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 = T 0 V x dt was positive; the ciliate moved forwards.
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 f x ds 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  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.

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.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 ν 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 Vx at various k-values are shown in figure 5(b); Vx is maximum when k 1.0, suggesting that an antiplectic metachronal wave with k = 1.0 is optimal in this cellular configuration.
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.

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 α v and α l .The swimming behaviours of an oblate ciliate (α v = 0.5) and a prolate ciliate (α 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 α v and α 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 H. Ito, T. Omori and T. Ishikawa a rigid spheroid became the minimum when α v ∼ 2, as shown in figure 7(b).Thus, large velocities were observed around α 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 α 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 α l in the moderate-aspect-ratio regime, as shown in figure 8.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 µm and the swimming velocity was 440 µm s −1 (i.e.6.8 body lengths per second).In figure 7 In phase (k = 0) In phase (k = 0) Symplectic (k = -1) average velocity when α 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.
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 ν 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 3.7 when N = 320 at in-phase beating, and was V x T/L 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.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) F ∼ c N ωL 2 , (3.1) where ω is the angular velocity of ciliary beat, L is the length of the cilium, c N = 4πµ/(ln(2L/r) − 0.5) is the drag coefficient of a slender body (Blake 1974), µ 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 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 θand φ-directions on a spherical cellular surface.By setting ν = 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.

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 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 The velocity on the surface of sphere r = a is given by The swimming velocity U and the coefficients A 1 and B 1 satisfy the following relation: 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 θ at r = a 0 + L + , where a 0 is the cell radius and L is the length of a cilium.Here is an offset used to avoid the singularity and to increase the accuracy of numerical integration; /L was set to 0.1.We averaged v r and v θ 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 θ -direction to estimate A 1 and B 1 , and the swimming velocity was then calculated employing (4.5).
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.

Energy dissipation
In the squirmer model, the instantaneous rate of working of the stress at the squirming surface a is given by where σ 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.The powers associated with movement of the cell body and cilia are given by where r c is the domain truncation.As the stress σ is proportional to r −3 in the far field, the energy σ : ∇v is proportional to r −6 (which can be numerically confirmed (data not shown)).We then assumed that σ : ∇v ∝ r −6 applied in the far field.In (4.11) P far out 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 near out 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.
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)

FIGURE 1 .
FIGURE 1. Problem setting.(a) Schematics of the body frames and the angles θ and φ.(b) Beat pattern of each cilium described by the parameters in table 1.

FIGURE 2 .
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 f x ds.A movie can be seen in the supplemental material available at https://doi.org/10.1017/jfm.2019.490.
FIGURE 3. (Colour online) (a) Swimming velocity of the ciliate as a function of time.(b) Propulsion force of individual cilia in three regions calculated by i f x ds i .Regions A, B and C are posterior, middle and anterior regions, as shown in the inset.

FIGURE 4 .
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.
FIGURE 5. (Colour online) 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 7 .FIGURE 8 .
FIGURE 7. Effect of aspect ratio α 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 α v .Here U is the fluid velocity relative to the rigid spheroid.

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

FIGURE 10 .
FIGURE 10.Average swimming velocity with different number of cilia N.
FIGURE 11. (Colour online) Trajectory of the centre of the spherical ciliate with dexioplectic metachronal wave (k = 1.0 and ν = 1.0).Black dot indicates starting point, while red dot is the end point.A movie can be seen in the supplementary material.

H
FIGURE 12. (Colour online) 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.
Brumley et al. (2012) 776 H. Ito, T. Omori and T. Ishikawa 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.