Finding the ciliary beating pattern with optimal efficiency

Finding the ciliary beating pattern with optimal efficiency

Natan Osterman J. Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia Department of Physics, Ludwig-Maximilians University Munich, Amalienstrasse 54, 80799 Munich, Germany    Andrej Vilfan J. Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 29, 1000 Ljubljana, Slovenia

We introduce a measure for energetic efficiency of biological cilia acting individually or collectively and numerically determine the optimal beating patterns according to this criterion. Maximizing the efficiency of a single cilium leads to curly, often symmetric and somewhat counterintuitive patterns. But when looking at a densely ciliated surface, the optimal patterns become remarkably similar to what is observed in microorganisms like Paramecium. The optimal beating pattern then consists of a fast effective stroke and a slow sweeping recovery stroke. Metachronal coordination is essential for efficient pumping and the highest efficiency is achieved with antiplectic waves. Efficiency also increases with an increasing density of cilia up to the point where crowding becomes a problem. We finally relate the pumping efficiency of cilia to the swimming efficiency of a spherical microorganism and show that the experimentally estimated efficiency of Paramecium is surprisingly close to the theoretically possible optimum.

Many biological systems have evolved to work with a very high energetic efficiency. For example, muscle can convert the free energy of ATP hydrolysis to mechanical work with more than 50% efficiency Kushmerick and Davies (1969), the F1-F0 ATP synthase converts electrochemical energy of protons to chemical energy stored in ATP molecules with even higher efficiency Yoshida et al. (2001), etc. At first glance, the beating of cilia and flagella does not fall into the category of processes with such a high efficiency. Cilia are hair-like protrusions that beat in an asymmetric fashion in order to pump the fluid in the direction of their effective stroke Sleigh (1974). They propel certain protozoa, such as Paramecium, and also fulfill a number of functions in mammals, including mucous clearance from airways, L-R asymmetry determination and transport of an egg cell in Fallopian tubes. Lighthill Lighthill (1952) defines the efficiency of a swimming microorganism as the power that would be needed to drag an object of the same size with the same speed through viscous fluid, divided by the actually dissipated power. Although the efficiency defined in this way could theoretically even exceed 100% Michelin and Lauga (2010), the actual swimming efficiencies are of the order of Purcell (1997); Chattopadhyay et al. (2006). In his legendary paper on life at low Reynolds number Purcell (1977) Purcell stated that swimming microorganisms have a poor efficiency, but that the energy expenditure for swimming is so small that it is of no relevance for them (he uses the analogy of “driving a Datsun [a fuel-efficient car of the period] in Saudi Arabia”). Nevertheless, later studies show that swimming efficiency is important in microorganisms. In Paramecium, more than half of the total energy consumption is needed for ciliary propulsion Katsu-Kimura et al. (2009).

When applied to ciliary propulsion, Lighthill’s efficiency Lighthill (1952) has some drawbacks. For one, it is not a direct criterion for the hydrodynamic efficiency of cilia as it also depends on the size and shape of the whole swimmer. Besides that it is, naturally, only applicable for swimmers and not for other systems involving ciliary fluid transport with a variety of functions, like L-R asymmetry determination Supatto et al. (2008). We therefore propose a different criterion for efficiency at the level of a single cilium or a carpet of cilia. A first thought might be to define it as the volume flow rate of the transported fluid, divided by the dissipated power. However, as the flow rate scales linearly with the velocity, but the dissipation quadratically, this criterion would yield the highest efficiency for infinitesimally slow cilia, just like optimizing the fuel consumption of a road vehicle alone might lead to fitting it with an infinitesimally weak engine. Instead, like engineers try to optimize the fuel consumption at a given speed, the well-posed question is which beating pattern of a cilium will achieve a certain flow rate with the smallest possible dissipation.

The problem of finding the optimal strokes of hypothetical microswimmers has drawn a lot of attention in recent years. Problems that have been solved include the optimal stroke pattern of Purcell’s three link swimmer Tam and Hosoi (2007), an ideal elastic flagellum Spagnolie and Lauga (2010), a shape-changing body Avron et al. (2004), a two- and a three-sphere swimmer Alouges et al. (2009) and a spherical squirmer Michelin and Lauga (2010). Most recently, Tam and Hosoi optimized the stroke patterns of Chlamydomonas’ flagella Tam and Hosoi (2011). But all these studies are still far from the complexity of a ciliary beat with an arbitrary three-dimensional shape, let alone from an infinite field of interacting cilia. In addition, they were all performed for the swimming efficiency of the whole microorganism, while our goal is to optimize the pumping efficiency at the level of a single cilium, which can be applicable to a much greater variety of ciliary systems.

So we propose a cilium embedded in an infinite plane (at ) and pumping fluid in the direction of the positive -axis. We define the volume flow rate as the average flux through a half-plane perpendicular to the direction of pumping Smith et al. (2008). With we denote the average power with which the cilium acts on the fluid, which is identical to the total dissipated power in the fluid filled half-space. We then define the efficiency in a way that is independent of the beating frequency as


As we show in Appendix 1, minimizing the dissipated power for a constant volume flow rate is equivalent to maximizing at a constant frequency. A similar argument for swimming efficiency has already been brought forward by Avron et al. Avron et al. (2004).

Furthermore, a general consequence of low Reynolds number hydrodynamics is that the volume flow only depends on the shape of the stroke and on the frequency, but not on the actual time dependence of the motion within a cycle. This is the basis of Purcell’s scallop theorem Purcell (1977). As a consequence, the optimum stroke always has a dissipation rate constant in time. We show this in Appendix 2.

We can make the efficiency completely dimensionless if we factor out the effects of the ciliary length , the beating frequency and the fluid viscosity . The velocity with which a point on the cilium moves scales with and the linear force density (force per unit length) with . The total dissipated power , obtained by integration of the product of the velocity and linear force density over the length, then scales with . The volume flow rate scales with . Finally, the efficiency scales with . The dimensionless efficiency can therefore be defined as


When optimizing the efficiency of ciliary carpets, we have to use the measures of volume flow and dissipation per unit area, rather than per cilium. We introduce the surface density of cilia , which is on a square lattice. In the following we show that the volume flow generated per unit area, , is also equivalent to the flow velocity above the ciliary layer. The fluid velocity above an infinite ciliated surface namely becomes homogeneous at a distance sufficiently larger than the ciliary length and metachronal wavelength. The far field of the flow induced by a single cilium located at the origin and pumping fluid in direction of the axis has the form Vilfan and Jülicher (2006)


with an arbitrary amplitude . For this field the volume flow rate is


and the velocity above an infinite field of such cilia is


which is independent of . In this regime, one can simplify the description of cilia by replacing them with a surface slip term with velocity Jülicher and Prost (2009).

We now define the collective efficiency as and in dimensionless form as


is a function of the beat shape, the dimensionless density and the metachronal coordination, which will be explained later. Additionally, for a single cilium or for collective cilia the efficiency also depends on the dimensionless radius of the cilium, , but this dependence is rather weak, of logarithmic order.

At this point we note that our definition of efficiency is different from that used by Gueron and Levit-Gurevich Gueron and Levit-Gurevich (1999). They define efficiency as volume flux through a specifically chosen rectangle above the group of cilia divided by the dissipated power. While this measure is useful for studying the effect of coupling and metachronal coordination (they show that the collective efficiency of a group of cilia increases with its size), it lacks the scale invariance discussed above. Gauger et al. Gauger et al. (2009) studied a model for individual and collective magnetically driven artificial cilia. Rather than introducing a single measure for the efficiency, they studied the pumping performance (which is the more relevant quantity in artificial systems) and dissipation separately. They showed that the pumping performance per cilium can be improved with the proper choice of the metachronal wave vector, while the dissipation per cilium remains largely constant. Both studies were limited to two-dimensional geometry (planar cilia arranged in a linear row) and neither of them uses a scale-invariant efficiency criterion proposed here. On the other hand, Lighthill’s criterion for swimming organisms shares the same scaling properties as ours (it scales with the square of the swimming velocity, divided by dissipation), but differs in definition because it measures the swimming and not the pumping efficiency. At the end we will show how the two measures are related to each other for a spherical swimmer.

Our goal is to find the beating patterns that have the highest possible efficiency for a single cilium, as well as the beating pattern, combined with the density and the wave vector that give the highest efficiency of a ciliated surface.

I The model

We describe the cilium as a chain of touching beads with radii . The first bead of a cilium has the center position , and each next bead in the chain is located at . The maximum curvature of the cilium is limited by the condition


Naturally, beads cannot overlap with the surface () or with each other .

We describe the hydrodynamics using the mobility matrix formalism. If the force acting on bead is , the resulting velocities are


In this formalism, each element is itself a matrix, corresponding to 3 spatial dimensions. In general, the above equation should also include angular velocities and torques, but they are negligible for small beads when the surface speeds due to rotational motion are much smaller than those due to translational motion. The mobility matrix is symmetric and positive-definite Happel and Brenner (1983). Therefore, one can always invert it to determine the friction matrix , which determines the forces on particles moving with known velocities


If the particles were at large distances relative to their sizes, the elements of the mobility matrix would be determined by Blake’s tensor Blake (1971), which represents the Green function of the Stokes flow in the presence of a no-slip boundary. In our case the condition of large interparticle distances is not fulfilled and we use the next higher approximation, which is the Rotne-Prager tensor in the presence of a boundary, as described in a previous paper Vilfan et al. (2010).

The volume flow rate in direction, averaged over one beat period , depends on -components of forces acting on particles and their heights above the boundary Smith et al. (2008):


The dissipation rate is simply the total power needed to move the beads against viscous drag,


We numerically maximized the quantity for a set of angles and different numbers of beads. We used the sequential quadratic programming algorithm (SQP) from NAG numerical libraries (Numerical Algorithms Group). The full details of the numerical procedure are given in Appendix 3.

To study the collective efficiency and metachronal coordination, we studied an array of cilia (unit cell) on a square lattice with a lattice constant . We introduced periodic boundary conditions by adding hydrodynamic interactions between particles and the representations beyond lattice boundaries. So if a certain element in the mobility matrix describing interaction between particles at and is , we replace it by . Here denotes the size of the unit cell. For the sake of numerical efficiency, we used the full Rotne-Prager form for the first instances () and approximated the interaction with its long range limit, independent of the actual particle positions, for the rest (SI).

We expect the optimal solution to have the form of metachronal waves with a wave vector . In order to satisfy the periodic boundary conditions, and have to be integer numbers, e.g., between and .

Figure 1: Optimal trajectories of the 1-particle model. A) Idealized case of a small particle restricted to . The solution consists of piecewise circular arcs, determined by geometric parameters and . B) Numerical solutions for finite-sized particles, plotted for different ratios . C) Optimal path for a particle at a constant distance from the origin, , with . The transparent hemisphere symbolizes the surface on which the particle can move. D) Dimensionless efficiency as a function of the dimensionless particle radius . The black line shows the model with variable distance and the red line with a fixed distance from the origin. The dashed line shows the limit of small radii (), .

Ii Results

ii.1 Single-particle model

We first start with some simple models that are not necessarily feasible in practice, but allow important insight into how the optimum is achieved. We will follow the spirit of the model used to study the synchronization of cilia Vilfan and Jülicher (2006), where we replace the cilium by a small spherical particle. There are many swimmer models building on similar assumptions, for example the three-sphere-swimmer Najafi and Golestanian (2004), and they all have in common that they assume the connections between spheres to be very thin and neglect any hydrodynamic forces acting on them.

So the first hypothetical model we study is a single sphere of radius that can move along an arbitrary path in the half space above the boundary, but in order to mimic the tip of a cilium it has to stay within the distance of the origin, . In order to simplify the calculation we also assume that the sphere is small, . In this limit, we can neglect the effect of the boundary on the hydrodynamic drag, which is then always . The dissipated power is then simply . Because it has to be constant in time, we can also write it as


with denoting the total distance traveled within one cycle and its period. The average volume flow follows from Eq. (10) as


where is the area of the particle trajectory, projected onto the plane. The resulting efficiency is (1)


To find the optimal path, we thus have to maximize the area-to-circumference ratio of the path, while fulfilling the constraints and . Obviously, there is no benefit in going out of the plane, but there is cost associated with it. Therefore, the optimum trajectories will be planar. As any curve that minimizes its circumference at a fixed surface area, the unconstrained segments of the trajectory have to be circle arcs. The curve has the shape shown in Fig. 1A. A numerical solution shows that the area-to-circumference ratio is maximal if the angle defined in Fig. 1A has the value . The resulting maximal efficiency in the limit is , or, in dimensionless form, .

Solutions for finite values of are shown in Fig. 1B and their efficiencies in Fig. 1D. The highest possible numerical efficiency of this model is , which is achieved at .

Another version of the single-particle model is one in which the particle has to maintain a fixed distance () from the origin, while it is free to move along the surface of a sphere (Fig. 1C). This is an additional constraint and can therefore only reduce the achievable efficiency. As shown by the red line in Fig. 1D, the efficiency indeed lies somewhat below that of the model with a variable distance and reaches a maximum value of .

