Functional geometry of ciliated tentacular arrays in active suspension feeders

Ciliary pumping with arrays of ciliated tentacles is one of the most universal mechanisms of water movement used by aquatic animals in generating suspension-feeding and respiratory flows. Examples of such arrays are found in bivalves, ascidians, brachiopods, bryozoans, phoronids, polychaetes and larval echinoderms (Strathmann, 1973; Jørgensen, 1989; Hart, 1991, 1994). The effectiveness of ciliary pumping with tentacular arrays in this large and diverse set of aquatic animals has far-reaching ecological and evolutionary consequences. Morphological details of the tentacle arrays vary across these taxonomic groups, but a basic geometrical design occurs repeatedly: a series of roughly cylindrical tentacular structures, arranged approximately parallel to each other in a plane normal to the flow, with propulsive cilia on the lateral surfaces where the gaps between neighboring tentacles are narrowest.

Repeated evolutionarily distinct occurrences of parallel tentacular structures with lateral cilia suggest that this geometry may be one of only a few effective design solutions to a common evolutionary problem: how to transport and exploit large volumes of water using cilia that are constrained to be small and that operate over relatively small distances and low speeds. The convergent tentacular geometry further suggests that similar hydromechanical principles may govern the functional morphology of ciliary pumping across taxa. In each example, the volume of water moved per unit time is determined largely by a balance between viscous traction forces transmitted from the cilia, losses due to flow resistance through the tentacular array and, in some cases, a pressure drop across the array opposing the flow, imposed by flow resistance in other components of the water-processing system (such as siphons in bivalves and ascidians).

However, the quantitative effects of variations in morphology in ciliated tentacle arrays are difficult to anticipate without detailed empirical and theoretical analysis. For example, increasing the spacing between adjacent tentacles has competing effects, both increasing the amount of fluid ‘dragged’ through the array by the cilia and increasing the ‘leakage’ through the array caused by an opposing pressure gradient. Thus, the effects of spacing on performance criteria such as flow rate can only be determined in the context of the water-processing system as a functional unit (LaBarbera, 1990). A large literature reporting observations of natural tentacle arrays and a smaller amount of hydrodynamic theory (Jørgensen et al. 1988; Silvester, 1988; Sleigh, 1989; Nielsen et al. 1993) and physical modeling (Emlet, 1991; Hansen and Tiselius, 1992) have addressed the functional morphology of ciliated tentacular arrays, how the biomechanics of pumping are reflected at the organismal level in clearance and respiration rates (Meyhöfer, 1985; Silverman et al. 1995) and, ultimately, how filtration mechanics help to determine the fluxes of energy and nutrients in communities of suspension feeders (Eckman and Duggins, 1993; Butman et al. 1994; Loo and Rosenberg, 1996). These studies are particularly well-developed in bivalves, especially the blue mussel Mytilus edulis (Jones et al. 1990, 1992; Beninger et al. 1992; Nielsen et al. 1993). However, no theoretical or experimental study has systematically examined how pumping performance varies with geometrical and hydrodynamic variables in a way that specifies design trade-offs from an evolutionary and ecological perspective. Such a systematic study is a requirement for making quantitative predictions about performance-maximizing geometries.

In this study, we investigate one aspect of performance in ciliary pumping systems by asking the following question: for a given set of morphological constraints and hydrodynamic conditions, how should a tentacle array be designed so as to maximize flow rate? We address this question using a computational hydrodynamic model that calculates flow rates through ciliated tentacle arrays. Using this model, we systematically examine how pumping performance varies with geometrical variables, such as gaps between adjacent ciliary bands and the pressure differential across the array, and look for general principles of favorable design for tentacular arrays under different hydrodynamic conditions. As a case study, we apply the model to ciliary pumping in Mytilus edulis and compare the predictions of the model for geometries that maximize flow rates with experimental observations from the literature.

We describe the geometry of the ciliated tentacular array and discuss some basic allometric and scaling properties of flow through such arrays. We then present simulation results illustrating the variation of pumping performance under various geometries and operating conditions, and present a scheme based on scaling laws and interpolation of our numerical results to predict the flow rates through a wide range of tentacle morphologies. We calculate the flow rates in M. edulis and explore performance trade-offs in geometrical variables for this case study. We discuss some implications of our results for the functional morphology of suspension feeders in general. Finally, we give a brief description of the fluid mechanical simulation methodology in the Appendix.

