## Abstract

A boundary element implementation of the regularized Stokeslet method of Cortez is applied to cilia and flagella-driven flows in biology. Previously published approaches implicitly combine the force discretization and the numerical quadrature used to evaluate boundary integrals. By contrast, a boundary element method can be implemented by discretizing the force using basis functions, and calculating integrals using accurate numerical or analytic integration. This substantially weakens the coupling of the mesh size for the force and the regularization parameter, and greatly reduces the number of degrees of freedom required. When modelling a cilium or flagellum as a one-dimensional filament, the regularization parameter can be considered a proxy for the body radius, as opposed to being a parameter used to minimize numerical errors. Modelling a patch of cilia, it is found that: (i) for a fixed number of cilia, reducing cilia spacing reduces transport, (ii) for fixed patch dimension, increasing cilia number increases the transport, up to a plateau at 9×9 cilia. Modelling a choanoflagellate cell, it is found that the presence of a lorica structure significantly affects transport and flow outside the lorica, but does not significantly alter the force experienced by the flagellum.

## 1. Introduction

Fluid dynamic phenomena on microscopic scales are fundamental to life, from the feeding and swimming of bacteria and single-celled eukaryotes to complex roles in the organ systems of mammals. Examples include sperm swimming, ovum and embryo transport, respiratory defence, auditory perception and embryonic symmetry-breaking. Cilia, which range from a few to tens of microns in length, are hair-like organelles that project from cell surfaces either in dense mats or one per cell, and through specialized side-to-side or rotational beating patterns and viscous interaction with the cell surface, propel fluid in a parallel direction. Flagella may be tens to hundreds of microns in length, typically occurring singly or two per cell, and generally produce propulsion tangential to the flagellum axis through the propagation of a bending wave, resulting in swimming as in sperm, or the generation of feeding currents as in *choanoflagellates*. These phenomena have since the work of Taylor, and Gray and Hancock in the 1950s, motivated the study of methods for solving the zero Reynolds number linear ‘Stokes flow’ equations,
1.1
where ** u**,

*p*,

*μ*and

**are the velocity, pressure, kinematic viscosity and force per unit volume, respectively. These equations model flow strongly dominated by viscous effects with negligible inertia, valid for cilia and flagella since the Reynolds number can be estimated as**

*f**Re*=

*μ*

*U*

*L*/

*ρ*<10

^{−2}. Hancock (1953), working with Sir James Lighthill, introduced the singular ‘Stokeslet’ solution, corresponding to a concentrated point force, or equivalently the purely viscous component of the flow driven by a translating sphere. The Stokeslet is the solution of the Stokes flow equations for unit force acting in the

*j*-direction and concentrated at

**, corresponding to taking**

*ξ***(**

*f***)=**

*x**δ*(

**−**

*x***)**

*ξ*

*e*_{j}, where

*δ*(

**−**

*x***) denotes the Dirac delta distribution and**

*ξ*

*e*_{j}is the appropriate basis vector. The

*i*-component of the velocity field driven by this force is written as

*S*

_{ij}(

**,**

*x***), and in an infinite fluid takes the form 1.2 where**

*ξ**r*

_{i}=

*x*

_{i}−

*ξ*

_{i}and , with

*δ*

_{ij}denoting Kronecker delta tensor. The flow due to a force

**concentrated at the point**

*F***corresponds to taking**

*ξ***(**

*f***)=**

*x**δ*(

**−**

*x***)**

*ξ***in equation (1.1), the solution being given by**

*F**u*

_{i}(

**)=(1/8**

*x**π*

*μ*)

*S*

_{ij}(

**,**

*x***)**

*ξ**F*

_{j}. We assume the summation convention over repeated indices throughout.

The Stokeslet forms the basis for slender body theory for Stokes flow (Burgers 1938; Hancock 1953; Tuck 1964; Batchelor 1970; Johnson 1980), its local simplification ‘resistive force theory’, introduced by Gray & Hancock (1955) and the single-layer boundary integral method for Stokes flow (e.g. Pozrikidis 1992, 2002). These techniques have proved very useful in the modelling of cilia-driven flows (Blake 1972; Liron & Mochon 1976; Gueron & Liron 1992; Smith *et al.* 2007), sperm motility (Dresdner *et al.* 1980; Smith *et al.* 2009*b*), bacterial swimming (Higdon 1979*b*; Phan-Thien *et al.* 1987; Ramia et al. 1993) and atomic force microscopy (Clarke *et al.* 2006), allowing three-dimensional flow problems with moving boundaries to be solved with greatly reduced computational cost compared with direct numerical solutions of the differential equations (1.1). The fluid velocity field is modelled with an integral equation of the form
1.3
where *S* is a collection of surfaces or lines, e.g. representing boundaries, cell surfaces, cilia or flagella, with the symbol denoting line or surface integrals as appropriate, and ** f**(

**) denoting force per unit area or length, respectively. The force exerted by the body on the fluid is given by**

*ξ***(**

*f***)**

*ξ**d*

*S*

_{ξ}, and the force exerted by the fluid on the body is given by−

**(**

*f***)**

