Fluid Mechanics of Mosaic Ciliated Tissues

Francesco Boselli<sup>1,*</sup> Jerome Jullien<sup>2,3,4,†</sup> Eric Lauga<sup>1,‡</sup> and Raymond E. Goldstein<sup>1,§</sup>

<sup>1</sup>Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge,
Cambridge CB3 0WA, United Kingdom

<sup>2</sup>Wellcome Trust/Cancer Research UK Gurdon Institute, Tennis Court Road, Cambridge CB2 1QN, United Kingdom

<sup>3</sup>Department of Zoology, University of Cambridge, Cambridge CB2 1QN, United Kingdom

<sup>4</sup>Inserm, Nantes Université, CHU Nantes, CRTI-UMR 1064, F-44000 Nantes, France

<figure><img src="image_1.png" /></figure>

(Received 30 March 2021; accepted 7 October 2021; published 2 November 2021)

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.

DOI: 10.1103/PhysRevLett.127.198102

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–7] and in mature physiology in areas from the reproductive system [8] to the brain [9]. The two extremes of this organismal 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–16]. Here we address the fluid mechanics of mosaic tissues using embryos of the amphibian *Xenopus laevis*. The cilia covering the epidermis drive a flow that sweeps away pathogens, preventing infection, and possibly contributing to oxygenation (Fig. 1) [17]. To date, the flow has served as a readout of cilia beating in the study

of tissue patterning and cilia disorders [18–20]; 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 [21], 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, we show that long range contributions of
other bundles slows the decay of this *endogenous* flow and
determines the shear stress at nonciliated cells. From
measurements of beating changes induced by *exogeneous*

<figure><img src="image_2.png" /><figcaption>FIG. 1. Ectoderm of embryonic *Xenopus laevis* at tailbud stages. (a) Schematic side view of MCCs (red) intermixed with secreting cells. (b) Location of MCCs across the embryo (adapted from Refs. [22,23]) and cilia-driven flow (blue arrows). (c) Confocal image of cell membranes (stained by membrane-RFP), with MCCs segmented in red, in ventral region of skin.</figcaption></figure>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Flow fields. (a) Lateral view of MCC showing (dashed) path of cilia tips and force F. (b) Experimental velocity field and vorticity in plane normal to skin near several MCCs. (c) Near an MCC, as in (b), with direction of cilia tip motion (black arrows) on $\partial\Omega_c$. (d) Estimated flow field $u_0$ for an isolated MCC (blue arrows): Point forces (red arrows) are used to fit velocity near cilia tips. Lateral velocity at (e) $x, y = (0, 0)$ and (f) $(\pm 40 \ \mu\text{m}, 0)$ in experiment (exp) and theory, with $u_0$ driven by an isolated bundle and $u_c$ by a bundle exposed to endogenous flow $u_a$ (see also Figs. S2–S4 [29]).</figcaption></figure>

flows, we determine the linear self-consistent relation describing the dissipative coupling between forces applied by bundles and the flows they generate; we find that bundles, as a collective, most efficiently shear the inter-vening tissue at the observed low area coverage. These results thereby suggest an explanation for the low area coverages observed in nature.

The *Xenopus* epidermis has strong similarities with human mucociliary epithelia [17]. The sparsely located MCCs whence emanate hundreds of cilia that drive a homogeneous anterior-to-posterior (A-P) flow (Fig. 1) are surrounded by nonciliated cells: “goblet cells” that cover most of the tissue secreting mucus-like material [24], mosaically scattered small cells [25,26] secreting serotonin vesicles that modulate the ciliary beat frequency [27], and ionocytes transporting ions important for homeostasis. Having evolved much earlier, *Xenopus* epidermis is simpler and more accessible than most vertebrate mucociliary epithelia. Unlike in airways, the fluid domain extends far above the cilia, and thus it is a model for fundamentals of mucociliary systems and their evolution.

