Hydrodynamic model of directional ciliary-beat organization in human airways

In the lung, the airway surface is protected by mucus, whose transport and evacuation is ensured through active ciliary beating. The mechanisms governing the long-range directional organization of ciliary beats, required for effective mucus transport, are much debated. Here, we experimentally show on human bronchial epithelium reconstituted in-vitro that the dynamics of ciliary-beat orientation is closely connected to hydrodynamic effects. To examine the fundamental mechanisms of this self-organization process, we build a two-dimensional model in which the hydrodynamic coupling between cilia is provided by a streamwise-alignment rule governing the local orientation of the ciliary forcing. The model reproduces the emergence of the mucus swirls observed in the experiments. The predicted swirl sizes, which scale with the ciliary density and mucus viscosity, are in agreement with in-vitro measurements. A transition from the swirly regime to a long-range unidirectional mucus flow allowing effective clearance occurs at high ciliary density and high mucus viscosity. In the latter case, the mucus flow tends to spontaneously align with the bronchus axis due to hydrodynamic effects.


Results
Ciliary beat reorientation in human bronchial epithelium cultures. The experimental system consists of well-differentiated human airway cultures at the air-liquid interface (ALI). These cultures exhibit the pseudo-stratified structure typical of bronchial epithelium, with basal cells, goblet cells that produce mucus, and ciliated cells with beating cilia 7 . As described by Khelloufi et al. 7 , the ciliogenesis is accompanied by the emergence of swirly mucus flows. The ciliary density on mature epithelium reaches around 70% of the epithelial surface. At high ciliary densities, the mucus is transported over macroscopic distances (see Movie S1 for a typical mucus flow) thanks to the continuous beating of cilia. The direction of the mucus flow and ciliary beats could be quantified as illustrated in Fig. 1 and Movies S1 and S2. Locally, as shown in Fig. 1(a), the cilia do not necessarily beat along the mucus flow direction but they present an average direction at large scale. Yet, we observed in several cultures that the ciliary beating directions tend to align along the mucus flow direction ( Fig. 1(a,b) and Movies S1 and S2) in a slow process of a maximum 35-40° angle reorientation within 24 hours (see angle distributions in Fig. 1(d,e)). In contrast, when mucus is washed out from the surface of the cell culture and replaced with culture medium, we observe within three days a loss of the coordination of beating directions ( Fig. 1(c)). In mature epithelia where mucus was regularly washed out, ciliary-beat directions remain stable overtime and do not exhibit particular temporal reorientation (see Fig. S1). Moreover, it should be mentioned that at this stage of the development, the tissue is jammed, namely the cells are not motile anymore 21 , and the turnover of epithelial cells is very slow (average half-life larger than 6 months for ciliated cells, see Rawlins and Hogan 22 ). The dynamics described in Fig. 1 is thus disconnected from these effects. In the disorganized state depicted in Fig. 1(c), only small patches of cilia attached to neighboring ciliated cells present the same beating direction. The characteristic diameter of these patches is around 10 to 20 μm, corresponding to groups of 5 to 15 cells. These patches can be seen as unit elements of beating direction, likely constituted by ciliated cells having a common cell polarity. There is no transport of fluid at large scale anymore. These observations disclose the existence of a strong mucus-cilia interplay and that mucus flow triggers an active response of cilia driving their reorientation.
Hydrodynamic model of ciliary-beat organization. We assume that the hydrodynamic forces generated by the mucus flow drives the observed long-range self-organization of the directions of ciliary beats, and propose the following two-dimensional model to examine this hypothesis. It is schematized in Fig. 2(a).
The main principles of the model are the following. Locally, ciliary beats exert a force on the mucus. Considering a simple modeling of this force, the overall mucus flow resulting from the joint forcing of all ciliated cells is computed by solving the Navier-Stokes equations. The model is based on a streamwise-orientation hypothesis: ciliary beats, and the associated local forces, reorient in response to the computed mucus flow and gradually align with the local flow direction. Such system, initialized with random orientations of the ciliary beats, thus evolves towards a self-organized solution. Distant cilia interact with each other through their individual contribution to the overall forcing field driving the global mucus flow.
In practice, the plane epithelium surface is discretized using hexagonal unit elements, which represent cells or groups of cells (see Fig. 2(a)). Their side length is denoted by D. Two types of unit elements are considered, namely multi-ciliated cells, or groups of neighboring multi-ciliated cells that exhibit a joint ciliary-beat direction, and passive elements that represent club and goblet cells and, possibly, dysfunctional ciliated cells. We chose a random spatial distribution of these elements on the epithelium, in agreement with observations 11 .
The surface fluid consists of two layers, a low-viscosity periciliary layer (PCL), in which cilia beat, and a high-viscosity mucus layer 1 . The PCL thickness is about the cilia length. The PCL is a shear region; the fluid velocity, in the plane parallel to the epithelium, is zero on the epithelial surface and is equal to the mucus velocity at the PCL/mucus interface 23,24 (see Fig. 2(a)). The mucus layer is modeled as a transport zone. Its velocity is mainly parallel to the epithelial plane and is uniform along the z axis, as schematized in Fig. 2(a). Overall, the dynamics of the mucus layer is modeled as a two-dimensional flow developing over a slip layer. This flow is driven by local propelling forces arising from ciliary beating. The forcing is spatially averaged over one or several neighboring multi-ciliated cells (unit element), and temporally over a time scale larger than the ciliary beating period.
Mucus is modeled as a Newtonian fluid, characterized by its density m ρ and dynamic viscosity m µ . The flow is governed by the two-dimensional incompressible Navier-Stokes equations. The flow momentum reads where t Δ is the model time step, Ω is a fixed angular velocity, θ 0 is an angle threshold, I n is the noise intensity characterizing rotational diffusion and η is a normalized Gaussian noise. The angle 0 θ , which is kept to very small values in this study, is introduced in order to allow the emergence of steady solutions (θ θ In presence of a transported mucus above the cilia, we observe a dynamic reorientation of the ciliary beat directions. Over 24 h the two peaks of the distribution get more narrow and are shifted toward the direction of the mucus flow measured from Movie S1. www.nature.com/scientificreports www.nature.com/scientificreports/ absence of noise. This alignment rule provides the hydrodynamic interaction mechanism between ciliated elements. It is introduced to drive the coupled cilium-flow system to a steady solution where ciliary forcing and mucus flow are locally aligned, namely when Eq. (2a) is satisfied. More details on the implemented alignment rule are reported in Methods.
Solutions of Eq. (1) are found using numerical simulations (see Methods). Periodic conditions are set at the boundaries of the computational domain. The ciliated elements are randomly placed and their total number is set according to the prescribed ciliary density S S c φ = , where S c designates the ciliated area and S is the total area. All computations are initialized with random orientations of the ciliated elements. Numerical and physical parameters. As noise is not expected to be a dominant mechanism in the real system, the focus is first placed on the model behavior when = I 0 n . In this case, Eq. (2) can evolve to a steady solution where all ciliary elements are close to aligned with the local flow, as defined by Eq. (2b). These steady organizations, once reached, do not depend on the angular velocity Ω. The latter parameter can thus be considered as a numerical parameter. The threshold angle 0 θ is kept to a small value, so that its influence on model solutions can be neglected.
The PCL friction coefficient κ is assumed to relate to PCL properties. From dimensional analysis, it reads κ µ δ ∼ / p p 2 , where p µ and δ p are the dynamic viscosity and thickness of the periciliary layer. In this model, mucus transport results from the competition between the ciliary forcing and the friction force. Accordingly, a typical mucus flow velocity can be derived from f c and κ, namely This indicates that the mucus velocity is assumed to mostly depend on the ciliary forcing and PCL properties.
Overall, the model mainly relies on six physical parameters, namely ρ m , m µ , φ, D, f c and κ. Following the Buckingham theorem, the number of parameters can be reduced to three non-dimensional and independent physical parameters, (i) the ciliary density φ, (ii) the Reynolds number which represents the ratio between inertial and viscous forces and (iii) the non-dimensional interaction length λ, defined as The λ parameter represents the typical range of influence of a ciliated element, resulting from the competition between viscous momentum diffusion over the epithelial plane and dissipation due to the friction force f v . Indeed, if λ is large (i.e. high mucus viscosity), the momentum transferred to the fluid through ciliary beating is rapidly diffused over a large fluid region; the limit case is the solid behavior, where the localized ciliary forcing is instantaneously transferred all over the material. In contrast, the range of influence decreases when λ is small, since the flow momentum tends to be damped faster than it is transported over the epithelial plane. Ciliated cells interact with each other only if their separation distance is low enough compared to their individual range of influence. Therefore, λ is referred to as hydrodynamic interaction length. Considering the above-mentioned connections between κ and the PCL properties, the interaction length is closely related to the mucus/PCL viscosity ratio, namely The role of λ is illustrated in Fig. 2(b,c), which show the typical flow in the vicinity of two ciliated elements, denoted by i and j. The non-dimensional flow vorticity is defined as U In the figure, high values of the vorticity (either positive or negative) emphasize fluid regions that are sheared along the epithelial plane, as an effect of the ciliary forcing. In this illustrative example, the ciliary force directions are fixed, i.e. the alignment rule Eq. (2) has been ignored. For low values of λ ( Fig. 2(b)), shear regions remain localized in the vicinity of the ciliated cells, and streamlines crossing the ciliated elements are mostly parallel to ciliary forces. Therefore, 0 i j θ θ Δ ≈ Δ ≈ , and Eq. (2b) is satisfied even in the absence of realignment process. For larger values of λ (Fig. 2(c)), large shear regions are observed, even though of lower magnitude than in the previous case, and flows induced by both ciliated cells tend to interact with each other. Consequently, streamlines are not parallel to ciliary forces, as emphasized by the locally-averaged flow velocities indicated in the figure. In this case, a trigger of the alignment rule Eq. (2) induces a re-alignment of the cells. This illustrates that λ is expected to control hydrodynamic interactions between ciliated elements.
To determine a relevant range of φ and λ for the simulations, one has to estimate their physiological range. During the ciliogenesis, the cilia density increases from 0 up to a final value of 70% which corresponds to the cilia density in the bronchus 25,26 . The relevant values of λ can be evaluated from Eq. (4). The size D for the unit elements is estimated from the size of the ciliated patches exhibiting a joint beating direction, which, as seen in Fig. 1, is of the order of ≈ × − D 2 10 m 5 . The thickness of the PCL has been measured by Button et al. 27 and is of the order of δ ≈ − 10 m p 5 . It is accepted that the PCL viscosity is close to that of water 27 (5), are then expected in the range [1 5, 150] .
. Finally, it should be mentioned that the flow behavior is expected to be governed by viscous effects, i.e.  Re 1.
Three different mucus flow regimes: localized, swirly, long-range aligned. A square domain composed of approximately 10 4 elements is first considered. Rotational diffusion is neglected ( = I 0 n ), and the Reynolds number is set to 0.1. The flow solution is computed over ranges of λ and φ, namely [1,150] . [0 1, 0 7]. The system dynamics is computed until a steady state is reached, when the ciliary-beat organization and the related mucus flow remain constant in time (see Movies S3, S4 and S5 for examples of transient dynamics). In this state, the ciliary forces and the flow are aligned locally. Three regimes associated with distinct flow patterns and ciliary-beat organizations are identified, as illustrated in Fig. 3.
When the ciliary density φ and the interaction length λ are small, typically 0 1 φ = . and 1 λ = , a poorly organized regime is observed (Fig. 3(a)). Only local ciliary-beat alignments are observed in regions of high local ciliary densities. No large-scale flow pattern emerges. The flow length scale is of the order of the typical cell size, as indicated by the localized and short-range shear regions close to ciliated elements.
When φ and λ are larger, typically φ = . 0 4 and 2 λ = (Fig. 3(b)), a remarkable swirly pattern, associated with local mucus recirculations, emerges. The ciliary beat orientation exhibits circular alignments over length scales much larger than the typical cell size, emphasizing the emergence of large-scale collective behavior. These large vortices are associated with large-scale vorticity regions. Both clockwise and counter-clockwise flow circulations are observed, as shown by the negative (blue) and positive (yellow) values of the vorticity.
When ciliary density and interaction length are further increased (Fig. 3(c)), typically φ = . 0 55 and 4 λ = , ciliary beat directions tend to fully align over the whole domain. In this case, the flow is unidirectional and almost uniform. As a result, the flow is hardly sheared, as indicated by the low magnitude of the vorticity. The fully aligned regime is the only regime allowing an overall fluid transport. In this case, the average flow velocity can be computed analytically from the momentum balance, expressed as c where the symbol U , in this equation and in the following, denotes the space-averaging operator. The mucus velocity is thus expected to linearly increase as a function of the ciliary density. In contrast, the mucus velocity appears to be independent of the mucus viscosity. This is an expected feature of the present model, which involves periodic boundary conditions, thus avoiding any energy dissipation due to boundary layers. Such dissipation is however expected to occur in real systems involving solid boundaries. In order to quantitatively characterize the three regimes, two physical quantities are defined. The first one is the polarization P, which is a spatial averaging of the unitary velocity vectors and quantifies the orientational order of the overall flow. 1 ≈ P is a signature of a unidirectional flow. The second quantity, Λ, is defined as, are auto-correlation functions of the vorticity along the x and y directions respectively.
Λ is a non-dimensional integral length. It determines the typical length scale of flow structures. In Eq. (8), it is normalized by / m µ κ, the dimensional interaction length. Therefore, Λ = 1 indicates that the flow length scale corresponds to the range of influence of one ciliated element, i.e. no long-range dynamics occurs. In contrast, Λ > 1 indicates that large-scale vorticity regions are present, as in Fig. 3(b).
The evolutions of P and Λ in the ( , ) φ λ plane are plotted in Fig. 4(a,b), for λ ranging from 1 to 5, i.e. the range where regime transitions occur. From these variations, we identify three regions in the φ λ ( , ) space: the low-λ/ low-φ region corresponds to a poorly organized regime, the low-λ/high-φ region corresponds to a swirly regime, and the high-λ/high-φ region corresponds to a fully aligned regime. We establish the flow regime phase diagram (Fig. 4(c)) by combining the two parameters P and Λ as following: The emerging regime may vary from one simulation to the other due to the random initialization of the computations. In Fig. 4(c), each point of the phase diagram represents the most frequent regime emerging over a set of ten simulations. The greyscale used in the figure illustrates the sensitivity of the regimes to initial conditions. As expected, flow regimes are sensitive to initial conditions close to the transition regions.  10); at each point, the most frequent regime over a set of ten randomly initialized runs is indicated by a triangle (poorly aligned), a circle (swirly) or a square (fully aligned) symbol. Symbols are colored using the occurrence frequency, ranging from 1/3 (white) to 1 (black). Dashed lines delimitate the three identified regions.

Scientific RepoRtS |
(2020) 10:8405 | https://doi.org/10.1038/s41598-020-64695-w www.nature.com/scientificreports www.nature.com/scientificreports/ The observed regimes remain robust in the presence of noise. In particular, the swirly regime is almost unaffected by the noise source even when the noise intensity I n is larger than the alignment velocity Ω (see Eqs. 2), although ciliary disorganization can be achieved for I n Ω  . Details on the effect of noise on the swirly regime are provided in the Supplementary information.

Swirl size scales with ciliary density and interaction length. The size of the vortices in the swirly
regime is determined using the auto-correlation of the flow vorticity, defined by Eq. (9). Additional series of simulations are performed in the region where the swirly regime typically emerges, namely [0 1, 0 7] φ ∈ .
. and [1, 2 5] λ ∈ . . For each point of the φ λ ( , ) space, 25 randomly-initialized simulations are carried out. Cases where the number of occurrences of the swirly regime is too small (lower than 10) are discarded. Simulations eventually converging to the swirly regime are selected for the computation of the swirl size. For each case, the average auto-correlation function τ R ( ) is computed, and the swirl size is defined as the minimum length scale τ satisfying τ = R ( ) 0. The swirl size is normalized by D and denoted by Λ ′ . The evolution of Λ ′ as a function of the ciliary density and interaction length is depicted in Fig. 5(a). Swirl sizes determined following this procedure are in agreement with qualitative visualizations of the flow. For λ = 2 and φ = .
0 4, the measured swirl size Λ ≈ ′ 20 is consistent with the swirls depicted in Fig. 3(b). At low ciliary densities, the globally-disorganized solution may include localized swirls, as observed in Fig. 3(a) for λ = 1 and φ = . 0 1. The size of these small vortices is consistent with the value 5 Λ ≈ ′ plotted in Fig. 5(a). Overall, the swirl size increases as a function of φ and λ. This qualitative trend can be understood based on simple physical arguments. First, it is expected that an increase of the hydrodynamic length scale λ results in a growth of the emerging self-organized patterns, a process that eventually leads to the overall alignment of the system when the interaction length is larger than the typical distance between ciliated elements. In addition, this distance naturally decreases as a function of the ciliary density, promoting the propagation of the alignment rule and the emergence of large-scale coherent structures.
The swirl sizes measured in the present in vitro experiments, as well as those reported by Khelloufi et al. 7 , have been included in 5(a) for comparison purpose. Only unconfined results, when the swirls are not impacted by the walls of the culture chamber, are reported here. In the figure, the swirl size has been normalized by D 2 10 m 5 = × − , namely the typical size of ciliary patches observed in the experiments (see Fig. 1). Both experimental data sets exhibit a consistent trend that is in agreement with the evolution predicted by the present model, in particular for the case λ = 2.
The evolution of Λ ′ in Fig. 5(a) follows a trend that is close to linear as a function of φ, with an increasing slope as a function of λ. Assuming a trend λ φ Λ = α ′ , a manual fit suggests 1 5 α ≈ . . The evolution / 1 5 λ Λ ′ . as a function of φ is plotted in Fig. 5(b). A striking collapse of the experimental and numerical data is noted, emphasizing the linear evolution the swirl size as a function of the ciliary density. However, the detailed physical mechanisms involved in this scaling remain to be understood.

Ciliary beats align with the axis of a virtual airway.
To get deeper insight in the role of the domain geometry on the flow direction in the fully aligned regime, we performed a series of 96 randomly initialized computations in a rectangular box with different aspect ratios. The aspect ratio is varied by increasing the domain size in the x direction (L x ); the domain size in the y direction (L y ) is kept constant. We use a box of the order of 2500 unit elements, with 5 λ = and φ = .
0 4 for which a fully aligned flow regime is predicted. For a square domain, no preferential flow direction is observed from one computation to the other (Fig. 6(a)). When the aspect ratio of the domain is equal to 2 ( Fig. 6(b)), a preferred flow direction emerges along the longest side of the domain. Upon Figure 5. Vortex size scales with φ and λ. The swirl size Λ ′ is plotted as a function of the ciliary density for different values of the interaction length in (a), and compared to experimental measurements issued from the present study and from the work of Khelloufi et al. 7 . In (b), all the data collapse on a linear curve when Λ ′ is normalized by 1 5 λ . , illustrating the scaling of Λ ′ with λ and φ. Experimental data are normalized using λ = 2.
This alignment mechanism relates to an anisotropic confinement effect. Here, the term confinement is employed as the size of the periodic domain (L D / 50 y = ) has an order of magnitude comparable to the range of influence of a ciliated element, determined by λ. Therefore, the periodicity of the domain has an influence on the ciliary organization. When the periodicity is anisotropic, the isotropy of the solution is lost and a preferential direction spontaneously emerges, even though this direction can hardly be predicted a priori. The present simulations indicate that the preferred direction is aligned with the longest side of the rectangular domain.
The rectangular geometry is particularly interesting as it is a Cartesian representation of the surface of a circular tube aligned with the x axis. The y periodicity of the rectangle relates to the azimuthal periodicity of the tube, and thus depends on its diameter. In the x direction, the periodic boundary condition ensures that the flow is not constrained in this direction. The possible connections with a natural airway are further discussed hereafter.

Discussion
We proposed a hydrodynamic model to investigate the emergence of orientational order on a bronchial epithelium. It is based on the description of the two-dimensional mucus flow over an array of self-oriented momentum sources. The orientation of these driving forces is controlled by a streamwise-alignment rule, providing a hydrodynamic coupling between the flow and the ciliary forcing. The model is governed by two independent parameters, the ciliary density and the interaction length, related to the friction over the epithelium and depending on the mucus/PCL viscosity ratio.
The first important result of the model is that it predicts the three main flow regimes experimentally observed: (i) the local flow regime (Fig. 1(c)), (ii) the swirly regime 7 and (iii) the aligned flow ( Fig. 1(a,b)). Experimentally, ciliary disorganization is observed when the mucus is washed using the culture medium, whose mechanical properties are expected to be close to that of water (see Fig. 1(c)). This configuration can be addressed with the present model using m p µ µ = , leading to 0 5 λ ≈ . . Even at high ciliary densities, only small mucus recirculations are predicted by the model in these conditions (see Fig. S5, which shows the self-organized pattern for λ = . 0 5 and φ = . . The fully aligned regime allows global mucus transport at the average velocity given by Eq. (6), depending on the ciliary forcing f c , the friction coefficient κ and the ciliary density φ. The typical force that can be developed by a single cilium on bronchial epithelium has been measured by Hill et al. 29 and is equal to ≈ × − F 6 10 11 N. Assuming that this force is only converted into hydrodynamic forcing during the power stoke of the ciliary beating cycle, the time-averaged propulsion force exerted by one cilium can be roughly estimated as F/2. Considering that a ciliated cell of radius = . × − R 2 5 10 m 6 carries around = n 250 c cilia 30 , the force density, equivalent to f c , can be approximated by , the average mucus velocity on a mature epithelium (φ = . 0 7) is | | ≈ . × − U 2 7 10 5 m.s −1 , a value that is close to the clearance velocity experimentally measured by Foster et al. 31 approximately equal to × − 4 10 5 m.s −1 . Interestingly, the model predicts a beneficial effect of the mucus viscosity on the clearance, by promoting long-range interactions between ciliated cells. Once the aligned regime is achieved, the overall transport is not www.nature.com/scientificreports www.nature.com/scientificreports/ affected by the mucus viscosity anymore, since the predicted flow velocity only depends on the ciliary forcing, ciliary density and friction coefficient (see Eq. (6)). In contrast, very large mucus viscosities are generally expected to result in poor mucus clearance in a physiological context. This effect might be due to additional local viscous dissipation mechanisms, related to the presence of obstacles or boundaries for instance, that have not been included in the present model. Moreover, non-hydrodynamic effects may emerge in the biological system when the mucus velocity decreases, such as bacterial proliferation within the mucus in the airways, local adhesion, dehydration of the mucus leading to an increase in mucin concentration and to cilia collapse, local increase in the amount of mucus leading to mucus plugs, etc. These antagonistic effects may suggest the existence of an optimal viscosity allowing long-range hydrodynamic interaction while ensuring a minimum mucus velocity, a concept that might be useful for future investigations on mucus clearance.
The role of the ciliary density on long-range mucus flow is also quantitatively described by the model. It could therefore be of interest for diseases affecting ciliary motility, as primary dyskinesia. Such a model can open the way to the prediction of the minimal density and spatial distribution of active cilia required to ensure a long-range mucus transport, i.e. the number of ciliated cells that has to be repaired on the airway.
A remarkable result of the model concerns the flow in rectangular domains with periodic boundary conditions. As stated above, this geometry can represent the 2D projection of a circular tube, similar to a natural airway. This comparison with a virtual human bronchus can be further analyzed by considering typical geometrical quantities. Indeed, an effect of the azimuthal periodicity on the directional organization is expected to occur only if the periodicity length and the interaction length have comparable orders of magnitude. Based on the morphometric model of Weibel et al. 32 , the typical cross-section perimeter of human airways is expected to range from L 10 4 ≈ − m to ≈ − L 10 2 m 33 . By using ≈ × − D 2 10 5 m, the non-dimensional periodicity length is therefore expected to range from ≈ L D / 5 to ≈ L D / 500. This range is of the same order of magnitude as the expected range of interaction lengths in natural airway, namely [1 5, 150] λ ∈ .
as discussed above, and interaction between both length scales is thus possible. The results of the model therefore suggest that, in addition to the biologically driven phenomena occurring during airway development, the distal-proximal alignment of ciliary beats is favored by hydrodynamic interactions in a tube geometry.

Methods
Experimental methods. In-vitro bronchial epithelium. Mature commercial cultures (MucilAir) of reconstituted human bronchial epithelium from primary cells, were bought from Epithelix (Switzerland). The tissue is cultured at the Air Liquid Interface in a transwell of 6 mm in diameter. The apical side of the tissue is at the air interface and the basal side is in contact with the culture medium (Epithelix MucilAir culture medium) via a porous membrane. The culture medium (700 μl) is replaced every two or three days. For experiments in absence of mucus, we removed the mucus by performing an apical washing. We added 200 μl of culture medium on the apical side for 10 minutes, followed by three rinses with the culture medium. After an apical washing, a thin layer of medium culture (~10-40 μm) remains on top of the cilia. Experiments were repeated on three different cultures.
Imaging. We imaged the cultures in bright field on an inverted Nikon Eclipse Ti microscope with a x20 objective and a Luminera infinity camera. The temperature is maintained at 37 °C and a humidified air flow of 5% CO2 is applied to maintain a physiological pH. To image the culture at the same place on different days, we used a motorized stage calibrated in the xy plane and saved the different positions.
Image processing. Beating directions of cilia were quantified using an in-house image processing routine applied on videomicroscopy acquisitions (40 fps). We first performed a standard deviation projection over 40 frames to reveal the trajectory of the tip of the cilia. Then we implemented in python the algorithm described in ref. 34 based on the structure tensor computation. numerical methods. The mucus flow, governed by Eq. (1), is simulated using a D Q 2 9 Bhatnagar-Gross-Krook lattice-Boltzmann method, which is second-order accurate in time and space 35 . The effect of the body force f n is ensured by the external forcing scheme proposed by Guo et al. 36 . The physical domain is discretized on a uniform and Cartesian numerical grid. The distance between two neighboring nodes, in the x or y direction, is denoted by Δn; it is set to Δn = D/5, where D denotes the side length of the hexagonal elements (see Fig. S2). An hexagonal element thus typically contains 65 fluid nodes. At each time step, the flow velocity is averaged in each ciliated elements, and the angle of the ciliary forcing is updated following Eq. (2). The parameter Ω, which drives the transient re-orientation of the ciliated elements but does not affect the final steady solution in the absence of noise, is set to Ω = U D / 0 , unless otherwise stated. It is recalled that U 0 designates the typical mucus velocity. A small threshold angle is employed, namely t 2 0004 0 θ = ΩΔ = . radians using the present numerical parameters, so that its influence on model solutions can be neglected. The Gaussian noise component η is updated at each time step using a Box-Muller transform. More details on the numerical methods are provided in the Supplementary information.