Three-Dimensional Numerical Analysis of Periciliary Liquid Layer: Ciliary Abnormalities in Respiratory Diseases

: Human pulmonary epithelial cells are protected by two layers of ﬂuid—the outer watery periciliary liquid layer (PCL) and the uppermost non-Newtonian mucus layer (ML). Aerosols and inhaled toxic particles are trapped by the ML which must then be removed swiftly to avoid adverse health implications. Epithelial cells are covered with cilia that beat rapidly within the PCL. Such ciliary motion drives the mucus transport. Although cilia can penetrate slightly inside the mucus to assist mucus movement, the motion of the underlying PCL layer within the airway surface liquid (ASL) is signiﬁcant in mucus and pathogens transport. As such, a detailed parametric study of the inﬂuence of different abnormal cilia characteristics, such as low beating frequency, short length, abnormal beating pattern, reduced ciliary density, and epithelium patchiness due to missing cilia on the PCL transport, is carried out numerically. Such abnormalities are found in various chronic respiratory diseases. In addition, the shear stress at the epithelium is assessed due to the importance of shear stress on the epithelial function. Using the immersed boundary (IB) method combined with the ﬁnite-difference projection method, we found that the PCL, under standard healthy conditions, has net forward motion but that different diseased conditions decrease the forward motion of the PCL, as is expected based on clinical understanding. the PCL locomotion for abnormal ciliary length. It was found cilia are, the more the PCL ﬂow rate is hindered. Compared to the standard length L = 6 µ m, the ciliary lengths of 5 µ m and 4 µ m led to 25% and 48% reduction in the ﬂuid ﬂow rate, respectively. However, increasing the length to 7 µ m could not alter the result signiﬁcantly—only 2% increase was observed. This may be related to the cilia tips getting closer to the PCL surface and resulting interactions, but the exact cause is left for future investigation.


Introduction
Airway surface liquid (ASL) covers the respiratory epithelia that is mainly comprised of two distinct layers: mucus layer (ML)-a non-homogenous, non-Newtonian, viscoelastic fluid [1]; and periciliary liquid layer (PCL)-treated as a watery lubricating (nearly Newtonian) fluid layer with much less viscosity [2]. Transport of the ASL is of vital importance to keep the respiratory system clear from inhaled particulate matters deposited on the ML. Such clearance mechanism is facilitated by periodic ciliary actions, with the ASL being swept forward during the effective stroke and backward during the recovery stroke [3].
Ineffective and poor mucus clearance is a major factor for lung health. There are generally two main reasons for muco-ciliary dysfunction in the respiratory system: airway mucosal diseases (such as cystic fibrosis (CF), chronic obstructive pulmonary disease (COPD), asthma, etc.) and ciliary abnormalities. Mucosal diseases occur in patients as a result of the conversion from healthy to pathologic mucus, where the mucus hydration and its biochemical constituents change by multiple mechanisms, leading to higher viscosity and elasticity of the mucus (see Fahy and Dickey [4], and the references therein). Defects in ciliary structure and function can cause abnormal ciliary beating. Primary ciliary dyskinesia (PCD) is known as the most significant disorder causing ciliary abnormalities [5].
There is fairly a limited number of numerical studies in the literature with the focus of modelling of diseased lungs. Guo and Kanso [6] emulated diseased environments in their 2D simulations by varying the viscoelastic properties and the thickness of the ML. Their findings revealed that the effects of the mucus elasticity depends very much on the ML thickness. Specifically, for the diseased states (thinner PCL), the higher elasticity led to a decrease in the mucus transport, whereas an opposite trend was observed for the healthy states. Chatelin et al. [7] attempted to simulate the high viscous effect of pathologic mucus in CF three-dimensionally by taking into account the highly shear-thinning nature of mucus and thus the non-Newtonian effects on the muco-ciliary clearance. Based on the mucus samples collected from patients with and without CF disorder, they concluded that shear-thinning effects can be of the same order as the viscoelastic effects for the mucus transport. Jayathilake et al. [8] studied ASL flow physics two-dimensionally in the presence of cilia beating under various ciliary abnormalities, such as ciliary beat pattern, length of cilia, immotile cilia, beating amplitude, and uncoordinated beating of cilia. According to the results obtained, each of the above ciliary abnormalities could greatly reduce the process of mucus clearance, except for uncoordinated beating.
The ultimate goal of this paper is to obtain a clear understanding of how the epithelial cilia-driven fluid flows perform under various ciliary abnormalities in respiratory diseases. Such abnormalities, whether inherited (genetic disorder) or acquired, are directly linked to lung infection and chronic respiratory diseases in humans. Table 1 summarises the most common lung disorders along with their damaging effects on the ciliary function. More details on varioues types of ciliary abnormalities, such as PCD, CF, asthma, and COPD, can be found in [5,9,10]. Table 1. Ciliary abnormalities due to inherited and acquired disorders.

