CILIUM HEIGHT DIFFERENCE BETWEEN STROKES IS MORE EFFECTIVE IN DRIVING FLUID TRANSPORT IN MUCOCILIARY CLEARANCE: A NUMERICAL STUDY

Mucociliary clearance is the first line of defense in our airway. The purpose of this study is to identify and study key factors in the cilia motion that influence the transport ability of the mucociliary system. Using a rodpropel-fluid model, we examine the effects of cilia density, beating frequency, metachronal wavelength, and the extending height of the beating cilia. We first verify that asymmetry in the cilia motion is key to developing transport in the mucus flow. Next, two types of asymmetries between the effective and recovery strokes of the cilia motion are considered, the cilium beating velocity difference and the cilium height difference. We show that the cilium height difference is more efficient in driving the transport, and the more bend the cilium during the recovery stroke is, the more effective the transport would be. It is found that the transport capacity of the mucociliary system increases with cilia density and cilia beating frequency, but saturates above by a threshold value in both density and frequency. The metachronal wave that results from the phase lag among cilia does not contribute much to the mucus transport, which is consistent with the experimental observation of Sleigh (1989). We also test the effect of mucus viscosity, whose value is found to be inversely proportional to the transport ability. While multiple parts have to interplay and coordinate to allow for most effective mucociliary clearance, our findings from a simple model move us closer to understanding the effects of the cilia motion on the efficiency of this clearance system.


1.
Introduction.With breathing, the human respiratory tract is constantly in contact with potentially infectious microorganisms and noxious substances in air.Lung epithelium is our first line of defense against these infectious or damaging particles.On the other hand, aerosol inhalation is an effective means of drug delivery for respiratory problems.In order to better protect ourselves from these external harmful agents, or to design a more effective delivery device for aerosol drugs, we will need to understand how these particles are transported.Where particles land in the airway through the airflow is one aspect of the picture; how these particles are redistributed via the mucociliary clearance (MCC) is the other important 1108 LING XU AND YI JIANG aspect of the system.The epithelium lining in the airway consists of several cell types, including the goblet cells that produce mucus and the ciliated cells that have cilia projected from the cell surface into the airway.Mucus traps foreign particles.The cyclic beating of the cilia drives a directional movement of the mucus layer, transporting these particles, and eventually moves them out of the airway.
The physiology of MCC has been well documented, including the composition and rheological properties of mucus, the structure, density and coordinated movement of cilia [1] [11] [38] [37] [39] [41] [42].Direct measurements of these factors are obtained from various experimental systems.In addition, because the compositions of cell and mucus differ at varying parts of the lung [32], one should expect the measurements being highly variable along the ciliated airway.Table 1 lists some of these factors, as well as their normal ranges in human and rabbit.These references in Table 1 contain just a small portion of all experimental studies of MCC.It is observed that the effective height change of the cilium during its cyclic motion is important to the mucus propulsion.Impaired ciliary motion typifies a number of human diseases [10] [24] [31] [46].Cystic fibrosis [28] and primary ciliary dyskinesia [34] are two primary examples of ciliary function failure, either in the cilia motion or in the cilia-mucus interaction.However, the detailed mechanism underlying cilia beating and mucus transport is not well understood, for example, which aspect of a beating cilium or multiple cilia contributes most to the fluid transport?We aim to address this question with a mathematical model.
Many mathematical models treat MCC as a problem of fluid-body interaction, which have shed insights into the aggregated effects of the mucociliary interactions.Barton & Raynor [3] modeled the cilium as a rigid rod, which shortened automatically during the recovery stroke.They used a resistance coefficient to approximate the impacts between the cilium and mucus flow.Blake ([4], [5]) developed the slender body theory, in which the individual cilium was modeled as an array of force singularities along the cilium centerline.Stemming from the slender body theory, Cortez [12] introduced the method of regularized stokeslets in which the singular force was smoothed.The method is widely applied in fluid-flexible body interactions in stokes flow, e.g.helical swimming [13] and bacterial flagella bundling [19].
More recently, direct numerical simulations of MCC have taken into account more detailed cilia structures and mucus hydrodynamics.One focus of these works has been the formation of the metachronal wave in MCC.Yang, Dillon & Fauci [47] considered the axonemal structure of the cilium and used elastic springs to mimic the cilium microtubule filaments as well as the dyneins and nexin links.The mucus flow is governed by the incompressible Navier-Stokes equations in the velocity-pressure form.Their compuation demonstrates the formation of synchrony and metachrony of cilia due to hydrodynamic couplings, and a net fluid transport is observed.A similar model was adopted by Dillon et al. [14] to simulate sperm flagella in a viscoelastic fluid flow.Mitran [33] used a 3D finite volume method to study the formation of the metachronal wave in rows of pulmonary cilia.He considered two fluid layers: mucus layer as a viscoelastic gel and the periciliary layer as a viscous fluid.The internal microtubule structure of an individual cilium was modeled using large-deflection, curved, finite element beams.The model has been applied to simulate motions up to 256 cilia.Elgeti & Gompper [16] simulated 2D cilia arrays using chains of beads in a 3D fluid flow.They found that the metachronal wave increases fluid propulsion significantly.In these studies, the metachronal waves emerge as a result of mucus-cilia interactions.The metachronal wave itself cannot be varied independently.The question "what is the key factor driving the fluid transport" has not got a clear answer.
In this work, we propose a rod-propel-fluid model to simulate the MCC process.This model enables us to isolate the parameters of the cilia motion and examine their effects individually.We focus our attention on the effects of asymmetry in the cilia motion on the transport ability of the mucus flow.At the cilium level, we study the cilium shape change, as well as the velocity difference, between the effective and recovery strokes.At the tissue level, we examine the cilia density and the phase lag among cilia that is associated with the metachronal wavelength.Our numerical results show that, contrary to some previous results, details of cilia shape change and the metachronal wave do not matter, and that the maximum cilium height difference between strokes is the key factor deciding the transport capability.In the following sessions, we will describe our model and demonstrate our results leading to this conclusion.

Methods.
2.1.Model setup.We start building our model with the following assumptions: 1. Our simulation is 2D, assuming that the motion of the cilia is in a plane.
In cultured human airway epithelium, the ciliary beat is found to be clearly planar [10] [9] [41].2. The fluid surrounding the cilium is homogeneous, incompressible and Newtonian.The fluid motion governed by the incompressible Navier-Stokes equations.This treatment effectively immerses the cilium in the periciliary fluid, simulating the setting of many in vitro cilia flow experiments of the liquid-air interface [41].3. The cilium is modeled as a rigid rod.The extending height of the cilium varies during effective and recovery strokes, which is modeled via changes in the rod length.This is motivated by [3]. 4. The cilium motion follows Equation (1).The numerical details are provided together with descriptions of Figure 2.
Those four assumptions provide a basic computational frame in our simulations.We realize that these simplifying assumptions might limit the applicability of the results.Our intension is to start with a simplest possible model, gaining insight into the basic principles, and then add the complexity gradually only when necessary, including increasing the number of cilia, modifying the cilium height, prescribing a different phase shift among cilia, and modulating the beating frequency.
Treating the cilium as a rigid rod is a significant simplification in the model.Sanderson & Sleight [38] depicted a typical motion of the cilium in one beating cycle (Figure 1a), based on their observations of the rabbit tracheal cilia.During the effective stroke, the height of the cilium barely changes; during the recovery stroke, the bending of the cilium shortens the height.We want to investigate this asymmetry between strokes.In our model, we mimic the bending of the cilium as the shortening of the rod length, similar to Barton & Raynor [3].The length of the rod is a function of its orientation θ.The lower and upper bounds of the orientation are approximately 0.8 to 2.4 radians (40 and 140 degrees in [38]).
We measure the effective length of the cilium, L c , as the straight distance from its tip to the base, shown as arrow-lines in Figure 1a, and fit all measured values using a third order polynomial, yielding the expression below L c = { 0.6 effective stroke, −0.20993θ 3 + 1.21485θ 2 − 2.14088θ + 1.64268 recovery stroke. ( (a) Numbers from 1 to 9 show consecutive motions of the cilium during the recovery stroke and numbers from 9 to 12 show that during the effective stroke.Arrow-lines measure the distance between the cilium tip and its base.(b) The rod length L c is plotted as a function of its orientation θ in our model, Equation (1).The rod has the same length as that of the arrow-line in (a).The solid line in (b) is for the rod length during effective stroke and the dash line is for the recovery stroke.
Figure 1b plots L c vs. θ in (1).All values in the figure are dimensionless.The angle θ is in radian, and L c = 0.6 corresponds to 6µm with a typical cilium length 10µm.The maximum length of the rod is L c = L max = 0.6 and the minimum is L min = 0.45.The motion of rods is prescribed with a given beating frequency and phase shift.Thus, the numerical model is in fact a moving boundary problem, where the beating rods propel the fluid flow.
2.2.Governing equation and numerical method.The fluid is assumed to be incompressible with a constant density.The Navier-Stokes equations in the formulation of fluid velocity u and the pressure P are used Here f is the force exerted in the fluid flow, and µ is the hydrodynamic fluid viscosity.We define the characteristic length L as the trace of the cilium tip in one beating cycle which is about L = 32µm, since a typical length of the cilium is about 10µm and the cilium rotates 3.2 radians in one cycle (from 0.8 to 2.4 radians, forth and back).The characteristic time is chosen to be T = 1/15sec as 15 Hz is a typical beating frequency (Table 1).The Reynolds number is Re = ρL 2 µT = 0.0172 with ρ = 1000kg/m 3 and µ = 8.94 × 10 −4 P a • sec.Here we use Reynolds number Re = 0.01 in computations.The effect of the Reynolds number is discussed in the Supporting Materials.
A schematic plot of the computational domain is shown in Figure 2. The domain is a 2D box, Domain = {(x, y), x min ≤ x ≤ x max , y min ≤ y ≤ y max }. ( The bottom of the box is the epithelial cell surface where the mucus flow is assumed no-slip; the top is a free surface with zero traction; the left and right sides are periodic.The rigid rods projecting from the bottom of the box are our idealized cilia.The rod has a Lagrangian representation and is parameterized by its arc length s.The velocity of the rod is u p (s, θ) = (sω cos θ, sω sin θ), where θ is the rod orientation and ω is the rod angular velocity.For multiple rods, a uniform phase shift φ is assigned among rods to describe the metachronal wave.The motion of all rods are synchronized if φ = 0.In our model, φ has an upper bound set by the number of the rods N rod in order to avoid rods entanglement.We consider both symmetric and non-symmetric rod motions.Two cases of non-symmetric motions are studied.One is that the rod beating frequency varies between the effective and recovery strokes, which can be considered as varying the cilium beating velocity while the cilium length is fixed.The other is to fix the rod beating frequency and vary the rod length between these two strokes.These two cases are most obvious asymmetries observed in experiments.Other asymmetries in 2D, for example, the cilium beating velocity is a function of its orientation, could be thought as a variation/combination of these two cases above and will not considered here.A finite difference grid is established.The grid point is denoted by (x i , y j ).Initially, the fluid flow is quiescent with a velocity of u(z, 0) = 0, where z = (x, y).Each rod consists of m discrete points where m = 200 for all simulations reported in this paper.A singular force value F(s, t) is present at each point z p (s, t) along the rod, z p (s, t) = (x p (s, t), y p (s, t)).The total force density of one rod is We apply the projection method [8] in combination with the Marker-And-Cell (MAC) staggered grid technique [23] to solve the Navier-Stokes equations ( 2).This combination helps to avoid the numerical instability that usually arises when evaluating velocity and pressure at the same grid point.A formally second order immersed boundary method [25] is implemented to compute the impact of the moving rod to the fluid flow.The immersed boundary method was first introduced by Peskin [35] to simulate aspects of cardiac dynamics, which has since become a major technique in biological fluid dynamics.We place particles with no mass in the fluid, and they are passively transported at the fluid velocity u.At any time, these particles in the flow serve to visualize streaklines as may be observed in laboratory experiments.The position of the fluid particle z is computed using dz/dt = u(z, t).In the simulation, the fluid particle is injected at (0, 0.75).The second order explicit Adam-Bashforth scheme [27] is applied to solve fluid particle positions based on velocities at the current and previous time steps.The verification of our numerical methods, including the convergence study, is in the Supporting Materials.Table 2 gives parameters for the problem setup and Table 3 lists parameter sets that are used to describe the rod motion, including ranges of the rod number, rod density and beating frequency.The angular velocity ω = 3.2f .When we simulate more than one rod, all rods are equally spaced in x = [−1, 1].

3.1.
No net fluid transport with symmetric rod motion.We start with a symmetric motion with one rod, and then systematically modify one parameter at a time.By symmetric motion, we mean that the cilium shape is the same between effective and recovery strokes, which corresponds to a constant rod length in the model, and the beating frequency is the same between stokes.
Figure 3 groups all results of symmetric rod motion.The computational parameters are listed in the caption.Figure 3a shows streamlines induced by five rods during t = [0.002,0.066].The angular velocity of each rod is 50rad/sec, corresponding to a beating frequency of 15.625Hz.The rod motion is periodic, leading to a periodic flow motion.
As the rod beating starts, it drives the flow into motion, immediately forming a recirculation region near the rod tips.Figure 3a shows that the instantaneous streamlines are dense near the rod tip, indicating a large flow speed.At the center of the recirculation region, the fluid velocity is zero.As the rods swing from the right to the left (t = 0.002 to 0.026), the rods drive the recirculation region in the same direction.The center of the rotating flow initially appears close to the rod, gradually shifts up, then sheds off from the rod tip, and eventually vanishes at the free fluid surface (t = 0.034).Around t = 0.034, the rods reverse the beating direction.The rods induce another small recirculation region near the rod tip, followed by a similar process as before, with streamlines now having an opposite sign.This new recirculation region grows in size as it moves to the right in accordance with the rod motion, shifts upward and vanishes at the top surface.As the rod motion repeats its cycle, the fluid motion repeats as well, with the same period as the beating rods.In figure 3a, during t = 0 to 0.026, the rod is moving from right to left, and during t = 0.034 to 0.066 the rod is moving from left to right.Note that the fluid close to the rod tip moves in the same direction of the rod owing to the no slip condition at the rod-fluid interface.When the rod number is decreased to one, similar periodic fluid motions remain.The only difference is that the recirculation region is smaller.
The plots of streamlines change only slightly when we add a phase shift of φ = 0.2 among these five rods (not shown).Such an out-of-phase rod motion forms a metachronal wave of a wavelength 1.However, since the rod motion is symmetric, the reverse strokes bring the fluid back by exactly the same amount, resulting in a net zero displacement of the fluid particle in one complete period.Figure 3b plots the normalized horizontal u x /u max and vertical u y /u max velocities along the vertical line of x = 0 in the computational domain.u max is the maximum flow speed, and therefore, the normalized velocities range from −1 to 1.The velocity profiles are statistically equally distributed on two sides of the zero line, indicating that the net velocity is almost zero.A more precise measure of the material transport is on the displacement of material points, as indicated in Figure 3(cd).
While the streamlines cannot reveal the exact period of the flow motion and its transport ability, the streaklines of fluid particles can.Particles are released into the fluid continuously, and t p refers to the time when the particle is released.Figure 3d are snapshots of streaklines and they show that these particles follow periodic paths as long as the rod motion is symmetric.The period in the motion of these particles is the same as the beating rod.At each time snapshot, we display all fluid particles that have been released in the fluid flow previously.The left column of Figure 3d shows streaklines from t = 0 to t = 0.084 which contains a whole period of the rod beating cycle, 0.067.Streaklines form a closed loop, suggesting that the fluid particle comes back to its original location.The observed period of the fluid particle motion is about 0.07, close to the period of the rod motion.The right column plots streaklines at later times, from t = 0.1 to 1.3, showing that all fluid particles follow periodic motions all the time, with a zero net displacement.It is worth mentioning that the periodicity of the computational domain (periodic boundary condition in x) means the simulation domain repeats indefinitely, eliminating the finite size effect of the simulation domain.
We denote by D the total displacement of the particle, Here x D , y D are the displacement of the fluid particle in the x, y direction, respectively.Figure 3c shows the total displacement D of all the fluid particles present at t = 2 as a function of the t p , using one rod, five synchronized rods (φ = 0), and five non-synchronized rods (φ = 0.2), from left to right, respectively.Results at other times show qualitatively similar behaviors.It is seen that increasing the rod number from N rod = 1 to 5 helps increase the particle displacements by about four-fold.With a phase shift of φ =0.2 among rods, the maximum D is larger than that of φ = 0. Proportionality between the maximum value of D and φ will be investigated in the next section.

3.2.
Net fluid transport with asymmetric rod motion.Two non-symmetric rod motions are considered in this section.One is that the cilium beating velocity changes between the effective and recovery stokes while fixing the rod length.The other one is to fix the beating frequency, but take into account the shape changes in the cilium between strokes by shortening the rod length during the recovery stroke.The net displacement of particles is in the direction of the effective stroke, or toward the positive x axis in our simulations.
3.2.1.Cilium beating velocity difference affects the transport ability.The cilium beating velocity is often slower in the recovery stroke [41], which is a means of asymmetry in the cilia motion.We consider a constant cilium velocity in either the effective or recovery stroke, but apply a velocity difference between strokes.Considering a fixed cilium length of L = 0.6 and the cilium tip velocity is always 3.2f L, we report the results in terms of the beating frequency f for convenience.Three cases of beating frequency difference are computed.The frequency of the effective stroke is doubled in each case from 7.8125 to 15.625, and to 31.25Hz.The frequency of the recovery stroke is always 7.8125Hz.These values of frequency are within the physiologically allowable range (Table 1).A more detailed range can also be obtained from Sears et al. , where velocity profiles of recovery and effective strokes are measured.We calculate the frequency differences from these measurements as follows: use the maximum effective tip velocity as 375µm/sec and the minimum recovery tip velocity as 10µm/sec, thus the beating frequency would have a higher bound of about 25Hz in the effective stroke and a lower bound of about 0.5Hz in the recovery stroke.Figure 4a shows streaklines of fluid particles at t = 2 in these three cases of cilium velocity asymmetry, and Figure 4b plots the corresponding displacements of all fluid particles.A larger transport distance is seen at a higher beating velocity difference.However, the key information of this figure is that, an exponential increase in velocity difference does not guarantee an equal increase in the transport.One might also notice that the particles have slightly displacements in the y direction as well.The displacement in y is far less than that in the x direction, thus one could expect that D is approximately the same as x D .Nevertheless, one will see soon in the next figure that the transport distance induced by the cilium beating velocity asymmetry is lower than those induced by cilium height differences between strokes.

3.2.2.
Cilium height difference between strokes is key to transport capability.The cilium bends during the recovery stroke (Figure 1a), which is mimicked by shortening the rod length in the simulation.Figure 5a shows streamlines with five beating rods with the rod length changing as in (1) or Figure 1b.In Figure 5a, from t = 0.26 and 0.66, two rotating regions are present and the small recirculation region at the left is well defined; at t = 0.34, there is only one dominant recirculation region.
Figure 5b shows the resulting streaklines at a sequence of times, which reveal marked differences from the symmetric case (Figure 3d).In Figure 5b, a clear preferred direction of these particles is seen, which is along the positive x axis.We also observe a slight displacement in the y direction, suggesting that particles move toward the epithelial surface.This downward movement is due to the different boundary conditions at above (free surface) and below (no slip).
As the asymmetry in the cilium height between effective and recovery strokes can have such a drastic effect on the particle transport, we carefully examine the rod length function that depicts such an asymmetry.Figure 5c compares different rod length functions.For a better illustration, we compare in Figure 6 the maximum particle displacements due to the cilium beating velocity difference and the cilium height difference between strokes.The comparison is meaningful since parameters in these computations are around the physiological values.The particle displacement resulting from the cilium height asymmetry is a few folds larger than that of the cilium beating velocity asymmetry, suggesting that the former is a more effective mechanism for driving the mucus transport.

3.2.3.
Transport capability increases and saturates as rod density or beating frequency increases.In order to examine the transport capability of the MCC as a function of the tissue level parameters, we systematically vary three parameters in our model, the rod number N rod , the rod beating frequency f , and the phase shift φ among rods, with their values listed in Table 3.The rod length follows (1) for all computations in the following sections.Figure 7 shows results of max(D) at t = 2 as functions of the rod number N rod at four beating frequencies, as indicated in the plot.The phase shift is fixed to be φ = 0 in order to avoid rod tangling.Figure 7(ab) visualize the same data in two different ways, providing different aspects of the results.Figure 7a plots max(D) vs. N rod , and Figure 7b plots max(D) vs. the beating frequency f .In Figure 7a, max(D) increases as the rod number increases at all beating frequencies and then saturates at N rod = 9.Meanwhile, when the rod number is fixed, the faster they beat the more they move the fluid.Alternatively, in Figure 7b, it is interesting to see the saturation in the maximum displacement max(D), when the value of frequency f is greater than 25Hz and the rod number is less than 5.

Phase shift between rods has no effect on transport capability.
To study the effect of the phase shift in the metachronal wave on the transport efficiency, we fix the rod number N rod = 5, and examine the joint results of the beating frequency and the phase shift.Figure 8 shows max(D) vs. f at various values of the phase shift,  4. Discussion.From arsenical cobalt that destroyed miners' lung in sixteenth century Europe, to automobile fumes that hurt us in our everyday life, and to anthrax spores as potential bio-weapons, environmental aerosol particles have long been recognized as a potential threat to our health.Meanwhile, medication in the aerosol form can be delivered painlessly from the airway using inhalers.Examples include the most commonly used pressurized metered-dose inhaler for asthma medication and the dry powder inhaler for delivering insulin to diabetic patients.Where particles come in contact with the epithelial cells and how they spread on the tissue surface are crucial to understanding both the initiation of damage and the efficiency of these inhalers.The mucociliary system on the lung epithelial cells is the last barrier before these particles reach the tissue surface.It is important for both protecting our airway health and treating ailments to understand how the cilia beating drives the fluid motion and distributes particles.
Our study attempts to isolate the effects of major influencing parameters of the beating cilia in order to provide an insight into basic bio-mechanical mechanisms of particle transport in mucus.Such a methodology of isolating each parameter is difficult, if not impossible, in the wet-lab experiments because of the inter-dependence among parameters.For example, as one attempts to control the cilia beating frequency via changing the calcium concentration, the coordinated motion of the cilia changes as well [40].Experimental images have captured the height difference of a beating cilium and the cilia synchronization [17] [38].However, only a few studies have tried to quantify their effects [30] [41] mainly because of the intrinsic messiness of cilia motions in the airway.The facts that the cilia are short, densely packed, and fast moving, all make it difficult to quantify the exact shape of a beating cilium and "convert these ciliary shape changes into forms that are useful for testing theoretical models of ciliary function" [41].The fluid in reality, however, is more complex.Mucus is a gel-like mesh-work of large molecules and behaves viscoelastic.With mucus sitting on top of the periciliary fluid, the multiphase flow is not as simple as the Stokes flow.A recent gel-on-brush proposal [6] suggests that mucin macromolecules form dense mesh-work near cilia tips with each cilium acting like a brush, which further complicates the fluid flow picture.
We use a rod-propel-fluid model to mimic the the beating cilia driving the mucus flow.Our model simplifications are motivated by in vitro human airway epithelial cultures, where the cilium motion is seen to be planar and the cilia are at the airliquid interface [41].With such simplifications, this model enables us to take full control of the cilia motion and capture the qualitative behaviors of the system to a large extend.Our model confirms that the asymmetry is the key to achieve the net transport in a fluid flow at a low Reynolds number, which is a well-received fact since Purcell [36].This test also serves as a verification that our model works properly.
Conventional belief has been that "Airway cilia depend on precise changes in shape to transport the mucus gel overlying mucosal surfaces" [41].Numerically, we have applied different functions to mimic the orbit of the cilium tip.It is interesting that our results suggest the key factor contributing to the magnitude of the net transport is the maximum cilium height difference between strokes, despite all the complexity involved in the MCC.In other words, the precise shape of the cilium does not matter; only the height difference, from the base to the tip of the cilium, between effective and recovery strokes, affects the transport efficiency.Similar conclusions are observed in [26] which studied shortened airway cilia due to smoking.This finding is surprising, but is easy to understand if we consider the oarlike motion [2]: because the cilium does not bend much during the effective stroke, it pushes the fluid motion more effectively than during the recovery stroke, while the bending of the cilium reduces its interaction with the fluid.On the other hand, Sears et al. [41] showed that the beating cilium of in vitro human airway epithelium does not change height by much between strokes; they suggested that the cilium velocity difference between strokes be the main mechanism for transporting fluid.Our simulations agree that the rate difference between strokes alone can generate a net fluid transport.We further showed that the cilium height asymmetry is more effective in fluid transport than cilium velocity asymmetry.
If we carry on the oar analogy, the effects of cilium beating frequency and cilia density are easier to understand.It is expected that the faster the oar rolls or more oars there are, the faster the boat will travel.Thus it comes with no surprise that our numerical study suggests that higher cilia beating frequency and denser cilia array both help to increase the mucus transport.However, it is found that there exists an upper bound in the cilia density and frequency above which the transport ability ceases to increase.This observation implies a tolerance in cilia density which helps explain the fact that the transport capacity does not differ at varying parts of the airway where cilia density varies at these places.It also suggests that minor damages of the ciliated cells will not significantly affect the transport capacity, unless too many ciliated cells are killed such that the density drops significantly.
Braiman & Priel [7] expected the propelling velocity of the mucus is linearly related to the cilia beating frequency.Numerically, we also observe the positive relation between the cilia beating frequency and the maximum fluid particle displacement.However, it is worth pointing out that we found the relation is exponential (Figure 8).This implies that in order to increase the transport efficiency substantially, one does not have to give the beating frequency an equally substantial rise.A scaling of 3/2 is observed in the maximum displacement of fluid particles as a function of the beating frequency.So far, we do not know of an appropriate way to explain such a scaling.However, we did not exclude that this scaling might be an artifact owing to the simplified assumptions in our numerical method.
The formation of metachronal waves via hydrodynamic interactions between the cilia and mucus has been the focus of some recent studies [33] [47].In particular, Elgeti & Gompper [16] showed that the metachronal wave could drastically increase the propulsion velocity and transport efficiency.However, in their models, it was hard to vary the metachronal wave alone while fixing other parameters, thus a quantitative examination of the metachronism is difficult.In our model, the metachronal wavelength is prescribed, and this enables us to study its effects on the mucus transport directly.Our numerical results show that, on the contrary, the magnitude of fluid particle displacements does not depend on the phase lag among cilia or the existence of metachronal wave in the cilia array.This finding agrees with Sleigh [43] who found that the propulsion of the mucus does not require strong ciliary coordination.In addition, the empirical observations that metachronal waves are much shorter and less well organized in human explants than e.g. in vitro frog [17] but still transport mucus fine [41], also suggest that the details of the metachronal waves play a less significant role in the transport.
In cystic fibrosis, the mucus viscosity increases due to the airway surface dehydration, leading to significant malfunction of mucociliary clearance (Donaldson et al. 2006 [15]).The mucus viscosity is related to the Reynolds number in our problem.If the typical length and velocity of the problem are fixed, a larger value of mucus viscosity implies a smaller value of Reynolds number.In computations, as we reduce the Reynolds number from 0.01 to 0.005 (results are shown in Supporting Materials), the magnitude of the net transport D is brought down by approximately one fold.
In summary, our simple rod-propel-fluid model reveals some new findings regarding the essential effects of cilia motion on the mucocilairy clearance.These new findings moves us closer to understanding the complex dynamical system, which may indicate ways to improve biological or physiological experiments, and potentially help the design of artificial fluid-particle transport systems.U is the mean traveling speed, computed by averaging the horizontal velocity of the swimming sheet in the domain over one time period.Similar comparisons are used by Fauci and Peskin [18] in the study of aquatic animal locomotion.The asymptotic solution of Taylor [44] is and the solution of Tuck [45] is The Reynolds number is calculated as Re = ρω/µκ 2 .Parameters in the verification are chosen to be the same as that in Fauci and Peskin.For convenience, we list their values in the table below.

Table 4. Values of Parameters
Figure 9a shows the computation domain and location of the swimming sheet.Figure 9b shows the convergence of both normalized horizontal and vertical velocities along x = 0.1 in the domain.It is seen that velocities converge.Figure 10 compares U/V among our simulation results (labeled as Computed ), and Taylor and Tuck's asymptotic analytic solutions.The problem setup of the mucociliary model is described in great details in the paper and will not be reiterated here.The convergence study of the velocity field is performed at t = 1 with five beating cilia (rods).Figure 11a plots the streamlines at t = 1.The length of the rod is always L c = 0.6, the beating frequency is 15.625Hz, and the phase shift among cilia is φ = 0.2. Figure 11b is the convergence study of the normalized horizontal and vertical velocities along x = 0 in the domain.Figure 12 shows the total displacements D of all the fluid particles present at t = 2 as a function of t p , It is seen that increasing the rod number from N rod = 1 to 5 helps increase the particle displacements almost by four-fold.Similar results are observed at Re = 0.01 (Figure 3 in the paper).
Figure 13a plots two functions of rod length L c against the rod orientation, and Figure 13b shows the displacements of all fluid particles present at t = 2 using these two functions.It is seen that the transport ability is still identical as long as the difference in the maximum cilium height difference stays the same during the effective and recovery stokes.E-mail address: xulingsue@gmail.comE-mail address: yjiang12@gsu.edu

Figure 1 .
Figure 1.(a) A schematic plot of the cilium shape in one beating cycle (modified with permission from Sanderson & Sleigh, 1981).Numbers from 1 to 9 show consecutive motions of the cilium during the recovery stroke and numbers from 9 to 12 show that during the effective stroke.Arrow-lines measure the distance between the cilium tip and its base.(b) The rod length L c is plotted as a function of its orientation θ in our model, Equation (1).The rod has the same length as that of the arrow-line in (a).The solid line in (b) is for the rod length during effective stroke and the dash line is for the recovery stroke.

Figure 2 .
Figure 2. A schematic plot of the 2D domain.Two rods are shown, and θ 1 and θ 2 are their orientations, respectively.The phase shift is φ = θ 2 − θ 1 .The angular velocity is ω = dθ/dt.The rod length is L c .The rod is parameterized by its arc length s with s = 0 at the rod base and s = L c at the rod tip.

2 Figure 3 .
Figure 3. (a) Streamlines induced by five rods with time indicated in the plots.(b) Normalized horizontal and vertical velocities (u x , u y )/u max during t = [0 : 0.002 : 0.1] along the vertical line of x = 0 in the computational domain.(c) Plots of displacement D of all fluid particles present at t = 2 as functions of t p , the time when the particle is released into fluid.(d) Snapshots of streaklines induced by five rods at varying times.The rod number for Figure 3(abd) are five, the rod length is L c = 0.6, and the rod beating frequency is f = 15.625Hz between stokes.The phase shift among rods is φ = 0 in (a), and φ = 0.2 in (b) and (d).

Figure 4 .
Figure 4. (a) Streaklines of fluid particles induced by five rods at t = 2.The frequency of the effective stroke is f = 31.25 (top), 15.625 (middle), 7.8125 (bottom).The frequency of the recovery stroke is f = 7.8125 for all.The phase shift among rods is φ = 0.The rod length is fixed to be L c = 0.6.(b) Displacement D of all fluid particles present at t = 2 as functions of t p , the time when the particle is released into fluid, for the three cases in (a), respectively.

Figure 5 (Figure 5 .
Figure 5. (a) Streamlines induced by five rods at a sequence of times, as indicated in the plots.The rod length shortens, following (1) or Figure 1b.(b) Streaklines of fluid particles at a sequence of times.(c1) Seven cases of the rod length L c are plotted as a function of its orientation θ.The solid line is the length for the effective stroke, and the dash lines are for the recovery stroke.(c2) Compare the fluid particle displacements D of curves 1 to 5 in (c1).(c3) Compare the fluid particle displacements D of curves 3, 6 and 7 in (c1).The displacement D is plotted as a function of t p , the time when the particle is released into fluid.For all plots of Figure 5, the number of rods N rod = 5, the beating frequency is the same between strokes, f = 15.625Hz, and the phase shift among rods is φ = 0.2.

Figure 5 (
c2) compares the particle displacement D of curves 1 to 5, elucidating the effects of the cilium height difference, L max − L min .In the figure, five different values of L min are considered, 0.4, 0.425, 0.45, 0.475, 0.5, and L max is 0.6 for all.These five curves are expressed in terms of third-order polynomials.Recall that four conditions are required to determine their function expressions, and they are L c = 0.6 at θ = 0.8 and 2.4, L c = L min and dL c /dθ = 0 at θ = 1.36.The arrow line in the figure indicates results with a decreasing L min .The observed trend shows that a decreasing value of L min contributes to an increasing value of D.

Figure 5 (Figure 6 .
Figure 6.Comparison of the maximum displacement max(D) due to cilium beating velocity and height differences between the effective and recovery strokes, respectively.In the figure, three data points (a1, a2, a3) at L max − L min = 0 correspond to those in Figure 4b; the data points (1, 2, 3, • • • , 7) correspond to the 7 cases of cilium height difference in Figure 5(c2,c3).

Figure 7 .
Figure 7.The maximum displacement max(D) of fluid particles present at t = 2 is shown as functions of (a) the rod number N rod in the semi-log scale and (b) the beating frequency f in the loglog scale at four beating frequencies, f = 7.8125, 15.625, 25 and 31.25Hz.The phase shift among rods is φ = 0.The rod length changes during recovery stroke, following (1) or Figure 1b.

Figure 8 .
Figure 8.The maximum displacement of fluid particles present at t = 2, max(D) is shown as a function of the beating frequency f at various phase shifts, φ = 0, 0.05, 0.1, 0.15, 0.2, 0.25 and 0.3.A dash-dot line of slope 3/2 is for visual guide.The number of rods is always N rod = 5.Four beating frequencies are considered, f = 7.8125, 15.625, 25 and 31.25Hz.The rod length changes during recovery stroke, following (1) or Figure 1b.

Figure 9 .Figure 10 .
Figure 9. (a) Computational domain and location of the sheet at t = 0. (b) Convergence study of the relative horizontal and vertical velocities, u x /u max and u y /u max , respectively, along the vertical line of x = 0.1.The time is at t = 20 at which the result is already steady state.The mesh size h is indicated in the figure legend.u max is the maximum velocity magnitude over the whole computational domain.

Figure 11 . 5 . 2 .
Figure 11.(a) Streamlines at t = 1.(b) Convergence study of the relative velocity u x /u max and u y /u max in the x direction and y direction respectively, along a vertical line of x = 0.The mesh size in (a) is h = 0.1.

Figure 12 .Figure 13 .
Figure 12.Plots of displacement D of all fluid particles present at t = 2 as functions of t p .The results are for symmetric motions of (a) one rod, (b) five synchronized rods.The beating frequency for all is f = 15.625Hz.The rod length is fixed to be L c = 0.6.The Reynolds number is Re = 0.005.

Table 1 :
Experimental measurements of the cilium length, cilia density, beating frequency and tracheal mucus velocity in rabbit or human.

Table 2 :
Parameters used in computations: starting/ending times [0, T end ], the time step ∆t, the size of the computational domain [x min , x max ] × [y min , y max ], and the number of grid points N x × N y in x and y directions.

Table 3 .
Parameters used to control the rod motion: the number of rods N rod , the phase shift φ between rods, and the beating frequency of the rod f .