Mathematical modelling of cilia-driven transport of biological fluids

D. J. Smith, E. A. Gaffney, J. R. Blake


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. (2008b) and Sleigh et al. (1988) for the airway surface liquid, Fauci & Dillon (2006) for the reproductive tract and Smith et al. (2008a) 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 Embedded Image 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. 2009a; and for recent work, see Clarke et al. 2006; Smith et al. 2007a). 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: Embedded Image 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. 2009a), which represents the slender body as a curved ellipsoid, gives the fluid velocity u at any point in the fluid x as Embedded Image 1.3 The slender body has centreline Embedded Image, parametrized by scaled arclength Embedded Image and the vector Embedded Image 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 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 ξ, u and f in this section. The Green function G is given by Embedded Image 1.4 where a0 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 Sjk(x;ξ)=(1/8π)(δjk/r+rjrk/r3) and Kjk(x;ξ)=(1/4π)(δjk/r3rjrk/3r5), where rj=xj−ξj and r2=|xξ|2. The term Sjk(x;ξ) is the jth component of the velocity field at x due to a unit point force at ξ in the k-direction, while Kjk(x;ξ) is the jth 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 Pk(x;ξ)=rj/4πr3, 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 Embedded Image must be solved. The kernel is nearly singular when x is a point on the body surface close to ξ, 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 Embedded Image, where Embedded Image 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 Embedded Image 2.1 where δ(x,t;ξ,η) is a delta distribution over three space dimensions and time, so that δ(x,t;ξ,η)Fi is a concentrated force per unit volume located at (ξ,η). We therefore seek Green’s functions Embedded Image and Embedded Image, singular at x=ξ and t=η, satisfying Embedded Image 2.2 Using the fact that the Newtonian Green’s functions Pj and Sjk satisfy −∂iPj(x;ξ)+∂kkSij(x;ξ)+8πδ(x;ξij=0, we can verify that Embedded Image 2.3 are the required solutions of equation (2.2). As Embedded Image 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, f and ξ 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 GMaxwell. We assume as a starting condition trivial dynamics for t<0 and Stokes dynamics for t=0. For t>0, we have Embedded Image 2.4 The notation Embedded Image denotes Embedded Image. 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 (A 1) with Embedded Image, we have Embedded Image 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 x3=0, on which we have the condition u(x3=0)=0. To take account of this condition in Newtonian slender-body theory, Blake (1971) derived the following stokeslet and image system: Embedded Image 2.6 The image system is located at (ξ12,−ξ3), so that R1=r1, R2=r2, R3=x33 and Embedded Image, 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 u3(x3=0)=0 will therefore be satisfied. Implementations of slender-body theory (e.g. Smith et al. 2007a, 2009a) typically do not add image corrections for the higher-order potential source dipole because of its rapid O(1/r3) decay.

(d) Numerical implementation

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

In order to discretize Embedded Image spatially, Smith et al. (2007a, 2009a) 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 Embedded Image is approximated by Embedded Image 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 Embedded Image, equivalent to the values of the discretized force at the nodes sq=(q−1/2)/N. Hence we now have Embedded Image 2.8

To solve for t=tn, we impose a prescribed cilium velocity boundary condition Embedded Image at the cilium surface x=Xn(s) as opposed to the centreline ξn(s). As shown by Smith et al. (2007a, 2009a), 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 2N collocation points X±(si)=ξn(sia(si)b(si), where b is the unit binormal—typically (0,1,0) for motion in the x1x3-plane— and a(si) is the cilium radius. Taking the sum of the equations formed by collocating at X+(si) and X(si), we have the following equation for i=1,…,N: Embedded Image 2.9 The boundary condition velocity satisfies un(X+(si))=un(X(si)), which we denote un(si). Rearranging, we have Embedded Image 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 Embedded Image. The second term on the left-hand side can be characterized as a ‘memory velocity’, which we denote Embedded Image. In the numerical implementation, we represent Embedded Image 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 Embedded Image 2.11 which may be evaluated by an adaptive quadrature routine, e.g. the NAG library D01AJF, we have the linear algebra problem Embedded Image for the unknown force strengths Embedded Image at the points sq. Each Aiq is a 3×3 matrix and so the matrix problem has 3N×3N scalar coefficients. The system is well conditioned due to the fact that when q=i, the surface points X±(sq) are close to ξ(si) and hence the values of the Green’s function Gjj(X±(sq),ξ(si)) 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 a0=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.

Figure 1.

The effective/propulsive stroke (a) and recovery/preparatory stroke (b) of the cilium beat from the experimental data of Sanderson & Sleigh (1981) and the mathematical description of Fulford & Blake (1986). The beat is divided into 120 timesteps, the effective stroke spanning steps 1 to 39 and the recovery stroke steps 40 to 120.

(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.

Figure 2.

Newtonian fluid. Force components per unit length, (a) f1(s,t), (b) f2(s,t), (c) f3(s,t), scaled with respect to μσL, where σ is the frequency in cycles per unit time and L the cilium length, for approximately Newtonian liquid. Non-dimensional relaxation time λ=10−6, equivalent to 1.6×10−8s for frequency σ=10 Hz.

Figure 3.

Viscoelastic fluid. Force components per unit length, (a) f1(s,t), (b) f2(s,t), (c) f3(s,t), scaled with respect to μσL, where σ is the frequency in cycles per unit time and L the cilium length, for the most elastic liquid simulated. Non-dimensional relaxation time λ=0.2, equivalent to 0.0032 s for frequency σ=10 Hz.

Figure 4.

Normal component of force density, n(s,t)⋅f(s,t) scaled with respect to μσL, sampled at the midpoint of the cilium (s=0.5) at time t=0.733, scaled with respect to inverse frequency, versus non-dimensional relaxation time.

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.

Figure 5.

Bending moment density m(s,t), as defined in appendix E over four cycles of the cilium beat, scaled with respect to μσL2, where σ is the radian frequency and L the cilium length. (a) Approximately Newtonian liquid, with non-dimensional relaxation time λ=10−6, equivalent to 1.6×10−8 s for frequency σ=10 Hz. (b) Non-dimensional relaxation time λ=0.2, equivalent to 0.0032 s for σ=10 Hz.

(b) Power dissipation

Power dissipation was calculated from the force density using equation (E 3). 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.

Figure 6.

Time-averaged power consumption Embedded Image, scaled with respect to μσ2L3.

(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 x2-direction. Figure 8 shows an instantaneous flow velocity profile in the plane x2=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 x3=3 is similar. In the example shown, the horizontal velocity decays with approximately Embedded Image in the Newtonian case and Embedded Image in the viscoelastic case, both values approaching the Embedded Image far-field decay of the Stokeslet and image system of Blake (1971).

Figure 7.

Example results, showing the magnitude of the velocity field |u(x,t)|, at x3=0.2,0.4,0.6,0.8 at four points during the beat cycle, with length scaled with respect to cilium length L and time scaled with respect to inverse radian frequency (ad). The velocity field is calculated from the solutions Embedded Image using equation (2.8).

Figure 8.

Instantaneous u1 velocity profile from simulations of Newtonian (λ=10−6) and viscoelastic (λ=0.2) liquids, calculated at x1=0=x2 and 0.15<x3<3.0 in units scaled with respect to cilium length L. The profile is shown in non-dimensional units, scaled with respect to σL and multiplied by a factor of 5 for visual clarity (dashed line: u1, scaled velocity profile (Newtonian); dotted line: u1, scaled velocity profile (λ=0.2); straight line: cilium centreline).

(d) Particle transport

Using a simple Eulerian particle tracking scheme, Embedded Image, 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 x1-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 (R2=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.

Figure 9.

Transport of a test particle in the x1-direction, in Newtonian (a,b:λ=10−6) and viscoelastic (b,d:λ=0.2) liquids over four cycles of the cilium beat. The initial location is marked ‘a’, the final location ‘b’. (b,d) Dotted lines show particle trajectory during the effective stroke, and solid lines during the recovery stroke.

Figure 10.

Transport of a test particle in the x1-direction, initially at position (0,0,1.05L) over one cycle of cilium beating, versus non-dimensional relaxation time, with empirical best-fit curve (open circle: simulation results; straight line: empirical fit Embedded Image).

(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 Embedded Image versus relaxation time is plotted in figure 11. Newtonian flow results in the highest value of Embedded Image, whereas for λ=0.2 the value of Embedded Image is negative, implying a mean backward flow of fluid in the domain Embedded Image. Despite the small negative flux averaged over Embedded Image, in certain regions such as near the cilium tip, particles may be transported forward as shown in the example result in figure 9.

Figure 11.

Mean dimensionless volume flow rate versus relaxation time λ. For λ=0.2, equivalent to 0.0032 s for σ=10 Hz, the mean volume flow is negative.

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. 2008b). 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. 2008b). 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. 2006b).

(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. 2008b). 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. 2007b). 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.


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, Embedded Image A 1 The result is derived as follows. Making the change of variable ζ=t−η, we have Embedded Image A 2 For any test function ϕ, the identities Embedded Image and Embedded Image 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 (A 2) is equal to Embedded Image A 3 for all T1 and T2 such that T1<0<t<T2. Hence this gives the limiting value as T1↑0 and T2t, proving equation (A 1).

Appendix B. Cubic basis functions

The basis functions ϕq(s) are chosen so that the discretized force fqϕq(s) (summation over q=1,…,N implied) is constant on intervals I1 and IN+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 Iq:=(sq−1,sq), and the nodal value of the force f(sq)=fq. The nodes are given by sq=(q−1/2)/N for q=1,…N, while at the ends of the cilium, s0=0 and sN+1=1.

The explicit expressions for the first two basis functions are Embedded Image B 1 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 Embedded Image B 2

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 9a,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 Embedded Image C 1 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 fn(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. (2007a), 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 u3 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+(si) and X(si). The alternative is to use collocation points on only one side of the cilium, corresponding to just the points X+(si), which has proved accurate in Newtonian slender-body theory (Smith et al. 2007a, 2009a). However, in the viscoelastic formulation, unphysical ‘lateral’ forces f2(s,t) are introduced by this, so we employ the averaged approach.

Figure 12.

Post hoc test of the numerical solution, comparing the fluid velocity field on the surface of the cilium X+(s,t) for 0<s<1 with the velocity boundary condition at the collocation points u(si,t) for i=1,…,N (dots: fluid velocity field; open circle: velocity boundary condition at collocation point).

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 Embedded Image for Embedded Image and the distal end of the cilium at Embedded Image 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: Embedded Image E 1 As ∂ξ/∂s is tangent to ξ, the vector product b∧∂ξ/∂s is equal to the normal n(s,t). Differentiating equation (E 1) with respect to s leads to the following expression for the bending moment density m(s,t):=∂M/∂s: Embedded Image E 2 During the effective stroke, f(s,t) will in general be positive, indicating a force on the fluid in the positive x1-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: Embedded Image E 3 The time-averaged power consumption Embedded Image is estimated from four beat cycles as (1/(8π)) Embedded Image, where t is scaled with respect to inverse radian frequency.

Appendix F. Volume flow rate

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

The volume flow rate due to a Stokeslet singularity S1j near to a no-slip surface was given by Liron (1978) as Embedded Image F 2 where x1, ξ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 K1j is zero for j=1,2,3. For j=2,3, this can be determined by inspection: choosing ξ123=0, we have K1j=−3x1xj/r5. For any x1≠0, the integrand is a regular odd function of xj, hence the integral with respect to xj is zero. For j=1, we have Embedded Image F 3 The indefinite integral in the right-hand side with respect to x3 is Embedded Image and hence the definite integral over the infinite domain Embedded Image is zero. In our model, we have restricted attention to the half-space x3>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 Embedded Image compared with the flow due to the Stokeslet, and so we neglect it.

Using the time-discretized equation (2.7) and the expression (F 2) for the Stokeslet, we have the following result for the volume flow rate at the nth timestep Q(tn) in terms of the cilium force density: Embedded Image F 4


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


View Abstract