Our simulations used a numerical technique known as the immersed boundary method, which has been used successfully in a number of biological fluid dynamics problems that involve unusual geometries and material properties. We use this method to calculate the two-dimensional flows generated through tentacular arrays by the movement of cilia at low to moderate Reynolds numbers (Re). A brief description of this method is given in the Appendix; see the references cited there for additional details.

The geometry of the tentacular array in our simulations is shown schematically in Fig. 1. The model assumes an infinite array of identical tentacles, each having three cilia on each lateral surface. The cilia undergo an angular displacement during their power strokes of 90 ° at a constant tip speed, propelling the surrounding fluid. We do not model the ciliary recovery stroke, during which relatively little water movement occurs; instead, cilia at the end of their power strokes are instantaneously moved to their starting positions and begin another power stroke. Cilia are impermeable and have no-slip boundary conditions (meaning that water in contact with a solid surface has the same velocity as that surface). Examples of flows generated by the model are shown in Fig. 2. In our simulations, both the cilia and the tentacles have compliance, that is, they bend and stretch in response to hydrodynamic forces. This distortion is visible in the cilia (Fig. 2B) but is negligible in the tentacles. Compliance may be an essential element of particle detection in suspension feeding flows; however, we do not consider this form of hydromechanical signaling in this paper.

where L is the tentacular spacing and s is the tip speed of the cilia. [Note that all flow rates cited in this paper are per unit length axially along the tentacles (i.e. out of the page in Fig. 1). Thus, Q has the dimensions m² s⁻¹. A summary of dimensional and non-dimensional variables is given in Table 1.] By using Q_n, the effects of isometric size and speed changes on Q are factored out.

We present our results in three stages. First, we show some numerical results of our simulations, focusing on allometry and the effects of geometrical variations. We then outline a scaling and interpolation scheme that allows the reader to estimate flows through any array geometry within the ranges covered by our simulations (this covers most natural examples of pumping with ciliated tentacular arrays of which we are aware). Finally, we use this scheme to assess pumping performance in the blue mussel Mytilus edulis.

Some basic allometric properties of flow through the ciliated arrays are shown in Fig. 3. Fig. 3A presents the dimensional (unscaled) flow rates between each pair of tentacles, for a standard tentacle geometry (variables are given in the figure legend) with variation in the speed of ciliary motion and in the diameter D of the tentacle (expressed as the tentacle diameter ratio r=D/c). When cilia speed and length and the gap between adjacent cilia tips are fixed, the tentacle diameter has a relatively small effect on the flow. Interestingly, the effect of increasing tentacle diameter is to increase flow rate. However, when the flow is normalized by the length of the array (as in Q_n), this slight increase in flow per tentacular unit is more than offset by the larger size of each unit (Fig. 3B). Thus, an array of fixed total width pumps less water with increasing tentacle diameter, because fewer tentacles fit into it. The scaled flow, Q_s, shows that the flow through these arrays slightly exceeds that predicted by the trapezoidal profile, a result that we found to be characteristic of relatively small cilia gaps and relatively high Re values.

Of particular interest in Fig. 3B is the weak variation of Q_s and Q_n over a relatively large range of Re. These results suggest that, in the absence of an external pressure gradient, the hydrodynamics of ciliated tentacles may not undergo any sudden fundamental changes in the transition from Re1 (at which inertia is strictly negligible) to moderate Re values (at which inertial effects could begin to be important). This is consistent with results from some other flows characterized by small length and velocity scales (e.g. flow around a sphere), in which experiments show that low-Re approximations may be useful predictors even for moderate Re values (Batchelor, 1983). If this conclusion is supported by further work, it may imply that the potential hydrodynamic operating range of ciliated tentacle geometries extends well into moderate Re values. That would in turn suggest that morphological and physiological constraints may place important limits on the sizes and speeds of these arrays, rather than abrupt changes in the hydrodynamic regime, at least in the absence of strong external pressure gradients.

Fig. 4 shows non-dimensional flow rates for a standard geometry (variables are given in the legend) with variation in the cilia gap (g) and with pressure drop (Δp) for three Re values that span much of the biologically relevant range (there is negligible variation for Re<0.1, so the curves in Fig. 4A apply to lower Re values as well). Two interesting points emerge from these results. First, in the absence of a pressure drop, the flow rates are concave downwards, meaning that there are diminishing returns for increasing cilia gap beyond a flow-maximizing value. The decrease of flow for large cilia gaps is especially pronounced at high Re values, but the low-Re curves also bend downwards at large gaps (not shown on these plots). The scaled flows (normalized by the ‘shape factor’ Q_e) are presented in Fig. 5. This figure shows that, at low Re, the trapezoid profile is a slight overestimate of the simulated flows for small cilia gaps and a more substantial overestimate for large cilia gaps. At the higher Re values, the trapezoid profile is an underestimate for small cilia gaps, but a substantial overestimate for large cilia gaps.

The second interesting feature of Fig. 4 is that, in each plot, equal increments in Δp_n result in equal decrements in flow, i.e. the flow deficit is a linear function of pressure drop throughout this range of Re values. This useful result is reflected in the plots of flow deficits for the three Re values in Figs 6–8. Note that while the flow deficits vary differently at the different Re values, in each case they collapse to a single curve characteristic of that Re when normalized as the scaled flow deficit, (Figs 6B, 7B, 8B). That is, is a function of Re and g_n, but is invariant in Δp_n in the simulated flows.

The results presented in Figs 3–8 (and many other simulations) can be reconstructed exactly using a simple scaling and interpolation scheme, together with the data in Tables 2 and 3. These tables are compilations of Q₀ and values from our simulations (these data are also available as an electronic file from the authors). Using equations 1–6, this scheme allows the reader to estimate the flow through any tentacle geometry within the scope of our simulations, i.e. within the variable ranges Re ⩽10 and 0.5 ⩽g_n ⩽10. We used two-dimensional bicubic interpolation (available in a number of standard numerical packages such as Matlab and Maple), but a simpler linear interpolation will also work well.

To recover the unscaled, physical flows for a given set of variables (g, c, s, L, D, Δp, ρ, μ, where ρ is fluid density) from these tables, one interpolates values of Q₀ and ⁠, using g_n, c_n, Re_nom=sL_nomρ/μ and Δp_n, where L_nom=L−g+g_nom, g_nom=12.208c, and other non-dimensional variables are as shown in Table 1. The nominal values of length and Reynolds number (L_nom and Re_nom) are to facilitate interpolation by preserving the rectangular layout of simulation results. Geometries with Re_nom<0.1 are treated as if Re_nom=0.1, because there is no Re-dependence below this Reynolds number. True values of Re can be calculated from the flows after interpolation and rescaling.

The blue mussel M. edulis is one of the best-studied suspension feeders that actively pump using ciliated tentacular arrays. We used variables from the literature to test whether the geometry of the mussel’s feeding structures correspond to flow-maximizing geometries predicted by our computational model. There is disagreement between authors about the precise values of many variables (Jones and Allen, 1986; Jørgensen et al. 1988; Silvester, 1988; Jones et al. 1992; Nielsen et al. 1993). The variables we used were c=15 μm, D=40 μm, g=10 μm and s=2 mm s⁻¹. For sea water at 10 °C, ρ=1023 kg m⁻³ and μ=1.4×10⁻³ Ns m⁻² (Roberson and Crowe, 1980; Denny, 1988).

Fig. 9 shows characteristic performance curves, i.e. plots of non-dimensional flow rate against pressure drop (or ‘head loss’ in engineering terminology), for M. edulis for several values of cilia gap. These plots show that there are strong differences in how arrays with different gaps respond to adverse pressure gradients. Arrays with large gaps (e.g. g=40 μm) have high flow rates when pressure gradients are absent or extremely weak. However, this geometry is not robust: flow rates fall dramatically in response to relatively moderate adverse pressure drops. In contrast, arrays with small gaps (e.g. g=5 μm) have lower flow rates when pressure drops are absent, but their flow rates decrease much less in response to an adverse pressure gradient. The linear relationship between back pressure and flow rate in M. edulis predicted by our model is similar to that observed experimentally by Jørgensen (1989). The cilia gap in a tentacular array functions much like a gear ratio in mechanical systems, with smaller cilia gaps corresponding to lower gear ratios. In many mechanical systems, each gear ratio has a range over which it maximizes power output. A similar effect is reflected in Fig. 9: when Δp is near zero, the wide gap results in the highest flow rates. As Δp increases, the next highest gap (g=20 μm) becomes the most effective. This is followed with increasing pressure by ranges in which gaps of 10 μm and 5 μm have the highest flow rates. From the literature cited above, mussels appear to operate at pressure drops across their tentacular arrays of approximately 3–5 mmH₂O (29–49 Pa) and at a cilia gap in the region of 10 μm. These ranges are consistent with the model’s predictions of the gaps that maximize flow rates under the prevailing hydrodynamic conditions.