*ξ**d*

*S*

_{ξ}. The force density is calculated by combining equation (1.3) with either prescribed motions of the bodies, or a coupled elastic/active force model, through the no-slip, no-penetration condition .

## 2. The method of regularized Stokeslets

Line distributions of Stokeslets, e.g. equation (1.3) with the boundary *S* parametrized as ** ξ**(

*s*) for scaled arclength parameter 0<

*s*<1, representing objects such as cilia or flagella, have the principal disadvantage that the flow field is singular at any point

**=**

*x***(**

*ξ**s*). To calculate the force per unit length, it is necessary to perform collocation at points on the ‘surface’ of the slender body, displaced a small distance from the centreline

**(**

*X**s*

_{q})=

**(**

*ξ**s*

_{q})+

*a*(

*s*

_{q})

**(**

*n**s*

_{q}), where

*a*(

*s*

_{q}) is equivalent to the slender body radius and

**(**

*n**s*

_{q}) is a unit normal vector. The velocity field is regular in the region designated outside the body, with the singular line distribution lying inside the notional surface of the cilium.

Point distributions of Stokeslets , e.g. representing a swarm of cells, are also singular at any ** x**=

*ξ*_{q}. Surface distributions of Stokeslets, as used in the single-layer boundary integral method for Stokes flow, do not result in singular velocities, but nevertheless require careful numerical implementation to ensure that the surface integrals are calculated correctly (Pozrikidis 1992, 2002).

To circumvent these difficulties and ensure a regular flow field throughout the flow domain, Cortez (2001) introduced the ‘regularized Stokeslet’. This is defined as the exact solution to the Stokes flow equations with smoothed point forces,
2.1
The symbol *ψ*_{ϵ} denotes a cutoff function or ‘blob’ with regularization parameter *ϵ*, satisfying . Cortez *et al.* (2005) showed that with the choice , the regularized Stokeslet velocity tensor is given by
2.2
Here and in the rest of the paper, we use the compact notation . We follow Cortez *et al.* (2005) in using the solution , and its counterpart near a no-slip boundary given in equation (4.1), as the basis for our computational study.

Cortez *et al.* (2005) derived the equivalent Lorentz reciprocal relation and hence a boundary integral equation for the regularized Stokeslet method (RSM). This leads to the following equation for the fluid velocity at location ** x**, where the surfaces and lines representing the cells and beating appendages are denoted

*S*, 2.3 As above, the unknown

**(**

*f***) is a force per unit area or length depending on whether**

*ξ***is on a surface or line comprising**

*ξ**S*. A significant advantage of the RSM is that the kernel remains regular for

**=**

*x***, even when**

*ξ**S*consists of lines or points rather than surfaces, and so collocation can be performed at such points without the need to make a finite displacement. This method has proved effective in modelling a number of biological flow problems, including flow due to an individual cilium and the swarming of bacteria (Ainley

*et al.*2008), the bundling of bacterial flagella (Flores

*et al.*2005; Cisneros

*et al.*2008), hydrodynamic interaction of swimming cells (Cisneros

*et al.*2007) and the flagellar motility of human sperm (Gillies

*et al.*2009).

## 3. Alternative implementation of the RSM as a boundary element method

The mathematical problem considered in this paper is the calculation of the force density ** f**(

**) from known boundary velocity, given by the time derivative , through applying the no-slip condition to equation (2.3). The boundary, representing a cell surface, a flagellum, an array of cilia, a network of fibres or a collection of particles, is discretized as a set of points {**

*ξ*

*ξ*_{q}:

*q*=1,…,

*N*}. Applying collocation at these points gives the following equation: 3.1

The implementation of the regularized Stokeslet method used in Cortez *et al.* (2005), Ainley *et al.* (2008), Gillies *et al.* (2009) and Cisneros *et al.* (2007) is equivalent to approximating the force density by constant values, ** f**(

*ξ*_{r})≈

*f*_{r}and replacing the integral of the Stokeslet with a low-order quadrature, giving 3.2 where {

*w*

_{1},…,

*w*

_{N}} are quadrature weights. Hence the matrix equation

*A*

*X*=

*B*can be formed, where 3.3

This discretization can be viewed as a constant-force implementation of a boundary integral method, using a low-order quadrature in which the abscissae are identified with the collocation points. However, while the force density ** f**(

**) varies relatively slowly, the kernel**

*ξ***S**

^{ϵ}(

*ξ*_{q},

**) will vary rapidly in the neighbourhood of**

*ξ*

*ξ*_{q}. Hence if the surface is discretized as patches

*S*

_{1},…,

*S*

_{N}, the term

*w*

_{r}

*f*_{r}⋅

**S**

^{ϵ}(

*ξ*_{q},

*ξ*_{r}) in equation (3.2) may not be an accurate approximation to , particularly when

*r*=

*q*.

Consequently, to obtain accurate solutions, it is necessary to use a relatively large number of nodes *ξ*_{q}, many degrees of freedom, and hence a large dense matrix system to calculate accurate solutions, as discussed in detail in appendix C*b*. This is an unnecessary computational effort, since even with a constant-force discretization, far fewer degrees of freedom for ** f** are required than quadrature nodes for the evaluation of the kernel integrals. It is also necessary to employ an empirical rule of the form

