Fluid Mechanics of Mosaic Ciliated Tissues

In tissues as diverse as amphibian skin and the human airway, the cilia that propel fluid are grouped in sparsely distributed multiciliated cells (MCCs). We investigate fluid transport in this “mosaic” architecture, with emphasis on the trade-offs that may have been responsible for its evolutionary selection. Live imaging of MCCs in embryos of the frog Xenopus laevis shows that cilia bundles behave as active vortices that produce a flow field accurately represented by a local force applied to the fluid. A coarse-grained model that self-consistently couples bundles to the ambient flow reveals that hydrodynamic interactions between MCCs limit their rate of work so that they best shear the tissue at a finite but low area coverage, a result that mirrors findings for other sparse distributions such as cell receptors and leaf stomata.

An indication of the importance of fluid mechanics in biology is the remarkable degree to which the structure of eukaryotic cilia has been conserved over the past billion years [1,2].These hairlike appendages provide motility to microorganisms [3,4] but also direct fluid flow inside animals during development [5][6][7] and in mature physiology in areas from the reproductive system [8] to the brain [9].The two extremes of this organisimal spectrum have a fundamental distinction.In unicellulars like Paramecium, cilia are uniformly and closely spaced on the cell surface [10], while in animals they are often grouped together in dense bundles on multiciliated cells (MCCs) [11] that are sparsely distributed on large epithelia, as in the trachea and kidney [12,13].This difference reflects the need in animal tissues to share surface area with cell types having other roles, such as mucus secretion.
The workings of cilia bundles and the significance of their sparse "mosaic" pattern for fluid transport have only begun to be investigated, primarily limited to in vitro or ex vivo studies [14][15][16].Here we address the fluid mechanics of mosaic tissues using embryos of the amphibian Xenopus laevis in which, by analogy to human airways, cilia driven flow sweeps away mucus and trapped pathogens (Fig. 1).To date, the flow has served as a readout of cilia beating in the study of tissue patterning and cilia disorders [17][18][19]; here we take advantage of the geometry of Xenopus embryos to obtain side views of cilia bundles and quantify the flows they drive.As those cilia collectively sweep through cycles consisting of an extended "power" stroke and compact "recovery" stroke close to the surface [20], the flow within each bundle appears as an active vortex.While the flow driven by a single such vortex decays quickly with distance from the skin, a coarse-grained model shows that long range contributions of other bundles slows the decay of this endogenous flow and determines the shear stress at non ciliated cells.From measurements of beating changes induced by exogeneous flows, we determine linear response coefficients describing the coupling between forces applied by bundles and the flows they generate; we find that hydrodynamic interactions between MCCs lead to maximization at low area coverage of shear at the intervening tissue.These results thereby suggest an explanation for the low area coverages observed in nature.
The epidermis of Xenopus has strong similarities with human mucociliary epithelia.
The sparsely located MCCs whence emanate hundreds of cilia that drive a homogeneous anterior-to-posterior, A-P or head-to-tail, flow (Fig. 1) are surrounded by non-ciliated cells secreting mucus-like material [23], including "goblet cells" that cover most of the tissue, mosaically scattered small cells [24,25] secreting serotonin vesicles that modulate the cil- iary beat frequency [26], and ionocytes transporting ions important for homeostasis.Wild-type Xenopus laevis embryos were obtained via in vitro fertilization [27,28], and grown in 0.1× Modified Barth's Saline at room temperature (or 15 • C to reduce the growth rate, if required).They were imaged at stage 28 [22] after treatment with a minimal dose of anaesthetic (∼ 0.01% Tricaine) to avoid twitching (without affecting cilia dynamics [19]).Embryos lie on one of their flat flanks at this stage, providing a side view of cilia bundles of ventral MCCs (Fig. 1), whose power strokes are in the A-P direction (left to right in figures) so cilia and the flows stay mostly within the focal plane.
In flow chamber experiments, embryos were perfused with a peristaltic pump while in a Warner Instruments chamber (RC-31A): a 4 mm × 37 mm channel cut into a 350 µm thick silicon gasket sandwiched between two coverslips that keep the embryo in place by pressing against its sides.The dorsal part of the embryo was positioned closer to the chamber wall, the anterior region of interest was > 2 mm away, and the A-P axis parallel to the channel axis, the main direction of the perfusing flow.Brightfield images of cilia and 0.2 − 0.5 µm tracers (mass fraction ∼ 0.01%) were acquired at 2000 frames/s for ≥ 1 s by a high speed camera (Photron Fastcam SA3) on an inverted microscope (Zeiss Axio Observer) with a long distance 63× objective (Zeiss LD C-Apochromat).Images were filtered by subtracting their moving average.Flow fields u = (u, v, w) were estimated by Particle Image Velocimetry (PIVlab) and averaged over time.
We set the stage by summarizing the important length and time scales.MCCs are spaced apart by 40 − 80 µm and uniformly distributed with average density P ≈ 2.6 × 10 −4 µm −2 , giving an average spacing d = 1/P ∼ 62 µm.With ∼ 15 µm ( ¯ = 14.52 ± 0.21µm) the cilia length, the average cellular area ∼ 287 ± 11 µm 2 is ∼ 2 , which gives a coverage fraction φ = ( /d) 2 ∼ 0.07.The cilia on MCCs beat at a frequency f ∼ 20 − 30 Hz and during a power stroke their tips move a distance ∼ 2 [Fig.2(a)] in half a period, reaching speeds V c ∼ 4f ∼ 1 mm/s, so the Reynolds number ρV c /µ (with ρ the density and µ the viscosity of water) is ∼ 0.01, well in the Stokesian regime.Using the fluid speed u v ∼ 0.5 mm/s between vortices as typical of the periciliary region, the Péclet number u v /D > 1 even for small molecules.
Cilia within an MCC are not synchronized; their tips move in a tank-treading manner [28], generating vorticity ω e y perpendicular to the beating plane.Each MCC is thus an active vortex, as seen in Figs.2(b,c).The vorticity can exceed ∼ 150 s −1 ∼ 2V c / , is colocalized with the cilia, and rapidly diffuses at larger z as the flow becomes parallel to the skin.Above non-ciliated cells between MCCs, there is a shear flow for z < , while further away (z d), the discreteness of the MCCs is washed out by viscosity and the horizontal velocity u is independent of x and falls off slowly with z [Figs.2(e,f)].
The first step toward understanding the coupling between cilia beating and fluid flow involves quantifying the contribution of a single MCC.We introduce a boundary ∂Ω c enclosing the volume Ω c of the active vortex [Fig.2(a)], and extrude it in y ∈ [−10 µm, 10 µm], the measured size of the vortex.The Stokes equations in the complement Ω c , are solved using an envelope approach [29,30] in which the dynamics of the cilia tips determine the flow, noting that the flow from an isolated cilium is well-approximated in the far field as that of a point force (Stokeslet) [31].We position N Stokeslets at s n = (x n , y n , h n ) in Ω c and find their strengths f n = (f n,x , 0, f n,z ) by fitting the velocity on ∂Ω c and a no-slip boundary at z = 0 [Fig.2(d)].This gives an estimate of the flow u 0 driven by a bundle in an otherwise quiescent fluid.The x-components are of interest, as they alone contribute to net AP flow.For z > 2l, the component u 0 falls off for z h like the flow u s (z) ∼ 3F h 2 /4πµz 3 above a single Stokeslet parallel to and a distance h above a no-slip wall [32] (Fig. S1 [28]).The data in Figs.2(e,f) show a much weaker fall-off than z −3 for z > 2l.This slow decay arises from the long range contribution of more distant MCCs, as we now show.
The flow u s due to a force f x can be approximated by its far field limit.In cylindrical coordinates (ρ, z) centered at the bundle, (1) The Stokeslet contributions can be lumped into an effective force F c = −1 n h n f n,x , which, applied at z = , matches the far field u s ≈ F c S11 .This also applies for an effective local moment (rotlet) Γ c = 2 F c [33].
The observed velocity u(z)e x in the region z > d is independent of x (Fig. S2 [28]) and can be described most simply in a model of the skin as a uniform distribution of x-directed Stokeslets with area density P. Summing the contributions of all bundles up to a cutoff radius Λ that incorporates finite embryo size, we obtain which has the scaling form u far = V G(z/Λ), with and For any fixed z, as the organism size Λ → ∞, χ → 0, giving a flow independent of z with speed V [34], while for any fixed lateral scale Λ, the asymptotic flow field vanishes as z/Λ → ∞.For the fitted parameters V 0.64 mm/s and Λ = 300 µm, u f ar provides an almost perfect fit of the data in Figs.2(e,f) for z > 2d (also Fig. S2 [28]).Direct summation of discrete Stokeslets on a lattice produces nearly identical results, validating the approximation of a continuous distribution for the far-field flow (Fig. S3 [28]).For this value of V and the observed density P, we have the farfield estimate F c 159 pN.A slightly larger effective force F 0 200 pN is obtained by solely fitting the velocity at the cilia tips (Fig. 2d), which is to be expected as this approach does not account for the endogenous flow u a e x to which a bundle is self-consistently exposed.
We estimate the endogeneous ambient flow u a by subtracting from u far the contribution from a single cilia bundle, taken as a distribution of radius d/2, As u c = u − u a represents the contribution of a single bundle, we test the far-field estimate with a near-field fit of u c in Ω c by means of Stokeslets within Ω c .Doing so for the volume (−15 < x < 15, −10 < y < 10, 0 < z < d), the superposition u c + u a gives an excellent fit to the data in Figs.2(e,f), above and between bundles, with F c 160 pN, nearly identical to the far-field estimate.
For comparison, the average lateral force generation over one cycle of beating (assuming only the power stroke contributes) can be estimated from resistive force theory [35] 3 pN [28], where ζ ⊥ = 4πµ/| ln( √ eε)| is the transverse drag coefficient for a slender filament of aspect ratio ε (for cilia, ε ∼ 75).We infer from the estimated F c that the effective number of cilia contributing to the Stokeslet is ∼ 53, about half the typically ∼100 cilia in an MCC, reflecting force cancellations from phase shifts between cilia.
The fact that F c /F 0 < 1 shows that it is necessary to incorporate the response of a cilia bundle to an ambient flow.We probed this experimentally by exposing the bundle to an exogenous shear flow γe ze x in a flow chamber (Fig. 3).When γe = 0, the cilia tips move with velocity V 0 and drag the fluid, generating a negative shear rate γ0 ≈ −23 s −1 .Pumping fluid in the same direction, the shear rate at the tips γc decreases linearly with the hydrodynamic load γe .The corresponding velocity V c tends to increase, but at a much slower rate.The rate of work above the cilia tip envelope ∝ − γc V c thus decreases almost at the same rate as γc , and for γe /V 0 > 0.3 becomes negative, consistent with a dissipative bundle.
The lateral velocity u γ e x just above the bundle (z < 2.5 ) is well fitted by the linear combination u γ ≈ C û0 + γ e z (Fig. 3c) with û0 ( F0 ) the profile at γe = 0. Note that C ≈ F/ F0 , as confirmed for the above calculations where u c ≈ (F c /F 0 )u 0 for z > (Fig. S2(b) [28]).The slope of F/ F0 versus γe /V 0 in Fig. 3, would be unity if the bundle dynamics were fully preserved (V c = V 0 ), and zero if the bundle's force were constant.The measured slope 0.76±0.06confirms the resistive behavior and allows us to parameterize the coupling of the cilia to the ambient flow by the linear relation F ≈ F0 −α γe , with α = 0.76 F0 /V 0 .
To close the loop on a self-consistent coupling of the bundles to the flow, we now replace γe with the endogenous shear γa = ∂u a /∂z| z=0 , giving F ≈ F 0 − α γa , with for large tissues, we have γa ≈ 3F /µd 3 , and thus where the λ = 3α/ µ is the key parameter of the selfconsistent theory.Using typical values ( F0 = 160 pN, V 0 = 1.3 mm/s, = 15 µm), we find λ ≈ 18.6.
The relation ( 5) can be used to address several aspects of sparse cilia distributions [28].The force applied to the wall per bundle is ∼ µ γa d 2 = 3F 0 a/(1 + λa 3 ), with a = /d, and has a maximum at d max = (2λ) 1/3 ≈ 50 µm as does the contribution F V of a single bundle to the rate of work.The force F w ≈ µ γa (d 2 − 2 ) applied to the non-ciliated cells is maximal for d ≈ 54 µm.Both values are in excellent agreement with those observed.Using the relation φ = ( /d) 2 we can express the wall force as a function of the coverage fraction φ, The contour plot of F w /F 0 in the φ−λ parameter space in Fig. 4 shows that the optimum area fraction is a strongly decreasing function of λ and can reach values far below unity for λ ∼ 20, as in the present study.Extra endogenous loads γe < 0, as expected for internal tissues, appear as an additional contribution to the ambient flow γa + γe .These loads will contribute to (5) as lower values of λ, and indeed, consistently larger, yet still low coverage fractions of the airways of several animals have been estimated to be 0.4 − 0.5 [14], qualitatively consistent with Fig. 4. We infer that the observed mosaic patterns are close to optimal in terms of the clearing force applied to non-ciliated cells.We close with comments on connections to other systems with sparse distributions of active elements.The force applied to the outer fluid by the cilia tips on the envelope ∂Ω c is equal and opposite to that applied to the skin, and we can simplify our results on the shearing of non-ciliated cells by reconsidering (2) as the flow of a patch of activity with given slip velocity V of radius Λ.The shear stress driving the flow is τ Λ ≈ 3µV /2Λ, which we assume constant over a bundle.Setting Λ = and integrating over a tissue with N bundles we obtain J ∼ N 2 τ .By contrast, if we set Λ = R, the local shear stress is τ R = 3µV /2R and the force over the entire surface is J R ∼ πR 2 τ R .With N = πR 2 φ/ 2 , the ratio measures how well a distribution of non-interacting MCCs shears the surface relative to the collection.The linear scaling of ( 7) with φ is expected, but the large prefactor R/ ∼ 20 (system size/MCC size) implies that J/J R can approach unity for area fractions as low as φ ∼ /R ∼ 5%.The form of this result mirrors one found by Jeffreys [36] for the evaporation rate from sparsely distributed leaf stomata, rediscovered years later [37] in the context of ligand binding to sparse cell receptors [38].
The results presented here suggest that long range hydrodynamic interactions between multiciliated cells allow efficient peri-ciliary transport at relatively low coverage, favoring the coexistence of multiple cell types in large tissues.This is likely just one aspect of more general mechanisms that maintain efficient transport in the upscaling events marking the evolutionary transition from unicellular to larger multicellular systems.
This This file contains additional experimental and calculational details.
Embryo Culture: Xenopus embryos were prepared as described previously [S1].Briefly, mature Xenopus laevis males and females were obtained from Nasco [S2].Females were injected with 50 units of pregnant mare serum gonadotropin 3 days in advance and 500 units human chorionic gonadotropin 1 day in advance in the dorsal lymph sack to induce natural ovulation.Eggs were laid in a 1× MMR buffer (5 mM HEPES pH 7.8, 100 mM NaCl, 2 mM KCl, 1 mM MgSO 4 , 2 mM CaCl 2 , 0.1 mM EDTA).Xenopus embryos were cultured at room temperature or 15 • C in the 0.1× MMR until they reached stage 27/28.Experiments with embryos were performed at the late tailbud stages (stages 28-30, as describe in Faber and Nieuwkoop [S3]).Embryos were terminated humanely immediately following the experiments.
Our work with Xenopus laevis is covered under the Home Office Project License PPL 70/8591 and frog husbandry and all experiments were performed according to the relevant regulatory standard.All experimental procedures involving animals were carried out in accordance with the UK Animals (Scientific Procedures) Act 1986.Moreover, we only used surplus embryos for this study, to conform with the NC3Rs guidance to exploit the possibility to minimise the use of animals by sharing embryos with collaborators.
Statistics: To fit any quantity y measured at a given hydrodynamic load x, we assume a linear relation y = a + bx and find 95% confidence intervals for the averages ā and b of the parameters a i and b i given by a leastsquare fit of the measured values (x i,m , y i,m ) acquired for the i-th MCC.We have [S4]: The standard errors for the mean parameters are the square roots of while the parameter variances for the i-th MCC are with The prime superscript indicates that the variance of the measurements σ 2 i,m = σ 2 i for the i-th MCC was assumed to be constant.It was estimated as the variance s 2 of the sample population: We normally imaged M = 4 conditions per MCC.For exceptions with M = 2, s 2 could not be computed directly from (S5) and was assumed to equal the largest value from the other experiments.
Fitting the near flow-field by the singularity method: The flow u c driven by the cilia in Ω c is modelled as the superposition of the flows arising from local point forces (Stokeslets) f n applied at s n ∈ Ω c : The tensor S is the well-known, exact solution for a Stokelset next to a no-slip plane at z = 0 [S5].The values of f n are found by fitting u c at M collocation points b i , with M > 2N to avoid numerical instabilities [S6].
As no-slip boundary conditions u c = 0 at z = 0 are implicitly satisfied, walls do not need to be discretized.The linear system (S6) is then simply recast in its matrix form Af = u b , with the 3M × 2N matrix We then solve for f using the backslash operator of Matlab.We set f 2 = 0, assuming the solution to be symmetric in y.Once f n are known, (S6) can be used to evaluate the fitted solution at any x.The flow field measured in the y plane is extruded by replications at 13 planes evenly spaced between −10 µm < y < 10 µm.We used 15 Stokeslets for each plane about the fictitious boundary ∂Ω c .
Coarse-graining the bundle: The flow u c driven by the Stokeslet in the bundle can be coarse-grained further, with a smaller number of Stokes flow singularities, moving away from the bundle.We compare the fitted flow u c , made up of N Stokeslets as discussed above, with the flow driven by the effective Stokeslet F c e x applied at (0, 0, ), and the effective rotlet 2 F c e y applied at (0, 0, /2).They share the same far-field, reflecting the fact it is driven by an active vortex.The flow driven by the entire bundle decays as 1/z 3 , as for a single singularity, for z > 2 (Fig. S1).
Far Field fitting: The flow given by (2), is used to fit the PIV measurements for z > 2d (Fig. S2).The velocity u f ar (z; Λ) depends linearly on V , but not on R. For a given value of R, we find V by a linear least-squares fit of the data.We then simply repeat this linear fit for candidate values in the range 30 µm < R < 1 mm, with increment ∆R = 10 µm, and select the value of R that minimizes the L 2 fitting error.Using the effective force F c estimated by the near field fitting, and summing up the exact contribution of each MCC, we retrieve the slow decay rate observed in vivo [Figs.2(e,f)] for Λ ∼ 300 µm (Fig. S3).This is the same result found by fitting the far-field flow with Eq. (2).
Resistive force theory estimate of the effective force applied by a single cilium: We adopt a simplified view of the power stroke of a cilium as a straight rod that pivots around its base.Let s ∈ [0, ] be arclength along a cilium, with s = 0 at the base and s = at the tip, and let φ be the angle between the cilium and the wall.The lateral component of the RFT force density at s is f ∼ (s/ )ζ ⊥ V c sin φ, and the resulting far field velocity, given by Eq. (1), is proportional to hf ds with h = s sin φ.Accordingly, the effective force f φ = −1 0 h(s, φ)f (s, φ)ds, matches the overall far field f φ S11 when applied at .We obtain f φ = sin 2 φ ζ ⊥ V c /3. Through the entire stroke, a cilium cycles through an angle ∆φ = 2π, and we assume that the recovery stroke does not contribute to the force, so f φ = 0 for π < θ < 2π.Averaging f φ gives the effective force f = ζ ⊥ V c (6π) −1 π 0 sin 2 φ dφ, which gives the expression f = ζ ⊥ V c /12 used in the main text.
Additional results from the self-consistent model: The self-consistent model in Eq. ( 5) can be

FIG. 3 .
FIG. 3. Response of a cilia bundle to shear flow.(a,b) Vorticity contours and velocity vectors before and during perfusion.(c) Velocities û0 and uγ, exogenous shear flow γez, and linear combination (F/ F0) û0 + γez fitting uγ.(d) Linear fits of the variation of estimated force F , velocity Vc and shear rate γc = ∂u/∂z measured above cilia tips, and rate of work W ∝ γcVc (overlaping γc), normalized by values at γe = 0. Shaded regions are 95% confidence intervals of averages over 10 samples.

9 FIG. 4 .
FIG. 4. Force on non-ciliated cells in the self-consistent model.Contour plot of (6) in parameter space.Dashed line traces optimization ridge.

FIG. S1 .
FIG.S1.The contribution of a single bundle of cilia decays as a single effective force Fc or moment 2 Fc, while the measured flow profile decays much more slowly due to the contributions from other MCCs (exp).The lateral velocity u0, with reference to Fig.2(d), is shown (a) above the bundle as a function of z, and (b) between bundles at z = , as a function of x.
FIG. S3.The lateral velocity of a two-dimensional array of Stokeslets of radius Λ ∼ 300, falls off as in experiments (exp), as shown for: (a,b) above and between bundles, respectively, as a function of z; (c) between bundles as a function of x. Results are obtained by direct summation of the exact solution of each Stokeslet, all of strength Fcex and z-offset .

FIG. S4 .
FIG. S4.Supplement of Fig. 4. Contour plot of (a) the effective force F [Eq. (5)], and (b) the limit velocity V , in parameter space.Dashed red line in (b) traces ridge of optimization for V .
Two-dimensional array of Stokeslets: Results similar to those presented in Fig. 2 for a uniform distribution of Stokelets can be obtained by positioning Stokeslets F c e x on a lattice with cut-off radius Λ.Each element s ij of the lattice is position at (x ij = id 11 + jd 12 , y ij = jd 22 , z ij = ).From confocal imaging of the closest neighboring cells of the bundle in Fig. 2, we estimate d 11 ∼ 70 µm, d 12 = 40.5 µm, d 22 = 53µm.