Free Interfaces at the Tips of the Cilia in the One-Dimensional Periciliary Layer

: Cilia on the surface of ciliated cells in the respiratory system are organelles that beat forward and backward to generate metachronal waves to propel mucus out of lungs. The layer that contains the cilia, coating the interior epithelial surface of the bronchi and bronchiolesis, is called the periciliary layer (PCL). With ﬂuid nourishment, cilia can move efﬁciently. The ﬂuid in this region is named the PCL ﬂuid and is considered to be an incompressible, viscous, Newtonian ﬂuid. We propose there to be a free boundary at the tips of cilia underlining a gas phase while the cilia are moving forward. The Brinkman equation on a macroscopic scale, in which bundles of cilia are considered rather than individuals, with the Stefan condition was used in the PCL to determine the velocity of the PCL ﬂuid and the height/shape of the free boundary. Regarding the numerical methods, the boundary immobilization technique was applied to immobilize the moving boundaries using coordinate transformation (working with a ﬁxed domain). A ﬁnite element method was employed to discretize the mathematical model and a ﬁnite difference approach was applied to the Stefan problem to determine the free interface. In this study, an effective stroke is assumed to start when the cilia make a 140 ◦ angle to the horizontal plane and the velocitiesof cilia increase until the cilia are perpendicular to the horizontal plane. Then, the velocities of the cilia decrease until the cilia make a 40 ◦ angle with the horizontal plane. From the numerical results, we can see that although the velocities of the cilia increase and then decrease, the free interface at the tips of the cilia continues increasing for the full forward phase. The numerical results are veriﬁed and compared with an exact solution and experimental data from the literature. Regarding the ﬁxed boundary, the numerical results converge to the exact solution. Regarding the free interface, the numerical solutions were compared with the average height of the PCL in non-cystic ﬁbrosis (CF) human tissues and were in excellent agreement. This research also proposes possible values of parameters in the mathematical model in order to determine the free interface. Applications of these ﬂuid ﬂows include animal hair, ﬁbers and ﬁlter pads, and rice ﬁelds.