*ϵ*=

*C*(

*δ*

*s*)

^{m}relating the regularization parameter to the mesh size

*δ*

*s*to obtain the expected result, as discussed in Cortez

*et al.*(2005) and Ainley

*et al.*2008.

For this reason, we instead implement a constant-force boundary element method where the integration and force discretization are ‘decoupled’, which we then apply to a number of example problems. On each patch *S*_{r}, the force is approximated by *f*_{r}. Then we have
3.4
This problem can again be written as a matrix equation *A**X*=*B*, with *A* now taking the form
3.5
The *ϵ*–*δ**s* coupling of the original method is now substantially weakened, for the following reason: the discretization chosen for the force no longer has any bearing on the evaluation of the kernel integrals in equation (3.5). If numerical integration is used, it is still necessary to choose a sufficiently refined method for the value of *ϵ* chosen, for the ‘near-singular’ integrals. However, for the values of *ϵ* typically used, this is far less computationally costly than increasing the number of degrees of freedom, as shown in appendix C, and moreover may be accomplished in principle with an automatic integration routine or in some cases analytic integration.

The choice of constant-force elements is purely for simplicity: higher-order methods may also be used, generalizing the above so that ** f**(

*s*) is approximated as , where {

*ϕ*

_{1}(

*s*),…,

*ϕ*

_{N}(

*s*)} is a set of basis functions. Cubic Lagrange interpolating polynomials were used by Smith

*et al.*(2009

*a*) in the context of the non-regularized viscoelastic Stokeslet distributions. The integral may be evaluated using any appropriate means: Gauss–Legendre numerical quadrature, a more advanced adaptive quadrature, or in the case of problems where the boundary is reduced to a line segment, the integrals may be performed analytically (see appendix B). In appendix C, we compare the two approaches (3.2) and (3.4).

## 4. The Ainley *et al.* image system for a no-slip plane boundary

To model fluid propulsion by cilia protruding from a cell surface, Blake (1971) derived the Stokeslet image system satisfying the no-slip, no-penetration condition ** u**=0 on

*x*

_{ 3}=0. The technique is equally useful in modelling the interaction of a glass surface and a sperm cell and removes the need for the cell surface boundary to be included in the boundary integral equation. We shall not repeat the solution here, however it consists of a Stokeslet in the fluid, an equal and opposite image Stokeslet, a Stokes-dipole and potential dipole. As an important step in generalizing the regularized Stokeslet method to the modelling of biological flows, Ainley

*et al.*(2008) derived the equivalent regularized image system, which we denote for flow near a no-slip boundary. We rewrite the solution in index notation for direct comparison with the result of Blake (1971): 4.1 The tensor

*Δ*

_{jk}takes value +1 for

*j*=

*k*=1,2, value −1 for

*j*=

*k*=3 and zero otherwise, originally written as (

*δ*

_{jα}

*δ*

_{αk}−

*δ*

_{j3}

*δ*

_{3k}) by Blake. The last line of equation (4.1) is not precisely equivalent to that in Ainley

*et al.*(2008), the latter having a typographical sign error. The term is a more slowly decaying blob than

*ψ*

_{ϵ}(

*R*), discovered by Ainley

*et al.*(2008) as the function generating the potential dipole appropriate to the no-slip image system. The additional term on the final line is the difference between rotlets arising from the two different blobs

*ψ*

_{ϵ}(

*R*) and

*ϕ*

_{ϵ}(

*R*), and is O(

*ϵ*

^{2}/

*R*

^{4}).

## 5. Application to biological flows

The numerical schemes given by equations (3.3) and (3.5) are applied to test problems involving rod and sphere translation in an infinite fluid in appendix C. In this section, we apply the method to two biologically inspired problems: determining the flow generated by an array of cilia protruding from a flat surface and by a choanoflagellate cell in an infinite fluid surrounded by a silica lorica structure. We non-dimensionalize with scalings *L*, 1/*σ* and *σ**L* for length, time and velocity, respectively, where *L* is the cilium or flagellum length and *σ* is the radian beat frequency. The appropriate scaling for force per unit length on a filament is *μ**σ**L*, and for force per unit area on a surface is *μ**σ*.

### (a) Modelling a patch of cilia

Cilia lining the airway epithelia and the surfaces of micro-organisms typically appear in dense ‘fields’, while cilia within the female reproductive tract generally appear in patches, with ciliated cells being interspersed between secretory cells. The density of cilia in such systems can be very high, approximately 6–8 *μ*m^{−2} (Sleigh *et al.* 1988), which for a square lattice is equivalent to a spacing of 0.35–0.41 *μ*m, or approximately 0.06–0.07*L* for a cilium length *L*=6 *μ*m. This contrasts with embryonic nodal primary cilia, for which there is only one cilium per cell.

Slender body theory does not, in general, produce exact solutions, except in certain special cases such as that of a straight ellipsoids in a uniform flow, or uniform shear (Chwang & Wu 1975). When modelling single bodies with non-uniform velocity or curvature, or multiple bodies, the flow field ** u**(

**,**

*x**t*

_{n}) generated by a numerical solution for

**(**

*f**s*,

*t*

_{n}) typically exhibits errors when comparing

**(**

*u***(**

*ξ**s*,

*t*

_{n})) with the prescribed boundary movement for values of

*s*that were not used as collocation points. These errors can be reduced in magnitude through mesh refinement, using higher-order singularities (Smith

*et al.*2007), or higher-order basis functions (Smith

*et al.*2009

*a*), however they cannot be eliminated entirely at the ends of the body without resorting to replacing the line distribution by a surface distribution, which incurs considerably higher computational expense. In an earlier computational study using singular Stokeslet distributions, Smith

*et al.*(2007) found that cilia could not be arranged in arrays with a spacing of less than 0.3

*L*, equivalent to approximately 1.8

*μ*m, without significant numerical errors occurring in the computational solution—a principal difficulty being that collocation points on the cilia may be close together, particularly during the recovery stroke. This limited the applicability of the technique in determining the effect of cilia density at realistic concentrations.

In this study, we test the RSM without the use of higher-order singularities and without the use of surface distributions, to model a densely packed patch of cilia. The array is modelled as a square lattice of *M*=*M*_{1}×*M*_{2} cilia, for 1≤*l*≤*M*_{1} and 1≤*m*≤*M*_{2}; at time *t*, the velocity field is given by
5.1
where the vector denotes the force on the *r*th segment of the cilium in the (*l*,*m*)-position in the array, at time *t*.

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). The mathematical specification has the following form:
5.2
where the functions *α*_{p}(*s*) and **β**_{p}(*s*) are cubic polynomials, the coefficients being given in Fulford & Blake (1986). Rather than imposing a metachronal wave on the cilia, we specify that they beat in synchrony, since fully developed metachronal waves are generally more typical of larger cilia arrays. Hence the position vector of the cilium in the (*l*,*m*)-position is simply
5.3
where the parameter *a* is the cilia spacing.

In the numerical implementation, the flow field is calculated at *t*_{n}=*n**δ**t*, where *n*=1,…,*N*_{T} and *δ**t*=2*π*/*N*_{T}. We find that results for the volume flow rate converge to three significant figures using *N*_{T}=40. As discussed in appendix B and appendix C*a*, we choose the regularization parameter based on the cilium radius-to-length ratio, so that *ϵ*=0.1/6=0.0167. Results are calculated using the University of Birmingham’s BlueBEAR Opteron cluster, however the most computationally expensive simulation of 169 cilia requires only 8 h 20 min of CPU time, and so could be performed on a desktop PC.

To monitor the accuracy of the numerical solution for closely approaching cilia, we use a similar technique to Gillies *et al.* (2009) and evaluate the post hoc collocation error at each timestep *t*=*t*_{n}, denoted . This is given by
5.4
where the integral is calculated numerically using the points *s*_{p}=(*p*−1/2)/120 for 1≤*p*≤120. For all results presented, we used *N*=20 elements per cilium and verified that the error did not exceed on any cilium (*l*,*m*) or any timestep *n*. Example results showing the post hoc comparison are given in figure 1. The discretized matrix *A* has similar favourable properties to the matrix resulting from the discretization using a singular kernel, with the diagonal entries being relatively large. The condition numbers calculated for the multi-cilia results were approximately 1.3×10^{−3}.

Figure 2 shows examples of the three-dimensional flow field computed with the model at four points in the beat cycle. The rapid decay of the flow magnitude in both horizontal and vertical directions is evident. The tendency of a patch of cilia to draw fluid ‘in’ from each side to the rear, and push fluid ‘out’ to the side in front is evident, as may be anticipated from the far-field form of the image system (Blake 1971).

An important function of cilia is the transport of liquid, and it is of interest to determine how the variation in cilia density observed in different biological systems might affect transport. To quantify this, we calculated the mean volume flow rate , using the formula described in appendix A. To examine the different fluid dynamic mechanisms at work, we simulate different cilia densities for flow driven by a finite patch in two ways: (a) keeping the number of cilia constant and varying the spacing through the size of the patch, (b) keeping the size of the patch constant, and varying the spacing through the number of cilia. Results are shown in figure 3.

For a 3×3 cilia patch, as the cilia are brought closer together, they hydrodynamically interact, reducing the fluid transport (figure 3*a*). Because of this interaction effect, for a patch of fixed dimensions, increasing the cilia number from 3×3 to 13×13 results only in a ≃50 per cent increase in fluid transport, and moreover 9×9 cilia give approximately the same transport as a 13×13 array (figure 3*b*). The common point in figure 3*a*,*b* corresponds to cilia spacing 0.5*L*, for which the non-dimensional mean volume flow rate is approximately 1.45. These results were computed with a prescribed beat pattern and frequency—in the biological system, these will vary depending on the flow field produced by the other cilia, as discussed in §6.

### (b) The modelling of a choanoflagellate

