Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.13(7); 2017 July**|**PMC5510810

Formats

Article sections

- Abstract
- Author summary
- Introduction
- Materials and methods
- Results
- Discussion
- Supporting information
- References

Authors

Related links

PLoS Comput Biol. 2017 July; 13(7): e1005552.

Published online 2017 July 14. doi: 10.1371/journal.pcbi.1005552

PMCID: PMC5510810

Mathieu Bottier, Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing,^{1,}^{2,}^{3} Marta Peña Fernández, Conceptualization, Formal analysis, Investigation, Software, Visualization,^{1,}^{2,}^{3} Gabriel Pelle, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing,^{1,}^{2,}^{3} Daniel Isabey, Funding acquisition, Investigation, Supervision, Writing – review & editing,^{1,}^{2,}^{3} Bruno Louis, Formal analysis, Funding acquisition, Investigation, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing,^{1,}^{2,}^{3} James B. Grotberg, Conceptualization, Formal analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing,^{4} and Marcel Filoche, Conceptualization, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing^{1,}^{2,}^{3,}^{5,}^{*}

Mark Alber, Editor^{}

University of California Riverside, UNITED STATES,

The authors have declared that no competing interests exist.

* E-mail: ude.euqinhcetylop@ehcolif.lecram

Received 2017 January 18; Accepted 2017 May 4.

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.

See "A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part I, experimental analysis" in volume 13, e1005605.

This article has been cited by other articles in PMC.

Mucociliary clearance is one of the major lines of defense of the human respiratory system. The mucus layer coating the airways is constantly moved along and out of the lung by the activity of motile cilia, expelling at the same time particles trapped in it. The efficiency of the cilia motion can experimentally be assessed by measuring the velocity of micro-beads traveling through the fluid surrounding the cilia. Here we present a mathematical model of the fluid flow and of the micro-beads motion. The coordinated movement of the ciliated edge is represented as a continuous envelope imposing a periodic moving velocity boundary condition on the surrounding fluid. Vanishing velocity and vanishing shear stress boundary conditions are applied to the fluid at a finite distance above the ciliated edge. The flow field is expanded in powers of the amplitude of the individual cilium movement. It is found that the continuous component of the horizontal velocity at the ciliated edge generates a 2D fluid velocity field with a parabolic profile in the vertical direction, in agreement with the experimental measurements. Conversely, we show than this model can be used to extract microscopic properties of the cilia motion by extrapolating the micro-bead velocity measurement at the ciliated edge. Finally, we derive from these measurements a scalar index providing a direct assessment of the cilia beating efficiency. This index can easily be measured in patients without any modification of the current clinical procedures.

Mucociliary clearance is the first line of defense mechanisms of the human airways. The 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 a mathematical and numerical model aiming at defining a global index for assessing the efficiency of this beating. Numerical simulations show that the bead velocity parallel to the wall varies according a parabolic profile with the distance to the wall. The velocity extrapolated at the wall is demonstrated to be a measurement of the momentum transfer between cilia and the surrounding fluid. This model allows us to interpret experimental measurements performed in a companion article and to propose a universal index characterizing the beating efficiency, which can be extracted in the current clinical setting.

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 (Fig 1, left). 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^{2} [6, 9]. A 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 several mm/min [10].

This work intends to build a mathematical model of the experiments described in a companion paper [11], and to derive from this model an efficient and reliable method to assess the cilia beating efficiency in a clinical setting. In the experiments, ciliated cell clusters issued from nasal brushing are immersed in a cell survival medium devoid of mucus, pushing forward the medium as they beat. One has to stress here that the absence of mucus is an important ingredient of the experimental setting and the developed mathematical model. Polystyrene micro-beads are then used as massless tracers to visualize and quantify the fluid velocity field around the cilia. The theoretical and numerical work presented here therefore aim at a quantitative modeling of the ciliary beating, then of the fluid flow generated by it.

In the literature, two main types of ciliary beating models can be found: *discrete-cilia* models and *volume-force* models [12]. In *discrete-cilia* model, each cilium is modeled as a discrete body and its shape is parametrized along its stroke period [1, 13, 14]. *Discrete-cilia* models are themselves divided into two types: in the *prescribed beating* models, the cilium motion is imposed as an input to the simulation [15]; in the *couple-internal mechanics/fluid-structure interaction* models [16, 17], cilia motions originate from the coupling between the internal structure of cilia and the external viscous forces. In contrast, in *volume-force* models cilia are modeled through a phenomenological continuous force distribution, varying in space and time as the cilia beat [18–20]. In this second type of model, the envelope modeling approach accounts for the formation of metachronal waves above the cilia layer [21]. The cilia tips are seen from the fluid as a “wavy wall”, hereby ignoring the details of the sub-layer dynamics [22, 23] (see Fig 1, right).

Many studies have addressed experimentally [24–26] and numerically [27–31] the effect of fluid visco-elasticity on transport and locomotion. We want to stress here that in our model, no finite layer of viscous mucus sits on top of the cilia, since they are surrounded by a semi-infinite layer of watery fluid.

In the following, we first compute the wave envelope boundary condition from the cilia motions, based on the work of Ross [22]. We then derive the non linear equations for a fluid flow periodic in the direction of the metachronal wave. These equations are expanded in *ε*, which is the ratio of the cilium amplitude to the fluid layer thickness, then solved using a Fourier decomposition. The steady contribution to the flow field in the vertical direction above the cilia is shown to exhibit a parabolic profile, to a very good approximation. We finally show that measuring microbead velocities as a function of the distance to the ciliated edge enables us to compute a scalar index which accounts for the transfer of momentum between the cilia and the fluid, and therefore assesses the cilia beating efficiency.

We present here a two-dimensional model of cilia, fluid, and micro-beads motion.

Each individual cilium is assumed to undergo a periodic motion in which its tip follows an elliptic trajectory, see Fig 1, 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 (Fig 1, right) [22, 23]. The *x** axis is chosen parallel to the ciliated edge, each ciliary tip being located around *y** = 0 on average (the ‘*’ notation stands for dimensioned quantities, and will be removed once we switch to dimensionless quantities). The tip of a cilium located at the horizontal coordinate *ξ** is assumed to follow a periodic elliptic trajectory centered in (*ξ**, 0) during each elementary beat (Fig 2). At time *t**, the tip coordinates $({X}_{w}^{*},{Y}_{w}^{*})$ are

$$\{\begin{array}{cc}{X}_{w}^{*}\hfill & ={\xi}^{*}-acos\left(\omega {t}^{*}\right)\hfill \\ {Y}_{w}^{*}\hfill & =\beta asin\left(\omega {t}^{*}\right)\hfill \end{array}$$

(1)

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

To reproduce the backwards traveling metachronal wave of wavelength λ, we introduce in the periodic motion of each cilium a phase shift $\frac{2\pi {\xi}^{*}}{\lambda}$ which linearly depends on the cilium position:

$$\{\begin{array}{cc}{X}_{w}^{*}\hfill & ={\xi}^{*}-acos\phantom{\rule{1pt}{0ex}}\left({\displaystyle \omega {t}^{*}+\frac{2\pi {\xi}^{*}}{\lambda}}\right)\hfill \\ {Y}_{w}^{*}\hfill & =\beta asin\phantom{\rule{1pt}{0ex}}\left({\displaystyle \omega {t}^{*}+\frac{2\pi {\xi}^{*}}{\lambda}}\right)\hfill \end{array}$$

(2)

The corresponding wave frequency is *f* = 2*πω* and its speed is *c* = *f*λ. By setting the space and time units respectively to be *h* and *ω*^{−1}, dimensionless parameters *ε* and *k*, and variables (*X*_{w}, *Y*_{w}) are introduced:

$$\epsilon =\frac{a}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}k=\frac{2\pi h}{\lambda}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}{X}_{w}=\frac{{X}_{w}^{*}}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}{Y}_{w}=\frac{{Y}_{w}^{*}}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}x=\frac{{x}^{*}}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}y=\frac{{y}^{*}}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}\xi =\frac{{\xi}^{*}}{h}\phantom{\rule{3.33333pt}{0ex}},\phantom{\rule{3.33333pt}{0ex}}t=\omega {t}^{*}$$

(3)

The motion equations of each tip are thus rewritten in a dimensionless form:

