Biomimetic cilia arrays generate simultaneous pumping and mixing regimes

Pumping Above the Cilia Tips.

Above the tips of the beating cilia, the time-averaged tracer motion is uniform, directed transport (Fig. 3A and Movies S1 and S2). This directed flow is analogous to the flow produced by nodal cilia (13, 35). At low magnification we have observed that this net transport is roughly uniform across the 1-mm width of the fluidic cell (Movie S3). Consistent with in vivo observations of the nodal flow (13, 35) and with computational work (36), the direction of the flow is determined by the direction of the motion of each cilium when its tip is farthest from the specimen floor (Movie S1). We also observe a slow recirculation of fluid along the chamber’s upper boundary due to the enclosed nature of the flow cell. This recirculatory flow has also been reported in the embryonic node (7, 13).

Our system offers a large degree of control over the directed fluid transport. The relative positioning of the actuating magnet and the cilia array determines the cilia tilt direction and thus the fluid transport direction (Fig. 2 A and B). By repositioning the magnet we can arbitrarily choose the pumping direction. Further, for a given magnet position the transport direction can be switched 180° by reversing the motor drive. We can also control the flow velocity by adjusting the cilia beat frequency and find a linear relationship between the two (Fig. S1). This flexibility in directionality and velocity of pumping would be a distinct advantage of biomimetic cilia devices in a microfluidics setting.

Cilia-Driven Pumping Is Poiseuille–Couette Flow.

A coarse-grained description of cilia-generated fluid flows is an essential component of large-scale computational modeling of ciliated systems. To construct a velocity flow profile we averaged tracer velocities in the direction of fluid transport at a number of focal planes above and below the biomimetic cilia tips (Fig. 4A). Considering only the portion of the velocity profile above the cilia tips (from z = 25 μm → 225 μm), we find that the directed flow can be explained by Poiseuille–Couette (PC) flow, a linear superposition of shear- and pressure-driven flow in a parallel plate geometry. PC flow is described by a velocity profile [2]where z^′ = z - 25 μm and the bottom plate at z^′ = 0 (at the cilia tips) slides with velocity u₀ relative to the stationary upper plate at z^′ = h, exerting a shear stress on a fluid of viscosity η under a pressure gradient ∇p. In Fig. 4A the solid curve was obtained as a least-squares fit of the biomimetic cilia-driven velocity profile above the cilia tips to Eq. 2. Thus, the collective motion of thousands of cilia generates a velocity profile that can be coarse-grained as an effective shear stress produced at the cilia tips (Couette flow) and an opposing induced pressure gradient generated by the interaction of the fluid with the closed flow cell boundaries (Poiseuille flow).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

The velocity flow profile and effective relative diffusivity of biomimetic cilia-driven flow. The area below the cilia tips is denoted by the shaded regions. (A) Red points are the tracers’ average velocity, in the transport direction, at various heights above the floor. The solid blue curve is a least-squares fit to the Poiseuille–Couette flow profile of Eq. 2 with fit parameters u₀ = 8.7 μm/s and ∇p = 1.05 Pa/m. The maximum average velocity occurs just above the cilia tips, demonstrating that the collectively driven transport at the tips acts analogously to the velocity u₀ of the sliding plane in PC flow. Because the chamber is enclosed, there is a recirculation of fluid near the upper boundary. (Inset) We fit the PC flow profile to mouse nodal flow velocity data at three heights taken from Okada et al. (13), with parameters u₀ = 4 μm/s and ∇p = 79 Pa/m. The discrepancy in ∇p between the node and our flow cell is likely due to the overall size differences of the two systems (see Appendix). (B) The biomimetic cilia motion induces mixing by producing a maximum relative diffusivity of 42 μm²/s, an enhancement of 25 in relation to the expected relative diffusivity of (vertical magenta line) from the Stokes–Einstein relation. Above the cilia tips, the relative diffusivity rapidly decays back to the expected value. Horizontal error bars are on the order of the dot size. (Inset) Tracer MSDs and least-squares linear fits to the log of Eq. 5, from above (orange) and below (green) the cilia tips. The solid magenta line is the theoretical MSD with the expected relative diffusivity D₀ and slope of 1.

Previous measurements of fluid flow within the node of a mouse embryo documented the variation of flow away from the plane of the nodal cilia tips (13). In Fig. 4A Inset, we have taken this data from Okada et al. (13) and fit it to the PC flow profile assuming the approximate geometry stated in their work. We can see that Poiseuille–Couette flow is consistent with the flow in our engineered fluid cell and that of the mouse node, suggesting that the three-dimensional nature of the flows in these systems is that of a “driven cavity” (see Appendix) (37). We note that, more recently, Okada et al. performed a numerical simulation of flow in an open pit that is driven by a shear plane, as in Couette flow, in order to confirm the persistence of the recirculatory flow despite their removal of the membranous covering of the nodal cavity (38).

We also note that neither fit to Eq. 2 results in a net zero flux as expected for a two-dimensional, closed system. As pointed out by Liron, this behavior is due to a portion of the return flow that, in three dimensions, is recirculated around the periphery of the system (39). This peripheral recirculation is also evident in the embryonic node in supplemental videos from Okada et al. (13).

Volume Flow Rate.