Wild-type *Xenopus laevis* embryos were obtained via *in vitro* fertilization [28,29], 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 [23] after treatment with 0.01% Tricaine to avoid twitching without affecting cilia dynamics [20]. Stage 28 embryos lie on one of their flat flanks, providing a side view of cilia in 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. The anterior region of interest was > 2 mm away from the chamber wall, and the A-P axis parallel to 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 sub-tracting their moving average. Flow fields $\mathbf{u} = (u, v, w)$ were estimated by particle image velocimetry (PIVlab) and averaged over time. Pulsatility of the flow is low (fluctua-tions < 10% of the mean [29]), so we treat the flow as steady.

We first summarize important length scales and time-scales. MCCs are spaced 40–80 µm apart and uniformly distributed with density $P \approx 2.7 \pm 0.73 \times 10^{-4} \ \mu\text{m}^{-2}$, from which we define a typical spacing $d = \sqrt{1/P} \sim 61 \ \mu\text{m}$. With $\ell = 14.52 \pm 0.21 \ \mu\text{m}$ the cilia length, the average cellular area of $287 \pm 11 \ \mu\text{m}^2$ is $\sim \ell^2$, and the coverage fraction $\phi \approx 0.075 \pm 0.02$ is $\sim (\ell/d)^2$. Cilia beat at $f \sim 20–30 \ \text{Hz}$; during a power stroke their tips move $\sim 2\ell$ in half a period, reaching speeds $V \sim 4f\ell \sim 1 \ \text{mm/s}$. The Reynolds number $V\ell/\nu \sim 0.01$ (with $\nu$ the kinematic viscosity of water) is in the Stokes regime.

Despite the tendency of nearby cilia to synchronize [4], those in an MCC are not in phase; a snapshot [Figs. 1(a) and 2(a)] shows cilia throughout the beat cycle [29], generating vorticity $\omega \parallel \mathbf{e}_y$ perpendicular to the beating
plane. Each MCC is thus an active vortex [Figs. 2(b),2(c)]: $\omega$ can exceed ~150 s⁻¹ ~2V/ℓ, is colocalized with the cilia, and rapidly diffuses at larger z as the flow becomes parallel to the skin. Above nonciliated cells, there is a shear flow for $z < \ell$, while for $z > d$ the discreteness of the MCCs is washed out and u is independent of x (Fig. S3 [29]) and falls off slowly with z [Figs. 2(e),2(f)].

A first step toward understanding the coupling between cilia and flow involves quantifying the contribution of a single MCC. We introduce a boundary $\partial\Omega_c$ enclosing the volume $\Omega_c$ of the active vortex [Fig. 2(a)], and extrude it in $y \in [-10 \ \mu m, 10 \ \mu m]$, the measured size of the vortex. The far-field flow from an isolated cilium is well approximated by that of a point force [34], and the velocity $\mathbf{u}_s$ driven by a point force $\mathbf{f}$ at s next to a no-slip wall is $\mathbf{u}_s(\mathbf{x}) = \mathbf{f} \cdot \mathbf{S}(\mathbf{x}, \mathbf{s})$ [29], with $\mathbf{S}$ the Blake tensor [35]. We model the bundle with N point forces $\mathbf{f}_n = (f_{n,x}, 0, f_{n,z})$ at $\mathbf{s}_n = (x_n, y_n, h_n)$ in $\Omega_c$, from which we define a single effective lateral force $F = \ell^{-1} \sum_n h_n f_{n,x}$, applied at $z = \ell$ that drives a lateral flow $u_s$ matching the far field of the entire bundle. This also equals the far-field velocity for a local moment (rotlet) $\Gamma = 2\ell F$ [36].