Marine choanoflagellates are the nearest unicellular relatives of the animal kingdom and exhibit a very diverse array of structures and behaviours that have yet to be fully explained biologically. Their feeding and swimming behaviour involves the beating of a single flagellum and has been the subject of previous fluid mechanical modelling using slender body theory (Orme *et al*. 2003). For recent review of the biology, see Leadbeater *et al*. (2009), and for details on flagellar movement and fluid transport, see Pettitt *et al*. (2002). Certain species exhibit a silica surrounding basket-like structure, known as a *lorica*, made up of a network of silica strips, termed *costae*. The lorica has been hypothesized to perform a number of roles including protecting the cell, reducing the translational velocity of free-swimming cells and directing feeding currents to its mouth. Simulating these behaviours is a subject for considerable future study; in this paper, we simply demonstrate that with the boundary element regularized Stokeslet method it is computationally inexpensive to simulate such flows and we make an initial study of the effect of surrounding the cell with a simple lorica structure.

Both the flagellum and the 16 costae are modelled as filaments, the cell body being modelled as a sphere with 6×4×4 elements, with *N*=30 constant-force elements for the flagellum and *N*=8 constant-force elements for each costa. The costae were modelled less precisely than the flagellum, since their role is simply to provide a hydrodynamic drag. For these simulations, we choose *ϵ*=0.01, corresponding to a flagellum aspect ratio of 0.1 *μ*m:10 *μ*m. The same regularization parameter was also used for the sphere surface, since this has been verified to give accurate results, as shown in appendix C, and for the costae. Simulations for a single timestep were performed on a desktop PC, requiring under 10 min of CPU time.

Our simplified model for the choanoflagellate flagellum, based on that in Higdon (1979*b*), consists of a sphere of radius *α*=0.15, located at (0,0,*h*_{0}), where *h*_{0}=0.4. The flagellum is described by
5.5
where *t* is scaled with respect to inverse radian frequency and is scaled with respect to flagellum extent rather than length. The wavenumber *κ*=2*π*, so that the wavelength is unity. The envelope parameter , where *κ*_{E}=*κ*. The coordinate is not an arclength parametrization, and so it is necessary to make a change of variable to arclength parameter *s* to calculate the force per unit length shown in figure 6.

The silica lorica is constructed using a simplified model, based on the more complex family of geometric models of Leadbeater *et al*. (2009), in which it was shown that there are two sets of costae that run vertically or spiral around the axis of the cell, creating a basket-like structure. We defined a family of *N*_{C} vertical strips, given by for *m*=1,…*N*_{C},
5.6

The family of *N*_{C} spiralling strips, denoted for *m*=1,…*N*_{C}, was then generated by changing the argument of the trigonometric functions in equation (5.6) to . For this study, we chose the envelope function , and a total of 2*N*_{C}=16 costae.

Figure 4 shows the instantaneous fluid velocity field, computed on horizontal and vertical surfaces, both with and without the lorica. While the flow close to the flagellum is relatively unchanged by the presence of the lorica, it is apparent from figure 4*c*,*d* that the remainder of the flow field is strongly suppressed in magnitude. Figure 5 shows the vertical component of the velocity field *u*_{3} only, on a horizontal cross section. The presence of the lorica reduces this component in magnitude by nearly 50 per cent within the lorica interior, and again almost completely abolishes any external flow. By contrast, the hydrodynamic interaction of the flagellum and lorica is relatively weak, with figure 6 showing that the force density components are almost unchanged by its presence.

## 6. Discussion and future work

A simple modification of the regularized Stokeslet method was presented, based on decoupling the force discretization and boundary integration. This reduced the number of degrees of freedom required to compute accurate solutions and removed the need for an empirical parameter relating regularization parameter and the discretization mesh for the force. The implementation was designed to be the most elementary possible, to emphasize its ease of use in biological fluid mechanics problems. For the Stokes’ law test example, a considerable reduction in computational cost was obtained, with approximately 1000-fold fewer kernel evaluations being required to obtain accuracy to within 1 per cent, and a 144-fold reduction in the number of degrees of freedom *N*. For problems involving *M* bodies, the storage and linear solver costs are at least O(*N*^{2}*M*^{2}), hence a reduction in the number of degrees of freedom for many-body problems will be very useful for the modelling of complex biological flows.

The regularization parameter *ϵ* has two distinct roles in models that employ surface distributions and line distributions. For surface distributions, *ϵ* plays the role of a tunable parameter, with smaller values giving more accurate solutions at the expense of more abscissae being required if numerical quadrature is used. It is sufficient to choose *ϵ* small enough that the regularization error is negligible; Cortez *et al*. (2005) showed for the blob function used in this paper that this error is O(*ϵ*^{2}), and in this study we found that for the test case of a translating sphere, a value of *ϵ*=0.01 gave errors of less than 1 per cent with 162 degrees of freedom. For line distributions, *ϵ* is instead equivalent to a physical parameter of the filament—the slenderness ratio, as shown in appendix B. A detailed asymptotic analysis of the slender body theory shall follow this paper.

The technique is capable of providing accurate results for denser arrays of cilia than have previously been possible, without recourse to higher-order solutions or basis functions. Arrays of up to 169 cilia, with spacing 0.08 times their length, could be simulated for a 40-timestep cycle in approximately 8 h with a single CPU. It was found that for prescribed cilia movement, cilia density of less than 0.125 times their length does not produce significant increases in transport, due to hydrodynamic interaction obliterating any gains resulting from there being more cilia—it is likely that the role of realistically high cilia densities will only be fully understood through more complex fluid–structure interaction models which take into account more biophysical details of the system in question.

For discussion of the role of cilia in reproduction, see Fauci & Dillon (2006). An additional complication when considering airway cilia is the effect of the serous–mucous bilayer, and possible pressure gradients that may exist due to surface tension, osmosis or both. Pressure gradients will have important effects on fluid transport (Mitran 2007; Smith *et al*. 2008*b*), and cilia penetration of the mucous layer may allow continued transport when many cilia are inactive (Fulford & Blake 1986).

The technique was also used to simulate the effect of the complex silica lorica structure of a marine choanoflagellate on the flow field and force density, with minimal computational cost. The lorica does not significantly affect the force distribution on the beating flagellum, but does significantly reduce the far-field flow produced. Future models that consider force and torque balance for a free-swimming cell, in addition to ambient currents, and other features of the cell geometry such as the collar, may provide insight into the role and function of loricae, and the diversity that exists in choanoflagellate species.

In line distribution representations, the force density depends on the regularization parameter—in particular, it can be estimated as as . For free-swimming problems, e.g. the work of Gillies *et al*. (2009), the swimming speed depends on force and torque balance on the force density, and hence on the value chosen for the regularization parameter. As shown in appendix C*a*, this parameter may be chosen based on the slenderness ratio. This observation is relevant only to line distribution models of flagella and not to surface distributions.

The calculations using line distributions were less accurate than those for translation of a sphere, however the method performed with similar accuracy to our previously published studies for cilia motion, and superior accuracy for dense arrays. Future work may involve the use of higher-order solutions, higher-order basis functions and non-uniform meshes to reduce the end error for slender bodies. Recently, Smith *et al*. (2009*a*) presented a slender body theory for modelling linear viscoelastic flow, which exploits computational methods developed using singular solutions of the Stokes flow equations. A limitation is that errors in the slender body velocity field are propagated through the solution, which can lead to numerical instability if the Deborah number is larger than approximately 0.06 for the problem considered. The RSM, possibly combined with higher-order Green’s functions, may provide a more accurate solution, particularly for problems with multiple cilia, and hence may extend the range of Deborah number that may be modelled.

Other researchers have regularized the slender body theory equations differently: the ‘Lighthill theorem’ (Lighthill 1976), later exploited by Dresdner *et al*. (1980) and by Gueron & Liron (1992); Gueron & Levit-Gurevich (1999) in the form of the Lighthill–Gueron–Liron theorem, are based on singular Stokeslet distributions, but employ integration to remove the near-field singularity, so that the local part of the integral equation is replaced by resistance coefficients. A related but different approach is evident in the work of Tornberg & Shelley (2004) on the modelling of complex suspensions of flexible fibres. These techniques have proved very successful in providing an efficient fluid mechanics solver as part of a model of cilia or fibre interactions with the ambient fluid (Gueron & Levit-Gurevich 1999, 2001; Tornberg & Shelley 2004). An advantage of the regularized Stokeslet method is that the modelling of particles, filaments and surfaces is unified within a single mathematical framework. This study is intended to expand the scope of the technique to problems that may previously have been computationally difficult, while retaining its ease of implementation. It should be noted that with sophisticated modern parallel computation, direct solution of the fluid flow and solid mechanics equations for large arrays of cilia is possible (e.g. Mitran 2007).

In this study, we considered the flow driven by prescribed cilia and flagellar motions. To understand and interpret phenomena such as the effects of viscosity, cilia coordination and cooperation, axoneme ultrastructure, accessory fibres and their effect on beat pattern and frequency, it is necessary to couple the fluid mechanics to models of elasticity and active bending moment generation. Important work has already been done on this by a number of authors, including those discussed above. The framework presented here may assist with the associated fluid mechanics computations, particularly for problems with many bodies, or where there is a complex geometry to be discretized, e.g. the surface of the ciliated ventral node of the developing embryo (Smith *et al*. 2008*a*; Hirokawa *et al*. 2009), or the interaction of sperm cells with the female reproductive tract.

## Acknowledgements

I thank my colleagues, particularly Professor John Blake, who introduced me to the biological fluid mechanics of cilia and flagella, and along with Dr Eamonn Gaffney continue to provide invaluable research training, mentoring and guidance. I also thank Dr Jackson Kirkman-Brown (Centre for Human Reproductive Science, Birmingham Women’s Hospital and University of Birmingham) for valuable guidance, and for the opportunity to work with real cilia and flagella first-hand, and Mr Hermes Gadêlha and Mr Henry Shum (University of Oxford) for fascinating conversations. I also thank Dr Dov Stekel and Dr Barry Leadbeater (University of Birmingham) for discussions on choanoflagellates. This work was funded by the MRC (Special Training Fellowship G0600178).