ii.2 N particles, stiff cilium

The next minimalistic model we will study is a stiff cilium: a straight chain of beads with radius and a total length of that can rotate freely around the center of the first bead. The problem is related to artificial cilia driven by a magnetic field Vilfan et al. (2010); Gauger et al. (2009); Downton and Stark (2009) in which the orientation of the cilium largely (although not completely) follows the direction of the magnetic field. A related optimization has been performed by Smith et al. Smith et al. (2008), but with two important differences. First, Smith et al. optimize the volume flow alone and not the efficiency. Their optimal stroke therefore touches the surface during the recovery stroke, while ours has to keep some distance in order to limit the dissipation. Second, they restrict themselves to cilia beating along tilted cones, whereas we allow any arbitrary pattern.

The motion of a stiff cilium on its optimal trajectory is shown in Figure 2A. The path of its tip closely resembles that of a single sphere at a fixed radius. The resulting dimensionless efficiency for beads is .


[ poster=Fig2A.jpg, label=figure2a.u3d, 3Daac=50.000000, 3Droll=0.000000, 3Dc2c=0 -3 5, 3Droo=300.092304, 3Dcoo=0.0 0.0 0.0, 3Dlights=CAD, ]6cm5cmFig2A.u3d


[ poster=Fig2B.jpg, label=figure2b.u3d, 3Daac=50.000000, 3Droll=0.000000, 3Dc2c=0 -3 5, 3Droo=300.092304, 3Dcoo=0.0 0.0 0.0, 3Dlights=CAD, ]6cm5cmFig2B.u3d


[ poster=Fig2C.jpg, label=figure2c.u3d, 3Daac=50.000000, 3Droll=0.000000, 3Dc2c=0 -3 5, 3Droo=300.092304, 3Dcoo=0.0 0.0 0.0, 3Dlights=CAD, ]6cm5cmFig2C.u3d

Figure 2: Optimal beating patterns of a cilium consisting of particles with different allowed bending angles . The gray surface shows the projection of the tip trajectory on the plane. A) A stiff cilium, . B) A flexible cilium, . The optimal stroke is symmetric in direction. C) Flexible cilium, . The symmetry in direction is broken. D) Dimensionless efficiency as a function of . The black line shows the optimal symmetric solution and the red line the asymmetric solution in cases where it is more efficient. The dashed line shows the maximum efficiency for and double (corresponding to the same curvature).

ii.3 N particles, flexible cilium

As the next level of complexity, and at the same time the first realistic description of biological cilia, we now study a flexible cilium consisting of spherical particles (we use and ). The bending angle per particle is restricted to (7). Such a constraint is necessary for two reasons. For one, the curvature of a biological cilium is restricted by the bending rigidity of the axoneme. In addition, our -particle model veritably represents a continuous cilium only if the curvature radius is sufficiently larger than the size of a sphere. Examples of beating patterns obtained by numerical optimization are shown in Figure 2B,C. Figure 2D shows the dimensionless efficiency as a function of .

It is instructive to look at fundamental symmetries of the problem at this point. First, as any of the problems studied here, it is symmetric upon reflection . For every clockwise beat, there is an equivalent counterclockwise beat with the same efficiency. All cycles discussed here spontaneously break the symmetry. In addition, the equations are invariant upon reflection with simultaneous time reversal, . This symmetry can be broken or not at the efficiency maximum. Interestingly, this depends on the allowed bending between adjacent elements . For example, the solution shown in Fig. 2B is -symmetric, while the one in Fig. 2C is not.

Figure 3: Optimal solutions at fixed wave vectors for interciliary distance , , and . A-F) Optimal solutions for wave vectors (A), (B), (C), (D), (E) and (F). The blue arrow ( axis) denotes the direction of pumping and the red arrow the wavelength and direction of metachronal wave propagation. G) Efficiency (red color represents high efficiency) as a function of the wave vector . The maximum efficiency is in this case achieved for and antiplectic waves are generally more efficient than symplectic. The synchronous solution represents the global minimum of efficiency.

ii.4 Multiple cilia and metachronal waves

We solve the optimization problem of cilia () with periodic boundary conditions by imposing a wave vector , finding the optimal solution for that vector and repeating the procedure for wave vectors. Examples of optimal solutions for 6 different wave vectors are shown in Fig. 3A-F. The efficiency as a function of the wave vector is shown in Fig. 3G. Note that all solutions determined in this section are for counterclockwise beating (viewed from above). For clockwise strokes the component of the wave vector would have the opposite sign. Optimal solutions for four different values of the interciliary distance are shown in Figure 4A-D and the optimal efficiency as a function of in Figure 4E.

