## ABSTRACT

Parallel tentacular structures with lateral cilia that produce suspension-feeding and respiratory flows occur repeatedly in many diverse taxonomic groups. We use a computational hydrodynamic model of flow through ciliated tentacles to simulate flow rates through ciliated tentacle arrays. We examine the functional relationship of one performance measure, flow rate per unit length of array, to geometrical variables, such as cilia length, cilia tip speed and the gap between adjacent tentacles, and to hydrodynamic operating conditions, such as adverse pressure drops across the array. We present a scaling and interpolation scheme to estimate flow rates for a wide range of geometries that span many taxa. Our estimates of flow rate can be coupled with the hydrodynamic characteristics of biological piping systems to understand design trade-offs between components of these systems. As a case study, we apply the model to the blue mussel *Mytilus edulis* by investigating the effect on performance of changes in the gap between neighboring tentacles. Our model suggests that the observed gaps between tentacles in *M. edulis* reflect flow-maximizing geometries. Even relatively weak adverse pressure drops have strong effects on flow-maximizing geometries and flow rates. One consequence is that an intermediate range of pressure drops may be unfavorable, suggesting that animals may specialize into high-pressure and low-pressure piping systems associated with differences in organism size and with their strategy for eliminating depleted water.

## Introduction

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.

## Materials and methods

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.

### Scaling of flow rates

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^{2} s^{−1}. 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.

*p*, against which the cilia are working. The importance of pressure forces depends in large part on their magnitude compared with other fluid forces, primarily viscous forces. A scaling that expresses the relative magnitude of pressure forces is the non-dimensional pressure, Δ

*p*

_{n}: where μ is the dynamic viscosity of the water and

*c*is the length of the cilia.

*Q*

_{0}to be the non-dimensional flow with no pressure drop (Δ

*p*

_{n}=0). A useful baseline estimate of this flow,

*Q*

_{e}, is obtained by assuming that the velocity profile between the tentacles has a trapezoidal shape, as in Fig. 1. The equation for the estimated trapezoidal flow

*Q*

_{e}is: where

*g*is the gap between the tips of the cilia on adjacent tentacles.

*Q*

_{e}is in a sense a ‘shape factor’ that discounts for the ‘porosity’ of the array (Jørgensen

*et al*. 1988). A scaled flow rate

*Q*

_{s}, where: should be approximately equal to unity if the true profile is similar to the trapezoidal profile shown schematically in Fig. 1.

*p*

_{n}>0) can be expressed as a flow deficit,

*Q*

_{d}, where:

*Q*

_{d}reflects the ‘leakage’ of water driven through the gap between the cilia by the drop in pressure (Fig. 1). As suggested by fluid dynamic theory and confirmed by the simulations of Nielsen

*et al*. (1993), at the

*Re*values characteristic of tentacle arrays, the leakage has an approximately parabolic profile (as in Poiseuille flow). The magnitude of this flow is roughly proportional to the pressure drop, to the square of the gap between cilia and to the inverse of the viscosity. These trends are factored out by considering a scaled flow deficit,

*Q*

_{d}

^{*}, where: Note that

*g*

_{n}is the non-dimensional cilia gap (see Table 1).

## Results

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

### Allometry and geometrical variation

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 *Re*1 (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.

### Interpolation scheme

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*_{0} 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*_{0} and , using *g*_{n}, *c*_{n}, *Re*_{nom}=*sL*_{nom}ρ/μ and Δ*p*_{n}, where *L*_{nom}=*L*−*g*+*g*_{nom}, *g*_{nom}=12.208*c*, 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.

### Ciliary pumping in Mytilus edulis

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^{−1}. For sea water at 10 °C, ρ=1023 kg m^{−3} and μ=1.4×10^{−3} Ns m^{−2} (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_{2}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_{2}O (19.6 Pa), flow rate decreases exponentially, while for Δ*p*<2 mmH_{2}O, the decrease is faster than exponential.

*k*

_{1}and

*k*

_{2}values indicate higher resistance to flow through the piping system, as would occur, for instance, with narrower or longer pipe segments. Jørgensen

*et al*. (1988) suggest that in

*M. edulis*there are two primary components of resistance. One is the tentacular array itself, which was considered separately from the ciliated region in their model but is part of the flow calculation in our model (and therefore is already accounted for in our flow rate estimates). The other is the kinetic energy carried away by the water as it exits the exhalant siphon, which was estimated by Jørgensen

*et al*. (1988) to have a much larger magnitude than losses at the inlet and outlet and in other conduit elements. For a kinetic energy loss term, the square root characteristic curve in equation 12 is most appropriate.

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*_{2}. 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*_{2}=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*_{2}=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.

## Discussion

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.

### Model simplifications

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.

### Implications for ciliary pumping

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_{2}O than at Δ*p*;:,1 mmH_{2}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.

### Appendix

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.

**X**(α,

*t*). Here, as the variable α varies between 0 and 1,

**X**(α,

*t*) sweeps through the points on the surface of the tentacle or along the length of the cilium at time

*t*. Each value of α refers to a particular material point on the tentacle or cilium, and the function

**X**(α,

*t*) with α held fixed describes the trajectory of this point in time. (This is a Lagrangian description of the tentacles and cilia.) Because the tentacle surface and the plane of our calculations intersect in a closed curve, the condition

**X**(0,

*t*)=

**X**(1,

*t*) holds for each tentacle and for all values of

*t*. Each point on a tentacle surface or cilium is in contact with the surrounding fluid and so its velocity must be consistent with the no-slip boundary condition. This gives us the equation of motion for the point as: where δ represents the Dirac delta function, and the integral of the velocity

**u**times δ[

**x**−

**X**(α,

*t*)] picks out the value of

**u**at the point

**X**(α,

*t*).

**F**[

**X**(α,

*t*),

*t*] along each tentacle surface and cilium. These forces are transmitted to the surrounding fluid by integrating

**F**[

**X**(α,

*t*),

*t*] times δ[

**x**−

**X**(α,

*t*)] over the relevant surfaces. Since each tentacle and cilium contributes to the total force density driving the fluid, the fluid force density is given by the equation: Here, the sum is over the

*N*tentacle surfaces and cilia in the calculation. This is the force density

**f**that appears in the Navier–Stokes equations (equation A1).

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.

## List of abbreviations

*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*_{1},*k*_{2}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*_{0}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

## ACKNOWLEDGEMENTS

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.

## References

*An Introduction to Fluid Dynamics*

*Pecten maximus*(Bivalvia: Pectinidae)

*Mar. Biol*

*Placopecten magellanicus*(Mollusca: Bivalvia) as revealed using video endoscopy

*Mar. Biol*

*J. comput. Phys*

*Mytilus edulis*L. as a function of boundary-layer flow

*Limnol. Oceangr*

*Math. Comput*

*Mar. Biol*

*Biology and the Mechanics of the Wave-Swept Environment*

*J. comput. Phys*

*Biol. Bull. mar. biol. Lab., Woods Hole*

*Am. Zool*

*Contemp. Math*

*J. comput. Phys*

*Calyptogena magnifica*, with a discussion of its nutrition

*Mar. Biol.*

*J. comput. Phys*

*J. theor. Biol*

*Membranipora membranacea*

*Limnol. Oceangr*

*J. Plank. Res*

*Biol. Bull. mar. biol. Lab., Woods Hole*

*Biol. Bull. mar. biol. Lab., Woods Hole*

*Mytilus edulis*L. and

*Cerastoderma edule*(L

*J. exp. mar. Biol. Ecol*

*Mytilus edulis*L

*J. exp. mar. Biol. Ecol.*

*Mytilus edulis*L

*J. exp. mar. Biol. Ecol.*

*Comp. Biochem. Physiol.*

*Mar. Ecol. Prog. Ser.*

*Science*

*Mytilus edulis, Cerastoderma edule, Mya arenaria*and

*Amphiura filiformis*

*J. Sea Res*

*Bryozoan Evolution*

*Mar. Biol*

*J. mar. Res.*

*Mytilus edulis*: video recordings and numerical modelling

*Mar. Biol*

*J. comput. Phys*

*Contemp. Math*

*Engineering Fluid Mechanics*

*Reproduction of Marine Invertebrates*

*Molluscs: Pelecypods and Lesser Classes*

*Dreissena polymorpha, Corbicula fluminea*and

*Carunculina texasensis*

*Biol. Bull. mar. biol. Lab., Woods Hole*

*Mytilus*gills

*J. exp. mar. Biol. Ecol*

*Comp. Biochem. Physiol*

*Acta protozool*

*Mar. Biol.*

*Spirobranchus giganteus*(Pallas) breaks a rule for suspension feeders

*J. exp. mar. Biol. Ecol*

*The Biology of Bryozoans*

*Bull. mar. Sci.*

*Advances in Bryozoology*