We first use an envelope approach [37,38] in which the cilia tips determine the flow, and find $\mathbf{f}_n$ by fitting the measured velocity on $\partial\Omega_c$ and a no-slip boundary at $z=0$ [Fig. 2(d)]. The flow driven by these forces, $\mathbf{u}_0(\mathbf{x}) \approx \sum_n \mathbf{f}_n \cdot \mathbf{S}(\mathbf{x}, \mathbf{s}_n)$, is in an otherwise quiescent fluid. This fit yields a force $F_0 \approx 209 \pm 14$ pN. For $z > 2\ell$ and $x=0$, the component $u_0$ falls off like the flow $u_s(z) \sim 3F_0\ell^2/4\pi\mu z^3$ above a single point force [35] (Fig. S2 [29]). The data [Fig. 2(e)] show a much weaker decay, which we now show arises from contributions of distant MCCs.

The simplest model for the ambient flow $u_a$ for $z > d$ assumes a uniform distribution of x-directed forces $F$, one for each MCC, with area density $\mathcal{P}$. In cylindrical coordinates $(\rho, \theta, z)$ centered on a bundle, the lateral flow due to any one force is $u_s(x) \approx F\tilde{S}_{xx}$, where $\tilde{S}_{xx}$ is the far-field limit of $S_{xx}$,

$$ \tilde{S}_{xx}(\rho, \theta, z) = \frac{3\ell}{2\pi\mu} \frac{z\rho^2\cos^2\theta}{(\rho^2 + z^2)^{5/2}}. \quad (1) $$

Integrating up to a radius $\Lambda$ that represents finite embryo size and subtracting that from the bundle at the origin, taken as a distribution of radius $d/2$, we obtain

$$ \begin{aligned} u_a(z) &= \mathcal{P}F \int_0^{2\pi} \int_{d/2}^{\Lambda} \rho d\rho d\theta \tilde{S}_{xx}(\rho, \theta, z) \\ &= U[G(z/\Lambda) - G(2z/d)], \end{aligned} \quad (2) $$

where $G(\chi) = 1 - (3\chi + 2\chi^3)/2(1 + \chi^2)^{3/2}$ and $U = \mathcal{P}Fl/\mu$. G decreases monotonically from $G(0) = 1$ to $G(\infty) = 0$. For any fixed $z \gg d$, as the organism size

$\Lambda \to \infty$, $\chi \to 0$, giving a z independent flow with speed $U$ [39], while for any fixed $\Lambda$, $u_a$ vanishes as $z/\Lambda \to \infty$.

A fit of $u_a$ to the data in Fig. 2 yields $U = 0.63 \pm 0.05$ mm/s and $\Lambda = 308 \pm 30$ µm [29], and using the observed $\mathcal{P} \approx 2.6 \pm 0.5 \times 10^{-4}$ µm⁻² we obtain the far-field estimate $F \approx 163 \pm 35$ pN. The fact that $u_a \sim u$ for $z > d$ confirms the slow decay above a bundle is established by distant MCCs [Figs. 2(e),2(f) and Fig. S3 in the Supplemental Material [29]). Direct summation of a lattice of point forces yields nearly identical results, validating the continuum approximation (Fig. S4 [29]).

Thus, we may model the flow as $\mathbf{u}(\mathbf{x}) \approx \mathbf{u}_c(\mathbf{x}) + u_a\mathbf{e}_x$, where $\mathbf{u}_c(\mathbf{x}) \approx \sum_n \mathbf{f}_n \cdot \mathbf{S}(\mathbf{x}, \mathbf{s}_n)$ is the local contribution, with the $\mathbf{f}_n$ found by fitting $u - u_a$. We obtain $F \approx 162 \pm 30$ pN, confirming the far-field estimate. Figures 2(e),2(f) show that $u_c + u_a$ gives an excellent fit to the data for $z > \ell$, above and between bundles.

The magnitude of $F$ can be compared to the average lateral force generated by the power stroke of a cilium estimated from resistive force theory [40] as $f \sim \zeta_{\perp}\ell V/12 \approx 3$ pN [29], where $\zeta_{\perp} = 4\pi\mu/|\ln(\sqrt{e}\epsilon)|$ is the transverse drag coefficient and $\epsilon \sim 75$ is the cilium aspect ratio. Thus, a bundle contributes with ~50–60 cilia, about half the typically ~100 cilia in an MCC, reflecting force cancellations from phase shifts between cilia.

To probe the ciliary response to flow, we exposed the epidermis to an exogenous shear flow $\dot{\gamma}_e z \mathbf{e}_x$ (Fig. 3). When $\dot{\gamma}_e = 0$, the cilia drag the fluid, generating a negative shear $\dot{\theta}(0) = \partial u(z)/\partial z|_{\ell} \sim -23$ s⁻¹ at the tips. Pumping fluid in the same direction, $\dot{\theta}$ decreases linearly with $\dot{\gamma}_e$, while the velocity $V$ at the tips increases, but much more slowly. The rate of work above the cilia tips $W \propto -\dot{\theta}V$ thus decreases almost at the same rate as $\dot{\theta}$, and for $\dot{\gamma}_e \ell/V(0) > 0.3$ becomes negative, consistent with a dissipative bundle. This reflects, in part, the dense packing of cilia in MCCs; the weak coupling is consistent with weak entrainment of dense bundles of ependymal cilia by exogenous flow [41].

Figure 3(c) shows that the lateral velocity $u(z; \dot{\gamma}_e)$ above the bundle ($z < 2.5\ell$) is well fit by $u(z; \dot{\gamma}_e) \approx Cu(z; 0) + \dot{\gamma}_e z$. Here, $C \approx F/F(0)$, as confirmed by the above calculations, where $u_c \approx (F/F_0)u_0$ for $z > \ell$ [Fig. S3(b) in Ref. [29]]. It follows that $V$ and force are coupled as $V = [F/F(0)]V(0) + \dot{\gamma}_e \ell$, and the slopes of the curves $F/F(0)$ and $V/V(0)$ versus $\dot{\gamma}_e \ell/V(0)$ are $-\alpha$, and $1 - \alpha$, with $\alpha$ an empirical parameter. If $\alpha = 1$, the bundle compensates for $\dot{\gamma}_e$ and the velocity at the tips is preserved, whereas if the force were constant ($\alpha = 0$), the bundle supplies no resistance, and the ambient flow simply adds to that of the bundle. The measured slope $\alpha \approx 0.76 \pm 0.06$ (Fig. 3) confirms resistive behavior.

The response $F-F(0)$ of the bundle to $\dot{\gamma}_e$ can be parametrized by a drag coefficient (viscosity × length), analogous to the Stokes drag factor $6\pi\mu R$ for a sphere of radius $R$, via the relation $\partial F/\partial(\dot{\gamma}_e \ell) = -\lambda_e \mu \ell$, with $\lambda_e = \alpha F(0)/\mu \ell V(0)$. For a bundle applying a force $F_0$ in
<figure><img src="image_4.png" /><figcaption>FIG. 3. Response of cilia bundles to exogenous shear. (a),(b) Vorticity and velocity vectors before and during perfusion. (c) Velocities $u(z; 0)$ and $u(z; \dot{\gamma}_e)$, shear flow $\dot{\gamma}_e z$, and sum $Cu(z; 0) + \dot{\gamma}_e z$ fitting $u(z; \dot{\gamma}_e)$, where $C = F(\dot{\gamma}_e)/F(0)$. (d) Variation with shear of estimated force $F$, velocity $V$ and shear rate $\dot{\theta} = \partial u(z)/\partial z|_z$ measured above cilia tips, and rate of work $W \propto \dot{\theta}V$ (overlapping $\dot{\theta}$), normalized by values at $\dot{\gamma}_e = 0$. Shaded regions are 95% confidence intervals of averages over 10 samples.</figcaption></figure>

an otherwise quiescent fluid, we infer an equivalent response $F - F_0$ to a general shear rate $\dot{\gamma}$, and drag coefficient $\partial F/\partial(\dot{\gamma}^\ell) = -\lambda\mu\ell$. Bundles experience the sum $\dot{\gamma} = \dot{\gamma}_e + \dot{\gamma}_a(F)$, where $\dot{\gamma}_a(F) = \partial u_a/\partial z|_{z=0}$. The constant $\lambda$ can be related to the empirical parameter $\lambda_e$ by matching derivatives of $F - F(0)$ and $F - F_0$ with respect to $\dot{\gamma}_e$. As $d/\Lambda \ll 1$, we have $\dot{\gamma}_a \approx 3F\ell/\mu d^3$, and thus $\lambda = [1/\lambda_e - 3(\ell/d)^3]^{-1} \approx 8.6 \pm 2.5$ for Xenopus; this is within a factor of 2 of the value $6\pi$ for a sphere.

To obtain the self-consistent model, we set $\dot{\gamma}_e = 0$, and find $F = F_0 - \mu\lambda\ell^2\dot{\gamma}_a(F, d)$. For $\phi = (\ell/d)^2 \ll 1$ we find

$$F(d) \approx \frac{F_0}{1 + 3\lambda(\ell/d)^3}, \quad (3)$$

so the force applied by a bundle decreases with coverage fraction $\phi$, as does the derivative $\partial U/\partial \phi$ of the limit velocity $U \propto F\phi$ in the outer region (2). With $U$ bounded by the cilia tip velocity $V$, the observation that $U \sim V/2$ at $\phi < 0.1$, is indeed consistent with lower rates $\partial U/\partial \phi$ at larger $\phi$. This implies a trade-off between the dissipated force $F_0 - F$, and the benefits of the collective contribution $u_a$, as seen in the extra shear force per bundle $F_w \approx \mu\dot{\gamma}_a(d^2 - \ell^2)$ applied to non-ciliated cells and the extra force per bundle $F_a = F\ell/d$. The latter relates to the power $FU = F_a^2/\mu\ell$ contributing to long-range transport [29]. Periciliary and outer transport are coupled for sparse tissues, as $F_w^2 \sim (\mu\dot{\gamma}_a d^2)^2 = 9\mu\ell FU$. The latter is maximal at $d_{max} = (6\lambda)^{1/3}\ell \approx 56$ µm, suggesting that further increase in coverage is inefficient.

In Fig. 4, we estimate these forces for all $\phi$ using the exact value of $\dot{\gamma}_a$ for a distribution of point forces [(S8,9), [29]]. The optimum of $F_w$ sharply decreases with $\lambda$ and for Xenopus is at $\phi \sim 0.055$, very close to the in vivo value. Increasing $\phi$ primarily increases transport in the outer region, where the collective power $FU$ is maximal at $\phi \sim 0.1$. At $\lambda = 0$, $F_w$ is maximal at $\phi = 0.17$, primarily

because the nonciliated area goes to zero with $\phi$, while $FU$ is maximal at $\phi = 1$. The observed mosaic patterns are nearly optimal based on the clearing force $F_w$.

Similar results arise from the envelope approach to flows driven by arrays of finite-size bundles [29]. For $\phi > 0.33$ these models lose realism as the vortices of size $\sim 2\ell$ start to overlap, but we expect an optimum shearing configuration at $\phi < 0.3$, as transport will be confined to the region above the cilia for larger $\phi$. Similar dissipative phenomena will contribute to systems working at larger coverage fractions, as in airways [14], though synchrony of bundles as well as mucus between cilia and the air above add complexity.

We close with connections to other sparse distributions of active elements. The force applied to the fluid by the cilia tips on the envelope $\partial\Omega_c$ is equal and opposite to that applied to the skin, and we can simplify our results on the shearing of nonciliated cells by reconsidering the flow $u(z) = UG(z/\Lambda)$ in Eq. (2) above a patch of activity of radius $\Lambda$ with slip velocity $U$. The shear stress driving the

<figure><img src="image_5.png" /><figcaption>FIG. 4. Collective efficiency of a distribution of point forces, in the space of coverage $\phi$ and dissipative coupling constant $\lambda$. Dashed lines trace optimization ridges of the extra forces per bundle $F_a$ and $F_w$ driving (a) flow in the outer region and (b) shear above nonciliated cells.</figcaption></figure>
flow is $\tau_{\Lambda} \approx -3\mu U/2\Lambda$, which we assume constant over a bundle. Setting $\Lambda = \ell$ and integrating over a tissue with $N$ bundles we obtain $J \sim N\ell^2\tau_{\ell}$. By contrast, if we set $\Lambda = R$, the local shear stress is $\tau_{R} = -3\mu U/2R$ and the force over the entire surface is $J_{R} \sim \pi R^{2}\tau_{R}$. With $N = \pi R^{2}\phi/\ell^{2}$, the ratio

$$ \frac{J}{J_R} \sim \frac{R}{\ell} \phi \qquad (4) $$

measures how a distribution of noninteracting MCCs shears the surface relative to the collection. The linear scaling with $\phi$ is expected, but the large prefactor $R/\ell \sim 20$ (system size/MCC size) implies that $J/J_R$ can approach unity for $\phi$ as low as $\sim \ell/R \sim 5\%$. This mirrors Jeffreys' result [42] for the evaporation rate from sparse leaf stomata, rediscovered years later [43] in the context of ligand binding to cell receptors [44].

The results presented here show that long-range hydrodynamic interactions between multiciliated cells allow efficient periciliary transport at low coverage, favoring coexistence of multiple cell types. This is likely but one example of the mechanisms that maintain efficient transport in the upscaling events marking evolutionary transitions from unicellular to larger multicellular systems.

This research was funded in whole, or in part, by Wellcome Trust Grants No. 101050/Z/13/Z (J. J.) and No. 207510/Z/17/Z (R. E. G.). For the purpose of open access, the authors have applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. Support is also acknowledged from Medical Research Council Grant No. MR/P000479/1 (J. J.), ERC Consolidator Grant No. 682754 (E. L.), Established Career Fellowship EP/M017982/1 from the Engineering and Physical Sciences Research Council, and the Schlumberger Chair Fund (R. E. G.).

[1] G. J. Pazour, N. Agrin, J. Leszyk, and G. B. Witman, Proteomic analysis of a eukaryotic cilium, *J. Cell Biol.* **170**, 103 (2005).

[2] C. Ainsworth, Tails of the unexpected, *Nature* (London) **448**, 638 (2007).

[3] D. R. Mitchell, The evolution of eukaryotic cilia and flagella as motile and sensory organelles, in *Eukaryotic Membranes and Cytoskeleton*, Advances in Experimental Medicine and Biology Vol. 607 (Springer, New York, NY, 2007), p. 130.

[4] R. E. Goldstein, Green algae as model organisms for biological fluid dynamics, *Annu. Rev. Fluid Mech.* **47**, 343 (2015).

[5] T. Nakamura and H. Hamada, Left-right patterning: Conserved and divergent mechanisms, *Development* (Cambridge, U.K.) **139**, 3257 (2012).

[6] O. Thouvenin, L. Keiser, Y. Cantaut-Belarif, M. Carbo-Tano, F. Verweij, N. Jurisch-Yaksi, P. L. Bardet, G. van Niel, and F. Gallaire, Origin and role of the cerebrospinal fluid bidirectional flow in the central canal, *eLife* **9**, e47699 (2020).

[7] R. R. Ferreira, A. Vilfan, F. Jülicher, W. Supatto, and J. Vermot, Physical limits of flow sensing in the left-right organizer, *eLife* **6**, e25078 (2017).

[8] L. J. Fauci and R. Dillon, Biofluidmechanics of reproduction, *Annu. Rev. Fluid Mech.* **38**, 371 (2006).

[9] R. Faubel, C. Westendorf, E. Bodenschatz, and G. Eichele, Cilia-based flow network in the brain ventricles, *Science* **353**, 176 (2016).

[10] H. Soares, B. Carmona, S. Nolasco, and L. V. Melo, Polarity in ciliate models: From cilia to cell architecture, *Front. Cell Dev. Biol.* **7**, 240 (2019).

[11] E. R. Brooks and J. B. Wallingford, Multiciliated cells: A review, *Curr. Biol.* **24**, R973 (2014).

[12] Y. Liu, N. Pathak, A. Kramer-Zucker, and I. A. Drummond, Notch signalling controls the differentiation of transporting epithelia and multiciliated cells in the zebrafish pronephros, *Development* (Cambridge, U.K.) **134**, 1111 (2007).

[13] A. Vasilyev, Y. Liu, S. Mudumana, S. Mangos, P.-Y. Lam, A. Majumdar, J. Zhao, K.-L. Poon, I. Kondrychyn, V. Korzh, and I. A. Drummond, Collective cell migration drives morphogenesis of the kidney nephron, *PLoS Biol.* **7**, e1000009 (2009).

[14] G. R. Ramirez-San Juan, A. J. T. M. Mathijssen, M. He, L. Jan, W. Marshall, and M. Prakash, Multi-scale spatial heterogeneity enhances particle clearance in airway ciliary arrays, *Nat. Phys.* **16**, 958 (2020).

[15] M.-K. Khelloufi, E. Loiseau, M. Jaeger, N. Molinari, P. Chanez, D. Gras, and A. Viallat, Spatiotemporal organization of cilia drives multiscale mucus swirls in model human bronchial epithelium, *Sci. Rep.* **8**, 2447 (2018).

[16] E. Loiseau, S. Gsell, A. Nommick, C. Jomard, D. Gras, P. Chanez, U. D'Ortona, L. Kodjabachian, J. Favier, and A. Viallat, Active mucus-cilia hydrodynamic coupling drives self-organization of human bronchial epithelium, *Nat. Phys.* **16**, 1158 (2020).

[17] P. Walentek and I. K. Quigley, What we can learn from a tadpole about ciliopathies and airway diseases: Using systems biology in *Xenopus* to study cilia and mucociliary epithelia, *Genesis* **55**, e23001 (2017).

[18] V. C. Twitty, Experimental studies on the ciliary action of amphibian embryos, *J. Exp. Zool.* **50**, 319 (1928).

[19] G. A. Deblandre, D. A. Wettstein, N. Koyano-Nakagawa, and C. Kintner, A two-step mechanism generates the spacing pattern of the ciliated cells in the skin of *Xenopus* embryos, *Development* (Cambridge, U.K.) **126**, 4715 (1999).

[20] M. E. Werner and B. J. Mitchell, Using *Xenopus* skin to study cilia development and function, *Methods Enzymol.* **525**, 191 (2013).

[21] C. Brennen and H. Winet, Fluid mechanics of propulsion by cilia and flagella, *Annu. Rev. Fluid Mech.* **9**, 339 (1977).

[22] K. Karimi, J. D. Fortriede, V. S. Lotay, K. A. Burns, D. Z. Wang et al., *Xenbase: A genomic, epigenomic and*

<footnote>* fb448@cam.ac.uk</footnote>

<footnote>† jerome.jullien@inserm.fr</footnote>

<footnote>e.lauga@damtp.cam.ac.uk</footnote>

<footnote>§ R.E.Goldstein@damtp.cam.ac.uk</footnote>
transcriptomic model organism database, Nucleic Acids Res. **46**, D861 (2018).

[23] J. Faber and P. D. Nieuwkoop, *Normal Table of Xenopus laevis* (Daudin) (Garland Publishing Inc., New York, 1994).

[24] S. Nagata, Isolation, characterization, and extra-embryonic secretion of the *Xenopus laevis* embryonic epidermal lectin, XEEL, *Glycobiology* **15**, 281 (2005).

[25] E. Dubaissi and N. Papalopulu, Embryonic frog epidermis: A model for the study of cell-cell interactions in the development of mucociliary disease, *Disease Mod. Mech.* **4**, 179 (2011).

[26] E. Dubaissi, K. Rousseau, R. Lea, X. Soto, S. Nardeosingh, A. Schweickert, E. Amaya, D. J. Thornton, and N. Papalopulu, A secretory cell type develops alongside multiciliated cells, ionocytes and goblet cells, and provides a protective, anti-infective function in the frog embryonic mucociliary epidermis, *Development* (Cambridge, U.K.) **141**, 1514 (2014).

[27] P. Walentek, S. Bogusch, T. Thumberger, P. Vick, E. Dubaissi, T. Beyer, M. Blum, and A. Schweickert, A novel serotonin-secreting cell type regulates ciliary motility in the mucociliary epidermis of *Xenopus* tadpoles, *Development* (Cambridge, U.K.) **141**, 1526 (2014).

[28] E. Hörmanseder, A. Simeone, G. E. Allen, C. R. Bradshaw, M. Figlmüller, J. Gurdon, and J. Jullien, $H_3K_4$ methylation-dependent memory of somatic cell identity inhibits reprogramming and development of nuclear transfer embryos, *Cell Stem Cell* **21**, 135 (2017).

[29] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.127.198102 for further details and results, which includes Refs. [30–33].

[30] www.enasco.com.

[31] P. R. Bevington and D. K. Robinson, *Data Reduction and Error Analysis for the Physical Sciences*, 3rd ed. (McGraw-Hill, New York, 1993).

[32] T. J. Park, B. J. Mitchell, P. B. Abitua, C. Kintner, and J. B. Wallingford, Dishevelled controls apical docking and planar

polarization of basal bodies in ciliated epithelial cells, *Nat. Genet.* **40**, 871 (2008).

[33] F. Boselli, D. Obrist, and L. Kleiser, A multilayer method of fundamental solutions for Stokes flow problems, *J. Comput. Phys.* **231**, 6139 (2012).

[34] D. R. Brumley, K. Y. Wan, M. Polin, and R. E. Goldstein, Flagellar synchronization through direct hydrodynamic interactions, *eLife* **3**, e02750 (2014).

[35] J. R. Blake, Note on the image system for a stokeslet in a no-slip boundary, *Math. Proc. Cambridge Philos. Soc.* **70**, 303 (1971).

[36] J. R. Blake and A. T. Chwang, Fundamental singularities of viscous flow. Part I: The image systems in the vicinity of a stationary no-slip boundary, *J. Eng. Math.* **8**, 23 (1974).

[37] J. R. Blake, Infinite models for ciliary propulsion, *J. Fluid Mech.* **49**, 209 (1971).

[38] T. J. Pedley, D. R. Brumley, and R. E. Goldstein, Squirmers with swirl: A model for Volvox swimming, *J. Fluid Mech.* **798**, 165 (2016).

[39] N. Osterman and A. Vilfan, Finding the ciliary beating pattern with optimal efficiency, *Proc. Natl. Acad. Sci. U.S.A.* **108**, 15727 (2011).

[40] J. Gray and G. J. Hancock, The propulsion of sea-urchin spermatozoa, *J. Exp. Biol.* **32**, 802 (1955).

[41] N. Pellicciotta, E. Hamilton, J. Kotar, M. Faucourt, N. Delgehyr, N. Spassky, and P. Cicuta, Entrainment of mammalian motile cilia in the brain with hydrodynamic forces, *Proc. Natl. Acad. Sci. U.S.A.* **117**, 8315 (2020).

[42] H. Jeffreys, XXX. Some problems of evaporation, *Philos. Mag.* **35**, 270 (1918).

[43] H. C. Berg and E. M. Purcell, Physics of chemoreception, *Biophys. J.* **20**, 193 (1977).

[44] R. E. Goldstein, Coffee stains, cell receptors, and time crystals: Lessons from the old literature, *Phys. Today* **71**, No. 9, 32 (2018).