Introduction
Breathing in polluted particles, such as dust and bacteria, is often unavoidable for people. At the same time, the innate immune system plays an important role in protecting the body by secreting mucus from goblet cells lining the airway epithelium to trap particles and remove them through the beat cycle of the cilia. This system is called muco-ciliary transport and has been studied by several authors [1][2][3][4][5][6][7][8][9]. Muco-ciliary clearance is fundamental for the health of the human respiratory system. Muco-ciliary dysfunction causes chronic airway diseases in humans, such as cystic fibrosis, primary ciliary dyskinesia, asthma, and chronic bronchitis. Moreover, deficient mucus clearance is not completely understood. For example, cystic fibrosis is a disease which depletes the volume of fluid or the height of the fluid in the respiratory system. Further investigations are vitally needed to investigate the effect of changing the volume or height of the fluid to improve muco-ciliary transport. In the respiratory system, cilia are contained in a layer called the periciliary layer (PCL), which consists of PCL fluid. Above the PCL is a layer of viscoelastic mucus. Figure 1 shows a portion of the patchwork of muco-ciliary clearance in the respiratory system.
At the epithelial cells, there are goblet cells, scattered among ciliated cells, which secrete mucus, a viscous fluid, to trap particulate matter and micro-organisms. This clearance prevents them from reaching the lungs. After this, the mucus forms a mucus blanket floating on a lower viscous fluid layer, the PCL, which contains the cilia covering the surface of the epithelial cells. Cilia beat back and forth about 12 times per second to propel mucus at 1 millimeter per minute [10], beating along the same plane without sweeping sideways during the effective stroke [11,12]. Above the mucous layer is air moving in and out of the human body. Because the height of the PCL contributes a significant role in the muco-ciliary clearance system, we first focus on the free surface or the height of the PCL fluid at the tips of cilia. Several studies modeled cilia movement in both experimentally and numerically. For example, Quek et al. [13] simulated the ciliary flow in a three-dimensional domain using the immersed boundary method. They concluded that the stiffness and the cilia density effected the slip velocity at the tips of cilia. Cilia in both high and low densities had a low slip velocity and the wave propagation was not supported by long stiff cilia. Sears et al. [9] used video microscopy to record ciliary motion. They converted planar ciliary motions into an empirical quantitative model, which was easy to use for examining ciliary function. If the tip of a cilium was modeled to penetrate the mucous layer, an integral equation could be used to express the mean velocity field in the mucous layer involving the tips of the cilia [2]. Fine-tuned methods of this kind can be found in [14][15][16][17].
Cilia were detected by several researchers. Some considered cilia and flagella to have the bending properties of an axoneme [18][19][20]. Den Toonder and Onck [21] studied and reviewed artificial cilia in many aspects, such as the bead-spring model for artificial cilia, a superparamagnetic filament for fluid transport, modeling the interaction of active cilia with species, electrostatic artificial cilia, microwalkers, artificial flagellar microswimmers, fluid manipulation by artificial cilia, and the measurement of fluid flow generated by artificial cilia.
Because the problem considered in this study is the one-dimensional prototype problem of finding the free surface at the tips of the cilia while the cilia move and drive the PCL fluid, instead of typical pressure gradient doing so, we first assume there to be a fluid-gas phase at the free surface. Because of the inference of the phases, the height of our domain, in this study, is proposed to change according to the angle θ that the cilia make with the horizontal plane, as shown in Figure 2. Figure 2 demonstrates the domain when (left figure) the cilia are perpendicular and (right figure) make a θ angle with the horizontal plane. The variable ζ is the length of a cilium and z is the axis. The heights of the domains in both figures depend on the tips of the cilia and the free surface that occurs at the distal ends of the cilia. To model the problem, since the PCL contains both PCL fluid and solid phases, i.e., cilia, it can be formed as a porous medium, see Figure 2. In this study, we consider the movement of bundles of cilia effecting the PCL fluid. Therefore, we use the Brinkman equation to predict the velocity of the PCL fluid in the porous region. Numerical and analytical studies can be found in the literature that focus on the fluid flow through the domain using the Brinkman equations [22][23][24][25][26]. For example, Morandotti [25] proved the existence and the uniqueness of the Brinkman model for the movement of a self-propelled microswimmer. Skrzypacz and Wei [26] proved the existence and uniqueness of a weak solution of the nonlinear Brinkman-Forchheimer-Darcy equation for the flow in chemical reactors of the packed-bed type. Liou and Lu [23] applied the Brinkman equation to flow over rough walls using the volume-averaging method. Cortez et al. [24] found the exact solution of the three-dimensional Brinkman equation where the flow was driven by regularized forces. Shavit et al. [22] suggested a modified Brinkman equation was a good fit with a free flow problem. There are a variety of numerical methods that study free boundary problems, such as the level-set method (LSM), the immersed boundary method, the simple line interface calculation (SLIC), the flux line-segment model for advection and interface reconstruction (FLAIR), and the volume-of-fluid method (VOF). For example, Ashgriz and Poo [27] initially developed the FLAIR method to predict the free surfaces of problems, originally taking the idea from Youngs [28] of using a sloped line in the interface reconstruction. Hirt and Nichols [29] established the volume-of-fluid method to extract the interface. Mashayek and Ashgriz [30] claimed that the volume-of-fluid method could handle large surface deformations. Although, there are several free boundary methods, to author's knowledge, there have been no studies conducted on the cilia movement in the PCL with a free boundary method. Because we assume there to be a fluid-gas phase at the free interface, we apply the Stefan problem [31,32] at the tips of cilia where the free surface is changed with respect to the angle θ, instead of time t. Because cilia tips touch the mucous layer only on the effective stroke [33], our objective was to investigate the free interface at the tips of the cilia for the forward stroke, provide feasible numerical results, and compare the results with the experimental data available in literature. The free interface obtained in this work becomes the lower boundary of the mucus. Our latter goal was to study changes in the sensitivity of the free interface of the PCL when the free surface of the mucus changes. If it is responsive, this indicates that we have to repeat our protocol to find the velocity of the mucus.
The mathematical model is provided in Section 2. To avoid the computational grid adjustment for each angle, we applied the boundary immobilization technique to the model in Section 3. Model discretization is presented in Section 4. The validation of the numerical results and numerical solutions is provided in Sections 5 and 6, respectively. The conclusion is drawn in Section 7.

