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 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/μ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]. 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.
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 {X w * = ξ * - a cos (ω t *) Y w * = β a sin (ω t *) (1) 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. 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: {X w * = ξ * - a cos (ω t * + 2 π ξ * λ) Y w * = β a sin (ω t * + 2 π ξ * λ) (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 (Xw, Yw) are introduced: ε = a h, k = 2 π h λ, X w = X w * h, Y w = Y w * h, x = x * h, y = y * h, ξ = ξ * h, t = ω t * (3) The motion equations of each tip are thus rewritten in a dimensionless form: {X w (ξ, t) = ξ - ε cos (k ξ + t) Y w (ξ, t) = β ε sin (k ξ + t) (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: x = ξ - ε cos (k ξ + t) (5) 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: y w (x, t) = Y w (ξ, t) = Y w (x, t) + (ξ - x) ∂ Y w ∂ ξ | ξ = x + (ξ - x) 2 2 ∂ 2 Y w ∂ ξ 2 | ξ = x + ⋯ (6) Since both Yw and (ξ − x) are first order in ε (Eqs 4 and 5), the expansion to the second order in ε of yw is y w (x, t) = ε β sin (k x + t) + ε 2 β k cos 2 (k x + t) (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: { uw (x, t) = (∂Xw∂t) ξ=εsin (kξ+t) vw (x, t) =∂yw∂t=εβcos (kx+t) −ε2βksin (2 (kx+t) ) (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 u w (x, t) = ∂ X w ∂ t | ξ = x + (ξ - x) ∂ u w ∂ ξ | ξ = x + … = ε sin (k x + t) + [ ε cos (k x + t) ] [ ε k cos (k x + t) ] = ε sin (k x + t) + ε 2 k cos 2 (k x + t) (9) 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: y w (x, t) = ε β sin θ + ε 2 β k 2 (1 + cos (2 θ) ) (10) u w (x, t) = ε sin θ + ε 2 k 2 (1 + cos (2 θ) ) (11) v w (x, t) = ε β cos θ - ε 2 β k sin (2 θ) (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 U → w, 1 = (u w, 1, v w, 1) and U → w, 2 = (u w, 2, v w, 2), respectively, defined such that: U → w = ε U → w, 1 + ε 2 2 U → w, 2 (13) In summary, the first and second orders in ε of the location and velocities of the cilia wall are: {y w, 1 (θ) = β sin θ = β 2 i e i θ - β 2 i e - i θ y w, 2 (θ) = β k (1 + cos (2 θ) ) = β k + β k 2 e 2 i θ + β k 2 e - 2 i θ u w, 1 (θ) = sin θ = 1 2 i e i θ - 1 2 i e - i θ u w, 2 (θ) = k (1 + cos (2 θ) ) = k + k 2 e 2 i θ + k 2 e - 2 i θ v w, 1 (θ) = β cos θ = β 2 e i θ + β 2 e - i θ v w, 2 (θ) = - 2 β k sin (2 θ) = i β k e 2 i θ - i β k e - 2 i θ (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 uw, 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 = 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 hω and μω, respectively. Hence, the dimensionless Navier-Stokes equation reads α 2 (∂ U → ∂ t + (U → · ∇ →) U →) = - ∇ → p + Δ U →, (15) where U → = (u, v) is the dimensionless velocity field, p is the dimensionless pressure, and α is the Womersley number defined by α 2 = ρ h 2 ω μ. 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. Since the flow is two-dimensional and incompressible, a natural formulation of the problem is obtained by introducing the stream function ψ such that: u = ∂ ψ ∂ y and v = - ∂ ψ ∂ x (23) Such a solution automatically fulfills the continuity equation div (U →) = 0. The dimensionless Navier-Stokes equation now reads: {α 2 (ψ y t + ψ y ψ y x - ψ x ψ y y) = - p x + ψ y x x + ψ y y y α 2 (- ψ x t + ψ y ψ x x + ψ x ψ x y) = - p y - ψ x x x - ψ x y y (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: α 2 {ψ x x t + ψ y y t + ψ y ψ y y x - ψ x ψ y y y + ψ y ψ x x x - ψ x ψ x x y} = ψ y y y y + 2 ψ x x ψ y y + ψ x x x x (25) This equation can finally be rewritten in the following compact form: α 2 {Δ ψ t + ψ y Δ ψ x - ψ x Δ ψ y} = Δ 2 ψ (26) 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: m d U → b * d t * = F → drag, (67) where m is the individual mass of the micro-bead, U → b * its velocity, and F → 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 F → drag = - 6 π R μ (U → b * - U → *), (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: d U → b d t = U → - U → b S t k with S t k = m ω 6 π R μ (69) 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: S t k = 4 3 π R 3 ρ b ω 6 π R μ = 2 9 R 2 ρ b ω μ (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 y0, the effective speed is computed as: where τ (y0) is the crossing time of the micro-bead entering at (0, y0), and T y 0 is the trajectory followed by this micro-bead. During each elementary step of this trajectory, the infinitesimal duration is d t = d s / ‖ U → (s) ‖, 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 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). The horizontal component of the steady pressure gradient p x (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 ∂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 Ox direction and dimensionless thickness H/h in the Oz direction thus reads: F w (s) = - H h 2 π k ∂ u (s) ∂ y | (y = 0) = - λ H h 2 ∂ u (s) ∂ 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: τ w = (μ ω h 2) F w (s) λ H = μ ω h 2 λ H λ H h 2 2 u (0) = 2 μ U 0 h (73) 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: U 0 = h ω u (0) = h ω a 2 2 h 2 (1 + 2 β - β 2 1 + 2 ϕ) 2 π h λ = π a 2 ω λ (1 + 2 β - β 2 1 + 2 ϕ * h) (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.
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).
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.
velocity
medicine and health sciences
body fluids
classical mechanics
fluid mechanics
mechanical stress
cellular structures and organelles
mucus
flow field
fluid dynamics
continuum mechanics
fluid flow
physics
shear stresses
cell biology
anatomy
cilia
momentum
physiology
biology and life sciences
physical sciences
motion
A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part II, modeling