The analogy with gear ratios is developed more fully in Fig. 10. Fig. 10A shows the flow-maximizing gap (g_max) as a function of pressure drop (Δp). To give an indication of how sensitive performance is to deviations from g_max, this plot also shows the ranges of g over which flow rate is at least 95 % and 90 % of maximum. These simulation results suggest that performance is good even for deviations from g_max as large as 30 %. Furthermore, the curve suggests that conservative strategies are possible: a gap of 10 μm is within 5 % of the maximum flow rate for the entire range 1.5 ⩽Δp ⩽7, which probably spans the biologically relevant range. In contrast, a gap of 20 μm, which is very close to g_max at Δp=1.5, is within 5 % of the maximum over a much smaller range of approximately 0.5 ⩽Δp ⩽2.2. Many suspension feeders experience rapid variations in pressure, for example those due to oscillatory external flows (waves and tides) and wind-driven currents or to intermittent interference from epizoites. Our results suggest that an animal that experiences a range of pressure drops but has a fixed cilia gap would perform best overall with a gap that is close to the g_max for the higher range of Δp values it is likely to encounter.

Fig. 10B shows the non-dimensional flow rates, Q_max, at the flow-maximizing cilia gap, g_max. Assuming that an animal is able to adjust its tentacle spacing to the best possible cilia gap, the cost (in terms of decreased flow rate) of increasing Δp is most pronounced at low pressure drops. For Δp ⩾2 mmH₂O (19.6 Pa), flow rate decreases exponentially, while for Δp<2 mmH₂O, the decrease is faster than exponential.

Because the pressure drop across the tentacular pump depends on the kind of piping system to which it is attached, while the flow rate through the piping system varies with the pressure drop across the pump, the real design trade-offs facing the animal are frequently not as simple as maximizing flow at fixed Δp. Equations such as 11 and 12 can be used to determine graphically the operating conditions of a piping system consisting of a tentacular pump and a set of associated conduit elements. This technique is discussed in many engineering texts on fluid mechanics (see Roberson and Crowe, 1980, for a good discussion). It stems from the requirement for consistency in both variables (Δp and Q) between the pump and the rest of the piping system. In Fig. 10B, flow rates of the form of equation 12 are shown for several values of k₂. Assuming, as before, that the animal adjusts its tentacle spacing to g_max for any Δp, the pump always operates somewhere along the solid curve in Fig. 10B. Suppose, for example, that the flow through the piping system as a function of Δp is characterized by equation 12 with k₂=0.2. Then, the piping system always operates along the corresponding dotted line. There is only one flow rate and associated pressure drop that satisfy both curves simultaneously, marked by their point of intersection. This point represents the values of Q and Δp at which the combined pump and piping system operate. The effect of varying the properties of the piping system can easily be determined using this method. For instance, halving the conductivity to k₂=0.1 (i.e. doubling the resistance to flow through the piping system) gives a new point of intersection. The new operating condition shows that doubling the resistance of the pipe decreases the flow rate by less than 25 %, far less than might be expected if the characteristics of the ciliary pump were not considered. Interestingly, ciliary gaps in M. edulis are observed to be smaller in smaller individuals (Jones et al. 1992), possibly reflecting adjustments for increased fluid resistance through narrower apertures, which is consistent with predictions based on Fig. 10B.

In both suspension feeding and respiration, many mechanisms and constraints may be involved in the functioning of ciliated tentacle arrays (Jones et al. 1990; Beninger et al. 1992). Nonetheless, given the large investment in terms of energy, biomass and body volume made by many of these animals in producing water currents, it is reasonable in many cases to hypothesize that their pumping structures have evolved towards geometries that maximize flow per unit cost. In this study, we used a hydrodynamic model of flows through arrays of ciliated tentacles to identify these favorable geometries. We explored the consequences of varying the operating characteristics of biological piping systems and compared them in a case study based on M. edulis with experimentally observed geometries and operating conditions. We discussed these results primarily in terms of flow rates per unit length of tentacle array. However, if the energy consumption per unit length of cilia band is assumed to be approximately constant over a range of operating conditions, these results may translate directly into relative measures of energetic efficiency (Clemmesen and Jørgensen, 1987). We also presented an interpolation and scaling scheme, together with tabulated flow rates from our simulations, that provide the reader with the means to estimate the flows for a wide range of geometrical and hydrodynamic variables relevant to many animals that we did not consider explicitly.

What are the limits to the volume of water that can be processed by an array of ciliated tentacles? Our simulations suggest that, when pressure drops are absent or small, a crude estimate such as the trapezoidal profile (equations 3, 4) is approximately correct for non-dimensional cilia gaps below g_n=3 and is an overestimate by approximately a factor of two for large gaps (Fig. 7). This type of estimate may give a preliminary indication of how physical constraints on cilia length and speed translate into limitations on water-processing rates, particularly if it is interpreted as a theoretical upper limit. However, this conclusion is complicated considerably when significant adverse pressure gradients are present. The degree to which flow is reduced by pressure is strongly affected by the gap between cilia on adjacent tentacles. The large gaps that perform exceptionally well with no pressure drop are also extremely susceptible to adverse pressures. Because pressure drops are potentially present in most animals that use ciliated tentacle arrays, simple analytical estimates such as equations 3, 4 are of limited usefulness in understanding the functional morphology of this pumping mechanism.

Both the limitations of our computational model and differences in array morphology among taxa imposed some simplifications on our model’s representation of geometrical detail, in particular with respect to the cilia. In most natural examples of tentacular arrays, large numbers of cilia are present in dense bands on the lateral tentacle surfaces (Silvester, 1988). Flows around individual cilia are complicated by a recovery phase in the stroke, by possible lateral plectic movements during the power phase and by three-dimensional ‘slippage’ around cilia (Sleigh and Aiello, 1972; Sleigh, 1989; Emlet, 1991; Hansen and Tiselius, 1992). The hydrodynamic effects of these complications can be approximated for isolated cilia; however, it is not yet possible to do so for dense ciliary arrays (Silvester, 1988; Nielsen et al. 1993). Because our primary interest is in the simpler flows induced in the far-field above the ciliary tips, and because this flow is observed to be a substantial fraction of the cilia power-stroke velocity in the vicinity of the cilia tips (Silvester, 1988; Jones et al. 1992), we modeled bands of large numbers of cilia with complex strokes and slippage as bands of relatively fewer cilia with no slippage and with simplified strokes. We note that each cilium in our two-dimensional simulation would translate into an impermeable flat plate in three dimensions. This corresponds to assuming that, in the real tentacular arrays, cilia are packed in dense rows axially along the tentacle and move approximately in phase along a row so that no fluid ‘leaks’ between them. These simplifications need to be evaluated experimentally – a more detailed theoretical or experimental analysis may suggest circumstances under which ciliary arrays that are less densely packed or out of phase have an important degree of permeability or an effective cilium tip speed lower than the nominal tip speed. However, we believe that our simulations provide a good first approximation to the aspects of the flows with which we are most concerned, i.e. the balance between viscous traction forces imparted by the cilia, resistance to flow through the tentacles and pressure forces that oppose the flow.

Our simulation results emphasize that even small pressure drops across tentacle arrays can have important functional consequences. In the case of M. edulis, for example, flow rates per unit length of tentacle array are predicted to be substantially higher at Δp=0 mmH₂O than at Δp;:,1 mmH₂O (9.8 Pa). Is it possible for animals to operate very near to zero pressure drop to take advantage of this difference?

One of the most important functions of piping system geometries in many suspension-feeding and respiratory flows is to prevent recirculation of depleted water. In many cases, such as exhalant siphons in bivalves and montecules (or ‘chimneys’) in highly integrated bryozoans, depleted water is eliminated by imparting it with sufficient kinetic energy to carry it beyond the incurrent. Because cilia operate through viscous traction forces, they are ineffective at directly producing flows with high kinetic energy. Instead, highly energetic flows must be produced in animals that use ciliary pumping by some sort of flow constriction in the piping system. Such a constriction exchanges potential energy (in the form of a pressure drop) for kinetic energy. Thus, animals using high-kinetic-energy excurrents to eliminate depleted water cannot operate with zero or near-zero pressure drops.

This constraint invites the speculation that animals moving water with ciliated tentacle arrays may commonly specialize into one of two strategies: (i) relatively small cilia gaps coupled with internal piping systems and energetic excurrent flows; and (ii) relatively large cilia gaps not coupled to internal piping systems and without energetic excurrent flows. According to this logic, animals should avoid the intermediate pressure drops that are large enough to reduce flow rates but too small to eliminate depleted water kinetically. Bryozoans may be an interesting group in which to test this prediction (Winston, 1977, 1978, 1979; McKinney and Jackson, 1989; Grünbaum, 1995, 1997), because closely related species may be either highly integrated (effectively with internal piping systems) or poorly integrated (effectively without internal piping systems). Interesting exceptions to the prediction might include species in which overall geometry or external flows prevent recirculation without energetic excurrent flows, such as the serpulid Spirobranchus giganteus (Strathmann et al. 1984), species that switch between suspension feeding modes depending on hydrodynamic conditions (Miller et al. 1992) or species in which the functions of ciliary pumping have changed over evolutionary time, such as hydrothermal vent fauna (Fiala-Médioni and Métivier, 1986).

The same reasoning might apply more broadly to suspension feeders in general: small suspension feeders, whose flow resistances are high and whose total flow rates are too low to produce energetic flows through constrictions, should be more successful with external tentacle arrays. In contrast, large suspension feeders, whose total flow rates are high enough to produce energetic flows through constrictions, should be more successful with internal tentacle arrays. Many suspension feeders with internal piping systems that are large as adults must pass through a juvenile stage in which they have internal piping systems but are small (Beninger et al. 1994). The allometric characteristics of ciliary pumping in our simulations suggest that small juvenile stages face biomechanical constraints on pumping that may significantly reduce growth and survivorship.

Many suspension feeders with internal piping systems as adults and juveniles utilize external filtration organs as planktonic larvae. For example, inarticulate brachiopods shift from external to internal flows in the course of their development (Barnes, 1986). Our results suggest that this might be a strategy to circumvent pumping limitations on very small juvenile stages. In keeping with this analysis, the size at which veligers (with external ciliary flow) of many bivalve species metamorphose into juveniles with internal piping systems is fairly consistently in the range 200–300 mm (Sastrys, 1979). Our results suggest that one of the factors affecting the timing of metamorphosis may be the hydromechanical disadvantages of having internal ciliary pumps smaller than this threshold size.

The simulations in this paper were carried out using the IBIS software package (Fogelson and Eyre, 1997). IBIS allows biologists and others without expertise in mathematics or computations to make use of a powerful computational tool for biological fluid mechanics known as the immersed boundary method. This method has been used successfully on a wide variety of problems including blood flow in a beating heart (Peskin, 1977; Peskin and McQueen, 1993), platelet aggregation during hemostasis and thrombosis (Fogelson, 1984), swimming of aquatic animals (Fauci and Peskin, 1988; Fauci, 1993), biofilm development (Dillon et al. 1996) and the dynamics of the inner ear (Beyer, 1992). What these problems have in common are interactions between a viscous fluid and one or more deformable, perhaps active, objects immersed in the fluid. The geometry of these objects can be quite complex, and IBIS was designed to allow the user to build up complex geometries and motions from a library of simpler ones. Information about the IBIS package, including examples of its use, is available from the authors or at web-site http://www.math.utah.edu/IBIS. Here, we give a brief description of the immersed boundary method as applied to the tentacle problem; much more detailed discussions about the method can be found in the references cited above.

For the ciliated tentacle simulations, we solve the equations governing the coupled motion of a viscous incompressible fluid and ciliated tentacles immersed in that fluid. The equations that describe the fluid motion are the fully time-dependent non-linear incompressible Navier–Stokes equations. We do not assume that fluid inertia is negligible, and so we determine flows with non-zero Reynolds number. Because the flow is driven in part by beating cilia, the flow is unsteady (i.e. not constant in time) and so we solve the fully time-dependent equations. We impose ‘no slip’ boundary conditions at points where the fluid abuts a solid object; this means that, at such points, the tangential and normal components of the velocity of the fluid are the same as those of the object.

Here, u(x,t) is the fluid velocity vector at location x and time t, p(x,t) is the fluid pressure, and ρ and μ are the constant fluid density and viscosity, respectively. The force density f(x,t) which drives the fluid motion is generated in large part by forces designed to maintain the shape and location of the tentacles and to drive the beating of the cilia. The Navier–Stokes equations (equations A1 and A2) hold throughout the domain in which our computations are carried out. The presence of the tentacles and the beating cilia is ‘felt’ by the fluid only through the force density f. This is described further below.

