## Abstract

Cilia-driven flow occurs in the airway surface liquid, in the female and male reproductive tracts and enables symmetry-breaking in the embryonic node. Viscoelastic rheology is found in healthy states in some systems, whereas in others may characterise disease, motivating the development of mathematical models that take this effect into account. We derive the fundamental solution for linear viscoelastic flow, which is subsequently used as a basis for slender-body theory. Our numerical algorithm allows efficient computation of three-dimensional time-dependent flow, bending moments, power and particle transport. We apply the model to the large-amplitude motion of a single cilium in a linear Maxwell liquid. A relatively short relaxation time of just 0.032 times the beat period significantly reduces forces, bending moments, power and particle transport, the last variable exhibiting exponential decay with relaxation time. A test particle is propelled approximately one-fifth as quickly along the direction of cilia beating for scaled relaxation time 0.032 as in the Newtonian case, and mean volume flow is abolished, emphasizing the sensitivity of cilia function to fluid rheology. These results may have implications for flow in the airways, where the transition from Newtonian to viscoelastic rheology in the peri-ciliary fluid may reduce clearance.

## 1. Introduction

### (a) Cilia-driven fluid mechanics

Cilia are microscopic hair-like organelles, from two to tens of microns in length, found in eukaryotic species from protozoa to mammals. Through their rapid ‘beating’ at 10–20 Hz and the resulting fluid flow generation, they perform a variety of functions, including the clearance of the protective airway surface liquid, the establishment of left/right asymmetry in the developing embryo, protozoa swimming and feeding, and the transport of gametes and embryos in the reproductive tracts. For review and research findings in cilia-driven flow, see Smith *et al.* (2008*b*) and Sleigh *et al.* (1988) for the airway surface liquid, Fauci & Dillon (2006) for the reproductive tract and Smith *et al.* (2008*a*) and Cartwright *et al.* (2007) for embryonic left/right symmetry breaking.

Cilia-driven flow involves Reynolds numbers much smaller than unity. In Newtonian liquids such as water or saline, the fluid flow is therefore accurately represented by the Stokes flow equations
1.1
where ** F** is the force per unit volume, due to drag of the cilium moving in the fluid, and μ the viscosity. Starting with the pioneering work of Hancock (1953), a number of authors have developed and applied slender-body theory for Stokes flow in order to model the propulsion of liquids by cilia and the swimming of sperm, in addition to other phenomena such as the generation of feeding currents by aquatic micro-organisms (for detailed review, see Sleigh

*et al.*1988; Smith

*et al.*2009

*a*; and for recent work, see Clarke

*et al.*2006; Smith

*et al.*2007

*a*). Slender-body theory, used to represent the cilium or flagellum, exploits the linearity of equation (1.1), in addition to the high aspect ratio or ‘slenderness’, expressing the flow velocity as line integral sums of known singular solutions.

Internal cilia, however, typically interact with liquids such as mucus which have non-Newtonian properties such as viscoelasticity. A variant on slender-body theory developed by Fulford *et al.* (1998) has been applied to the modelling of flagellar motility in a linear viscoelastic Maxwell liquid, and several recent studies have applied nonlinear viscoelastic analysis to modelling sperm motility (e.g. Lauga, 2007; Fu *et al.* 2007, 2008) and other forms of low-Reynolds-number propulsion (Normand & Lauga 2008).

These studies have provided a number of important insights regarding the generation of bending waves, effects of viscoelasticity on efficiency and the breakdown of the ‘scallop theorem’ under viscoelastic effects. All of these analyses have used the approximation of low-amplitude beating; in this study, we develop a large-amplitude theory using slender-body theory, subject to the restriction of linear viscoelasticity. We shall solve for zero-Reynolds-number flow, with rheology described by the linearised Maxwell equations:
1.2
where **τ** and **e** are deviatoric stress and rate-of-strain tensors, respectively, λ the relaxation time and μ the equivalent viscosity (e.g. Phan-Thien 2002). Generalization to other linear viscoelastic constitutive relations is straightforward within the framework we shall present.

### (b) Singularity solutions

Slender-body theory in Stokes flow involves the use of line integrals of Stokeslet singularities, corresponding to point forces and higher-order singularities such as the potential flow source dipole. A formulation used recently (Smith *et al.* 2009*a*), which represents the slender body as a curved ellipsoid, gives the fluid velocity ** u** at any point in the fluid

**as 1.3 The slender body has centreline , parametrized by scaled arclength and the vector is the force per unit length on the slender body or ‘force density’. Here and in what follows, we assume that the variables have been non-dimensionalized with scalings**

*x**L*, 1/σ, σ

*L*, μσ

*L*and μσ for length, time, velocity, force density and pressure, respectively, where

*L*is the length of the cilium and σ is the beat frequency in radians per second. The Stokes flow equations (1.1) have no explicit time dependence, so we suppress the time dependence in the variables

**ξ**,

**and**

*u***in this section. The Green function**

*f***G**is given by 1.4 where

*a*

_{0}is the slender-body radius at its midpoint. The singularities

**S**and

**K**are the Stokeslet and source dipole, respectively, which in infinite three-dimensional space have the form

*S*

_{jk}(

**;**

*x***ξ**)=(1/8π)(δ

_{jk}/

*r*+

*r*

_{j}

*r*

_{k}/

*r*

^{3}) and

*K*

_{jk}(

**;**

*x***ξ**)=(1/4π)(δ

_{jk}/

*r*

^{3}−

*r*

_{j}

*r*

_{k}/3

*r*

^{5}), where

*r*

_{j}=

*x*

_{j}−ξ

_{j}and

*r*

^{2}=|

**−**

*x***ξ**|

^{2}. The term

*S*

_{jk}(

**;**

*x***ξ**) is the

*j*th component of the velocity field at

**due to a unit point force at**

*x***ξ**in the

*k*-direction, while

*K*

_{jk}(

**;**

*x***ξ**) is the

*j*th component of the velocity field due to a unit source dipole, aligned in the

*k*-direction. The pressure field associated with a Stokeslet in the

*k*-direction is

*P*

_{k}(

**;**

*x***ξ**)=

*r*

_{j}/4π

*r*

^{3}, whereas the source dipole is associated with a constant pressure.

For the problem with prescribed cilium or flagellum movement, the surface velocity ** u**(

**ξ**) is known and hence a first-kind integral equation for the unknown force density must be solved. The kernel is nearly singular when

**is a point on the body surface close to**

*x***ξ**, which results in a well-conditioned linear system. This formulation of slender-body theory is analogous to a first-kind boundary integral method, involving line distributions of Stokeslets as opposed to surface distributions.

## 2. A mathematical formulation for the modelling of cilia-driven flow in a linear Maxwell liquid

### (a) The ‘viscoelastic Stokeslet’

A Maxwell liquid with linear time derivative is defined by equation (1.2). This can be written as , where is the invertible linear operator (1+λ∂/∂*t*)^{−1}, and indeed the theory below applies for any such operator that commutes with spatial derivatives. In index notation with the summation convention, the momentum equation in zero inertia flow with point force is
2.1
where δ(** x**,

*t*;

**ξ**,η) is a delta distribution over three space dimensions and time, so that δ(

*x*,

*t*;

**ξ**,η)

*F*

_{i}is a concentrated force per unit volume located at (

**ξ**,η). We therefore seek Green’s functions and , singular at

**=**

*x***ξ**and

*t*=η, satisfying 2.2 Using the fact that the Newtonian Green’s functions

*P*

_{j}and

*S*

_{jk}satisfy −∂

_{i}

*P*

_{j}(

**;**

*x***ξ**)+∂

_{kk}

*S*

_{ij}(

**;**

*x***ξ**)+8πδ(

**;**

*x***ξ**)δ

_{ij}=0, we can verify that 2.3 are the required solutions of equation (2.2). As commutes with spatial derivatives, the mass conservation equation is also satisfied. A similar result holds for the potential dipole and also for higher-order singularities such as the Stokes doublet.

### (b) Viscoelastic slender-body theory

For the intrinsically time-dependent viscoelastic problem, the variables ** u**,

**and**

*f***ξ**are now functions of both space and time. To formulate an analogous slender-body theory, we shall represent the solution formally as a sum in both time and space of the singularity solutions

**G**

^{Maxwell}. We assume as a starting condition trivial dynamics for

*t*<0 and Stokes dynamics for

*t*=0. For

*t*>0, we have 2.4 The notation denotes . Note that the Maxwell kernel is now replaced by the Newtonian kernel, with the additional complication of the delta distribution and temporal derivative. Using equation (A1) with , we have 2.5 Hence the delta distribution and temporal integral are removed, and the problem is reduced to expressions from the Newtonian slender-body theory, with the addition of a leading time derivative.

### (c) Image systems representing a no-slip boundary

Cilia extend from a cell surface, which acts as a no-slip, no-penetration boundary for fluid flow. It is convenient to model this mathematically as a planar surface *x*_{3}=0, on which we have the condition ** u**(

*x*

_{3}=0)=0. To take account of this condition in Newtonian slender-body theory, Blake (1971) derived the following stokeslet and image system: 2.6 The image system is located at (ξ

_{1},ξ

_{2},−ξ

_{3}), so that

*R*

_{1}=

*r*

_{1},

*R*

_{2}=

*r*

_{2},

*R*

_{3}=

*x*

_{3}+ξ

_{3}and , and in the above, the Arabic indices take values 1, 2 and 3 and the index α takes values 1 and 2 only. The image system may be added directly in the viscoelastic slender-body theory because the equations are linear, augmenting the kernel

**G**and the condition

*u*_{3}(

*x*

_{3}=0)=0 will therefore be satisfied. Implementations of slender-body theory (e.g. Smith

*et al.*2007

*a*, 2009

*a*) typically do not add image corrections for the higher-order potential source dipole because of its rapid

*O*(1/

*r*

^{3}) decay.

### (d) Numerical implementation

The initial force distribution at *t*=0 can be found using the methods of Smith *et al.* (2007*a*), as the dynamics are equivalent to that of the Newtonian problem. In what follows, we discuss the calculation at subsequent timesteps *t*^{n}>0. To discretize the time derivative, we use the backward Euler finite difference, with *t*^{n}:=(*n*−1)δ*t*. For compactness, we use the notation *u*^{n}(** x**):=

**(**

*u***,**

*x**t*

^{n}), and similar notation for and . The time-discretized equations are 2.7

In order to discretize spatially, Smith *et al.* (2007*a*, 2009*a*) used piecewise constant basis functions. For this problem, greater accuracy is required, so we choose piecewise cubic Lagrange basis functions which we denote ϕ_{q}, so that is approximated by with summation over *q*=1,…,*N* implied. The explicit form of the basis functions is given in appendix B. There are *N* degrees of freedom for the coefficients , equivalent to the values of the discretized force at the nodes *s*_{q}=(*q*−1/2)/*N*. Hence we now have
2.8

To solve for *t*=*t*^{n}, we impose a prescribed cilium velocity boundary condition at the cilium surface ** x**=

*X*^{n}(

*s*) as opposed to the centreline

**ξ**

^{n}(

*s*). As shown by Smith

*et al.*(2007

*a*, 2009

*a*), the kernel given in equation (1.4) reduces, but does not remove completely, the error in the fluid velocity around the cross section of a curved slender body in a non-uniform flow. Because we wish to minimize any error propagating in time, we use the following approach, which averages the velocity between points on opposite sides of the cross section. We take 2

*N*collocation points

*X*^{±}(

*s*

_{i})=

**ξ**

^{n}(

*s*

_{i})±

*a*(

*s*

_{i})

**(**

*b**s*

_{i}), where

**is the unit binormal—typically (0,1,0) for motion in the**

*b**x*

_{1}

*x*

_{3}-plane— and

*a*(

*s*

_{i}) is the cilium radius. Taking the sum of the equations formed by collocating at

*X*^{+}(

*s*

_{i}) and

*X*^{−}(

*s*

_{i}), we have the following equation for

*i*=1,…,

*N*: 2.9 The boundary condition velocity satisfies

*u*^{n}(

*X*^{+}(

*s*

_{i}))=

*u*^{n}(

*X*^{−}(

*s*

_{i})), which we denote

*u*^{n}(

*s*

_{i}). Rearranging, we have 2.10 For the purpose of time-stepping, the left-hand side of the equation contains known values and the right-hand side contains the unknown force distribution . The second term on the left-hand side can be characterized as a ‘memory velocity’, which we denote . In the numerical implementation, we represent using a cubic smoothing-spline representation (NAG library E02BEF), with up to 16 points that accurately represents this term. The purpose of this is to ensure that small perturbations do not grow unstably and this ensures that results may be calculated for a wider range of values of the relaxation time λ.

Defining the following coefficients
2.11
which may be evaluated by an adaptive quadrature routine, e.g. the NAG library D01AJF, we have the linear algebra problem for the unknown force strengths at the points *s*_{q}. Each **A**_{iq} is a 3×3 matrix and so the matrix problem has 3*N*×3*N* scalar coefficients. The system is well conditioned due to the fact that when *q*=*i*, the surface points *X*^{±}(*s*_{q}) are close to **ξ**(*s*_{i}) and hence the values of the Green’s function *G*_{jj}(*X*^{±}(*s*_{q}),**ξ**(*s*_{i})) are large but finite. This gives a diagonally dominant problem, as first noted by Liron & Mochon (1976).

## 3. Results

All results are calculated using the ‘9+2’ planar cilium beat cycle given by Fulford & Blake (1986), as described in appendix C and shown in figure 1. This beat pattern, involving asymmetric ‘effective’ and ‘propulsive’ strokes, was obtained from analysis of tracheal electron micrographs. We also chose the parameter *a*_{0}=0.01, corresponding to the cilium having radius 0.1 μm and length 10 μm. The cilium was discretized with *N*=12 collocation points. Convergence and accuracy testing for the numerical algorithm are discussed in appendix D.

### (a) Force and bending moment density

Figures 2 and 3 show the computed solution ** f**(

*s*,

*t*) for four beat cycles in approximately Newtonian and viscoelastic cases, respectively. The four effective strokes are visualized as red streaks in panel

*a*in each figure, whereas the zero lateral force computed with our algorithm is evident in panel

*b*. The force density is greatest at or near the cilium tip, due to it having maximum velocity, with nearly zero force close to the base. Comparing the figures, it is evident that the essential ‘shape’ of the force distribution is similar in the viscoelastic and Newtonian cases, although the magnitude is reduced. In order to investigate this further, we plot in figure 4 the calculated normal component of the force density

**(**

*n**s*,

*t*)⋅

**(**

*f**s*,

*t*) at the midpoint of the cilium, at

*t*=0.733 early in the effective stroke. A non-monotonic variation was found, as for small values of λ (0–0.02) the force density slightly increases, but thereafter decreases almost linearly as λ increases to 0.2. The instantaneous midpoint normal force density at

*t*=0.733 with λ=0.2 is around 76 per cent of the Newtonian value. This ratio varies depending on time

*t*, from around 50 per cent to more than 130 per cent, with average approximately 73 per cent.

The bending moment density resulting from the force distribution, calculated as in appendix E, is shown in figure 5, comparing the Newtonian fluid (*a*) and viscoelastic fluid (*b*). Again, the shape of the results in both fluids is similar, the magnitude being reduced in the viscoelastic case. The bending moment density differs from force density in that the largest values occur near the base of the cilium, occurring because of the combined effect of all of the force exerted by the cilium on the fluid along its length.

### (b) Power dissipation

Power dissipation was calculated from the force density using equation (E3). Figure 6 shows how increasing the relaxation time from nearly zero to a dimensionless value of 0.2, equivalent to 0.032 of the cilium beat period, results in an approximately exponential decay in the power consumption of the cilium beat, occurring due to the elastic storage and release of energy by the viscoelastic liquid. With relaxation time λ=0.2, approximately one-tenth of the effective stroke duration, the power consumption is halved.

### (c) Velocity field

Example results for the three-dimensional variation in the fluid velocity magnitude |** u**(

**,**

*x**t*)|, at four points during the beat cycle, are shown in figure 7, for relaxation time λ=0.2. The velocity field decays rapidly in the

*x*

_{2}-direction. Figure 8 shows an instantaneous flow velocity profile in the plane

*x*

_{2}=0 for values of relaxation time 10

^{−6}and 0.2. The flow velocity is smaller in the viscoelastic case, although the decay rate in the ‘far-field’, evaluated at

*x*

_{3}=3 is similar. In the example shown, the horizontal velocity decays with approximately in the Newtonian case and in the viscoelastic case, both values approaching the far-field decay of the Stokeslet and image system of Blake (1971).

### (d) Particle transport

Using a simple Eulerian particle tracking scheme, , we simulated the transport of a test particle initially at (0,0,1.05), just above the cilium tip. Example results are shown in figure 9, with initial position labelled ‘a’ and final position after four cycles of beating labelled ‘b’. The Newtonian and viscoelastic cases differ significantly, with the introduction of elastic properties resulting in reduced particle transport during the effective stroke, and additionally proportionately greater recoil during the recovery stroke, greatly reducing the overall particle clearance. The variation of transport for this test particle in the *x*_{1}-direction with relaxation time is shown in figure 10. Transport for λ=0.2 is just 21 per cent of transport for λ=10^{−6}. The results show an excellent agreement (*R*^{2}=0.99) with an exponential best-fit curve. Simulations were also carried out for initial positions (0,0,1.20), and again the results were fitted very closely by an exponential best-fit curve, although with differing coefficients.

### (e) Volume flow rate

The volume flow rate *Q*(*t*) is a measure of the instantaneous pumping effect of the cilium movement. Calculation of *Q*(*t*) exploiting Newtonian Stokes flow theory is discussed in appendix F. Mean volume flow rate versus relaxation time is plotted in figure 11. Newtonian flow results in the highest value of , whereas for λ=0.2 the value of is negative, implying a mean backward flow of fluid in the domain . Despite the small negative flux averaged over , in certain regions such as near the cilium tip, particles may be transported forward as shown in the example result in figure 9.

## 4. Discussion

### (a) The viscoelastic Stokeslet and slender-body theory

The viscoelastic Stokeslet derived in this article may be useful in mathematical modelling in a number of microscopic flow problems, not restricted to cilia and flagella movement. The framework presented with the Maxwellian or other viscoelastic singularity solutions may also form the basis for a novel boundary integral formulation of zero-Reynolds-number linear viscoelastic problems. In particular, the use of viscoelastic singularities may circumvent the inhomogeneous equations, and associated substantial calculations, which result from the viscoelastic contribution to the deviatoric stress tensor in boundary integral approaches that can accommodate creeping linear viscoelastic flows (e.g. Phan-Thien and Fan 2002). The ‘Maxwell Stokeslet’ or higher-order singularities such as the Maxwell Stokes dipole may also prove useful in modelling bacterial populations as self-propelling particles in viscoelastic media. Viscoelasticity may have important effects on bacterial fluid dynamics and may, for example, be an important factor in the establishment of bacterial colonies (Matsui *et al.* 2006). The development of an equivalent viscoelastic ‘Stokesian dynamics’ may be possible for the modelling of the movement of large numbers of particles in viscoelastic media.

### (b) Rheological properties of mucus and other polymers

Fitting a Maxwell model to oscillatory rheometry data leads to estimates of the relaxation time of 0.027 s for respiratory mucus (Smith *et al.* 2008*b*). However, the airway surface liquid in the healthy state is differentiated into an upper layer of mucus and a lower layer of ‘peri-ciliary liquid’, the latter being the liquid in which the cilia move. Peri-ciliary liquid is generally considered to be more ‘watery’, typically having a much lower polymer concentration and hence shorter elastic relaxation time, at least in the healthy state (for further references, see Smith *et al.* 2008*b*). There does not appear to be any published data on the liquid in the embryonic node; however, it is known that signalling proteins are present (for review, see Raya & Izpisúa Belmonte 2006) and these may cause increased viscoelasticity. Methylcellulose solutions are often used in biological imaging experiments on ciliated epithelium (Gheber & Priel 1998) and sperm migration (Ivic *et al.* 2002; Smith *et al.* 2009b). From rheometry experiments on 1 per cent solutions we estimate the relaxation time to be approximately 0.006 s (Smith *et al.* 2006*b*).

### (c) Transport of particles and fluid pumping in a viscoelastic medium

In contrast to the results of Fulford *et al.* (1998), in which it was found that cell propulsion is not significantly affected by relaxation time, we found that particle transport by cilia is significantly reduced by viscoelasticity. Even over the relatively modest range of relaxation time λ explored in our study, transport varied fivefold. For fixed initial position and varying relaxation time, the predicted transport was approximated remarkably well by an exponential best-fit curve. The underlying physics do not appear to be obvious, although we note that integral formulations of viscoelastic constitutive equations involve an exponentially decaying factor representing the elastic ‘memory effect’ of the fluid. It should also be noted that force density, bending moment density and power of cilia beating did not show such clear exponential variation with relaxation time, the force density, in particular, showing non-monotone behaviour for small values of λ.

The mean volume flow rate at λ=0.2 was negative, implying that despite the cilium sweeping closer to the epithelium during the recovery stroke than the effective stroke, the recovery stroke transported fluid in the negative direction more effectively. Zero-Reynolds-number viscous flow is well known to be reversible, in the sense that symmetric effective and recovery strokes will produce zero net propulsion. However, asymmetric effective and recovery strokes produce a positive net flow. This contrasts with a solid-like material exhibiting a purely elastic response, which will not exhibit flow. The viscous or elastic response of a viscoelastic material depends on the timescale of the driving force relative to the relaxation time of the material. As the recovery stroke occurs over a longer timescale than the effective stroke, for sufficiently long relaxation time, viscous flow during the effective stroke becomes relatively smaller in magnitude than the viscous flow during the recovery stroke. This may explain why negative mean flow occurs for sufficiently large λ.

### (d) Biological implications

Our results suggest that a relatively short relaxation time of a fraction of the ciliary beat, less than 0.032, can significantly alter the mechanics of cilia movement and resulting particle transport. At a beat frequency of 10 Hz, this corresponds to a dimensional relaxation time of 0.0032 s. Hence in the airways, even slightly increased viscoelasticity of the peri-ciliary layer, in which cilia beat, may be harmful to particle transport. It is thought that the rheological properties of this layer may be a factor in the effective functioning of the muco-ciliary system (Knowles & Boucher 2002). Increased peri-ciliary liquid viscoelasticity may occur because of the presence of DNA and albumin, as found in the sputum of asthmatic patients, and/or because of mucus ‘tethering’ to the epithelium, as observed in the asthmatic lung (see Rogers 2004 for further details). Rheology may alter transport in the embryonic node, in which cilia propel morphogen-containing vesicles by performing a circular ‘whirling’ motion rather than the nearly planar beat studied here. Effective transport of morphogen-containing vesicles in this cavity is likely to be the mechanism by which left/right asymmetry is established in the developing embryo, and the rheological properties of the extraembryonic liquid may affect both the motion and the particle transport by these cilia.

In the airway muco-ciliary system, a high-viscosity mucous layer overlies the peri-ciliary liquid and this leads to complex fluid dynamic interactions (Smith *et al.* 2008*b*). Two important effects are (i) cilia generally engage with this liquid only during the effective stroke, providing an increased propulsive force and (ii) because the mucous layer remains relatively flat, pressure gradients are predicted to ‘assist’ with liquid propulsion during the effective stroke (Smith *et al.* 2007*b*). These effects may compensate for the abolition of positive flow caused by viscoelasticity predicted in our study.

### (e) Future work

It may be instructive to develop a formulation of this method for the modelling of flagellar motion and the resulting motility of sperm, bacteria and other cells in viscoelastic liquids with relatively short relaxation time. Such motion differs from that of cilia in that movement of the slender body tangential to its axis is important, in addition to motion that is primarily normal. The modelling of the sperm head or cell body will also prove an important challenge. When data are available on the contents and rheological properties of the very small volumes of liquid in the embryonic node, it will be possible to determine how such properties may affect morphogen transport.

The linearized viscoelastic theory can only be rigorously justified for small-amplitude movements, whereas in this study we considered large-amplitude movements of a cilium. Neglecting nonlinear effects is necessary to formulate slender-body theory; however, future comparison with computational solutions of the fully nonlinear equations will be important in evaluating the accuracy of the theory. Our model does not take into account the cilium internal mechanics, but rather specifies the beat pattern and frequency. However, in the biological system, these are not prescribed, but rather emerge from the interaction of an internal bending moment and external forces from the fluid. How the beat pattern and frequency are affected by viscoelasticity is an important subject for biophysical modelling, in the same way that the effect of viscous load has previously been addressed (Gueron & Liron 1998; Gueron & Levit-Gurevich 2001; Brokaw 2002; Dillon *et al.* 2006). In the embryonic node, cilia are relatively widely distributed (Nonaka *et al.* 1998); however, in the reproductive tracts and the airway, cilia density may be very high, approximately 6–10 μm^{−2} (Sleigh *et al.* 1988), and so modelling the effect of multiple cilia beating in a viscoelastic fluid is an important subject for future work. Finally, the modelling of bacteria in gels such as distributions of singularity solutions and even the development of a viscoelastic ‘Stokesian dynamics’ may be possible.

## Acknowledgements

The authors thank Dr Jackson Kirkman-Brown of the Centre for Human Reproductive Science, Birmingham Womens Hospital, Mr Henry Shum and Mr Hermes Gadêlha, Mathematical Institute, University of Oxford, and members of the Reproductive Biology and Genetics Group, University of Birmingham, for invaluable discussions. D.J.S. was supported by MRC Special Training Fellowship G0600178 in Computational Biology; part of this research was performed while D.J.S. was supported by a Wellcome Trust ‘Value in People’ Fellowship.

## Appendix A. Required lemma to formulate the slender-body theory

To derive the slender-body theory (2.5), we use the following lemma: for any differentiable function *H*,
A1
The result is derived as follows. Making the change of variable ζ=*t*−η, we have
A2
For any test function ϕ, the identities and hold (e.g. Lighthill 1964). Taking *H*(*t*−ζ) as the test function and noting that (d/dζ)*H*(*t*−ζ)=−*H*′(*t*−ζ), it follows that the right-hand side of equation (A2) is equal to
A3
for all *T*_{1} and *T*_{2} such that *T*_{1}<0<*t*<*T*_{2}. Hence this gives the limiting value as *T*_{1}↑0 and *T*_{2}↓*t*, proving equation (A1).

## Appendix B. Cubic basis functions

The basis functions ϕ_{q}(*s*) are chosen so that the discretized force *f*_{q}ϕ_{q}(*s*) (summation over *q*=1,…,*N* implied) is constant on intervals *I*_{1} and *I*_{N+1} and is a Lagrange cubic interpolating function elsewhere. To give the explicit form for the basis functions, we use the following notation: each interval *I*_{q}:=(*s*_{q−1},*s*_{q}), and the nodal value of the force ** f**(

*s*

_{q})=

*f*_{q}. The nodes are given by

*s*

_{q}=(

*q*−1/2)/

*N*for

*q*=1,…

*N*, while at the ends of the cilium,

*s*

_{0}=0 and

*s*

_{N+1}=1.

The explicit expressions for the first two basis functions are
B1
The expressions for ϕ_{N−1} and ϕ_{N} are analogous to those for ϕ_{2} and ϕ_{1} respectively. For 2<*q*<*N*−1, the basis functions are given by
B2

## Appendix C. Cilia beat cycle specification

We use the mathematical formulation of the tracheal cilium beat cycle given by Fulford & Blake (1986), based on the micrograph observations of Sanderson & Sleigh (1981). Visualizations of the cilium beat are given in figures 1 and 9*a*,*c*. The beat cycle consists of two phases, an effective/propulsive stroke, during which the cilium is relatively straight, and a recovery/preparatory stroke, during which the cilium returns to its initial position but bends closer to the cell surface. The recovery stroke is approximately twice the duration of the effective stroke. The mathematical specification is
C1
where *s* is scaled with respect to cilium length, *t* is scaled with respect to inverse radian frequency, so that *t*=2π corresponds to one complete beat cycle, and the functions **α**_{p}(*s*) and **β**_{p}(*s*) are cubic polynomials, the coefficients being given in Fulford & Blake (1986).

## Appendix D. Convergence of the algorithm and boundary condition testing

Recall that time is non-dimensionalized with respect to inverse radian frequency, so that 0<*t*<2π corresponds to one complete beat cycle. Results were calculated with *N*=12 degrees of freedom for the force density *f*^{n}(*s*) at each timestep. This was verified to be accurate by doubling the number of degrees of freedom and noting that the root mean square error was within 1.3 per cent. We found that temporally converged results for the force density ** f**(

*s*,

*t*) were obtained with timestep 120th of the cilium beat cycle, δ

*t*=2π/120. Results were calculated for a range of non-dimensional relaxation times from λ=10

^{−6}to 0.2. For a beat frequency of 10 Hz, the parameter λ=0.2 corresponds to a dimensional relaxation time of 0.2/(2π×10)=0.0032 s.

As described in Smith *et al.* (2007*a*), the accuracy of the numerical algorithm can be tested post hoc by comparing the boundary condition velocity ** U**(

*s*,

*t*) with the fluid velocity

**(**

*u***ξ**(

*s*,

*t*)), as shown in figure 12. Results are satisfactory along the great majority of the length of the cilium, with only small deviations from the expected values near the base and tip of the cilium, largely occurring in the

*u*

_{3}component. There are relatively small oscillations between collocation points, which are introduced by the ‘averaged’ collocation scheme of equation (2.9) that sums the velocity fields at

*X*^{+}(

*s*

_{i}) and

*X*^{−}(

*s*

_{i}). The alternative is to use collocation points on only one side of the cilium, corresponding to just the points

*X*^{+}(

*s*

_{i}), which has proved accurate in Newtonian slender-body theory (Smith

*et al.*2007

*a*, 2009

*a*). However, in the viscoelastic formulation, unphysical ‘lateral’ forces

*f*

_{2}(

*s*,

*t*) are introduced by this, so we employ the averaged approach.

## Appendix E. Bending moment density and power dissipation

Bending of the cilium will be produced by a balance of viscous and elastic moments due to the fluid and passive structures of the cilium, and active moments due to dynein molecular motor activity (e.g. Hines and Blum 1978; Brokaw 2002). It is instructive therefore to consider the bending moment density along the cilium, and now it varies with the rheological properties of the liquid. Following Hines & Blum (1978), we may derive an expression for the bending moment density as follows: suppose that the cilium is instantaneously cut at arclength *s* and time *t*. In order to maintain the motion of the distal end of the cilium, a force ** F**(

*s*,

*t*) and a moment

*M*(

*s*,

*t*) must be applied at the cut end, replacing the effect of the proximal end of the cilium. Using the fact that the viscous force per unit length on the cilium is for and the distal end of the cilium at is not subject to any forces or moments other than the cilium and surrounding fluid, we may form the following moment balance equation at

*s*: E1 As ∂

**ξ**/∂