Figure 4: Optimal solutions for various interciliary distances: (A), (B), (C) and (D). For reasons of clarity the front rows of cilia are omitted and instead of individual spheres used in the calculation a tube connecting them is shown. E) Highest efficiency as a function of the interciliary distance for cilia consisting of (circles) and (squares) spheres. For , we set in order to allow the same maximum curvature.

Figure 3G shows that the efficiency depends more strongly on the longitudinal () component of the wave vector than on the lateral (). This could partly be due to the nature of the hydrodynamic interaction, which is stronger in longitudinal direction, and partly because the cilia exert larger motion in longitudinal direction and therefore come closer to their neighbors along the axis. Antiplectic metachronal waves (waves propagate in the opposite direction from the fluid pumping) generally have a higher efficiency than symplectic, but the fine structure is much more complex. For low ciliary densities, the optimal solution is found for waves propagating along the axis. However, for higher densities solutions with a positive are more efficient, which means that the waves are dexio-antiplectic. Efficiency also grows with increasing density. But when the interciliary distance reaches crowding becomes a problem and the efficiency drops again. At even higher densities the solution becomes increasingly difficult to find because of the complicated topology of densely packed cilia. Another problem is that the wavelength of the optimal solution, relative to the lattice constant, becomes increasingly long at high densities, which would require a unit cell larger than used in our calculations.

ii.5 Swimming efficiency of a ciliated microorganism

We can finally use these results to estimate the maximum possible swimming efficiency of a ciliated microorganism. For the sake of simplicity, we assume that the swimmer has a spherical shape with radius . According to Lighthill’s definition, the swimming efficiency is defined as


where is the velocity and the total dissipated power Lighthill (1952). Assuming that the layer of cilia is thin in comparison with the size of the organism (), the swimming velocity can be calculated as Stone and Samuel (1996); Jülicher and Prost (2009)


Here is the propulsion velocity above the ciliated layer. The dissipation can be expressed as


In the second equality we used the definition (6), as well as the relationship between and (5). In order to obtain the maximum at a given , the two functional derivatives should be related through a Lagrange multiplier


which is fulfilled if the angular dependence has the form . This leads to




Together, these equations give Lighthill’s efficiency


With a maximum and a typical ratio , we obtain .

For comparison, an optimized envelope model yields an efficiency around with the same parameters (if we translate the ciliary length into maximal displacement of the envelope) Michelin and Lauga (2010), which shows that the latter is a relatively good approximation for cilia if they operate close to the optimal way.

Iii Discussion

We introduced a scale invariant efficiency measure for the fluid pumping by cilia and started with optimizing three simple instructive systems: a free sphere, a sphere at a constant distance from the origin and a stiff cilium. Of those three, the free sphere can reach the highest efficiency. But they are all topped by a flexible cilium. The thickness of a cilium has only a small effect on its efficiency, which strengthens our choice to describe the cilium as a chain of beads. Depending on the allowed bending, the flexible cilium can have different shapes of the optimal beating pattern. In most cases the cilium curls up during the recovery stroke, rather than sweeping along the surface. Such shapes appear “unnatural” if we compare them with those observed in microorganisms Brennen and Winet (1977).

But the collective optimization of ciliary carpets leads to beating patterns that are strikingly similar to what is observed in many ciliated microorganisms. Unlike isolated cilia, they contain a recovery stroke during which they sweep along the surface. This is primarily due to the fact that beating patterns that are optimal for a single cilium (e.g., as shown in Fig. 2B) are not possible on a dense grid due to steric hindrance. The sweeping recovery stroke, on the other hand, allows dense stacking of cilia (best seen in Fig. 4D) which further reduces drag as well as backward flow. The optimal effective stroke becomes significantly faster than the recovery stroke. While a single cilium reaches its highest efficiency if the effective stroke takes about of the cycle, the optimum is around for densely packed cilia. A similar ratio has been observed in Paramecium Gueron et al. (1997). The distance between adjacent cilia in Paramecium is between and Sleigh (1969), consistent with the predicted optimum around . It is interesting that the efficiency of any other wave vector is higher than the efficiency of the synchronous solution (wave vector 0). This is in agreement with some previous simpler, one dimensional models Gauger et al. (2009), but has not yet been shown on a 2-D lattice. We also find that antiplectic waves are generally more efficient than symplectic, although symplectic solutions with a relatively high efficiency exist, too. For high densities and cilia beating counterclockwise, the waves become almost dexioplectic (meaning that the effective stroke points to the right of the wave propagation) and the wavelength becomes similar to the cilium length - both findings are in agreement with observations on Paramecium Machemer (1972); Sleigh (1974). For cilia beating clockwise, laeoplectic waves would be more efficient, which is indeed observed Gheber and Priel (1990). Although the effect of thickness is small, it is interesting to note that thicker cilia have a slightly higher efficiency when isolated or at low surface densities, but are outperformed by thinner cilia at high densities.

The total energetic efficiency of swimming in Paramecium has been measured as Katsu-Kimura et al. (2009). This figure includes losses in metabolism and force generation – the hydrodynamic swimming efficiency alone has been estimated as . This comes close (by a factor of 2) to our result for the maximally possible Lighthill efficiency of a spherical ciliated swimmer . A biflagellate swimmer like Chlamydomonas has a lower theoretical efficiency of Tam and Hosoi (2011), but it is still within the same order of magnitude.

Although efficiencies below seem low, we have shown that Paramecium still works remarkably close to the maximum efficiency that can be achieved with its length of cilia. While longer cilia might have a higher swimming efficiency, there are other considerations that are not included in this purely hydrodynamic study. For example, the bending moments and the power output per ciliary length can be limiting Sleigh and Blake (1977). Thus, our study shows that at least for ciliates like Paramecium, Purcell’s view that efficiency is irrelevant for ciliary propulsion has to be revisited. Efficiency of swimming does matter for them, and in their own world they have well evolved to swim remarkably close to the optimal way.

We have benefited from fruitful discussions with Frank Jülicher. This work was supported by the Slovenian Office of Science (Grants P1-0099, P1-0192, J1-2209, J1-2200 and J1-0908).

Appendix A Appendix 1. Justification of the efficiency criterion

In this section we prove that instead of minimizing the dissipated power for a constant volume flow rate we can maximize the numerical efficiency at a constant beating frequency . Let the volume flow rate and dissipation be functionals of the trajectory shape and functions of . The efficiency , defined as , is only a functional of , but independent of . The relationship


proves immediately that minimizing the dissipated power while keeping the volume flow constant is equivalent to maximizing at a constant .

Appendix B Appendix 2. Optimal solutions have a constant dissipation

In the following we demonstrate that a trajectory with optimal efficiency always has a dissipation that is constant in time. The pumping performance only depends on the shape of the trajectory and the period , but not on the velocity along the trajectory. We have to determine the latter in a way that minimizes the average dissipation. We parameterize the trajectory as where the phase is a function of time that fulfills and . The dissipation at any moment is quadratic in velocity, , and can be written as


We obtain the average dissipation by integrating over one period,


with and . Minimization of the average dissipation while keeping the period constant leads to the following Euler-Lagrange equation (note the swapped roles of and )


This proves that the optimal solution has a dissipation that is constant in time.

Appendix C Appendix 3. Numerical optimization procedure

c.1 Single cilium

For a single cilium, modeled as a chain of beads, we parameterize the stroke as a sequence of equally spaced time steps with duration . Let the vector represent the position of the bead at time step . Because the trajectory is periodic in time, we have . For a given stroke the volume flow rate is calculated according to Eq. 10 as


denotes the viscosity of the surrounding fluid and the height of -th bead above the surface. The dissipated power, defined in Eq. 11, can be written as


The dimensionless efficiency, which we want to maximize, follows as with (2). The friction matrix that appears in above expressions is obtained by numerical inversion of the mobility matrix , which is calculated using the Rotne-Prager approximation in the presence of a boundary, as described in Vilfan et al. (2010).

We computed the optimal strokes using the NAG (The Numerical Algorithms Group Ltd., Oxford, UK) Fortran Library E04UGF routine, which maximizes an arbitrary smooth function subject to constraints using a sequential quadratic programming method. At each time step we parameterize the bead positions by a set of pairs of polar () and azimuthal () angles, as described in the main text (Section “The Model”). Therefore, the parameter space in which we search for the for the efficiency maximum is dimensional. For most calculations presented here (, ), this means 3192 dimensions. We still consider optimization in such high-dimensional space preferable to parameterizing the beat shape in terms of Fourier modes (see, e.g. Tam and Hosoi (2011)), which would have less variables, but therefore a more complex landscape with more local maxima.

The stroke shapes are subject to two types of constraints: (i) beads are modeled as hard core objects that allow no overlapping of two beads or a bead with the surface, (ii) the maximum bending angle defined between three adjacent beads is limited to (7).

We initialized the optimization procedure with a simple beating pattern (a tilted cone) that had the correct handedness (clockwise) and checked that the result was otherwise independent of the initial state.

c.2 Carpets of cilia

For an array of cilia the size of parameter space would grow with the square of the system size and the demand to calculate the hydrodynamic mobility matrix with its fourth power. The problem would become numerically unsolvable even for a relatively small field of cilia. However, in an infinite system (or a system with periodic boundary conditions) we expect the optimal solution to have the form of metachronal waves with an unknown wave vector. The optimization problem at a fixed wave vector then requires the same number of parameters as the optimization of a single ciliary beat. The globally optimal solution can be found by repeating the calculation for all wave vectors that fulfill the periodic boundary conditions.

For a wave vector , the -th sphere of the cilium at lattice position follows the path


where is the path of the cilium at lattice position , which we are optimizing. Similarly, the force acting on the same bead has the following time dependence


The equations of motion of -th bead belonging to the cilium read


In this expression represents the element of the mobility matrix that links the response of -th bead of the cilium to the force acting on -th bead of the cilium at time .

Figure S1: The elements of the mobility matrix describe the response of the particles forming the central cilium () to forces acting on all other particles. Cilia are positioned on a square lattice with lattice constant , while the unit cell comprising cilia has the edge length . To facilitate the calculation we use the full Rotne-Prager form for unit cells and use the long-range limit expressions for more distant cilia.

Periodic boundary conditions are fulfilled if and are multiples of , with denoting the size of the field (Fig. S1). We can therefore use integer wave vectors instead. We can make use of the periodicity by carrying out the summation in the above equation




The indices and run over one unit cell, centered around the origin. For odd , this would be from to . For an even , they run from to , but the boundary terms are given half weight each. This distribution is necessary in order to preserve the symmetry of the mobility/friction matrix.

For the sake of numerical efficiency, we use the full Rotne-Prager form for the first instances () and approximate the interaction with


for larger indices. We used for the optimization and checked the final efficiency with (the correction was not significant). The contribution of distant cilia towards can easily be calculated numerically to quadratic order in and and has the form


Numerical values of the coefficients are listed in table S1.

0 4.51681 14.0613 -2.60821 -12.8518
1 1.80961 0.680011 0.182076 -0.210573
2 1.1135 0.161671 0.0474549 -0.0445071
3 0.801452 0.0608975 0.0181854 -0.0163511
Table S1: Coefficients describing the contributions of cilia beyond periods, for different orders .

With discretized time steps the expression for velocity (S10) becomes


with denoting the generalized mobility matrix that now couples the forces and velocities at different time steps. The brackets denote a function which is if the condition is fulfilled and otherwise. We also set the number of time steps to be a multiple of the lattice size, . is in total a dimensional matrix, which can be re-ordered into a block-diagonal form with blocks for greater numerical efficiency.

We calculate corresponding generalized friction matrix by inverting the generalized mobility matrix,


The friction matrix now allows us to calculate the forces if we know the velocities of all beads at all times,


Within the finite difference approximation this leads to equations


Now we can numerically optimize the quantity as a function of the angles that parameterize the beat shape (in total variables). We initialize the system in a way that only solutions with clockwise rotation (as seen from above) are considered. The constraints are similar to those on a single cilium, except that hard-core repulsion between neighboring cilia has to be taken into account, too.

In conclusion, the combination of an efficient description of hydrodynamics, the usage of periodic boundary conditions and the translational invariance of the metachronal wave allowed us to reduce the optimization problem to one with the same number of variables as for a single cilium. The computational demand to calculate the efficiency for a certain set of variables scales as , rather than the fourth power which would be the case if all cilia were treated independently. The optimization for one wave vector typically takes a few days on one core of the Intel Xeon E5520 CPU. To find the global maximum we have to repeat the calculation with all wave vectors, in total.


  • Kushmerick and Davies (1969) M. J. Kushmerick and R. E. Davies, Proc. R. Soc. Lond. B Biol. Sci. 174, 315 (1969).
  • Yoshida et al. (2001) M. Yoshida, E. Muneyuki, and T. Hisabori, Nat. Rev. Mol. Cell. Biol. 2, 669 (2001).
  • Sleigh (1974) M. A. Sleigh, ed., Cilia and Flagella (Academic Press, London, 1974).
  • Lighthill (1952) M. J. Lighthill, Comm. Pure Appl. Math. 5, 109 (1952).
  • Michelin and Lauga (2010) S. Michelin and E. Lauga, Phys. Fluids 22, 111901 (2010).
  • Purcell (1997) E. M. Purcell, Proc. Natl. Acad. Sci. USA 94, 11307 (1997).
  • Chattopadhyay et al. (2006) S. Chattopadhyay, R. Moldovan, C. Yeung, and X. L. Wu, Proc. Natl. Acad. Sci. USA 103, 13712 (2006).
  • Purcell (1977) E. M. Purcell, Am. J. Phys. 45, 3 (1977).
  • Katsu-Kimura et al. (2009) Y. Katsu-Kimura, F. Nakaya, S. A. Baba, and Y. Mogami, J. Exp. Biol. 212, 1819 (2009).
  • Supatto et al. (2008) W. Supatto, S. E. Fraser, and J. Vermot, Biophys. J. 95, L29 (2008).
  • Tam and Hosoi (2007) D. Tam and A. E. Hosoi, Phys. Rev. Lett. 98, 068105 (2007).
  • Spagnolie and Lauga (2010) S. E. Spagnolie and E. Lauga, Phys. Fluids 22, 031901 (2010).
  • Avron et al. (2004) J. E. Avron, O. Gat, and O. Kenneth, Phys. Rev. Lett. 93, 186001 (2004).
  • Alouges et al. (2009) F. Alouges, A. DeSimone, and A. Lefebvre, Eur. Phys. J. E Soft Matter 28, 279 (2009).
  • Tam and Hosoi (2011) D. Tam and A. E. Hosoi, Proc. Natl. Acad. Sci. USA 108, 1001 (2011).
  • Smith et al. (2008) D. J. Smith, J. R. Blake, and E. A. Gaffney, J. R. Soc. Interface 5, 567 (2008).
  • Vilfan and Jülicher (2006) A. Vilfan and F. Jülicher, Phys. Rev. Lett. 96, 058102 (2006).
  • Jülicher and Prost (2009) F. Jülicher and J. Prost, Eur. Phys. J. E Soft Matter 29, 27 (2009).
  • Gueron and Levit-Gurevich (1999) S. Gueron and K. Levit-Gurevich, Proc. Natl. Acad. Sci. USA 96, 12240 (1999).
  • Gauger et al. (2009) E. M. Gauger, M. T. Downton, and H. Stark, Eur. Phys. J. E Soft Matter 28, 231 (2009).
  • Happel and Brenner (1983) J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Kluwer, Dodrecht, 1983).
  • Blake (1971) J. R. Blake, Proc. Camb. Phil. Soc. 70, 303 (1971).
  • Vilfan et al. (2010) M. Vilfan, A. Potočnik, B. Kavčič, N. Osterman, I. Poberaj, A. Vilfan, and D. Babič, Proc. Natl. Acad. Sci. USA 107, 1844 (2010).
  • Najafi and Golestanian (2004) A. Najafi and R. Golestanian, Phys. Rev. E 69, 062901 (2004).
  • Downton and Stark (2009) M. T. Downton and H. Stark, Europhys. Lett. 85, 44002 (2009).
  • Stone and Samuel (1996) H. A. Stone and A. D. Samuel, Phys. Rev. Lett. 77, 4102 (1996).
  • Brennen and Winet (1977) C. Brennen and H. Winet, Ann. Rev. Fluid Mech. 9, 339 (1977).
  • Gueron et al. (1997) S. Gueron, K. Levit-Gurevich, N. Liron, and J. J. Blum, Proc. Natl. Acad. Sci. USA 94, 6001 (1997).
  • Sleigh (1969) M. A. Sleigh, Int. Rev. Cytol. 25, 31 (1969).
  • Machemer (1972) H. Machemer, J. Exp. Biol. 57, 239 (1972).
  • Gheber and Priel (1990) L. Gheber and Z. Priel, Cell Motil. Cytoskeleton 16, 167 (1990).
  • Sleigh and Blake (1977) M. A. Sleigh and J. R. Blake, in Scale Effects in Animal Locomotion, edited by T. J. Pedley (Academic Press, London, 1977), pp. 243–256.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description