## Appendix A.

### The volume flow rate due to a point force in the vicinity of a no-slip boundary

The pumping effect of a cilium in the *x*_{1}-direction can be quantified using the volume flow rate *Q*,
A1
the choice of *x*_{1} being arbitrary due to conservation of mass. The mean volume flow rate over one beat cycle is denoted by .

The volume flow rate due to a Stokeslet singularity *B*_{1j} near to a no-slip surface was given by Liron (1978). Our definition of the Stokeslet differs in that we do not include the leading 1/8*π**μ* factor, so in our notation we have
A2
where *x*_{1}, *ξ*_{1} and *ξ*_{2} are arbitrary.

By conservation of mass, the volume flow rate does not depend on *x*_{1}. Taking *x*_{1} and hence *r* to be large, we note that the regularization error , and the regularization error for the other images decays at least as rapidly as O(*ϵ*^{2}/*r*^{4}). Hence the regularization error , the double integral of which tends to zero as , and it cannot make any contribution to the volume flow rate. This proves that result (A2) holds exactly for the regularized solution , and so combining equations (5.1), (A1) and (A2) we have
A3
In §5*a*, we report results in non-dimensional units, with *Q* scaled with respect to *σ**L*^{3}.

## Appendix B.

### Straight-line integrals of the regularized Stokeslet and their relation to slender body theory

For filaments composed of straight-line segments, e.g. the three-link swimmer (Purcell 1977; Becker *et al*. 2003; Tam & Hosoi 2007), equation (3.5) reduces to the calculation of straight-line integrals. Similarly, provided the discretization is sufficiently fine, curved line integrals can be approximated by a sum of straight-line integrals. These integrals can be calculated explicitly and used in a similar manner to the singular Stokeslet-dipole integrals employed by Higdon (1979*a*). While this is not required to implement our computational method, in certain situations it is more efficient than using numerical quadrature, e.g. the results shown in appendix C*a* and §5*a*. As the radius of curvature of the relative filament length decreases, it is necessary to subdivide the elements into multiple straight-line segments, and the computational savings are less. For this reason, we use numerical integrals for the choanoflagellate example of §6.

For an element centred on ** ξ**(

*s*

_{r}), having tangent

**′(**

*ξ**s*

_{r}), the approximating straight-line segment is

**(**

*ξ**s*)≈(

*s*−

*s*

_{r})

**′(**

*ξ**s*

_{r})+

**(**

*ξ**s*

_{r}). The regularized Stokeslet integrals are of the form B1 We perform the coordinate transformation , where

**Θ**is a rotation matrix. The first row (Θ

_{11},Θ

_{12},Θ

_{13}) is given by the tangent

**′(**

*ξ**s*

_{r}), with the second and third rows completing an orthogonal basis. The integral becomes, in ‘local’ coordinates, B2 The integral in the original coordinate system is given by . The coordinate translation does not affect

**D**since the integrand

**S**depends on

**and**

*x***only through their difference.**

*ξ*The definite integrals in the local coordinate system are , where are the corresponding indefinite integrals. In concise notation (*x*,*y*,*z*) for , these can be written as
B3
with .

The diagonal matrix entries, occurring when ** x**=

**(**

*ξ**s*

_{r}), require higher-order numerical quadrature than the other entries since the kernel is nearly singular at this point, and hence are particularly suited to analytic evaluation. These ‘local’ integrals are , the components of which do not depend on

**or**

*x**s*

_{r}. They have explicit form B4 The local integrals also provide a link from the numerical implementation to classical slender body theory. Taking

*ϵ*as a proxy for slender body radius and assuming that the length scale on which the force density varies, and the radius of curvature, are small in comparison with

*ϵ*, then it is reasonable to take

*ϵ*≪

*δ*

*s*≪1. The ‘local’ integrals then give the dominant contribution of the Stokeslet distribution in the calculation of the force density, corresponding to the resistance coefficients of Gray and Hancock theory. Moreover, rescaling , we have the asymptotic dependence on for , as found by Batchelor (1970).

## Appendix C.

### Test problems

Two test problems are considered: the translation of a fibre normal to its axis and the translation of a sphere, both in an infinite fluid. The latter problem was previously considered extensively by Cortez *et al*. (2005).

#### (a) Test case 1: force on an axially translating slender body

A classical slender body theory result (e.g. Batchelor 1970, noting the slightly different notation) is that the magnitude of the force on a rod of length *l* and radius *a*, translating with velocity *u* parallel to its axis is, to leading order . Non-dimensionalizing with force scale *μ**u**l*, we calculate a numerical solution using a line distribution of regularized Stokeslets, with the regularization parameter chosen as *ϵ*=*a*/*l*. We compare the two numerical implementations given in equations (3.2) and (3.4) with *a*/*l*=0.01 and examine the force at the midpoint of the rod. For the latter case, integrals are calculated analytically, as described in appendix B. The ‘boundary element’ implementation (3.4) gives a relative error of less than 3.5 per cent compared with the asymptotic estimate for *N*=5,10,20,40,80,160. The original implementation (3.2) gives relative error of over 35 per cent for *N*=10, and it is necessary to use *N*=80 degrees of freedom to give similar accuracy to the boundary element version.

Benchmarking the two algorithms on a desktop PC with *N*=20 for the boundary element version and *N*=80 for the original implementation, the former requires less than half the CPU time, using analytic kernel integration. If we repeat the calculation with *M*=5 equally spaced parallel rods, the boundary element version requires approximately one-sixtieth of the CPU time compared with the original implementation. The reason for the significant saving is that the matrix setup cost is O(*M*^{2}*N*^{2}) and the direct solution cost is O(*M*^{3}*N*^{3}).

#### (b) Test case 2: force on a translating sphere

To test the method with a surface distribution, we compute the drag force on a translating sphere in an infinite fluid. Stokes’ law gives this as exactly *F*=6*π**μ**u**a*, where *a* is the sphere radius. Non-dimensionalizing with force scale *μ**u**a*, we compare the boundary element method (3.4) with the results of Cortez *et al*. (2005). We use the same geometric discretization, projecting the faces of a cube onto its surface and subdividing these faces with a square mesh, as shown in figure 7. The square mesh provides a natural parametrization of each element, which we denote (*α*,β). For the boundary element method, numerical integration is performed using two-dimensional Gauss–Legendre quadrature with the parametrized coordinates, taking 12×12 points for the ‘near-singular’ cases *q*=*r* in equation (3.4) and 4×4 points otherwise. The surface metric is calculated as the magnitude of the vector product |*∂*** ξ**/

*∂*

*α*∧

*∂*

**/**

*ξ**∂*β|, as described in Pozrikidis (2002). All of the collocation and quadrature points lie on the curved surface of the sphere.

For problems involving surface integrals, refining the mesh by halving the element width results in four times the number of degrees of freedom, 16 times the matrix setup and storage cost, and at least the same cost increase for an iterative linear solver. With a conventional direct solver, the cost increases with O(*N*^{3}), i.e. 64 times. For this reason, a method that gives accurate results with a relatively coarse mesh is beneficial. Consistent with Cortez *et al.* (2005), we compute solutions using the generalized minimum residual method with zero initial guess, although for the cases *N*≤6×12×12; solution by *L**U* factorization is possible in a few seconds on a desktop PC.

Table 1 is based on the results with regularization parameter *ϵ*=0.01 reported by Cortez *et al.* (2005) using equation (3.2), although we include relative percentage error and the number of scalar kernel evaluations necessary to construct the matrix *A*. To obtain a relative error of less than 1 per cent, it is necessary to use a mesh with 6×36×36 patches, giving 23 328 scalar degrees of freedom for the unknown force components, the matrix setup requiring 5.44×10^{8} kernel evaluations.

Table 1 shows results calculated using equation (3.4), with coarser meshes and fewer degrees of freedom. To obtain similar accuracy, a 6×3×3 mesh with 162 degrees of freedom is sufficient. While the matrix setup cost per DOF is greater due to the higher-order quadrature, the greatly reduced degrees of freedom results in the number of kernel evaluations being three orders of magnitude smaller. The matrix solution cost is reduced similarly.

N | DOF | matrix entries | no. kernel evaluations | rel. % error in drag |
---|---|---|---|---|

Original method | ||||

6×12×12 | 2592 | 6.72×10^{6} | 6.72×10^{6} | 12.6 |

6×24×24 | 10368 | 1.07×10^{8} | 1.07×10^{8} | 2.76 |

6×36×36 | 23328 | 5.44×10^{8} | 5.44×10^{8} | 0.849 |

6×48×48 | 41472 | 1.72×10^{9} | 1.72×10^{9} | 0.265 |

Boundary element method | ||||

6×3×3 | 162 | 26244 | 482112 | 0.827 |

6×4×4 | 288 | 82944 | 1.44×10^{6} | 0.626 |

6×6×6 | 648 | 419904 | 6.97×10^{6} | 0.431 |

6×9×9 | 1458 | 2.12×10^{6} | 3.48×10^{7} | 0.320 |

6×12×12 | 2592 | 6.72×10^{6} | 1.08×10^{8} | 0.279 |

The effect of varying the regularization parameter *ϵ* on the predicted drag was also examined. Using *N*=6×4×4 and constant force discretization, the relative error for *ϵ*=0.05 was 1.4 per cent, reducing to 0.63 per cent for *ϵ*=0.01. The error remains at approximately 0.62 per cent for *ϵ*=0.005 and 0.0025, implying that this is a satisfactory value—the remaining error arises from the force discretization. Cortez *et al.* (2005) showed that the error associated with the regularization parameter can be estimated as O(*ϵ*^{2}). Hence a value of *ϵ* that gives accurate results for the sphere problem will, in general, give accurate results for any solid surface for which the length scales are significantly larger than *ϵ*. Moreover, once *ϵ* has been chosen sufficiently small, accurate results may be computed for any values of *ϵ* that are smaller, provided that the kernel quadrature is computed accurately.

## Footnotes

- Received June 2, 2009.
- Accepted August 26, 2009.

- © 2009 The Royal Society