*s*is tangent to

**ξ**, the vector product

**∧∂**

*b***ξ**/∂

*s*is equal to the normal

**(**

*n**s*,

*t*). Differentiating equation (E1) with respect to

*s*leads to the following expression for the bending moment density

*m*(

*s*,

*t*):=∂

*M*/∂

*s*: E2 During the effective stroke,

**(**

*f**s*,

*t*) will in general be positive, indicating a force on the fluid in the positive

*x*

_{1}-direction, whereas

*m*(

*s*,

*t*) will in general be negative, indicating a ‘clockwise’ bending moment.

Instantaneous fluid-dynamic power dissipation *P*(*t*) is calculated by integrating force density multiplied by cilium centreline velocity as follows:
E3
The time-averaged power consumption is estimated from four beat cycles as (1/(8π)) , where *t* is scaled with respect to inverse radian frequency.

## Appendix F. Volume flow rate

The pumping effect of a cilium in the *x*_{1}-direction can be quantified using the volume flow rate *Q*(*t*),
F1
the choice of *x*_{1} being arbitrary due to conservation of mass. In this study, the mean volume flow rate over four beat cycles is denoted .

The volume flow rate due to a Stokeslet singularity *S*_{1j} near to a no-slip surface was given by Liron (1978) as
F2
where *x*_{1}, ξ_{1} and ξ_{2} are arbitrary. In non-dimensional variables, the factor μ in the denominator is absent.

The volume flow rate due to the infinite domain potential dipole *K*_{1j} is zero for *j*=1,2,3. For *j*=2,3, this can be determined by inspection: choosing ξ_{1}=ξ_{2}=ξ_{3}=0, we have *K*_{1j}=−3*x*_{1}*x*_{j}/*r*^{5}. For any *x*_{1}≠0, the integrand is a regular odd function of *x*_{j}, hence the integral with respect to *x*_{j} is zero. For *j*=1, we have
F3
The indefinite integral in the right-hand side with respect to *x*_{3} is and hence the definite integral over the infinite domain is zero. In our model, we have restricted attention to the half-space *x*_{3}>0; however, it can be seen from the above and equation (1.4) that the error in the volume flow rate resulting from this will be compared with the flow due to the Stokeslet, and so we neglect it.

Using the time-discretized equation (2.7) and the expression (F2) for the Stokeslet, we have the following result for the volume flow rate at the *n*th timestep *Q*(*t*^{n}) in terms of the cilium force density:
F4

## Footnotes

- Received January 12, 2009.
- Accepted April 28, 2009.

- © 2009 The Royal Society