Simulation Study on the Mass Transport Based on the Ciliated Dynamic System of the Respiratory Tract

To study the mass transport of mucociliary clearance of the human upper respiratory tract, a two-dimensional mass transport model based on the ciliated movement was established by using the immersed boundary-lattice Boltzmann method (IB-LBM). In this model, different characteristics of the mucus layer (ML) and the periciliary liquid (PCL) were taken into account. A virtual elastic membrane was introduced to divide the two layers dynamically. All moving boundaries that were involved in the present simulation were modeled with the immersed boundary. The Newtonian fluid was used to model the flow in PCL, and the viscoelastic fluid based on the Oldroyd-B model was used for the flow in ML; the two types of flow were both solved by the LBM framework. Based on the model, the ML thickness, the cilia density, and the phase difference of adjacent cilia were regulated, respectively, to study the transport velocity of the ML. In addition, the motion law of solid particles in PCL was also studied. According to the results, four primary conclusions were drawn. (1) At a given beating pattern, the increase of the ML thickness will decrease its transport velocity. (2) Increasing the cilia density can promote the mean transport velocity of the ML. (3) By raising the phase difference of adjacent cilia to a certain scope, the transport of ML can be accelerated. (4) In PCL, particles initially located on the upper part of the cilia tend to migrate upward and then get close to the ML. The above study can provide some reasonable explanations for the mechanism of the mucociliary clearance system, which is also helpful to the further understanding of the mass transport principle of the human upper respiratory tract.


Introduction
Mucociliary clearance (MCC) system is an important part of the human immune system, which is responsible for the removal of bacteria, viruses, dust, and various harmful particles in the inhaled air. It plays an important role in protecting human health [1]. e MCC mainly consists of cilia and mucus attached to the respiratory tract surface. e mucus is made up of two layers, i.e., the mucus layer (ML) and the periciliary liquid (PCL). ML is a non-Newtonian fluid layer of high viscosity. It can capture and filter harmful particles in the air. PCL is a Newtonian fluid layer of less viscosity, which has similar physical properties to the water. e PCL is important to lubricate the surface of the respiratory tract and facilitate the swing of the cilia [2]. By performing a coordinated spontaneous beating, cilia promote the movement of ML to the throat and then transport foreign materials out of the body [3].
ere are many factors affecting the mass transport of MCC, such as human activity, drug stimulation, pH level variation, and external environment [4]. Other studies reveal that the particulate matter 2.5 (PM2.5) can change the secretion of respiratory mucus and even destroy the structure of cilia, resulting in abnormal cilia function and dysfunction of ML transport [5]. In addition, some artificial traumas and congenital diseases, such as surgery, primary ciliary dysplasia (PCD), and immobile cilia syndrome (ICS), can affect the structure or the movement state of cilia and then affect the work efficiency of MCC [6,7]. erefore, the study of cilia movement mechanism is practically significant. However, up to now, the relevant experiment study is insufficient for two major reasons. e first is that the respiratory structures are complex, and the cilium contains various microscopic motion patterns. e other is that the appropriate noninvasive experimental equipment is still lacking in application.
In recent years, with the rapid development of computer science and numerical technology, computer simulation of MCC has attracted more attention and the mass transport of MCC has become an attractive topic in the fields of both physiology and physics [8]. To model the beating cilia, the immersed boundary method (IBM) is one of the most popular research methods. Peskin introduced this method for the first time to study the cardiac blood flow around the heart valve [9]. Now, it has been widely used in solving biological problems related to moving boundaries [10]. By combining IBM with the finite difference method, Lee et al. [11] established a two-dimensional (2D) model to study the underlying mechanisms of some diseases in the MCC system. ey concluded that the beat frequency and density of cilia and the depth of PCL were the key factors affecting the transport velocity of the ML. Jayathilake et al. [12] studied three cilia beating patterns: windscreen wiper motion, rigid planar motion, and normal motion. It was found that abnormal cilia beating patterns, such as windscreen wiper motion and rigid plane motion, greatly reduced or prevented the transportation of ML. Moreover, Jayathilake et al. [13] developed a 3D model to simulate the movement of human lung cilia in a single layer of mucus. e results showed that thinner PCL or higher cilia beating frequency could increase the mean flow velocity of PCL. By applying the immersed boundary-lattice Boltzmann method (IB-LBM) framework, Shahmardan et al. [14] studied the effect of mucus thickness on its transport efficiency. e results showed that an increase in mucus thickness caused by air pollution and smoking reduces the mean mucus velocity. In particular, for the first time, Sedaghat et al. [15,16] used IB-LBM to solve the two-layer mucus model. e viscoelastic flow was established with the Oldroyd-B model in the upper layer ML, and the Newtonian fluid was set in the lower layer PCL.
eir study revealed that the viscosity had a significant effect on mucus transport efficiency.
Although great progress has been made in the study of the MCC system, there are still some principles in MCC that need to be discovered. On one hand, Hoffmann and Cortez [17] studied the motion of two adjacent cilia; they concluded that the speed of ML tended to increase in the out-of-phase beating. We think that, in this topic, the phase difference of adjacent cilia on ML transportation deserves further study. On the other hand, besides the cilia and flow, some floating rigid particles are usually found in MCC [18], which may be the foreign matter or come from the respiratory tract surface. Moreover, some other topics, such as the decontamination mechanism of cilia [19] and drug effect on the motion of cilia [20], are also worth further exploring.
In this study, the Oldroyd-B model was used to model the viscoelastic fluid in ML. A 2D IB-LBM was applied to simulate the transport velocity of ML in the MCC system. A virtual membrane was used to separate the two layers of ML and PCL dynamically. Compared to the model in References [15,16], we improved the motion control of cilia. A moving virtual cilium was introduced to lead the motion of the actual cilia through elastic springs. To do this, we can obtain anticipant beating patterns of the cilia. Based on the model, by changing the thickness of ML, the density of cilia, and the phase difference of adjacent cilia, the transport efficiency of the ML was studied in detail. Moreover, we also conducted a study on the particle motion principle in PCL, where the particle size and initial position on the particle migration were discussed. e rest of this paper is organized as follows. Section 2 gives the physical model and the numerical method. Section 3 verifies our models. Section 4 exhibits and discusses the numerical results. And the final conclusions are listed in Section 5.