Another goal in the coarse graining of cilia-driven fluid flows is to be able to predict flow rates and fluid velocities of a system of cilia from knowledge of the beat shape of a single cilium in that system. A result from resistive force theory predicts the volume flow rate for a single cilium based on the parameters of the tilted conical beat, [3]where ω is the angular frequency and g = 4/3[1 + 2 ln(2q/a)]^-1 is a geometric factor in which a is the diameter and q is a length parameter of order (36). In contrast to our experiment in an enclosed flow cell, this result is derived in a semiinfinite space, and so we take this value to be the maximum possible flow rate a single cilium can produce in a system of negligible resistance, or . Thus, from the average beat parameters of the cilia driving our velocity profile (see Materials and Methods) we have Q_max = 1.3 ± 0.3 pL/s per cilium.

Experimentally, we can measure the volume flow rate produced by the entire system by noting that for flow in a cell whose width w is larger than its height h (as is the case in both the mouse node and our flow cell) the total volume flow rate can be approximated as [4]where u(z^′) is the velocity flow profile. Note that we use |u(z^′)| because both the primary flow and the recirculation are attributed to the cilia beat. Performing the integration on the velocity profile of Fig. 4A we get an experimental value of Q^total = 460 ± 50 pL/s.

In order to compare this experimental measurement with the theoretical maximum flow rate per cilium of Q_max = 1.3 pL/s, we would need to know how the flow driven by each cilium in our system contributes to the fluid flow of the entire system. In general, however, this requires knowledge of the interactions between neighboring cilia for which there is no current theory.

From an analogy with classical pumps we can provide limits for the comparison with experiment. In the case of noninteracting cilia, which experimentally would only be realized in a sparse array of cilia, the flow in the system is simply a superposition of each cilium, and thus the N cilia in the array should act like classical pumps that are configured in parallel such that the maximum possible total flow rate would be , or . This upper bound is well above the experimental value of 460 pL/s; however, the actual operating point of the model would depend on losses in the system.

In contrast, for strongly interacting cilia we assume that only cilia arrayed across the width of the fluid cell, perpendicular to the transport direction, would behave as classical pumps in parallel. The rest of the cilia, arrayed down the channel, would then act as pumps in series, which would not contribute to . Thus, for n_p cilia arrayed across the width of the cell we have . In our experiment we have n_p ≈ 55, giving a lower bound of . In this case, the discrepancy with our experimental measurement implies that cilia arrayed down the length of the channel are still able to make an effective contribution to .

Thus, our experimental value falls within upper and lower bounds considered for non- and strongly interacting cilia. We are optimistic that future experiments performed by varying the number of cilia in an open-ended channel may shed further light on the nature of cilia-cilia interactions.

Enhanced Diffusion Below the Cilia Tips.

Below the cilia tips, the average tracer motion is strikingly different from the directed flow above (Fig. 3B and Movie S2). Here, the fluid motion is dominated by local vortices associated with each beating cilium and tracer trajectories that briefly orbit within then jump to nearby vortices. While similar vortical motion in close proximity to a single cilium has been qualitatively described in experiments (13, 29, 40) and simulations (10, 36), no quantitative analysis of transport phenomena has been undertaken in any ciliated system.

In contrast to the directed transport regime, the flow below the cilia tips shows little overall (time-averaged) directionality, with the standard deviations of the average tracer velocities being at least as large as the mean velocities, except over the 5-μm transition just beneath the cilia tips (Fig. 4A). Tracers exhibit large fluctuations in speed, with the magnitude of individual tracer velocities being largest just below the tips where tracers reach average speeds as high as 30–40 μm/s, as measured over a 0.25-s window.

The nature of this flow regime is suggestive that cilia motion could mix fluids at the microscale. At low Reynolds number, fluids cannot support turbulence and so mixing is often diffusion-limited. Thus, it has been well-documented that mixing can be an obstacle to the efficacy of many microfluidic devices (41). In order to homogenize a fluid over a length l, diffusion requires a time scale t_m = l²D₀. A driven process that generates an enhanced effective diffusivity D_eff will therefore lead to shorter mixing times if the intrinsic diffusivity of the species D₀ is less than D_eff.

To investigate particle transport in this flow regime we analyzed the relative dispersion of tracer pairs with separations given by R(t) = x₂(t) - x₁(t), which has the advantage of suppressing correlated motion of pairs due to uniform advection (42). Mass transport can be characterized by the mean square displacement (MSD) of an ensemble of pair separations, which for motion in two dimensions is given by [5]where τ is a lag time, γ is the exponent of the time dependence, and is the effective relative diffusivity. The MSD of a diffusive process has a linear time dependence, or γ = 1, while ballistic motion results in γ = 2. By a linear least-squares fit to log[〈R²(τ)〉] vs. log(τ), we find that 〈R²(τ)〉 grows in time with γ ≈ 1 at all heights above the sample floor (Fig. 4B Inset). Thus, despite the dominance of advection in the individual tracer paths as represented by the Péclet number, the temporal evolution of the relative quantity R²(τ) can be characterized as an effective diffusion process with a transport rate given by .

Below the tips, the biomimetic cilia motion produces a maximum enhancement of 25 in in relation to the expected relative diffusivity (42) of 500 nm tracers (Fig. 4B). This enhancement would produce a commensurate decrease in the mixing time. We note that the relative diffusivity of the directed transport regime above the cilia tips is equal to that expected for the particles’ intrinsic diffusive motion because the uniform transport is suppressed in the relative dispersion calculation. As an additional control, the expected relative diffusivity was accurately obtained for diffusing tracers tracked while the biomimetic cilia were motionless (as in Fig. 3C).