Search tips
Search criteria 


Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS Comput Biol. 2017 July; 13(7): e1005552.
Published online 2017 July 14. doi:  10.1371/journal.pcbi.1005552
PMCID: PMC5510810

A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part II, modeling

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 & editing1,2,3,5,*
Mark Alber, Editor


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.

Author summary

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) [14]. 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 [57].

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/μm2 [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].

Fig 1
Schematic representation of the stroke of an individual cilium and the envelope model.

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 [1820]. 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 [2426] and numerically [2731] 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.

Materials and methods

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

From the individual cilium motion to the metachronal wave

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 (Xw*,Yw*) are


where β is the ellipse eccentricity, 2a 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.

Fig 2
Schematic elliptic motion of an individual ciliary tip.

To reproduce the backwards traveling metachronal wave of wavelength λ, we introduce in the periodic motion of each cilium a phase shift 2πξ*λ which linearly depends on the cilium position:


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 (Xw, Yw) are introduced:


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


In the following, we assume that the ciliary beating amplitude is much smaller than the film thickness, i.e. ε [double less-than sign] 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:


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


Since both Yw and (ξx) are first order in ε (Eqs 4 and 5), the expansion to the second order in ε of yw is

yw(xt) = εβsin(kxt) + ε2βkcos2(kxt)

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:


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


In summary, at location x in the Eulerian frame, the second order expansion in ε of the vertical position yw and of the velocity vector (uw, vw) of the ciliated wall are:



vw(xt) = εβcosθε2βksin(2θ)

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 Uw,1=(uw,1,vw,1) and Uw,2=(uw,2,vw,2), respectively, defined such that:


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


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 uw,2. This contribution is proportional to k and will be the origin of the steady horizontal motion of the fluid above the cilia.

Computing the flow field

We now compute the oscillatory flow field of a fluid of density ρ and viscosity μ in the channel comprised between y = yw(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 and μω, respectively. Hence, the dimensionless Navier-Stokes equation reads


where U=(u,v) is the dimensionless velocity field, p is the dimensionless pressure, and α is the Womersley number defined by α2=ρh2ωμ. 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−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

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):


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:


Next to the cilia envelope (y = yw), 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]:


If expressed in dimensional quantities, the parameter ϕ becomes a length ϕ* = . 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 U in powers of ε, U=εU1+ε22U2, the boundary conditions can be expressed at all orders in ε. Since yw is first order in ε (see Eq 10), the boundary condition can be expanded around y = 0 both for the velocity U and its normal derivative U/y:



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:



The stream function

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


Such a solution automatically fulfills the continuity equation div(U)=0. The dimensionless Navier-Stokes equation now reads:


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:

α2 {ψxxtψyytψyψyyxψxψyyyψyψxxxψxψxxy} = ψyyyy + 2ψxxψyyψxxxx

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

α2ψtψyΔψxψxΔψy} = Δ2ψ

Asymptotic expansion of the stream function

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


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:

α2Δψ1,t = Δ2ψ1

while the second order is:

α2ψ2,t + 2ψ1,yΔψ1,x - 2ψ1,xΔψ1,y} = Δ2ψ2

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:

ψ1,y|(x,0)ϕ ψ1,yy = sinθ

ψ1,x|(x,0)ϕ ψ1,xy =  - βcosθ

ψ2,y|(x,0)ϕ ψ2,yyk(1 + cos(2θ)) - 2βsinθ (ψ1,yyϕ ψ1,yyy)

ψ2,x|(x,0)ϕ ψ2,xy = 2βksin(2θ) - 2βsinθ (ψ1,xyϕ ψ1,xyy)

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.

First order of the stream function

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


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


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


this linear ODE is solved using its characteristic equation:

δ4 - (2k2n2inα2) δ2k2n2 (k2n2inα2) = 0

whose 4 solutions are:


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 and e (sin θ and cos θ), we only retain the coefficients corresponding to n = 1 and n = −1, which leads to the following expression for ψ1:

ψ1(xyt) = [A1e-kyB1e-γ1y] eiθ + [A-1e-kyB-1e-γ-1y] e-iθ

Since γ-1=γ¯1, imposing that ψ1 is real valued implies that the coefficients also satisfy A-1=A¯1 and B-1=B¯1. Therefore, ψ1 can be expressed as:

ψ1(xyt) = [A1e-kyB1e-γ1y] eiθc.c. , 

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


A1 and B1 are determined from the two boundary conditions at y = 0 in Eq (34):


which finally gives:


Second order of the stream function: The steady contribution

ψ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:


Since the first order ψ1 of the stream function contains only terms in e and e, 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 e2 and e−2. 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 ψ2(s), solves the following equation


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


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 ψ2(s), leaving only the first parenthesis:


Integrating along y yields:


where 2A2 is an integration constant. Note that ε2/2 × 2A2 = ε2 A2 can actually be interpreted as the steady component of the pressure gradient px(s) 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:

ψ1(xyt) = a1(y) eiθc.c. with a1(y) = A1e-kyB1e-γ1y

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

(ψ1,yψ1,yxψ1,xψ1,yy)(s)=ik|a1|2ika1a¯1+c.c.=2k Im(a1a¯1)

(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 ψ2(s) (the second boundary condition deals with x-derivatives that vanish in the steady flow due to the translation invariance of the problem):


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:


Therefore the constant contribution C can be expressed from the constants A1 and B1 computed in the previous section:


Using the values of A1 and B1 obtained in the previous section gives immediately:


We now turn to Im(a1a¯1) in order to determine ψ2(s). Using Eq 40, this term can be rewritten:


Using the relation γ¯12=k2-iα2 allows us to simplify the above expression to:


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


where B2 and C2 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 A2, B2, and C2 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:


where G1, G2, and G3 are defined such that:

{G14α4k=Re(A1B¯1 e(k+γ¯1)(k+γ¯1)2)|B1|2 e(γ1+γ¯1)(γ1+γ¯1)2G24α4k=Re(A1B¯1 e(k+γ¯1)(k+γ¯1))+|B1|2 e(γ1+γ¯1)(γ1+γ¯1)G3+k+C4α4k=Re(A1B¯1(k+γ¯1)2)ϕ Re(A1B¯1(k+γ¯1))|B1|2(γ1+γ¯1)2ϕ|B1|2(γ1+γ¯1)

The coefficients of the steady parabolic profile are finally:


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 ϕ.

The steady velocity profile

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 ek [double less-than sign] 1 and α [double less-than sign] 1. In this situation, G1 and G2 are negligible compared to G3, and the coefficients of the parabolic profile are:


This leads to the following steady flow:


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α4k in G3 can be neglected, leading to:

G3 ≈  - kC =  - kβkβ(1 - β) Re(γ1) ≈ (β2 - 2β - 1)k

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


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, ε2k/(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:


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:


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 βc=1+22.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.

Modeling the micro-bead motion

Once the oscillatory velocity field 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:


where m is the individual mass of the micro-bead, Ub* its velocity, and Fdrag 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


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:


Stk 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:


The micro-beads are about 4.5 μm diameter, and made out of polystyrene of density ρb of order 1−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 y0, the effective speed is computed as:


where τ(y0) is the crossing time of the micro-bead entering at (0, y0), and 𝒯y0 is the trajectory followed by this micro-bead. During each elementary step of this trajectory, the infinitesimal duration is dt=ds/U(s), U being the fluid velocity at curvilinear abscissa s of the trajectory. The effective speed Veff corresponds to the quantity measured in our experiments.


Simulation of micro-bead motions

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 3
Numerical simulation of the velocity field and micro-bead trajectories.

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 A2, B2, and C2 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).

Fig 4
Parabolic velocity profile.

The ciliary beating efficiency index

The horizontal component of the steady pressure gradient px(s) (see Eq 48), applies a force on the fluid directed negatively along Ox. This force is exactly compensated by the force exerted by the cilia on the fluid, proportional to the dimensionless shear stress [partial differential]u/[partial differential]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 Ox direction and dimensionless thickness H/h in the Oz direction thus reads:


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:


where U0 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:


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 (Ox, Oy), which is perpendicular to the optical axis Oz 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 Uw = ε2k/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 ε2k/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 ϕ* = (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.

Model limitations

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

Supporting information

S1 File

Force balance on a fluid volume.

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


Funding Statement

This work has been supported by a grant from Fondation pour la Recherche Médicale (, 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 ( 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

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