Mathematical Model
In order to consider the motion of cilia collectively, rather than that of a single cilium, the governing equation employed in this work is the Brinkman equation on the macroscopic scale. They were upscaled using the hybrid mixture theory (HMT) [34][35][36][37], i.e., by applying the averaging theorem to average the microscale field equations to obtain the macroscale equations. In this work, we consider the model in steady state. That is [37,38] where l is the porosity; v l and v s are the velocities of the liquid and solid phases, respectively; µ is a dynamic viscosity; k −1 is the inverse of the permeability tensor; p is pressure; ρ is fluid density; g is gravity;˙ l is the material time derivative of the porosity with respect to the solid phase, Therefore, the system of Equations (1) and (2) can be rewritten as Figure 3 illustrates the fluid flow at the tips of the cilia, where the domain containing cilia is considered as a porous medium. In this work, we assume that the pressure gradient is zero in the z direction and a constant in the x direction. The velocity and the other variables depend on the z direction. Substituting (5) into (4), we obtain (4) For the boundary condition between the free fluid and porous medium interface, we were motivated by [39], where they used Beaver and Joseph boundary condition at the top of a porous region. We then use the boundary condition where r is a continuous function depending on the z−axis and the angle θ; q is a function of the angle θ at the free surface; and c 3 is a constant. We also employ the Stefan problem [31] at the free-fluid/porous-medium interface with the rate of change of the interface with respect to the angle θ.
That is with the initial condition where c is a constant and θ 0 is the initial angle.

Boundary Immobilization Technique
Because of the free boundary at the free-fluid/porous-medium interface, we employ a boundary immobilization technique [32] to the problem. Define the spatial coordinate transformation as to fix the moving surface to a fixed coordinate system in space, where η is the dimensionless variable, z is the spatial coordinate, and q(θ) is the function of the angle θ at the free-fluid/porous-medium interface. Therefore, at z = q(θ), η = 1. Thus, Applying partial differential relations to the velocity v, we have Substituting the second derivative (13) into (6), we obtain Similarly, the boundary and initial conditions become To calculate the numerical results, model discretization is provided in the next section.

Model Discretization
In order to find the numerical solutions, the discretization of the Brinkman equation is provided in this section. We apply a finite element method to the system of Equations (14)- (18). Using the Galerkin finite element method, where the test function w comes from a chosen trial function [40], the weak formulation of the one-dimensional governing Equation (14) is to find v ∈ H 1 (Ω), such that where Ω is our computational domain; H 1 (Ω) is the Hilbert space; and w ∈ H 1 0 (Ω) is the test function. The n element discretization of (19) is Define the finite-dimensional subspaces of the Hilbert space H 1 (Ω) as follows where T h is a single line segment of the domain. The approximate solution where V is the vector of the velocity; ψ m is called a basis function; Ψ is its vector form; and the integer M is determined by the interpolation function. For a linear shape function, M = 2. In this study, we use linear isoparametric element, with where ξ is the variable in natural coordinates as shown in Figure 4, which illustrates the mapping between the physical (η) and natural (ξ) coordinate systems using the shape functions (23) and the interpolation (24). This is effective for numerical integration because the linear element is normalized between −1 and 1 in the natural coordinate system. Therefore, where h = η i − η i−1 is the distance between adjacent points of a uniform mesh. Hence, the vector of the velocity and the vector form of the basis function are respectively. Substituting the vector form of the basis functions, Ψ in (26) into w in (20), and the approximate velocity v in (22) into (20), we obtain the element matrix form of (20), which is where the element matrices, in the natural coordinate, are where the variable j = dη dξ is the Jacobian, 0 ≤ η i ≤ 1 and is the vector of the solid velocity in each element. After assembling the element matrices, we have the global matrix form, where we use the same notation V for both element and global vectors of the velocities for the sake of notation simplicity and the boundary term in (20) is added into the last component of the vectorF. Let Because the velocity is zero at the bottom of the domain, Since ψ 1 = 0 and ψ 2 = 1 at η = 1, Therefore, where the second and third equalities come from (15) and (32), respectively. In this study, we employ the porosity l and the velocity of cilia from [41], where they obtained the solid velocities from [9], which are the eighth-order polynomial functions as shown in Table 1 in [41]. Therefore, and where ζ = z/ sin θ as shown in Figure 2. To discretize (17), we first change the variable q to be u using (32) and find its derivative with respect to θ. That is Substituting (40) into (17), we obtain A suitable discretization formula of (41) is [32] where u at an initial angle is where ∆θ is the difference between two adjacent angles θ and N is the number of nodes of the numerical porous-medium domain. The validation of the numerical solutions and numerical results are presented in the next sections.

Numerical Validation
The numerical result of the discretized model is validated in this section. The exact solution of (14) is provided in [42] with a linear velocity of cilia, i.e., v s = z, and a fixed boundary conditions v(0) = 0 and v(1) = 100, which is where The values of variables are µ = 3 × 10 −6 , g = −1, ρ = 992.2 × 10 −15 , k = 0.0017, l = 0.74867, and dp/dx = −9.7335 × 10 −15 , where k and l are set to be the values that the cilia are perpendicular to the horizontal plane. Figure 5 shows the exact and numerical solutions with 10, 20, 100, and 400 elements. The l 2 −norm errors are presents in Table 1. They illustrate that the errors decrease when the number of elements increases, showing the convergence of the numerical results to the exact solution. Notice that the velocity rapidly changes at the upper part of the domain, which is close to z = 1. The numerical result is different from the analytical result, but the difference is not observable from the graph when the number of elements is 400. This example problem is well posed [38] and is independent of time, so the numerical results were expected to be very accurate.

Numerical Results
The velocities of the PCL fluid due to the forward stroke of the cilia are illustrated in this section. We assume that the cilia start at θ = 140 • [9] and move forward with a decreasing angle, and then the cilia stop at θ = 40 • . In this study, we calculate the PCL velocity when the angle decreases by 10. The velocities of the cilia, produced from data presented in [9], at each angle θ are shown in Table 1 in [41] which are plotted in Figure 5 in [41]. In this work, we arrange the velocities of the cilia at the angles θ = 130 • , 120 • , 110 • , 100 • to equal θ = 50 • , 60 • , 70 • , 80 • , respectively. To find the numerical results of the Brinkman equation, the initial velocity is assumed to occur at θ = 140 • , which is defined to be an increasing function from 0 to a value which is less than 10 µm/s along the length of cilia.
We place the initial velocity in (42) and then calculate u θ+∆θ . Then, u at the angle θ + ∆θ is placed in (27). After assembling the element matrix, we obtain the velocity of the PCL fluid. Then, it becomes an initial velocity in (42). The iterative process is continued until we obtain the velocity of the PCL fluid for every angle θ.
The values of variables used in the program are the same as in Section 5, i.e., µ = 3 × 10 −6 g/(µm·s), ρ = 992.2 × 10 −15 g/µm 3 , dp/dx = −9.7335 × 10 −15 g/(µm 2 s 2 ). Regarding the other variables, the gravity g = −9.81 × 10 6 µm/s 2 , c = 0.035, r = e θ , c 3 = d(1 − z 3 )/U 0 , where d = 7.5 is the greatest length of cilia; U 0 is the maximum velocity of the cilia; and z 3 = e z s , where z s is the greatest vertical length when cilia make an angle θ with the horizontal plane. Because of the dimensionless coordinate, the values of z s are less than 1 and change depending on the angle θ. The values of permeability k and the porosity l were obtained from [41,43], respectively. Although Chamsri and Bennethum [43] found the permeability in a three-dimensional domain, in this work, we only adopt the permeability in the x direction, which is k(1, 1) in [43]. This is because from Figure 13 in [41], the velocity in the x direction affects the average velocity of all directions the most, which has an essential role in the fluid flow above the cilia. Moreover, the diagonal element of the permeability highly influences the fluid flow.
With 200 elements, Figure 6 shows the velocities of the PCL fluid from the angle 130 • to 90 • . Notice that the velocities of the PCL fluid increase when the angles decreases. Figure 7 shows the velocities of the PCL fluid from the angle 80 • to 40 • . In this forward movement, the velocities of the cilia are slowing. Therefore, the velocities of the PCL fluid decrease as the angles decrease. Figure 8 shows the combination of the velocities of the PCL fluid in Figures 6 and 7. This illustrates that the velocities increase until the 80 • angle, then they decrease from the 70 • to 40 • angles.    Figure 10 shows the free surface, q values, for the forward stroke at the tips of the cilia from the 140 • to 40 • angle. Although the angle becomes smaller, the free surface exhibits an increasing behavior. The numerical results are compared with the results in [45]. They measured the height of the PCL in human lower airway biopsies and in mice. The average height of the PCL in non-cystic fibrosis (CF) tissues was 5.60 ± 0.28 µm, and in CF tissue, it was 4.52 ± 0.47 µm. In mice, the non-CF PCL height was 4.70 ± 0.13 µm and the CF height was 4.10 ± 0.09 µm. They stated that the PCL fluid was reduced in patients with cystic fibrosis (CF). From Figure 10, we can see that the average of the free surface is 5.5858, which is close to the height of the PCL in non-CF human tissues. Figures 11 and 12 present the free surfaces for different variables r, which are e θ and e −θ , respectively, with values of c = 0.035, 0.045, ..., 0.085, where the angles are the same as those shown in Figure 10. The shapes of the free surfaces for different functions r are different but they are increasing functions of the angle θ for the forward stroke, and for each r, they have the same shape for different values c. The height of the free surface moves up when the initial values c increase.