In this paper, we consider situations in which the tentacle geometry and the flow variables depend only on two spatial variables, those in a plane that cuts perpendicular to the main tentacle axes (see Fig. 11A). Nothing in our methodology depends essentially on the simulations being two-dimensional; however, the computational requirements would increase considerably for three-dimensional calculations.

The force density F for the power stroke of a beating cilium is illustrated in Fig. 11B. In the model, a target cilium sweeps through a prescribed angle at a prescribed angular velocity during the power stroke. Each point of the actual cilium is ‘tethered’ by an elastic spring to the corresponding point of the target cilium. The spring generates a force at the cilium point directed towards the corresponding point on the target cilium. How closely the actual cilium tracks the prescribed motion of the target cilium depends on the stiffness of the springs and on the properties of the local flow. The latter matters because the action of the spring forces on the cilium is mediated through the flow, as described above. At the end of the power stroke, the locations of the target cilium and the actual cilium are re-initialized to their starting locations in preparation for the next power stroke. During this recovery motion, the cilium applies no force to the fluid.

The contributions to the force density F at points on the tentacle surface are depicted in Fig. 11C. Again, there is a target tentacle surface and each point of the actual tentacle surface is tethered by an elastic spring to the corresponding point of the target tentacle surface. Since we specify that the target tentacle points are stationary, these springs serve to anchor the tentacle in an approximately fixed location. A second contribution to the tentacle force density F comes from elastic springs within the tentacle surface ring itself that respond to motions that perturb the size or shape of the ring away from its prescribed equilibrium configuration.

For actual simulations, the model system given by equations A1–A4 is approximated by a discrete system of algebraic equations. Approximate solutions are obtained at discrete instances of time t_n=nΔt separated by a time step Δt. For the Navier–Stokes equations, a uniform computation grid is placed over the domain and a finite-difference scheme due to Chorin (1968) is used to find approximate values of the velocity u at points of the grid. Each function X(α,t_n) is also represented by a discrete set of points each of which moves according to a discrete analogue of equation A3. The integrals that appear in equations A3 and A4 are approximated by sums, and the Dirac delta function is replaced by an approximate δ function (see Peskin, 1977) whose influence is spread over a small but non-vanishing portion of the domain. For the simulations, the numerical variables (e.g. time step and grid spacing) were reduced until the results became insensitive to further reduction. We re-emphasize that users of the IBIS software need not be concerned with these technical details of the immersed boundary method. IBIS provides simple user interfaces for setting up the problem and for displaying the simulation results graphically.

 
• c

  length of the cilia

 
• c_n

  non-dimensional length of the cilia

 
• D

  tentacle diameter

 
• g

  gap between tips of cilia on adjacent tentacles

 
• g_max

  flow-maximizing cilia gap

 
• g_n

  non-dimensional cilia gap

 
• g_nom

  nominal cilia gap

 
• k₁, k₂

  conductivity coefficients

 
• L

  tentacle spacing

 
• L_nom

  nominal tentacle spacing (L−g+g_nom)

 
• Q

  flow rate

 
• Q_d

  non-dimensional flow deficit scaled flow deficit

 
• Q_e

  non-dimensional trapezoidal ‘shape factor’

 
• Q_max

  non-dimensional flow rate at flow-maximizing cilia gap, g_max

 
• Q_n

  non-dimensional flow rate

 
• Q_s

  scaled flow rate

 
• Q₀

  non-dimensional flow with no pressure drop (Δp_n=0) tentacle diameter ratio,

 
• r

  tentacle diameter ratio,D/c

 
• Re

  Reynolds number

 
• Re_nom

  nominal Reynolds number

 
• S

  tip speed of the cilia

 
• V

  flow velocity

 
• Δp

  pressure drop

 
• Δp_n

  non-dimensional pressure drop

 
• μ

  dynamic viscosity

 
• ν

  kinematic viscosity

 
• ρ

  fluid density

The authors thank R. Emlet, C. Jordan, M. LaBarbera, J. Levinton, R. Self, R. Strathmann and J. West for helpful conversations and comments on the manuscript. D.G. was supported by the Special Year in Mathematical Biology at the University of Utah, NSF Grant 9503478, to H. Othmer, F. Adler, M. Lewis, A. Fogelson and J. Keener, and by Office of Naval Research Grant N00014-92-J-1527 to H. Caswell. The work of A.F. and D.E. was supported in part by NSF Grant DMS-9307643. The authors are indebted to the staff and scientists of Friday Harbor Laboratories and its director, A. O. Dennis Willows.