# A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part I, experimental analysis

Mathieu Bottier, Sylvain Blanchon, Gabriel Pelle, Emilie Bequignon, Daniel Isabey, André Coste, Estelle Escudier, James B. Grotberg, Jean-François Papon, Marcel Filoche, et al.

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

## To cite this version:

Mathieu Bottier, Sylvain Blanchon, Gabriel Pelle, Emilie Bequignon, Daniel Isabey, et al. A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part I, experimental analysis. PLoS Computational Biology, 2017, 13 (7), pp.e1005605. 10.1371/journal.pcbi.1005605. hal-01579454

HAL Id: hal-01579454

https://hal.sorbonne-universite.fr/hal-01579454

Submitted on 31 Aug 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

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

# A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part I, experimental analysis

Mathieu Bottier¹,²,³, Sylvain Blanchon¹,²,³,⁴, Gabriel Pelle¹,²,³, Emilie Bequignon¹,²,³,⁵, Daniel Isabey¹,²,³, André Coste¹,²,³,⁵,⁶, Estelle Escudier⁷,⁸,⁹, James B. Grotberg¹⁰, Jean-François Papon¹,³,¹¹, Marcel Filoche¹,²,³,¹², Bruno Louis¹,²,³*

1 Inserm U955, Equipe 13, Créteil, France, 2 Université Paris-Est, Faculté de médecine, Créteil, France, 3 CNRS, ERL 7240, Créteil, France, 4 CHU Toulouse, Hôpital des Enfants, Service de pneumologie-allergologie pédiatrique, Toulouse, France, 5 AP-HP, Hôpital H. Mondor-A. Chenevier, Service d'ORL et de chirurgie cervico-faciale, Créteil, France, 6 Hôpital intercommunal, Service d'ORL et de chirurgie cervico-faciale, Créteil, France, 7 Inserm, U933, Paris, France, 8 Université Pierre et Marie Curie, U933, Paris, France, 9 AP-HP, Hôpital Armand-Trousseau, Service de génétique et d'embryologie médicale, Paris, France, 10 Department of Biomedical Engineering, University of Michigan, Ann Arbor, MI, USA, 11 AP-HP, Hôpital Bicêtre, Service d'ORL et de chirurgie cervico-faciale, Le Kremlin-Bicêtre, France, 12 Physique de la Matière Condensée, Ecole Polytechnique, CNRS, Université Paris-Saclay, Palaiseau, France

* bruno.louis@inserm.fr

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

OPEN ACCESS

**Citation:** Bottier M, Blanchon S, Pelle G, Bequignon E, Isabey D, Coste A, et al. (2017) A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part I, experimental analysis. PLoS Comput Biol 13(7): e1005605. https://doi.org/10.1371/journal.pcbi.1005605

**Editor:** Mark Alber, University of California Riverside, UNITED STATES

**Received:** January 18, 2017

**Accepted:** May 30, 2017

**Published:** July 14, 2017

**Copyright:** © 2017 Bottier et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work has been supported by a grant from Fondation pour la Recherche Médicale (https://www.frm.org/, DI being the PI): FRM programme Bio-Ingénierie pour la Santé 2014, DBS 20140930771. The PhD Scholarship of MB has been provided by CNRS INSIS 2013 (http://www.cnrs.fr/). Some authors of this article are participants in BEAT-PCD (COST Action 1407). The

## Abstract

Mucociliary clearance is one of the major lines of defense of the respiratory system. The mucus layer coating the pulmonary airways is moved along and out of the lung by the activity of motile cilia, thus expelling the particles trapped in it. Here we compare *ex vivo* measurements of a Newtonian flow induced by cilia beating (using micro-beads as tracers) and a mathematical model of this fluid flow, presented in greater detail in a second companion article. Samples of nasal epithelial cells placed in water are recorded by high-speed video-microscopy and ciliary beat pattern is inferred. Automatic tracking of micro-beads, used as markers of the flow generated by cilia motion, enables us also to assess the velocity profile as a function of the distance above the cilia. This profile is shown to be essentially parabolic. The obtained experimental data are used to feed a 2D mathematical and numerical model of the coupling between cilia, fluid, and micro-bead motion. From the model and the experimental measurements, the shear stress exerted by the cilia is deduced. Finally, this shear stress, which can easily be measured in the clinical setting, is proposed as a new index for characterizing the efficiency of ciliary beating.

## Author summary

Mucociliary clearance is the first line of defense of the human pulmonary airways. Mucus transporting debris, particles, microorganisms and pollutants is carried away by the coordinated motion of cilia beating at the surface of the airway epithelium. We present here an experimental, mathematical and numerical study aiming at defining a global index for assessing the efficiency of this beating. We measure experimentally the ciliary beat
funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

frequency, ciliary beat amplitude, and metachronal wavelength on ciliated edges obtained from nasal brushing. Properties of fluid motion are simultaneously extracted from micro-bead tracking next to the ciliated edge. A mathematical and numerical model is developed to describe the fluid motion induced by the cilia tips considered as a moving wall. Experimental and numerical results show that the bead velocity is a parabolic function of the distance to the wall. It allows us to infer the shear stress exerted by the cilia on fluid from micro-bead tracking. This quantity is proposed as a universal index characterizing the beating efficiency, which can be extracted in the current clinical setting.

## Introduction

Mucociliary clearance is one of the major defense mechanisms of the respiratory airway system. The mucus layer coating the epithelial surface of the airways filters the inhaled air by trapping potentially harmful material (fungi, bacteria and other particles) [1–4]. This mucus layer is continuously carried away and out of the airways by the activity of motile cilia. Neighboring cilia beat in an organized manner with a small phase lag, their tips creating an undulating surface on top of the cilia layer which deforms in a wave-like fashion called the metachronal wave [5–7].

The beat pattern of an individual cilium displays a two-stroke effective-recovery motion [8]. During the effective stroke, cilia beat forwards and engage with the mucous layer, propelling it forward. In contrast, during the recovery stroke, they return to their initial position in the underlying periciliary fluid, minimizing thereby the drag on the mucus in the opposite direction. This asymmetry in the beat pattern is responsible for a net fluid flow in the direction of the effective stroke. In the airways, each mature ciliated cell may be covered with up to 200 cilia, with a surface density around 5–8 cilia/µm² [6, 9]. Each cilium, approximately 6 µm long and of diameter around 0.2 µm, beats 12 to 15 times per second, resulting in a velocity of the mucus layer of approximately 1 mm per minute [10].

Defects in mucociliary clearance may result in chronic airway inflammation and infections causing injury and structural changes to the airway epithelium, leading to a variety of diseases like bronchiectasis and chronic sinusitis. Two main reasons may lead to impaired mucociliary clearance. The first one is related to alterations of the mucus properties, as in cystic fibrosis where the mucus become too thick and sticky to be moved properly by the cilia. The second one is linked to a dysfunction in ciliary motion or ciliary coordination. Such dysfunctions may be either inherited as in Primary Ciliary Dyskinesia (PCD) or acquired. Diagnosis of dyskinesia is difficult. For instance, in the diagnosis of PCD, a combination of transmission electronic microscopy, nasal nitric NO, genetic testing and use of high-speed video-microscopy is recommended [11].

*In situ* observation of ciliary beating and mucociliary clearance is almost impossible in patients at present stage, and one currently lacks a reliable and general method for evaluating mucociliary clearance in the clinical field. In the past, integrated assessment of mucociliary clearance was achieved through techniques using saccharine [12], a drop of blue marker [13], or clearance of radioactive tracers [14–17]. More recently, micro-optical coherence tomography was also proposed [18]. However, due to their various requirements (patient cooperation for the saccharine test, endoscopic examination, inhalation of radiopharmaceutical, ...) these different methods are rarely used in the current clinical practice. Light microscopy observation of ciliated edge obtained by nasal or bronchial brushing is the most common method used to evaluate ciliary beating. This observation is often associated with High-Speed Video-
Microscopy (HSVM) analysis [19–23] that quantifies only cilia motion and ciliary beat pattern. However, none of these tests provides information about the global efficiency of ciliary beating regarding mucociliary clearance.

This present study aims at a deeper understanding of the relationship between ciliary beat pattern, measured by ciliary beating analysis, and the induced motion of the surrounding fluid in order to assess the global efficiency of ciliary beating. An original method of micro-bead tracking (MBT) is presented where micro-beads act as markers of the flow generated by the cilia [24]. A numerical model is then developed (presented in detail in a companion paper [25]) based on envelope modeling approach, allowing us to simulate the flow generated by the cilia as in the MBT experiment. We finally propose a new and global index for characterizing the efficiency of the ciliary beating. One very appealing aspect of this index is that it does not require any modification of the present clinical practice of data collection (nasal or bronchial brushing).

# Materials and methods

## Experimental analysis

**Patients.** Studies of ciliary beating and MBT were performed in 11 consecutive patients referred to our diagnostic center. All patients were investigated because of chronic upper and/or lower respiratory tract infections, i.e., bronchitis and/or bronchiectasis and sinusitis. Patients were selected only after ruling out genetic disorder as cystic fibrosis or PCD. Informed consent was obtained from all patients, and this study was approved by the local Ethics Committee (Comité de Protection des Personnes Ile-de-France XI).

**Digital high-speed video-microscopy.** Ciliated samples were obtained by brushing the middle part of the inferior turbinate with a 2 mm cytology brush. Cells were suspended in a survival medium and examined within three hours. 20 µL of 4.5 µm polystyrene micro-beads at the concentration of 0.125%w/v were added to 80 µL of the medium. All observations were performed within 20 min, at controlled room temperature (20–25°C). We used an inverted microscope (with a LD condenser 0.35 in H position, i.e., without any phase contrast or differential interference contrast) in brightfield conditions associated with a ×40 objective. 100 µL of medium containing beating ciliated edges and microbeads in suspension are comprised between a microscope slide and a cover slide. The medium is delimited in the horizontal plane by a circle of grease of diameter approximately equal to 1 cm (Fig 1). Recorded ciliated edges are located in the horizontal focal plane of the microscope (so that observed cilia are horizontal). The “local plane” of the cell cluster, which corresponds to the cell-medium interface, is therefore essentially perpendicular to the substrate in the focal plane. The observed region lies at the intersection of the horizontal focal plane and the locally vertical cluster plane.

Cilia were recorded with a digital camera at a rate of 358 frames per second (brightness 0%, gamma 220, shutter 720 s, gain 1.160 dB). Each movie was composed of 1,800 frames with a definition of 256 × 192 pixels, each individual pixel being (0.32 × 0.32) µm². All areas containing intact undisrupted ciliated epithelial edges larger than 50 µm, devoid of mucus, beating in the plane of the camera and close to micro-beads were recorded. As recommended in [26], isolated ciliated cells were excluded. Basically and except for the presence of microbeads, the high-speed video-microscopy procedure used here is equivalent (setup requirement and protocol) to the ones used by different groups for PCD diagnosis [19–21, 23].

**Micro-bead tracking method.** Locations of ciliated edges were determined using in-house software which reports positions of computer mouse clicks on paused video frames. Practically, the operator chooses 5 points delimiting 4 line segments. These segments define the location of the cilia wall (Fig 2). Micro-bead motions are tracked using an automatic image
<figure><img src="image_6.png" /><figcaption>**Fig 1. Scheme of experimental setup.**</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g001

processing method programmed in Python language. First, the PINK Image Processing Library is used to threshold the frames, thus separating the micro-beads from the background. Secondly, micro-beads are labeled and their trajectories are plotted. Instantaneous micro-bead velocities, as well as their distance to the ciliated edge, are measured and saved in a database.

**Analysis of ciliary beating.** Ciliary beat frequency (CBF) is determined from movies in two different ways. In both methods, the operator defines a line by placing five points on the ciliated edge. In the first method, the variation of the mean grey level of this line during the movie is calculated. A Fast Fourier Transform (FFT) is performed on this signal, and the frequency of highest amplitude in the FFT power spectrum is taken as CBF. In the second method, video-kymography is used by stacking all lines obtained from each frame of a movie on the same graph. The beating period is measured on this graph, from which the ciliary beat frequency is deduced.

To characterize the ciliary beating amplitude (CBA), we used a simplified quantitative analysis of ciliary beat pattern previously described by Papon et al. [20]. In short, one individual cilium that can be followed during a complete beating cycle (excluding cilia whose tip ran out of the focal plane) is selected. Video sequences are played back frame by frame in order to determine three points characterizing the complete cycle of the cilium. These three points correspond to the position of the base of the cilium (P0) and the positions of the tip before the active and recovery strokes (P1 and P2, respectively). These measurements are used to determine the cilium length and the distance traveled by the cilium tip (i.e., the ciliary beating amplitude).

The metachronal wavelength $\lambda$ is obtained by inferring the phase variation of the main beat frequency along the ciliated edge. To this end, 104 line segments perpendicular to the cilia wall are defined (see Fig 3, left), each segment being made of 24 pixels. In each segment, the measured mean grey level oscillates cyclically in time. The phase of this variation, $\theta$, is computed by a Fast Fourier Transform. In the case where this phase $\theta(x)$ depends linearly on the position x along the ciliated edge, a linear regression is used to determine the average slope $p$, hence the wavelength of the metachronal wave given by $\lambda = 2\pi/p$ (Fig 3, right).

In addition, cilia density is estimated for each ciliated edge by comparing the mean grey level of the movie background with the grey level of the cilia area. In the image, cilia appear
<figure><img src="image_7.png" /><figcaption>**Fig 2. Determining the location of the ciliated edges.** The operator places 5 points (blue crosses) which are the ends of 4 line segments (red dashed lines) defining the location of the cilia wall.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g002

darker than the background (Fig 4, left). Cilia density is therefore calculated by an in-house software as the percentage of pixels belonging to the ciliated edge darker than the mean grey level of the background. This procedure allows us to define a characteristic length of the distance between two cilia ($d_c$) and the fraction of area covered by the space without cilia ($f_c$).

Let us consider an homogeneous distribution of cilia. In this case, the cilia density is given by:

$$ 
\begin{aligned}
\rho_c &= \left( \frac{\text{number of pixel} \le \text{Mean Background}}{\text{number of pixel}} \right) \Big|_{\text{cilia zone}} \\
&= \frac{(e_c/\text{pixel size}) \times \text{number of pixel}}{[(e_c+d_c)/\text{pixel size}] \times \text{number of pixel}} = \frac{e_c}{e_c+d_c},
\end{aligned} 
\quad (1) $$

where $e_c$ is the cilium diameter (0.2 µm). As such, what is called “cilia density” here is a relative measure of the cilia density obtained by averaging over a sufficient large number of individual
<figure><img src="image_8.png" /><figcaption>**Fig 3.** Example of wavelength determination. (Left) 104 line segments are placed along the cilia wall, perpendicular to it. (Right) The mean grey level of each segment oscillates in time. The corresponding phase each segment is plotted as a function of the curvilinear abscissa of the center of each line. A linear regression is performed on the phase-abscissa relationship, the slope of this regression directly giving the metachronal wavelength (blue dashed line).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g003

pixels in the region of the ciliated edge. Let us note that $(1 - \rho_c)$ is in this case the fraction of area without cilia.

**Data selection.** Only the movies where the inputs required for ciliary beating simulation could be measured were retained (ciliary beat frequency, ciliary beating amplitude, metachronal wavelength and cilia density). These movies are the ones in which the cluster surface is essentially vertical so that the cilia edge located at the intersection of the cluster plane and the focal plane can clearly be observed. To remain within the assumptions of the model, we only selected micro-bead following straight trajectories close to a ciliated edge (i.e., not influenced by the external environment or other distant clusters). Finally, all analyses were performed by one single operator.

<figure><img src="image_9.png" /><figcaption>**Fig 4.** Measuring the cilia density. (Left) Two different regions are delimited on the real system, the green rectangle corresponds to the cilia region while the red ellipse corresponds to the background. (Right) Schematic representation of the same system: black rectangles represent the cilia, $e_c$ is the cilium diameter and $d_c$ is the typical distance between two cilia.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g004
<figure><img src="image_10.png" /><figcaption>Fig 5. Schematic representation of the stroke of an individual cilium and the envelope model. (Left) The trajectory of the stroke cycle is assumed to follow an elliptic motion. (Right) Representation of the envelope model covering the cilia layer and the propagation of the metachronal wave. (Inspired by Velez-Cordero et al. [28]).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g005

## Numerical model

We present here a brief overview of the two-dimensional model describing the system composed of the cilia, the surrounding fluid, and the micro-beads. The entire model is extensively detailed in [25]. Inspired by [27], we propose an approach in which the momentum transfer from discrete cilia to the fluid (induced by the ciliary beating) is modeled through an appropriate continuous boundary condition. The ciliated edge is chosen parallel to the x axis, the ciliary tip being located around y = 0 on average. Each cilium is assumed to undergo a periodic elliptic motion (see Fig 5, left). Taking the limit of a continuous cilia distribution, the cilia array is simplified as an undulating surface that covers the cilia layer, ignoring the details of the sub-layer dynamics, see Fig 5, right (inspired by Velez-Cordero et al. [28]).

The tip of a cilium located at the horizontal coordinate $\xi^*$ is assumed to follow a periodic elliptic trajectory centered in ($\xi^*$, 0) during each elementary beat (Fig 6). The '*' notation is used here to represent dimensional quantities, for consistency reason with [25]. At time $t^*$, the tip coordinates ($X_w^*$, $Y_w^*$) satisfy

$$ \begin{cases} X_w^* = \xi^* - a \cos(\omega t^*) \\ Y_w^* = \beta a \sin(\omega t^*) \end{cases} \qquad (2) $$

where $\beta$ is a function of the ellipse eccentricity, $2\beta a$ is its minor axis in the $y^*$ direction, and $2a$ its major axis in the $x^*$ direction. For $\beta > 0$, the tip orbits clockwise, while for $\beta < 0$, the tip orbits counterclockwise.

In this picture, the metachronal wave is materialized by prescribing the motion of the envelope. This envelope is then used as a boundary condition to compute the fluid velocity field dominated by viscous forces, inside a channel of height $h$ (see Fig 6). The fluid is assumed to be stagnant above this height. Indeed, the parameter $h$ essentially summarizes in our 2D model (see companion paper) the effect of the external environment of the clump, this environment imposing a limit to the spatial extension of the cilia-induced flow field. This fact
<figure><img src="image_11.png" /><figcaption>Fig 6. Schematic elliptic motion of an individual ciliary tip.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g006

explains why this parameter *h* cannot be estimated a priori and has to be measured as an external free parameter.

The trajectories of micro-beads in this velocity field are calculated by solving the equation of motion with Stokes drag. This theoretical model provides us with a prediction of the effective speed of micro-beads as a function of their average altitude above the ciliated edge. This velocity profile is found to follow essentially a parabolic profile as a function of the height above the cilia wall. Finally, the velocity of the micro-beads extrapolated at the cilia wall gives a direct estimate of the shear stress exerted on the fluid, this index being proposed as an index of the ciliary beating efficiency [25].

**Simulation of micro-beads effective velocity.** The bead inertia in the fluid flow is characterized by their Stokes number, $St_k$. For particles of mass *m* and radius *R* traveling through a fluid of viscosity *μ*, of typical velocity *U* and of typical length *L*, this number is defined as:

$$St_k = \frac{mU}{6\pi R\mu L} \qquad (3)$$

For spherical particles of radius *R* and density $\rho_b$ in an oscillatory flow of pulsation *ω*, this number also reads:

$$St_k = \frac{\frac{4}{3}\pi R^3 \rho_b \omega}{6\pi R \mu} = \frac{2 R^2 \rho_b \omega}{9 \mu} \qquad (4)$$

The micro-beads are about 4.5 µm diameter, and made out of polystyrene of density $\rho_b$ of order $1 \text{ g.cm}^{-3}$. At 10 Hz in water, the corresponding Stokes number is therefore about $10^{-4}$,
which means that the beads can be considered as massless tracers. Their velocity can be assumed to be permanently equal to the fluid velocity at the same location.

For each micro-bead entering the simulation window at $x = 0$ and a given altitude $y_0$, the effective speed is computed as:

$$V_{\text{eff}}(y_0) = \frac{L_w}{\tau(y_0)} = \frac{L_w}{\oint_{T_{y_0}} dt} = \frac{L_w}{\oint_{T_{y_0}} \frac{ds}{\| \vec{u}(s) \|}}, \quad (5)$$

where $L_w$ is the length of the observation window, $\tau(y_0)$ is the crossing time of the micro-bead entering at $(0, y_0)$, and $T_{y_0}$ is the trajectory followed by this micro-bead. During each elementary step of this trajectory, the infinitesimal duration is $dt = ds/\|\vec{u}(s)\|$, $\vec{u}$ being the fluid velocity at curvilinear abscissa $s$ of the trajectory. The effective speed $V_{\text{eff}}$ corresponds to the quantity measured in our experiments.

# Results

## Micro-bead tracking method

78 movies were recorded, corresponding to a total of 24 ciliated edges. From these movies, the trajectories of 195 micro-beads were retained. Three examples of MBT movie are displayed online (see supporting information).

Micro-beads velocities are essentially oriented along the $x^*$ direction (parallel to the ciliated edge), as can be seen in Fig 7. Only 4% of the micro-beads exhibit a vertical velocity component larger than 25% of the horizontal component. Moreover, this vertical component is about equally distributed among positive and negative values (98 vs 97). This result suggests that the bead velocities can be modeled in a good approximation as parallel to the ciliated edge.

Micro-bead velocities range from 0.0 to 253.8 µm.s⁻¹ (mean = 42.2 µm.s⁻¹, std = 35.0 µm.s⁻¹). This wide variation is explained by the spread of the distances of the micro-beads to the ciliated edge which range from 0.3 to 70.9 µm (mean = 12.7 µm; std = 11.1 µm). Indeed, as observed in Fig 8, left which presents micro-bead velocity measurements on three different ciliated edges, velocities appear to be strongly correlated to these distances, the fastest micro-beads being the ones closest to the ciliated edge. To confirm this observation, velocities of all micro-beads were separated into two groups, respectively above and below to the median value of the distance (Fig 8, right). The comparison between the two groups was performed with a statistical software package using non-parametric tests (Mann-Whitney U-test). A $p$ value <0.05 was considered significant. Here again, one can distinctly observe a clear link between the value of the micro-bead velocity and its distance to the ciliated edge.

## Correlation between micro-bead velocity and other ciliary beating measurements

In order to assess the parameters that influence significantly the micro-bead velocity, a multiple linear regression analysis was performed between the values of the micro-bead velocities (the dependent variable in the regression) and five exploratory variables, 4 of which are determined from measurements of ciliary beating (CBF, CBA, $\lambda$, $\rho_c$), and the fifth being the distance $y_0$ to the cilia. The resulting multiple regression equation of micro-bead velocities is given by:

$$V_{\text{eff}}(y_0) = (6.96 \times \text{CBF}) + (230.6 \times \rho_c) + (11.95 \times \text{CBA}) - (1.22 \times y_0) - (0.83 \times \lambda) - 243.35, \quad (6)$$

with a coefficient of determination equal to 0.68 and a probability smaller than $10^{-6}$ (F-test).
<figure><img src="image_12.png" /><figcaption>**Fig 7.** Vertical velocity component versus horizontal velocity component. Each square corresponds to one of the 195 micro-beads measured. Red, green and purple dashed lines correspond respectively to $|u_x| = |u_y|/10$, $|u_x| = |u_x|/4$ and $|u_x| = |u_x|/2$.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g007

Moreover, all probabilities associated to regression coefficients (t-test) are smaller than $10^{-6}$, except for the wavelength ($p = 0.048$). One can remark that this regression implicitly assumes a linear dependency on the exploratory variables, while a further and finer analysis using our mathematical model reveals a more complex behavior, as for instance the micro-bead velocity profile above the edge which is shown in fact to be essentially parabolic.

## Comparison between experimental measurement and numerical model

**Input of the numerical model.** To run simulations and compare their results with our MBT experiments, the numerical model presented in [25] requires 5 inputs: CBF, CBA, meta-chronal wavelength $\lambda$, distance $h$, and the slip length $\phi$. CBF, CBA, and $\lambda$ are directly measured by a microscopic analysis of the ciliary beating (see above). The distance $h$ has to be determined by fitting the parabolic velocity profile. Finally, the slip length $\phi$ is shown to be directly correlated to the cilia density.
<figure><img src="image_13.png" /><figcaption>**Fig 8.** Micro-bead velocities versus distance to the ciliated edge. (Left) Examples of mean velocity obtained in 25 micro-beads in 3 ciliated edges (blue, green and red). Each square corresponds to one micro-bead. Horizontal and vertical error bars display for each bead the standard deviation of the distance to cilia and velocity, respectively. Dashed lines are the parabolic regression on each ciliated edges. (Right) The 195 measured micro-beads are presented into 2 groups of equal sizes, according to their distance to the ciliated edge. Bars exhibit a significant velocity difference (p < 0.05). Error bars represented the standard error of the mean (SEM).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g008

*Distance h.* The mathematical model predicts that the micro-bead velocity exhibits a parabolic profile along *y* [25]. The distance *h* between the cilia wall and the region of stagnant fluid is read as the *x*-intercept of the parabolic fit of the measured velocity profile.

*Slip length φ.* The fifth and last parameter required as input to the numerical model is the slip length *φ* which cannot be directly retrieved from a microscopic measurement at the cilia wall. This parameter accounts for the partial momentum transfer between the wall and the fluid, due to the non continuous coverage of the cilia, through the following boundary condition:

$$ \vec{u} - \phi \frac{\partial \vec{u}}{\partial y} \Big|_{(x,y,t)} = \vec{u}_w(x,t), \qquad (7) $$

where $\vec{u}$ is the velocity and $\vec{u}_w(x, t)$ is the wall velocity. This condition is analogous to that of a fluid flow next to a porous wall, where the presence of pores reduces the transfer of momentum between the wall and the fluid [29, 30]. In this situation, *φ* can be understood an *effective slip length* [31]. The parameter *φ* can be retrieved from MBT experiments by fitting the measured micro-bead velocity profile above the cilia wall with a parabola. If *φ* = 0, one recovers a non slip boundary condition, while *φ* → +∞ corresponds to a perfect sliding condition. Sbragalia, et al. [32] showed that the effective slip length produced by a solid no-slip boundary with a distribution of free slip area was proportional to the product of the length scale of free slip by the fraction of surface covered by them. In our problem, if we assume that the tips of cilia are “no-slip” while the spaces between cilia are “free slip”, we expect the value of *φ* (extracted from the aforementioned fitting procedure) to be directly proportional to $(1 - \rho_c)^2/\rho_c$ where $\rho_c$ is the cilia density directly estimated from a grey level measurement. As shown in Fig 9, left, *φ* and $\rho_c$ appear closely related. This correlation shows that a sparse density of cilia is a main determinant for a poor momentum transfer to the fluid. A simple regression method leads to the
<figure><img src="image_14.png" /><figcaption>**Fig 9. Effective porosity parameter φ.** (Left) The effective slip length φ, obtained in each one of the 24 ciliated edges by fitting the measured micro-bead velocity with a parabola, plotted against cilia density. (Right) Effective slip length φ fitted from micro-bead measurements plotted against the effective slip length computed from the cilia density using Eq 8. The red line corresponds to the linear regression.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g009

following analytic relationship between these two parameters φ and ρc:

$$ \phi \approx \phi_{\text{calc}}(\rho_c) = \Lambda \frac{(1 - \rho_c)^2}{\rho_c} \quad \text{with} \quad \Lambda = 750 \ \mu\text{m} \qquad (8) $$

Fig 9, right, displays a plot of the effective slip length φ extracted from the fitting procedure and the calculated slip length φcalc obtained from the formula in Eq 8. Except for a couple outliers, most values appear to be very close.

**Micro-bead velocities.** We now compare micro-bead velocities measured by MBT with numerical simulations using the 4 input parameters deduced from the analysis of ciliary beating: CBF, CBA, metachronal wavelength λ, and φcalc (this last value being deduced from the measured cilia density). We add to these measured parameters the distance h fitted from MBT experiment. In Fig 10, left, the velocities Vb measured by MBT are plotted against the velocities Veff predicted by the numerical model. Each point corresponds to one micro-bead. The correlation coefficient R² between cilia density and slip length is found to be equal to 0.71. In Fig 10, right, the Bland-Altman plot displays the difference between the two velocities against their average, i.e., each point corresponds to $(\frac{V_b + V_{eff}}{2}, V_b - V_{eff})$. Again, experimental and numerical velocities are found to be in close agreement, proving that our model faithfully captures the relationship between fluid velocity and the microscopic measurements of cilia motion.

## Ciliary beating efficiency index

The shear stress τw exerted by the cilia wall on the fluid is given by (see [25]):

$$ \tau_w = \frac{2\mu}{h} U_w^* \qquad (9) $$
<figure><img src="image_15.png" /><figcaption>**Fig 10.** Comparison between measured and simulated micro-bead velocities. (Left) Experimental velocities against simulated velocities in the 195 measured microbeads. The red line corresponds to $y = x$. (Right) Bland-Altman plot. The blue dashed line corresponds to the mean value while the red dashed lines correspond to the mean value ± 1.96 standard deviation.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005605.g010

where $U_w^*$ is the extrapolated velocity at the ciliated wall, $h$ is the horizontal intercept (vanishing velocity) of the fitted parabolic profile, and $\mu$ is the dynamic viscosity of the fluid. We propose to use this force $\tau_w$ as an index for assessing the efficiency of the ciliary beating. MBT measurements enable us to retrieve $U_w^*$ and $h$, and therefore to compute the shear stress.

In Table 1 we report the results of the analysis performed on 11 patients. For each patient, one or several ciliated edges have been analyzed. For each edge, the number of tracked microbeads is reported, followed by measured values of the input parameters of the model. The last 3 columns show the values of $h$ and $U_w^*$, retrieved from the parabolic profiles, and the value of the deduced shear stress.

Discussion

The experimental MBT results (see Fig 10, left and right), allowed us to validate the numerical model. The model predicts that the micro-bead velocity is an increasing function of CBA and CBF while it decreases with $\lambda$ and $\phi$ (the last dependency means that micro-bead velocity increases with $\rho_c$). These behaviors are consistent with the various signs of the regression coefficients of Eq 6 deduced from MBT experiments.

The good agreement observed between the numerical model and experimental MBT results led us to propose the steady component of the shear stress exerted by the cilia wall on the fluid as a new index of the ciliary beating efficiency. To our knowledge, there exist very few evaluations of the shear stress exerted by the respiratory cilia in the literature. If we assume a cilia density per unit area of 5 cilia/µm² [9], we find a force per cilium equal to 0.013±0.009 pN. Such a value is much lower (at least 3 orders of magnitude) than the ones reported in a study on human bronchial epithelial cell culture [9], or in a study on culture grown from frog esophagus [33]. However, in these two last cases, the measured quantity was an oscillating force and not the steady component resulting from the entire beating as in our study. This suggests that these two types of measurement should not be directly compared. Moreover, the shear stress we find is of the same order of magnitude than the shear stress induced by the so-called
<table><caption>Table 1. Value of efficiency index for all ciliated edges.</caption><thead><tr><th>Patient<br>N<sub>o</sub></th><th>n</th><th>CBF<br>(Hz)</th><th>CBA<br>(µm)</th><th>λ<br>(µm)</th><th>ρ<sub>c</sub><br>(%)</th><th>h<br>(µm)</th><th>U<sub>w</sub><br>(µm.s<sup>-1</sup>)</th><th>τ<sub>w</sub><br>(mPa)</th><th>mean ± std dev.</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>8</td><td>14.0</td><td>7.5</td><td>13.0</td><td>84</td><td>87.5</td><td>138.3</td><td>3.2</td><td rowspan="4">5.8±4.1</td></tr><tr><td>4</td><td>15.0</td><td>6.8</td><td>14.6</td><td>86</td><td>52.6</td><td>240.9</td><td>9.2</td></tr><tr><td>5</td><td>13.5</td><td>7.5</td><td>19.8</td><td>89</td><td>47.4</td><td>224.8</td><td>9.5</td></tr><tr><td>5</td><td>13.6</td><td>4.0</td><td>28.0</td><td>82</td><td>70.9</td><td>49.4</td><td>1.4</td></tr><tr><td rowspan="3">2</td><td>3</td><td>5.3</td><td>7.0</td><td>10.8</td><td>73</td><td>7.0</td><td>27.6</td><td>0.7</td><td rowspan="3">2.5±1.6</td></tr><tr><td>3</td><td>5.6</td><td>6.6</td><td>11.3</td><td>80</td><td>32.9</td><td>46.9</td><td>2.9</td></tr><tr><td>9</td><td>6.3</td><td>7.3</td><td>11.0</td><td>68</td><td>28.9</td><td>55.2</td><td>3.8</td></tr><tr><td>3</td><td>6</td><td>14.0</td><td>5.4</td><td>10.5</td><td>69</td><td>53.0</td><td>51.3</td><td>1.9</td><td>1.9</td></tr><tr><td>4</td><td>4</td><td>11.5</td><td>7.5</td><td>12.0</td><td>74</td><td>66.3</td><td>76.8</td><td>2.3</td><td>2.3</td></tr><tr><td rowspan="5">5</td><td>3</td><td>7.3</td><td>7.8</td><td>16.0</td><td>58</td><td>82.1</td><td>24.0</td><td>0.6</td><td rowspan="5">2.0±1.3</td></tr><tr><td>3</td><td>9.4</td><td>7.0</td><td>15.5</td><td>80</td><td>91.5</td><td>80.1</td><td>1.8</td></tr><tr><td>3</td><td>4.9</td><td>8.4</td><td>17.3</td><td>88</td><td>53.5</td><td>64.7</td><td>2.4</td></tr><tr><td>4</td><td>5.6</td><td>6.9</td><td>13.3</td><td>95</td><td>52.8</td><td>106.5</td><td>4.0</td></tr><tr><td>3</td><td>8.1</td><td>6.7</td><td>15.4</td><td>82</td><td>139.7</td><td>85.4</td><td>1.2</td></tr><tr><td>6</td><td>4</td><td>3.7</td><td>6.2</td><td>16.7</td><td>89</td><td>34.2</td><td>33.6</td><td>2.0</td><td>2.0</td></tr><tr><td>7</td><td>4</td><td>7.7</td><td>7.2</td><td>12.0</td><td>70</td><td>42.6</td><td>49.3</td><td>2.3</td><td>2.3</td></tr><tr><td>8</td><td>3</td><td>3.6</td><td>6.3</td><td>17.7</td><td>98</td><td>40.1</td><td>32.3</td><td>1.6</td><td>1.6</td></tr><tr><td rowspan="2">9</td><td>12</td><td>7.7</td><td>8.8</td><td>12.1</td><td>62</td><td>60.0</td><td>51.7</td><td>1.7</td><td rowspan="2">1.9±0.3</td></tr><tr><td>5</td><td>6.6</td><td>7.0</td><td>10.3</td><td>73</td><td>52.3</td><td>54.4</td><td>2.1</td></tr><tr><td>10</td><td>3</td><td>6.7</td><td>6.5</td><td>10.4</td><td>87</td><td>28.6</td><td>60.2</td><td>4.2</td><td>4.2</td></tr><tr><td rowspan="4">11</td><td>56</td><td>3.6</td><td>4.7</td><td>10.2</td><td>96</td><td>51.2</td><td>58.9</td><td>2.3</td><td rowspan="4">2.1±0.3</td></tr><tr><td>11</td><td>3.5</td><td>5.3</td><td>12.6</td><td>89</td><td>48.6</td><td>45.6</td><td>1.9</td></tr><tr><td>6</td><td>4.9</td><td>4.8</td><td>18.5</td><td>89</td><td>50.5</td><td>44.3</td><td>1.8</td></tr><tr><td>28</td><td>6.0</td><td>5.2</td><td>16.7</td><td>86</td><td>45.4</td><td>56.4</td><td>2.5</td></tr></tbody></table>

N<sub>o</sub> corresponds to the patient number; n is the number of micro-beads per ciliated edges; CBF is the ciliary beat frequency; CBA is the ciliary beat amplitude; λ is the metachronal wavelength; h is the distance fitted from the micro-bead velocity profile; ρ<sub>c</sub> is the cilia density; U<sub>w</sub><sup>*</sup> is the extrapolated velocity at the ciliated wall through MBT; τ<sub>w</sub> is the shear stress exerted by the cilia wall on the fluid (see Eq 9). When several ciliated edges are measured for the same patient, the standard deviation of τ<sub>w</sub> is computed.

https://doi.org/10.1371/journal.pcbi.1005605.t001

Couette flow between two parallel flat plates:

$$ \tau_{\text{Couette}} = \frac{\mu}{e} V, \qquad (10) $$

where e and V are the distance and the relative velocity between the two plates, respectively. If we put on a par e with h, and V with U<sub>w</sub><sup>*</sup>, respectively, one finds the formula for the two stresses to be very similar. The factor 2 observed in our mode originates from the difference of velocity profiles (linear for a Couette flow, parabolic in our case).

The parameter β, i.e., the ratio between minor and major axis is very difficult to measure on a large part of our movies. We decided to set β at 0.14 for all simulations as it was the mean value observed on our movies. The measured values of all other parameters obtained by ciliary analysis in this study were consistent with the literature. As an example, our values of metachronal wavelength fall in the same range than the few numbers that can be found in the literature for other cellular models (paramecium [34, 35], frog oesophagus [36] or rabbit trachea [5]).

In contrast with other analyses of the ciliary beating [20], or evaluations of the beat pattern [19], τ<sub>w</sub> provides a direct estimation of the force that is potentially applied by the ciliated
epithelium on the surrounding fluid, hence on the mucus. As such, it appears as a very good
candidate for a global index of the potential ciliary beat efficiency. It encompasses the beating
of the ciliated edge as a whole rather than focusing on an individual cilium. In comparable
experimental conditions (here the cell survival medium at room temperature), τw allows us to
compare the ciliary beat efficiency of different patients with each other. Moreover, this index
that can be easily obtained from a MBT experiment via a simple parabolic fitting requires nei-
ther a subjective interpretation by an operator, nor a specific procedure in terms of human
data sampling, but only a classical nasal (or bronchial) brushing.

Low values of τ<sub>w</sub>, corresponding to a low ciliary beat efficiency can be explained by several distinct causes: ciliary beating parameters alterations (CBF, CBA, metachronal wavelength), reduced cilia density, or a loss of coordination between cilia. Indeed our model assumes, via the metachronal wave, a perfect cilia coordination as well as a pure rectilinear geometry of the edge. As a consequence, high values of τ<sub>w</sub> are expected to correspond to well coordinated cilia without degradation of ciliary beating parameters (CBF, CBA, metachronal wavelength). However, for low values of τ<sub>w</sub>, a ciliary beat pattern analysis may be required to find the reason of a degraded efficiency (ciliary beating parameters, loss of coordination, ...).

Table 1 displays values of τ<sub>w</sub> measured in several patients. These values exhibit a relatively high intra-patient heterogeneity which is not surprising as it was already observed in other characteristics of the ciliary beating, such as the “the distance traveled by the cilium tip in one second weighted by the percentage of beating ciliated edges” [20]. Despite its intra-patient variability, this distance is able to discriminate non-PCD from PCD patients with a specificity and sensitivity above 0.95 [20]. Similarly, using τ<sub>w</sub> as a screening index in clinical studies would probably require analyzing several edges per patient.

Interestingly, the patient exhibiting the highest values of τ<sub>w</sub> (patient №1, see Table 1) seems to be also the patient with the most normal clinical state. To this date, it is difficult to pinpoint a precise threshold of τ<sub>w</sub> that would allow one to discriminate between clinically healthy and pathological edges. Determining this threshold would require a clinical study including a larger cohort of control patients compared to patients with well defined pathological ciliary beating phenotype, a study far beyond the scope of the present study.

Ciliary beating parameters (CBF, CBA, pattern, ...) depend on the characteristics of the fluid at least in part. In our MBT experiments, cilia are beating in a fluid whose physical properties are very similar to water. Such conditions are far from the real conditions of airways coated with mucus. One may wonder how the modification of the surrounding fluid influences cilia beating. Several numerical studies have investigated the way viscosity influences ciliary beating and the difference between a Newtonian vs. non Newtonian surrounding fluid. Jayathilake *et al.* have explored numerically the effect of an increased PCL viscosity on the motion of cilia embedded, and have found that, for a given beating frequency (set around 10 Hz, consistent with our measurements), the velocities of the induced PCL flow were almost unaffected by a 5-fold increase of the viscosity [37]. Other models taking into account the thermodynamic characteristic of the mucus (generalized Newtonian fluid, ...) can be found in the literature [28]. More recently, Sedaghat *et al.* have studied in detail the two layer structure (PCL and mucus), PCL being modeled as a Newtonian fluid and mucus as a viscoelastic fluid [38]. Their study showed firstly that mucus viscosity has a very limited effect on the mucus flow, due to the dominating part of the elastic part in the viscosity, and secondly that ciliary beating frequency plays a major role in the mucus velocity. These results were confirmed by very recent 3D numerical simulations [39]. In these works, the beating frequency was identified as a major determinant of the mucus flow, before the mucus viscosity itself. However, in all simulations, this frequency is always a prescribed parameter whereas in reality it results from a delicate force and momentum balance in the fluid-structure interaction. As a
consequence, it is very complicated to really assess the net influence of mucus rheological
properties on cilia beating, considering the extreme difficulty to reproduce *in vitro* the *in vivo*
situation. The essential point here lies in the value of the beating frequency, which appears to
be very similar between *in vivo* data from the literature and our ex vivo experiments.

Replacing the surrounding fluid with water may also induce difference in biochemical
interactions. The mucociliary clearance process is generally described with three main compo-
nents: mucins secreted by goblet that will give mucus, cilia that move the mucus and an ion
transport process allowing maintaining an adequate aqueous environment on the airway epi-
thelium [40]. These components can be, at least partially, controlled by local agonist, extracel-
lular nucleotides, and nucleosides released from the epithelium. For example, it is known that
ATP, UTP, and adenosine increase the ciliary beat frequency while ATP and UTP stimulate
the secretion of mucins.

However, to our knowledge, the relation between these local agonists and the cilia beat pat-
tern remains to be established. In 2000 Chilvers *et al*. proposed to use digital high speed video
microscopy to visualize the beat pattern of human nasal cilia in the absence of mucus [19].
Since 2009, the European Respiratory Society strongly recommended to combine this tech-
nique (i.e., in the absence of mucus) with classical tests in order to ensure primary ciliary dys-
kinesia diagnosis [23, 41]. At this time, studying cilia beat pattern in the presence of mucus for
patient specific evaluation would be very difficult to carry out because it would require tracheal
or bronchial explants which are ethically problematic in clinical practice.

Clearly, the experimental procedure described in this study is not intended at evaluating
the effect of possible alterations of the cilia environment (change in rheological properties of
the mucus, modification of the nucleotides in the periciliary layer...). Our MBT experiments,
the developed numerical model and the proposed index aim essentially at evaluating patholo-
gies impairing mucociliary clearance resulting from a defect in ciliary motion (ciliopathies). In
future studies, one could even envision adding exogenous drug to the survival medium to test
the potential effect of a drug treatment. For now, the goal is to simulate the conditions of our
experiment corresponding to a measurement that can be realized in a clinical setting.

Conclusion

We have presented here a first experimental validation of the numerical model introduced in
[25], describing the Newtonian fluid flow induced by cilia motion above a ciliated edge. This
study shows that the computational model satisfactorily predicts the profile of velocities of
micro-beads in the Newtonian flow, this profile being essentially parabolic. Recovering this
profile allows us to assess the shear stress locally applied by the ciliated edge onto the fluid.
Our model suggests that this shear stress at the cilia wall characterizes the momentum transfer
between the cilia and the fluid, and thus the efficiency of the ciliary beating. Interestingly, the
estimation of this index does not require any modification of the present clinical practice of
data collection (nasal or bronchial brushing). This study opens the broad perspective of using
this index in the future to characterize ciliary function of transport during normal and patho-
logical conditions of either congenital origin such as PCD or acquired dyskinesia secondary to
pathogenic invasion and/or occupational exposures.

**Supporting information**

**S1 Video.** Example n°1 of MBT movie. Patient №1, third edge presenting a good coordinated edge and good cila density.
(MP4)
**S2 Video.** *Example n°2 of MBT movie.* Patient №4 presenting a good coordinated edge but relatively low cilia density.
(MP4)

**S3 Video.** *Example n°3 of MBT movie.* Patient №10 presenting miscoordinated edge with non measurable metachronal wavelength.
(MP4)

## Author Contributions

**Conceptualization:** Mathieu Bottier, Gabriel Pelle, Daniel Isabey, James B. Grotberg, Jean-François Papon, Bruno Louis.

**Data curation:** Mathieu Bottier, Gabriel Pelle, Bruno Louis.

**Formal analysis:** Mathieu Bottier, Gabriel Pelle, Bruno Louis.

**Funding acquisition:** Daniel Isabey, Bruno Louis.

**Investigation:** Mathieu Bottier, Sylvain Blanchon, Gabriel Pelle, Emilie Bequignon, Daniel Isabey, André Coste, Estelle Escudier, James B. Grotberg, Jean-François Papon, Marcel Filoche, Bruno Louis.

**Methodology:** Mathieu Bottier, Gabriel Pelle, Daniel Isabey, James B. Grotberg, Marcel Filoche, Bruno Louis.

**Project administration:** Daniel Isabey, Bruno Louis.

**Software:** Mathieu Bottier, Gabriel Pelle, Bruno Louis.

**Supervision:** Gabriel Pelle, Daniel Isabey, James B. Grotberg, Jean-François Papon, Marcel Filoche, Bruno Louis.

**Validation:** Mathieu Bottier, Sylvain Blanchon, Gabriel Pelle, Emilie Bequignon, André Coste, Estelle Escudier, James B. Grotberg, Jean-François Papon, Marcel Filoche, Bruno Louis.

**Visualization:** Mathieu Bottier, Gabriel Pelle, Bruno Louis.

**Writing – original draft:** Mathieu Bottier, Marcel Filoche, Bruno Louis.

**Writing – review & editing:** Mathieu Bottier, Gabriel Pelle, Emilie Bequignon, Daniel Isabey, André Coste, Estelle Escudier, James B. Grotberg, Jean-François Papon, Marcel Filoche, Bruno Louis.

## References

1. Fulford GR, Blake JR. Muco-ciliary Transport in the Lung. J Theor Biol. 1986; 121:381–402. https://doi.org/10.1016/S0022-5193(86)80098-4 PMID: 3796001

2. Wanner A, Salathe M, O'Riordan TG. Mucociliary clearance in the airways. Am J Respir Crit Care Med. 1996; 154:1868–1902. https://doi.org/10.1164/ajrccm.154.6.8970383 PMID: 8970383

3. Knowles MR, Boucher RC. Mucus clearance as a primary innate defense mechanism for mammalian airways. J Clin Invest. 2002; 109:571–577. https://doi.org/10.1172/JCI15217 PMID: 11877463

4. Fahy JV, Dickey BF. Airway mucus function and dysfunction. N Engl J Med. 2010; 363(23):2233–47. https://doi.org/10.1056/NEJMra0910061 PMID: 21121836

5. Sanderson MJ, Sleigh MA. Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and meta-chrony. J Cell Sci. 1981; 47:331–47. PMID: 7263784

6. Sleigh MA, Blake JR, Liron N. The propulsion of mucus by cilia. Am Rev Respir Dis. 1988; 137(13): 726–41. https://doi.org/10.1164/ajrccm/137.3.726 PMID: 3278666
7. Gheber L, Priel Z. On metachronism in ciliary systems: a model describing the dependence of the meta-chronal wave properties on the intrinsic ciliary parameters. Cell Motil Cytoskeleton. 1990; 16(3):167-81. https://doi.org/10.1002/cm.970160304 PMID: 2364445

8. Gheber L, Priel Z. Extraction of cilium beat parameters by the combined application of photoelectric measurements and computer simulation. Biophys J. 1997; 72:449-462. https://doi.org/10.1016/S0006-3495(97)78686-7 PMID: 8994632

9. Hill DB, Swaminathan V, Estes A, Cribb J, O'Brien ET, Davis CW, et al. Force Generation and Dynamics of Individual Cilia under External Loading. Biophys J. 2010; 98:57-66. https://doi.org/10.1016/j.bpj.2009.09.048 PMID: 20085719

10. Salathe M. Regulation of mammalian ciliary beating. Annu Rev Physiol. 2007; 69:401-22. https://doi.org/10.1146/annurev.physiol.69.040705.141253 PMID: 16945069

11. Lucas JS, Barbato A, Collins SA, Goutaki M, Behan L, Caudri D, et al. European Respiratory Society guidelines for the diagnosis of primary ciliary dyskinesia. European Respiratory Journal. 2017; 49(1). https://doi.org/10.1183/13993003.01090-2016

12. Canciani M, Barlocco EG, Mastella G, De-Santi MM, Gardi C, Lungarella G. The saccharin method for testing mucociliary function in patients suspected of having primary ciliary dyskinesia. Pediatr Pulmonol. 1988; 5(4):210-4. https://doi.org/10.1002/ppul.1950050406 PMID: 3237448

13. Braat J, Ainge G, Bowles J, Richards D, Van-Riessen D, Visser W, et al. The lack of effect of benzalkonium chloride on the cilia of the nasal mucosa in patients with perennial allergic rhinitis: a combined functional, light, scanning and transmission electron microscopy study. Clin Exp Allergy. 1995; 25(10):957-65. https://doi.org/10.1111/j.1365-2222.1995.tb00398.x PMID: 8556567

14. Yeates DB, Pitt BR, Spektor DM, Karron GA, Albert RE. Coordination of mucociliary transport in human trachea and intrapulmonary airways. J Appl Physiol. 1981; 51(5):1057-64. PMID: 7298443

15. Del-Donno M, Pavia D, Agnew JE, Lopez-Vidriero MT, Clarke SW. Variability and reproducibility in the measurement of tracheobronchial clearance in healthy subjects and patients with different obstructive lung diseases. Eur Respir J. 1988; 1(7):613-20. PMID: 3181408

16. Mortensen J. Bronchoscintigraphy and pulmonary clearance of 99mTc-albumin colloid in the study of mucociliary clearance. Dan Med Bull. 1998; 45(2):195-209. PMID: 9587703

17. Munkholm M, Nielsen KG, Mortensen J. Clinical value of measurement of pulmonary radioaerosol mucociliary clearance in the work up of primary ciliary dyskinesia. EJNMMI Res. 2015; 5:39. https://doi.org/10.1186/s13550-015-0118-y

18. Liu L, Chu K, Houser G, Diephuis B, Li Y, Wilstreman E, et al. Method for Quantitative Study of Airway Functional Microanatomy Using Micro-Optical Coherence Tomography. Plos One. 2013; 8(1):e54473. https://doi.org/10.1371/journal.pone.0054473 PMID: 23372732

19. Chilvers MA, O'Callaghan C. Analysis of ciliary beat pattern and beat frequency using digital high speed imaging: comparison with the photomultiplier and photodiode methods. Thorax. 2000; 55:314-317. https://doi.org/10.1136/thorax.55.4.314 PMID: 10722772

20. Papon JF, Bassinet L, Cariou-Patron G, Zerah-Lancner F, Vojtek AM, Blanchon S, et al. Quantitative analysis of ciliary beating in primary ciliary dyskinesia: a pilot study. Orphanet J Rare Dis. 2012; 7:78. https://doi.org/10.1186/1750-1172-7-78 PMID: 23057704

21. Amirav I, Mussaffi H, Roth Y, Schmidts M, Omran H, Werner C. A reach-out system for video microscopy analysis of ciliary motions aiding PCD diagnosis. BMC Research Notes. 2015; 8:71. https://doi.org/10.1186/s13104-015-0999-x PMID: 25869032

22. Quinn SP, Zahid MJ, Durkin JR, Francis RJ, Lo CX, Chennubhotla SC. Automated identification of abnormal respiratory ciliary motion in nasal biopsies. Sci Transl Med. 2015; 7(299):299ra124. https://doi.org/10.1126/scitranslmed.aaa1233 PMID: 26246169

23. Jackson CL, Behan L, Collins SA, Goggin PM, Adam EC, Coles JL, et al. Accuracy of diagnostic testing in primary ciliary dyskinesia. Eur Respir J. 2016; 47(3):837-848. https://doi.org/10.1183/13993003.00749-2015 PMID: 26647444

24. Bottier M, Blanchon S, Filoche M, Isabey D, Coste A, Escudier E, et al. Characterization of upper airway ciliary beat by coupling isolated and collective cilia motion analysis. Cilia. 2015; 4(Suppl 1):P86.

25. Bottier M, Pena-Fernández M, Pelle G, Isabey D, Grotberg JB, Louis B, et al. A New Index for Characterizing Micro-bead Motion in a Flow Induced Ciliary Beating: part II, Modeling. PLoS Comput Biol. 2017.

26. Thomas B, Rutman A, O'Callaghan C. Disrupted ciliated epithelium shows slower ciliary beat frequency and increased dyskinesia. Eur Respir J. 2009; 34:401-404. https://doi.org/10.1183/09031936.00153308 PMID: 19648518

27. Smith DJ, Gaffney EA, Blake JR. Modelling mucociliary clearance. Resp Physiol Neurobiol. 2008; 163:178-188. https://doi.org/10.1016/j resp.2008.03.006
28. Velez-Cordero JR, Lauga E. Waving transport and propulsion in a generalized Newtonian fluid. J Non-Newton Fluid. 2013; 199:37–50. https://doi.org/10.1016/j.jnnfm.2013.05.006

29. Beavers GS, Joseph DD. Boundary conditions at a naturally permeable wall. J Fluid Mech. 1967; 30(1): 197–207. https://doi.org/10.1017/S0022112067001375

30. Verma VS, Tripathee SM. A planar model for muco-ciliary transport in the human lung: effects of mucus visco-elasticity, cilia beating and porosity. IJMMPS. 2013; 1(1).

31. Zhang XP, Lund NJ, Hendy SC. Effective slip length: some analytical and numerical results. ANZIAM J. 2015; 7:79–88. https://doi.org/10.1017/S1446181115000127

32. Sbragaglia M, Prosperetti A. Effective velocity boundary condition at a mixed slip surface. J Fluid Mech. 2007; 578:435–451. https://doi.org/10.1017/S0022112007005149

33. Teff Z, Priel Z, Gheber LA. Forces applied by cilia measured on explants from mucociliary tissue. Bio-phys J. 2007; 92:1813–1823. https://doi.org/10.1529/biophysj.106.094698 PMID: 17142280

34. Machemer H. Ciliary activity and the origin of metachrnoy in paramecium: effects of incresead viscosity. J Exp Biol. 1972; 57:239–259. PMID: 5075893

35. Okamoto K, Nakaka Y. Reconstitution of metachronal waves in ciliated cortical sheets of paramecium. J Exp Biol. 1994; 192:73–81. PMID: 9317363

36. Gheber L, Priel Z. Metachronal activity of cultured mucociliary epithelium under normal and stimulated conditions. Cell Motil Cytoskeleton. 1994; 28(4):333–45. https://doi.org/10.1002/cm.970280407 PMID: 7954860

37. Jayathilake PG, Tan Z, Le DV, Lee HP, Khoo BC. Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method. Comp Fluids. 2012; 67:130–137. https://doi.org/10.1016/j.compfluid.2012.07.016

38. Sedaghat MH, Shahmardan MM, Norouzi M, Jayathilake PG, Nazari M. Numerical simulation of muco-ciliary clearance: immersed boundary-lattice Boltzmann method. Comp Fluids. 2016; 131:91–101. https://doi.org/10.1016/j.compfluid.2016.03.015

39. Chatelin R, Poncet P. A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus. J Biomech. 2017; 49:1772–1780. https://doi.org/10.1016/j.jbiomech.2016.04.009

40. Davis CW, Lazarowski E. Coupling of airway ciliary activity and mucin secretion to mechanical stresses by purinergic signaling. Respiratory Physiology and Neurobiology. 2008; 163:208–213. https://doi.org/10.1016/j resp.2008.05.015 PMID: 18635403

41. Barbato A, Frischer T, Kuehni CE, Snijders D, Azevedo I, Baktai G, et al. Primary ciliary dyskinesia: a consensus statement on diagnostic and treatment approaches in children. Eur Respir J. 2009; 34:1264–1276. https://doi.org/10.1183/09031936.00176608 PMID: 19948909