Physical Model.
e diagrammatic sketch of the physical model is displayed in Figure 1.
e left and right boundaries of the channel are set as periodic.
e top boundary is a slip wall. e bottom boundary is a nonslip wall, on which the root ends of cilia are fixed. Similar to the model in References [15,16], the mucus is divided into two layers, ML and PCL. ey are divided by a virtual elastic membrane. e fluid is viscoelastic in ML and is Newtonian in PCL.
In order to measure the transport efficiency of the ML, as shown in Figure 1, we set a group of tracepoints to float with the flow. e mean velocity of ML is estimated by calculating the mean velocity of these tracepoints within a time scope of 3T. Worth noting is that, in the chosen time window, the flow in ML has been fully developed.

Mathematical Models of Flow and
Structure. In this model, the cilia are immersed in a viscous fluid; their autonomous beating motion drives the flow of the surrounding mucus. e flow is governed by the Navier-Stokes equation: where u is the velocity, μ is the kinetic viscosity of the fluid, p is the pressure, and F is the resultant force on the cilia. To divide ML and PCL, the cilia and the membrane consist of a set of nodes connected by springs in a consecutive way [21,22]. e force on cilia nodes through the spring is made up of three components, the stretching force F s , the bending force F b , and the driving force F c from the virtual cilia. en we get the resultant force F � F s + F b + F c . For the membrane, each node involves two force components of F s and F b . en there is F � F s + F b . e stretching force F s satisfies Hooke's law of elasticity in the tangent direction. It maintains the inherent length of the boundary, which can be written as where K s is the extensional coefficient. e bending force F b is the bending moment of the boundary in the normal direction, which is derived from the Frechet derivative of the bending energy formula based on the virtual work principle [23]. It is calculated by where K b represents the bending coefficient. In addition, the moment and transverse stress are zero for nodes at both ends of the cilia [24]: In this model, in order to realize the autonomous beat of cilia, a virtual cilium with a fixed beat posture is introduced to lead the moving of the actual cilium boundary. A group of virtual springs are used to connect the virtual and the actual cilium. ey are displayed in Figure 2. e elastic force F c is defined as F c � K sl (Y − X), where K s1 is the corresponding elastic coefficient and Y represents the boundary of the virtual cilium. A hard particle model was used in the present study, which follows Newton's second law. e velocity of the particle is made up of two components, i.e., the translational velocity U p and the angular velocity Ω p ; the corresponding forces on the particle are formulated as [25] where the subscripts p and f, respectively, stand for particle and fluid. m p is the mass of the particle. ρ is the density. I p is the moment of inertia. X s − X c stands for the position vectors of a Lagrangian point and the center of particle, respectively. CS stands for the surface of the particle, and n is the normal vectors of triangles. σ s is the fluid stress tensor, and g is the gravitational acceleration. Moreover, the repulsive effects between two boundaries also need to be considered, including the cases of particle surface to particle surface, particle surface to cilium boundary, and particle surface to wall boundary. e repulsive force is marked with F rep ; it is calculated according to the model suggested in Reference [26]. en the total force on the particle is In the Oldroyd-B model for the viscoelastic fluid in ML, the pressure tensor σ is decomposed into two parts as [27] where σ M,N is the part of the elastic tensor in the Newtonian fluid, and σ M,E is the part of the elastic tensor in the non-Newtonian fluid, where σ M,N is given as where D is the strain rate tensor and it is defined as For σ M,E , it is derived from the Upper Convected Maxwell model (UCM), which is given as where λ is the relaxation time, which accounts for the magnitude of hydroelastic properties. σ In addition, the viscosity in this model can be divided into the Newtonian part (η M,N ) and the elastic part (η M,E ). e total viscosity η M of the fluid is η M � η M,N + η M,E . ere are several methods to solve the model to get the tensor σ. In the present study, the method based on LBM in Reference [28] is used to solve the Oldroyd-B model.

e IB-LBM Framework.
In our study, the D2Q9 model of the lattice Boltzmann method is employed to solve the flow. e single relaxation time lattice Boltzmann equation is [29,30] where g i (x, t) is the distribution function of particles, e i is the particle velocity, ∆t is the time step, g eq i (x, t) is the equilibrium distribution function, τ is the nondimensional relaxation time, and G i shows the effect of volume force on the distribution function.
In equation (11), g eq i (x, t) and G i are given as [31] where f is the vector of the body force density, and ω i is the weight defined by ω 0 � 4/9, ω i � 1/9 for i � 1 to 4, and ω i � 1/ 36 for i � 5 to 8.
Δt is the sound speed. In this study, the immersed boundary method is applied to deal with the interaction between cilia and fluid. e equation of IBM is given by where F(s, t) is a Lagrangian force that acts on the fluid by cilia, and Dirac's delta function . e velocity on cilia can be obtained by interpolation of fluid points nearby. To update its location, the velocity of cilia U(s, t) is given as

The Verification of the Oldroyd-B Model
In this section, a 2D Poiseuille flow is used to verify the numerical method for the Oldroyd-B model. e model is shown in Figure 3, in which the channel consists of two parallel planar no-slip walls. e grid-scale L × H is 120 × 40. A parabolic velocity profile is set at the inlet, and the free flow condition is set at the outlet. e nonequilibrium extrapolation scheme is utilized to model the microboundary conditions in LBM.
As the viscoelastic Poiseuille flow is fully developed, the analytical results of velocity and configuration tensors can be expressed by the following formulas: where u max is the maximum speed at the entrance, and L y is the characteristic length of the geometric model, which is the height of y-direction. e Weissenberg number W e is set as 0.01 and 1. e Reynolds number Re � 4, and viscosity ratio Rv � η M,N /(η M,N + η M,E ) � 0.9. e maximum velocity at the inlet is U max � 0.1.
From Figure 4, it can be concluded that numerical and analytical results are very close at W i � 0.01 and 1.
is indicates that our model can simulate the viscoelastic flow correctly.

Results and Discussion
e relevant parameters are set as follows. e side length of the grid of LBM is ∆x � ∆y � 1, corresponding to 0.1 μm in the physical unit. e size of the channel is 600∆x × 102∆y. e Reynolds number Re � 0.05. e phase difference of the adjacent cilia is 0.02T (T is the beating cycle of cilium). e depth of PCL is fixed at 6.2 μm [32]. e viscosity ratio of ML to PCL is 40. For the flow in ML, the Weissenberg number W e � 0.01, and the viscosity ratio of the Newtonian part to the total viscosity is η M,N /(η M,N + η M,E ) � 0.975 [32]. e cilia beat frequency is 20 Hz [33]. e total simulation time is set as 30T, where T is the beating cycle. e cilium is structured by the node-spring model, and the default spacing of the adjacent nodes is ∆s � 1. e length of the cilium is L � 55∆s (5.5 μm); its beat frequency is set to be 20 Hz [33].
e extensional coefficient K s in equation (2) is set as 2.0 in nondimensional unit. is can restrict the stretching rate within the range of ±1% in simulation. e bending coefficient K b in equation (3) is set as 1.0; this can make the cilium have certain bending rigidity and then perform a stable beating. e extensional coefficient K sl for the elastic force F c is set as 1.0, which can make the actual cilium move with the virtual cilium synchronously. For the membrane to divide the ML and PCL, it has the same settings of K s and K b to the cilium.
Based on the above settings, the beating amplitude was regulated at the case that the distance between two adjacent cilia is 0.55L [16], where we obtained the mean velocity of ML to be about 42 μm/s, which is in good agreement with the results of other studies in the literature [32,34].

Effect of the ML ickness and Cilia Density.
For the single PCL layer model in Reference [13], the effect of PCL depth on the transport efficiency has been discussed. It is known that if the depth of PCL increases, the transport velocities of the stream will increase at a decreasing rate. In the present study, a double-layer (PCL and ML) model is applied; it has been found that the transport efficiency is different for each layer. Here, we fix the depth of PCL as 6.2 μm and then change the thickness of ML and the density of the cilia, based on which the transport velocity of ML is studied. To conduct the simulation, a thickness sequence of ML is set from 2 to 10 μm with an interval of 1 μm. Meanwhile, the density of the cilia is regulated by controlling the distance of two adjacent cilia, where three patterns are set as 0.55L, 0.73L, and 0.91L. e results are shown in Figure 5.
As shown in Figure 5, it is found that, in the three patterns of cilia density, the mean mucus velocity decreases as the ML thickness increases. However, in the single-layer model, the transport velocity will also increase if the depth of PCL is increased [13]. We think that the primary reason for this phenomenon is the effect of the viscoelastic fluid in ML. Two detailed reasons are analyzed below. First, the viscosity of the flow in ML is much higher than that of PCL, and the flow has less flowability. Second, the elastic property of the fluid can also affect the flow state in ML, which makes the fluid have some characteristics of the solid. is can hinder the transport of ML.
In addition, it is also found that, by increasing the cilia density (decreasing the distance between two adjacent cilia), the transport of ML will accelerate. at is, a higher density can result in faster transport of ML. is is qualitatively in accord with the results from Reference [11]. In order to exhibit more details about the transport velocity under different cilia density, a distance sequence between two adjacent cilia is set as 1.09L, 0.91L, 0.78L, 0.69L, 0.6L, and 0.55L. And three patterns of ML thickness, 4 μm, 5 μm, and 6 μm, are also considered. e results are shown in Figure 6.
ese simulation results show that if the ciliary density is decreased, the velocity of the ML will also decrease. Moreover, the cilia density has an almost linear relationship with the mean velocity of the ML. is implies that smokinginduced ciliary dysfunction and even ciliary death may hinder the excretion of foreign materials by the mucociliary clearance system.

Effect of the Phase Difference.
We know, in most cases, the cilia beat asynchronously, where there is a phase difference between the adjacent cilia. However, how this phase difference works on the transport efficiency of ML is not clear now; it deserves a further study. In this section, a sequence of phase difference between adjacent cilia is set to study the transport efficiency of ML, 0.03T, 0.05T, 0.075T, 0.1T, 0.133T, 0.167T, 0.2T, and 0.25T. Meanwhile, three cases of cilia density are also taken into account, with a distance sequence of 0.55L, 0.73L, and 0.91L between two adjacent cilia. Worth noting is that our channel is periodic; the first beating cilium at the left inlet must keep a proper phase difference with the one at the right outlet. To meet this condition, the length of the channel must be lengthened or shortened to some extent. is operation will not influence the transport efficiency of ML because the number of cilia within a unit length of the channel is constant. e results are shown in Figure 7.
From Figure 7, we know that, in the scope of 0 to 0.25T, the velocity of the ML tends to increase by raising the phase difference. In the ranges of 0.03T to 0.05T and 0.133T to 0.25T, the increasing trend is relatively slow, and from 0.05T to 0.133T, the increasing trend is relatively fast. is indicates that the synchronous beating of the cilia is not advantageous to the transport of ML, and certain phase differences can promote transport efficiency. Actually, it can be deduced that if the phase difference is continuously raised to more than 0.25T, the transport efficiency will not be enhanced anymore.

Particle Migration in PCL.
In general, some floating particles can be found in MCC. ey may be the foreign matter from the air or some fragments from the respiratory tract surface. We know that the flow in ML has less flowability than that in PCL, so the movement of the particle in the ML is relatively simple, where the particle tends to float with the flow of the ML. In contrast, in PCL, the flow has lower viscosity and more flowability; furthermore, driven by beating cilia, this can generate a more complex flow, and the particles tend to move in complicated patterns. In such a case, to study the movement law of these particles can help us understand the clearance mechanism of the impurity in the respiratory tract.
In the present study, 20 hard particles were placed in PCL, where the diameter of the particle was set as 0.4 μm.
eir initial distribution is shown in Figure 8(a). ey are at a y-coordinate of 4.5 μm away from the bottom. With the time passing on, the snapshots of the particle distribution at t � 2T and 20T are picked up and displayed, respectively, in Figures 8(b) and 8(c). As a result, it has been found that almost all particles move to the top area of these cilia and get close to the ML. is indicates that the beating cilia can lead to a targeted motion of these particles, in which the particles may be captured by ML.
Furthermore, a more detailed study was conducted by varying the size and the initial position of these particles. Here, three particle diameters of 0.4 μm, 0.6 μm, and 0.8 μm are considered for simulation. Meanwhile, the initial y-coordinates are set, respectively, as 2 μm, 2.5 μm, 3 μm, 3.5 μm, 4 μm, 4.5 μm, 5 μm, and 5.5 μm. In each case, 20 hard particles are placed with the distribution like that in Figure 8(a), and the mean y-coordinate of all these particles is set to express the overall trajectory of particle motion. e total simulation time is set as 50T.
e results are shown in Figures 9(a)-9(c), from which we know that, for the three types of particles in the different sizes, if the initial y-coordinate is larger than 4.5 μm (including 4.5 μm), the particles tend to migrate to the top area of the beating cilia and get close to ML within 5T. Once the particles reach the top area of the beating cilia, they tend to move along with ML stably. In contrast, particles with other initial positions are not able to move to the top area of the Computational and Mathematical Methods in Medicine cilia; they just pace up and down in a confined area between two adjacent cilia.

Conclusion
In this paper, we employed the IB-LBM framework to model a 2D human respiratory mucociliary clearance system. Based on the model, the ML thickness, cilia density, and phase difference between adjacent cilia were varied to study the transport velocity of ML. And the migration law of particles with different sizes in PCL was also investigated. According to the simulation results, two conclusions can be summarized as follows. On one hand, for the same beating pattern in a fixed depth of PCL, the mean velocity of ML will increase by reducing the ML thickness, raising the cilia density, or enlarging the phase difference between adjacent cilia within a certain range. On the other hand, in PCL, the migration law of hard particles with different sizes is related closely to their initial position. Particles initially located on the upper part of the cilia tend to migrate upward and get close to the ML. In summary, the simulation results are reasonable, and our study can help people further understand the physical and physiologic function of the human respiratory system.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.  Computational and Mathematical Methods in Medicine