Inherited Disorders Ciliary Abnormalities
We therefore performed a parametric analysis to investigate the effects of ciliary abnormalities, like low beating frequency, short length, abnormal beating pattern, and reduced ciliary density, as well as epithelium patchiness due to missing cilia on the PCL transport as an extension to previous work reported in Reference [23]. To achieve this, the three-dimensional equations of fluid motion are solved using the IB method coupled with finite-difference projection method. Some theoretical studies have concluded that the PCL remains stationary during the process of clearance [24][25][26][27][28][29]. However, experimental findings [30] argued that if PCL was stationary, inhaled hydrophilic toxicants might remain in the lung for extended periods of time. Hence, axial PCL transport is needed in order to efficiently clear the lung from this class of toxic substances. Overall, by analysing the transport behaviour of PCL, the present 3D model not only could provide valuable insights into muco-ciliary clearance, the results of this parametric study could be useful for designing lab-on-a-chip microfluidic devices which are used for pumping and/or mixing purposes [31][32][33].
The paper is organized as follows: Section 2-details of the numerical method are given; Section 3-the results and relevant discussions are presented in terms of PLC flow rate, velocity changes above the tips of cilia, and alterations in shear stress at the epithelium; and Section 4-conclusions for the results are provided, and directions for future research are suggested.

Numerical Method
In this work, the time-dependent Newtonian PCL flow is solved by a projection method and the force exerted by cilia on the fluid is computed via the IB method. Figure 1a shows where u = (u, v, w) is the fluid velocity, p is the fluid pressure, and ρ and µ are constant density and dynamic viscosity of the fluid, respectively.  [34], which shows the changes in shape of each cilium X = (X 1 , X 2 , X 3 ) at intervals of T/12 through one beat cycle (0 < t < T). Cilia move in a planner motion during their rapid effective stroke (phase No. ∼ 1-7 ), but the recovering cilia (phase No. ∼8-12 ) are not performing a planner beat. When viewed from above, the ciliary tips follow a anticlockwise oval path. (d) The cilium velocity components. (e) Interpretation of ciliary activity on cultured rabbit tracheal epithelium from the model of Sanderson and Sleigh [3]. Arrows e and m indicate the directions of effective stroke and metachronal wave, respectively. (f) Adopted ciliary pack for all the simulation cases unless stated otherwise.
Since one portion of an infinite problem is simulated in the present work, periodic boundary conditions are applied to the side walls normal to the x and y axes. Herawati et al. [35] observed that the fluid velocity at the epithelium boundaries remain static due to the slow movements of the basal bodies (a protein structure found at the base of a cilium) at the cell boundaries of cultured mouse tracheal multiciliated cells. Hence, the bottom boundary at the epithelium is treated as no-slip wall in our simulations. Given that the boundary condition at the PCL-mucus interface is influenced by the surface tension force between the ML and PCL, as well as the complex fluid rheology of the mucus, it is reasonable to assume the top boundary as the free slip wall in the present mucus-free model. The mentioned boundary conditions can be presented in the mathematical form as: , u| y=−24µm = u| y=+24µm periodic boundaries in the xand y-directions. ( Note that in this work the geometry is described using terms such as top and bottom, but within the respiratory tracts, the orientation may be in any direction. Gravity is not considered in the present model. To study the flow-structure interaction in biological systems, the immersed boundary (IB) methods have been applied successfully. Employing this method [36][37][38][39] (see Figure 2b), each cilium can be represented by a set of N control points X l , l = 1, . . . , N as functions of Lagrangian coordinate S (0 < S < L) and time t (0 < t < T) embedding in the Eulerian fluid domain. Thus, the cilium length can be given by L = ∆S.N where ∆S is the Lagrangian grid spacing along the IB, as depicted in Figure 2a. Here, T is the period of cilium beat cycle. The control points can be treated as force generators to the fluid, which are computed by using the following direct forcing method at Cartesian grid point (i, j, k): where ∆t is the marching time step between the time levels n and n + 1, and − → U n+1 i,j,k is the distributed cilium velocity which can be obtained as Since the cilium control points − → X l do not coincide with the fluid grid points x i,j,k , the filter function is employed to transfer the information between those two grids [40], which is defined by where h is the Eulerian grid spacing (see Figure 2a). It should be noted that ∆S/h < 1 to avoid penetration of streamlines into the Lagrangian grid spacing. Three-dimensional kinematics of the cilia U(X l ) are prescribed using the ciliary beat profiles shown in Figure 1b,c, which is introduced as the 'standard beat' throughout this work. It is assumed that cilia beat with the frequency of CBF = 1/T. It can be clearly seen in Figure 1d that, through one beat cycle, the magnitude of the cilium velocity components in the x-direction U is significantly higher compared with the other two velocity components-V in the y-direction and W in the z-direction. This means that the driving force imparted by cilia is dominant in the x-direction. The packing among the cilia is also depicted in Figure 1f. The metachronal wave produced by actions of cilia travels approximately at an angle of 135 • to the effective stroke, as reported in the experiment on ciliary activity of cultured rabbit trachea epithelium by Sanderson and Sleigh [3] (see Figure 1e). Hence, considering one full wave over the computational domain using the same propagation angle, the initial arrangements of the cilia are chosen such that they are in similar phase along the wave fronts, with a fixed phase difference of ∆Φ = π/6 and the interciliary spacing being 3 √ 2 µm. The velocity of the cilia U(X l ) is prescribed accordingly and the position of them at every time step ∆t is calculated as X n+1 = X n + ∆t × U(X l ). Once the force exerted to the fluid F was computed explicitly, Equations (1) and (2) are then discretised and solved for the velocity field and pressure using the pressure-increment projection algorithm, analogous to that described in Reference [41]. The staggered marker-and-cell (MAC) Eulerian mesh [42] is used for the spatial discretization (as shown in Figure 2c), and the intermediate velocity field u * is calculated as follows: where λ = 2ρ/µ.∆t and Since the u * (∼ u + ∇φ) does not satisfy the divergence free condition in Equation (2), a pressure increment φ n+1 is computed as To calculate the Helmholtz equation for u * in (7) and the Poisson equation for φ n+1 in (9), we employ FISHPACK [43] fast solvers. Having solved for the u * and φ n+1 , the velocity and pressure fields are then updated as The standard projection method has been successfully applied in many muco-ciliary studies [6,8,23,44,45], and all the authors have reported satisfactory values for the fluid transport when compared with the available experimental data. Nonetheless, the accuracy of this method at low Reynolds numbers has been questioned because the projection on divergence-free fields has been observed to induce tangential spurious velocities [46]. Although this is an important limitation to consider, the actual degree of loss of accuracy that might result is unknown for muco-ciliary systems. What makes the projection method attractive is that, at each time step, one only needs to solve a sequence of decoupled elliptic equations for the velocity and the pressure, which makes it more computationally efficient. However, it is necessary to evaluate the accuracy of the algorithm which is tested for the Couette flow, as well as for Stokes flow, in a lid-driven cavity, for Re = 0.01, which is a creeping flow condition of analogous flow regime to muco-ciliary transport. Figure 3a provides results for the Couette flow at a range of times compared with the analytical solution of the Navier-Stokes equations, which can be expressed as [47] where L y is the width of the channel, ν is the fluid kinematic viscosity, and u(y, t) is the fluid velocity calculated at position y and time t. Comparable to the results obtained by the analytical solution, the present results deviated by no more than 0.12% which is acceptable within the current context. As for the second validation (Figure 3b), the streamline results of the present work are in good agreement with those obtained by an analytical method in the works of Shankar [48] and Pan and Acrivos [49], supporting the use of the numerical algorithm to simulate Stokes flows for our purposes. Detailed discussions on the validity and robustness of the implemented IB method are found in Reference [23,44,50] and hence not repeated here. Table 2 summarizes the computational steps used in the present work. For the numerical simulations, we use the parameter set displayed in Table 3. Table 2. The numerical procedure for each time step of the present method.
Step 5: Update the velocity and pressure fields by Equations (10) and (11).

Cilia-Driven Fluid Flow
To assess the numerical convergence of the discretization scheme, we evaluate the effects of time step, fluid grid, and number of control points on the time-averaged fluid velocity in the x-direction. This is computed over the total PLC volume using u = 1/(TV f ) u(x, y, z, t)dV f dt where V f is the fluid volume and T = 0.3 s. The results show that the solution changes little under spatial and temporal mesh refinement and the discretization scheme reaches convergence over the third cycle, as seen in Figure 4. Therefore, the fluid grid size    Figure 5a. In addition to this, since the cilia beat most effectively in the x-direction, it is expected to have regions with higher velocities in this direction. The largest values of the mean velocity magnitude u are observed on the z = 8 µm plane (Figure 5b) which can facilitate the mobility of the ML. Figure 5c,d depict the xand y-component of velocity field (i.e., u and v ) at two representative planes: y = 0 and x = 0, respectively. It is observed that the maximum contour value of u is approximately four times larger than that of v . To quantify these flow changes, the time-averaged velocity profiles as functions of z are constructed in the middle of the computational domain (intersection of (x, z) and (y, z) planes), as seen in Figure 5e. These findings clearly confirm that a considerable amount of flow in the direction parallel to effective stroke is generated. Referring to Figure 5c, one can also observe that the main direction of the fluid transport above the cilia tips is towards the positive x-direction which agrees well with the actual clearance mechanism; unidirectional transport of the ASL over time. Furthermore, our calculations for the time-averaged fluid velocity which is computed over the total PLC volume in y-direction v = 1/(TV f ) v(x, y, z, t)dV f dt led to a very small value: v = 0.145 µm/s vs. u = 3.68 µm/s. Henceforth, we shall present our findings for flow parallel to the direction of cilia effective stroke.  shown with vertical dotted line for −24 ≤ x ≤ 24 µm (right panel). The cilia are labelled 1-6 from left to right. We now seek to understand the interaction between the fluid and cilia at various time instants. Viewing from above, many small size vortices are formed above the cilia tips due to the packed distribution of cilia, as well as the complex nature of cilia-driven fluid field. Such turbulence-like streamlines are due to the flow separations in the vicinity of ciliary actions. However, the overall fluid regime is still laminar. Outside the ciliary domain (i.e., |(x, y)|≥ 16.5 µm), it is more likely to have bigger recirculation zones. The elevation view shows a row of six cilia which are at different beating phase, depending on the simulation time. Case (a) shows a snapshot of the results at the beginning of the beating cycle, i.e., t = 0 s . In case (b), cilium 1 is seen to stand up and perform its effective stroke. The action of this cilium is pushing the fluid forward in the positive x-direction; the corresponding velocity profile shows this positive fluid motion. The cilium continues to move forward (case (c)), but compared to the previous case, the corresponding fluid velocity profile shows smaller values. Case (d) depicts the recovery stroke for this cilium where an opposing vortex is set up in the fluid above the cilium tip, resulting in a reverse flow with a negative velocity profile. One can note that this relatively strong vortex is formed above the other cilia with similar beating phase in the other snapshots. Eventually, at t = 0.08 s, which corresponds to the case (e), the cilium 1 start reaching to its previous initial position to begin the next cycle. This biological fluid-structure interaction is very similar to what is observed in nature. For example, birds extend their wings to overcome the viscous resistance and then bend them closer to their body to avoid cancelling out their effective forward motion.
To provide information on the PCL transport, the flow rate is computed as The surface mean velocity where the ML sits on top is also obtained as The fluid-flow induced wall shear stress at the epithelia are also another important factor of interest that, to our knowledge, has not been numerically investigated, which is calculated as It has been experimentally observed that the shear stress can alter the epithelial paracellular permeability (i.e., diffusion of ions or molecules via the intercellular space) significantly [53]. The magnitude of fluid vorticity is the last metric of interest which is averaged over the PCL volume and computed as Figure 7 shows the transients of flow rate, PCL surface average velocity, wall shear stress, and volume-averaged vorticity magnitude obtained from Equations (13)- (16). With 72 cilia beating at 10 Hz, the average of PCL flow rate, surface mean velocity, wall shear stress, and vorticity magnitude are obtained as follows: Q = 1533 µm 3 /s , u sur f = 8.19 µm/s, τ w = 0.08 dynes/cm 2 , and |ω| = 0.95 s −1 , respectively. Our surface velocity results are qualitatively consistent with those reported in Reference [54], but different quantitatively, with the time-averaged PCL mean velocity being 3.18 µm/s. The reasons for this discrepancy are due to differences in the adopted ciliary density and prescribed beating pattern. In addition, they considered a synchronized cilia beating (no methachronal wave) and a fluid with variable viscosity, both of which can contribute to the further difference.

Parametric Investigations
To achieve the goal of this paper, we perform thorough parametric investigations to evaluate the effects of parameters including ciliary beat frequency CBF, length of cilia L, ciliary beat pattern (CBP), ciliary pack, ciliary density, and missing ciliary patch on the simulation results. It will also be explained how changes in these parameters could represent some of the abnormalities listed in Table 1.
The results of our parametric study are assessed in terms of three metrics: (i) Time-averaged flow rate as (ii) time-averaged surface mean velocity as and (iii) time-averaged wall shear stress as It should be noted that in all of the following results, the numerical setting is the same as the standard setting given in Table 3, except when stated differently. Figure 8 shows the variation of time-averaged flow rate Q generated by the ciliary action for six different parameters. Parameter (I): Ciliary beat frequency (CBF). In the healthy human lung, cilia beat rapidly at about 8-15 Hz [55]. However, experimental findings in PCD shows a considerable CBF variability from 1.4-24.9 Hz due to gene mutations [56]. Further, in acquired disorders, such as cigarette smoking, environmental pollutants, and asthma, CBF tends to decrease [9,10,19]. To study the effect of CBF on the flow rate, the beat frequency is varied in the range of 5-25 Hz. The flow rate varies linearly with CBF, as shown in Figure 8a. The effect of CBF on the results is significant-a decrease in CBF from 10 to 5 Hz, due to diseased conditions for example, reduces the flow rate Q by 50%. Note that 10 Hz is considered as the standard or healthy beat frequency throughout this work.
Parameter (II): ciliary length (L). Clinical findings reveal that cigarette smokers, people exposed to environmental pollutants and patients suffering from COPD have shorter cilia [15,16,21]. As such, it is important to investigate the PCL locomotion for abnormal ciliary length. It was found that the shorter the cilia are, the more the PCL flow rate is hindered. Compared to the standard length L = 6 µm, the ciliary lengths of 5 µm and 4 µm led to 25% and 48% reduction in the fluid flow rate, respectively. However, increasing the length to 7 µm could not alter the result significantly-only 2% increase was observed. This may be related to the cilia tips getting closer to the PCL surface and resulting interactions, but the exact cause is left for future investigation.
Parameter (III): CBP. Studies of CBP in individuals with asthma and COPD showed impaired CBPs in bronchial epithelial [20] and nasal cilia [22], respectively. Two distinctive beat patters have been reported for cilia from PCD disease subjects: stiff and wiper planar motions [12,13]. The earlier is constructed by assuming a stiff cilium beats within a range of 30 • ∼150 • with twelve equally spaced phase numbers over one cycle. The later motion is constructed in the same manner but the range of motion is assumed 50 • ∼130 • , where a straight line pointing from the cilium base to the cilium tip makes the angle with respect to the bottom wall. In addition to these ciliary abnormalities, two-dimensional beating kinematics (planar) is also incorporated into the model and compared against the three-dimensional ciliary motion (standard beat). This motion is constructed by setting the y-component of 3D beat profile equal to zero.
Comparing the standard CBP to planar CBP, the flow rate Q was almost identical. This agrees well with the results of Ding et al. [57], who concluded that both 2D and 3D ciliary kinematics show qualitatively similar transport due to the metachronal waves. In the stiff planar and wiper cases, the fluid flow rate drops significantly: by 70% in planar and by 94% in wiper. Parameter (IV): ciliary pack. The spatial relationships of cilia on the epithelial surface is not very clear because cilia are tightly packed on the epithelia which makes it difficult to visualise them by means of scanning electron microscopy (SEM). Sanderson and Sleigh [3] interpreted an area of ciliary activity on cultured rabbit tracheal epithelium by a model in which cilia had equilateral triangular packing (Pack 1; Figure 9a). Here, we suggest two other possible packing types (see Figure 9a): isosceles triangle (Pack 2)-which is initially employed in the present work for the cilia-driven fluid flow analyses and then for the first three parametric studies; and square pack (Pack 3)-which has been commonly used by researchers [7,45,54,57]. To find out which packing type could be more effective in the fluid transport, the three mentioned packings are compared with one another. Note that the ciliary density (i.e., cilium/area taken up for each cilium) is kept constant at 0.055 cilium/µm 2 for this comparison. That is to say, only 64, 68, and 72 cilia could be fitted within the ciliary domain using Pack 3, Pack 1, and Pack 2, respectively. The results show that Pack 1 provides the highest flow rate which is slightly larger than that of Pack 1, revealing that although pack 1 has a fewer number of driving ciliary forces, it is capable to push the fluid more efficiently. However, the square pack (Pack 3) decreases the forward motion of fluid up to nearly 40%. These findings could play an important role in the optimal design of cilia-on-chip systems. Parameter (V): missing ciliary patch. Examination of the airway ciliated cells of smokers and CF patients shows loss of cilia [10,14]. Exposure to airborne pollutants is another reason of loss of ciliated cells [17]. In order to study the regions of missing cilia, three patches of different sizes are considered using Pack 3, as depicted in Figure 9b. As seen in Figure 8e, Q is inversely proportional to the patch size.
Parameter (VI): ciliary density. Decrease of ciliary density on the nasal epithelium is another respiratory impairment in COPD subjects [58]. To explore the effect of ciliary density, 36-576 equally spaced cilia are distributed (using Pack3) over the ciliary domain (33 µm 2 ) in the xand y-directions which correspond to inter-ciliary spacing range of 1.43-6.6 µm within the same computational domain. Given that the reported value for the spacing between cilia is 0.35 µm [59], further decreasing the spacing between two adjacent cilia can cause spurious numerical effects since the neighbouring cilia may overlap, and as a result, two Lagrangian points with opposite velocities can occupy the same fluid grid cell. It is also important to emphasise that the wavelength formed by motile cilia is kept constant for all of the ciliary densities. Thus, the phase numbers of cilia are interpolated according to their positions within the ciliary domain. The linear dependency of Q on frequency (Figure 8a) and number of cilia (Figure 8e) can be expected as the driving force linearly increases with frequency and number of cilia in the creep flow regime and therefore the average PCL flow would lineally depend on these two parameters. Figure 8f shows the relation between the ciliary density and Q . This simulation results show that although there is a decrease for the case with 64 cilia, the flow rate initially increases with increasing the number of cilia because the total driving force increases. However, when the cilia are too close to each other the flow generated by one cilium would be partially cancelled out by the surrounding cilia and one can therefore see a plateau in Q as the number of cilia approaches 400. In all of the simulations presented, results obtained for the 2nd and 3rd cycles have been almost identical, regardless of packing type, patch size, and ciliary density. Hence, the time-averaged results presented in this section are obtained over the cycle 3 with the period of T = 0.1 s.
It is now interesting to investigate the variation of time-averaged surface mean velocity u sur f for all of the six above-mentioned parameters (see Figure 10). One can clearly see that the curves exhibit relatively similar trends as the time-averaged flow rates Q , particularly in parameters (I), (IV), (V), and (VI); Figure 10a,d-f. This similarity is to be expected since the flow rate is the product of area and velocity. The linear relationship between u sur f and CBF in Figure 10a is also to be expected due to the negligible inertia forces (very low Reynolds number) for this problem which, consequently, results in a linear momentum Equation (1). Additionally, the experimental findings also indicate that ciliary forces depend linearly on the CBF [60]. Figure 10b shows that the surface mean velocity depends highly on the ciliary length L. Decreasing L from 7 to 4 µm, for example, could reduce u sur f by 206%, indicating that smoke-induced shortening of cilia could decrease the muco-ciliary clearance significantly. Interestingly, it was found that the standard beat is able to slightly improve the surface mean velocity compared with the planar beat (see Figure 10c), although it had less fluid transport capability. On the other hand, u sur f experienced a reduction of over 80% when using the stiff planar beat, and the surface remained almost stationary for the wiper motion. As mentioned in above, u sur f exhibits relatively similar trend as Q when comparing different ciliary pack(parameter (IV)). That is to say, Pack 1 has the largest value, followed by Pack 2 and Pack 3, respectively (see Figure 10d). Hussong et al. [61] measured the in-plane fluid velocities above the epithelium surface of mucus-free mice tracheae using the micro Particle Image Velocimetry technique and observed that regions where no beating cilia are located had smaller velocities compared to the regions with ciliated cells. This indicates that the cilia-induced PCL velocity slows down when patchiness occurs within the ciliated cells, thus hindering the effective mucus transport (Figure 10e); see Table 1 for disorders associated with cilia loss. As shown in Figure 10f, for 400 cilia and above, the surface mean velocity seems to level off and is not very sensitive to further changes in cilia number.
As mentioned earlier, little is known regarding the effects of fluid-induced shear stress on the epithelial function. The relationship between parameters (I)-(VI) and shear stress τ w at the epithelium is presented in Figure 11. The results show that the time-averaged wall shear stress τ w increases almost linearly with CBF (Figure 11a), while the increase of τ w with L is non-linear (Figure 11b). It is clear that the velocity gradients at the epithelium intensifies with the increment of beating frequency and ciliary length. The shear stress for different CBPs is compared in Figure 11c. The shear stress of the planar, stiff planar, and wiper CBPs is 20%, 16%, and 38% smaller than that of the standard beat pattern, respectively. As can be expected, if the number of cilia is decreased, there would be smaller fluid velocity gradients from the ciliary actions, thus decreasing τ w at the epithelium. As such, since Pack1 (with 68 cilia) and Pack 3 (with 64 cilia) have smaller number of cilia, they have lower levels of shear stress than Pack2 (with 72 cilia) (see Figure 11d). Similarly, the shear stress becomes smaller with decreasing ciliary density and increasing the size of missing ciliary patch, as depicted in Figure 11e,f. It is worth mentioning that exploring the effects of shear force on the epithelial function can be of very complex task biologically. Nevertheless, the experimental work of Sidhayewhich [53] showed that the shear stress can regulate barrier function, alter paracellular permeability, and provide insight into the homeostatic mechanism in the respiratory epithelia.

Conclusions
By using a finite-difference fluid solver coupled to an IB method, we first examined cilia-driven fluid flows, and then thorough parametric studies were performed to better understand the fluid flow characteristics of PCL in humans with ciliary disorders. Our results can be summarized as follows: (1) A complex region was generated by the ciliary carpet, and a substantial amount of flow is transferred in the direction parallel to effective stroke. Further, the visualization of the velocity field, particularly above the cilia tips, showed unidirectional transport of the ASL over time. This is in accord with the actual muco-ciliary mechanism. It was also found that vortex-like structures are formed above the cilia tips once they are moving backward during their recovery stroke. However, when the cilia move forward during their effective stroke, small vortices appear below the cilia tip.
(2) Our flow rate and surface velocity results imply the significant role of ciliary beat frequency, suggesting that the discovery of novel pharmaceuticals that can target cilia and increase their frequency of beating could enhance the PCL transport, and as a result, better transport of the ML under diseased conditions, such as asthma, or in cigarette smokers. It should be noted that increasing the beat frequency results in a high shear stress at the epithelium which may adversely affect the epithelial functions.
(3) It was shown that decreasing the ciliary length, which is the case in COPD patients and smokers, hinders the PCL flow rate and the surface mean velocity. It is also found that if the length of cilia is increased, it is less likely to see any significant changes in the PCL flow rate, while the time-averaged surface velocity would increase.
(4) In PCD patients with stiff planar and wiper ciliary beat patterns, pathogens are not effectively cleared from the respiratory tracts due to the large reduction of PCL flow rate, as well as the surface velocity. Nevertheless, the impaired ciliary beat patterns would not significantly change the shear stress.
(5) Generally, the triangular ciliary packs could outperform the square pack, thus increasing the pumping performance of engineered ciliated tissues in microfluidic chips. Note that the case of the equilateral triangle pack showed a higher fluid transport rate than the isosceles triangle pack, but less wall shear stress.
(6) If the patch size is increased, the flow rate, surface velocity, and the shear stress of fluid will be reduced accordingly. Thus, CF patients and cigarette smokers are more prone to infection from the loss of the cilia.
(7) It was found that the flow rate and surface velocity first increase sharply with increasing ciliary density (or number of cilia within the same computational domain) and then reach a plateau. On the other hand, the changes of shear stress depend highly on the number of cilia and did not level off at the higher ciliary density.
In the future work, the present model can be extended to study how drugs dissolved into the PCL would be absorbed by the epithelial cells, or how chemical species secreted by epithelial cells would be dispersed by PCL flows.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: