Multi-ciliated microswimmers–metachronal coordination and helical swimming

Abstract The dynamics and motion of multi-ciliated microswimmers with a spherical body and a small number N (with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5< N < 60$$\end{document}5<N<60) of cilia with length comparable to the body radius, is investigated by mesoscale hydrodynamics simulations. A metachronal wave is imposed for the cilia beat, for which the wave vector has both a longitudinal and a latitudinal component. The dynamics and motion is characterized by the swimming velocity, its variation over the beat cycle, the spinning velocity around the main body axis, as well as the parameters of the helical trajectory. Our simulation results show that the microswimmer motion strongly depends on the latitudinal wave number and the longitudinal phase lag. The microswimmers are found to swim smoothly and usually spin around their own axis. Chirality of the metachronal beat pattern generically generates helical trajectories. In most cases, the helices are thin and stretched, i.e., the helix radius is about an order of magnitude smaller than the pitch. The rotational diffusion of the microswimmer is significantly smaller than the passive rotational diffusion of the body alone, which indicates that the extended cilia contribute strongly to the hydrodynamic radius. The swimming velocity is found to increase with the cilia number N with a slightly sublinear power law, consistent with the behavior expected from the dependence of the transport velocity of planar cilia arrays on the cilia separation. Graphic abstract Supplementary Information The online version contains supplementary material available at 10.1140/epje/s10189-021-00078-x.


Introduction
Cilia and flagella are the ubiquitous machinery in eukaryotic cells and organisms to generate fluid flow and to propel cells and microorganisms in a fluid environment [1,2].While flagella have the beat pattern of a sinusoidal traveling wave, and are usually employed to propel single cells like sperm [3], cilia have two distinct phases in their beat cycle -the power and the recovery stroke -, and often work together in pairs like in Chlamydomonas reinhardtii [4], or in large cilia carpets.Examples for the concerted action of cilia in carpets are the transport of mucus in the airways [5,6], the flow generation of the cerebrospinal fluid in brain ventricles [7][8][9], and the swimming of multi-ciliated microorganisms, such as the green alga Volvox [10], the protozoan Paramecium [11], and the placidozoan Opalina [12].
A remarkable feature of cilia carpets is that the beat is highly coordinated in the form of metachronal waves [12][13][14][15][16][17], where the beat cycles of neighboring cilia have a fixed phase shift.The powerstroke direction can be parallel or antiparallel to the propagation direction of the metachronal wave, which is denoted as symplectic or antiplectic wave, but can also point right-wise or left-wise from the wave direction, which is denoted dexioplectic or laeoplectic wave.The later wave form clearly requires some chirality in the system, which can either be in the aplanarity of the ciliary beat [18], or result from the spatial arrangement of cilia.
The origin and effect on the transport efficiency of the ciliary beat and of the metachronal wave is fascinating, and has thus been investigated intensively.
It is now well established that hydrodynamic interactions are strong enough to cause coordination between neighboring cilia [19][20][21][22][23], and suffice to explain the formation of metachronal waves in cilia arrays [17].It has also been shown that the transport efficiency can be much higher than for perfect beat synchronization, related to the fact that always a fraction of the cilia is in the power stroke, thus avoiding a forward-backward motion of the fluid.Additionally, flow generated by the power stroke at the wave crest experiences only small resistance from cilia in the neighbouring wave trough, which can all remain close to the anchoring surface during the recovery stroke -both without much steric hindrance [14,17,18].A further interesting issue is the coordination of the beat directions in cilia carpets, which is now believed to be a self-organized process mediated by the fluid flow [6,7,14,22,23].On the other hand, the synchronization of the beat of the two flagella of Chlamydomonas reinhardtii arise from an elastic mechanical coupling at their basal foot [24,25].
More complex is the cilia coordination in multiciliated spherical or spheroidal microorganisms.One reason is the well-know "hairy-ball" theorem, which states that there is no non-vanishing continuous tangent vector field on a surface of spherical topology [26].This implies that the power-stroke directions of neighboring cilia cannot be parallel everywhere on a spherical surface, but there have to be at least two defects, which can either be of hedgehog or of swirl type.A second reason is that plane metachronal waves are not possible on curved surfaces.
Volvox is a perfect model system for experimental studies of swimming and cilia synchronization [27][28][29].
These studies reveal the existence of a symplectic metachronal wave [27], and that the average metachronal coordination is punctuated by periodic phase defects during which synchrony is partial and limited to specific groups of cells [28].Under conditions of decreasing nutrient concentration, Volvox colonies were found to grow larger and increase their flagellar length, separating the somatic cells further, with the opposing effects of increasing beating force and flagellar spacing balance not significantly affecting the fluid speed at the colony surface [29].
The theoretical description of the swimming of multi-ciliated microorganisms has lead to the early development of the squirmer model [30,31], in which the effect of the ciliary beating is mimicked by a prescribed surface velocity.This model has been generalized more recently to spheriodal shapes [32], and is employed nowadays to describe the collective swimming behavior of many types of microswimmers [1,33].The squirmer model has also been generalized to capture the effect of metachronal waves and to include azimuthal swirl on the continuum level [34].This model predicts mean swimming speeds and angular velocities as a function of the colony radius qualitatively correct, but underestimates both velocities quantitatively [34].
The squirmer model applies in the limit that the cilia length and the separation of their anchoring points on the surface is much smaller than the body size.The dynamics and flow generation of individual cilia becomes more important in the opposite limit.A model with explicit cilia and a prescribed metachronal wave has been introduced recently [35].The model facilitates the calculation of hydrodynamic interactions between cilia and the cell body under free-swimming conditions.
An antiplectic metachronal wave is predicted to be optimal in the swimming speed with various cell-body aspect ratios, which is consistent with former theoretical studies [15,36].The swimming velocity of model ciliates was well represented by the squirmer model.The effect of oblique wave propagation is also briefly touched, and is found to lead to a helical swimming trajectory [35].
Also, there has been significant progress recently on the experimental techniques for cilia characterization [37], as well as the construction of artificial, externally actuated cilia carpets [38], which relates to the goal of the construction of soft microbots [39].Soft robots with antiplectic waves have been shown to exhibit much higher locomotion speed than those with symplectic waves [38].
We employ a similar model of ciliated microorganisms as studied in Ref. [35], but focus on the regime of a smaller number of longer cilia -comparable in length to the body size.Volvocalean algae, which falls into this parameter range are, for example, the 16-celled Pandorina and the 32-celled Eudorina [29,40].We focus on the efficiency of metachronal waves for translational and rotational motion of the swimmer.In particular, we investigate the influence of the wave direction, the phase lag between neighboring cilia, and the cilia density, on the swimming efficiency and the persistence of the swimming trajectory.In particular, we identify parameter combinations for which helical swimming trajectories emerge.

Ciliated microswimmer and hydrodynamic simulation
The ciliated microswimmer is modeled as a spherical body of radius R, to which several cilia are attached (see Fig. 1).The body consists of N b point particles covering the surface of a sphere, with an additional particle at its center.The particles on the surface are connected to their neighbors by stiff harmonic springs to form a triangular network, as well as to the center particle to form an essentially rigid sphere.Following Refs.[17,41], we model a cilium of length L by three semi-flexible polymers, which are arranged on the surface of a (hypothetical) cylinder with triangular cross section.Each semi-flexible polymer is described by a bead-spring model, consisting of N c beads connected by harmonic springs with rest length c and spring constant k c .The three polymers are interconnected by springs in order to retain the cylindrical shape over time.The cilia are anchored in the body with a "clamped" boundary condition, which is achieved by extending the filament by a segment of length R/4 inside the body.This anchoring is implemented by stiff springs which connect the sub-surface part of the cilium to the sphere surface, mimicking the embedding of the basal bodies of the cilia in the cell body.In detail, the first and fourth particle of each of the three polymers constituting the cilium are connected to the closest particle on the sphere surface as well as its next-nearest neighbors.This construction anchors the cilium tightly to the body without a noticeable body deformation -even for the largest applied ciliary forces.
The cilia beat pattern consists of a power and a recovery stroke (see Fig. 2).The ciliary beat is generated by varying the equilibrium spring lengths c of one selected polymer, both spatially along the cilium and periodically in time, which creates a spatially and temporarily varying cilium curvature.The selection of the active polymer defines the beat plane.The beat pattern of power and recovery stoke is obtained by prescribing an analytic function for the desired local curvature of the cilium.For details, see Appendix A. This is inspired by the molecular mechanism which drives ciliary beating, where molecular motors apply torques along the flagellum, but each motor has a maximum force it can generate.To mimic the stall force, and to avoid artificial cilia shapes due to unnaturally large local torques, we limit the change of equilibrium bond length such that a maximum energy of 1.0 k B T per MPC time step (see below) can be inserted into the system, where k B T is the thermal energy.The dynamics of the beat is not only determined by the time-dependence of the internal torques, but is also affected by the flow field around the swimmer and the elastic properties of the cilium.To model the hydrodynamics of the embedding fluid, we employ multi-particle collision dynamics (MPC), a mesoscale simulation technique, which is ideally suited for simulations with a particle-based model of an active microswimmer.In this approach, the fluid consists of point particles, each characterized by its location r i (t) and velocity v i (t).These particles move ballistically during the streaming step for a time interval h.In the subsequent collision step, all particles are sorted into the cells of a simple cubic lattice with lattice constant a. Particles in each cell interact by exchanging momentum, but in such a way to conserve the mass and linear momentum within each cell.A celllevel canonical thermostat (with Maxwell-Boltzmann scaling) is applied after every collision step to maintain a constant temperature T [42].A detailed description of the MPC technique and a review of its application to many systems in soft, active and living matter is provided in Refs.[43][44][45].
For the investigation of the self-organization of beating cilia arrays into metachronal waves in Ref. [17], the switch from power to recovery stroke, and back, needs triggers, which are reached faster or slower depending on the environmental conditions ("geometrical clutch" hypothesis [46]).In Ref. [17], a critical value curvature of the cilium, and a maximum angle between the extended cilium and the normal vector to the cell surface, were employed to trigger these switching events.We use here a variant of this model, in which the force generation is modified such that the switch between power and recovery stroke becomes deterministic in time, with a power-stroke time τ p and a recovery-stroke time τ r , with τ r > τ b , which results in a constant beat period τ b = τ p + τ r .In this case, a phase lag between the beats of neighboring cilia is imposed to induce a metachronal wave (see below).
We consider a ciliated microswimmer, where a total of N cilia are placed equidistantly on three latitudinal rings.In spherical coordinates (ϕ, θ), one ring is located at the equator, θ = 0, and two rings at θ = ±π/4, see Fig. 1.In order to have a roughly constant cilia density, the number of cilia at the "polar circles" is reduced by a factor close to 1/ cos θ = √ 2 compared to their number N eq at the equator.The "first" cilium on each ring is initially placed at ϕ = 0; subsequently, the cilia positions on both polar rings are shifted azimuthally in the same direction by a small ∆ϕ = π/N eq to avoid perfect registry for one particular longitude.Several examples of microswimmer with various numbers of cilia are shown in Fig. 1).
The power-stroke direction of the cilia can be along e θ , i.e. be parallel to the circles of longitude, or deviate from this highly symmetric case.This is modeled by rotating the beat plane of each cilium around the radial axis e r by an angle θ r (with θ r = 0 corresponding to the longitudinal beat direction).A beat-plane orientation with θ r = 0 introduces chirality in the swimmer propulsion pattern.
In the case of an imposed metachronal wave, we define a local phase Ψ(ϕ, θ) for each cilium on the surface of the sphere, with −π < θ < π and 0 ≤ ϕ < 2π, where the direction of the wave is determined by the wave vector k = (k ϕ , k θ ), which has longitudinal and latitudinal components, k θ and k ϕ .Since we want to have continuous wave solution traveling around the circles of constant latitude, k ϕ has to be an integer number.For N eq equally-spaced cilia on the equator, k ϕ = N eq /2 results in a phase lag between neighboring cilia of ∆Ψ = π.Since k ϕ ∈ N, we limit it to k ϕ = 0, 1, ..., (N eq //2 + 1), where // denotes integer division.The longitudinal wave vector k θ ∈ R is not restricted by any physical boundary conditions.We employ the phase lag between neighboring rings to quantify k θ , where ∆θ = π/4 is the fixed latitudinal angular distance between successive rings.For zero longitudinal wave component, k ϕ = 0, a positive phase lag χ corresponds to a symplectic metachronal wave, which travels in the direction of the power stroke, a negative phase lag χ to an antiplectic metachronal wave, which travels in the direction of the recovery stroke.For non-zero longitudinal component  k ϕ , the metachronal wave becomes either dexioplectic (k θ > 0) or laeoplectic (k θ < 0).Because the microswimmer is constructed essentially symmetric to the main axis, dexioplectic waves show the same effect on the propulsion as laeoplectic, except for an opposite direction of axial rotation Ω n .Therefore, we restrict our analysis to symplectic, antiplectic (k ϕ = 0) and dexioplectic (k ϕ > 0) metachronal waves with varying phase lags χ.
For simulations, we employ the following parameters.The spherical body consists of N b = 643 mesh points, the cilia are constructed with N c = 26 beads for each polymer strand, they have length (outside the body) of L = 3R, and the body radius in terms of MPC collision box size a is R = 8a, which guarantees a good resolution of the hydrodynamic flow fields.The MPC simulations employ collisions described by stochastic rotation dynamics, with a time step h = 0.05, cell size a = 1, rotation angle α 0 = 130 • , and an overall simu-lation box size of (100a) 3 .

The 5-7-5 Swimmer with Longitudinal Beat Direction
We focus on the swimming properties of a spherical swimmer with 7 cilia on the equator and 5 cilia on the two polar rings see Fig. 1, with power stroke direction along the main body axis, θ r = 0. Examples for the beating dynamics with phase lag χ = −77 • with latitudinal wave numbers k ϕ = 0 and k ϕ = 1 are shown in Fig. 3 and Fig. 4, respectively.
We consider the propulsion velocity v n , the velocity fluctuations (v n − v n ) 2 around the average, and the rotational velocity Ω n around the main body axis, see Fig. 5.All these quantities depend on both the phase lag χ in the longitudinal direction and the wave number k ϕ in the latitudinal direction.The propulsion velocity v n in the direction of the main body axis shows a pronounced "sinusoidal" dependence on the phase lag χ for all wave numbers k ϕ , see Fig. 5a.Here, the strongest variation is found for a metachronal wave with k ϕ = 0.The swimming velocity v n is maximal for a negative phase lag of χ −50 • (see Fig. 5a) and reaches almost a body length per beat cycle, v n = 0.9R/τ b .Such a negative phase lag corresponds to an antiplectic metachronal wave, with cilia of one ring lagging behind those of the subsequent ring by a little bit more than 1/4 of a beat cycle, where the wave travels against the direction of the power stroke.The smallest velocity of v n = 0.2R/τ b corresponds to symplectic wave with a positive phase lag of χ 50 • , where the wave travels with the direction of the power stroke.These results are consistent with former theoretical studies of efficiency optimization [15,16].
The variation (v n − v n ) 2 of the swimming velocity is shown in Fig. 5b.The swimming velocity v n fluctuates mostly for swimmers with k ϕ = 0, whereas it is nearly independent of the phase lag χ for latitudinal wave numbers with k ϕ > 0. In particular for the synchronous case, with χ = 0, the swimmer moves quickly forward during the power stroke, but reverses its direction of motion during the recovery stroke, which in sum leads to a relative slow average velocity with high fluctuations.
For k ϕ ≥ 1, the latitudinal component of the wave leads to an additional rotation velocity Ω n around its main body axis (see Fig. 5c).The rotation is most pronounced for k ϕ = 1 and −50 • < χ < 50 • , and is very small for k ϕ ≥ 2 and all phase lags χ.This behavior can be understood by considering two contributions.(i) For k ϕ ≥ 1, rotation is enhanced when the metachronal wave travels in the latitudinal direction, i.e. χ 0. (ii) Latitudinal wave numbers k ϕ ≥ 2 imply short metachronal wave lengths, which cannot propel fluid effectively, because the opposing beat of neighboring cilia just generates local swirls.For example, k ϕ = 3 corresponds to a phase lag between neighboring cilia on the equator of ∆Ψ = k ϕ 2π/N eq , which implies a large phase lag of ∆Ψ 150 • (for N eq = 7).The rotational component of the propulsion also slightly reduces the component that contributes to the swimming velocity, compare Fig. 5a.
For k ϕ ≥ 1, the chirality of the wave pattern not only induces a body rotation, but also implies a helical swimming trajectory [41,47].Consider a swimmer, for which the main body axis n rotates around a fixed axis e in the lab reference frame with an opening angle α.In the absence of translational and rotational noise, this corresponds to the trajectory of a perfect helix with axis e and azimuthal direction e ⊥ (t) where e ⊥ (t) = (cos(Ω c t), sin(Ω c t)) in Cartesian coordinates in the plane with normal vector e .The opening angle α, the rotation frequency Ω c , and swim velocity v n , are related to the helix parameters -helix radius R h , pitch length P h , an helix angle α h -as cos(α) and For α = 0 the microswimmer moves on a straight lines, whereas for α = 90 • it moves on a circle.In general, the directional auto-correlation function of a swimmer consists of two factors, the correlation due to the helical motion, which depends on the inclination angle α, and an exponential decay due to thermal or active noise, Correlation functions for fixed χ = −50 • and various k ϕ are shown in Fig. 5d.Results for the fitted parameters of the helical motion and the decay time 1/κ in Eq. ( 6) are displayed in Fig. 6.The calculation of these parameters requires very long simulation times of several hundred beat periods τ b or more.Even in this case, trajectories are sometimes too short for a reliable parameter estimation, in particular for large decorrelation times 1/κ.
For most cases, the helix angle α < π/4, which means that the constant term in Eq. (3) dominates over the oscillatory term, and thus the helix is thin and elongated.Only in a few cases, like k ϕ = 3, we find a nearly circular motion and a tightly wound helix.The circling frequency Ω c is pronounced for k ϕ = 1, and small phase lags |χ|, which leads to a pronounced helical swimming trajectory, see Fig. 7.This is closely related to the large internal spinning frequency Ω n (see Fig. 5c).For k ϕ ≥ 2, the circling frequency is typically very small, which is again related to the inefficiency of propulsion for short metachronal wave lengths.The data in Fig. 6 sometimes appear to have a larger "scatter" for different wave numbers k ϕ and different phase lags χ.However, two points should be noticed.(i) We expect smooth curves for fixed k ϕ as a function of χ when χ is varied in small steps; however, we vary χ in rather large discrete steps of about 50 • , which results in significantly different metachronal waves.(ii) Similarly, different k ϕ generate very different wave patterns, compare Figs. 3 and 4 for k ϕ = 0 and k ϕ = 1, respectively.
The ratio Ω c /|Ω n | is displayed in Fig. 8; it demonstrates that the circling and spinning frequencies are closely related.In many cases, Ω c /|Ω n | 1, which  corresponds to a "twisted-ribbon-like" motion (or to a "tidal-locking-like" motion, as for the moon, which always presents the same side to the earth).However, there are also cases where Ω c /|Ω n | is close to 0 or to around 2, which indicates that the two frequencies don't have to be locked always.An example for Ω c /|Ω n | = 0 is a microswimmer which spins around its body axis, but swims on a straight trajectory.
It is interesting to note that the microswimmers with k ϕ = 0 do not have a vanishing circling frequency.The reason is that they are not perfectly axisymmetric, and therefore have some inherent chirality, due to the non-symmetric location of the cilia on the body.For example, there is a particular longitude, where the cilia in the various rings are closest together, which generates a non-axisymmetric flow field.This leads to the helical trajectory displayed in Fig. 7, with a nonvanishing but very small Ω c .The spinning frequency Ω n around the main body axis nearly vanishes in this case (see Fig. 5c), because for θ r = 0 there is essentially no component of the cilia beat in the latitudinal direction.
Finally, the decay (or decorrelation) time 1/κ is determined by thermal fluctuations and Stokes friction, and possibly by active fluctuations generated by the beat.In the thermal case, we expect a rotational diffusion time τ rot = 1/(2D rot ) with rotational diffusion coefficient D rot = k B T /(8πηR 3 ).The scaled decorrelation times 1/(κτ rot ) are mostly larger than 4, significantly larger than unity.This implies that (i) activity has a small effect on the rotational diffusion, and (ii) that the cilia contribute significantly to reduce the rotational diffusion by increasing the effective hydrodynamic radius.A factor 4 reduction is equivalent to a hydrodynamic radius R hydro = 1.6R, which is not implausible as the geometric radius from body center to cilia tips is 4R.The least persistent motion is observed for k ϕ = 1 and χ = 100 • , which indicates a rather strong active contribution to the noise, as might be expected for the strong wiggling motion for k ϕ = 1 (compare movie SM2).
We can now use the expressions (5) to extract the helix radius R h from the data presented in Figs. 5 and  6.The results are shown in Fig. 6(d).They indicate that the helix radius R h is typically rather small, from essentially zero to just a few times the body radius R, with a few exceptions, like k ϕ = 0 and χ = 150 • .The helix pitch is usually a factor 10 larger than the  radius, as typically P h /(2πR h ) = cot(α) 1 (compare Eq. ( 4) and Fig. 6(a)).This is in good agreement with the trajectories displayed in Fig. 7.

The 5-7-5 Swimmer with Oblique Power-Stroke Direction
Except for the oblique propagation direction of the metachronal wave, there are other possibilities to achieve a chirality of the dynamic beat pattern of a ciliated microswimmer.One of these possibilities is to vary the power stroke direction away from the main body axis, and rotate it in the local tangent plane to the left by a tilt angle θ r .For simplicity, we consider in this case only a metachronal wave in the main body direction, i.e. k ϕ = 0.An example for the beating dynamics with phase lag χ = −77 • and tilt angle θ r = 22.5 • is shown in Fig. 9.
Results for the swimming properties with oblique power-stroke direction are displayed in Fig. 10.The results for the swimming velocity, velocity fluctuations, and rotational motion are qualitatively similar as for the oblique metachronal wave.There is again a "sinusoidal" dependence of the velocity v n on the phase lag χ, with a maximum for at χ = −50 • , i.e. for an antiplectic metachronal coordination, see Fig. 10a.However, there are also several pronounced qualitative and quantitative differences.The velocity decreases with increasing θ r , because an increasing fraction of the beat is employed for body rotation rather than forward propulsion.Velocity fluctuations (v n − v n ) 2 peak for synchronously beating cilia for all θ r , as all cilia beat in synchrony in this case (with k ϕ = 0), see Fig. 10b.The body rotation is much more pronounced for all θ r > 0, see Fig. 10c, compared to the case of oblique metachronal wave with θ r = 0. Interestingly, the body rotation becomes very slow for χ 50 • for all θ r > 0. This minimum of spinning frequency Ω n correlates with the minimum of propulsion velocity v n , i.e. weak propulsion is accompanied by slow spinning.Correlation functions for fixed χ = −50 • and various θ r are shown in Fig. 10d.
Results for the fitted parameters of the helical motion and the decay time 1/κ in Eq. ( 6) are displayed in Fig. 11.For most cases, the helix angles α are small, around 10 • to 25 • , which implies that the constant term in Eq. ( 3) dominates over the oscillatory term, and the helices are very thin and elongated.The circling frequencies Ω c , displayed in Fig. 11(b), are generally very small, which implies that cilia beat orientation with θ r > 0 results more in spinning than in circling.Scaled decorrelation times 1/κτ rot are typically larger than 4, which indicates a large hydrodynamic radius, in agreement with the conclusions of Sec.3.1.
We can use again the expression (5) to extract the helix radius R h from the data presented in Figs. 10  and 11.The results are shown in Fig. 11(d).They indicate that similarly as for the case θ r = 0 displayed in Fig. 6(d), the helix radius R h is typically rather small, from nearly zero to just a few times the body radius R. The helix pitch is typically about a factor 10 larger than the radius.

Variation of Cilia Number
We have considered so far a spherical microswimmer with a fixed number N = 5 + 7 + 5 = 17 cilia.It is now of course interesting to see how the swimming behavior, in particular the swimming velocity, depends on the number of cilia.For the swimmer with imposed metachronal wave and θ r = 0, as described in Sec.3.1, the results are shown in Fig. 12(a), both for the maximum and minimum velocity obtained for various phase lags χ.In general, the maximum velocity is found to increase with cilia number, consistent with the results of Ref. [35] for high cilia numbers.For all cilia numbers considered, the highest velocity is obtained for k ϕ = 0, i.e. for a wave direction along the main body axis, in agreement with our arguments in Sec.3.1 above.In contrast, the variation of the minimum velocity on cilia number is much less pronounced.
We want to compare this result with the case of a self-organized metachronal wave -similar to what has been studied in Ref. [17] for planar cilia arrays.The cilia are distributed randomly, but as homogeneously as possible on the body surface.We consider several different distributions, so that the variance of the results for the same cilia number can be taken as a measure for the sensitivity of the swimming velocity on the spatial distribution.Results are shown in Fig. 12(b).Again, the velocity increases with cilia number, but seems to level off at about 40 cilia.This saturation is not too surprising, as with increasing cilia density, the effect of each individual cilium on the propulsion diminishes due to the interaction with the neighbors.For N cilia on the surface, the average distance d between them is approximately determined by d (4πR 2 /N ) 1/2 , which yields d/L 1/6 for N = 40 and cilia length L = 3R.
The results for the self-organized metachronal wave on a spherical body can now be compared with those on a planar substrate as a function of d/L.For the planar case, the fluid transport velocity was found to scale as v f luid ∼ (d/L) −γ with γ 1.4 [17].This implies a dependence of the swimming velocity v swim on cilia number as v swim ∼ N γ/2 , i.e. a behavior somewhere between linear and square root, which seems not inconsistent with the numerical results of Fig. 12(b).The simulation results of Ref. [35] for microswimmers with large cilia numbers (in the range N = 20 to N = 320) also show a sublinear dependence of v swim on N , which a nearly linear dependence for intermediate values of the phase lag χ, and a strongly sublinear dependence for χ = 0.It is also interesting to note that the results for cilia arrays on a planar substrate indicate that d/L = 1/6 is in the regime where transport with metachronal coordination is much more efficient than synchronous beating [17] -in qualitative agreement with the results displayed in Fig. 12(a).

Summary and Conclusion
The dynamics and motion of multi-ciliated microswimmers with a spherical body and a small number N , in the range 5 < N < 50, of long cilia, with length L three times the body radius, has been investigated by mesoscale hydrodynamics simulations.A metachronal wave is imposed for the cilia beat, for which the wave vector has both a longitudinal, k θ , and a latitudinal, k ϕ component.The dynamics and motion is characterized by the swimming velocity v n along the main body axis, the variance of the velocity averaged over a full beat cycle, the spinning velocity Ω n around the main body axis, as well as the parameters of the helical trajectory, which are the circling velocity Ω c , the helix angle α, the helix radius R h and pitch P h .
Our simulation results show that the microswimmer motion strongly depends on the latitudinal wave number k ϕ and the longitudinal phase lag χ = k θ (π/4).We find, not unexpectedly, that the microswimmers usually spin around their own axis, and swim on helical trajectories.It is important to notice that spinning and helical motion are not necessarily correlated, as a spinning particle can move on a perfectly straight trajectory.However, chirality in the metachronal beat pattern generically generates helical trajectories.In most cases, the helices are found to be thin and stretched, i.e. the helix radius R h is about an order of magnitude smaller than the pitch.An interesting result is also that the rotational diffusion of the microswimmer is significantly smaller than the passive rotational diffusion of the body alone.This indicates that active contributions to rotational diffusion are small, and that the extended cilia make a pronounced contribution to the hydrodynamic radius.
The swimming velocity v swim increases with the number N of cilia on the body.Our simulation results indicate a slightly sublinear dependence on N .The comparison with the transport velocity of planar cilia arrays on the cilia separation predicts a dependence v swim ∼ N 0.7 .Our simulation results for self-organized metachronal waves are not inconsistent with such a relation.
Finally, it is interesting to note that already a relatively small number of about ten cilia, beating with a phase lag in the form of a metachronal wave, are sufficient to generate a steady propulsion and smooth swimming motion -in contrast to the strongly oscillatory motion of Chlamydomonas [48,49] with its two cilia beating in synchrony.Such a smooth swimming motion can be of significant benefit for marine microorganisms, because large disturbances can be exploited by predators to locate their prey [50].

Acknowledgments
Support by the Deutsche Forschungsgemeinschaft (DFG) through the priority program on "Microswimmers -from single particle motion to collective behavior" (SPP 1726) is gratefully acknowledged.Computing time has been granted through JARA-HPC on the supercomputer JURECA [51] at Forschungszentrum Jülich.The cilium is considered as a semi-flexible filament, which is modeled by point particles (beads) that are interconnect by springs.Beads are arranged in three linear chains that form a bundle with crosssection of an equilateral triangle (see Fig. A1).Each chain consist of N c = 26 individual beads equally spaced at a distance of 0.5 a, where a is the linear MPC collision-cell size.Springs with equilibrium length l b connect each bead to its neighboring beads along the line.Two diagonal connection on each face provide lateral stability (see Fig. A1).The forces along the active chain are generated by changing the equilibrium length of the springs along the chain [17].

Appendix A. Ciliary Beat
The beat pattern is generated by a heuristic model for the time evolution of bond forces along the chain.In this description, the beat can be controlled by a few key parameters, and allows adaptation to external flows.The beat is regulated by making the equilibrium curvature along the cilium depend on a pivot point i 0 .Thus, the bond length along the active filament varies as where i is the segment number which varies from n 0 to L/l b .The first three beads of each linear chain are passive and stay at the bond-length value l i<n0 = l b .Therefore, we set n 0 = 3.The difference between power and recovery stroke depends on the pivot-point position i 0 .During the power stroke, the pivot point is set to the first point of active beating n 0 .It stays at this position, until the power stroke time τ p has passed.Then, during the recovery stroke, the pivot-point i 0 moves along the cilium with a constant velocity v rec until the recoverystroke time τ r has passed.Together, this results in a beat as displayed in Fig. 2, with a constant beat period τ b = τ p + τ r .
In addition to the global beat dynamic, the dynamics of the beat pattern is limited by the maximum energy E max each motor can inject into the system.If the bond-length change δl i compared to the current distance between the beads exceeds the maximal energy E max a motor can exert, the motor stall and δl i is reduced such that a maximum energy of E max is not exceeded.Due to the thermal noise in the system, the actual distance between two beads fluctuates as well.Although we do allow for "helpful" fluctuations, we prohibit the "obstructive" fluctuations by clipping the bond length at the previous value δl i (t+1) = δl(t).This update scheme ensures step-wise directional movement similar to an active Brownian ratchet.
The parameters of the cilia beat pattern are specified in Tab.A1.

Appendix B. Parameter Extraction from Correlation Functions
Due to the long time scales and the inherent noise of the auto-correlation function, it is difficult to extract parameters reliably.The signal-to-noise ratio is particularly high in the limit where the time scales 1/Ω c and 1/κ approach the total simulation time T s .Due to the approximate rotational symmetry of the microswimmer, we assume that the tangent vector to its center-of-mass trajectory is along the main body axis n.This allows us to access the correlation functions, which determine the trajectory in two different ways.Two examples for the quality of the fit results are shown in Fig. B1.In order to estimate the quality of the fits, we employ the difference between two different -although related -auto-correlation functions which characterize the helical swimming motion.On the one hand, we use the correlation function n(t) n(t + τ ) of the main body axis, as discussed in the main text; on the other hand, the correlate the tangent vectors of the smoothed center-of-mass trajectory with a window of τ .
For a fixed total simulation time T s , the error of the auto-correlation function increases proportional to τ /T s , which is used as a weight factor for the fit.As a second error estimate, we employ the difference between a fit of the full trajectory with τ ≤ T s , and a fit for τ < 0.75 T s .Parameters, where the fit does not converge, or where the difference between the parameter obtained from the full fit and a partial fit exceeds the error threshold displayed in Tab.B1, are omitted and are not shown in the plots.Table B1.Error thresholds for the fits of the auto-correlation functions.
Furthermore, the microswimmer can perform motions which more complex than a simple helix, e.g. a nutation-like motion.This is not captured by our analytic formula, Eq. ( 6).We thus omit all results from fits with a too large mean squared deviation.

Figure 2 .
Figure 2. Sketch of the ciliary beat pattern -The color indicates time.The elongated conformation during the power stroke (orange to green) is followed by the buckled conformation during the recovery stroke (green to blue).

Figure 3 .
Figure 3. Temporal beat pattern of a 4-6-4 swimmer with cilia beat in the longitudinal direction (θr = 0), with kϕ = 0 and phase lag χ = −77 • .The varying cilium color indicates the instantaneous stage in the beat cycle (compare Fig.2).The progression of the beat is indicated by the time t, given in units of the beat period τ b .The motion of the swimmer is from left to right.The translational motion is not to scale.See also movie SM1 for an illustration of the swimming behavior.

Figure 4 .
Figure 4. Temporal beat pattern of a 4-6-4 swimmer with cilia beat in the longitudinal direction (θr = 0), with kϕ = 1 and phase lag χ = −77 • .The varying cilium color indicates the instantaneous stage in the beat cycle (compare Fig.2).The progression of the beat is indicated by the time t, given in units of the beat period τ b .The motion of the swimmer is from left to right.The translational motion is not to scale.See also movie SM2 for an illustration of the swimming behavior.

Figure 6 .
Figure 6.Swimming properties of a 5-7-5 swimmer with cilia beat in the longitudinal direction (θr = 0), for various phase lags χ and latitudinal modes kϕ.a) Alignment angle α; b) Rotation frequency Ωc, normalized by the beat period τ ; c) Inverse correlation time κ, normalized by the rotational diffusion time τrot = 1/(2Drot).d) Helix radius R h in units of the body radius R, as obtained from Eq. (5).

Figure 9 . 45 Figure 10 .Figure 11 .
Figure 9. Temporal beat pattern of a 4-6-4 swimmer with cilia beat in an oblique direction (θr = 22.5 • ), with latitudinal wave vector kϕ = 0 and phase lag χ = −77 • .The varying cilium color indicates the instantaneous stage in the beat cycle (compare Fig.2).The progression of the beat is indicated by the time t, given in units of the beat period τ .The motion of the swimmer is from left to right.The translational motion is not to scale.See also movie SM3 for an illustration of the swimming behavior.

Figure 12 .
Figure 12.Swimming velocity for microswimmers with various numbers of cilia.(a) With a predefined metachronal wave with various latitudinal wave numbers kϕ, as indicated.Full (dashed) lines show that maximum (minimum) velocity obtained by variation of the phase lag χ.(b) With self-organized metachronal waves and random spatial arrangement.The red line is a spline fit to different cilia arrangements.See also movie SM4 for an illustration of the beat and swimming behavior.

1 k
b T 13 a 20'000 k b T /a 2 0.19 l b /h TableA1.Model parameters (given in MPC units of thermal energy k B T , collision cell size a, and collision time h) used to generate the beat pattern of the cilium in the simulations.

Figure B1 .
Figure B1.Fit results of the temporal auto-correlation function which characterizes the helical trajectory.Either the autocorrelation function of the unit normal vector n of the main body axis is employed, or of the tangent vector along the smoothed center-of-mass trajectory, both with a temporal delay τ .Simulation data are shown as solid lines, fits to Eq. (6) as dashed lines.