Simulation study on the mass transport in PCL based on the ciliated dynamic system of the respiratory tract

In this paper, a two-dimensional dynamic model of human respiratory cilia is established by immersed boundary-lattice Boltzmann method. In the model, periciliary liquid and mucus layer with different flow characteristics are considered, in which the non-Newtonian fluid in the viscous layer is described by Oldroyd-B model; the area between two layers is dynamically divided by a virtual membrane. By changing the cilium beat amplitude and the depth of periciliary liquid, this paper focuses on the removal of particle impurities by cilium beat. The results show that when the initial position is near the top of cilium, the particle tends to migrate upward to approach the mucus layer. However, if the initial position is near the middle or at lower part of the cilia, most of the particles will not be transferred to the vicinity of mucus layer by the swinging cilium. It is also found that by increasing the beat amplitude of cilia, particles are more likely to approach the mucus layer. By conducting this study, the migration of the particle in periciliary liquid is discussed in detail. The results are significant to help us to know the material transport law of human respiratory tract.


Introduction
Muco-ciliary clearance (MCC) system plays an important role in protecting human health. It is responsible for the removal of bacteria, viruses, dust and various harmful particles in the air from the lung [1]. The MCC consists of cilia and mucus attached on the epithelial surface of the airway. Upper respiratory mucus can be divided into upper mucus layer (ML) and periciliary liquid (PCL). Cilia promotes the movement of ML to the throat through the coordinated spontaneous beat and then transports foreign materials out of the body through swallowing or coughing [2,3].
With the rapid development of computer science and technology, computer numerical simulation of MCC has become a hot topic. Previous research can be found from Smith et al, 2008 [4]. Since then, researchers have used different numerical methods to research different aspects of MCC. W.L. Lee [5] established a 2D, two-layer Newtonian fluid model using a mix of immersion-boundary method(IBM) and finite difference methods to understand the underlying mechanisms of some respiratory diseases, such as cystic fibrosis. They concluded that the cilia beat frequency, the density of cilia, and the depth of PCL were the key factors affecting the transport speed of the mucus ML. Sedaghat [6] used Lattice Boltzmann method (LBM) combined with the finite difference method and IBM to solve the two-layer mucus model. For the first time, they took Oldroyd-B model as the constitutive equation of mucus. The results showed that compared with cilia beat frequency, the viscosity had a more significant effect on mucus transport efficiency.
Franz Hoffmann [7] placed two cilia in a 3D long channel and found that the speed of ML increased when the phase difference between the two cilia was generated. Alapati [8] confirmed that when there were large particles in the mucus, the position of other particles and the rotation of particles would affect the motion state of particles. Tripathi [9] explored the decontamination ability of cilia. Passive cilia have the ability of self-cleaning, and it can use the oscillation in the flow to prevent the attachment of particles and biological cells.
In this study, Oldroyd-B model was used as the constitutive equation for mucus, and a 2D model was established by using IB-LBM to simulate the transport velocity of mucus with given cilia beat pattern in two-layer mucus. A virtual membrane is used to separate the two layers of mucus with different viscosity, and the oscillation of the actual cilia is directed by a virtual cilium with a fixed motion trajectory. By changing the depth of PCL and oscillation amplitude of cilia, the effects on the transport of particle in different positions of PCL was simulated.

Physical model
In order to establish the muco-ciliary transport module in a long channel, the diagrammatic sketch of this model is displayed in Fig. 1. The structure of a cilium has two parts: one part is 55 nodes (which means L = 55∆x), and the other part is the virtual spring between the nodes.
For standard parameter sets, the number of cilia is set at 20. The cilia length is set to 5.5μm, and the distance between adjacent cilia was 3μm [10]. The standard height of PCL is 6.7 μm and the depth of ML is 4 μm [11]; the viscosity of PCL η p =0.275, the viscosity of ML η M =11. We assume that the standard viscosity of Newtonian part of the mucus is the same as that of PCL( η M,N =0.275). The standard viscosity of the elastic part is η M,E =10.725, and the viscosity ratio is β=η M,E /η M =0.975. The cilia beat frequency is 20Hz/s [12].

Numerical Method
In this model, the cilia are immersed in a viscous fluid. Cilia's autonomous beat drives the flow of surrounding mucus. The flow is governed by Navier-Stokes equations [13].
The cilia consist of a set of nodes connected by springs in a consecutive way [14,15]. The resultant force F on each cilia node consists of three forces, including stretching force F s , bending force F b [16], and tension of virtual cilia F c (s, t) (s, t) (s, t) (s, t) The physical cilia and the virtual cilia are associated by a virtual spring with stiffness, which is showed in Fig. 2. Apply F p to the nodes on the surface of the particles when building PCL particles. F p is the force which the particles are subjected to, it is calculated by: Where F ps and F pb are the stretchings and bending forces of the adjacent nodes on the particle surface, F pp is the elastic force between the particles, and F pc is the interaction force between the particles and the cilia.
In this study, Oldroyd-B model will be used as the constitutive equation of mucus [17]. Fluid is assumed to be discrete lattice in LBM [18,19]. In addition, IBM is applied to deal with the interaction between cilia and fluid [20].

Results and discussion
(a) (b) (c) Figure 3. The trajectory of particles in the vertical direction of PCL.(a) particle initial positions on Y direction is 4.5μm, (b) particle initial positions on Y direction is 5μm, (c) particle initial positions on Y direction is 5.5μm In this section, we change the depth of PCL and oscillation amplitude of cilia to study the effects on the transport of particle in different positions of PCL layer. For the first time, the depth of PCL was taken as the research object to study its influence on the movement law of particles in PCL. Numerical setting is the same as standard setting except the depth of PC whose value is set as 6.7, 7.2, 7.7μm at three different particle initial positions on Y direction (4.5, 5, 5.5μm from the bottom of the channel). The specific trajectory of this study is shown in Figure 3. According to the experimental results, if the depth of PCL increased, the final equilibrium position of the particles in Y direction increases. In addition, the final equilibrium position of the particles will increase as the initial position in Y direction increases.
With different oscillation amplitude of cilia, particle's movement feature in different PCL positions also studied. It is worth mentioning that in order to increase the amplitude of cilia beat, the distance between cilia is increased to 0.93L (L is the length of the cilia). Figure 4 shows the motion pattern of  figure 4(d), when the initial position in Y direction is 5μm, particles move upward in both beat modes. For particles with a lower initial position, such as 3.5μm (figure 4(e)~ figure 4(h)), its position decreases slightly or remians unchanged with small amplitude of cilia, while if ciliary oscillation amplitude is about 120°, most of the particles move upward. Figure 5(a) shows the result of the ciliary oscillation amplitude is about 90°, and Figure 5(b) is about 120°. According to the experimental results, if the ciliary oscillation amplitude is about 90°, when the initial position of the particle in Y direction is greater than 4.5 μm, the particle will move up to the junction of PCL and ML under the action of cilium beat. When the initial position of the particle is less than or equal to 4.5 μm, the position of the particle in Y direction will decrease slightly or remain unchanged. When the ciliary oscillation amplitude is about 120°, all particles tend to move upward. In addition, the average position in Y direction of uppermost particles with the beat amplitude of 120° is lower compared with 90°. What's more, it seems that particles aggregation extent is higher when ciliary beat amplitude is about 120°.

Conclusion
In the present study, a 2D numerical model is proposed to explore the migration law of particle in the PCL of human respiratory tract. It is found that driven by cilium beat, the migration of particle is closely related to its initial position and the beat amplitude of cilia. Firstly, when putting the particles initially near the top of cilium, they tend to migrate to the mucus layer. Secondly, if the initial positions are located near the middle and lower part of the cilia, most of the particles will not migrate to the vicinity of the mucus layer. Thirdly, if increasing the beat amplitude of cilia, particles will have more opportunities to approach the mucus layer. Therefore, this indicates that larger oscillation amplitude of cilia has a bigger impact on removing the solid particles in PCL.