Conclusions
We determine the free surface at the tips of the cilia due to the movement of cilia using the Brinkman equation on a macroscopic scale with the Stefan problem. In this work, we consider the free boundary between the phases, which moves with the change in the angle θ rather than time. We summarize the process and the results as follows.

•
To avoid the numerical grid adjustment for every change in the the angle θ, we adopted the boundary immobilization technique to the mathematical model; • The one-dimensional isotropic finite element method was applied to the Brinkman equation and a finite different scheme was employed to the Stefan problem to calculate the numerical solutions; • The numerical results were verified by comparing them with an exact solution with a fixed boundary condition at η = 1 when θ = 90 • . The results converge to the exact solution, where the l 2 -norm errors are provided in Table 1; • The numerical results of each angle θ are provided in Figures 6-8 for the forward stroke from the angles 130 • to 40 • . Even if the velocity of the cilia is highest at θ = 90 • ( Figure 5 in [41]), the highest velocity of the PCL fluid is found at θ = 80 • . That is, the velocities of the PCL fluid increase from the starting angle and decrease from the angles 70 • to 40 • ; • The free surfaces due to the movement of cilia are shown in Figures 10-12. We discovered that the height of the surface increases for the full forward stroke although the velocities of the PCL fluid increase in the first half of the effective stroke and decrease for the second half. Moreover, we found that the initial value c does not change the shape of the free surface, but the different functions r cause the change in the shape of the surface as illustrated in Figures 11 and 12; • To the author's knowledge, no experimental data in the literature provide the shape of the free surface, only the average value is given. The numerical results were compared with the available experimental data. For the non-CF tissues, the average heights of the PCL from our results and the experimental data were 5.58 and 5.60 ± 0.28 µm, respectively. The difference between 5.60 ± 0.28 and 5.58 is 0.02 ± 0.28, which is small. This is one of the justifications that shows that our mathematical model is suitable for this problem.
Even though the mathematical model was developed in [32,37,38], in this study, we determined the free interfaces at the tips of cilia during the full forward stroke. From this, we hope to be able to further assess the natural velocity of mucus in future work. This study is a benchmark in terms of assessing the height and shape of the free boundary at the tips of cilia. In forthcoming work, we hope to determine the free boundary of the mucous layer, in which the result from this study will become the lower boundary of the mucous region. Consequently, we hope to extend the problem into the twoand then three-dimensional domains.