$$\{\begin{array}{cc}{X}_{w}(\xi ,t)\hfill & =\xi -\epsilon cos(k\xi +t)\hfill \\ {Y}_{w}(\xi ,t)\hfill & =\beta \epsilon sin(k\xi +t)\hfill \end{array}$$

(4)

In the following, we assume that the ciliary beating amplitude is much smaller than the film thickness, i.e. *ε* 1. In the limit of continuously distributed cilia along the *x* axis, the envelope of their tip motions forms a continuous boundary that generates a forcing on the fluid layer above the cilia. Starting from Eq 4, we now express the position of a particle of this envelope (or wall) in the Eulerian frame of reference (*x*, *y*, *t*). A tip located at position *x* at time *t* corresponds to a cilium centered in (*ξ*, 0) such that:

(5)

This equation shows that (*x* − *ξ*) is of order *ε*. Therefore the vertical location of the cilia wall *y*_{w}(*x*, *t*), which is also the y-coordinate of the corresponding cilium, *Y*_{w}(*ξ*, *t*), can be developed around *ξ* = *x* using a Taylor expansion:

$${y}_{w}(x,t)={Y}_{w}(\xi ,t)={Y}_{w}(x,t)+(\xi -x)\frac{\partial {Y}_{w}}{\partial \xi}\phantom{\rule{1pt}{0ex}}{|}_{\xi =x}+\frac{{(\xi -x)}^{2}}{2}\frac{{\partial}^{2}{Y}_{w}}{\partial {\xi}^{2}}{|}_{\xi =x}+\cdots $$

(6)

Since both *Y*_{w} and (*ξ* − *x*) are first order in *ε* (Eqs 4 and 5), the expansion to the second order in *ε* of *y*_{w} is

(7)

Now that we have determined the location of the ciliated wall at time *t*, we can derive its velocity. The horizontal component of the wall velocity is obtained as the time derivative of the Lagrangian velocity of the tip at location *x* at time *t*, while its vertical component is calculated from the time derivative of the vertical coordinate of the wall at position *x* and time *t*:

$$\{\begin{array}{cc}{u}_{w}\left(x,t\right)& ={\left(\frac{\partial {X}_{w}}{\partial t}\right)}_{\xi}=\epsilon \text{sin}\left(k\xi +t\right)\hfill \\ {v}_{w}\left(x,t\right)& =\frac{\partial {y}_{w}}{\partial t}=\epsilon \beta \text{cos}\left(kx+t\right)-{\epsilon}^{2}\beta k\text{sin}\left(2\left(kx+t\right)\right)\hfill \end{array}$$

(8)

The value of *ξ* to insert in the above horizontal velocity is given by Eq 5. Expanding the horizontal component of the wall velocity in Taylor series to the second order in *ε* yields

$$\begin{array}{ccc}\hfill {u}_{w}(x,t)& =\hfill & \frac{\partial {X}_{w}}{\partial t}\phantom{\rule{1pt}{0ex}}{|}_{\xi =x}+(\xi -x)\frac{\partial {u}_{w}}{\partial \xi}{|}_{\xi =x}+\phantom{\rule{3.33333pt}{0ex}}\dots \hfill \\ & =& \epsilon sin(kx+t)+[\epsilon cos(kx+t)][\epsilon kcos(kx+t)]\hfill \\ & =& \epsilon sin(kx+t)+{\epsilon}^{2}k{cos}^{2}(kx+t)\hfill \end{array}$$

(9)

In summary, at location *x* in the Eulerian frame, the second order expansion in *ε* of the vertical position *y*_{w} and of the velocity vector (*u*_{w}, *v*_{w}) of the ciliated wall are:

$${y}_{w}(x,t)=\epsilon \beta sin\theta +\frac{{\epsilon}^{2}\beta k}{2}\phantom{\rule{1pt}{0ex}}(1+cos(2\theta \left)\right)$$

(10)

$${u}_{w}(x,t)=\epsilon sin\theta +\frac{{\epsilon}^{2}k}{2}\phantom{\rule{1pt}{0ex}}(1+cos(2\theta \left)\right)$$

(11)

(12)

where *θ* = *kx* + *t* is the local phase of the metachronal wave. We now introduce the first and second orders of the *ε*-expansion of the wall velocity, called ${\overrightarrow{U}}_{w,1}=({u}_{w,1},{v}_{w,1})$ and ${\overrightarrow{U}}_{w,2}=({u}_{w,2},{v}_{w,2})$, respectively, defined such that:

$${\overrightarrow{U}}_{w}=\epsilon \phantom{\rule{3.33333pt}{0ex}}{\overrightarrow{U}}_{w,1}+\frac{{\epsilon}^{2}}{2}\phantom{\rule{3.33333pt}{0ex}}{\overrightarrow{U}}_{w,2}$$

(13)

In summary, the first and second orders in *ε* of the location and velocities of the cilia wall are:

$$\{\begin{array}{cc}{y}_{w,1}\left(\theta \right)\hfill & {\displaystyle =\beta sin\theta =\frac{\beta}{2i}{e}^{i\theta}-\frac{\beta}{2i}{e}^{-i\theta}}\hfill \\ {y}_{w,2}\left(\theta \right)\hfill & {\displaystyle =\beta k\phantom{\rule{1pt}{0ex}}(1+cos(2\theta \left)\right)=\beta k+\frac{\beta k}{2}{e}^{2i\theta}+\frac{\beta k}{2}{e}^{-2i\theta}}\hfill \\ {u}_{w,1}\left(\theta \right)\hfill & {\displaystyle =sin\theta =\frac{1}{2i}{e}^{i\theta}-\frac{1}{2i}{e}^{-i\theta}}\hfill \\ {u}_{w,2}\left(\theta \right)\hfill & {\displaystyle =k\phantom{\rule{1pt}{0ex}}(1+cos(2\theta \left)\right)=k+\frac{k}{2}{e}^{2i\theta}+\frac{k}{2}{e}^{-2i\theta}}\hfill \\ {v}_{w,1}\left(\theta \right)\hfill & {\displaystyle =\beta cos\theta =\frac{\beta}{2}{e}^{i\theta}+\frac{\beta}{2}{e}^{-i\theta}}\hfill \\ {v}_{w,2}\left(\theta \right)\hfill & =-2\beta ksin\left(2\theta \right)=i\beta k\phantom{\rule{3.33333pt}{0ex}}{e}^{2i\theta}-i\beta k\phantom{\rule{3.33333pt}{0ex}}{e}^{-2i\theta}\hfill \end{array}$$

(14)

We can observe therefore that all velocity terms at the ciliated wall are oscillatory, except one steady contribution to the horizontal velocity at the second order in *ε* that appears in *u*_{w,2}. This contribution is proportional to *k* and will be the origin of the steady horizontal motion of the fluid above the cilia.

We now compute the oscillatory flow field of a fluid of density *ρ* and viscosity *μ* in the channel comprised between *y* = *y*_{w}(*x*, *t*) and *y* = *h* in the vertical direction. Due to the periodic nature of the forcing from the wall, we will consider periodic solution in the *x* direction. Given the normalization chosen in Eq 3, the normalization factors for velocity and pressure are *hω* and *μω*, respectively. Hence, the dimensionless Navier-Stokes equation reads

$${\alpha}^{2}\phantom{\rule{1pt}{0ex}}(\frac{\partial \overrightarrow{U}}{\partial t}+(\overrightarrow{U}\xb7\overrightarrow{\nabla})\phantom{\rule{3.33333pt}{0ex}}\overrightarrow{U})=-\overrightarrow{\nabla}p+\Delta \overrightarrow{U}\phantom{\rule{3.33333pt}{0ex}},$$

(15)

where $\overrightarrow{U}=(u,v)$ is the dimensionless velocity field, *p* is the dimensionless pressure, and *α* is the Womersley number defined by $\alpha}^{2}=\frac{\rho {h}^{2}\omega}{\mu$. In the limit *α* → 0, one recovers the classical stationary Stokes equation. In our case, typical values for the channel thickness and beating frequency are *h* = 50 *μ*m and *f* = 10 Hz. In water (*ρ* = 1 g.cm^{−3} and *μ* = 1 cP), the above expression leads to a Womersley number *α* ≈ 0.4, which means *α*^{2} of the order of 0.1.

Boundary conditions are periodic in the horizontal direction with a period equal to the wavelength of the metachronal wave (corresponding here to λ/*h* = 2*π*/*k*):

$$\overrightarrow{U}\phantom{\rule{1pt}{0ex}}\left({\displaystyle \frac{2\pi}{k},y,t}\right)=\overrightarrow{U}(0,y,t)$$

(16)

At the upper side of the domain (*y* = 1), one assumes vanishing vertical velocity and vanishing shear stress due to the presence of stagnant fluid above the channel:

$$v(x,1,t)=0\phantom{\rule{2.em}{0ex}}\text{and}\phantom{\rule{2.em}{0ex}}\frac{\partial u}{\partial y}{|}_{(x,1,t)}=0$$

(17)

Next to the cilia envelope (*y* = *y*_{w}), we introduce a Robin-type boundary condition, which involves both the velocity and its normal derivative (Eq 18). The factor in front of the normal derivative, called *ϕ*, accounts for the partial momentum transfer between the wall and the fluid due to the non continuous coverage of the cilia. *ϕ* = 0 corresponds to a no slip boundary condition, while *ϕ* → +∞ corresponds to a vanishing shear stress at the wall (perfect sliding). 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 [32, 33]:

$$\left({\displaystyle \overrightarrow{U}-\varphi \frac{\partial \overrightarrow{U}}{\partial y}}\right){|}_{(x,{y}_{w},t)}={\overrightarrow{U}}_{w}(x,t)$$

(18)

If expressed in dimensional quantities, the parameter *ϕ* becomes a length *ϕ** = *hϕ*. For a porous wall, this length is proportional to the square root of the permeability. In our case, it is interpreted as a characteristic sliding length and will be related to the cilia surface density.

Expanding the flow field $\overrightarrow{U}$ in powers of *ε*, $\overrightarrow{U}=\epsilon \phantom{\rule{3.33333pt}{0ex}}{\overrightarrow{U}}_{1}+\frac{{\epsilon}^{2}}{2}\phantom{\rule{3.33333pt}{0ex}}{\overrightarrow{U}}_{2}$, the boundary conditions can be expressed at all orders in *ε*. Since *y*_{w} is first order in *ε* (see Eq 10), the boundary condition can be expanded around *y* = 0 both for the velocity $\overrightarrow{U}$ and its normal derivative $\partial \overrightarrow{U}/\partial y$:

$$\begin{array}{ccc}\hfill \overrightarrow{U}\left({y}_{w}\right)& =\hfill & \overrightarrow{U}\left(0\right)+{y}_{w}\frac{\partial \overrightarrow{U}}{\partial y}{|}_{y=0}+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \\ & =& \epsilon {\overrightarrow{U}}_{1}\left(0\right)+\frac{{\epsilon}^{2}}{2}{\overrightarrow{U}}_{2}\left(0\right)+[\epsilon {y}_{w,1}+\frac{{\epsilon}^{2}}{2}{y}_{w,2}]\phantom{\rule{1pt}{0ex}}[\epsilon \frac{\partial {\overrightarrow{U}}_{1}}{\partial y}{|}_{0}+\frac{{\epsilon}^{2}}{2}\frac{\partial {\overrightarrow{U}}_{2}}{\partial y}{|}_{0}]+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \\ & =& \epsilon \phantom{\rule{3.33333pt}{0ex}}{\overrightarrow{U}}_{1}\left(0\right)+\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}[{\overrightarrow{U}}_{2}\left(0\right)+2{y}_{w,1}\frac{\partial {\overrightarrow{U}}_{1}}{\partial y}]+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill \frac{\partial \overrightarrow{U}}{\partial y}{|}_{{y}_{w}}& =\hfill & \frac{\partial \overrightarrow{U}}{\partial y}{|}_{0}+{y}_{w}\frac{{\partial}^{2}\overrightarrow{U}}{\partial {y}^{2}}{|}_{0}+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \\ & =& \epsilon \frac{\partial {\overrightarrow{U}}_{1}}{\partial y}{|}_{0}+\frac{{\epsilon}^{2}}{2}\frac{\partial {\overrightarrow{U}}_{2}}{\partial y}{|}_{0}+[\epsilon {y}_{w,1}+\frac{{\epsilon}^{2}}{2}{y}_{w,2}]\phantom{\rule{1pt}{0ex}}[\epsilon \frac{{\partial}^{2}{\overrightarrow{U}}_{1}}{\partial {y}^{2}}{|}_{0}+\frac{{\epsilon}^{2}}{2}\frac{{\partial}^{2}{\overrightarrow{U}}_{2}}{\partial {y}^{2}}{|}_{0}]+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \\ & =& \epsilon \phantom{\rule{3.33333pt}{0ex}}\frac{\partial {\overrightarrow{U}}_{1}}{\partial y}{|}_{0}+\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}[\frac{\partial {\overrightarrow{U}}_{2}}{\partial y}{|}_{0}+2{y}_{w,1}\frac{{\partial}^{2}{\overrightarrow{U}}_{1}}{\partial {y}^{2}}{|}_{0}]+\mathcal{O}\left({\epsilon}^{3}\right)\hfill \end{array}$$

(20)

Inserting Eqs 19 and 20 into Eq 18 finally gives the first two orders of the *ε*-expansion of the boundary condition at the ciliated wall:

$${\overrightarrow{U}}_{1}\left(0\right)-\varphi \frac{\partial {\overrightarrow{U}}_{1}}{\partial y}{|}_{0}={\overrightarrow{U}}_{w,1}$$

(21)

$${\overrightarrow{U}}_{2}\left(0\right)-\varphi \frac{\partial {\overrightarrow{U}}_{2}}{\partial y}{|}_{0}={\overrightarrow{U}}_{w,2}-2{y}_{w,1}(\frac{\partial {\overrightarrow{U}}_{1}}{\partial y}{|}_{0}-\varphi \frac{{\partial}^{2}{\overrightarrow{U}}_{1}}{\partial {y}^{2}}{|}_{0})$$

(22)

Since the flow is two-dimensional and incompressible, a natural formulation of the problem is obtained by introducing the stream function *ψ* such that:

$$u=\frac{\partial \psi}{\partial y}\phantom{\rule{2.em}{0ex}}\text{and}\phantom{\rule{2.em}{0ex}}v=-\frac{\partial \psi}{\partial x}$$

(23)

Such a solution automatically fulfills the continuity equation $\text{div}\left(\overrightarrow{U}\right)=0$. The dimensionless Navier-Stokes equation now reads:

$$\{\begin{array}{cc}{\alpha}^{2}({\psi}_{yt}+{\psi}_{y}{\psi}_{yx}-{\psi}_{x}{\psi}_{yy})\hfill & =-{p}_{x}+{\psi}_{yxx}+{\psi}_{yyy}\hfill \\ {\alpha}^{2}(-{\psi}_{xt}+{\psi}_{y}{\psi}_{xx}+{\psi}_{x}{\psi}_{xy})\hfill & =-{p}_{y}-{\psi}_{xxx}-{\psi}_{xyy}\hfill \end{array}$$

(24)

We derive a Partial Differential Equation (PDE) in *ψ* only by taking the y-derivative of the first equation and subtracting it to it the x-derivative of the second equation:

(25)

This equation can finally be rewritten in the following compact form:

(26)

We now expand the stream function to the second order in *ε*, as it was done for the flow field:

$$\psi =\epsilon {\psi}_{1}+\frac{{\epsilon}^{2}}{2}{\psi}_{2}$$

(27)

The zero order term *ψ*_{0} vanishes since the fluid is set into motion only by the cilia beating (*ε* > 0). The two orders in *ε* are solved in sequence. The first order of Eq 26 is:

(28)

while the second order is:

(29)

The boundary condition at the lower wall (the cilia envelope) couples the first and second orders *ψ*_{1} and *ψ*_{2} through the Robin condition of Eq 18. The two velocity components and the two orders in *ε* translate into 4 equations at the wall boundary:

(30)

(31)

(32)

(33)

The stream function can thus be computed in first approximation as if the domain were a straight channel comprised between *y* = 0 and *y* = 1.

*ψ*_{1}(*x*, *y*, *t*) is a periodic function of period 2*π*/*k* along *x* that solves the following system:

$$\{\begin{array}{c}{\Delta}^{2}{\psi}_{1}-{\alpha}^{2}\Delta {\psi}_{1,t}=0\hfill \\ {\psi}_{1,y}-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,yy}=sin\theta \phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=0\hfill \\ {\psi}_{1,x}-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,xy}=-\beta cos\theta \phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=0\hfill \\ {\psi}_{1,yy}=0\phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=1\hfill \\ {\psi}_{1,x}=0\phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=1\hfill \end{array}$$

(34)

We expand *ψ*_{1} in Fourier series along the *x* direction, each term corresponding to a wave traveling forward or backward at dimensionless speed 1.

$${\psi}_{1}(x,y,t)=\sum _{n=-\infty}^{+\infty}\phantom{\rule{1pt}{0ex}}{a}_{n}\left(y\right)\phantom{\rule{3.33333pt}{0ex}}{e}^{in(kx+t)}$$

(35)

Due to the linearity of the PDE in Eq 34, each term of this Fourier series satisfies the same equation, which implies:

$$({a}_{n}^{\prime \prime \prime \prime}-2{k}^{2}{n}^{2}{a}_{n}^{\prime \prime}+{k}^{4}{n}^{4}{a}_{n})-in{\alpha}^{2}\phantom{\rule{1pt}{0ex}}({a}_{n}^{\prime \prime}-{k}^{2}{n}^{2}{a}_{n})=0\phantom{\rule{3.33333pt}{0ex}}.$$

(36)

this linear ODE is solved using its characteristic equation:

(37)

whose 4 solutions are:

$$\delta =\pm k\left|n\right|\phantom{\rule{2.em}{0ex}}\text{and}\phantom{\rule{2.em}{0ex}}\delta =\pm {\gamma}_{n}\phantom{\rule{1.em}{0ex}}\text{with}\phantom{\rule{3.33333pt}{0ex}}{\gamma}_{n}=\sqrt{{k}^{2}{n}^{2}+in{\alpha}^{2}}$$

(38)

The coefficients of the Fourier decomposition of *ψ*_{1} are determined by the boundary conditions. In real experiments carried out on fluid flow next to ciliated edges [11], the value of *k* is much larger than 1 (which means that the channel width along y is much larger than the metachronal wavelength). This enables us to translate the boundary condition at *y* = 1 (the last two of Eq 34) into the cancellation of the coefficients attached to the roots with positive real part in Eq 37. Finally, since the only non homogeneous terms appearing in these boundary conditions are in *e*^{iθ} and *e*^{−iθ} (sin *θ* and cos *θ*), we only retain the coefficients corresponding to *n* = 1 and *n* = −1, which leads to the following expression for *ψ*_{1}:

(39)

Since ${\gamma}_{-1}={\overline{\gamma}}_{1}$, imposing that *ψ*_{1} is real valued implies that the coefficients also satisfy ${A}_{-1}={\overline{A}}_{1}$ and ${B}_{-1}={\overline{B}}_{1}$. Therefore, *ψ*_{1} can be expressed as:

(40)

where c.c. stands for complex conjugate. The flow field at first order in *ε* is then:

$$\{\begin{array}{c}{u}_{1}={\psi}_{1,y}=-2\phantom{\rule{3.33333pt}{0ex}}\text{Re}\left\{\right[k{A}_{1}{e}^{-ky}+{\gamma}_{1}{B}_{1}{e}^{-{\gamma}_{1}y}\left]\phantom{\rule{3.33333pt}{0ex}}{e}^{i\theta}\right\}\hfill \\ {v}_{1}=-{\psi}_{1,x}=-2k\phantom{\rule{3.33333pt}{0ex}}\text{Im}\left\{\right[{A}_{1}{e}^{-ky}+{B}_{1}{e}^{-{\gamma}_{1}y}\left]\phantom{\rule{3.33333pt}{0ex}}{e}^{i\theta}\right\}\hfill \end{array}$$

(41)

*A*_{1} and *B*_{1} are determined from the two boundary conditions at *y* = 0 in Eq (34):

$$\{\begin{array}{cc}-k{A}_{1}-{\gamma}_{1}{B}_{1}-\varphi \phantom{\rule{1pt}{0ex}}({k}^{2}{A}_{1}+{\gamma}_{1}^{2}{B}_{1})\hfill & {\displaystyle =\frac{1}{2i}}\hfill \\ ik[{A}_{1}+{B}_{1}-\varphi \phantom{\rule{1pt}{0ex}}(-k{A}_{1}-{\gamma}_{1}{B}_{1}\left)\right]\hfill & {\displaystyle =-\frac{\beta}{2}}\hfill \end{array}$$

(42)

which finally gives:

$${A}_{1}=\frac{i}{2(1+k\varphi )}\phantom{\rule{1pt}{0ex}}\left(\frac{{\displaystyle \frac{\beta {\gamma}_{1}}{k}-1}}{{\gamma}_{1}-k}\right)\phantom{\rule{2.em}{0ex}}\text{and}\phantom{\rule{2.em}{0ex}}{B}_{1}=\frac{i}{2(1+{\gamma}_{1}\varphi )}\phantom{\rule{1pt}{0ex}}\left(\frac{1-\beta}{{\gamma}_{1}-k}\right)$$

(43)

*ψ*_{2}(*x*, *y*, *t*), the second order in *ε* of the stream function, is a periodic function of period 2*π*/*k* along *x* that solves the following system:

$$\{\begin{array}{c}{\Delta}^{2}{\psi}_{2}-{\alpha}^{2}\Delta {\psi}_{2,t}=2{\alpha}^{2}\{{\psi}_{1,y}\Delta {\psi}_{1,x}-{\psi}_{1,x}\Delta {\psi}_{1,y}\}\hfill \\ {\psi}_{2,y}-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{2,yy}=k(1+cos(2\theta \left)\right)-2\beta sin\theta \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,yy}+2\beta \varphi sin\theta \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,yyy}\phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=0\hfill \\ {\psi}_{2,x}-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{2,xy}=2\beta ksin\left(2\theta \right)-2\beta sin\theta \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,xy}+2\beta \varphi sin\theta \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,xyy}\phantom{\rule{2.em}{0ex}}\hfill & \text{at}\phantom{\rule{3.33333pt}{0ex}}y=0\hfill \end{array}$$

(44)

Since the first order *ψ*_{1} of the stream function contains only terms in *e*^{iθ} and *e*^{−iθ}, as seen in the previous section, a rapid analysis of the above boundary conditions shows that the second order *ψ*_{2} will only contain terms either constant or proportional to *e*^{2iθ} and *e*^{−2iθ}. We are interested here only in the constant term, which corresponds to the first non vanishing contribution to the steady flow field in the *ε*-expansion. From Eq 29, we know that the steady part of *ψ*_{2}, here called ${\psi}_{2}^{\left(s\right)}$, solves the following equation

$${\Delta}^{2}{\psi}_{2}^{\left(s\right)}=2{\alpha}^{2}\phantom{\rule{1pt}{0ex}}\phantom{\rule{1pt}{0ex}}{\left\{{\psi}_{1,y}\Delta {\psi}_{1,x}-{\psi}_{1,x}\Delta {\psi}_{1,y}\right\}}^{\left(s\right)}$$

(45)

where the superscript ‘(s)’ on the right hand side stands for the steady part. More precisely, the bracketed term can be rewritten as:

$$\begin{array}{ccc}{\psi}_{1,y}\Delta {\psi}_{1,x}-{\psi}_{1,x}\Delta {\psi}_{1,y}& =& {\psi}_{1,y}\left({\psi}_{1,xxx}+{\psi}_{1,xyy}\right)-{\psi}_{1,x}\left({\psi}_{1,xxy}+{\psi}_{1,yyy}\right)\hfill \\ & =& \left({\psi}_{1,y}{\psi}_{1,xyy}-{\psi}_{1,x}{\psi}_{1,yyy}\right)+\left({\psi}_{1,y}{\psi}_{1,xxx}-{\psi}_{1,x}{\psi}_{1,xxy}\right)\hfill \\ & =& {\left({\psi}_{1,y}{\psi}_{1,xy}-{\psi}_{1,x}{\psi}_{1,yy}\right)}_{y}+{\left({\psi}_{1,y}{\psi}_{1,xx}-{\psi}_{1,x}{\psi}_{1,xy}\right)}_{x}\hfill \end{array}$$

(46)

Let us examine the *x*-dependency of the parentheses appearing in the last line. All products contain terms that are either constant or oscillating in *x*. Therefore the contribution of the second parenthesis to the steady stream function can only be constant. Its *x*-derivative then vanishes and can be removed from the equation satisfied by ${\psi}_{2}^{\left(s\right)}$, leaving only the first parenthesis:

$${\psi}_{2,yyyy}^{(s)}=2{\alpha}^{2}\phantom{\rule{1pt}{0ex}}{({\psi}_{1,y}{\psi}_{1,yx}-{\psi}_{1,x}{\psi}_{1,yy})}_{y}^{(s)}$$

(47)

Integrating along *y* yields:

$${\psi}_{2,yyy}^{(s)}=2{\alpha}^{2}\phantom{\rule{1pt}{0ex}}{({\psi}_{1,y}{\psi}_{1,yx}-{\psi}_{1,x}{\psi}_{1,yy})}^{(s)}+2{A}_{2}$$

(48)

where 2*A*_{2} is an integration constant. Note that *ε*^{2}/2 × 2*A*_{2} = *ε*^{2}
*A*_{2} can actually be interpreted as the steady component of the pressure gradient ${p}_{x}^{\left(s\right)}$ which appears in the right hand side of the Navier-Stokes equation (Eq 15). The model therefore predicts a constant pressure gradient along *x*. From the previous section, we know that *ψ*_{1} can be written as:

(49)

Inserting this expression of *ψ*_{1} into the right hand side parenthesis of Eq 48 leaves:

$${({\psi}_{1,y}{\psi}_{1,yx}-{\psi}_{1,x}{\psi}_{1,yy})}^{(s)}=-ik|{{a}^{\prime}}_{1}{|}^{2}\phantom{\rule{1pt}{0ex}}-\phantom{\rule{1pt}{0ex}}ik{a}_{1}{{\overline{a}}^{\u2033}}_{1}+c.c.=2k\text{Im}({a}_{1}{{\overline{a}}^{\u2033}}_{1})$$

(50)

(The first term of the above right hand side being imaginary, its contribution vanishes when adding its complex conjugate). Eq 48 and the first boundary condition of Eq 44 thus provide the following simpler system for ${\psi}_{2}^{\left(s\right)}$ (the second boundary condition deals with *x*-derivatives that vanish in the steady flow due to the translation invariance of the problem):

$$\{\begin{array}{c}{\psi}_{2,yyy}^{\left(s\right)}=4{\alpha}^{2}k\phantom{\rule{3.33333pt}{0ex}}\text{Im}\left({a}_{1}{\overline{a}}_{1}^{\prime \prime}\right)+2{A}_{2}\hfill \\ {\psi}_{2,y}^{\left(s\right)}(x,0)-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{2,yy}^{\left(s\right)}(x,0)=k+C\hfill \end{array}$$

(51)

where *C* stands for the constant contribution of −2*β* sin *θ* (*ψ*_{1, yy} − *ϕ*
*ψ*_{1, yyy})|_{y = 0}. This constant *C* is determined using the expression of *ψ*_{1} in Eq 49:

$${\psi}_{1,yy}-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{1,yyy}=({a}_{1}^{\prime \prime}-\varphi \phantom{\rule{3.33333pt}{0ex}}{a}_{1}^{\prime \prime \prime})\phantom{\rule{3.33333pt}{0ex}}{e}^{i\theta}+c.c.$$

(52)

Therefore the constant contribution *C* can be expressed from the constants *A*_{1} and *B*_{1} computed in the previous section:

$$C=2\beta \phantom{\rule{3.33333pt}{0ex}}\text{Im}({a}_{1}^{\prime \prime}\left(0\right)-\varphi \phantom{\rule{3.33333pt}{0ex}}{a}_{1}^{\prime \prime \prime}\left(0\right))=2\beta \phantom{\rule{3.33333pt}{0ex}}\text{Im}({k}^{2}{A}_{1}+{\gamma}_{1}^{2}{B}_{1}+{k}^{3}\varphi {A}_{1}+{\gamma}_{1}^{3}\varphi {B}_{1})$$

(53)

Using the values of *A*_{1} and *B*_{1} obtained in the previous section gives immediately:

$$C=2\beta \phantom{\rule{3.33333pt}{0ex}}\text{Im}\phantom{\rule{1pt}{0ex}}\left[\frac{i{k}^{2}}{2}\phantom{\rule{1pt}{0ex}}\right(\frac{{\displaystyle \frac{\beta {\gamma}_{1}}{k}-1}}{{\gamma}_{1}-k})+\frac{i{\gamma}_{1}^{2}}{2}\phantom{\rule{1pt}{0ex}}(\frac{1-\beta}{{\gamma}_{1}-k}\left)\right]=\beta k+\beta (1-\beta )\phantom{\rule{3.33333pt}{0ex}}\text{Re}\left({\gamma}_{1}\right)$$

(54)

We now turn to Im(${a}_{1}{\overline{a}}_{1}^{\prime \prime}$) in order to determine ${\psi}_{2}^{\left(s\right)}$. Using Eq 40, this term can be rewritten:

$$\begin{array}{ccc}\hfill \text{Im}\left({a}_{1}{\overline{a}}_{1}^{\prime \prime}\right)& =\hfill & \text{Im}\phantom{\rule{1pt}{0ex}}\left\{\right({A}_{1}{e}^{-ky}+{B}_{1}{e}^{-{\gamma}_{1}y}\left)\right({k}^{2}{\overline{A}}_{1}{e}^{-ky}+{\overline{\gamma}}_{1}^{2}\overline{{B}_{1}}{e}^{-\overline{{\gamma}_{1}}y}\left)\right\}\hfill \\ & =& \text{Im}\phantom{\rule{1pt}{0ex}}\{{\overline{\gamma}}_{1}^{2}{A}_{1}\overline{{B}_{1}}{e}^{-(k+{\overline{\gamma}}_{1})y}+{k}^{2}{\overline{A}}_{1}{B}_{1}{e}^{-(k+{\gamma}_{1})y}+{\overline{\gamma}}_{1}^{2}{B}_{1}\overline{{B}_{1}}{e}^{-({\gamma}_{1}+{\overline{\gamma}}_{1})y}\}\hfill \end{array}$$

(55)

Using the relation ${\overline{\gamma}}_{1}^{2}={k}^{2}-i{\alpha}^{2}$ allows us to simplify the above expression to:

$$\begin{array}{ccc}\hfill \text{Im}\left({a}_{1}{\overline{a}}_{1}^{\prime \prime}\right)& =\hfill & \text{Im}\phantom{\rule{1pt}{0ex}}\{-i{\alpha}^{2}{A}_{1}\overline{{B}_{1}}{e}^{-(k+{\overline{\gamma}}_{1})y}-i{\alpha}^{2}{\left|{B}_{1}\right|}^{2}{e}^{-({\gamma}_{1}+{\overline{\gamma}}_{1})y}\}\hfill \\ & =& -{\alpha}^{2}\phantom{\rule{1pt}{0ex}}\left[\text{Re}\right({A}_{1}\overline{{B}_{1}}{e}^{-(k+{\overline{\gamma}}_{1})y})+{\left|{B}_{1}\right|}^{2}{e}^{-({\gamma}_{1}+{\overline{\gamma}}_{1})y}]\hfill \end{array}$$

(56)

This expression is now inserted into the first equation of Eq 51, which gives after 2 integrations:

$${\psi}_{2,y}^{\left(s\right)}={A}_{2}\phantom{\rule{1pt}{0ex}}{y}^{2}+{B}_{2}\phantom{\rule{1pt}{0ex}}y+{C}_{2}-4{\alpha}^{4}k\phantom{\rule{1pt}{0ex}}\phantom{\rule{1pt}{0ex}}\left[\text{Re}\phantom{\rule{1pt}{0ex}}\phantom{\rule{1pt}{0ex}}\left({A}_{1}{\overline{B}}_{1}\phantom{\rule{1pt}{0ex}}\frac{{e}^{-\left(k+{\overline{\gamma}}_{1}\right)y}}{{\left(k+{\overline{\gamma}}_{1}\right)}^{2}}\right)+{\left|{B}_{1}\right|}^{2}\phantom{\rule{1pt}{0ex}}\frac{{e}^{-\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)y}}{{\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)}^{2}}\right],$$

(57)

where *B*_{2} and *C*_{2} are integration constants. We observe here that the steady contribution to the horizontal fluid velocity appears as the sum of one parabolic and two exponential profiles. The parabolic profile depends on integration constants *A*_{2}, *B*_{2}, and *C*_{2} that are determined using the boundaries conditions. Velocity and shear stress vanish at *y* = 1, so that the fluid in the flowing layer matches the stagnant fluid above. The boundary conditions take the following form:

$$\{\begin{array}{cc}{\psi}_{2,y}^{\left(s\right)}\left(1\right)\hfill & =0={A}_{2}+{B}_{2}+{C}_{2}+{G}_{1}\hfill \\ {\psi}_{2,yy}^{\left(s\right)}\left(1\right)\hfill & =0=2{A}_{2}+{B}_{2}+{G}_{2}\hfill \\ {\psi}_{2,y}^{\left(s\right)}\left(0\right)-\varphi \phantom{\rule{3.33333pt}{0ex}}{\psi}_{2,yy}^{\left(s\right)}\left(0\right)\hfill & =k+C={C}_{2}-\varphi {B}_{2}+{G}_{3}+k+C\hfill \end{array}$$

(58)

where *G*_{1}, *G*_{2}, and *G*_{3} are defined such that:

$$\{\begin{array}{ccc}\frac{{G}_{1}}{4{\alpha}^{4}k}\hfill & =& -\text{Re}\phantom{\rule{1pt}{0ex}}\left({A}_{1}{\overline{B}}_{1}\frac{{e}^{-(k+{\overline{\gamma}}_{1})}}{{\left(k+{\overline{\gamma}}_{1}\right)}^{2}}\right)-|{B}_{1}{|}^{2}\frac{{e}^{-({\gamma}_{1}+{\overline{\gamma}}_{1})}}{{\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)}^{2}}\hfill \\ \frac{{G}_{2}}{4{\alpha}^{4}k}\hfill & =& \text{Re}\phantom{\rule{1pt}{0ex}}\left({A}_{1}{\overline{B}}_{1}\frac{{e}^{-(k+{\overline{\gamma}}_{1})}}{\left(k+{\overline{\gamma}}_{1}\right)}\right)+|{B}_{1}{|}^{2}\frac{{e}^{-({\gamma}_{1}+{\overline{\gamma}}_{1})}}{\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)}\hfill \\ \frac{{G}_{3}+k+C}{4{\alpha}^{4}k}\hfill & =& -\text{Re}\phantom{\rule{1pt}{0ex}}\left(\frac{{A}_{1}{\overline{B}}_{1}}{{\left(k+{\overline{\gamma}}_{1}\right)}^{2}}\right)-\varphi \text{Re}\phantom{\rule{1pt}{0ex}}\left(\frac{{A}_{1}{\overline{B}}_{1}}{\left(k+{\overline{\gamma}}_{1}\right)}\right)-\frac{|{B}_{1}{|}^{2}}{{\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)}^{2}}-\varphi \frac{|{B}_{1}{|}^{2}}{\left({\gamma}_{1}+{\overline{\gamma}}_{1}\right)}\hfill \end{array}$$

(59)

The coefficients of the steady parabolic profile are finally:

$$\{\begin{array}{cc}{A}_{2}\hfill & {\displaystyle =\frac{{G}_{1}-(1+\varphi ){G}_{2}-{G}_{3}}{1+2\varphi}}\hfill \\ {B}_{2}\hfill & {\displaystyle =\frac{-2{G}_{1}+{G}_{2}+2{G}_{3}}{1+2\varphi}}\hfill \\ {C}_{2}\hfill & {\displaystyle =\frac{-2\varphi \phantom{\rule{3.33333pt}{0ex}}{G}_{1}+\varphi \phantom{\rule{3.33333pt}{0ex}}{G}_{2}-{G}_{3}}{1+2\varphi}}\hfill \end{array}$$

(60)

In summary, the velocity field of the fluid flow can be expanded in powers of *ε*. The first non vanishing order in *ε* for the oscillatory part of the fluid flow is of order 1 while at the second order, a steady contribution appears. This steady flow is the one responsible for the micro-bead motion observed in the experiments. It has essentially a parabolic profile whose coefficients are directly determined from the values of the quantities *k*, *α*, *β*, and *ϕ*, which are in turn deduced from the fundamental parameters of the cilia motion and the surrounding fluid (frequency, amplitude, density, viscosity). In the experiments, fitting the measured micro-bead velocities with a parabolic profile will allow us to extract quantities that are not easily measurable, such as the cilium amplitude *a* or the transfer parameter *ϕ*.

In experiments carried out on real samples, the values of *k* are larger than 10 (in most cases much larger), and *α*^{2} is about 0.1. We can therefore safely assume that *e*^{−k} 1 and *α* 1. In this situation, *G*_{1} and *G*_{2} are negligible compared to *G*_{3}, and the coefficients of the parabolic profile are:

$${A}_{2}=-\frac{{G}_{3}}{1+2\varphi}\phantom{\rule{1.em}{0ex}},\phantom{\rule{1.em}{0ex}}{B}_{2}=\frac{2{G}_{3}}{1+2\varphi}\phantom{\rule{1.em}{0ex}},\phantom{\rule{1.em}{0ex}}{C}_{2}=-\frac{{G}_{3}}{1+2\varphi}$$

(61)

This leads to the following steady flow:

$$u\left(y\right)=-\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}(\frac{{G}_{3}}{1+2\varphi}\phantom{\rule{3.33333pt}{0ex}}{y}^{2}+\frac{2{G}_{3}}{1+2\varphi}\phantom{\rule{3.33333pt}{0ex}}y-\frac{{G}_{3}}{1+2\varphi})=-\frac{{\epsilon}^{2}}{2}\frac{{G}_{3}}{1+2\varphi}\phantom{\rule{3.33333pt}{0ex}}{(y-1)}^{2}$$

(62)

One recovers here a parabolic profile with vanishing velocity and vanishing shear stress at *y* = 1. Moreover, from Eq 59, the term proportional to 4*α*^{4}*k* in *G*_{3} can be neglected, leading to:

(63)

The dimensionless extrapolated velocity at the ciliated wall in this case is then:

$$u\left(0\right)=\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}\left(\frac{1+2\beta -{\beta}^{2}}{1+2\varphi}\right)k=\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}\left(\frac{2-{(\beta -1)}^{2}}{1+2\varphi}\right)\phantom{\rule{1pt}{0ex}}k$$

(64)

One can see here that the extrapolation of the flow velocity (hence of the microbead) at the ciliated wall depends on *β*, *ϕ*, and *k*. It has a maximum value, *ε*^{2}*k*/(1 + 2*ϕ*), which is achieved for *β* = 1, i.e., for a circular motion of the cilia tips. The case *β* = 0, together with a finite value of *ε* = *a*/*h* corresponds to the limit of the flat horizontal motion of the cilia tips. The dimensionless fluid velocity at the cilia wall in this situation is:

$$u\left(0\right)=\frac{{\epsilon}^{2}k}{2(1+2\varphi )}$$

(65)

One would think that the limit case *β* = 0 would lead to a vanishing steady velocity of the fluid, since all cilia undergo the same oscillatory motion, with a left/right symmetry. What we see here is that the metachronal wave (antiplectic, i.e. going backwards for *k* > 0) breaks this symmetry and leads to a positive steady contribution to the fluid velocity. When *k* = 0 (no wave), one recovers a pure oscillatory motion of the fluid, with no net fluid motion.

Finally, the limit *β* → +∞ corresponds to a flat vertical motion of the cilia tips. Since *ε* = *a*/*h* is the adimensioned horizontal beating amplitude, the limit has to be taken with *ε* → 0, *βε* remaining constant as the adimensioned vertical amplitude *b*/*h*. The horizontal component of the steady velocity is then:

$$u\left(0\right)=\underset{\beta \to +\infty}{lim}\frac{{\epsilon}^{2}}{2}\phantom{\rule{1pt}{0ex}}\left(\frac{1+2\beta -{\beta}^{2}}{1+2\varphi}\right)\phantom{\rule{1pt}{0ex}}k=-\frac{{b}^{2}k}{2{h}^{2}(1+2\varphi )}$$

(66)

Interestingly, the metachronal wave also breaks the left/right symmetry in this case, but in the opposite direction. In summary, the existence of an antiplectic metachronal wave triggers a net fluid motion along *Ox*, in the positive direction for smaller values of *β* and in the opposite direction for larger values of *β*, the critical value *β*_{c} being obtained when (*β*_{c} − 1)^{2} = 2 (according to Eq 64), i.e., for ${\beta}_{c}=1+\sqrt{2}\approx 2.41$. In realistic situations, the trajectories of cilia tips are very flat, so that one is always in the case of small values of *β*: the fluid steady motion occurs always in the direction opposite to the metachronal wave propagation.

Once the oscillatory velocity field $\overrightarrow{U}=(u,v)$ is computed, the trajectories and the crossing times of micro-beads in this field are calculated by solving their equation of motion:

$$m\phantom{\rule{3.33333pt}{0ex}}\frac{d{\overrightarrow{U}}_{b}^{*}}{d{t}^{*}}={\overrightarrow{F}}_{\text{drag}}\phantom{\rule{3.33333pt}{0ex}},$$

(67)

where *m* is the individual mass of the micro-bead, ${\overrightarrow{U}}_{b}^{*}$ its velocity, and ${\overrightarrow{F}}_{\text{drag}}$ is the drag force applied on the bead by the fluid flow. Assuming a spherical shape for the bead and small speed differences between the bead and the fluid, the drag force takes the Stokes expression

$${\overrightarrow{F}}_{\text{drag}}=-6\pi R\mu \phantom{\rule{1pt}{0ex}}({\overrightarrow{U}}_{b}^{*}-{\overrightarrow{U}}^{*})\phantom{\rule{3.33333pt}{0ex}},$$

(68)

where *R* stands for the bead diameter. Using the previously defined space and time units *h* and *ω*^{−1}, the adimensioned version of Eq 67 reads:

$$\frac{d{\overrightarrow{U}}_{b}}{dt}=\frac{\overrightarrow{U}-{\overrightarrow{U}}_{b}}{S{t}_{k}}\phantom{\rule{2.em}{0ex}}\text{with}\phantom{\rule{4.pt}{0ex}}\phantom{\rule{1.em}{0ex}}S{t}_{k}=\frac{m\omega}{6\pi R\mu}$$

(69)

*St*_{k} is the bead Stokes number, which characterizes the effective inertia of the bead in the fluid flow. For spherical particles of radius *R* and density *ρ*_{b}, this number also reads:

$$S{t}_{k}=\frac{{\displaystyle \frac{4}{3}\pi {R}^{3}{\rho}_{b}\omega}}{6\pi R\mu}=\frac{2}{9}\frac{{R}^{2}{\rho}_{b}\omega}{\mu}$$

(70)

The micro-beads are about 4.5 μm diameter, and made out of polystyrene of density *ρ*_{b} of order 1 g.cm^{−3}. At 10 Hz in water, the corresponding Stokes number is about 10^{−4}, which means that the micro-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}}\left({y}_{0}\right)=\frac{L}{\tau \left({y}_{0}\right)}=\frac{L}{{\oint}_{{\mathcal{T}}_{{y}_{0}}}dt}=\frac{L}{{\displaystyle {\oint}_{{\mathcal{T}}_{{y}_{0}}}\frac{ds}{\Vert \overrightarrow{U}\left(s\right)\Vert}}}\phantom{\rule{3.33333pt}{0ex}},$$

(71)

where *τ*(*y*_{0}) is the crossing time of the micro-bead entering at (0, *y*_{0}), and 𝒯_{y0} is the trajectory followed by this micro-bead. During each elementary step of this trajectory, the infinitesimal duration is $dt=ds/\Vert \overrightarrow{U}\left(s\right)\Vert $, $\overrightarrow{U}$ being the fluid velocity at curvilinear abscissa *s* of the trajectory. The effective speed *V*_{eff} corresponds to the quantity measured in our experiments.

The fluid velocity field is solved using the Fourier transform decomposition exposed earlier. This field is periodic both in space and time. Fig 3A and 3B displays the horizontal and vertical components of the fluid velocity, respectively, while Fig 3D presents the micro-bead trajectories in this fluid, simulated by numerically solving Eq 69. One can observe that the beads follow slightly wavy trajectories when close to the cilia wall, but that these trajectories become almost perfect straight lines when moving away from the wall farther than a cilia length. This means that the transit time of a micro-bead across the simulation window is dominated by the steady part of the flow field.

Fig 4A displays the various contributions to the flow profile in the *y*-direction for typical values of the input parameters CBF, CBA, λ, *ϕ*, and *h*. One observes that the total bead velocity (blue line) is essentially dominated by the steady parabolic contribution (black line) determined by the coefficients *A*_{2}, *B*_{2}, and *C*_{2} in Eq 57. The oscillatory part coming from the first order in *ε* (green line) brings a significant contribution only very close to the ciliated edge, and the exponential contribution in the steady part is negligible everywhere. Consequently, measuring the bead velocity as function of the distance to the ciliated edge essentially amounts to measuring the steady parabolic profile of the flow. Fie 4B shows the excellent agreement between the model and the actual measurements performed on different ciliated edges (here 3 samples from 3 different subjects are presented).

The horizontal component of the steady pressure gradient ${p}_{x}^{\left(s\right)}$ (see Eq 48), applies a force on the fluid directed negatively along O*x*. This force is exactly compensated by the force exerted by the cilia on the fluid, proportional to the dimensionless shear stress *u*/*y* (see S1 File for the detailed force balance between the cilia and the fluid). The dimensionless steady force applied by the cilia to a volume fluid of length 2*π*/*k* (the dimensionless wavelength) in the O*x* direction and dimensionless thickness *H*/*h* in the O*z* direction thus reads:

$${F}_{w}^{\left(s\right)}=-\frac{H}{h}\frac{2\pi}{k}\frac{\partial {u}^{\left(s\right)}}{\partial y}{|}_{(y=0)}=-\frac{\lambda H}{{h}^{2}}\frac{\partial {u}^{\left(s\right)}}{\partial y}{|}_{(y=0)}$$

(72)

Using the parabolic velocity profile *u*^{(s)}(*y*) = *u*(0)(*y* − 1)^{2} found in Eq 62, and putting back dimensional quantities finally leads to the expression of the local shear stress *τ*_{w} applied to the fluid by the cilia:

$${\tau}_{w}=\frac{\left(\mu \omega {h}^{2}\right){F}_{w}^{\left(s\right)}}{\lambda H}=\frac{\mu \omega {h}^{2}}{\lambda H}\frac{\lambda H}{{h}^{2}}2u\left(0\right)=\frac{2\mu {U}_{0}}{h}$$

(73)

where *U*_{0} is the velocity extrapolated at *y* = 0 from the measurements of microbead velocities above the cilia. This shear stress characterizes the momentum transfer between cilia and the surrounding fluid. Since *h* is also directly measured by fitting the microbead velocity profile with a parabolic profile, it means that *τ*_{w} can be directly deduced using this microbead tracking technique. Consequently, we propose this shear stress as an index for assessing the efficiency of the ciliary beating. One has to stress that we do not intend to reproduce in this model the *in vivo* condition. The shear stress measured using this micro-bead velocity technique is not assumed to be identical to the shear stress experienced by mucus in the pulmonary airways. It is only a way to assess of the ability of the ciliated edge to transfer momentum into the surrounding fluid, thus defining a usable clinical index.

Finally, the extrapolated velocity at the wall can be compared to the one predicted from Eq 64:

$${U}_{0}=h\omega \phantom{\rule{3.33333pt}{0ex}}u\left(0\right)=h\omega \frac{{a}^{2}}{2{h}^{2}}\phantom{\rule{1pt}{0ex}}\left(\frac{1+2\beta -{\beta}^{2}}{1+2\varphi}\right)\phantom{\rule{1pt}{0ex}}\frac{2\pi h}{\lambda}=\frac{\pi {a}^{2}\omega}{\lambda}\phantom{\rule{1pt}{0ex}}\left(\frac{1+2\beta -{\beta}^{2}}{{\displaystyle 1+2\frac{{\varphi}^{*}}{h}}}\right)$$

(74)

The cilium beating amplitude *a* (CBF), the the cilium beating frequency *ω*/2*π* (CBF), and the metachronal wavelength λ are measured directly by microscopic measurements on the cilia, while *h* is extracted from the parabolic fitting of the microbead velocity profile. In the approximation of a flat beating (*β* = 0), this set of measurements therefore provides a way to assess *ϕ**, the equivalent sliding length introduced in the boundary condition of Eq 18.

In human airways, the coordinated motion of the cilia covering the epithelium induces a complex displacement of a double layer consisting in a bottom periciliary layer (PCL) and a top mucus layer. Both layers have similar thicknesses, about 5 to 10 μm, but very different viscosities. The situation presented in this paper is however very different from the *in vivo* condition. It reproduces in fact the *ex vivo* experiments performed on a few ciliated cells obtained from nasal brushing. In these experiments, clusters formed of a few ciliated cells are immersed in water, and observed between two microscope slides by high speed video-microscopy. The ciliated edges are therefore lying horizontally in the simulation plane (O*x*, O*y*), which is perpendicular to the optical axis O*z* of the microscope.

The first and foremost difference lies in the fact that the fluid surrounding the cilia is now essentially water, i.e., a homogeneous Newtonian fluid. Secondly, unlike the in vivo situation where mucus is propelled only by the forward trajectory of the cilia tips, here both forward and backward parts of the stroke cycle contribute to the fluid motion. Consequently, as expressed through the mathematical model developed here, the oscillatory motion of the cilia induces at first order in the amplitude *a* an oscillatory motion of the fluid, and at second order a net steady component that pushes the fluid in the direction of the cilia beating. The net motion of the fluid originates from the phase shift between neighboring cilia tips, which breaks the left/right symmetry. More precisely, it creates an unbalance on the envelope between the respective densities of tips moving forward and those moving backwards. This unbalance appears in the second term of Eq 9. One can also note that the direction of the net fluid motion does not depend on the sense of rotation of the tips (clockwise or counterclockwise), but only the direction of the metachronal wave: the fluid is propelled in the direction opposite to the metachronal wave. This property due to the antiplectic nature of the wave might have an interesting consequence: Although the forcing mechanism on the fluid is very different in the presence of mucus (the tips of the cilia enter the mucus layer and push it forward), the fluid covering the cilia would still be pushed in the same direction if the mucus viscosity was significantly lowered or if the mucus would disappear.

In the *ex vivo* experiments, the Womersley number *α* is close to 0.4 while the dimensionless wave vector *k* = 2*πh*/λ is at least about 10 and most of the time much larger than that (*h* is comprised between 30 μm and 140 μm in the experiments, see [11]). This means that the contribution of the convective acceleration to momentum conservation is very small in all cases, and that in first approximation, the system can be understood using a linear Stokes equation point of view. The steady and oscillatory contributions to the wall velocity decouple and induce independent and additive steady and oscillatory contributions to the flow field, respectively. Regarding micro-bead velocity far from the ciliated edge, the system then behaves as an equivalent stationary system with an effective wall velocity *U*_{w} = *ε*^{2}*k*/2.

At a distance larger than a fraction of the metachronal wavelength from the ciliated edge (the wavelength λ is of order 10–20 μm in all experiments, see [11]), the oscillatory behavior of the flow field vanishes and only the net steady horizontal flow remains. The microbead tracking technique therefore probes this steady part of the flow field, which exhibits a parabolic profile along *y*. The steady contribution to the dimensionless horizontal wall velocity is *ε*^{2}*k*/2 (see Eq 11). In the case of flat beating (*β* = 0) and no slip boundary condition at the cilia wall (*ϕ* = 0), this value is exactly the steady fluid velocity at the wall expressed in Eq 65. As a consequence, it also corresponds to the extrapolated micro-bead velocity at the wall. Our model shows that the synchronized elliptic motion of the cilia generates a correction to this velocity by a factor (1 + 2*β* − *β*^{2}) (Eq 64). In addition, an imperfect momentum transfer between the cilia and the fluid (*ϕ* > 0) reduces this velocity.

Introducing a partial slip at the cilia wall, through a Robin type boundary condition with an effective sliding length *ϕ** = *hϕ* (Eq 18), reduces the fluid velocity by a factor (1 + 2*ϕ*), a result that is also obtained in a purely steady model for a parabolic profile. As shown in the companion paper [11], this effective sliding length is directly correlated to the average density of cilia. This is very consistent with the use of this type of boundary condition in fluid flow next to porous flow, the pore openings corresponding to the empty space between cilia. In clinical application of this model, it is expected that *ϕ* will not be a fitting parameter but predicted by the mere measurement of the cilia density.

The model proposed in this article relies on several assumptions, hence has a few limitations that we examine now.

First, the motion of each individual cilium tip is assumed to be elliptic. Although this is a generally accepted hypothesis, microscopic imaging of the actual motion of the cilium shows a more complicated trajectory [34]. In particular, the path followed by the tip on the backward trajectory seems to be closer to the forward path than for an elliptic motion. However, we have seen that the net steady motion of the fluid is essentially generated by the metachronal wave. This discrepancy of the trajectory should be accounted for through a change of the parameter *β* representing the ellipse eccentricity. One also has to stress that, since our model reproduces an experimental setup in which cilia are surrounded only by water, one might expect a different motion from the one achieved in *in vivo* situations where cilia are beating in a periciliary liquid while their tips are entering the bottom of the mucus layer.

The model also assumes an exact synchronization of the individual cilia motions, generating a metachronal wave of constant velocity. This implies that the phase shift between two cilia is a linear function of the distance separating them along the direction of the metachronal wave propagation. *Ex vivo* experiments carried out in [11] show that this hypothesis is satisfactorily verified at the observational scale of the experiment, i.e., an edge of a cluster containing a few ciliated cells. A deviation from the linear phase shift would disrupt the translational periodicity of the system and alter our solution based on a Fourier decomposition along the horizontal axis. The fluid velocity at the cilia wall computed in our model can thus be considered as an “ideal velocity”, reached for a perfect metachronal wave. Any measured discrepancy with respect to this ideal velocity could be therefore be attributed either to an imperfect momentum transfer (through the effective slip length *ϕ*), or to a perturbed metachronal wave.

In the model, the ciliated edge is supposed to be flat. This assumption is generally valid in real experiments (see [11]). Cases in which cell clusters appear to have a rough, curved, or disrupted ciliated edge are removed from the measurement procedure.

An important assumption lies in the fact that the fluid surrounding the cilia is stagnant at a distance *h* above the edge. The stagnation of the fluid above the cilia is observed experimentally, and originates from the external environment of the measured cell cluster. Various explanations can account for it: it can be due to the presence of other cell clusters at a few hundred microns distance which create a complex fluid flow in the entire system, or by the friction of the fluid on the upper and lower slides, or the appearance of boundary layers. Determining this value would require modeling the full 3D geometry of the system whose geometry is not easily accessible by microscopic observation. The distance *h* therefore is treated as the only fitting parameter of the model (the sliding length *ϕ* being calculated from the measured cilia density).

Finally, the bead velocities are found in the model to be almost always parallel to the cilia wall. In real experiments, a vertical component (i.e., in the direction perpendicular to the wall) might appear, although much smaller than the horizontal component. This is in particular the case when the ciliated edge is not strictly flat, or when the observation window is such that the ciliated edge is on one side of the window and the bead passing on the other side of the same window. The micro-beads are in this situation not set into motion only by the considered ciliated edge, and external influence originating from outside the observation window interfere, contributing to create a more complex fluid flow pattern. Such situations should be discarded from the analysis, either through eye examination or by an automatized procedure.

We have developed a mathematical model of micro-bead velocity in a fluid set into motion by the periodic beating of a ciliated edge. The cilia wall is represented as a continuous envelope whose motion is calculated from the coordinated movement of all cilia. The boundary condition imposed by the effective wall induces a motion of the fluid above the wall. This motion has two components: one oscillatory, at the spatial and temporal periodicity of the metachronal wave and the cilia beating, respectively, and one steady, oriented along the direction of the wall. The oscillatory component vanishes at a distance of the order of the metachronal wavelength. The steady component extends on a much longer distance and exhibits a parabolic profile in the direction perpendicular to the wall. The parameters governing this profile are the fluid viscosity, the ciliary beating frequency (CBF), the ciliary beating amplitude (CBA), the metachronal wavelength (λ), the cilia density (through the sliding length *ϕ**), and the distance from the wall to the region of stagnant fluid (*h*). Polystyrene micro-beads immersed in the fluid act as massless tracers, allowing the measurement of the local fluid velocity, hence the measurement of the velocity profile. All aforementioned parameters can be determined through local measurement by high speed video-microscopy, except for the distance *h* which remains the only fitting parameter of the model extracted from the velocity profile. The velocity extrapolated from the profile at the wall is found to be linearly related to the shear stress exerted by the cilia on the fluid. This shear stress is proposed as a new index for assessing the efficiency of the ciliary beating. One has to stress that this index can be measured from nasal brushing in the clinical setting without any modification of the current clinical procedures. Preliminary tests (see [11]) have already shown that this index has the potential to be a powerful screening test, able to distinguish patients suffering from various alterations of the cilia beating such Primary Ciliary Dyskinesia (PCD).

Derivation of shear stress exerted by the cilia motion onto the surrounding fluid.

(PDF)

Click here for additional data file.^{(111K, pdf)}

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 funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

All relevant data are within the paper and its Supporting Information files.

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

2.
Wanner A, Salathe M, O’Riordan TG. Mucociliary clearance in the airways. Am J Respir Crit Care Med. 1996;154:1868–1902. doi: 10.1164/ajrccm.154.6.8970383
[PubMed]

3.
Knowles MR, Boucher RC. Mucus clearance as a primary innate defense mechanism for mammalian airways. J Clin Invest. 2002;109:571–577. doi: 10.1172/JCI15217
[PMC free article] [PubMed]

4.
Fahy JV, Dickey BF. Airway mucus function and dysfunction. N Engl J Med. 2010;363(23):2233–47. doi: 10.1056/NEJMra0910061
[PMC free article] [PubMed]

5.
Sanderson MJ, Sleigh MA. Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony. J Cell Sci. 1981;47:331–47. [PubMed]

6.
Sleigh MA, Blake JR, Liron N. The propulsion of mucus by cilia. Am Rev Respir Dis. 1988;137(13):726–41. doi: 10.1164/ajrccm/137.3.726
[PubMed]

7.
Gheber L, Priel Z. On metachronism in ciliary systems: a model describing the dependence of the metachronal wave properties on the intrinsic ciliary parameters. Cell Motil Cytoskeleton. 1990;16(3):167–81. doi: 10.1002/cm.970160304
[PubMed]

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. doi: 10.1016/S0006-3495(97)78686-7
[PubMed]

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. doi: 10.1016/j.bpj.2009.09.048
[PubMed]

10.
Salathe M. Regulation of mammalian ciliary beating. Annu Rev Physiol. 2007;69:401–22. doi: 10.1146/annurev.physiol.69.040705.141253
[PubMed]

11.
Bottier M, Blanchon S, Pelle G, Bequignon E, Isabey D, Coste A, et al.
A New Index for Characterizing Micro-bead Motion in a Flow Induced by Ciliary Beating: part I, Experimental Analysis. PLoS Comput Biol. 2017. [PMC free article] [PubMed]

12.
Smith DJ, Gaffney EA, Blake JR. Modelling mucociliary clearance. Resp Physiol Neurobi. 2008;163:
178–188. doi: 10.1016/j.resp.2008.03.006 [PubMed]

13.
Blake J. A model for the micro-structure in ciliated organisms. J Fluid Mech. 1972;55
doi: 10.1017/S0022112072001612

14.
Lauga E, Powers TR. The hydrodynamics of swimming microorganisms. Rep Prog Phys. 2009;72(9):096601
doi: 10.1088/0034-4885/72/9/096601

15.
Lee WL, Jayathilake PG, Tan Z, Le DV, Lee HP, Khoo BC. Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density. Comput Fluid. 2011;49(1):214–221. doi: 10.1016/j.compfluid.2011.05.016

16.
Mitran SM. Metachronal wave formation in a model of pulmonary cilia. Comput Struct. 2007;85:763–774. doi: 10.1016/j.compstruc.2007.01.015
[PMC free article] [PubMed]

17.
Yang X, Dillon RH, Fauci LJ. An Integrative Computational Model of Multiciliary Beating. Bull Math Biol. 2008;70(4):1192–1215. doi: 10.1007/s11538-008-9296-3
[PubMed]

18.
Keller SR. Fluid mechanical investigations of ciliary propulsion. California Institute of Technology; 1975.

19.
Blake JR, Winet H. On the mechanics of muco-ciliary transport. Biorheology. 1980;17:125–34. [PubMed]

20.
Smith DJ, Gaffney EA, Blake JR. Discrete Cilia Modelling with Singularity Distributions: Application to the Embryonic Node and the Airway Surface Liquid. Bull Math Biol. 2007;69(5):1477–1510. doi: 10.1007/s11538-006-9172-y
[PubMed]

21.
Childress S. Mechanics of Swimming and Flying. Cambridge University Press; 1981.
Cambridge Books

22.
Ross SM, Corrsin S. Results of an analytical model of mucociliary pumping. J Appl Physiol. 1974;34(3):333–40. [PubMed]

23.
Velez-Cordero JR, Lauga E. Waving transport and propulsion in a generalized Newtonian fluid. J Non-Newton Fluid. 2013;199:37–50. doi: 10.1016/j.jnnfm.2013.05.006

24.
Shen XN, Arratia PE. Undulatory Swimming in Viscoelastic Fluids. Phys Rev Lett. 2011;106:208101
doi: 10.1103/PhysRevLett.106.208101
[PubMed]

25.
Liu B, Powers TR, Breuer KS. Force-free swimming of a model helical flagellum in viscoelastic fluids. Proc Natl Acad Sci USA. 2011;108:19516–19520. doi: 10.1073/pnas.1113082108
[PubMed]

26.
Dasgupta M, Liu B, Fu HC, Berhanu M, Breuer KS, Powers TR, et al.
Speed of a swimming sheet in Newtonian and viscoelastic fluids. Phys Rev E. 2013;87:013015
doi: 10.1103/PhysRevE.87.013015 [PubMed]

27.
Lauga E. Propulsion in a viscoelastic fluid. Phys Fluids. 2007;19
doi: 10.1063/1.2751388

28.
Fu H, Wolgemuth W, Powers TR. Swimming speeds of filaments in nonlinearly viscoelastic fluids. Phys Fluids. 2009;21
doi: 10.1063/1.3086320 [PubMed]

29.
Teran J, Fauci L, Shelley M. Viscoelastic Fluid Response Can Increase the Speed and Efficiency of a Free Swimmer. Phys Rev Lett. 2010;104:038101
doi: 10.1103/PhysRevLett.104.038101
[PubMed]

30.
Zhu L, Do-Quang M, Lauga E, Brandt L. Locomotion by tangential deformation in a polymeric fluid. Phys Rev E. 2011;83:011901
doi: 10.1103/PhysRevE.83.011901 [PubMed]

31.
Zhu L, Lauga E, Brandt L. Self-propulsion in viscoelastic fluids: Pushers vs. pullers. Phys Fluids. 2012;24(5). doi: 10.1063/1.4718446

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

33.
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).

34.
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. doi: 10.1136/thorax.55.4.314
[PMC free article] [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |